最小二乘法拟合直线

2024-05-21

最小二乘法拟合直线（精选7篇）

最小二乘法拟合直线篇1

1 引言

MLS法是对传统加权最小二乘法的进一步推广, 克服了最小二乘法的不足, 突出了加权函数的紧支性, 据此建立的拟合函数, 具有局部近似的特点。MLS法的基本理论本文不再赘述详见参考文献[4]。本文针对具体二维和三维算例分析了MLS的局部近似性, 为其他研究领域中引入该方法解决复杂问题提供了一种新的思路。

2 算法流程

结合移动最小二乘法的特点, 给出了基于MATLAB软件编写的MLS法的曲线曲面拟合程序, 程序设计流程如图1所示。

3 权函数的影响分析

由图2所示:权函数为恒定常数, 样本点的权函数值均为1, 在求解过程中系数a (x) 不再是自变量x的函数, 而是在整个区域内a (x) 为常数, 可以看出系数a (x) 与权函数的选取有关, 因此, 曲线拟合就退化成最小二乘法的曲线拟合。这样权函数形成的近似为线性拟合。

由图3所示:权函数为紧支常数, 只有在中心点x (28) 5.0权函数存在非零区域, 非零区域尺寸即为节点权函数支持域尺寸, 这样在任意一点x的支持域内均有两个节点对其有贡献 (n (28) m (28) 2) , 权函数为1, 其他节点的权函数为零, 权函数构造的近似就变成了分段线性插值, 类似有限元法的函数逼近。由图4所示:权函数为紧支且光滑的函数 (三次样条函数) , 计算点x的支持域内包含的节点数大于基函数的个数 (n (29) m) , 为使局部近似误差最小, 加权平方和J (x) 对a (x) 求偏导, 最终得到MLS的近似位移, 这样的求解过程用于EFG方法中, 即使是线性基函数, 近似值也是连续光滑的, 因为MLS近似继承了权函数的连续光滑性。因此, 本部分的分析, 选取三次样条函数。

4 曲线拟合

采用常数基函数]1[, 拟合结果如图5所示。

采用一次基函数[1x], 拟合结果如图6所示。

采用二次基函数[1xx2]拟合结果如图7所示。

由图7可以看出, MLS拟合函数及其一阶导数拟合的效果比较好, 二阶导数拟合呈现稍微的震荡。为了更加精确的分析拟合结构, 引入相对误差范数公式:

精确解公。式中:N为样本点个数, yN和yE分别为对应的数值解和

在移动最小二乘法中基函数的选取也有一定的规律, 从以上分析, 结合表1可以看出:同一个函数采用不同基函数, 对计算结果有直接的影响, 线性基的基函数的次数越高, 逼近效果越好且数值计算的相对误差也越小, 二次基函数能够得到良好的拟合效果。

5 曲面拟合

采用MLS法拟合以下函数:

设定区域在区域W上均匀布置8 (8个节点, 图8显示了精确值的曲面分布图, 利用MATLAB中Griddata函数插值得到的。为进行曲面拟合, 采用具有C0 (W) 连续的高斯型权函数 (b=30.) 和线性基函数[1xy], 在区域上均匀划分25 (25个规则的节点网格, 由MLS法公式, 求出各节点的近似函数值, 绘制网格曲面, 如图10所示, 从图中可以看出, 用移动最小二乘法拟合的曲面比较光滑, 并与Griddata插值形成的曲面相比, 如图9所示, 曲面光滑准确, 从而验证该方法的有效性。

为了衡量拟合的质量, 定义拟合的相对误差:

其中:zEXACT (x, y) 是精确解, zMLS (x, y) 是使用移动最小二乘法计算值, 本文得到的拟合误差为:L=2.1526, 这说明MLS拟合具有比较高的精度, 虽然MLS拟合的计算量稍微大一些, 和其精确性相比, 是可以接受的。

6 结语

移动最小二乘法虽然不具备传统有限元或边界元形函数所具有的插值特征, 即不满足Kronecker函数的属性, 给边界条件处理带来了困难, 但是它采用了较低的基函数通过选取适当的权函数, 获得具有较高连续性和相容性的形函数, 数值精度较高, 数据简单等特点, 是其他数值近似方法无法比拟的, 对以后的数值近似的研究具有重要的参考价值。

参考文献

[1]王能超.数值分析简明教程[M].北京:高等教育出版社, 2003.

[2]常青, 余湘娟.软土次固结变形特性的试验研究[J].南水北调与水利科技, 2008 (1) :148-150.

[3]张荣.利用VB实现水泵性能曲线的自动绘制[J].水利科技与经济, 2006 (1) :62-64.

[4]程玉民.移动最小二乘法研究进展与述评[J].计算机辅助工程, 2009 (2) :5-11.

最小二乘法拟合直线篇2

关键词：最小二乘法,曲线拟合,Matlab

1.用MATLAB求解拟合问题

最小二乘法做数据拟合在理论上已经介绍,然而实际计算过程较为烦琐,目前经常用数学软件处理计算问题,常用的有Mathmatical,Matlab等,这里我们用Matlab软件做最小二乘拟合法处理实际计算问题.

1.1用Matlab作线性最小二乘拟合[1,2,3]

(2)多项式在x处的值y的计算命令:y=polyval(a,x)

1.2用Matlab作非线行最小二乘拟合

两个求非线性最小二乘拟合的函数:lsqcurvefit、lsqnonlin.

2.实际应用实例

下表为吐鲁番市某年1月1日的部分气象数据见表1,为找出干球温度及地表温度变化规律,需要找出相应的数据拟合曲线.

这里用Matlab做最小二乘拟合.作出干球温度与地表温度的逐时温度值的曲线图,并进行观察.容易发现,散点图分布情况接近多项式函数,可以结合Matlab编程,找出拟合曲线.

输出图像为:(加号为地表温度曲线,圆点为干球温度曲线)

结束语

由图1和图2进行对比,图3和图4进行对比,可以看出,在用最小二乘法的Matlab曲线拟合时,若拟合函数为多项式类,需要适当调整最高次数,才能找到最佳的逼近曲线.根据曲线图像与原始数据进行对比,可见,曲线并没有经过所有原始数据所在点,从这里可以看出,用最小二乘法做曲线拟合与函数插值的区别,插值法要求过所有实验数据点,最小二乘法只需求出最佳拟合曲线.这样就可以避免少数不符合规律的实验数据被使用的问题.

参考文献

[1]高照玲,周浩尚.蒋波VC++6.0实现计算方法中的曲线拟合[J].农业装备与车辆工程,2011(11).

[2]王可,毛志伋.基于Matlab实现最小二乘曲线拟合[J].北京广播学院学报(自然科学版),2005(02).

最小二乘法拟合直线篇3

工程施工中, 我们会经常取得一些相关的数据, 这些数据往往来自与施工密切相关的测量或试验中, 比如石灰剂量的滴定试验中, EDTA用量与灰土剂量存在某种关系;又如预应力千斤顶与油表的配套校验中, 油表读数与千斤顶实际张拉力又有一种关系, 这些原始数据一般是5组以上。

我们可以通过作图或多段插值取得变量之间的联系, 但作图和插值查图往往误差较大, 这时可采用最小二乘法先拟合出一个多项式 (也称经验公式) , 再根据此多项式求解任一自变量所对应的因变量较精确的结果, 据此绘图可得到较精确、较合理的曲线。

2 最小二乘法原理

现实中通过测量或试验取得的各组数据 (xi, yi) 其本身不可避免的地带有测试误差, 如果构造一个较为简单的插值法P (x) 来逼近真实函数f (x) , 当个别数对误差影响较大时就会引起插值函数发生严重波动, 从而影响逼近精度, 因为插值法要求插值函数在节点 (已知数对) 处满足标准条件, 即P (xj) =f (xj) , (j=0, 1, …, n) 。这时候, 为尽可能减小测试误差对逼近精度的影响, 我们可以用另一种方法构造一个经验公式, 使得该公式在每一个节点上所求得的结果与原测试结果的差的平方和最小, 即曲线拟合的误差最小, 精度最高, 这就是最小二乘法原理, 用定义表述为:

设有n对数据 (xi、yi) , (j=0, 1, …, n) , 通过这些数据找一个m次近似多项式P (x) =a0+a1x+…+amxm (m

使得为最小值,

则称P (x) 为最小二乘拟合多项式, 或称x、y之间的经验公式。

3 最小二乘原理应用的具体做法

根据定义, 我们最终目的是根据已知数据组求出合适的系数ai, 使多项式φ (ai) 取得最小值, 运用多元函数求极值的方法可导出如下方程组:

这是以ai为未知数的m+1阶线性方程组, 写成矩阵形式有

因此, 我们只要求出S0, S1, S2, …, Sm, …, S2m以及T0、T1、…、Tm, 建立正规方程组即可求出系数ai, 代入多项式P (x) 即得拟合曲线方程。

4 应用实例

4.1 灰剂量EDTA滴定法标准曲线绘制备5种试样共10个样品的灰土混合料, 其石灰剂量为已知, 在6%~18%范围内按3%剂量递增, 该范围已涵盖2:8灰土到3:7灰土 (V:V) 的石灰剂量范围。分别对样品进行滴定, 测得EDTA标准液耗量 (ml) 如下:

步骤1:作草图, 先将数据点画在坐标纸上, 连线分析其变化的粗糙趋向, 确定多项式的大致形式, 可以看出, 本例为一缓和抛物线。

步骤2:设拟合曲线方程为y=a+bx+cx2, 列表求和计算Sk、Tk (见下表)

步骤3:建立正规方程组

解得c=-0.0013b=0.382a=0.253

因此得方程y=0.253+0.382x-0.0013x2

该经验公式求得的结果与测试数据的最大离差Δ=-0.2, 离差的平方和最小为0.06, 据此可绘制石灰剂量 (%) 与EDTA耗量 (ml) 关系的标准曲线, 应用中当已知某灰土样品液滴定的EDTA耗量x时, 直接代入公式或查标准曲线便可得到较为精确的待测灰土的石灰剂量, 方便、快速而又准确。

4.2 求解张拉千斤顶与油表读数的回归方程预应力千斤顶与油表的配套校验中, 分级张拉数据可达到5~20组, 而张拉力与油表读数实际为线性关系, 一般只需两组数据便可确定其关系式, 但数据越多, 回归方程越真实, 越精确。此时采用最小二乘法可使每一组数据参与回归。

对线性方程y=ax+b

用最小二乘法公式求系数a、b (斜率和截距) , 公式简化为:

某千斤顶校验数据见下表 (n=7)

由于第一数据 (0, 0) 将可能使计算失真, 误差较大, 回归时舍去该组数据, 取n=6组, 计算得:

则有a=0.035 b=0.68,

故得经验公式y=0.035x+0.68

回归相关系数r=0.99993, 精度较高。应用时将所需张拉力的值代入公式便可得到相应张拉力下的油表读数, 以控制张拉力。

5 结语

最小二乘法作为函数逼近的一种重要方法, 在工程技术中的曲线拟合、求取经验公式等方面有着广泛的应用, 但一般来说拟合式不高于三次多项式, 否则在计算上将十分繁琐, 且拟合精度将受到影响甚至出现病态。而实际工程应用中拟合的曲线多为一次线性或二次抛物线型, 因此, 最小二乘法拟合曲线在工程应用中仍占有着重要的位置。

摘要：施工中常遇到一些相关数据的分析处理, 并要求拟合曲线以便反映数组规律和扩大应用范围, 本文介绍了目前应用较广、较为精确的曲线拟合方法, 以及回归方程的建立。

关键词：最小二乘法,拟合曲线,施工,应用

参考文献

最小二乘法拟合直线篇4

图像质量评价(Image Quality Assessment,IQA)属于信号处理学科中较为年轻的领域,分为主观评价和客观评价两种方式。主观评价方式讲究样本采样,以样本分布模拟总体分布,把参与评价的主体的平均打分作为图像质量主观评分(Mean Opinion Score,MOS)[1,2,3]。主观评价因为评价周期长、参与人数多、对主客观条件均要求严格,不利于在实践中推广。与之相反,客观评价方式借助高性能计算机和基于人类视觉模型的算法设计,不仅能快速给出始终一致地图像质量评分,而且方便嵌入其他工程应用[4]。近些年,研究者们对IQA客观算法的设计与改进取得了不少优秀成果。然而,针对客观算法的准确性评估一直都存在争议。

Richard Dosselmann[5]对IQA算法中获得广泛应用的结构相似度[6](Structural Similarity index,SSIM)算法与均方误差(Mean Squared Error,MSE)算法进行统计分析和主观比较研究之后发现,上述两种算法的性能有很多相近的地方。这与文献[6]和文献[7]的结论不同。此后,他与Xue Dong Yang合作[8]为SSIM算法和MSE算法的度量值建立了代数关系式,更进一步拉近了两种算法之间的关系。以上成果提供了审慎反思的新视角,它对IQA算法中仅凭少数个例,或通过缺乏可信度的统计分析而得到的结论提出了质疑。为了给不同的IQA算法做准确评估,Hamid Rahim Sheikh等[7]和Nikolay Ponomarenko等[9]采取了相同的策略,他们分别设计了LIVE database图片库和TID2008图片库来代表各种图像内容和失真类型,以对应图像的MOS值作为评估IQA算法的依据。不同的是,Sheikh[7]提出首先用1个含5个参数的Logistic方程做数据的非线性回归拟合,然后选用3种性能度量指标分别做评估;不足之处在于Logistic方程缺乏普遍意义,以及3种性能指标的统计学意义含重复的部分;Nikolay Ponomarenko[9]则直接选用Spearman相关系数和Kendall相关系数对各算法分别进行排序,这2个统计指标显得过于粗糙、简陋。

视频质量专家组(Video Quality Experts Group,VQEG)作为推进视频质量评价VQA(Video Quality Assessment),由于VQA与IQA的评估标准相同,本文中一律用IQA代替)领域的专家组织,前后推出3次针对视频质量评价算法评估的报告[1,2,3]。在统计分析策略方面,新报告相比之前版本均有变动。最新版本报告[1]的策略为:首先用1个保证单调递增的含4个参数的三次多项式方程做为数据的非线性回归拟合;然后选用3种性能度量指标分别从不同角度做评估。该报告未提供三次多项式方程形式的非线性回归拟合的具体实施过程。

针对VQEG提供的IQA算法评估策略中未作说明的,三次多项式形式的非线性回归拟合过程开展研究,力争探索一种客观评估IQA算法性能的方法。本文选用TID2008图片库[9]及其MOS数据作为实验素材,以峰值信噪比(PSNR,Peak Signal to Noise Radio)、SSIM[6]和多尺度SSIM[10](MSSIM,Multi-scale SSIM)3种算法为评估的对象举例,给出了基于有约束最小二乘法的曲线拟合过程。

1 对IQA算法的评估过程

主观评价被认为是评价图像质量最有效力、最可信赖的方式[11]。因此,本文认为对IQA算法做评估的过程就是衡量IQA算法度量值与主观评分的相似度和差异度的过程。本文仅对图像逼真度(Image Fidelity)感兴趣,因此只涉及全参考IQA算法的质量评估。其中所采用的三次多项式形式的回归映射函数和3个性能指标,均参照VQEG[1]的做法;考虑到TID2008数据库中MOS值的获取途径,原始图像(或参考图像)的MOS值存在缺失,因此,将文献[1]中出现DMOS值的位置用失真图像的MOS值代替。

1.1 回归映射函数

回归映射函数的方程形式为

MOSp=ax3+bx2+cx+d (1)

式中:要求式(1)在x定义域范围内保持单调递增。

1.2 评估的三个性能指标

1.2.1 Pearson相关系数(R)

$R = \frac{\sum_{i = 1}^{Ν} (X_{i} - X) \cdot (Y_{i} - Y)}{\sqrt{\sum (X_{i} - X)^{2}} \cdot \sqrt{\sum (Y_{i} - Y)^{2}}} (2)$

式中:Xi代表主观评分(MOS值);Yi代表客观评分(MOSp值);N代表参与评估的图像总数。Pearson相关系数R衡量了一种算法输出的预测值与主观数据之间的线性关系,值越大,表示算法的单调性越好。

1.2.2 均方根误差(RMSE)

$R Μ S E = \sqrt{{\frac{1}{Ν - d} \sum_{Ν} [Μ Ο S (i) - Μ Ο S p (i)]^{2}}} (3)$

式中:N代表参与评估的图像总数;d表示映射函数方程中自由度的个数,文中d=4。RMSE描述了IQA算法的准确度,值越小,算法的准确度越高。

1.2.3 离散率(Outlier Radio,OR)

$Ο R = \frac{Τ o t a l Ν o Ο u t l i e r s}{Ν} (4)$

式中:1个有效的离散点是满足式(5)约束条件的点。

$| Μ Ο S (i) - Μ Ο S p (i) | > Κ_{2} \cdot \frac{σ (m o s (i))}{\sqrt{Ν s u b j s}} (5)$

式中:K2=1.96, $σ (m o s (i))$ 表示与第i幅图像相关的标准差;Nsubjs=33,表示参与评价第i幅图像的人数[8]。OR表征了IQA算法评分与主观评分的一致程度,值越小,算法与主观评分的一致性越好。

2 有约束线性最小二乘法的数学模型

考虑到IQA算法的目的是为了模拟人眼视知觉的判断,本文选定将均方根误差最小化作为回归映射函数的直接目标。将IQA算法应该满足的基本假设为:图像质量越高,即MOS值越大,算法度量值(MOSp)应该越大。因此,在定义域内保持单调递增是回归映射函数需要满足的约束条件。因为均方根运算不影响函数的单调性,所以目标函数的形式可简化为

$\min f = \sum_{Ν} (Μ Ο S (i) - Μ Ο S p (i))^{2} (6)$

将公式(1)代入式(6),得

$\min f = \sum_{Ν} (Μ Ο S (i) - a \cdot x^{3} - b \cdot x^{2} - c \cdot x - d)^{2} (7)$

式中: f为因变量,对任何i值,MOS(i)均为常量。以上述形式充当目标函数的方法即为最小二乘法。

当将x作为自变量考虑时,x的定义域可归一化为[0,1]区间,显然公式(1)为连续函数,为了使MOSp(x)在x的定义域内满足单调递增,当且仅当MOSp(x)的一阶导数MOSp′(x)在[0,1]区间满足

MOSp′(x)=3ax2+2bx+c≥0 (8)

式中:目标函数和约束不等式十分复杂。参考文献[12,13,14,15]中提及该目标函数和约束条件既不属于多元线性回归问题,也不属于简单的曲线拟合和非线性回归的范畴。从回归分析的角度来看,其可定性为有约束的线性回归问题,此时a,b,c,d为待求变量。

当将a,b,c,d看作自变量重新考虑上述过程时,发现公式(8)即为线性约束条件。离散情况下,当x的数据量足够大时,若对每一个x值,均有公式(8)成立,那么可近似认为MOSp(x)在[0,1]区间上单调递增。当然,前者为后者的必要不充分条件,结果是否符合预期需要验证。简言之,对第i幅图像(i=1,2,…,N),使每一个MOSp(xi)均满足公式(8),同时最优化目标函数(7)的值最小,最后通过查看映射函数曲线来反过来验证函数MOSp(x)是否单调。上述即为本文采取的策略,可用数学模型表达为

${\begin{cases} \min_{x} ∥ C x - d ∥_{2}^{2} \\ A \cdot x \geq b \\ l b \leq x \leq u b \end{cases} (9)$

式中:C和A为矩阵;d,b,lb,ub和x为矢量,其中x=[a,b,c,d]为待求变量。

3 实验结果与分析

实验平台采用Matlab实现,其数学模型为

${\begin{cases} \min_{x} \frac{1}{2} ∥ C x - d ∥_{2}^{2} \\ A \cdot x \leq b \\ A e q \cdot x = b e q \\ l b \leq x \leq u b \end{cases} (10)$

式中:C,A和Aeq为矩阵;d,b,beq,lb,ub和x为矢量,其中x为待求变量。

对本文而言,各参数的含义如下:x=[d,c, b,a]T,A=[0,-I,-2y,-3y.^2]1 700×3,y为某一IQA算法1 700×1维的度量值,I=[1,1,…,1]T1 700×1,0=[0,0,…,0]T1 700×1,b=[0,0,…,0]T1 700×1,C=[I,y,y.^2,y.^3]1 700×3,d=[mos1,mos2,…,mos1 700]T。其中mosi表示对应第i幅图像的MOS值,Aeq=beq=[],lb=[-Inf,-Inf,-Inf,-Inf],ub=[Inf,Inf,Inf,Inf]。上述符号均按Matlab语言表述,由此解得的自变量x即为公式(1)中映射函数的系数。

根据PSNR、SSIM和MSSIM三种算法基于TID2008数据库的计算结果,使用Matlab优化工具箱拟合得到的回归映射函数分别为

MOSp(PSNR)=5.251 5-16.686 1·psnr_t+

38.142 1·psnr_t.^2-19.070 9·psnr_t.^3 (11)

MOSp(SSIM)=2.773 8-0.423 2·SSIM+

0.488 9·SSIM.^2+3.010 9·SSIM.^3 (12)

MOSp(MSSIM)=-46.031 6+195.219 8·MSSIM-

267.960 4·MSSIM.^2 124.763 6·MSSIM.^3 (13)

式中:psnr_t=psnr/50,使得psnr_t值中99.94%的比例落在区间[0,1]之间。

PSNR、SSIM和MSSIM 3种算法的散点图和回归映射函数曲线如图1所示。图1的横坐标分别表示图像的PSNR、SSIM和MSSIM算法归一化到[0,1]区间的度量值,这些度量值通过对TID2008图像库的失真图像及其参考图像应用各算法计算得到;纵坐标表示图像的MOS值,数据直接取自TID2008图像库[9],0值表示主观判断图像质量最差,9表示最好。如果存在理想算法,那么在图1的坐标系下其散点将全部分布在第1象限的某条直线上。观察图1中曲线可知,根据此方法得到的映射函数的确在数据定义域内单调递增。

实验所得的统计参量如表1,表2所示。表1为不同算法的统计参量绝对值。该组数据显示,3种算法与MOS值之间的线性相关系数R全部低于90%,RMSE平均相差0.7个等级(全部为0～9共10个等级),OR系数显示有至少94%以上的算法度量误差阈值大于人眼。表2列出了不同算法的相对差异程度,用0和1表示,0表示两种算法的差异不明显,1表示两种算法间存在显著的差异。由以上两组数据可以得到3种算法的性能排序,但是其结果与人眼判断的MOS值差距很大,均不能令人满意。这组性能指标与文献[6,7]的实验数据相比,数据意义简单、明确,对客观地理解和判断IQA算法的性能有所帮助。

上述实验结果表明,目前主流的IQA算法与主观评价值之间差异较大,图像的客观评价算法仍有待进一步提高性能。本文提供了一种衡量图像客观评价与主观评价之间差异的方案,为图像评级研究提供了帮助。应该指出,本文的方法还有待完善,如公式(1)不一定满足在[0,1]整个区间内单调递增,具体单调区间受算法的度量值x的上下限xmin、xmax影响,只能满足在[xmin, xmax]区间内的单调递增,下一步研究将对该曲线拟合方法继续完善。

4 结论

本文针对VQEG对IQA算法的评估策略,采用基于有约束最小二乘法的数学模型,对算法数据的三次多项式形式非线性回归拟合过程求解。建议的方法约束条件为线性,所以求解过程效率高,可以为研究IQA算法提供性能比较的实验平台。

最小二乘法拟合直线篇5

在现代高铁、地铁隧道矿山法施工中, 人们对工程成品结构的精度要求越来越高, 这就要求有与之相适应的精密测量技术。在浇筑成品混凝土所用钢模台车模板精密测量中, 圆或圆弧的几何特征参数的测量及评定是最基本、也是最重要的测量内容。

在施工中, 模板台车在没有任何线性参数的情况下, 也就是我们常说的没有平曲线和竖曲线参数, 需要进行检验台车拼装精度或变形状态, 其主要是进行模板台车几何参数评定, 对相对尺寸进行解析, 进而是对空间位置的圆心和半径的求取, 其方法就是测得表面上若干个点, 然后拟合成平面圆, 再求出圆心。本文用最小二乘法拟合圆应用于台车尺寸检查与传统方法进行比较, 具有操作简单, 速度快, 线性整体性好, 精度高等特点。

1 传统方法

首先对台车模板某一横断面进行数据采集, 假设测点在纵向影响为零, 将采集数据通过Auto CAD展点生成一个横断面, 任意连接半径相同的3个测点可确定一个圆, 由于测量误差、模板台车加工误差或模板台车本身变形, 各个圆心并不重合, 将再任意选取前部分得到的3个圆心点进行确定一个圆, 以此类推, 满足精度要求为止。这样确定圆心后根据标准断面设计半径画出圆, 将各个测量沿径向量取偏差值, 此偏差值作为台车的变形或检查成果。

2 用最小二乘原理拟合法

在台车的几何参数评定中, 主要是圆心和半径的求取, 其方法是在某一横断面上测取若干坐标, 然后拟合出平面圆, 从而求出圆心和半径, 再进行偏差值的对比。如图1:

其中Ri为测点半径, 是测点位置的相角。此时, 测点的直角坐标 (Xi, Yi) 为:

其中 (i=1, 2, 3…, m)

根据文献1可得出近似二乘估算值为:

上述算法测点需均匀分布于断面上, 且需要满足两个条件, 即小偏差假设和小误差假设。然而在实际中, 加工误差、拼装误差和测量误差使之不能实现, 处理数据也比较繁琐, 由此可用直接依赖测量数据求取任何平面圆的圆心和半径的方法。解法如下:

其中 (i=1, 2, 3…, m-1) 。

其中 (i=1, 2, 3…, m-1) 由文献2可得到圆心坐标:

式 (9) 、 (10) 、 (11) 为圆曲线拟合的圆心坐标和半径的通用算法。

3 实例分析

为了验证两种方法的有效性和适应性, 本题依托在建工程深圳地铁安托山停车场出入线矿山法隧道E型断面进行分析。如图2:

E型台车断面由5个半径的圆构成, 测点较均匀的分布在断面上。如图3。

假如O为原点, X方向为横轴, Y为纵轴, 下表为已换算为当前坐标系下的测量数据, 如表1:

通过将数据进行计算机处理, 得出两种方法对比成果, 如表2:

注:“/”为此半径圆上测点太少, 作图法不能得到相应值。

3 结语

评定钢模台车的相对几何参数时, 通过传统作图法与最小二乘原理拟合圆成果对比表明, 作图法在每个不同半径圆上测点必须大于等于3个, 测点图形强度要求高, 作图繁琐, 误差较大, 极易产生人为误差, 不能进行台车的整体线性判定;而后者在满足小误差和小偏差条件下, 可用 (5) (6) (7) 直接计算, 计算过程简便, 如果不满足小误差和小偏差条件, 用直接依赖测量数据用式 (9) 、 (10) 、 (11) 进行计算, 虽计算量大, 但借助计算机能较快实现。总之, 最小二乘原理拟合圆应用于隧道钢模台车尺寸检查具有计算过程严密, 借助计算机计算速度快, 效率高, 避免了不必要产生的误差, 为计算机编程奠定了理论基础, 能在现场实时测量计算偏差值, 更快捷指导现场施工。

摘要：在高铁、地铁隧道矿山法施工中, 常用钢模台车浇筑洞身二次衬砌, 在此过程中, 不时会发生台车变形, 运用最小二乘法拟合圆应用于台车尺寸检查比传统方法具有操作简便, 速度快, 线性整体性好, 精度高等特点。

关键词：变形,最小二乘法,操作简便,速度快,精度高

参考文献

[1]熊有伦.精密测量的数学方法, 中国计算出版社, 1989;

[2]田社平, 张守愚, 李定学等.平面圆圆心及半径的最小二乘拟合, 1995;

[3]徐国旺, 廖明潮.拟合圆的几种方法[J].武汉工业学院学报, 2002 (4) :104~105;

最小二乘法拟合直线篇6

强震后余震是大地震发生后的常见灾害之一,震后余震对救灾活动、难民安置以及灾后重建等工作有着严重危害。所以,对震后余震的分析及预测就显得尤为重要。

文献[1]在研究了我国20世纪30年代到70年代接近50年内发生的28次地震的强余震序列后,发现了“等待时间”与余震“发生时刻”的对数线性关系,并对其作了初步理论解释。文献[2]分析了广州邻近区域的几次典型的余震序列,发现余震“发生时刻”的对数与“等待时间”对数具有较好线性关系,进一步印证了文献[1]的理论。文献[3]系统地分析了1966~1996年间70余次Ms5.5以上的地震的余震序列,发现强余震“等待时间”模型对6级以上主余震型序列预测成功率较高。文献[4]对多次大震和中强震进行比较研究,发现“等待时间法”对大震的强余震预测效果较好。文献[5]利用混沌理论对汶川地震余震“等待时间”序列进行研究,发现“等待时间”序列存在混沌动力学特征。文献[6]以“等待时间”作为权值弥补双对数变化对余震预测误差的影响,发现经过加权改正后,模型的预测精度有所提升。然而,现有的利用“等待时间法”进行余震预测的研究都是基于文献[1]的理论,仅考虑了“等待时间”的随机误差,运用最小二乘法(LS)或加权最小二乘法对文献[1]的模型进行求解,未能顾及“等待时间”模型中由“发震时刻”构成的系数矩阵的随机误差。

本文根据“等待时间”与余震“发生时刻”的对数线性关系,利用总体最小二乘方法(TLS),同时考虑“等待时间”观测向量的随机误差和以“发震时刻”构成的系数矩阵的随机误差来进行拟合;并将该方法应用在尼泊尔地震的强余震序列预测和拟合中,得到了有益的结果。

1 利用总体最小二乘法拟合的“等待时间”强余震预测模型

1.1“等待时间”强余震预测模型

连续发生的两个强余震之间的时间间隔为后一个强余震的“等待时间”。即:

文献[1]发现“等待时间”与“发震时刻”存在以下对数关系:

在式[1]及式(2)中i为第i次强余震的序号,ti为第i次强余震“发生时刻”(距主震时间),Δti则表示第i次强余震的“等待时间”,以上时间单位均为小时。β0,β1为拟合参数。式(2)中的β0,β1值可通过总体最小二乘方法拟合已知的强余震序列求出。根据得到的β0,β1值可以得到拟合直线,再通过最近一次余震的“发震时刻”对数即可预测出下次余震的时间对数,进而得到下次余震的“发震时刻”。

顾及观测向量和系数矩阵同时存在的随机误差的情况,建立EIV模型如下:

式中y为观测向量值,ey为观测值误差:

A为系数矩阵,EA为系数矩阵误差:

ξ为未知参数:

总体最小二乘的目标函数:

式(8)中vec为矩阵拉直符号。

利用总体最小二乘方法求解的迭代算法可表示如下[7]:

[1]以最小二乘方法得到的为起始点ξ^[1]:

(2)计算出残差:

(3)将式(9)、(10)得到的结果代入式(11)求出:

(4)重复步骤(2)、(3),当:

时终止迭代,本文中取ε=10-10。

精度评定:

1.2 作穷举线求交点

根据总体最小二乘方法拟合得出的β0,β1作拟合直线,并将求出的直线同穷举线作交点,穷举线绘制是依照以下公式进行的:

穷举线与拟合线的交点即为预测的第N+1个强余震,交点的横坐标即为预测下次强余震“发生时刻”的对数值,经过简单变换即可以得到预测时间[6]。

2 以尼泊尔地震为例的强余震序列拟合和预测

2015年4月25日在尼泊尔博拉克(北纬28.2°,东经84.7°)发生的Mw7.8强震给当地人民造成了巨大的生命和财产损失,接踵而来的余震也对救援和安置工作造成了很大的困难。本文分别应用最小二乘与总体最小二乘方法对尼泊尔地震强余震序列进行拟合及预测,并对结果进行讨论分析。

2.1 数据选择

对于余震目录的筛选,在利用“等待时间”序列进行余震预测时,一般采用主震震级减2~3级为余震目录的震级下限,并把时间间隔小于前一强余震“等待时间”的1/20的余震与前一次余震当作同一次能量释放,可以把它从余震目录中剔除[8]。

本例采用的余震目录是由美国地质勘探局(USGS)在其官网提供的尼泊尔Mw7.8级地震5级以上余震序列,并剔除了时间间隔小于前一强余震“等待时间”的1/20余震序列,时间截至2015年4月27日12:35:53(UTC时)。余震目录见表1。

2.2 震例分析

为了验证总体最小二乘方法的拟合效果,本文首先应用最小二乘法和总体最小二乘法对表1中的前17个样本进行拟合,利用总体最小二乘法得到β0=-0.638,β1=1.087,σ02=0.054。利用最小二乘法得到β'0=-0.612,β'1=0.980,σ'02=0.132。以两种方法求得的单位权中误差σ02和σ'02作为拟合线残差。通过求得的参数即可绘出拟合直线,见图1。

根据式(14)作穷举线,穷举线与拟合直线的交点即为预测的第N+1个强余震。预测残差是下一次强余震的“发生时刻”对数值的预测值和实际值之差。分别用eLS和eTLS来表示最小二乘与总体最小二乘的预测残差。结果见图2。

为了体现样本个数对拟合和预测结果的影响,本文进行了如下实验:采用前n个样本预测第n+1次余震,分别利用最小二乘和总体最小二乘方法进行直线拟合并预测下一次余震,样本数从11个逐次增加到19个,进行9次实验,结果见表2。

通过表2及图3可以看出,在该震例中,总体最小二乘对于余震“等待时间”和“发震时刻”对数模型的拟合结果要优于采用最小二乘拟合方法,预测精度在总体上也比直接采用最小二乘有较好的结果。预测残差的绝对平均值由0.082降为0.071,降幅约为13%。拟合残差的绝对平均值由0.071降为0.059,降幅约为17%。图3中随着样本数的增加,两种方法的预测误差都向0值靠拢,而总体最小二乘的结果在整体上更加接近0值。

3 其它地震实例

为了进一步验证总体最小二乘方法应用于余震“等待时间”和发震时刻对数模型的拟合和预测效果,本文另选了几组国内外典型强震的强余震序列样本进行最小二乘方法与总体最小二乘方法对比实验,结果见表3及图4。

由表3和图4可以看到,利用总体最小二乘方法得到的结果比直接运用最小二乘方法得到的结果在拟合精度上有明显的优势,在预测精度方面也有着不同程度的提高。总体上,拟合残差绝对均值从0.133降低为0.053,降低幅度约为60%,预测残差绝对均值从0.170降低为0.148,降低幅度约为13%。其中日本311地震精度提高幅度最大,拟合残差绝对均值由0.301降低为0.078,预测残差绝对均值由0.101降低为0.066,降低幅度分别为74%和35%。

4 总结

本文提出了同时顾及“等待时间”观测向量的随机误差和以“发震时刻”构成的系数矩阵的随机误差,采用总体最小二乘方法来进行强余震“等待时间”序列的拟合,并将其应用到2015年尼泊尔地震的强余震预测中。再经过对10个典型震例的实验分析后,结果表明,总体最小二乘方法可以用于“等待时间”余震模型的拟合以及余震预测;在拟合方面,因为总体最小二乘方法顾及了系数矩阵的误差,所以较最小二乘而言有较好的提升;在余震预测精度方面也有一定的提高。但是由于本文所涉及的余震震例较少,可能未能全面反映总体最小二乘在该方面的优劣性,利用该方式进行实际的余震预测有待进一步研究。

摘要：本文针对“等待时间法”余震预测模型在拟合过程中未考虑系数矩阵的随机误差的问题,利用总体最小二乘方法进行该模型的拟合,并将该方法应用于尼泊尔地震及其它震例的强余震预测。结果表明,运用总体最小二乘方法进行的“等待时间”余震序列拟合和预测所得结果精度比传统的方法在总体上有一定的提升。

最小二乘法拟合直线篇7

传统的分段曲线拟合根据主观经验和绘制数据散点图来确定拟合的经验函数和分段点。文献[5] 提出分段区间重合的拟合方法, 由每4个数据点决定一个三次曲线, 但分段区间太密, 不适用于密集的数据拟合。文献[6]提出的多项式基函数的全局连续拟合方法, 只限于2个分段区间。文献[7]提出多分段区间全局连续的曲线拟合方法, 但基函数只限于一次多项式。文献[5—7]提出的方法都是根据主观经验来确定经验函数和分段区间, 然后进行拟合。文献[8]引入拟合误差限度 εmax, 在误差大于该限度的点处分段, 但该限度 εmax不容易界定。

提出一种自动分段多项式曲线拟合方法, 该方法在实际运用中具有如下有点: 1提供几种不同的经验函数, 根据不同经验函数拟合的数据和实测数据的误差的方差, 自动确定较优的经验函数; 2自动确定数据的分段区间, 在各个区间进行最小二乘拟合。通过数值实验, 证明该方法有较高的拟合精度。

1经验函数的确定及拟合步骤

1. 1经验函数的确定

数据拟合方法有很多, 例如对数曲线拟合, 反函数曲线拟合, 二次曲线拟合, 三次曲线拟合, 幂函数曲线拟合, 指数曲线拟合等。一般先观察散点图来确定曲线的类型, 不过散点图都是相关关系的粗略表示, 有时候散点图可能与几种曲线都很接近, 这时建立相应的经验函数可能都是合理的, 但由于选择不同的曲线, 得到同一个问题的多个不同经验函数, 怎样从这些经验函数中选择最优的一个。文献[1] 提出的, 用几种函数进行拟合, 计算历史数据点实测值和拟合值的误差平方和最小的为最优经验函数这一方法, 可能存在这样的问题: 误差平方和最小, 但误差波动较大, 即一些点误差很小, 一些点误差相对较大。针对这种情况, 本文提出一种新的确定经验函数的方法, 用几种不同的函数进行拟合, 从中选取最优的经验函数, 最优经验函数确定的条件如下:

1) 历史数据点误差为正和误差为负的个数之差小于适应性参数, 在本文试验中选用的适应性参数为3。

2) 计算误差的方差, 方差最小的为最优的拟合函数。

方差是各个数据与平均数之差的平方和的平均数。通俗点讲, 就是各点与平均数偏离的程度, 来衡量一批数据相较于平均数的波动的大小。方差的数学描述为

式 ( 1) 中把x－视为第一个历史点的误差, xi为后面各历史点的误差, 该公式即可理解为各点误差相较于第一个点的误差的波动大小, 波动越小, 说明各个点的误差相差越小, 分布较均匀, 那么该函数即为最优经验函数。经验函数确定的流程图如图1所示。

1. 2自动分段拟合的步骤

本文选取三种较为典型的回归方程类型

分段拟合步骤如下, 分段拟合流程图如图2所示。

1) 根据1. 1节提出的方法拟合所有的历史数据, 在上述的三种拟合函数中选择最优的拟合函数。

2) 由第1) 步选取的经验函数, 计算历史数据点拟合值和实测值的误差, 计算历史数据点误差绝对值的均值S。

3) 比较各历史点误差与误差均值S, 若连续3个点的误差绝对值大于均值S, 则从第一个误差大于均值的点处分段, 执行第4) 步, 若不存在拟合值和实测值误差大于均值的点, 或者这种点只有一个或两个, 则不进行分段, 执行第5) 步。

4) 从分段点处到历史数据最后一个点重新拟合, 选择最优的函数, 重复以上步骤。

5) 根据以上步骤得到的分段, 以及各段选出的最优经验函数对历史数据进行拟合。

2数值实验

2. 1自动分段拟合实验

为检验所论述的自动分段多项式曲线拟合方法, 采用两组数据, 厦门高崎机场1990 ~ 2012年的起降架次的数据［9］和1990 ~ 2004年中国人均GDP的数据来进行拟合。

提出的自动分段多项式曲线拟合方法对厦门高崎机场1990 ~ 2012年起降架次统计数据和1990 ~ 2004年中国人均GDP数据进行拟合, 用JAVA编程实现, JFree Chart绘制曲线。表1和图3是厦门高崎机场起降架次的数据与经过分段曲线拟合后的数据的比较, 表2和图4是中国人均GDP数据与拟合后的数据的比较。

对于厦门高崎机场起降架次的数据, 提出的该方法自动将数据分为3段, 即数据1990 ~ 1994, 1994 ~ 2003, 2003 ~ 2012, 这三段数据选取的最优经验函数都是指数函数形式。对于中国人均GDP的数据, 程序自动将数据分为三段, 即数据段1990 ~ 1995, 1995 ~ 1999, 1999 ~ 2004, 这三段数据选取的最优经验函数分别是多项式函数形式, 直线形式和多项式函数形式。由表1、表2和图3、图4可以看出, 提出的自动分段拟合方法拟合精度较高, 能够很好的拟合历史数据。

2. 2对比实验

本文选取1990 ~ 2004年中国人均GDP数据, 分别用本文提出的自动分段拟合方法和传统直线拟合、指数拟合、多项式拟合方法对该组数据进行实验, 对比实验各种方法的相对误差如表3所示。

由表3可以看出, 自动分段拟合方法相对误差比传统拟合方法相对误差小, 且相对误差波动比传统方法波动小。自动分段拟合方法不仅在拟合精度上比传统方法高, 而且分段方便, 不需用户根据主观经验和绘制散点图来分段, 方便用户使用。

3结论

提出的自动分段曲线拟合方法采取对历史数据进行拟合, 自动选取较优的经验函数, 自动进行分段, 使得在进行数据拟合时, 不需要人为的绘制出历史数据散点图来选取经验函数和分段。通过编程实现和数值实验验证, 该方法不仅便于在计算机上编程实现, 而且拟合效果较好。

摘要：针对传统的分段曲线拟合方法在选择拟合函数和确定分段区间时经验成分较多的不足, 提出一种自动分段多项式曲线拟合方法, 根据误差方差和误差均值, 自动确定经验函数和分段区间。通过实际数据的检验, 验证了该方法的拟合效果。

关键词：数据拟合,分段拟合,多项式曲线,最小二乘法

参考文献

[1] 吴宗敏.散乱数据拟合的模型、方法和理论.北京:科学出版社, 2007

[2] 杨维, 张晓明.利用最小二乘法进行回归分析及经验公式的确定.沈阳工业大学学报, 1991;13 (2) :1—6

[3] 侯超钧, 曾艳姗, 吴东庆, 等.全局连续的分段最小二乘曲线拟合方法.重庆师范大学学报 (自然科学版) , 2011;28 (6) :44—48

[4] Grama L, Rusu C.Phase approximation by divide-and-conquer piecewise linear fitting of gain.International Symposium on Signals, Circuits and Systems, 2009:1—4

[5] 蔡山, 张浩, 陈洪辉, 等.基于最小二乘法的分段三次曲线拟合方法研究.科学技术与工程, 2007;7 (3) :352—355

[6] 刘晓莉, 陈春梅.基于最小二乘原理的分段曲线拟合方法.伊犁教育学院学报, 2004;17 (3) :132—134

[7] Gluss B.An alternative method for continuous line segment curve-fitting.Information and Control, 1964;7 (2) :200—206

[8] 王成山, 黄碧斌, 李鹏, 等.燃料电池复杂非线性静态特性模型简化方法.电力系统自动化, 2011;35 (7) :64—69

【最小二乘法拟合直线】推荐阅读：