秩亏问题在大地测量数据处理中经常出现，由于法方程系数阵奇异，不能使用经典的平差算法进行解算，给计算带来难度，并导致未知参数的平差结果不唯一。目前，解决秩亏问题有许多方法，如附加条件法通过建立基准条件来消除法方程的秩亏^[3]，广义逆法以解算广义逆矩阵的方式获得秩亏问题的平差值^[2]；拟稳平差法根据最小范数准则使未知参数拟合于秩亏网中稳定点的初值，得到未知参数的解^[3]等。这些处理秩亏问题的方法大多是利用数学技巧来保证解的唯一性，虽然可以得到唯一解，但大多是加入了一些人为的平差准则。测绘工程中常存在一些关于未知参数的先验信息，它们能够对未知参数形成有效的约束，可以作为虚拟观测信息来弥补观测不充分导致的模型秩亏问题。

结合先验信息建立不等式约束，并同误差方程一起平差，可对平差结果起到约束作用^[4]，避免解的不唯一。附不等式约束的平差方法多年来深受国内外专家学者的关注，1976年Liew^[5]提出把不等式约束的最小二乘解转换成线性互补问题，用优化理论思想求解；1981年Schaffrin^[6]用Liew的方法研究大地控制网的优化，填补了不等式约束应用于大地测量的空白；1993年Remondi^[7]将高程的先验信息转换成不等式约束，改善了模糊度的初始化；2008年宋迎春等^[8]结合凸二次规划理论提出规划类算法，把优化理论应用到不等式约束平差算法中；2018年谢雪梅等^[9]将不等式约束看成是一个可行域，建立基于Wolf-Powell算法的非精确快速搜索算法，得到了不等式约束平差的新方法。然而这些方法主要研究的是如何对不等式约束平差模型进行解算，建立的模型都不是秩亏的。当模型秩亏时，上述算法都将失效。

在变形监测、控制网平差等实际工程应用问题中，秩亏问题常伴随着先验信息出现。充分利用这些信息，对未知参数进行有效的约束，可以补充观测信息，解决因为观测信息不足导致的秩亏问题。2008年王乐洋等^[10]利用不等式约束表达各种先验信息特点，解决了大地测量反演解非唯一性的难题。由于带有不等式约束的秩亏平差模型是一类新的平差模型，没有现成的有效算法，本文针对未知参数具有不等式约束的情形，提出利用Lemke算法解决这类平差问题，并给出参数求解的具体步骤。最后，通过算例验证该算法的可行性，拓展现有的误差理论和研究方法。

1 不等式约束秩亏平差模型

平差模型为：

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}{\rm{, }}权阵\mathit{\boldsymbol{P}} $

(1)

式中，L为m维观测向量，A为m×n维系数矩阵，e为m维随机误差向量，X =(x₀, x₁, …, x_n)^T为n维未知参数向量。最小二乘平差准则为：

$ \begin{array}{*{20}{c}} {{\rm{min}}}&{f\left( \mathit{\boldsymbol{X}} \right)} \end{array} = {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right) $

(2)

因在目标函数(L - AX)^T P (L - AX)= X^T(A^T PAX-2 L^T PAX + L^T PL中，L^T PL是一个常量，故可以把式(2)转化为一个二次规划问题：

$ \begin{array}{*{20}{c}} {{\rm{min}}}&{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}}} \end{array} $

(3)

式中，若A列满秩，则A^T PA是正定矩阵，f(X)是严格凸二次函数；若A秩亏，且附不等式约束GX ≤ h，X≥0，则得到具有不等式约束的秩亏平差模型：

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}{\rm{, }}权阵\mathit{\boldsymbol{P}} $

(4a)

${\rm{s}}{\rm{.t}}\;\;\;\mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{h}}, \mathit{\boldsymbol{X}} \ge 0$

(4b)

此处没有把式(4b)中的2个不等式合成1个是为了便于式(6)的转换。由上面的分析可知，式(4)能够解决如下的二次规划问题：

${\rm{min}}\;\frac{1}{2}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{NX}} + {\mathit{\boldsymbol{c}}^{\rm{T}}}\mathit{\boldsymbol{X}}$

(5a)

${\rm{s}}{\rm{.t}}\;\;\;\mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{h}}, \mathit{\boldsymbol{X}} \ge 0$

(5b)

式中，N = A^T PA为n×n维的对称半正定矩阵，c =-(L^T PA)^T为n维的列向量, G为s×n维的行满秩矩阵，h为s维的列向量。

从函数模型上来看，不等式约束秩亏平差模型属于非线性规划类的范畴。经上述推导被转化为二次规划形式(5)：目标函数是二次函数，约束条件是不等式。鉴于二次规划函数的特殊性质^[11]，可以采用一些更有效的方法求解。较为经典的是Wolfe算法和Lemke算法^[12]。这2种算法均是线性规划中单纯形算法改进形成的。其中，Lemke算法在设计上精巧，简明易懂，计算稳定、效率高，是解决二次规划问题的主流方法之一。基于Lemke算法的这些优势，本文提出将其应用于处理不等式约束秩亏平差问题。

2 线性互补问题与Lemke算法

使用Lemke算法时，需先将二次规划问题转化为线性互补问题，然后再结合该算法求解。由于A秩亏，其法方程系数阵N是一个半正定矩阵，f(X)为非严格凸二次函数，局部最小值也是全局最小值，Karush-Kuhn-Tucker点是问题(5)的最优解。根据Karush-Kuhn-Tucker条件：

$\mathit{\boldsymbol{NX}} + {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{Y}} + \mathit{\boldsymbol{c}} - \mathit{\boldsymbol{u}} = 0,\;\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{h}} - \mathit{\boldsymbol{GX}} \ge 0$

(6a)

${\mathit{\boldsymbol{v}}^{\rm{T}}}Y = 0,\;{\mathit{\boldsymbol{u}}^{\rm{T}}}\mathit{\boldsymbol{X}} = 0$

(6b)

$\mathit{\boldsymbol{u}} \ge 0, \;{\boldsymbol{v}}{\rm{}} \ge 0, \;\mathit{\boldsymbol{X}} \ge 0, \;\mathit{\boldsymbol{Y}} \ge 0$

(6c)

整理式(6)得：

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{u}}\\ \mathit{\boldsymbol{v}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{N}}&{{\mathit{\boldsymbol{G}}^{\rm{T}}}}\\ { - \mathit{\boldsymbol{G}}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{X}}\\ \mathit{\boldsymbol{Y}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{c}}\\ \mathit{\boldsymbol{h}} \end{array}} \right] \ge 0 $

(7)

$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}^{\rm{T}}}{\rm{}}}&{{\mathit{\boldsymbol{v}}^{\rm{T}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\boldsymbol{X}}\\ {\boldsymbol{Y}} \end{array}} \right] = 0 $

(8)

令$ \boldsymbol{Z}=\left[\begin{array}{l}\boldsymbol{X} \\ \boldsymbol{Y}\end{array}\right], \boldsymbol{M}=\left[\begin{array}{cc}\boldsymbol{N} & \boldsymbol{G}^{\mathrm{T}} \\ -\boldsymbol{G} & 0\end{array}\right], \boldsymbol{q}=\left[\begin{array}{l}\boldsymbol{c} \\ \boldsymbol{h}\end{array}\right]$，得：

$ \mathit{\boldsymbol{MZ}} + \mathit{\boldsymbol{q}} \ge 0\mathit{\boldsymbol{Z}} \ge 0,{\mathit{\boldsymbol{Z}}^{\rm{T}}}\left( {\mathit{\boldsymbol{MZ}} + \mathit{\boldsymbol{q}}} \right) = 0 $

(9)

其中，M ∈ R^(n+s)×(n+s)，q ∈ R^n+s，Z ∈ R^n+s。

问题(9)称为线性互补问题(linear complementary problem)。由式(6)到式(9)的推导过程可知，二次规划问题(5)可以转换为线性互补问题，其作用是把不等式约束平差模型归化到一个线性互补问题中求解。其等价形式为：

$ \mathit{\boldsymbol{w}} - \mathit{\boldsymbol{MZ}} = \mathit{\boldsymbol{q}},\mathit{\boldsymbol{w}} \ge 0,\mathit{\boldsymbol{Z}} \ge 0,{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{Z}} = 0 $

(10)

若q ≥0，则(w, Z)=(q, 0)是线性互补问题(10)的互补基本可行解；若q < 0，则引入单位列向量E和单位列向量Z₀，原式转化为：

$ \mathit{\boldsymbol{w}} - \mathit{\boldsymbol{MZ}} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{Z}}_0} = \mathit{\boldsymbol{q}},\mathit{\boldsymbol{w}} \ge 0,\mathit{\boldsymbol{Z}} \ge 0,{\mathit{\boldsymbol{Z}}_0} \ge 0, \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;且\;{\mathit{\boldsymbol{w}}^T}\mathit{\boldsymbol{Z}} = 0 $

(11)

显然，当式(11)满足下列条件：

$ \mathit{\boldsymbol{q}} + \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{Z}}_0} \ge 0, \;\mathit{\boldsymbol{Z}} = 0\mathit{\boldsymbol{w}} = \mathit{\boldsymbol{q}} + \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{Z}}_0},\\ \;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{Z}}_0} = {\mathit{\boldsymbol{Z}}_0} $

(12)

则(w, Z, Z₀)为式(11)的可行解。若Z₀=0，则(w, Z)是线性互补问题(10)的互补基本可行解，并可通过拟互补基本可行解(w_i, Z_i, Z₀)迭代得到^[13]，其中，w_i^T Z_i=0，Z₀>0，i=1, 2, …, n。

优化理论将大于0的分量且个数少于n+s的拟互补基本可行解称为退化的拟互补基本可行解，将所有分量全大于0的拟互补基本可行解称为无退化的拟互补基本可行解。无退化的拟互补基本可行解对应的线性互补问题是非退化的^[11]。

对于非退化的线性互补问题，Lemke算法可以有效地求解。其本质是对单纯形算法修改，求二次规划问题的Karush-Kuhn-Tucker点^[14]。首先，确定准互补基本可行解(w, Z, Z₀)，其中，Z₀>0，Z₀=max{- q_i|i=1, …, n+s}=-q_S，S表示行数。然后，从准互补基本可行解出发，依据负检验数最大原则，即Z₀=max{- q_i|i=1, …, n+s}时引进非基变量Z₀；依据最小比值原则^[15]换出基变量w_S，完成进基变量和离基变量选择过程。最后，利用主元消去法得到新相邻的拟互补基本可行解(w_i, Z_i, Z₀)。求解过程沿特定的搜索方向迭代，并遵循以下规则^[14]：

1) 保持可行性，按照最小比值规则确定离基变量；

2) 保持准互补性，若w_i(或Z_i)是离基变量，则Z_i(或w_i)是进基变量。

在此原则下进行迭代，直到得到Z₀=0的互补基本可行解(w, Z)。根据式(10)可知，求得线性互补问题的互补基本可行解(w, Z)就是二次规划问题的K-T点。因Z =(X, Y)，故可得附不等式约束秩亏平差问题的平差值${\mathit{\boldsymbol{\widehat X}}}$。

算法步骤^[11]如下：

1) 求解线性互补问题(10)，若q≥0，停止计算，(w, Z)=(q, 0)是它的一个互补基本可行解；若q < 0，则引入单位列向量E和单位列向量Z₀，原式转化为：

$ \begin{array}{l} \mathit{\boldsymbol{w}} - \mathit{\boldsymbol{MZ}} - \mathit{\boldsymbol{E}}{\mathit{\boldsymbol{Z}}_0} = \mathit{\boldsymbol{qw}} \ge 0, \\ \mathit{\boldsymbol{Z}} \ge 0, {\mathit{\boldsymbol{Z}}_0} \ge 0\;\;{\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{Z}} = 0 \end{array} $

(13)

其中，引入的人工变量Z₀是单位列向量。依据优化理论中单纯形算法的思想，将式(13)的信息体现在表 1中。

表 1全面反映了不等式约束秩亏平差模型的信息，是Lemke算法解决线性互补问题依据的准则。式(9)中的秩亏平差模型和不等式约束模型构成了向量M和向量q。设

$ - {\mathit{\boldsymbol{q}}_s} = {\rm{max}}\left\{ { - {\mathit{\boldsymbol{q}}_i}|i = 1, \ldots , n + s} \right\} $

(14)

以S行为主行，Z₀对应列为主列，进行主元消去，令y_S= z_S。

2) 设在现行表变量y_S下面的列d_S，若d_S≤0，则停止计算。否则，按最小比值规则确定相应的主行数r，使满足：

$ \frac{{{{\mathit{\boldsymbol{\bar q}}}_r}}}{{{\mathit{\boldsymbol{d}}_{rS}}}} = {\rm{min}}\left\{ {\frac{{{{\mathit{\boldsymbol{\bar q}}}_i}}}{{{\mathit{\boldsymbol{d}}_{iS}}}}{\rm{|}}{\mathit{\boldsymbol{d}}_{iS}} > 0} \right\} $

(15)

如果r行的基变量是Z₀，则转步骤4)；否则，进行步骤3)。

3) 设r行的基变量为w_l或Z_l(对于某个l≠S)，变量y_S进基，以r行为主行，y_S对应的列为主列，进行主元消去，如果离基变量是w_l，则令y_S= Z_l；如果离基变量是Z_l，则令y_S= w_l，转步骤2)。

4) 变量y_S进基，Z₀离基，以r行为主行，y_S对应的列为主列，进行主元消去，得到互补基本可行解，停止计算。

秩亏平差问题由于缺少必要的起算数据，引起法方程系数阵秩亏，导致结果有无穷多组解^[3]。附加条件法、广义逆法、拟稳平差法^[16]等经典算法处理秩亏问题在最小二乘的前提下增加等式基准约束条件，即假设最小范数准则$\hat{\boldsymbol{x}}_{P}^{\mathrm{T}} \boldsymbol{P}_{x} \hat{\boldsymbol{x}}_{P}$=min或基准约束条件$\boldsymbol{S}^{\mathrm{T}} \boldsymbol{P}_{x} \hat{\boldsymbol{x}}_{P}$=0成立^[3]，虽然克服了秩亏问题基准亏损的缺点，使秩亏问题得到唯一解，但未能充分利用先验信息。为了充分利用先验信息，并使平差结果唯一，本文构建不等式约束秩亏平差模型，并结合Karush-Kuhn-Tucker条件，将平差问题转化为求解线性互补的问题。由Lemke算法的步骤可知，本算法实现了准互补可行解到拟互补可行解的转换，且这种状态持续到Z₀离基才停止。最终，得到线性互补问题的互补基本可行解(w, Z)唯一。因此，Lemke算法在求解不等式约束秩亏平差的问题上具有一定的优势。下面结合模拟的GPS观测网算例，对本算法的有效性进行验证。

3 算例分析

图 1为模拟的GPS观测网，网中各条基线向量彼此独立，表 2给出了监测站点的真实坐标及近似坐标，表 3给出了基线向量和观测值的情况。

假设全程为等精度独立观测，则观测值的权阵P为单位权矩阵。A₁、A₂、A₃、A₄站的坐标改正数表示为${\widehat x_i}、{\widehat y_i}、{\widehat z_i}$(i=1, …, 4)。受外界条件的影响，经过一段时间观测，监测站点在3个方向上发生轻微移动：d_x≥0，d_y≥0，d_z≥0。移动的真值及与其他几种方法的比较见表 4。参数估值与真实值之间的残差平方和为$m=\|\hat{\boldsymbol{X}}-\tilde{\boldsymbol{X}}\|^{2}$，目标函数为$\boldsymbol{F}=(\boldsymbol{L}-\boldsymbol{A} \hat{\boldsymbol{X}})^{\mathrm{T}} \boldsymbol{P}(\boldsymbol{L}-\boldsymbol{A} \hat{\boldsymbol{X}})$。取$\mathit{\boldsymbol{\widehat X}}$ =[d_x, d_y, d_z]^T，误差方程为(单位m)：

$ \begin{aligned} &v_{1}=\hat{x}_{1}-\hat{x}_{4}+0.0107, v_{2}=\hat{y}_{1}-\hat{y}_{4}+0.0045, v_{3}=\hat{z}_{1}-\hat{z}_{4}-0.0298\\ &v_{4}=\hat{x}_{2}-\hat{x}_{4}-0.0137, v_{5}=\hat{y}_{2}-\hat{y}_{4}+0.0043, v_{6}=\hat{z}_{2}-\hat{z}_{4}+0.0155\\ &v_{7}=\hat{x}_{3}-\hat{x}_{4}+0.0411, v_{8}=\hat{y}_{3}-\hat{y}_{4}+0.0064, v_{9}=\hat{z}_{3}-\hat{z}_{4}-0.0047\\ &\begin{aligned} v_{10} &=\hat{x}_{2}-\hat{x}_{1}-0.0209, v_{11}=\hat{y}_{2}-\hat{y}_{1}-0.0162, v_{12}=\hat{z}_{2}-\hat{z}_{1}+0.0406 \\ v_{13} &=\hat{x}_{3}-\hat{x}_{1}-0.0054, v_{14}=\hat{y}_{3}-\hat{y}_{1}+0.0448, v_{15}=\hat{z}_{3}-\hat{z}_{1}+0.0260 \end{aligned}\\ &v_{16}=\hat{x}_{3}-\hat{x}_{2}+0.0325, v_{17}=\hat{y}_{3}-\hat{y}_{2}+0.0037, v_{18}=\hat{z}_{3}-\hat{z}_{2}-0.0259 \end{aligned} $

已知先验信息d_x≥0、d_y≥0，d_z≥0，根据GX ≤ h, X ≥0加入以下不等式约束条件：

$ \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&{ - 1}&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&{ - 1}&0&0&0&0&0\\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\ 0&0&0&0&0&0&0&0&0&0&0&{ - 1} \end{array}} \right], \mathit{\boldsymbol{h}} = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ 0\\ \ldots \\ 0 \end{array}} \right] $

在已知误差方程系数阵A、常数项L和不等式约束系数阵G、常数项h后，根据式(9)构建出Lemke算法需要的M和q。为了检验Lemke算法的可靠性和有效性，将本文方法与处理秩亏问题的附加条件法、广义逆法，拟稳平差法、高斯迭代法^[17]、罚函数法进行比较，结果见表 4(单位m)。

由表 4可以得出以下结论：

1) 未考虑未知参数先验信息的附加条件法和广义逆法偏离真值较大，效果不理想。

2) 因拟稳平差是用部分相对稳定点的未知参数范数最小代替所有未知参数范数最小的方法求解，故本文选取A₁和A₂站作为控制网的相对稳定点进行平差。从平差结果上看，拟稳平差法较前面2种方法的效率更高一些，但其结果不符合先验信息。

3) 高斯迭代法和罚函数法都考虑了未知参数的先验信息，平差结果十分接近。

4) 本文提出的不等式约束秩亏平差模型的Lemke算法充分利用了先验信息，对未知参数形成了有效约束，平差结果同真值最为接近且效果最好，显示出了其优越性。

4 结语

大地测量中，未知参数的先验信息反映了被测目标的几何或物理特征、内在相关性、测量精度等。在数据处理中，结合有效的先验信息对未知参数的等式或不等式约束进行平差，可以补充观测信息的不足，使得秩亏问题的解唯一。当前不等式约束平差的解算方法大多是针对非秩亏问题设计的，还不能解决平差模型的秩亏问题。利用先验信息解决秩亏问题，可以降低人的主观性，使平差结果更接近测绘工程的实际。本文就如何利用不等式约束先验信息解决秩亏问题，提出基于Lemke算法求解不等式约束秩亏平差问题的新算法，可以充分利用测绘工程中的先验信息，解决秩亏问题，降低平差模型的复杂性，丰富了平差算法理论。