病态问题广泛存在于大地测量数据处理过程中，主要是由观测信息不充分、对参数物理信息和先验信息不了解及计算方法选择不适当等原因所致^[1-2]，其研究方法可以大致分为3类：

1) 间接法。主要围绕平差模型构成的法方程进行研究，通过对法方程附加一定的约束条件进行正则化处理，需要选择合适的正则化参数，如正则化方法。正则化方法原理是通过附加全部或部分参数(或其改正数)加权平方和极小化条件来克服不适定性，其实质是通过增加约束或补充(先验)信息的方法，降低不适定性^[3]。在大地测量中，参数附加的约束不同，或者说参数的权阵选择不同，就有不同物理意义的解，许多学者对正则化参数的选择作了大量研究^[4-6]。利用正则化方法处理病态问题的不足主要表现在缺少对参数物理信息或先验信息的掌握，附加信息无从获取，对于一些复杂的先验信息无能为力，确定正则化参数具有主观随意性。

2) 直接解算法。通过观测方程设计矩阵分解的方式进行解算，不需要选择正则化参数，如截断奇异值法^[7]。截断奇异值法主要是通过数学手段截断相对较小的奇异值，用损失未知数的无偏性为代价换取均方根误差的减小，但截断参数的选取比较困难，对于条件数较大的病态问题效果不好。另外，实际工程技术中的法矩阵往往具有多个小特征值对模型求解构成危害，仅考虑克服最小特征值对解算结果的影响是不够的^[8]。

3) 迭代法。代表性的解法有谱修正迭代算法、共轭梯度法、增广方程组法、快速搜索法及其他算法^[8-11]，这些求解方法主要是基于非线性规划理论，求解过程十分复杂，计算量较大，只能解一些低阶及病态程度不高的平差问题，将非线性规划算法引入到测量数据处理中，还需要解决初值选定、收敛速度等关键问题^[12-13]。

本文通过分析近似计算中病态问题与局部最优解的关系，给出局部最优解的快速迭代方法，有效解决了大规模系数矩阵高病态的平差问题。

1 基于优化理论的病态原因分析

设平差模型L= AX+e，其中，A为m×n(m≥n)维列满秩矩阵，L为n维观测向量，e为n维随机误差向量，X=(x₁, x₁, …, x_n)^T为n维参数向量。最小二乘平差准则可以看作是无约束二次规划问题，即

$ {\rm{min}}\mathit{F}\left( \mathit{\boldsymbol{X}} \right) = {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right) $

(1)

式中，P为权矩阵。在目标函数F(X)=(L-AX)^TP(L-AX)=X^T(A^TPA)X-2L^TPAX+L^TPL中，L^TPL为一个常量，式(1)可以转化为二次规划问题，即

$ {\rm{min}}\ \mathit{ f}\left( \mathit{\boldsymbol{X}} \right) = \frac{1}{2}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{NX}} + {\mathit{\boldsymbol{c}}^{\rm{T}}}\mathit{\boldsymbol{X}} $

(2)

由于A为列满秩矩阵，N为正定矩阵，所以f(X)为严格凸二次函数。当法方程系数矩阵N=A^TPA为非病态时，式(1)的最优解为${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$=(A^TPA)^-1A^TPL。

对于N病态情形，现分析一个简单的例子：

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}} = \left( {\begin{array}{*{20}{c}} 1&3\\ 1&{3.000\;001}\\ {1.000\;001}&3 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}} \end{array}} \right) + \mathit{\boldsymbol{e}} $

(3)

其中，权矩阵P为单位矩阵，观测值L=AX_real+e由真值X_real=(1, 2)^T与随机误差e= (-0.037 161 457 061 89, 0.034 226 173 736 33, 0.028 947 447 163 43)^T生成。由于法方程系数矩阵N=A^TPA的条件数为cond(N)=3.458 9×10¹³，严重病态，利用最小二乘平差准则，计算出最小二乘平差解为${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$=(12 213.727 669 115 2, -4 068.908 663 853 5)^T，严重偏离真值，目标函数值为F(${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$)=0.003 076 176 855 38。下面3个解分别由截断奇异值法、修正奇异值法和快速搜索法得到：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_1} = }\\ {{{\left( {0.700\;867\;219\;175\;50, 2.102\;301\;190\;281\;92} \right)}^{\rm{T}}}, }\\ {{\rm{目标函数值为}}\mathit{F}\left( {{{\mathit{\boldsymbol{\hat X}}}_1}} \right) = 0.003\;164\;822\;221\;28;} \end{array} $

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_2} = }\\ {{{\left( {0.702\;087\;888\;090\;72, \;2.102\;187\;757\;056\;31} \right)}^{\rm{T}}}, }\\ {{\rm{目标函数值为}}\mathit{F}\left( {{{\mathit{\boldsymbol{\hat X}}}_2}} \right) = 0.003\;164\;885\;153\;43;} \end{array} $

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_3} = }\\ {{{\left( {0.703\;133\;226\;723\;54, \;2.101\;845\;854\;696\;11} \right)}^{\rm{T}}}, }\\ {{\rm{目标函数值为}}\mathit{F}\left( {{{\mathit{\boldsymbol{\hat X}}}_3}} \right) = 0.003\;164\;822\;167\;99。} \end{array} $

如果在计算中对于目标函数F(X)保留3位小数，上述4种方法得到的目标函数值都是0.003，此时，${{{\mathit{\boldsymbol{\hat X}}}_1}}$、${{{\mathit{\boldsymbol{\hat X}}}_2}}$和${{{\mathit{\boldsymbol{\hat X}}}_3}}$都是二次规划式(1)的最优解。因为N病态，其特征根分别为0.000 000 000 000 87和30.000 008 000 001 13，如果保留3位小数，说明N有一个特征根为0，不是一个正定矩阵，因此F(X)不是一个严格的凸函数，其最优解不唯一。当N是满秩矩阵时，理论上，只要计算过程能够保持足够的精度，最小二乘解${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$一定是全局最优解，其目标函数值一定是最小的。但在现实中，没有绝对精确的数据，计算机的字长也限制了计算过程中不能保留足够长的小数位数，所以寻找一种满足某种约束的最优解，或寻找一个满足某种条件的局部最优解，是解决病态问题最根本的办法。

2 无约束二次规划问题共轭梯度算法

迭代算法的基本思想是：给定一个初始解X₀，按照某一迭代规则产生一个点列{X_k}，使得当{X_k}是有穷点列时，其最后一个点是式(2)的最优解；当{X_k}是无穷点列时，有极限点，其极限点是式(2)的最优解。设X_k为第k次迭代点，d_k为第k次搜索方向，α_k为第k次步长因子，则第k次迭代为：

$ {\mathit{\boldsymbol{X}}_{k + 1}} = {\mathit{\boldsymbol{X}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_k} $

(4)

不同的步长因子和不同的搜索方向构成不同的迭代方法，在最优化方法中，搜索方向d_k是f(X)在点X_k处的下降方向，满足f(X_k+α_kd_k) < f(X_k)。如果搜索方向d_k已经确定，从X_k出发，沿d_k可求得α_k，使目标函数沿d_k达到极小值，即

$ f({\mathit{\boldsymbol{X}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_k}) = \mathop {{\rm{min}}}\limits_{\alpha > 0} F({\mathit{\boldsymbol{X}}_k} + \alpha {\mathit{\boldsymbol{d}}_k}) $

(5)

则称这样的搜索为最优一维搜索，或精确线搜索，α_k为最优步长因子。由于函数f(X)是严格凸二次函数，在精确线搜索下各种算法是等价的，但对于一般的非线性函数，各种确定α_k的算法结果不一定相同。共轭梯度法是一种常用的搜索算法，设f(X)在X处的梯度为g(X)= NX + c，搜索方向d_k由如下公式确定：

${\mathit{\boldsymbol{d}}_k} = \left\{ \begin{array}{l} - g\left( {{\mathit{\boldsymbol{X}}_k}} \right), k = 0\\ - g({\mathit{\boldsymbol{X}}_k}) + {\beta _k}{\mathit{\boldsymbol{d}}_k}, {\rm{ }}k \ge 1 \end{array} \right. $

(6)

对于β_k的确定，Fletcher等^[14]提出β_k=$\frac{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}$，Klessig等^[15]提出β_k=$\frac{{{{\left[ {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right) - g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right]}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}$，Murty等^[16-18]将以上2种方法结合起来，实现了自动正交校正及自动重新开始的原则，使一般共轭梯度法在求解的精度和速度方面有所改善，具体如下：

设${{\bar \beta }_k}=\frac{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}$，v_k= $ \frac{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}$，β_k=β_k(1-v_k)，如果梯度向量g(X_k)和g(X_k+1)相互正交，那么v_k=0，则β_k=β_k，即与Klessig等^[15]结论相同；如果v_k=1，那么β_k=0，此时搜索方向自动回到最速下降方向，即算法重新开始。此算法称为正交校正共轭梯度法，即CGM-OC算法，具体步骤为：

1) 令k=0，初始点为X₀，F₀=F(X₀)，g(X₀)= NX₀+c, d₀=-g(X₀)；

2) 线性搜索，令$\mathop {{\rm{min}}}\limits_{\alpha > 0} F$(X_k+αd_k)=F(X_k+α_kd_k)，求得α_k；

3) X_k+1=X_k+α_kd_k；

4) g(X_k+1)=NX_k+1+c，新的搜索方向为d_k=-g(X_k+1)+β_kd_k，设${{\bar \beta }_k}=\frac{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}$，v_k=$\frac{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{{\left( {g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)} \right)}^{\rm{T}}}g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}$，如果|v_k| < ε(ε=0.175)，则β_k=β_k；否则，β_k= _k(1-v_k)；

5) 若满足收敛条件‖X_k+1-X_k‖ < δ₁及‖F(X_k+1)-F(X_k)‖ < δ₂(δ₁和δ₂为事先给定的精度)，则停止计算；否则，令k=k+1，转至步骤2)。

文献[16-18]证明，由该算法得到的序列{X_k}对应的函数f(X_k)是严格单调递减序列，如果变量X没有约束，直接利用该算法可以求出最优解。

3 迭代初始解的选取

前文已经说明，当N病态时，近似计算中的N可能不正定，式(2)的目标函数f(X)不是一个严格的凸函数，从近似计算的角度来看，其最优解不唯一。从不同的初始值出发，利用CGM-OC算法得到的解就是不同的局部最优解，下一步就是选择一个初始值来求局部最优解。

CGM-OC算法利用目标函数的负梯度方向作为迭代的搜索方向，每步迭代都是沿着负梯度方向取最优步长。在计算过程中，目标函数开始时下降较快，接近极值点时收敛速度显著变慢。病态情形下，最速下降方向仅是算法的一个局部性质，最速下降法在开始的几步迭代计算中步长较大，自变量的改变和目标函数值的下降幅度也较大；但当接近收敛点时，步长很小，自变量的改变量也很小，目标函数值下降很慢，导致收敛速度十分缓慢。从近似的角度来看，算法设置的迭代终止条件导致截断误差的出现，实际上得到一个局部的最优解，选择不同的初始解迭代可得到不同的局部最优解。选择一个岭估计作为初始值，即

$ {{\mathit{\boldsymbol{\hat X}}}_0} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}} + \alpha \mathit{\boldsymbol{I}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $

(7)

利用CGM-OC算法求解无约束二次规划问题，因为初始值${{\mathit{\boldsymbol{\hat X}}}_0}$是一个岭估计，A^TPA+αI能够改善法矩阵的病态性，从${{\mathit{\boldsymbol{\hat X}}}_0}$出发求得的局部最优解具有岭估计的特征，当α足够小时，${{\mathit{\boldsymbol{\hat X}}}_0}$能够非常接近最小二乘解${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$。显然，减小α可提高迭代速度，但当α的取值太小时，法方程系数矩阵的病态性无法改善，使迭代算法失去意义。本文关于α的确定方法来源于文献[9]，其选取原则是既能改善矩阵的病态性，又能提高迭代速度。

α的初值根据(A^TPA+αI)^-1的谱半径ρ= $ \frac{\alpha }{{\alpha + \lambda }}$确定，若要平衡矩阵病态性与迭代收敛速率，需对谱半径ρ取折衷值，这时α指数为λ指数的一半加1，即

$ \alpha = {10^{{\rm{abs}}((\lg({\lambda _{{\rm{min}}}})/2) + 1)}} \times {\lambda _{{\rm{min}}}} $

(8)

其中，λ为矩阵A^TPA的特征值，abs为求绝对值的函数。对于高病态矩阵，λ非常接近于0，为避免出现α=0的情形，设λ为微小值eps，即2^-52。理论上，A^TPA非负定，其特征值也应非负，但在计算中，当λ从负方向接近0时，由于计算误差的存在，特征值可能为负数，因此对λ取绝对值计算是必要的。

4 算例与分析

为得到基于快速搜索算法的病态平差模型解算方法，利用2个算例进行分析：1)算例1，利用Hilbert矩阵建立一个病态平差模型，分别设置n=100、500、1 000、2 000的阶数，用以分析CGM-OC算法的迭代速度和精度；2)算例2，构建一个GPS单点定位实例，比较一些常用方法处理病态问题的效果。

4.1 算例1

以Hilbert矩阵为系数矩阵构造平差模型，即L=HX+e，其中，H为n阶Hilbert矩阵，其元素满足${h_{ij}} = \frac{1}{{i + j - 1}} $(i, j=1, 2, …, n)，X为n维参数向量，设其真值X_real=(1, 1, …, 1)^T，L为n维观测向量，e为n维随机误差向量。实例中，观测值由L= HX_real+e生成，e是由MATLAB随机函数生成区间为[-0.005，0.005]的随机误差值。表 1为不同阶数对应Hilbert矩阵的条件数、CGM-OC算法迭代的次数、目标函数值$F\left( {\mathit{\boldsymbol{\hat X}}} \right) = {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{H\hat X}}} \right)^{\rm{T}}}{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{H\hat X}}} \right)}$及算法结果${\mathit{\boldsymbol{\hat X}}}$的偏差平方和$m\left( {\mathit{\boldsymbol{\hat X}}} \right) = {\left( {\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{X}}_{{\rm{real}}}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{X}}_{{\rm{real}}}}} \right)$。

由系数矩阵构成的法方程系数矩阵严重病态，表 1中给出了相应的条件数，但现有的一些病态问题算法无法处理特别大的条件数。CGM-OC算法的初始值由式(7)和式(8)给出，算法中设定的精度为δ₁=10^-7和δ₂=10^-8。由于不涉及逆矩阵运算，且初始解非常接近最小二乘解，迭代次数很小，算法的精度(表 1中的偏差平方和)说明算法的效果很好。

4.2 算例2

取历元间隔为2 s，观测5颗卫星，使用4个历元解算整周模糊度，误差方程的系数阵见表 2。设参数的真值为X_real=(12.256 4, 5.148 9, 8.503 1, -1, -13, 4, 16)^T，观测值L=AX_real+e，e是由MATLAB随机函数生成区间为[-0.05，0.05]的随机误差值。

由系数矩阵构成的法方程系数矩阵N=A^TA的条件数cond(N)=8.928 2×10⁹，严重病态，N的最小特征根r_min=0.000 000 004 468 52，由式(7)可算出α=4.468 5×10^-4，CGM-OC算法的初始解${{\mathit{\boldsymbol{\hat X}}}_0}$=(A^TPA+αI)^-1A^TPL，再利用CGM-OC算法进行解算，迭代算法中给定的精度为δ₁=10^-5和δ₂=10^-8。为比较不同算法的效果，表 3给出基于最小二乘方法(${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$)、正则化迭代算法^[8](${{\mathit{\boldsymbol{\hat X}}}_{{\rm{RBC}}}}$)、截断奇异值法^[7](${{\mathit{\boldsymbol{\hat X}}}_{{\rm{TSVD}}}}$)、奇异值修正法^[7](${{\mathit{\boldsymbol{\hat X}}}_{{\rm{R}}}}$)及CGM-OC算法的快速搜索法(${{\mathit{\boldsymbol{\hat X}}}_{{\rm{CGM-OC}}}}$)计算得到的结果。

从表 3可以看出：

1) 正则化迭代算法中没有给定初值，当使用其他初值时，计算效果非常差，完全不能解决病态问题。实例中采用了初值设定方法，才能接近其他的解。

2) 正则化迭代算法中，干扰源向量的个数s需要根据实际情况进行调节，在数据处理中使用不方便，实例中是参考了真值才给出的。

3) 奇异值修正法中，正则化参数通过L曲线法给定，其求取不方便。

4) 截断奇异值法虽然有一定的通用性，但截断参数的确定无法自动进行，需要按照不同的情况确定，不能适应大规模的病态问题。

5) CGM-OC算法是一种快速搜索算法，不需要设定参数，得到的解与其他算法一致，适用于大规模病态问题的解算。

5 结语

病态平差模型可以转化为无约束二次规划问题。当法方程系数矩阵N病态时，N可能非正定，相应的二次规划问题目标函数f(X)可能不是一个严格的凸函数，从近似计算的角度来看，其最优解不唯一。从不同的初始值出发，利用CGM-OC算法得到的解只能是一个局部最优解，该算法是利用目标函数的负梯度方向作为迭代的搜索方向，每步迭代都是沿着负梯度方向取最优步长。在计算过程中，目标函数开始下降较快，接近极值点时收敛速度显著变慢，病态情形下，最速下降方向仅是算法的一个局部性质，已有的病态问题迭代算法很少考虑初始的选择。在选择一个初始值求局部最优解时，CGM-OC算法没有采用任何矩阵求逆运算，不需要选择正则化参数和奇异值截断参数，在处理法方程病态情形中有较好的优势。本文借鉴最优二次规划算法，为病态问题的处理提供了一种新的方法，该算法简单明了，可以很好地应用于测量数据处理，有效提高参数估计的准确率。