最小二乘多项式

最小二乘多项式（精选7篇）

最小二乘多项式 篇1

1 概述

多项式拟合在工程计算领域得到较广泛的应用, 如在强噪声地震资料中的应用, 短期潮位补缺中的应用, 机动目标运动补偿算法中的应用等等。在工程分析计算过程中, 多项式拟合的结果往往需要保存到文件中, 作为数据处理的中间数据。市场中存在很多可用于多项式拟合的软件 (如Grapher) , 虽说可以完成多项式的拟合, 并看到拟合曲线, 但大多数软件不能够一次性的查看多次拟合图形, 并形成比较, 且不能根据用户需要保存理想的拟合结果。因此很多企业或公司需要开发自己的软件已满足特定行业需求。那么, 编写程序来求解最小二乘多项式系数成为此类软件开发过程的关键技术, 青海师范大学陈桂秀老师在[1]中提供了程序法求解最小二乘多项式系数的方法。但其算法无论是空间复杂度还是时间复杂度, 都达到让人无法接受的结果。实验数据结果分析, 如果实验数据达到800条, 运行一次拟合就需要3-4分钟 (根据常用办公室计算机测试) 。这就为设计较小空间复杂度和时间复杂度的程序算法提供了契机。该文主要包括三个部分:1) 多项式拟合基本原理介绍;2) 程序算法;3) 通过实验验证结果的正确性。

2 多项式拟合的基本原理

假设给定数据点 (xi, yi) (i=0, 1, ..., m) , φ为所有次数不超过n (n≤m) 的多项式构成的函数类, 求, 使得

当拟合函数为多项式时, 称为多项式拟合, 满足式 (1) 的pn (x) 称为最小二乘拟合多项式。特别地, 当n=1时, 称为线性拟合或直线拟合。显然为a0, a1, ..., an的多元函数, 因此上述问题即为求I=I (a0, a1, ..., an) 的极值问题。由多元函数求极值的必要条件, 得

即

(3) 是关于a0, a1, Lan的线性方程组, 用矩阵表示为

式 (3) 或式 (4) 称为正规方程组。

可以证明, 方程组 (4) 的系数矩阵是一个对称正定矩阵, 故存在唯一解。从式 (4) 中解出ak (k=0, 1, …, n) , 从而可得多项式

所以, 求解最小二乘拟合多项式系数的步骤如下:

第一、计算正规方程组的系数矩阵和常数项元素

第二、利用解线性方程组的方法求出正规方程组的解

3 程序算法

求解系数矩阵的程序如下:

其中xs M为系数矩阵。时间复杂度为O (n*n) 。

通过比较可以看出, 此程序算法可以大幅度减少空间和时间复杂度。

使用C#语言实现的求增广矩阵的方法为:

4 实验

为了验证算法的正确性, 使用vs2005设计一个DEMO, 界面如图1所示。

原始数据有796条 (删除了严重偏差的4条) , 一次性完成从2次曲线拟合到10次曲线拟合, 并实现动态加载原始数据和保存拟合结果 (为了简化程序设计的难度, 坐标系并没有标明坐标值) 。程序运行后, 首先单击左上角“载入数据”按钮, 然后单击“绘图”按钮, 程序一次画出从2次拟合到10次拟合的原始数据点线图和拟合后的曲线图。其中黑色的为原始数据点连接成线的图形, 从图形可以看出, 数据点产生强烈震荡而不平滑。其中中间平滑的红线为拟合后的曲线图。很明显介于原始图形震荡的趋于中间的位置。说明了程序拟合结果的正确性, 实验中完成9次拟合共用时6秒钟 (包括绘制图形部分) , 完全在用户可以接收的时间范围之内。通过直观的图形, 用户可以轻松的得出结论——哪一次拟合的结果更有助于分析数据, 然后在相应的拟合图形上单击以保存拟合后的数据。

5 结论

该文借助矩阵压缩存储的思想, 通过数组法, 使用一维数组来保存系数矩阵各项值, 减少了时间复杂度和控件复杂度, 极大的提高了程序运行的速度, 为用户开发适合特定工程需要的软件产品提供思想上和方法上的指导。

参考文献

[1]陈桂秀.用程序求解最小二乘拟合多项式的系数[J].青海师范大学学报:自然科学版, 2010 (3) :15-17.

[2]钟伟, 杨宝俊, 张智.多项式拟合技术在强噪声地震资料中的应用研究[J].地球物理学进展, 2006, 21 (1) :184-189.

[3]陆伟, 刘杰, 孙艳军.多项式拟合在短期潮位补缺中的应用分析[J].科技信息, 2011 (4) .

[4]侯成宇, 沈一鹰.基于最小二乘拟合的机动目标运动补偿算法[J].现代雷达, 2009, 31 (1) :54-57.

最小二乘多项式 篇2

在实验数据的处理、分析时, 数据拟合是经常采用的一种方法。数据拟合的目标是找到能反映变量之间关系的一种表达式, 使其在某种准则下最佳地接近已知数据。其原理有最小二乘法、契比雪夫法等, 且以最小二乘法最为常见和重要。

随着人类认识能力的不断进步以及计算技术的快速发展, 对于变量之间的未知关系, 应用曲线拟合的方法揭示其内在规律具有重要的理论和现实意义。针对最小二乘曲线拟合的有关问题以及相应的MATLAB实现进行探讨。

1 最小二乘曲线拟合

曲线拟合是指:已知n个数据点 (xi, yi) , i=1, 2, ……, n, 其中xi不全相同, 寻求函数f (x;a1, a2, …, am) 的待定参数a1, a2, …, am的一组取值, 使得在这组取值之下, 函数f (x;a1, a2, …, am) 与已知n个数据点整体上最为接近。

最小二乘曲线拟合方法根据已知数据, 首先构造出能够反映含有待定参数的函数f (x;a1, a2, …, am) 与n个数据点 (xi, yi) , i=1, 2……, n偏离程度的函数:

然后应用数学方法求函数J (a1, a2, …, am) 的最小值a1, am2, i…n, am=J (a1, a2, …, am) , 此时a1, a2, …, am的取值就是所求的待定值。这样一组取值使得函数f (x;a1, a2, …, am) 与n个数据点在二次平方和意义下最为接近。

2 最小二乘曲线拟合的MATLAB实现

2.1 线性最小二乘曲线拟合的MATLAB实现

线性最小二乘曲线拟合的含有待定参数的函数形式为:

其中rk (x) 为事先选定的一组已知函数, ak为待定系数, k=1, 2, …, n, m

该方程组可简化为RTRA=RTY, 其中:

当r1 (x) , …, rm (x) 线性无关时, RTR可逆, 方程组RTRA=RTY有唯一解。

根据线性最小二乘拟合理论可得关于待求参数构成的矩阵A的解, 利用MATLAB下命令A=RY, 可以直接求得待求参数, 下面举例说明。

例1下表给出了一组实测数据

已知数据 (xi, yi) 的函数原型为:

试用已知数据进行曲线拟合, 求出待定参数的值。

在MATLAB中输入以下程序:

运行程序, 输出待定参数结果为:

下面画出拟合得到的曲线及已知的数据散点图:

由图1可见, 曲线拟合效果非常好。事实上, 所给数据就是由已知曲线:

2.2 非线性最小二乘曲线拟合

在最小二乘曲线拟合中, 如果待定参数不能全部以线性形式出现, 如指数拟合函数f (x) =a0+a1e-a2·x, 这便是非线性最小二乘曲线拟合。

非线性最小二乘曲线拟合与线性最小二乘曲线拟合的原理没有什么区别, 但是最小二乘的解常常难以通过人的手算实现, 从而制约了该方法的应用。随着计算机技术的进步、专业软件的不断涌现, 这一问题的求解已不再困难。但是, 非线性曲线拟合中初值的选取是一个重要的问题, 目前为止还没有固定的理论或方法给出一般性的结论。

在MATLAB中实现非线性最小二乘曲线拟合有三个常用的命令。

(1) lsqcurvefit () 命令, 其使用格式为x=lsqcurvefit (fun, x0, xdata, ydata) , 其中fun是要拟合的非线性函数, x0是初始参数, xdata, ydata是拟合点的数据, 该函数最终返回系数矩阵。

(2) nlinfit () 命令, 其应用格式为beta=nlinfit (x, y, fun, beta0) , 其中x和y是拟合点数据, fun是被回归 (拟合) 的函数, beat0是初始参数。

(3) lsqnonlin () 命令, 其应用格式为x=lsqnonlin (fun, x0) , 其中fun的定义与lsqnonlin () 函数中fun的定义有差别, x0仍为初始参数向量, 将输出的系数结果放在变量x中。

例2假设已知

并已知该函数原型为y (x) =a1exp (a2x) +0.54exp (a3x) .sin (a4x) , 其中ai为待定系数。

在MATLAB中输入以下命令:

%建立函数原型, 则可以根据他来进行下面的求取系数的计算

运行结果为:

所求得的系数与原式中的系数相近。

如果想进一步提高精度, 则需修改最优化的选项, 函数的调用格式也将随之改变。

在MATLAB中输入以下命令:

%修改精度, 即是修改其限制条件

%两个空矩阵表示系数向量的上下限a', res运行结果为:

下面画出由拟合得到的曲线及已知的数据散点图, 如图2。

3 结论

首先介绍了最小二乘法的原理, 其次对线性最小二乘曲线拟合和非线性最小二乘曲线拟合应用MATLAB进行了具体实现, 使得相关曲线拟合理论更加生动易懂, 这使我们可以在今后的研究和工作中建立曲线模型对相应对象进行曲线拟合, 并应用MATLAB来实现, 从而找到更好更形象的反映变量之间关系的曲线, 也就找到最合适的模型了。

参考文献

[1]薛定宇, 陈阳泉.高等应用数学问题的MATLAB求解[M].北京:清华大学出版社, 2004.

[2]刘琼荪, 何中市, 傅鹂等.数学实验[M].北京:高等教育出版社, 2000.

[3]梁国业, 廖健平.数学建模[M].北京:冶金工业出版社, 2004.

最小二乘多项式 篇3

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面。其中一种简单的情形是:根据相关影长数据确定参照物位置, 这一技术也称为太阳影子定位技术。在本文中, 根据2015 年全国大学生数学建模竞赛A题附件1 的数据[1], 确定出该地点的位置。

2 影长与时间的函数关系

2.1 影响影长的各个参数

太阳高度角是指太阳光的入射方向和地平面之间的夹角, 可以用太阳高度角来估计影子场电影。 太阳高度角h与实测地的地理纬度θ, 太阳直射角纬度 φ, 时角 Ω 满足以下函数关系式[2]:

其中, 关于各个位置参数确定如下:

太阳直射点每时都在向西移动, 每小时移过15 度经度, 在地理题的计算中可粗略取每天移动0.25 度纬度[3]。 据此求解出太阳直射角纬度与日期之间的关系, 进而通过MATLAB编程, 得到了每天太阳直射点的纬度, 部分表见表1, 在4 月18 日的太阳直射点纬度为 φ=7。

一个天体的时角表示该天体是否通过了当地的子午圈, 其数值则表示了该天体与当地子午圈的角距离, 以小时来计量 (1 时角=15 度) 。根据以上理论, 建立时角的关系表达式[4], 即:

其中Cr表示北京时;LC表示经度校正 (4min/度) , 经度大于120 的经度校正为正, 小于120 的经度校正为负;EQ为时差。

2.2 建立影长随时间变化函数

根据式 (1) , 我们得到了太阳高度角与时角、太阳直射点纬度以及实测地的地理纬度的函数关系式, 进而我们可以通过太阳高度角与影子长度之间的关系, 建立影长随时间变化的函数关系式, 其步骤如下:

根据太阳高度角的定义, 得到太阳高度角与影子长度之间满足以下关系:

其中L为影子长度;H为旗杆高度3m;h为太阳高度角;

根据式 (1) , 将各个参数值代入公式进行整理, 得到太阳高度角与各个参数之间的函数关系:

再根据式 (4) , 将 (5) 转化为:

其中:

3 利用最小二乘拟合确定直杆高度以及所处位置

3.1 模型的建立

已知对应时刻的影子长度 (ti, Li) , 影长与时间的函数关系, 见公式 (6) , 利用最小二乘原理, 建立如下模型:

其中c1, c2, c3的初始值设为1, 0.1, 10

3.2 参数及经度的确定

根据参数的初始值, 将21 组 (ti, Li) 代入拟合方程中进行拟合。 在拟合之前, 本文根据c1, c2, c3的表达式做了以下猜想:发现c1, c2关于其自身的自变量的表达式十分简单, 而c3表达式较为复杂, 并且c3在拟合方程中的表达式也较为复杂, 所以关于c3的拟合效果相对于c1, c2来说, 拟合效果较差, 可能得不到精确的解。

而在实际拟合过程也证明了我们的猜想。 在实际拟合过程中, 随着初始值的不断变化, c1, c2变化不大或者只在固定的值的范围内小幅度变化, 而c3的变化幅度较大, 且不能找到唯一确定的值, 因此通过求出的纬度以及直杆高度来确定经度。经求解:c1=2.1316m, c2=0.313或0.1287, 即H=2.1316m, φ=18.24N或7.3969N。

将H, φ 中的数据代入 (5) 中, 用MATLAB编程计算, 我们得到21个时间点所对应的经度位置, 我们发现, 这21 个经度变化都十分小, 故我们经度该位置的经度就选取求出的21 个经度的平均值, 得到4 个经度, 进而得到四个固定直杆可能所处的地理位置: (18.24N, 105.9528E) , (18.24N, 37.4972E) , (7.3969N, 107.1988E) , (7.3969N, 36.2512E) 。

4 地理位置的检验

根据测定的关于x, y的坐标, 我们发现随着时间的增加, 影子长度也在不断增加, 说明当地时间已经过了正午时间12:00, 而此时的北京时间为14:42-15:42, 结合北京时间以东八区120 度为标准, 每向西15 度, 当地时间往前1 小时的规律, 我们计算出在79.5 E时, 北京时间为14:52, 当地时间恰好为12:00, 因此我们所确定的经度必须大于79.5E。 进而确定可能的地理位置为: (18.24N, 105.9528E) 和 (7.3969N, 107.1988 E) 。

根据地图确定地理位置为 (18.24N, 105.9528E) 的地点为越南, 地理位置为 (7.3969N, 107.1988E) 的地点为泰国湾 (为靠近马来西亚的一个海峡) , 故这个位置不符合我们的标准, 故我们确定的直杆所处位置大致在越南。

5 模型的检验

我们将所确定的地理位置代入关于太阳高度角的式子 (5) 中, 画出14:42-15:42的影子长度的变化曲线l1, 同时拟合出附件1中影子长度随时间变化曲线l2, 两者之间的误差, 可以这21个时刻的误差分别为:0.0004, 0.0005, 0.0005, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0006, 0.0005, 0.0005, 0.0005, 0.0005, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004, 0.0004。两条影子长度随时间变化的曲线 (见图1) 。

根据图1, 结合21 个点的误差, 可以发现两条曲线误差最大的为0.0006, 误差非常小, 认为该模型合理。

参考文献

[1]http://www.shumo.com/2015cumcm.html[OL].

[2]http://baike.baidu.com/view/86609.htm[OL].

[3]http://baike.baidu.com/link?url=_a Uu Zt IGSo3DBZgj Wi SGC2Nq H9NNBj Ec XX37JYu5MWs F6M22qio Ndju Y81e Nje P31Llpv JKTBhf_e PWF9g5CXa#1[OL].

[4]王国安, 米鸿涛, 邓天宏, 李亚男, 李兰霞.太阳高度角和日出日落时刻太阳方位角一年变化范围的计算[J].气象与环境科学, 2007, 30:161-163.

广义逆的计算与最小二乘估计 篇4

在系统理论、信号处理和矩阵方程求解等许多应用领域中,往往需要给出形如:

的矩阵的某种类型的广义逆,其中表示样本数目的k通常都很大,而且随时间增长而增大,因而计算量也迅速增加,这在实时应用中是应该避免的。

本文给出了某些类型的广义逆的算法,并用所得结果,给出了求解由多个矩阵方程所构成的系统的最小二乘估计,从而推广了[1]的工作,并简化了[2]中有关求多个方程公共解的算法,文中使用广义逆的记号均引自[3]。

1广义逆的计算

设k≥1为一整数,为具有相同行数的矩阵且满足

再设的某种广义逆,我们将证明,存在着满足

的矩阵Mk和Gk,使得

为的同种广义逆。

1°,如果(·)-表{1}-逆;

2°,如果(·)-表{1,2}-逆;

1°的证明:注意到Qk-1为对称幂等矩阵且Ak=(AkTAk)(1)AkT=(AkTAk)(1)HkTQk-1[3],我们有

2°的证明:注意到对任意A(1,3)∈A{1,3},总存在着矩阵Z,使得A(1,3)=A-+(I-A-A)Z,因此AA-=AA+,于是由(1),(2)和(3)可知

为对称阵,另外由1°可知:

3°的证明:由2°,我们只需再进一步证明

引理3[5]如果Gk如下确定:

从引理1到引理3知,为获得的广义逆,必须求出

容易看出,以上两式的计算量将随着k的增大而迅速增大,为避免这一困难,注意到

或等价地,

显然,(5)～(8)的计算量与k无关.

综上所述,我们证明了以下两定理。

定理1设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

Q-1=I,Ak=Qk-1Hk,Qk=Qk-1-AkGk

则

定理2设Mk和Zk由(1),(2)和(3)所定义且对k=0,1,…,令

下面,我们将定理1与定理2应用于

其中和Hk为具有相同列数的矩阵,且注意到如下事实[6]

我们有

定理3设由(9)所定义,且令

则

定理4设和Zk由(9)和(10)所定义且对k=0,1,…,令

则

2最小二乘估计

考虑观测系统

其中Yk∈Rm×p,Hk∈Rm×n为已知,Vk∈Rm×p为零均值随机矩阵,X为待估未知矩阵,对任何k,称X赞k为X的最小二乘估计(LSE),如果X赞k为最优化问题

的最小范数解,而

及加权矩阵wi∈Rm×n皆正定,i=0,1,…,k.

首先考虑无加权情形,即的情形。

为使估计误差的方差最小[6,7],我们有

其中的Moore-Penrose逆.

为得到(14)的递推式,注意到定理4以及(9)、(10)和(13),我们有

定理5系统(11)的X的无加权的极小范数最小二乘估计可如下递推地给出

其中Gk由定理4递推而得。

证明:

如果m=p=1,则定理5便退化为[1]的结果。定理5还有如下简单的推论:

若(11)中的Vk=0,坌k,则定理5给出多个方程所构成的系统

的公共解,亦即我们有

则(16)的公共解可按

递推地求得。

该推论中要求计算Ak+,考虑广义逆的性质和定理3,如果Gk=Qk-1Ak(1),则(17)给出(16)的公共解,从而只需计算A(1)∈A{1}而使计算量大为减少。类似地,如果Gk=Ak(1,4),则(17)给出(16)的具有极小范数的公共解。

现在再考虑一般加权最小二乘估计,因

故可通过无加权问题(18)的右端来处理加权问题。其中

或等价地,

应用定理5于无加权问题(18),便有

定理6令P-1=0,Q-1=I,且对k=0,1,…令

H′kTCk,则系统(11)的加权极小范数二乘估计可如下递推地求得:

为使(11)的加权最小二乘估计是最优估计[8],其中加权矩阵应如何选择?答案如下:

设而其方差设为

则由[4]中的推论6.4.4知,当时,为(11)在给定Yo,…,Yk时的最优最小二乘估计。

关于定理3和定理4的更多应用,可进一步参考[6,9,10]等。

3结论

应用定理1、2、3的结果,将[1]中关于单变量情形下的最小二乘估计推广到了多变量的情形,给出了后者的最优最小二乘估计的算法,同时还简化了[2]中求多个相容矩阵方程公共解的算法。

摘要：给出了某些广义逆的计算方法,将这些结果应用于观测系统的参数估计,给出了系统无加权和加权最小二乘的最新估计结果,推广了[1]的工作,同时简化了[2]中求多个相容矩阵方程公共解的算法。

关键词：算法,最小二乘估计,广义逆,方差,加权

参考文献

[1]Albert A,Sittle R.W.A method for computing leastsquares estimators that keep up with data SIAM[J].Control,1965,3:384-417.

[2]Morris G.L,Odell P.L.Common solutions for a matrixequations with applications[J].Assoc comput.math,1968(15):272-274.

[3]A.B.Israel,T.N.E.Greville.Generalized inverses:Theory and Applications[M].New York:John Wiley,1974.

[4]Campbell S.L,Meyer,C.D.Jr.Generalized Inverses of lin-ear transformations[M].London:Pitman publishing Ltd,1979:104-115.

[5]Cline R.E.Representations for the generalized inverses ofa partitioned matrix[J].Soc.Indust.Appl.Math.1964(12):588-600.

[6]X.Liu(刘轩黄).Recursive computation of generalizedinverses with applications to optimal state estimation[J].Control theory and advanced technology.1995(10):1485-1497.

[7]Rao G.R.Linear statistical inference 2nd its applications(2nd edtion)[M].New York:John Wiley&Sons,1973.

[8]刘轩黄.无初始状态统计知识时多变量系统的最优平滑与预报[J].海南大学学报,1998,16(3):210-214.

[9]X.Liu(刘轩黄).Optimal and minimum energy optimaltracking problems[J].苏州科技学院学报,2005,22(1):1-9;22(2):9-18.

最小二乘多项式 篇5

空间直线拟合在工程中具有重要应用[1,2,3,4,5,6,7,8],如在计算机视觉[4]、三坐标测量问题[5]、机器人焊缝中心三维定位[6]和空间管线网络干线的确定[7]中。

由于最小二乘准则下问题求解比较容易,一些文献在该准则下给出确定空间直线的算法。文献[1]将直线理解为两个平面的交集,通过求取两个平面方程,从而确定直线。文献 [2,3]根据最佳平方逼近原理,证明直线必过数据中心,并用逐步迭代给定方向向量,从而确定点向式直线方程。文献 [5]则通过两次降维,先按全最小二乘准则拟合二维平面,将空间中的点投影到该平面上,再将该平面上数据利用全最小二乘准则拟合直线,从而确定空间直线。文献[6]先确定直线经过的点,再求导得到一个非线性方程组,解方程组确定方向向量,从而最终确定空间直线。这些算法比较复杂,本文给出一种基于主成分分析的简便算法,通过计算机仿真验证了算法的有效性。

由于最小一乘准则下确定空间直线是一个不可微优化问题,难度较大。文献 [7]提出了空间中按全加权最小一乘准则下确定直线问题,但没有给出具体的算法,文献[8,9]给出了平面上的求解方法,没有对空间中全加权最小一乘准则下确定直线问题进行讨论。本文给出一种基于无约束优化的搜索算法,并进行了计算机仿真。

1 空间直线拟合问题

空间直线拟合问题可以描述为:在三维空间给定n个点P1(x1,y1,z1),P2(x2,y2,z2),…,Pn(xn,yn,zn),求空间直线L:$\frac{x-x{}_0}{m}=\frac{y-y{}_0}{n}=\frac{z-z{}_0}{p}$,其中M(x,y,z)为直线上任意点,M0(x0,y0,z0)为直线上确定的一点,非零向量$\overrightarrow{s}=(m,n,p)$为直线的方向向量,使得目标函数:

${\displaystyle \sum _{i=1}^{n}d}{}^2({\rm P}{}_i,L)(1)$

为最小,其中$d({\rm P}{}_i,L)=\frac{\left|\overrightarrow{{\rm M}{}_0{\rm P}{}_i}\times \overrightarrow{s}\right|}{\left|\overrightarrow{s}\right|}$表示点Pi到直线L的距离[10] ,其中|·|表示向量的模。式(1)实际上是按全最小二乘拟合空间直线,等价于文献[1,2,3,4,5,6]中拟合空间直线的准则。

如在上述假设下,改为求目标函数:

${\displaystyle \sum _{i=1}^{n}d}({\rm P}{}_i,L)(2)$

最小,则是按全最小一乘拟合空间直线,这样得到的直线不易受奇异点(outliers)的影响。

如在上述假设下,改为求目标函数:

${\displaystyle \sum _{i=1}^{n}w}{}_{i}d({\rm P}{}_i,L)(3)$

最小,其中wi(i=1,2,…,n)为权重,则是按加权全最小一乘拟合空间直线,这样得到的直线可以用于空间中管道的主干线确定问题。当wi=1(i=1,2,…,n)时,式(3)变为式(2)。

2 拟合空间直线的算法

2.1 全最小二乘准则下的拟合算法

将空间中三维分布的点集拟合为一维的直线,其本质就是降维,因此可以按照主成分的思想,对三维分布的空间数据进行主成分分析,提取第一主成分向量作为直线方向向量。文献[3]证明在最小二乘准则下,空间直线必通过数据中心点,因此可以由点向式确定空间直线方程。按此思想给出求问题式(1)中直线的算法如下:

第1步 输入三维空间中点坐标的n个数据点矩阵Xn×3,求出均值$\overline{X}=(\overline{x},\overline{y},\overline{z})$,拟合直线必过数据中心$(\overline{x},\overline{y},\overline{z})$;

第2步 求协方差矩阵$S=\frac{1}{n}(X-\overrightarrow{1}{}_n\overline{X}){}^{{\rm T}}\cdot (X-\overrightarrow{1}{}_n\overline{X})$,其中1n为n行1列元素全为1的n维列向量;

第3步 求出S的最大特征值λ1和对应特征向量$\overrightarrow{v}{}_1$,以$\overrightarrow{s}=\overrightarrow{v}{}_1=(a,b,c)$作为直线方向向量;

第4步 由点向式确定直线$\frac{x-\overline{x}}{a}=\frac{y-\overline{y}}{b}=\frac{z-\overline{z}}{c}$,其中(x,y,z)为直线上任意一点。

由主成分意义,点到直线距离平方和为n(λ2+λ3),其中λ2与λ3为S的其余两个特征值,与由式(1)得到的结果相同。

2.2 加权全最小一乘准则下的拟合算法

由于问题式(2)是式(3)的特殊情形,因此讨论式(3)。在问题式(3)中目标函数中含有绝对值,是一个不可微函数,因而求解比较困难。

要在空间中确定一条直线,无论采用确定两个相交平面还是点向式方程都至少需要确定6个未知参数。

将问题式(3)中的目标函数看作是一个含有6个自变量的不可微无约束优化问题,可以采用文献[11]中不依赖于导数的Nelder-Mead的单纯形算法与Hooke-Jeeves的模式搜索法求解。MATLAB中有基于单纯形算法的fminsearch()命令可以直接调用。

本文采用点向式方程确定直线,约定6个自变量按下面方式排列$\overrightarrow{v}=(x{}_0,y{}_0,z{}_0,a,b,c)$,给出求问题式(3)中直线的算法如下:

第1步 输入三维空间中点坐标的n个数据点矩阵Xn×3,给出初始值,给出$\overrightarrow{v}=(x{}_0,y{}_0,z{}_0,a,b,c)$的初始值。初始值可以选为全最小二乘准则下给出的值。

第2步 利用初始值和MATLAB中fminsearch()对目标函数${\displaystyle \sum _{i=1}^{n}w}{}_{i}d({\rm P}{}_i,L)$进行迭代求解,给出收敛值,从而确定直线经过的点(x0,y0,z0)和方向向量(a,b,c)。收敛的最终结果(x0,y0,z0)不唯一,方向向量(a,b,c)经标准化并指定方向后唯一。

第3步 由点向式确定直线$\frac{x-x{}_0}{a}=\frac{y-y{}_0}{b}=\frac{z-z{}_0}{c}$,其中(x,y,z)为直线上任意一点。

2.3 空间点到拟合直线的投影

为了在计算机仿真中绘制空间中点到拟合直线的垂线,这里给出点到直线的垂足坐标确定方法。

由于空间拟合直线已经求出,故可以确定空间中每一个点Pi(xi,yi,zi)在直线上的垂直投影点P′i(x′i,y′i,z′i)。由空间解析几何知识[10]可知投影点P′i(x′i,y′i,z′i)为直线在经过点Pi(xi,yi,zi)且与直线垂直平面上的垂足,联立直线的参数式方程与直线垂直的平面方程。

$\{\begin{array}{l}x=\overline{x}+aty=\overline{y}+btz=\overline{z}+ct\\ a(x-x{}_i)+b(y-y{}_i)+c(z-z{}_i)=0\end{array}(4)$

可以解出$t=\frac{a(x{}_i-\overline{x})+b(y{}_i-\overline{y})+c(z{}_i-\overline{z})}{a{}^2+b{}^2+c{}^2}$,从而代入参数式方程解出(x,y,z),即为垂足P′i(x′i,y′i,z′i)。

3 计算机仿真

为了检验本文算法的正确性和有效性,给出下面例子。

仿真环境为 Intel酷睿2双核 P8400、2G内存、Windows Vista OS、MATLAB7.5环境。

例1 设直线方程为x=1+t,y=1+t,z=1+t,t为参数。

参数t取值区间为[-6,6],先将该区间等分为30个小区间,t取区间端点,从而可在直线上确定31个点,可以确定31×3的数据矩阵。对直线上每个点进行随机扰动,让x,y,z分别加上随机数种子为rand(‘twister’,10),rand(‘twister’,20),rand(‘twister’,30) 的服从[-1,1]上均匀分布的n维随机列向量。绘出这些点三维散点图,如图1所示。

表1给出不同准则和算法得到的数值结果。其中方向向量已经标准化并约定向上(c>0),加权最小一乘准则中取权重全为1,即最小一乘,优化使用MATLAB自带的fminsearch()。

按最小二乘准则,绘制出拟合直线和点到拟合直线的垂线,如图2所示。

按最小一乘准则绘制出拟合直线和点到拟合直线的垂线,如图3所示。

4 结 论

本文讨论了空间中的直线拟合问题,针对全最小二乘和全加权最小一乘准则下,分别给出了基于主成分分析和无约束优化的空间直线拟合方法,这些方法在计算机上实现起来比较简便,计算机仿真验证了算法的有效性。这对使用空间直线拟合方法的科研人员具有一定的参考价值。

空间中曲线拟合比较复杂,未见到统一的方法,然而特殊情况比较容易。如空间中拟合圆周,可以先拟合平面,将三维数据投影到拟合平面上(降维到2),再根据平面法向量的方向角对平面上投影数据旋转使得一个分量全部为0,将旋转后的数据拟合圆周,得到圆周对应的圆心和半径。最后进行旋转逆变换,将球面方程与平面方程联立即得空间圆周方程。

参考文献

[1]袭杨.空间直线拟合的一种方法[J].齐齐哈尔大学学报,2009,25(2):64-68.

[2]许海涛,高彩霞.直线拟合算法[J].电脑知识与技术,2009,5(4):864-867.

[3]杜明芳.空间直线拟合[J].北京印刷学院学报,1996,4(2):27-31.

[4]李畅,张剑清,邹松柏,等.自动识别街道立面凹凸边界[J].计算机工程与应用,2009(15):6-8.

[5]胡学军,王俊亮,王刚.全最小二乘法在三坐标测量中的应用[J].武汉工程大学学报,2008,30(4):112-113.

[6]郭学军,王喜平,刘贺平.机器人焊缝中心三维定位的解析方法[J].数学的实践与认识,2009,39(9):123-127.

[7]梁怡,吴可法.一类非光滑最优化问题的有限步解法[C].中国运筹学会第六届学术交流会论文集(章祥荪等主编),Global-Link出版公司,香港,2000,801-811.

[8]林诒勋,尚松蒲.平面上的点—线选址问题[J].运筹学学报,2002,6(3):61-68.

[9]冯守平,杨桂元.关于加权全最小一乘法[J].应用概率统计,2009,25(2):135-142.

[10]同济大学应用数学系.高等数学(第五版上册)[M].北京:高等教育出版社,2002.

最小二乘多项式 篇6

考虑下列二阶抛物问题

$\{\begin{array}{l}u{}_t(x,t)-\nabla \cdot (a(x)\nabla u(x,t))=f(x,t),\\ (x,t)\in \Omega \times J,\\ u(x,t)=0,(x,t)\in \Gamma \times J,\\ u(x,0)=u{}_0(x),x\in \Omega ‚\end{array}(1)$

(1)式中J=[0,T],Ω⊂Rn为有界凸域,边界Γ是Lipschitz连续的。假定a(x)有界,即存在常数α*和α*,使得

0≤α*≤a(x)≤α*,∀x∈Ω (2)

对上述问题人们提出了许多数值模拟格式,如有限元格式[1],混合有限元格式[2],扩展混合有限元格式[3],最小二乘混合有限元格式[4]等。现将最小二乘思想和扩展混合有限元方法相结合,提出并分析了最小二乘扩展混合有限元格式。该格式继承了最小二乘和扩展混合元方法的优点,即格式能同时高精度逼近未知函数,未知函数的梯度以及流体的通量,所形成的有限元方程组是对称正定的,有限元空间无须满足LBB相容性条件。数值算例说明了该格式的有效性。

1 最小二乘扩展混合元格式

定义空间Λ=(L2(Ω))n,V=H${}_0^1$(Ω),W=H(div;Ω)={q∈(L2(Ω))n,∇·q∈L2(Ω)},相应范数为$\parallel \mu \parallel {}_{0,\Omega }=(\mu ,\mu ){}^{\frac{1}{2}},\parallel v\parallel {}_{1,\Omega }=(\parallel v\parallel {}_{0,\Omega }^2+\parallel \nabla v\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}},\parallel q\parallel {}_{{\rm H}(\text{d}\text{i}\text{v};\Omega )}=(\parallel q\parallel {}_{0,\Omega }^2+\parallel \nabla \cdot q\parallel {}_{0,\Omega }^2){}^{\frac{1}{2}}$,其中(·,·)表示空间(L2(Ω))n中的内积。作时间剖分τ=T/N,N为正的常数,在时间tn=nτ,n=1,2,…,N处逼近方程的解。$\text{Ω}={\displaystyle \underset{{\rm K}\in \tau {}_h}{\cup }{\rm K}}$是Ω的一个剖分,h=max{diam(K):K∈Th}。设Vh,Λh和Wh分别是V,Λ和W的有限维子空间,具有通常逼近性质。Vh,Λh和Wh的一个标准的选法是k,r和s次的分片多项式空间。

记λ=-∇u,σ=aλ,并用一阶差分近似ut,则式(1)可变成离散形式I

$\{\begin{array}{l}\overline{\partial }{}_{t}u{}^n=\nabla \cdot \sigma {}^n=f{}^n+R{}_1^n,x\in \Omega ,\\ \lambda {}^n+\nabla u{}^n=0,x\in \Omega ,\\ a\lambda {}^n-\sigma {}^n=0,x\in \Omega ,\\ u{}^n(x)=0,x\in \Omega ,\\ u{}^0(x)=u{}_0(x)x\in \Omega \end{array}(3)$

(3)式中$\overline{\partial }{}_{t}u{}^n=(u{}^n-u{}^{n-1})/\tau ,R{}_1^n=\overline{\partial }{}_{t}u{}^n-u{}_t^n$。

对于(u,λ,σ)∈V×Λ×W和(v,u,q)∈V×Λ×W,定义双线性形式

$\begin{array}{l}a{}_{{\rm N}}(u,\lambda ,\sigma ;\upsilon ,\mu ,q)=(u+\tau \nabla \cdot \sigma ,\upsilon +\tau \nabla \cdot q)+\\ \tau (\lambda +\nabla u,\mu +\nabla \upsilon )+\tau (a{}^{-1}(a\lambda -\sigma ),a\mu -q)(4)\end{array}$

根据离散形式式(3),省略时间截断误差定义最小二乘扩展混合元格式

格式Ⅰ 给定初始值u${}_h^0$∈Vh,对于n=1,2,…,N,求(u,λ,σ)∈Vh×Λh×Wh,使得

$\begin{array}{l}a{}_{{\rm N}}(u{}_n^n,\lambda {}_n^n,\sigma {}_n^n;\upsilon {}_h,\mu {}_h,q{}_h)=(u{}_h^{n-1}+\tau f{}^n,\upsilon {}_h+\\ \tau \nabla \cdot q{}_h),\forall (\upsilon {}_h,\mu {}_h,q{}_h)\in V{}_h\times \Lambda {}_h\times W{}_h。(5)\end{array}$

2 解的存在唯一性

为证明格式Ⅰ解的唯一性,首先说明双线性形式aN(·;·)是强制的。

引理1 对任意的(v,μq)∈V×A×W存在c>0,使得aN(v,μ,q;v,μq)≥c(‖v‖${}_{1,n}^2$+‖μ‖${}_{0,n}^2$+‖q‖H(div n))。 (6)

证明 引入适当选取的常数β>0,式(4)可整理为

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)=(\tau \nabla \cdot q+(1-\beta )\upsilon ,\\ \tau \nabla \cdot q+(1-\beta )\upsilon )+\tau (q-a(\mu +\beta \upsilon ),\\ a{}^{-1}(q-a(\mu +\beta \upsilon )))+\\ 2\tau (\mu +(1-a\beta )\nabla \upsilon ,\mu +(1-a\beta )\nabla \upsilon )+\\ (2\beta -\beta {}^2)(\upsilon ,\upsilon )。(7)\end{array}$

由于υ(x)=0,x∈Γ,所以由Poincaré-Friedrichs不等式知,存在正常数CF,使得

‖v‖${}_{0,n}^2$≤CF‖∇v‖${}_{0,n}^2$ (8)

结合式(2)有

$\begin{array}{l}a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge \beta (2(\tau \alpha {}_{*}+C{}_F)-\\ \beta (\tau \alpha {}^{*2}+\tau \alpha {}^{*}+C{}_F))\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2。(9)\end{array}$

取β=τα*/(τα*2+τα*+CF),则

$a{}_{{\rm N}}(\upsilon ,\mu ,q;\upsilon ,\mu ,q)\ge C\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\ge C\parallel \upsilon \parallel {}_{1,\Omega }^2。(10)$

从而由式(10)和式(4),知

$\begin{array}{l}\parallel \mu \parallel {}_{0,\Omega }^2\le 2\parallel \mu +\nabla \upsilon \parallel {}_{0,\Omega }^2+2\parallel \nabla \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(11)\\ \parallel q\parallel {}_{0,\Omega }^2\le 2\parallel a\mu -q\parallel {}_{0,\Omega }^2+2\parallel a\mu \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(12)\\ \parallel \nabla \cdot q\parallel {}_{0,\Omega }^2\le \frac{2}{\tau }\parallel \upsilon +\tau \nabla \cdot q\parallel {}_{0,\Omega }^2+\frac{2}{\tau }\parallel \upsilon \parallel {}_{0,\Omega }^2\le Ca{}_{{\rm N}}\\ (\upsilon ,\mu ,q;\upsilon ,\mu ,q),(13)\end{array}$

结合式(10)—式(13)便证得式(6)成立。

定理1 设f∈L2(Ω),则格式Ⅰ存在唯一解(μ${}_h^n$,λ${}_h^n$,σ${}_h^n$)∈Vh×Λh×Wh。

证明 考虑空间V×Λ×W,其范数定义为‖υ‖1,Ω+‖μ‖0,Ω+‖Q‖H(div;Ω)。由引理1知双线民生形式αN(·;·)在空间Vh×Λh×Wh中是强制的;而且由Cauchy不等式易证αN(·;·)在Vh×Λh×Wh中是有界的。因此由Lax-Milgram引理知格式Ⅰ解存在唯一。

3 数值算列

为验证所提格式的有效性,在区域[0,1]×[0,1]上考虑方程

$\{\begin{array}{l}u{}_t-(au{}_x){}_x=f,(x,t)\in (0,1)\times (0,1],\\ u(0,t)=u(1,t)=0,t\in [0,1],\\ u(x,0)=\sin \pi x,x\in [0,1]\end{array}(14)$

(14)式中$a=\frac{1}{100},f=\left(1+\frac{\pi {}^2}{100}\right)\text{e}{}^t\sin \pi x$。

容易验证方程的真解为$u=\text{e}{}^t\sin \pi x,\lambda =-\pi \text{e}{}^t\cos \pi x,\sigma =-\frac{1}{100}\pi \text{e}{}^t\cos \pi x$。取Λh为分片常数多项式空间,Vh和Wh为分片线性多项式空间,用MATLAB编程计算近似解uh,λh和σh,绘出真解和最小二乘扩展混合元解的比较图,见图1—图3。

摘要：对二阶抛物问题提出了最小二乘扩展混合元数值模拟格式。数值算例验证了格式的有效性。

关键词：二阶抛物问题,最小二乘,扩展混合有限元方法

参考文献

[1]Thome V.Galerkin finite element methods for parabolic problems.Lecture Notes in Mathematics,Springer,1984

[2]Johnson C V.Thome,Error estimates for some mixed finite element methods for parabolic type problem.RAIRO Numer Anal,1981;15:41—78

[3] Guo Ling,Chen Huanzhen.An expanded characteristic mixed finite element method for a convection-dominated transport problem.Journal of Coputational Mthematics,2005;23:479—490

最小二乘多项式 篇7

近年来,加权总体最小二乘是测绘学科数据处理领域的热点研究问题。针对观测向量L和系数矩阵A同时存在误差的情况,在1980 年,Golub和Van Loan[1]从数值分析的观点首次对这种情况进行了全面的分析,并正式命名为总体最小二乘。通过对总体最小二乘问题的进一步研究,发现观测向量L和系数矩阵A并不可能都是独立等精度的,而在总体最小二乘算法的求解过程中是将两者当做独立等精度的情况进行求解的, 面对这一问题,Schaffrin和Wieser[2]提出了加权总体最小二乘方法。目前,比较新的加权总体最小二乘算法除了Schaffrin和Wieser[2]提出的Schaffrin-Wieser算法,还有就是2012 年Mahboub[3]提出的Mahboub算法,以及2012 年Amiri-Simkooei[4]提出的Amiri-Simkooei算法。本文主要研究了这三种算法,通过近似垂直、近似水平和普通直线拟合的模拟算例以及在六参数坐标转换模型中的应用对这三种算法进行了比较和分析。对今后在实际工程中的应用,具有一定的理论价值和现实意义。

1 加权总体最小二乘算法

1.1 Schaffrin-Wieser算法

误差方程式为:

式中,y ∈ Rn ×1为观测向量,A ∈ Rn × m为系数矩阵,∈ Rm ×1为未知参数,ey∈ Rn ×1、EA∈ Rn × m分别为观测向量和系数矩阵的误差向量,Qy∈ Rn × n、QA∈ Rmn × mn分别为观测向量和系数矩阵的协方差矩阵,其中QA被分解为Q0∈ Rm × m、Qx∈ Rn × n,Py∈ Rn × n、PA∈Rmn × mn分别为观测向量和系数矩阵的权阵,vec(·) 表示矩阵的拉直运算,表示kronecker积。

加权总体最小二乘的平差准则或目标函数为:

误差模型在该准则下的拉格朗日极值函数为:

通过推导,得Schaffrin-Wieser算法如下[2]:

(6) 当时,计算结束。

1. 2 Mahboub算法

2012 年Mahboub[3]提出了五条原则来构造系数矩阵的协方差矩阵[5]:

(1) 如果系数阵的某个元素重复出现,认为这两个元素100% 相关,因此这两个元素之间的协方差等于重复元素的自方差;

(2) 假如系数矩阵的某个元素以其相反数的形式重复,认定这两个元素100% 负相关,因此这两个元素之间的协方差等于重复元素的自方差的相反数;

(3) 如果系数阵的某个元素是常数,认为其方差为零;

(4) 系数阵中两个不同元素,若两者明显相关,他们的协方差即为其实际值,否则为零;

(5) 上述规则在同方差情况中同样适用,若元素是随机数只需用数字1 作为其方差,若是常数其自方差为零。

利用上述5 条规则,Mahboub在重新构造系数矩阵的协方差阵的基础上,提出了一种新的算法。具体解算步骤如下[3]:

(1) 根据上述五条规则生成QA;

(5) 当时,计算结束。

1. 3 Amiri-Simkooei算法

2012 年,Amiri-Simkooei根据经典最小二乘理论,并结合文献[3] 的迭代方法,提出了一种类似于最小二乘解形式的加权总体最小二乘算法。具体的解算步骤如下[4]:

(3) 当时,计算结束。

2 算例及分析

2. 1 直线拟合

为了比较各种加权总体最小二乘方法计算得到的待定系数估值的差别,本文利用MATLAB软件编程进行数值模拟实验,分别用上述三种算法,以及总体最小二乘方法( TLS)、最小二乘方法( LS)和混合总体最小二乘方法( LS-TLS) 进行求解系数,对直线的不同情况进行拟合,最后比较分析各种方法在不同情况下的差异。 本文将迭代的Frobenius范数差的阈值 δ0设为10- 8。

(1) 近似垂直的直线

模拟一条近似垂直于x轴的直线y=ax+b,其直线的设计模型为y=2000x-2000,其自变量x从0.1到10之间每隔0.1模拟一个观测点,因变量y的范围为(-1800,18000),共100组观测数据,然后在每组观测点的x和y的坐标上均加上均值0,方差为σ2i=1/i×diag(0.01,0.01)的随机误差,其中i∈(1,…,100)表示点号。作1000次随机试验,利用MATLAB软件计算,各种方法的参数估计结果见表1。图1为各方法求出的参数估计值与真值差值的范数分布图,图2为利用Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法和TLS算法对近似垂直直线拟合做100次随机试验,得到的参数估计值与真值的差值分布图。

从表1 的数据和图1 中可以看出,在近似垂直直线拟合中使用加权总体最小二乘方法( WTLS)进行求解,拟合的结果比TLS、LS以及LS-TLS方法拟合的结果都要好。而LS-TLS方法拟合的结果比TLS拟合的结果好。这是由于在LS-TLS方法中考虑到了在系数矩阵中存在常数项,不需要进行误差改正,从而使得得到的EA比较可靠。而在加权总体最小二乘方法中,不仅考虑到系数矩阵中某些元素是常数,不存在误差,而且考虑到各元素之间是不等精度的,以及各元素之间存在相互联系,所以采用WTLS解算的结果更加接近真值。从图2 中可以看出,在近似垂直直线拟合中Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法这3 种算法所求得的参数估值波动较小,而TLS算法求得的参数估值与真值的差值波动较大,这也直观反映出,在系数矩阵和观测向量中各元素精度不同时,采用加权总体最小二乘算法更好。

(2) 近似水平的直线

模拟一条近似水平的直线y = ax + b,其直线设计模型为y = x /500 + 1,自变量x从0. 1 到10 之间每隔0. 1 模拟一个观测点,因变量y的范围为(1. 0002,1. 02),共100 组观测数据,然后在每组观测点的x和y的坐标上均加上均值为0,方差为σi2= 1 / i × diag(0. 01,0. 01) 的随机误差,其中i ∈(1,…,100) 表示点号,作1000 次随机试验。表2为各种方法对近似水平直线模拟的估计值。图3 为各方法求出的参数估计值与真值差值的范数分布图,图4 为利用Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法和TLS算法对近似水平直线拟合做100 次随机试验,得到的参数估计值与真值的差值分布图。

从表2的数据和图3中可以看出,在近似水平直线拟合中使用加权总体最小二乘方法(WTLS)进行求解,拟合的结果比TLS、LS以及LS-TLS方法拟合的结果都要好。从图4中也能较直观地看出Schaffrin-Wieser算法、Mahboub算法、Amiri-Simkooei算法这3种算法所求得的参数估值波动较小,比较接近真值,且这3种算法所求得的参数估值的结果都较为接近。

(3) 普通直线

模拟普通直线y = ax + b,直线的设计斜率范围为( - 5 ≤ a ≤5,且a ≠0),b = 1 自变量x从0. 1 到10 之间每隔0. 1 模拟一个观测点,共100 组观测数据,然后在每组观测点的x和y的坐标上均加上均值为0,方差为 σi2= 1 / i × diag(0. 01,0. 01) 的随机误差,其中i ∈ (1,…,100) 表示点号。表3 代表不同斜率直线的截距的模拟结果,b的真值恒为1。表4 代表不同斜率直线的斜率参数模拟结果。

从表3 和表4 的数据可以看出,在普通直线拟合中使用加权总体最小二乘方法( WTLS) 进行求解,拟合的结果更为接近实际真值。而三种加权算法得到的参数估值也基本相同。从以上直线拟合的算例中可以看出,上述三种加权算法虽然表现形式不同,但实际上是等价的,它们都是通过构造Lagrange目标函数,在这一基础下,通过一定目的的公式变形得到参数估值的迭代求解公式。

2. 2 坐标转换算例及分析

六参数坐标转换模型可写成:

当控制点个数为n(n ≥ 3) 时,六参数转换模型可以表示为:

其中(Xs,Ys)、(Xt,Yt) 分别为原始坐标系和目标坐标系中控制点的坐标值;sX、sY为沿X轴和Y轴的缩放比例,β 为旋转参数,c1、c2 为平移参数,δ 为关于坐标系两个坐标轴旋转的尺度因子。

通过MATLAB软件进行数值模拟实验。选取13 个点的二维坐标,设转换参数a1 = 0. 9,b1 = -0. 8,c1 = 1,a2 = 0. 6,b2 = 0. 7,c2 = 5,按六参数转换模型得到新坐标系下的坐标值(x,y),实验数据如表5 所示。在两套坐标上均加上均值为0,标准差为 σi= i × diag(0. 001,0. 001) 的随机误差,其中i表示点号,以使各坐标点具有不同的精度。

分别用Mahboub算法、Amiri-Simkooei算法、TLS、LS、LS-TLS方法求解坐标转换参数, 由于Schaffrin-Wieser算法在系数矩阵协方差结构上的限制,无法应用于二维坐标转换模型,因此在本算例中未使用该方法。最后将各方法求解的参数估值列于表6。

从表6 的数据和图5 可以看出使用加权总体最小二乘方法(WTLS) 进行求解,拟合的结果更为接近实际真值。而WTLS方法得到的结果比LS-TLS得到的结果更好,则可以看出,由于在WTLS方法中引入了PA、Py权阵,达到部分修改系数矩阵A的作用,即只修改数据元素,而固定常数元素,以及使用了Mahboub算法中构造系数矩阵协因数阵的五条原则[3],使得到的系数矩阵的残差阵EA更为准确,从而较LS-TLS方法得到更合理的模型。

3 总结

上述三种算法有相同点,同时也存在着明显的不同点以及适用范围:

(1) Schaffrin-Wieser算法的特点在于系数矩阵的定权方面,提出了QA构造法QA= Q0Qx,其中QA被表示成两个矩阵的直积形式[2],能够使迭代计算更为简便。但是该方法的缺点在于拆分后的矩阵Q0可能是奇异的,从而导致QA也是奇异的,以及当QA中行向量和列向量没有固定形式时,无法使用,例如在二维坐标转换模型应用中;

(2) Mahboub算法是在Schaffrin-Wieser算法的基础上提出的,在算法中利用五条原则来重新构造系数矩阵的协方差矩阵,这使得在变量误差模型中系数矩阵的协方差阵结构得到更好的保持。这一算法也使得实施过程中协方差阵不必限制于其他条件,如Schaffrin-Wieser算法中的直积形式[2]; 但是该方法仍存在不足之处,由于在构造系数矩阵时需要人工计算或通过编写程序代码自动获得系数矩阵中各个元素的协因数值而使得该方法操作过程复杂,计算所需时间长;

(3) Amiri-Simkooei算法是在Schaffrin-Wieser算法和Mahboub算法的基础上,依据最小二乘理论提出的,因此它更容易被理解,应用起来更为方便、灵活,在高斯- 马尔科夫模型中,参数估值是根据线性迭代算法求得的, 因此应用AmiriSimkooei算法进行参数求解时能够避开需要将公式线性化这一麻烦[4],这样使得计算过程更为简便。而且利用最小二乘理论可以直接得到参数估值的协方差矩阵,而不必重新推导; 但是该算法在表示方面为了与最小二乘方法形式一致,在算法中引入了较多代换量,例如这也使得算法显得冗长。

在实际应用中,需要根据实际情况判断选用哪一种加权总体最小二乘方法更为合适。

摘要：介绍了求解加权总体最小二乘问题的Schaffrin-Wieser算法和Mahboub算法以及Amiri-Simkooei算法,通过MATLAB软件编程进行数值模拟实验比较了三种方法在近似垂直、近似水平和普通直线拟合中的应用以及在二维坐标转换参数求解中的应用,并分析了三种算法的区别与联系。