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

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

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

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

$ \boldsymbol{L}=\boldsymbol{A} \boldsymbol{X}+\boldsymbol{e} $

(1)

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

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

$ \left\{\begin{array}{l} \boldsymbol{L}=\boldsymbol{A X}+\boldsymbol{e} \\ \text {s. t. }\|\boldsymbol{X}\|^{2} \leqslant c \end{array}\right. $

(2)

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

$ \begin{aligned} &\min (\boldsymbol{L}-\boldsymbol{A} \boldsymbol{X})^{\mathrm{T}}(\boldsymbol{L}-\boldsymbol{A} \boldsymbol{X}) \\ &\text {s. t. }\|\boldsymbol{X}\|^{2}=c \end{aligned} $

(3)

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

$ \boldsymbol{X}(\lambda)=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{A}^{\mathrm{T}} \boldsymbol{L} $

(4)

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

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

$ \boldsymbol{A}=\boldsymbol{P}\left[\begin{array}{l} \boldsymbol{B} \\ 0 \end{array}\right] \boldsymbol{Q}^{\mathrm{T}} $

(5)

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

$ \boldsymbol{B}=\boldsymbol{C S D}^{\mathrm{T}} $

(6)

那么，

$ \boldsymbol{A}=\boldsymbol{P}\left[\begin{array}{c} \boldsymbol{C S D}^{\mathrm{T}} \\ 0 \end{array}\right] \boldsymbol{Q}^{\mathrm{T}}=\boldsymbol{U} \boldsymbol{\varSigma} \boldsymbol{V}^{\mathrm{T}} $

其中，$\boldsymbol{U}=\boldsymbol{P}\left[\begin{array}{ll} \boldsymbol{C} & 0 \\ 0 & \boldsymbol{I} \end{array}\right], \boldsymbol{\varSigma}=\left[\begin{array}{l} \boldsymbol{S} \\ 0 \end{array}\right], \boldsymbol{V}=\boldsymbol{Q} \boldsymbol{D}$。令$\boldsymbol{X}_{1}(\lambda)=\boldsymbol{Q}^{\mathrm{T}} \boldsymbol{X}(\lambda), \hat{\boldsymbol{L}}=\boldsymbol{P}^{\mathrm{T}} \boldsymbol{L}=\left(\begin{array}{l} \boldsymbol{L}_{1} \\ \boldsymbol{L}_{2} \end{array}\right)$，其中，L₁为n维列向量，L₂为(m－n)维列向量，则式(4)可以写为：

$ \boldsymbol{X}_{1}(\lambda)=\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B}+\lambda \boldsymbol{I}\right)^{-1} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{L}_{1} $

(7)

式(7)等价于求解：

$ \min \left\|\left(\begin{array}{c} \boldsymbol{B} \\ \sqrt{\lambda} \boldsymbol{I} \end{array}\right) \boldsymbol{X}_{1}(\lambda)-\left(\begin{array}{c} \boldsymbol{L}_{1} \\ 0 \end{array}\right) \right\|^{2} $

(8)

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

$ \boldsymbol{W}^{\mathrm{T}}\left(\begin{array}{c|c} \boldsymbol{B} & \boldsymbol{L}_{1} \\ \sqrt{\lambda} \boldsymbol{I} & 0 \end{array}\right)=\left(\begin{array}{c|c} \boldsymbol{B}_{\lambda} & \boldsymbol{Z}_{1} \\ 0 & \boldsymbol{Z}_{2} \end{array}\right) $

(9)

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

$ \boldsymbol{B}_{\lambda} \boldsymbol{X}_{1}(\lambda)=\boldsymbol{Z}_{1} $

(10)

由X₁(λ)=Q^TX(λ)得‖X(λ)‖²=‖X₁(λ)‖²。

设

$ \begin{gathered} \omega(\lambda)=\|\boldsymbol{X}(\lambda)\|^{2}=\left\|\boldsymbol{X}_{1}(\lambda)\right\|^{2}=\\ \boldsymbol{L}_{1}^{\mathrm{T}} \boldsymbol{B}\left(\boldsymbol{B}^{\mathrm{T}} \boldsymbol{B}+\lambda \boldsymbol{I}\right)^{-2} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{L}_{1} \end{gathered} $

(11)

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

$ \omega(\lambda)=\boldsymbol{L}_{1}^{\mathrm{T}} \boldsymbol{C} \boldsymbol{S}\left(\boldsymbol{S}^{2}+\lambda \boldsymbol{I}\right)^{-2} \boldsymbol{S C}^{\mathrm{T}} \boldsymbol{L}_{1}=\sum\limits_{i=1}^{n}\left(\frac{\sigma_{i} \widetilde{L}_{i}}{\sigma_{i}^{2}+\lambda}\right)^{2} $

(12)

式中，σ_i为B的奇异值，$\widetilde{\boldsymbol{L}}=\boldsymbol{C}^{\mathrm{T}} \boldsymbol{L}_{1}$。

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

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

$ \lambda_{i+1}=\lambda_{i}-\frac{\omega\left(\lambda_{i}\right)}{\omega^{\prime}\left(\lambda_{i}\right)-\left(\omega^{\prime \prime}\left(\lambda_{i}\right) \omega\left(\lambda_{i}\right) / 2 \omega^{\prime}\left(\lambda_{i}\right)\right)} $

(13)

式中，$\omega^{\prime}(\lambda)=-2 \boldsymbol{X}_{1}^{\mathrm{T}}(\lambda)\left(\boldsymbol{B}_{\lambda}^{\mathrm{T}} \boldsymbol{B}_{\lambda}\right)^{-1} \boldsymbol{X}_{1}(\lambda), \omega^{\prime \prime}(\lambda)=6 \boldsymbol{X}_{1}^{\mathrm{T}}(\lambda)\left(\boldsymbol{B}_{\lambda}^{\mathrm{T}} \boldsymbol{B}_{\lambda}\right)^{-2} \boldsymbol{X}_{1}(\lambda)$。

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

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

2 参数的确定方法

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

$ \sigma_{1} \geqslant \cdots \geqslant \sigma_{n-1} \geqslant \sigma_{n} $

对于测量平差数据处理中的一些病态问题，其系数矩阵中往往会存在一个接近于0的特征值，也就是说σ_n≪σ_n－1，当$ \widetilde{L}_{n} \neq 0$时，有：

$ c_{\mathrm{LS}} \approx\left(\frac{\tilde{L}_{n}}{\sigma_{n}}\right)^{2} \gg\left(\frac{\tilde{L}_{i}}{\sigma_{i}}\right)^{2}, 1 \leqslant i \leqslant n-1 $

式中，$\tilde{L}_i$与式(12)中相同，c_LS是λ=0处的ω(λ)值，因此对于足够小的λ，存在：

$ \frac{\sigma_{i}^{2} \widetilde{L}_{i}^{2}}{\left(\sigma_{i}^{2}+\lambda\right)^{2}} \ll \frac{\sigma_{n}^{2} \widetilde{L}_{n}^{2}}{\left(\sigma_{n}^{2}+\lambda\right)^{2}}, 1 \leqslant i \leqslant n-1 $

即

$ c \approx \frac{\sigma_{n}^{2} \tilde{L}_{n}^{2}}{\left(\sigma_{n}^{2}+\lambda\right)^{2}} \approx \frac{\sigma_{n}^{4} c_{\mathrm{LS}}}{\left(\sigma_{n}^{2}+\lambda\right)^{2}} $

由此可得到λ的估计值：

$ \hat{\lambda}=\sigma_{n}^{2}\left(\sqrt{\frac{c_{\mathrm{LS}}}{c}}-1\right) \approx \lambda^{*} $

(14)

下面分析式(14)中$\hat{\lambda}$作为迭代初始值的可行性。

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

定义一个新函数：

$ \varphi(\sigma)=(\sigma+\hat{\lambda} / \sigma)^{2} $

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

$ \varphi^{\prime}\left(\sigma^{*}\right)=2\left(\sigma^{*}+\hat{\lambda} / \sigma^{*}\right)\left(1-\hat{\lambda} /\left(\sigma^{*}\right)^{2}\right)=0 $

即

$ \sigma^{*}=\sqrt{\hat{\lambda}}, \varphi\left(\sigma^{*}\right)=4 \hat{\lambda} $

当σ=σ_n时，φ(σ_n)=σ_n²·c_LS/c。如果$\varphi(\hat{\sigma})$=φ(σ_n)，且$\hat{\sigma} \neq \sigma_{n}$，有：

$ \hat{\sigma}=\frac{\hat{\lambda}}{\sigma_{n}}=\sigma_{n}\left(\sqrt{\frac{c_{\mathrm{LS}}}{c}}-1\right) $

令$\theta=\sqrt{\frac{c_{\mathrm{LS}}}{c}}-1$，则$\sigma^{*}=\sigma_{n} \sqrt{\theta}, \hat{\sigma}=\sigma_{n} \theta$。为便于分析，根据以上变量，绘制φ(σ)在θ>1和θ≤1情况下的大致图形(图 1、2)。

1) 假设φ(σ_i)≥φ(σ_n), i < n，其等同于：

$ \left(\sigma_{i}-\sigma_{n}\right)\left(\sigma_{i}-\sigma_{n} \theta\right) \geqslant 0 $

(15)

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

$ \omega(\hat{\lambda})=\sum\limits_{i=1}^{n} \tilde{L}_{i}^{2} / \varphi\left(\sigma_{i}\right) \leqslant\|\tilde{\boldsymbol{L}}\|^{2} / \varphi\left(\sigma_{n}\right) $

(16)

令$\eta=\|\tilde{\boldsymbol{L}}\|^{2} / \tilde{L}_{n}^{2} \geqslant 1$，可以将式(16)简化为：

$ \begin{gathered} \omega(\hat{\lambda}) \leqslant \frac{\|\tilde{\boldsymbol{L}}\|^{2}}{\sigma_{n}^{2}} \cdot \frac{c}{c_{\mathrm{LS}}}=\frac{\|\widetilde{\boldsymbol{L}}\|^{2}}{\widetilde{L}_{n}^{2}} \cdot \frac{\widetilde{L}_{n}^{2}}{\sigma_{n}^{2}} \cdot \frac{c}{c_{\mathrm{LS}}} \leqslant\\ \frac{\left\|\tilde{\boldsymbol{L}}\right\|^{2}}{{\widetilde{L}}_{n}^{2}} c=\eta c \end{gathered} $

即

$ \left\|\boldsymbol{X}_{1}(\hat{\lambda})\right\|^{2} / c \leqslant \eta $

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

$ \varphi\left(\sigma_{i}\right) \geqslant \varphi\left(\sigma^{*}\right)=4 \sigma_{n}^{2} \theta, i \leqslant n $

因此

$ \begin{gathered} \frac{\left\|\boldsymbol{X}_{1}(\hat{\lambda})\right\|^{2}}{c} \leqslant \frac{\|\widetilde{\boldsymbol{L}}\|^{2}}{4 \sigma_{n}^{2} \theta c} \leqslant \\ (\eta / 4) \sqrt{\frac{c_{\mathrm{LS}}}{c}}\left(1-\sqrt{\frac{c}{c_{\mathrm{LS}}}}\right) \leqslant(\eta / 2) \sqrt{c_{\mathrm{LS}} / c} \end{gathered} $

综上所述，在θ≤1(c_LS≤4c)的情况下，$\left\|\boldsymbol{X}_{1}(\hat{\lambda})\right\|^{2} / c \leqslant \eta$。在θ>1(c_LS>4c)的情况下，当σ_i∉Y(i < n)时，上式仍成立；当存在一些B的奇异值在区间Y内时，$\left\|\boldsymbol{X}_{1}(\hat{\lambda})\right\|^{2} / c \leqslant(\eta / 2) \cdot \sqrt{c_{\mathrm{LS}} / c}$。说明$\left\|\boldsymbol{X}_{1}(\hat{\lambda})\right\|^{2} / c$是有限的，即$\omega(\hat{\lambda})$不会远大于c。

3 数值实验 3.1 数值模拟实验

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

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

3.2 病态测边网实验

病态测边网见图 4，P₁₀、P₁₁为未知点，假设其坐标真值分别为(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)。P₁，P₂，…，P₉为已知点，其坐标及到2个未知点的观测距离见表 1。P₁₀与P₁₁之间的观测距离为d_{10, 11}=13.107 9 m，各观测量精度相等，中误差为0.001 m。

计算发现，该法方程系数矩阵的条件数为4.592 0×10³，说明法方程病态。由于误差的随机性，导致先验信息存在很多不确定性，所以在实际应用中，大多采用前期数据处理的结果以及经验信息来确定参数的约束范围。本文利用岭估计法与奇异值分解法计算的未知点坐标改正数结果来确定约束值，即$c=\frac{\hat{\boldsymbol{X}}_{L}^{\mathrm{T}} \hat{\boldsymbol{X}}_{L}+\hat{\boldsymbol{X}}_{\mathrm{SVD}}^{\mathrm{T}} \hat{\boldsymbol{X}}_{\mathrm{SVD}}}{2}$=0.018 2 m²($\hat{\boldsymbol{X}}_{L}$表示岭估计计算的解，$\hat{\boldsymbol{X}}_{\mathrm{SVD}}$表示奇异值分解法计算的解)。分别用最小二乘法、岭估计法、奇异值分解法、球约束^[6]以及本文算法求解该算例，结果见表 2(单位m)。

表 2中，$\|\widetilde{\boldsymbol{X}}-\boldsymbol{X}\|^{2}$为计算值与真值的误差平方和，其中，$\widetilde{\boldsymbol{X}}$为待测点坐标的真值，X为其计算值。从计算结果可以看出，由于测边网的法方程组病态严重，导致误差在传播的过程中被放大，使得最小二乘法的误差平方和最大，求得的未知点坐标严重偏离真值；岭估计法和奇异值分解法虽然能够在一定程度上改善系统的病态性，但是其结果与加入参数先验信息的球约束以及本文的新岭估计算法相比还有一些差距。这说明将参数的先验约束信息融入到平差模型中可以有效改善平差模型的病态性。

本文算法在求解2个算例的过程中，因利用岭参数的估计式$\hat{\lambda}=\sigma_{n}^{2}\left(\sqrt{\frac{c_{\mathrm{LS}}}{c}}-1\right)$作为迭代的初始值与利用0作为迭代的初始值的‖X(λ)‖²/c的值不同，所以二者迭代的次数也不同，具体见表 3。

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

4 结语

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