非线性不等式约束平差的一种新岭估计算法

李文娜 宋迎春 邓 伟 谢雪梅

1 中南大学地球科学与信息物理学院,长沙市麓山南路932号，410083 2 中南林业科技大学土木工程学院,长沙市韶山南路498号，410004

处理测量数据病态问题时常用岭估计[1]方法，但该方法仅从数学角度对平差模型的病态性进行改善，缺乏对参数实际情况的考虑，难以保证显著提高平差精度。解决病态问题的另一个思路是充分利用一些基于前期观测和附加(或额外)的先验信息。这些附加的先验信息常表现为参数的约束条件。

近年来国内学者对等式约束[2]、不等式约束[3]、区间约束[4]的平差理论与算法的研究取得了许多新成果。但这些约束基本上都是线性约束，用来描述先验信息存在明显不足，对于解决数据处理中的问题也有一定的局限性。杨婷等[5]指出，当系数矩阵具有近似复共线性(模型病态)时，参数的最小二乘估计的长度将在向量空间中的某些方向上偏大，此时一般的线性约束条件不再适合。 最近，学者们提出一种新的非线性约束形式的“二次约束”，如范数约束、球约束和椭球约束，其“二次”特性更容易与平差准则融合，但直接将它们线性化有时会使约束条件失去意义。针对病态问题的非线性不等式约束平差模型的解算，左廷英等[6]提出一种带有球约束的迭代算法，该算法只对参数进行迭代。在实际应用中，通常采用参数的最小二乘解作为迭代的初始值。但当法方程系数矩阵严重病态时，参数的最小二乘解与参数真值相距甚远，此时不仅会延长迭代的时间，也会降低计算结果的精度。肖兆兵等[7]给出参数带椭球约束平差的具体算法，该算法获取高精度解的前提是需要选择一个合适的拉格朗日乘子迭代初始值。对于不同的观测数据，应该有不同的初始值，由于缺乏对平差模型的深入了解，该算法通常选取0作为迭代初始值，导致可能会出现迭代不收敛的现象。综合以上问题，本文基于岭参数与约束条件有关这一特性，提出一种新的岭估计算法求解非线性不等式约束平差模型，不仅可以保证迭代收敛，而且还能提高计算效率。

1 非线性不等式约束平差模型及其算法

一般的平差模型可以表示为：

L=AX+e

(1)

式中，L为m×1的线性化后的常数向量，A为m×n的观测值系数矩阵，设m≥n，rank(A)=n，X为n×1的未知参数向量，e为随机误差向量，期望为0(e～N(0,σ2I),Σ=σ2I为协因数矩阵)。

如果参数带有不确定性，对X加以先验约束，得到：

式中，‖·‖ 表示二范数，c≥0且已知。根据最小二乘准则，将式(2)转化为非线性不等式约束最小二乘问题。当X没有约束时，直接利用最小二乘法可求得参数解为XLS。若‖XLS‖2≤c，则该最小二乘解就是平差问题(2)的解；当系数矩阵病态时，往往有‖XLS‖2>c，式(2)中的非线性不等式约束问题可转化为非线性等式约束问题：

min(L-AX)T(L-AX)

s.t.‖X‖2=c

(3)

式(3)也可通过凸优化理论得到，因为(L-AX)T(L-AX)是一个凸函数，其最优值将在边界上取得。使用拉格朗日乘数法求解式(3)：

X(λ)=(ATA+λI)-1ATL

(4)

式(4)为通常的岭估计，其中，λ的取值与观测值相关，即与参数的先验信息有关。

为研究λ的确定方法，采用Householder变换将矩阵A进行双对角化[8]：

式中，P为m阶的正交矩阵，Q为n阶的正交矩阵，B为n×n的上双对角矩阵，且B和A的奇异值相同。因此，如果B的奇异值分解为:

B=CSDT

(6)

那么，

X1(λ)=(BTB+λI)-1BTL1

(7)

式(7)等价于求解：

为减少计算量，利用Givens变换简化式(8)。设W为Givens矩阵的乘积，即为2n阶的正交矩阵，则有：

式中，Bλ为n阶的上双对角矩阵。因此，X1(λ)可由式(10)计算：

BλX1(λ)=Z1

(10)

由X1(λ)=QTX(λ)得‖X(λ)‖2=‖X1(λ)‖2。

设

ω(λ)=‖X(λ)‖2=‖X1(λ)‖2=

将式(6)代入式(11)：

式中，σi为B的奇异值，=CTL1。

显然， 对于λ>0，ω(λ)是单调递减函数。由于ω(0)=‖XLS‖2>c，因此存在唯一的λ*>0，使ω(λ*)=c。对于非线性方程ω(λ)=c，很难求得其精确解。通常情况下，求出的λ只需要使‖X1(λ)‖2/c与1足够接近即可。

根据Halley迭代公式[9]可以得到λ的迭代公式为：

具体迭代步骤为：1)利用先验信息确定约束值c；2)设定λ的初始值λ0和一个非常小的限值ε；3) 将λn代入式(9)求出Bλn，根据式(10)得到X1(λn)，从而求出ω(λn)，计算|ω(λn)-c|<ε是否成立，若成立取出X1(λn)，则X(λn)=QX1(λn)作为最终解，迭代终止，反之则进行步骤4)；4)利用X1(λn)求出ω′(λn)、ω″(λn)，然后再利用式(13)求出λn+1，回到步骤3)。

在上述算法中，迭代初始值λ0的确定非常关键，只有选取适当的初始值，才能保证其迭代的收敛性，同时提高计算效率。根据上文可知，该算法中岭参数λ与参数的先验信息有关，基于此特性，下面推导岭参数的估计式，将其作为迭代的初始值，并阐述其合理性。

2 参数的确定方法

假设矩阵B的奇异值从大到小排列为：

σ1≥…≥σn-1≥σn

对于测量平差数据处理中的一些病态问题，其系数矩阵中往往会存在一个接近于0的特征值，也就是说σn≪σn-1，当n≠0时，有：

即

由此可得到λ的估计值：

λ*

(14)

ω()的值不会远大于c，否则令作为迭代初始值相比于0作为迭代初始值的优势不明显。

定义一个新函数：

φ(σ)=(σ+/σ)2

由此函数的性质可知，当σ>0时，φ(σ)有唯一的最小值，假设该最小值在σ=σ*处取得，则：

φ′(σ*)=2(σ*+/σ*)(1-/(σ*)2)=0

即

图1 θ>1时φ(σ)的图像Fig.1 Image of φ(σ) when θ>1

图2 θ≤1时φ(σ)的图像Fig.2 Image of φ(σ) when θ≤1

1)假设φ(σi)≥φ(σn),i

(σi-σn)(σi-σnθ)≥0

(15)

如果θ≤1，则式(15)成立；如果θ>1，且σi>σnθ，说明B的奇异值分布范围不在区间Y内，即σi∉Y，则式(15)仍成立。在这2种情况下有：

令η=‖‖2/，可以将式(16)简化为：

即

‖X1()‖2/c≤η

2)假设θ>1，且B的奇异值有一些分布在区间Y内，则存在：

因此

综上所述，在θ≤1(cLS≤4c)的情况下，‖X1()‖2/c≤η。在θ>1(cLS>4c)的情况下，当σi∉Y(i

3 数值实验

3.1 数值模拟实验

为验证本文新岭估计算法在处理病态问题中的有效性，采用文献[10]给出的第一类Fredholm积分方程的具体函数形式模拟数据进行计算，并将所得结果与利用L-曲线法确定岭参数的岭估计方法计算的结果进行比较，见图3。

图3 本文算法与岭估计法的比较Fig.3 Comparison of our algorithm and ridge estimation algorithm

从图3可以看出，用本文方法计算出的函数曲线与真值的函数曲线吻合程度更高，本文方法计算的结果与真值的误差平方和为0.080 1，而利用岭估计法计算的结果与真值的误差平方和达到0.713 6。虽然两者都能够对系数矩阵的病态性进行改善，但本文算法考虑了参数的先验信息，所以精度更高。

3.2 病态测边网实验

病态测边网见图4，P10、P11为未知点，假设其坐标真值分别为(0 m，0 m，0 m)和(7 m， 10 m，-5 m)，其坐标近似值分别取(0.01 m，-0.01 m，0.02 m)和(7.01 m，9.99 m，-5.01 m)。P1，P2，…，P9为已知点，其坐标及到2个未知点的观测距离见表1。P10与P11之间的观测距离为d10,11=13.107 9 m，各观测量精度相等，中误差为0.001 m。

图4 空间测边网[11]Fig.4 Spatial geodesic netwo

表1 已知点的坐标和边长观测值[11]Tab.1 Coordinates of known points and side length observations

表2 各种方法对病态问题解算结果的比较Tab.2 Comparison of solution results among various methods for ill-posed problem

表3 初始值对迭代次数的影响Tab.3 Influence of initial values on the number of iterations

从表3可以看出，将岭参数的估计式作为迭代的初始值是合理的，不仅能保证迭代收敛，同时还能提高计算效率。

4 结 语

大地测量中存在大量关于参数的先验信息，利用这些先验信息对参数加以约束，可以有效改善平差模型的病态性。本文针对参数带有非线性不等式约束的最小二乘问题，建立相应的平差模型，给出求解该平差模型的一种迭代算法。该算法在解算效率和解算精度上有较明显的优势，但由于误差的随机性和参数真值的未知性，导致参数的约束范围通常需要采用前期数据处理的结果来确定，所以此时解算结果的质量会受到前期数据处理提供的先验信息的影响。通常在实际应用中遇到的问题比函数模型所展示的问题要复杂得多，如何充分利用先验信息获取可靠的参数信息，进而提高算法的准确性，是接下来研究工作的重点。