一种新的整体最小二乘迭代解法.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 的解算问题这种算法的提出是建立在下面的数学定理之上的：定理：如果 ，且存在正交矩阵 ，mnCR1[.]mUuR，使得 1[.]Vvv 1(.),,in(,)T ppCVdiag m即 的 SVD 分解可以表示为 ，又若 C 的秩为 r，有 对于)Ssv Triiuv，则有：1TkiiCuv1） 122()mn.krakDC2） 2() 1i ,min(,).pkiFFrnk ik在上面的数学定理的基础上，对模型 ，进行变换得到 ，如果对矩Axb01TxAb阵 进行 svd 分解，按上面的数学定理，在以式 为平差准则的情Ab 2[;];min况下，可以得到下面的结论： 1, ,1, 1ˆ[.]ˆ[]][Tt tt ttxvvAbAbu（1）式中 t 表示待求参数的个数由于这种算法是建立在上面的既定数学定理之上，理论严密，自提出以来一直是解决 TLS 问题有效的方法之一。

但这种方法有个缺陷，就是Golub 是直接对矩阵 做 SVD 分解，认为 A、b 均含有误差，由于 A 中可能并不是所[;]Ab有的元素均含有误差，如果 A 中仅是部分元素含有误差，对其不加区分的做上面最小范数的约束，是不严密的，可能会造成解算结果偏差较大 [9]这一点会在下面的算例部分讨论1.2 迭代算法TLS 的另外一种实用解法是由 Burkhard Schaffrin 提出的迭代解法 [7]，其基本原理如下，设 TLS 平差模型如式（2）所示(2)210,~(,)yAeNP其中 y 是 n 维观测向量，A 是 维系数矩阵， 是 维参数向量，e 是 n 维随机观nmm测误差向量若令 A 矩阵的改正矩阵为 ，且令 则整体最AE20()~(,)AmvcENI小二乘平差的平差准则可以写成：(3)inTAe为了求得目标函数的最小值，联立式（2）与式（3）构造拉格朗日函数(4)(,)2(()TTTAA nAeyeIe逐项求导可得：(5)1ˆ02ˆ()102ˆAnATAeIyeE令 ，则整理式（5）可以得到,TTNAPcy(6)ˆˆ1TNcvy其中， ，在式（6）的基础上迭代过程可以按下面的过程进行：ˆˆTvyc1.由 ，计算(1)N(1)2.由式 和 ，计算ˆˆTvyc()i()ˆiv3.由式 ，计算(1)()i iv(1)i4. 计算 是否成立，否则重复 2、3 步，是则退出迭代。

)(ˆiiBurkhard Schaffrin 在给出了上述迭代解法之后，给出了一种参数精度评定的方法，但是给出的公式是在对原有方程近似的基础上给定的，得到的结果并不是整体最小二乘平差下参数的实际精度2．TLS 的新迭代算法由于实际应用中上述 TLS 解算方法存在一定的缺陷，本文提出一种新的迭代算法，并在此基础上给出精度评定的策略间接平差的误差方程式为： 整体最小二乘是考虑到了系数矩阵的误差lxBV这时的误差方程变为： ，现在方程中要求的未知数是观测值的改正数)（V，系数矩阵的误差 ，以及未知参数的改正数 x，常规最小二乘是直接极值的问题，首先目标函数是 ，又因为 V 是 x 的线性函数，因此可以直接求导，得出一阶minPT导数为零的极值点整体最小二乘不同，此时不仅目标函数不同，而且在误差方程式中右边有多出 项，这是一个非线性方程，这给方程的解算带来了困难B下面给出在给出既定平差准则情况下整体最小二乘迭代解法，整体最小二乘的观测方程为：(7)dB)XVL（这时系数矩阵是有误差的，设观测值的平差值是 ，系数矩阵的平差值是 ，未知参ˆLˆB数平差值是 ，引入平差准则：ˆX(8),2211,iˆˆ()()minjtinni ijiji jLB要求的 、 、 是在满足式(7)的基础上，使得式(8)取得最小值的最优解。

把式(7)带入ˆLBˆX式(8)中有：(9)21 1ˆˆ(dL)()innt ntijiiijijij jB把式(9)设为 F，F 分别对待求参数求导： t1111t2 211t111t22111ˆˆˆ0ˆFˆˆˆ0ˆFˆˆjjjjttj tjt jjnttnjtBXdLBXdLBX（ ） （ ）（ ） （ ）（ ） （ ）（ ） （ ）（ ） （ t11112221 1111112221ˆ0Fˆ ˆ. 0ˆˆˆˆˆˆˆ.ˆ ntjt t tj j njnnj j jt t tj tj tnjnntj j jt dLBXdLBXdLBxBXdL  ）（ ） （ ） （ ）（ ） （ ） （ ）(10)上面给出了 个方程，其中 和 是未知数，ntˆ(1,;,)ijnjt  ˆ(,)jxt共有 个待求参数，方程有唯一解但是方程是非线性的，而且方程的形式非常复杂，t如果通过线性化求解，则会面临方程不收敛的问题本文采用一种新的迭代方法计算，通过参数变换构造迭代方程首先对上面的式子分类，分为对 求导的式子以及对未知参数 求导的式子，对于对 求导的式子，以其中ˆijBˆjxˆijB的第一式为例化简：(11)t11112211ˆˆˆ0ˆ.()jj ttBXdLBXLdX（ ） （ ）将对 求导的 个式子整理如下：ˆijnt(12)2112111222121112 2ˆˆˆˆ().()ˆˆˆˆ.()()()ttt ttt tXBXBLdXXB（ +） 2112ˆˆˆˆ.()()tntntnttntXBLdX上式写成矩阵形式，可以表述为：(13)ˆˆTTbN其中， ， 。

由（13）式可以看出，获ˆTbNEX12( )nLdLd取了观测值和未知参数的平差值就可以得到系数矩阵的平差值对未知参数 求导的式子，仍以第一式为例化简：ˆjx1112221 112211122212112ˆˆˆˆˆ. 0ˆ.()ˆˆˆ.ˆˆt t tj j njnnj j jttnnBXdLBXdLBXdLBBX （ ） （ ） （ ） 11211ˆˆ.()0ˆˆ.()ntnnii itiiiBXdLB（14）将对 求导的 个式子整理如下：ˆjxt 211211121211ˆˆˆˆ.()ˆˆˆˆ.()nnnnii itiiinnnnitit itiitiBXBXdLB（15）上式写成矩阵形式，可以表述为：(16)ˆTTBXL由（16）式可以看出，获取了观测值和系数矩阵的平差值，就可以得到未知参数的平差值由此，参数 和 可通过联立方程（13）和（16）来求解，即 ，ˆBXˆˆTTbNBXL但这是非线性方程，是一个迭代的过程，具体操作可以描述如下：1） 获取未知参数的初值 。

02） 根据观测值信息以及未知参数初值 ，取 ，由式（16） ，求取未知参数的平0(0)ˆB差值 1)ˆX3） 根据 ，求 ，ˆTbNE(1)(1)ˆTbNEX4） 根据求得的 、未知参数的平差值和观测值信息，由式（13） ，求取设计矩阵平差值(1)1)ˆB5） 重复 2-4 步，直到两次计算的参数值之差小于一定的域值退出迭代，输出结果上述迭代方程的结构简单直观，和经典最小二乘方法的参数求解过程非常相似，编程实现简单3.精度评定3.1 TLS 单位权方差估计公式TLS 单位权方差的估计公式一直是 TLS 理论研究难以解决的问题，在国内外还没有一个定论的公式经典最小二乘中，估计公式都是用 除以自由度在 TLS 中对组成系TVP数矩阵的观测量也进行了改正，如何在这种情况确定平差自由度是解决单位权方差估计公式的关键但无论是何种情况，估计出来的单位权方差都应满足统计上的最优性质在无偏性的准则下可以得到下面的单位权方差估计公式 [10]：(17)T20VPˆtr()Q在式（17）基础上，确定 就成了解决这一问题的途径，但是完全可以避免繁tr()V琐的推导，只要采用实测数据，求得不同情况下的 的具体数值，从中就可以归纳tr()V出 的表达式，从而确定式(17)的具体形式。

的计算方法可以参见本文 3.2tr()VQ部分给出的精度评定策略1）实验一 验证 大小V为了测试 的值，在本文采用二维平面直线模型，实验数据取自参考文献[10]，Qk=0.4,b=1.2；分三对、四对、五对分别模拟了两组观测值，以验证 和观测值个数的关VQ系实验结果如下表一 不同情形下 的迹V实验一 实验二三对观测值 1.014 1.035四对观测值 2.186 2.175五对观测值 3.025 3.003（2）实验二：验证 与所求参数无关VQ为了验证 ，是否与所求的参数值大小有关，变更参数的值k=4; b=12;模拟观测值计算 ，在四对观。