一种新的整体最小二乘迭代解法.doc

一种新的整体最小二乘迭代解法
Jian Kong1 Yibin Yao1,2 Han Wu1 Jianqing Cai2
(1 School of Geodesy and Geomatics,Wuhan University,129 Luoyu road,430079,Wuhan,China)
(2 Institute of Geodesy，University of Stuttgart, Geschwister-Scholl-Str. 24D,D-70174 Stuttgart,Germany)
E-mail:ybyao@

[摘要] 整体最小二乘(TLS)问题首先由Golub和 Van Loan提出并给出了第一个数值稳定解法，经过二十多年的研究，已经从数学的角度给出了整体最小二乘有解的充分条件和解的形式。近年来，不少学者提出过许多TLS的新解法，其中较为实用的方法有SVD方法、迭代法。本文首先介绍现有常用的TLS解法，指出了这些方法实际应用中存在的缺陷，并在此基础上提出一种新的TLS迭代算法；TLS平差方法的精度评定一直是困扰测量数据领域的难题，本文在新迭代算法的基础上，提出了精度评定的策略，并通过算例归纳确定了TLS的自由度；最后本文通过工程算例验证了新算法的可行性。
[关键词] 整体最小二乘 非线性方程求解 TLS自由度 迭代计算 奇异值分解 坐标转换

0.引言
一直以来，最小二乘是测量数据处理理论中最基本、最常用的方法，但随着测量仪器精度的不断地提高，测量数据处理也更趋向于追求处理过程中估计理论的严密性。上个世纪末，有关学者提出在二维直线拟合问题中，由于观测点坐标在x、y方向均含有观测误差，如果将y坐标作为观测量，那么这时平差模型中不仅观测向量含有误差，由x坐标组成的系数矩阵也是含有误差的，而经典最小二乘数据处理过程中无法顾及这项误差[1,2,3]。
这类问题可以归结为系数矩阵含有误差的高斯-马尔科夫问题，在数据处理过程中，经常会遇到平差模型中系数矩阵也是有误差的情况，传统最小二乘处理过程中忽略掉这项误差，这样做显然是不合理的，估计出来的结果，从统计上来看是有偏的，而不是最优的[2]。整体最小二乘的提出正是为了解决这个问题，但是整体最小二乘解法的复杂性却制约了其自身的推广应用。为此，近年来不少学者给出了很多TLS的解算方法[3-8]，但这些方法在实际应用中存在缺陷，本文提出的迭代解法不仅解决了TLS算法复杂度的问题，而且在理论上完善了TLS算法。
1. TLS的SVD解法和迭代算法
1.1 SVD解法
TLS的SVD解法[3]由Golub和Van Loan提出，首次解决了TLS的解算问题。这种算法的提出是建立在下面的数学定理之上的：
定理：如果，且存在正交矩阵，，使得。即的SVD分解可以表示为，又若C的秩为r，有。对于，则有：
1）
2）
在上面的数学定理的基础上，对模型，进行变换得到，如果对矩阵进行svd分解，按上面的数学定理，在以式为平差准则的情况下，可以得到下面的结论：
（1）
式中t表示待求参数的个数。由于这种算法是建立在上面的既定数学定理之上，理论严密，自提出以来一直是解决TLS问题有效的方法之一。但这种方法有个缺陷，就是Golub是直接对矩阵做SVD分解，认为A、b均含有误差，由于A中可能并不是所有的元素均含有误差，如果A中仅是部分元素含有误差，对其不加区分的做上面最小范数的约束，是不严密的，可能会造成解算结果偏差较大[9]。这一点会在下面的算例部分讨论。
1.2 迭代算法
TLS的另外一种实用解法是由Burkhard Schaffrin提出的迭代解法[7]，其基本原理如下，设TLS平差模型如式（2）所示
(2)
其中y是n维观测向量，A是维系数矩阵，是维参数向量，e是n维随机观测误差向量。若令A矩阵的改正矩阵为，且令。则整体最小二乘平差的平差准则可以写成：
(3)
为了求得目标函数的最小值，联立式（2）与式（3）构造拉格朗日函数
(4)
逐项求导可得：
(5)
令，则整理式（5）可以得到
(6)
其中，，在式（6）的基础上迭代过程可以按下面的过程进行：
1.由，计算
2.由式和，计算
3.由式，计算
4. 计算是否成立，否则重复2、3步，是则退出迭代。
Burkhard Schaffrin在给出了上述迭代解法之后，给出了一种参数精度评定的方法，但是给出的公式是在对原有方程近似的基础上给定的，得到的结果并不是整体最小二乘平差下参数的实际精度。
2．TLS的新迭代算法
由于实际应用中上述TLS解算方法存在一定的缺陷，本文提出一种新的迭代算法，并在此基础上给出精度评定的策略。
间接平差的误差方程式为：。整体最小二乘是考虑到了系数矩阵的误差。这时的误差方程变为：，现在方程中要求的未知数是观测值的改正数V，系数矩阵的误差，以及未知参数的改正数x，常规最小二乘是直接极值的问题，首先目标函数是，又因为V是x的线性函数，因此可以直接求导，得出一阶导数为零的极值点。整体最小二乘不同，此时不仅目标函数不同，而且在误差方程式中右边有多出项，这是一个非线性方程，这给方程的解算带来了困难。
下面给出在给出既定平差准则情况下整体最小二乘迭代解法，整体最小二乘的观测方程为：
(7)
这时系数矩阵是有误差的，设观测值的平差值是，系数矩阵的平差值是，未知参数平差值是，引入平差准则：
(8)
要求的、、是在满足式(7)的基础上，使得式(8)取得最小值的最优解。把式(7)带入式(8)中有：
(9)
把式(9)设为F，F分别对待求参数求导：
(10)
上面给出了个方程，其中和是未知数，共有个待求参数，方程有唯一解。但是方程是非线性的，而且方程的形式非常复杂，如果通过线性化求解，则会面临方程不收敛的问题。
本文采用一种新的迭代方法计算，通过参数变换构造迭代方程。首先对上面的式子分类，分为对求导的式子以及对未知参数求导的式子，对于对求导的式子，以其中的第一式为例化简：
(11)
将对求导的个式子整理如下：
(12)
上式写成矩阵形式，可以表述为：
(13)
其中，，。由（13）式可以看出，获取了观测值和未知参数的平差值就可以得到系数矩阵的平差值。
对未知参数求导的式子，仍以第一式为例化简：

（14）
将对求导的个式子整理如下：
（15）

上式写成矩阵形式，可以表述为：
(16)
由（16）式可以看出，获取了观测值和系数矩阵的平差值，就可以得到未知参数的平差值。
由此，参数和可通过联立方程（13）和（16）来求解，即，但这是非线性方程，是一个迭代的过程，具体操作可以描述如下：
1） 获取未知参数的初值。
2） 根据观测值信息以及未知参数初值，取，由式（16），求取未知参数的平差值。
3） 根据，求，
4） 根据求得的、未知参数的平差值和观测值信息，由式（13），求取设计矩阵平差值。
5） 重复2-4步，直到两次计算的参数值之差小于一定的域值退出迭代，输出结果。
上述迭代方程的结构简单直观，和经典最小二乘方法的参数求解过程非常相似，编程实现简单。

3.精度评定
3.1 TLS单位权方差估计公式
TLS单位权方差的估计公式一直是TLS理论研究难以解决的问题，在国内外还没有一个定论的公式。经典最小二乘中，估计公式都是用除以自由度。在TLS中对组成系数矩阵的观测量也进行了改正，如何在这种情况确定平差自由度是解决单位权方差估计公式的关键。但无论是何种情况，估计出来的单位权方差都应满足统计上的最优性质。在无偏性的准则下可以得到下面的单位权方差估计公式[10]：
(17)
在式（17）基础上，确定就成了解决这一问题的途径，但是完全可以避免繁琐的推导，只要采用实测数据，求得不同情况下的的具体数值，从中就可以归纳出的表达式，从而确定式(17)的具体形式。的计算方法可以参见本文3.2部分给出的精度评定策略。
（1）实验一 验证大小
为了测试的值，在本文采用二维平面直线模型，实验数据取自参考文献[10]，k=0.4,b=1.2；分三对、四对、五对分别模拟了两组观测值，以验证和观测值个数的关系。实验结果如下。
表一 不同情形下的迹

实验一
实验二
三对观测值
1.014
1.035
四对观测值
2.186
2.175
五对观测值
3.025
3.003
（2）实验二：验证与所求参数无关
为了验证，是否与所求的参数值大小有关，变更参数的值。k=4; b=12;模拟观测值计算，在四对观测值下，得到结果为=2.000。
（3）实验结论
1. 的迹与观测值的个数有关。
2. 与所求参数的大小无关。
3. ，亦即采用TLS和采用LS具有相同的自由度。
4. TLS下的单位权方差估计公式：
（18）

3.2 TLS平差值的精度估计
利用式(13)、式(16)迭代计算可以得到系数矩阵以及待估参数的平差值，将结果带回观测方程中可得观测向量的平差值：
(19)
至此即完成了参数估计的工作。应该给出估计参数的精度信息，式(13)、式(16)是一个隐函数方程组，直接线性化很困难，但它内在的确定了如下两个函数关系：
(20)
、对观测信息、的线性信息可以用、和对、的一次偏导数来表示。即：
(21)
给出了上面的隐函数，可以运用多元函数的隐函数求导法则来确定上面的四项偏导数。
但如果要直接获取上面的四个偏导数，根据隐函数求导的知识，需要对上面推导得到的个式子同时求导，这时计算量是很大的，为了简化运算，根据全微分不变定理，可以仍然把上面的式子按原有的分类分开求导，从公式(16)中可以得到：
(22)
从公式(13)中可以得到：
(23)
这样求导可以减少解算各项偏导数的复杂度，根据全微分的不变性，公式(22)与公式(23)的表示与公式(21)是等价的。
在得到公式(22)与公式(23)之后，由于这两个式子都是线性的，B，L是测量数据，已知它们的误差信息，从而可以根据误差传播定律确定、的误差信息。在得到、协因数阵之后，对公式(19)线性化可以得到观测向量平差值的精度信息：
(24)
4.实例分析
下面以二维坐标系转换为例，分析验证本文所提出来的新的整体最小二乘迭代解法的正确性。
本文选取了某地3个坐标公共点数据，它们同时具有1980西安坐标系以及CGCS2000坐标系的坐标成果。用整体最小二乘的方法选用相似变换模型，编程计算得到了相应的结果。
表二 新迭代法与svd法参数平差值（单位 m）

新迭代法
svd法
经典最小二乘
Mcos(a)
1.000000444
1.0000018113
1.000000448
Msin(a)
-0.000000503
-0.0000004261
-0.000000469
disX
119.65814
115.6107
119.5321
disY
-8.3051
-12.6153
-8.2999
表三 迭代次数统计
实验
退出阈值
迭代次数
一
0.0001
2
二
0.00001
8
三
0.000001
10
表四 TLS与LS单位权方差

TLS
LS
单位权方差
0.0083
0.0116
从表二中可以看出，三种不同方法中，迭代法，经典最小二乘的结果符合，而与传统svd方法的结果相差较大，按照本文上面的分析，这是由于在相似变换模型中系数矩阵只是部分含有误差，而svd方法仍以作为平差准则是不合理的，另一方面，在系数矩阵中有的坐标观测量会出现两次，这时如果不加区分的采用作为平差准则，相当于对坐标观测值进行了加权，这样的处理方式也是不完善的[6-7]。同时为了考察本文提出的迭代算法的效率，表三中统计了不同迭代退出阈值时迭代的次数，可见新迭代法的计算工作量并不大。
5结论
在考虑到系数矩阵误差的情况下，整体最小二乘较最小二乘法而言，理论上更严密，计算结果具有明显的统计优势，数据分析结果也证明了这点。
整体最小二乘解算后的精度评定还未得到广泛的研究，本文中提出的精度评定方法，首先通过对隐函数求导提取估计量对观测量的线性信息，然后通过误差传播定律估计误差。这从数据处理的角度完善整体最小二乘的理论。同时，本文通过推演的方法确定了整体最小二乘下单位权方差估计公式。
由于传统整体最小二乘的解法十分复杂，实际应用中存在一定的缺陷，限制了方法的推广应用。本文给出的整体最小二乘的迭代解法，计算简便，易于编程，理论严密，并同时给出了相应的精度评定方法，对整体最小二乘的推广是有意义的。

参考文献
[1]俞锦成.关于整体最小二乘问题的可解性[J].南京师大学报(自然科学版).1996, 19(1):13-16
[2]邱卫宁,陶本藻,姚宜斌等.测量数据处理理论与方法[M].武汉大学出版社.2008-6
[3]Sabine Van Huffel, Hongyuan Zha.The Total Least Squares Problem. Handbook of Statistics: 377-407.
[4] 张青富,保铮,焦李成.一种求解整体最小二乘问题的神经网路.通信学报,15(4)1994: 79-85
[5] Musheng Wei.Algebraic relations between the total least squares and least square problems with more than one solution. Numer．Math, 1992, 62:123-148
[6] O.Akyilmaz,Total Least Squares Solution of Coordinate Transformation[J]. Survey Review.2007,1:68-80.
[7]Burkhard Schaffrin.Total Least-Sqares(TLS) for geodetic straight-line and plane adjustment.Anno lxv Bollettino Di Geodesiae Scienze Affinin.200,3:141-166.
[8] 丁克良,欧吉坤,赵春梅.正交最小二乘曲线拟合法[J].测绘科学 2007,32(3):17-19
[9] 董校洪.整体最小二乘法在工程测量上的应用[硕士论文].上海:同济大学,2009,3
[10] 武汉大学测绘学院测量平差学科组,误差理论与测量平差基础[M].武汉大学出版社.2003