非线性最小二乘平差.doc

o 6-2 非线性模型平...
o 6-3 非线性模型平...
o 6-4 非线性模型...
o 6-1 问题的提出
非线性最小二乘平差

6-1问题的提出
经典平差是基于线性模型的平差方法。然而在现实世界中，严格的线性模型并不多见。测量上大量的数学模型也是非线性模型。传统的线性模型平差中的很多理论在非线性模型平差中就不一定适用；线性模型平差中的很多结论在非线性模型平差中就不一定成立；线性模型平差中的很多优良统计性质在非线性模型平差中就不一定存在。例如，在线性模型平差中，当随机误差服从正态分布时，未知参数X的最小二乘估计具有一致无偏性和方差最小性。但在非线性模型平差中，即使随机误差严格服从正态分布，未知参数X的非线性最小二乘估计也是有偏的。其方差一般都不能达到最小值。
 对于测量中大量的非线性模型，在经典平差中总是进行线性近似（经典的测量平差中称之为线性化），即将其展开为台劳级数，并取至一次项，略去二次以上各项。如此线性近似，必然会引起模型误差。过去由于测量精度不高，线性近似所引起的模型误差往往小于观测误差，故可忽略不计。随着科学技术的不断发展，现在的观测精度已大大提高，致使因线性近似所产生的模型误差与观测误差相当，有些甚至还会大于观测误差。例如，GPS载波相位观测值的精度很高，往往小于因线性近似所产生的模型误差。因此，用近似的理论、模型、方法去处理具有很高精度的观测结果，从而导致精度的损失，这显然是不合理的。现代科学技术要求估计结果的精度尽可能高。这样，传统线性近似的方法就不一定能满足当今科学技术的要求。另外，有些非线性模型对参数的近似值十分敏感，若近似值精度较差，则线性化会产生较大的模型误差。由于线性近似后，没有顾及因线性近似所引起的模型误差，而用线性模型的精度评定理论去评定估计结果的精度，从而得到一些虚假的优良统计性质，人为地拔高了估计结果的精度。
 鉴于上述各种原因，对非线性模型平差进行深入的研究是很有必要的。非线性模型的平差和精度估计以及相应的误差理论研究也是当前国内外测绘界研究的前沿课题之一。

电子教材 > 第六章 非线性模型平差  > 6-2 非线性模型平差原理
一、非线性误差方程
测量中大量的观测方程是非线性方程。比如导线测量中，以待定点坐标为未知参数的角度观测方程和边长观测方程分别为：
（6-2-1）
式中：为待定点坐标的真值，分别为角度观测值和边长观测值的真误差。角度观测值和边长观测值的观测方程（6-2-1）式是待定点坐标真值（）的非线性函数。又如在GPS伪距测量中，第j颗卫星至测站k的几何距离的观测方程为：

也是测站点k的待定坐标真值（）的非线性函数。一般地，用L表示的观测向量，用表示的未知参数向量的真值，用△表示的真误差向量，则非线性观测方程可写为：
（6-2-2）
式中: ，是由n个的非线性函数组成的的向量;。（6-2-2）式就是我们所要讨论的一般的非线性模型。
在一般的非线性模型（6-2-2）式中，用未知参数向量和真误差向量的估计值代替其真值，得非线性误差方程如下：
（6-2-3）
式中：V为观测值的改正数向量（残差向量）；为参数向量的估值。 

 二、非线性模型平差
由非线性误差方程（6-2-3）式知，非线性误差方程（6-2-3）式中仅有n个方程，而有n + t 个未知数（n个观测值的改正数和t 个参数）。因此非线性误差方程（6-2-3）式是非线性不定方程组，有无穷组解。在这无穷组解中，必然有一组解能使
（6-2-4）
我们将满足（6-2-4）式的一组解作为最优解，并称（6-2-4）式所确定的为的一个非线性最小二乘[23]估计。本书中将求解非线性最小二乘估计的过程称为非线性模型平差。
可见，非线性模型平差与线性模型平差的是完全一致的。（6-2-4）式的几何意义就是观测空间至解空间的距离最短，或者说是解轨迹π上离观测值L最近的点（见图6-1）。L到π的距离就是‖V‖。

图6-1
在非线性模型（6-2-3）式中，若存在一阶连续偏导数，且的非线性最小二乘估计量存在，则残差向量V在处垂直于切空间T（见图6-1）[24]。
一、非线性最小二乘估计的近似解
当非线性模型（6-2-3）式的非线性强度[24]较弱时，可以将非线性模型在处线性近似，并用线性模型的求解理论和方法来近似地求解非线性模型（6-2-3）式。这也就是我们大家所熟悉的传统方法——线性化方法，即将非线性模型（6-2-3）式在处用台劳级数展开，取至一次项，得：
（6-3-1）
令 （6-3-2）
（6-3-3）
则（6-3-1）式可写为：
（6-3-4）
（6-3-4）式就是我们熟悉的间接平差的误差方程。
由间接平差知，根据最小二乘原理可解得：
（6-3-5）
于是参数X的非线性平差结果为：
（6-3-6）
例6-1（本例取自参考文献[24]）已知非线性模型为。其中参数和的真值为。的5个真值（用参数的真值X算得）和相应的5个同精度独立观测值列于表6-1。
表6-1 的真值和相应的观测值
i
1
2
3
4
5
真值
4.202834
3.258924
2.527006
1.959469
1.519394
观测值
4.20
3.25
2.52
1.95
1.51
观测方程为：

取参数X的近似值为。将观测方程在处线性近似，得误差方程：

由（6-3-5）式得：

于是，由（6-3-6）式得参数X的平差值为：

参数估值的真误差为：

其范数为：

二、非线性最小二乘平差的迭代解
当非线性模型的非线性强度很强时，线性近似可能产生大于观测误差的模型误差，所以对于非线性模型，一般采用迭代的方法求解。
求解非线性误差方程（6-2-3）式的最小二乘平差值，就是求参数X的估值，使
（6-3-7）
由于是一常量，所以（6-3-7）式等价于目标函数为
（6-3-8）
的非线性无约束最优化问题。
因为是的非线性函数，所以对（6-3-8）式求一阶偏导数，并令其为零，得不到的显表达式。故求不出的解析解。因此，我们只能设法寻找某一近似解，使
（6-3-9）
成立。寻找使（6-3-9）式成立的近似解，一般只有采用迭代的方法。为此，下面介绍几种常用的迭代方法。

1．牛顿法
设的极小值的一个近似值为，在附近将展为台劳级数，取至二次项得：
（6-3-10）
式中：
（6-3-11）
（6-3-12）
称为处的Hessian矩阵。
（6-3-13）
是在处的梯度方向。
由于是的一个已知的近似值，故（6-3-10）式只是的函数，为了求得使6-3-10）式成立的，将（6-3-10）式对求偏导，并令其为零，得：

移项后两边转置，顾及（6-3-12）式，得
（6-3-14）
当Gk非奇异时，由（6-3-14）式可解得使（6-3-10）式成立的：
（6-3-15）
当充分小时，能使（6-3-10）式成立。但由于未知，故不能充分小，需不断迭代，直至充分小，其迭代公式为：
（6-3-16）
（6-3-16）式就是牛顿迭代的基本公式迭代终止条件：
（6-3-17）
或 =0 （6-3-18）
由于是一个绝对值较大的数，而的各元素的绝对值都很小，因此，由于计算机有效数字的限制，以（6-3-17）式作为迭代收敛条件比（6-3-18）式作为迭代收敛条件收敛要快一些。
牛顿法的迭代步骤为：
（1）选取初值，并令k=0。
（2）按（6-3-11）式计算梯度方向，若=0则转至（7）。
（3）计算Hessian矩阵。
（4）解线性方程组（6-3-14）式，得。
（5）按（6-3-16）式计算新的近似值。
（6）计算目标函数值，若则转至（2）继续迭代。
（7）终止迭代，输出和，结束。
例6-2 在例6-1中，仍设，用牛顿法求例6-1中非线性模型的非线性最小二乘平差值。
解：由例6-1知P=I，故目标函数为：

将代入计算，G0后，按以上迭代程序迭代，结果列于表6-2。
表6-2 牛顿法迭代计算
k
1
2
3
4
5
6

-1.205024908
0.3991833382
0.02889398018
0.0001691624122
-3.9492×10-9
-2.4012×10-9

-17.1530503
7.037242713
0.4948407424
0.002938012264
2.4039×10-7
-1.5569×10-7

5.333013265
5.41719809
5.422708003
5.442744565
5.422744593
5.422744582

-0.2539145225
-0.2542573375
-02556634078
-0.2556720853
-0.2556720877
-0.2556720866

-40.21054702
-40.58524686
-40.63522342

-40.63549278
-40.63549281
-40.63549281
迭代6次后，有== -40.63549281，所以停止迭代，得X的非线性最小二乘解为

则，

由例6-1知，本迭代解与其真值的距离比线性近似解与其真值的距离要小一个数量级。
当初值取时，迭代发散，这说明牛顿法对初值很敏感。
2．信赖域法
牛顿法具有很快的收敛速度，但它总是局部收敛的。因为牛顿法的基本思想是用二次函数

去逼近。只有当充分小时，才能很好地逼近。
既然只有当充分小时，才能逼近，那么可以对dX加以限制，然后在限制条件下来寻求的极小值。这个思想相当于求解下列约束最优化问题：
目标函数：
（6-3-19）
约束条件：

式中：为一正数，它随迭代而变化。
约束条件 限制了，使的长度不大于，这样总在一个给定的小区域中活动。这个区域是可信赖的，所以称该方法为信赖域法。
常数取决于对的逼近程度。这个逼近程度可用下式来描述
（6-3-20）
越接近于1，对的逼近程度越好，于是
（6-3-21）
这样，可总结出信赖域法的迭代程序：
（1）选取初值，，。
（2）按（6-3-11）式和（6-3-12）式计算梯度方向和矩阵，若=0则转至（7）。
（3）按（6-3-15）式计算，并检查是否满足约束条件。若不满足，则采取适当方法对予以压缩。然后在区域内求使=min的。
（4）计算的新的近似值。
（5）按（6-3-20）式计算，并按（6-3-21）式确定。
（6）检查是成立。若不成立，则转（2）继续迭代。
（7）终止迭代，输出和，结束。
例6-3，设，=0.08，用信赖域法求解例6-1中非线性模型的非线性最小二乘平差值。
、g和Gk的表达式同例6-2，用信赖域法迭代计算的结果列于表6-3。
表6-3 信赖域法迭代计算
k
1
2
3
4
5
6

-1.1205024908
0.1247630234
0.00548626096
0.000123053314
7.605186×10-5
4.700803×10-5

-17.1530503
2.746427624
0.09609929786
0.002154167737
1.331179×10-3
8.228461×10-4

5.342193488
5.421482607
5.422716379
5.422727153
5.422733806
5.422737919

-0.2446971417
-2553841042

-0.2556656317
-0.2556680975
-0.255669209
-0.255670562

-40.21054702
-40.62555285
-40.63548246
-40.6254928
-40.63549279
-40.63549279
迭代六次，有== -40.63549279，所以停止迭代，得：
当初值时，和牛顿法一样发散。这说明信赖域法也与初值有关，仍然是局部收敛，并不像想象的那样全局收敛。
3．拟牛顿法
牛顿法是基于二次模型

的。当R的形式很复杂时，求R的二阶偏导数阵将非常困难。为了避免求二阶偏导数，我们考虑用一个仅包含一阶偏导数信息的对称矩阵去逼近，然后再按牛顿法予以迭代。可见拟牛顿法与牛顿法的差别就是用代替。
拟牛顿法的关键是寻找一个只包含一阶偏导数信息的矩阵。此处介绍按“数值法”确定方法[24]。
由（6-3-12）式并顾及（6-3-11）式有

根据多元函数偏导数的定义：

知：

令 （6-3-22）
则去掉极限后，得的近似矩阵：
（6-3-23）
用（6-3-23）式定义的对称矩阵既能较准确地逼近，又只包含R的一阶偏导数信息，不需要求二阶偏导数，有了后，一切迭代均按牛顿法进行。由于拟牛顿法一开始就要按（6-3-23）式计算阵，所以计算前除了给定X的初值外，还必须给定的初值。可以这样确定：当给定后，将减去一个很接近的向量。则差值就是。即
（6-3-24）
开始计算时，用和即可。
例6-4，用拟牛顿法求解例6-1中非线性模型的非线性最小二乘平差值。初值为。
由于给定，则取，由（6-3-24）式得。取定和后，按拟牛顿法迭代的结果列于表6-3-4。由表6-3-4可以看出，迭代6次后，有，停止迭代，得

与以上迭代解相同，这表明用数值法确定的阵，能很好地逼近阵。
表6-4 拟牛顿法迭代计算
k
1
2
3
4
5
6

-0.221245771
0.856623242
-0.03115185858
-0.00384076635
-0.00010347425
2.358884×10-7

-2.965959852
12.93742451
-0.4936312734
-0.0575278579
-0.0015617973
3.739449×10-6

5.451669534
5.420543661
5.422705929
5.422744651
5.422744561
5.422744560

-0.256609692
-0.255545476
-0.2556687702
-0.2556720931
-0.2556720852
-0.2556720849

-40.62536846
-40.48402317
-40.63523365
-40.63548925
-40.6354928
-40.6354928

2.5022628
2.753000625
2.765194851
2.763778092
2.76373428
2.76373428

26.20604801
33.01991882
31.12556489
31.13656682
31.1394768
31.1394768

398.1482354
533.4314266
472.3603199
469.4625468
469.7401401
469.7401401

2.311619915
3.003695448
2.780051726
2.763839956
2.763740298
2.763740270

21.2196711
35.54791966
31.443537
31.13982534
31.13799006
31.13798972

289.8392513
546.7992839
475.1482883
469.7863468
469.7540942
469.7540898
4、高斯——牛顿法
以上介绍的几种方法，都是求目标函数的非线性最优化算法。与我们在《误差理论与测量平差基础》中已掌握的平差方法相去甚远。而高斯-牛顿法则不同，几乎和我们已经掌握的平差方法相同。
高斯-牛顿法的基本出发点就是在初值处对非线性模型进行线性近似。并按传统的平差方法求出一次近似值，然后反复迭代，直至前后两次的值相等，即。迭代步骤如下：假设非线性模型（6-2-2）式存在一阶连续偏导数，且参数X之间相互独立，则在近似值处线性化，得误差方程：

式中： 为用按（6-3-2）式算得的误差方程系数矩阵。
根据最小二乘原理，有

求得后，再以为近似值继续迭代，其迭代公式为： （6-3-25）
终止迭代条件：。
高斯-牛顿法具有一定的合理性。因为若（6-2-2）式是线性模型，则有=B，=B。于是：

上式表明：若（6-2-2）式是线性模型，则由高斯-牛顿法从任意初值出发，经一次迭代就可得到最小二乘平差的精确解。当非线性模型（6-2-2）式的非线性强度
[1][1]较弱时，高斯——牛顿法是较好的方法。
例6-5，设，用高斯——牛顿法求解例6-1中非线性模型的非线性最小二乘平差值。
按（6-3-25）式迭代的结果列于表6-5
表6-5 高斯——牛顿法迭代计算
k
1
2
3
4
5

5.394141331
5.422298989
5.422744502
5.422744573
5.422744573

-0.250050
-0.255618
-0.255672
-0.255672086
-0.255672086

-39.78568664
-40.62829761
-40.63549238
-40.6354928
-40.6354928

当时，迭代发散。这说明虽然高斯——牛顿法有一定的合理性，但在具体执行时可能会产生一些问题。首先是对初值的依赖性较大。当初值较差时，会出现迭代发散现象，使迭代无法进行下去。好在我们在实际计算时，总是用观测值算出，即如此求得的初值与X的真值很接近，故一般可迭代收敛。
一、非线性最小二乘平差结果的统计性质
　　通过《误差理论与测量平差基础》的学习，我们知道在线性模型中，当服从正态分布时，最小二乘估计量和均为无偏估计。并且和均具有最小方差。即在线性模型中，当服从正态分布时，最小二乘估计量具有优良的统计性质。那么，在非线性模型中，当仍服从正态分布时，非线性最小二乘估计量和是否还有这些优良统计性质呢？回答是否定的（参见文献[24]）。非线性最小二乘估计量和为有偏估计，而且和的方差达不到最小值。

二、单位权中误差
　　文献[25]已推导出非线性模型平差中单位权方差的严密估计公式。由于非线性模型平差中单位权方差的严密估计公式非常复杂，建议在实际工作中仍用
（6-4-1）
去估计单位权方差，并称（6-4-1）式为非线性模型平差中单位权方差的近似估计公式。

三、非线性函数的误差传播
　　为简单起见，仅讨论独立观测的情况。设独立观测向量的真值为，观测值为，观测值的真误差为服从正态分布，即。其中

　　现有独立观测向量的非线性函数
（6-4-2）
式中为常数。现要求根据独立观测向量L的方差来求非线性函数的方差。这就是非线性函数的误差传播问题。为了根据观测向量L的方差求非线性函数的方差，将（6-4-2）式在观测值L处展为台劳级数，并取至二次项得：

式中：
（6-4-3）
令

对取数学期望，并顾及
（6-4-4）
得：

根据方差的定义知：

=
顾及（6-4-4）式，得：
（6-4-5）
例6-6 在GIS矢量数据库中得A、B两点在空间直角坐标系中的坐标观测值为
L的方差阵为：

试求空间直线

的方差
解：将观测值代入上式可算得该空间直线得近似值为：

于是有：

由公式（6-4-5）得：

若将y的展开式仅取至一次项，即按线性函数求y的方差，则有：

可见。这说明一个非线性函数若按线性函数计算其方差，则人为地将函数的精度过度拔高了。
　　以上就非线性函数展开至二次项的误差传播问题作了简要讨论。关于严格的非线性函数误差传播问题、参数估值及其函数的精度评定理论和方法，都是需要深入研究的课题。

参考文献：
　　[23] P. J. G. Teunissen. Nonlinear least squares. Manuscript Geodaetica, Vol.15,137-150,1990
　　[24]王新洲，非线性模型参数估计理论与应用，武汉大学出版社，2002年8月
　　[25]王新洲，非线性模型平差中单位权方差的估计。武汉测绘科技大学学报，2000（4），358-361