-
大地测量数据处理中，参数本身附加的先验约束信息和物理约束信息常常被用来解决因观测信息不充分导致的病态问题[1-2]。如大地测量反演中，由先验模型或实际观测资料构成的反演计算约束条件[3]，变形监测中根据经验确定变形的范围和变形的方向[1-2]。在保持参数先验值和验后估值的统计、几何或物理意义下，这些附加信息可以转变成约束方程来降低模型的不适定性，保证解的稳定性和唯一性[4-5]。许多学者研究了如何利用参数内在的相关性、精度、几何或物理约束信息来保证参数解的唯一性和稳定性[6-13]。然而，这些研究利用的附加信息主要还是协方差阵以及一些容易处理的简单信息，如等式约束和不等式约束[14]。大地测量中，还存在一些用区间、集合以及椭球形式描述的约束信息，这些约束信息的统计性质不清楚，只有一个模糊变化的范围，如参数的可行空间、噪声变化的范围等[10]。大地测量反演问题实质上是一个附有约束条件的最优化问题，即在约束范围内求目标函数极小的问题[3]。一些复杂的附加先验信息，如单调性、光滑型、边界条件或源数据的误差界等加入反演模式中，显然可以提高反演的准确性。然而，先验信息的建模非常困难，实际工程中，为了简化平差模型，常利用区间约束或集合方法来描述复杂先验信息。
在测量数据处理中，带有区间约束的平差模型还没有得到广泛的应用；数学上，学者们给出的算法非常复杂[15-17]，难以在测量中进行应用。对于区间约束的平差模型，通常是将其转化为不等式约束的平差模型，但这种转换会增加更多的不等式约束，因而解算方法变得更加困难[13-14]。本文针对参数带有区间约束的平差模型，把平差问题转化成一个简单的二次规划问题，通过矩阵正则分裂，利用积极集对二次规划问题进行划分与重组，结合无约束共轭梯度优化算法，给出了带有区间约束的平差算法，并同时给出了参数解的精度评估。
-
带有区间约束的平差模型是指参数向量带有区间约束的平差模型，即：
$$ \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}}\\ {\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{X}} \le \mathit{\boldsymbol{u}}} \end{array}} \right. $$ (1) 式中，$\mathit{\boldsymbol{X}} = {({x_1}, {x_2} \cdots {x_n})^{\rm{T}}}$为n维未知参数；L为m维观测向量；A是m×n维系数矩阵（m≥n）；e为m维随机误差向量；$\mathit{\boldsymbol{l}} = {({l_1}, {l_2} \cdots {l_n})^{\rm{T}}}$和$\mathit{\boldsymbol{u}} = {({u_1}, {u_2} \cdots {u_n})^{\rm{T}}}$分别为X的下界和上界向量。
利用最小二乘平差准则求式（1）的参数解，实际上就是求$\;\left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right\|_\mathit{\boldsymbol{P}}^2$的最小值问题，其中P为式（1）的权矩阵，即：
$$ \left\{ {\begin{array}{*{20}{l}} {{\rm{min}}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)}\\ {\begin{array}{*{20}{c}} {{\rm{s}}.{\rm{t}}.}&{\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{X}} \le \mathit{\boldsymbol{u}}} \end{array}} \end{array}} \right. $$ (2) 利用拉格朗日方法求式（2）的极小值，令：
$$ \begin{array}{l} F\left( \mathit{\boldsymbol{X}} \right) = \left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right\|_\mathit{\boldsymbol{P}}^2 = {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\\ \left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right) = {\mathit{\boldsymbol{X}}^{\rm{T}}}\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}}} \right)\mathit{\boldsymbol{X}} - 2{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{PAX}} + {\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{PL}} \end{array} $$ (3) 因为LTPL为一个常量，二次规划问题（2）可以转化为如下的二次规划问题：
$$ \left\{ {\begin{array}{*{20}{l}} {{\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}}}\\ {\begin{array}{*{20}{c}} {{\rm{s}}.{\rm{t}}.}&{\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{X}} \le \mathit{\boldsymbol{u}}} \end{array}} \end{array}} \right. $$ (4) 式中，N=ATPA；c=-ATPL。
令$g\left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{NX}} + \mathit{\boldsymbol{c}}$表示F(X)在X处的梯度，对于规划问题(4)不带约束时，共轭梯度法是一种常用的解算方法, 它的迭代形式如下：
$$ {\mathit{\boldsymbol{X}}_{k + 1}} = {\mathit{\boldsymbol{X}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_k} $$ (5) 式中，${\mathit{\boldsymbol{X}}_k} = {({x_{k, 1}}, {x_{k, 2}} \cdots {x_{k, n}})^{\rm{T}}}$；dk称为搜索方向：
$$ {\mathit{\boldsymbol{d}}_k} = \left\{ {\begin{array}{*{20}{l}} { - g\left( {{\mathit{\boldsymbol{X}}_k}} \right), k = 0}\\ { - g\left( {{\mathit{\boldsymbol{X}}_k}} \right) + {\beta _k}{\mathit{\boldsymbol{d}}_k}, k \ge 1} \end{array}} \right. $$ (6) 式中，βk为共轭梯度参数，不同的βk的确定方法对应着不同的共轭梯度法。如果函数F(X)是严格凸二次函数, 那么在精确线搜索下各种算法是等价的，但对于一般的非线性函数, 各种算法在理论和数值上的表现很不相同。式（6）中的αk称为步长，${\alpha _k} = \mathop {{\rm{min}}}\limits_{\alpha > 0} \;F\left( {{\mathit{\boldsymbol{X}}_k} + \alpha {\mathit{\boldsymbol{d}}_k}} \right)$，它可以通过某种搜索算法得到。
共轭梯度法CGM-OC（conjugate gradient method-orthogonality correction）[18-20]的解算步骤如下：
1) 令k=0，初始点为X0，F0=F(X0)，$g\left( {{\mathit{\boldsymbol{X}}_0}} \right) = \mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_0} + \mathit{\boldsymbol{c}}$, ${\mathit{\boldsymbol{d}}_0} = - g\left( {{\mathit{\boldsymbol{X}}_0}} \right)$。
2) 线性搜索$\mathop {{\rm{min}}}\limits_{\alpha > 0} \;F\left( {{\mathit{\boldsymbol{X}}_k} + \alpha {\mathit{\boldsymbol{d}}_k}} \right) = F\left( {{\mathit{\boldsymbol{X}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_k}} \right)$，求得αk。
3) ${\mathit{\boldsymbol{X}}_{k + 1}} = {\mathit{\boldsymbol{X}}_k} + {\alpha _k}{\mathit{\boldsymbol{d}}_k}$。
4) $g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right) = \mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_{k + 1}} + \mathit{\boldsymbol{c}}$，新的搜索方向为：${\mathit{\boldsymbol{d}}_k} = - g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right) + {\beta _k}{\mathit{\boldsymbol{d}}_k}$，设：
$$ {{\bar \beta }_k} = \frac{{{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)g\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}, {v_k} = \frac{{{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}}{{{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)g\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right)}} $$ (7) 如果$\left| {{v_k}} \right| < \varepsilon \left( {\varepsilon = 0.175} \right)$，则βk=βk；否则，βk=βk(1-vk)。
5) 若满足收敛条件：$\left| {{\mathit{\boldsymbol{X}}_{k + 1}} - {\mathit{\boldsymbol{X}}_k}} \right| < {\delta _1}$，$\left| {F\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right) - F\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right| < {\delta _2}$，δ1和δ2为事先给定的精度，则停止；否则, 令k=k+1，转步骤2）。
对于上述CGM-OC算法得到的序列{Xk}，文献[18]证明了它所对应的函数序列f(Xk)是严格单调递减序列。如果变量X没有约束，直接利用CGM-OC算法可以求出最优解。当X有区间约束限制时，在迭代过程中，还需作一些必要的改进，使迭代产生的解不至于越界。
-
积极集方法是求解问题(4)的一种比较好的方法，在点$\mathit{\boldsymbol{X}} = {({x_1}, {x_2} \cdots {x_n})^{\rm{T}}}$的积极集定义为$A\left( \mathit{\boldsymbol{X}} \right) = L\left( \mathit{\boldsymbol{X}} \right) \cup U\left( \mathit{\boldsymbol{X}} \right)$，其中$L\left( \mathit{\boldsymbol{X}} \right) = \left\{ {i:{x_i} = {l_i}} \right\}$，$U\left( \mathit{\boldsymbol{X}} \right) = \left\{ {i:{x_i} = {u_i}} \right\}$。设二次规划问题(4)的可行域$D = \left\{ {\mathit{\boldsymbol{X}} \in {{\bf{R}}^n}:\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{X}} \le \mathit{\boldsymbol{u}}} \right\}$，若$\mathit{\boldsymbol{X}} = {({x_1}, {x_2} \cdots {x_n})^{\rm{T}}} \in D$，满足：
$$ \left\{ {\begin{array}{*{20}{c}} {{x_i} = {l_i}, {g_i}\left( \mathit{\boldsymbol{X}} \right) > 0{\rm{}}}\\ {{l_i} < {x_i} < {u_i}, {g_i}\left( \mathit{\boldsymbol{X}} \right) = 0}\\ {{x_i} = {u_i}, {g_i}\left( \mathit{\boldsymbol{X}} \right) < 0} \end{array}} \right. $$ (8) 称X为问题(4)的一个稳定点或KKT（Karush-Kuhn-Tucher）点，其中gi(X)表示g(X)的第i个分量。X为问题(4)的最优解的充分必要条件是X为问题(4)的KKT点[21]。一般地，X是问题(4)的KKT点，它也是如下问题的KKT点[21]：
$$ \left\{ {\begin{array}{*{20}{l}} {{\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}}}\\ {{\rm{s}}.{\rm{t}}.{\rm{}}\;\;\;\;{x_i} = \left\{ {\begin{array}{*{20}{c}} {{l_i}, i \in L\left( {\mathit{\boldsymbol{\bar X}}} \right)}\\ {{u_i}, i \in U\left( {\mathit{\boldsymbol{\bar X}}} \right)} \end{array}} \right.} \end{array}} \right. $$ (9) 换言之，带有区间约束的优化问题(4)等价于带有等式约束的优化问题(9)，但是，由于问题(9)依赖于问题(4)的解X，因此不能直接通过求解问题(9)的方式来解算问题(4)。求解区间约束问题(4)的积极集法的基本思想是从问题(4)的初始可行解X0出发，如果X0不是问题(4)的最优解，则用L(X0)和U(X0)来代替问题中的L(X)和U(X)，从而得到X1，依次类推，直到获得最优解。定义指标集：
$$ \begin{array}{*{20}{c}} {I = \left\{ {i:{x_i} = {l_i}且{g_i}\left( \mathit{\boldsymbol{X}} \right) > 0} \right\} \cup \left\{ {i:{x_i} = {u_i}且{g_i}\left( \mathit{\boldsymbol{X}} \right) < 0} \right\}}\\ {{\rm{}}J = \left\{ {1, 2 \cdots n} \right\} - I} \end{array} $$ (10) 对应指标集I和J重新划分g = NX + c中的参数向量X、矩阵N、向量c、向量g，有：
$$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}_{II}}}&{\mathit{\boldsymbol{N}}_{JI}^{\rm{T}}}\\ {{\mathit{\boldsymbol{N}}_{JI}}}&{{\mathit{\boldsymbol{N}}_{JJ}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_I}}\\ {{\mathit{\boldsymbol{X}}_J}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{C}}_I}}\\ {{\mathit{\boldsymbol{C}}_J}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{g}}_I}}\\ {{\mathit{\boldsymbol{g}}_J}} \end{array}} \right] $$ (11) 式中，XI是X中对应下标$i \in I$的分量xi构成的向量，xi=li或xi=ui，XI的值保持固定，而让XJ进行迭代。令gJ=0，由式(11)有：
$$ {\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{X}}_J} = - {\mathit{\boldsymbol{C}}_J} - {\mathit{\boldsymbol{N}}_{JI}}{\mathit{\boldsymbol{X}}_I} $$ (12) 由于NJJ是N的主子阵，由N的正定可知，NJJ也是正定的，则有：
$$ \begin{array}{*{20}{c}} {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}} = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}}_I^{\rm{T}}}&{\mathit{\boldsymbol{X}}_J^{\rm{T}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}_{II}}}&{\mathit{\boldsymbol{N}}_{JI}^{\rm{T}}}\\ {{\mathit{\boldsymbol{N}}_{JI}}}&{{\mathit{\boldsymbol{N}}_{JJ}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_I}}\\ {{\mathit{\boldsymbol{X}}_J}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{C}}_I^{\rm{T}}}&{\mathit{\boldsymbol{C}}_J^{\rm{T}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_I}}\\ {{\mathit{\boldsymbol{X}}_J}} \end{array}} \right] = }\\ {\left[ {\frac{1}{2}\mathit{\boldsymbol{X}}_J^{\rm{T}}{\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{X}}_J} + {{\left( {{\mathit{\boldsymbol{N}}_{JI}}{\mathit{\boldsymbol{X}}_I} + {\mathit{\boldsymbol{C}}_J}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{X}}_J}} \right] + \left( {\frac{1}{2}\mathit{\boldsymbol{X}}_I^{\rm{T}}{\mathit{\boldsymbol{N}}_{II}}{\mathit{\boldsymbol{X}}_I} + \mathit{\boldsymbol{C}}_I^{\rm{T}}{\mathit{\boldsymbol{X}}_I}} \right)} \end{array} $$ (13) 对已选择的指标集I，式(13)中等号右边的第2项不变，第1项看作无约束的二次规划问题，可以利用CGM-OC来进行内部迭代求解。当X有区间约束限制时（$\mathit{\boldsymbol{l}} \le \mathit{\boldsymbol{X}} \le \mathit{\boldsymbol{u}}$），在迭代过程中，如果XS（$S \subseteq J$）达到或违反上、下界约束，便缩短搜索步长，使XS达到其上界或下界，然后把S加到I中，再对矩阵和变量进行重新划分，重新开始CGM-OC迭代，直到式(12)的迭代结果中没有变量违反区间约束条件，就完成了内部迭代。这时当${l_i} < {x_i} < {u_i}$($i \in J$)时，有gi(X)=0，Xk满足最优性条件式(8)。接着对所求的Xk开始新的外部迭代，一旦Xk的指标集Ik=Ik-1，那么Xk最优性条件完全满足，算法停止，得到最优解，否则开始新的内部迭代。
共轭梯度积极集算法为：选择初始点为X0，使得$\mathit{\boldsymbol{l}} \le {\mathit{\boldsymbol{X}}_0} \le \mathit{\boldsymbol{u}}$，令k=0，$I = \left\{ {1, 2 \cdots n} \right\}$，步骤如下：
1) 外部迭代：k=k+1，令Xk=Xk-1，g(Xk)=NXk +c，Ik-1=I，有：
$$ {I_k} = \left\{ {i:{x_{k, i}} = {l_i}且{g_i}\left( {{\mathit{\boldsymbol{X}}_k}} \right) > 0} \right\} \cup \left\{ {i:{x_{k, i}} = {u_i}且{g_i}\left( {{\mathit{\boldsymbol{X}}_k}} \right) < 0} \right\} $$ (14) 若Ik=Ik-1，则算法停止，找到最优解；否则，令I=Ik，开始内部迭代。
2) 内部迭代
ⅰ)$J = \left\{ {1, 2 \cdots n} \right\} - I$，按指标集重新划分${\mathit{\boldsymbol{X}}_k} \to \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_{{I_k}}}}\\ {{\mathit{\boldsymbol{X}}_{{J_k}}}} \end{array}} \right]$，$\mathit{\boldsymbol{c}} \to \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{c}}_I}}\\ {{\mathit{\boldsymbol{c}}_J}} \end{array}} \right]$，$g\left( {{\mathit{\boldsymbol{X}}_k}} \right) \to \left[ {\begin{array}{*{20}{c}} {{g_I}\left( {{\mathit{\boldsymbol{X}}_k}} \right)}\\ {{g_J}\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \end{array}} \right]$，N→$\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{N}}_{II}}}&{\mathit{\boldsymbol{N}}_{JI}^{\rm{T}}}\\ {{\mathit{\boldsymbol{N}}_{JI}}}&{{\mathit{\boldsymbol{N}}_{JJ}}} \end{array}} \right]$，NJJ是s×s阶正定矩阵，利用CGM-OC算法解方程组${\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{X}}_{{J_k}}} = - {\mathit{\boldsymbol{c}}_J} - {\mathit{\boldsymbol{N}}_{JI}}{\mathit{\boldsymbol{X}}_{{I_k}}}$，序列${\mathit{\boldsymbol{Z}}_q} = {({z_{q, 1}}, {z_{q, 2}} \cdots {z_{q, s}})^{\rm{T}}}$是解向量XJk的近似值，搜索方向${\mathit{\boldsymbol{p}}_q} = {({p_{q, 1}}, {p_{q, 2}} \cdots {p_{q, s}})^{\rm{T}}}$，方程组的残向量${\mathit{\boldsymbol{r}}_q} = {({r_{q, 1}}, {r_{q, 2}} \cdots {r_{q, s}})^{\rm{T}}}$，设q=0，${\mathit{\boldsymbol{Z}}_0} = {\mathit{\boldsymbol{X}}_{{J_k}}}$，${\mathit{\boldsymbol{p}}_0} = {\mathit{\boldsymbol{r}}_0} = - {\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{Z}}_0} - {\mathit{\boldsymbol{N}}_{JI}}{\mathit{\boldsymbol{X}}_{{I_k}}} - {\mathit{\boldsymbol{c}}_J}$。
ⅱ)计算新的迭代点和残量，acg是方向pq上的一维搜索步长：
$$ \begin{array}{*{20}{c}} {{a_{cg}} = \frac{{\mathit{\boldsymbol{r}}_q^{\rm{T}}{\mathit{\boldsymbol{r}}_q}}}{{\mathit{\boldsymbol{p}}_q^{\rm{T}}{\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{p}}_q}}}}\\ {{a_m} = {\rm{min}}\left\{ {\mathop {{\rm{min}}}\limits_{\begin{array}{*{20}{c}} {j = 1, 2 \cdots s}\\ {{p_{q, j}} < 0} \end{array}} \frac{{{l_j} - {z_{q, j}}}}{{{p_{q, j}}}}, \mathop {{\rm{min}}}\limits_{\begin{array}{*{20}{c}} {j = 1, 2 \cdots s}\\ {{p_{q, j}} > 0} \end{array}} \frac{{{u_j} - {z_{q, j}}}}{{{p_{q, j}}}}} \right\}} \end{array} $$ 所求步长${a_q} = {\rm{min}}\left\{ {{a_{cg}}, {a_m}} \right\}$，${\mathit{\boldsymbol{Z}}_{q + 1}} = {\mathit{\boldsymbol{Z}}_q} + {a_q}{\mathit{\boldsymbol{p}}_q}$，${\mathit{\boldsymbol{r}}_{q + 1}} = {\mathit{\boldsymbol{r}}_q} - {a_q}{\mathit{\boldsymbol{N}}_{JJ}}{\mathit{\boldsymbol{p}}_q}$。
ⅲ)内部迭代终止准则
① 如果rq+1=0，XJk=Zq+1，重新开始外部迭代。
② 如果$\left\{ {j:{z_{\left( {q + 1} \right), j}} = {l_j}或{u_j}} \right\} = \emptyset $，转④。
③ ${\mathit{\boldsymbol{X}}_{{J_k}}} = {\mathit{\boldsymbol{Z}}_{q + 1}}$，$I = \left\{ {i:{x_{k, i}} = {l_i}或{u_i}} \right\}$，如果$I = \left\{ {1, 2 \cdots n} \right\}$，重新开始外部迭代，否则重新开始内部迭代。
④ 计算新的搜索方向:${\mathit{\boldsymbol{p}}_{q + 1}} = {\mathit{\boldsymbol{r}}_{q + 1}} + {b_q}{\mathit{\boldsymbol{p}}_q}$，
$$ {\beta _q} = \frac{{\mathit{\boldsymbol{r}}_{q + 1}^{\rm{T}}{\mathit{\boldsymbol{r}}_{q + 1}}}}{{\mathit{\boldsymbol{r}}_q^{\rm{T}}{\mathit{\boldsymbol{r}}_q}}}, {v_q} = \frac{{\mathit{\boldsymbol{r}}_q^{\rm{T}}{\mathit{\boldsymbol{r}}_{q + 1}}}}{{\mathit{\boldsymbol{r}}_{q + 1}^{\rm{T}}{\mathit{\boldsymbol{r}}_{q + 1}}}} $$ (15) 如果$\left| {{v_q}} \right| < \varepsilon {\rm{}}\left( {\varepsilon = 0.175} \right)$，则bq=βq；否则，${b_q} = {\beta _q}\left( {1 - {v_k}} \right)$，q=q+1，转②。
在上述算法中，式（13）等号右边的第1项是利用CGM-OC来求解的，其目标函数值严格单调递减，因为在每次进行内部迭代时，至少要进行一次CGM-OC迭代，故在内部迭代中求解式(13)时，不会出现重复现象。而指标集$I$的选取是有限的，则对应的方程(13)的个数有限，故算法必在有限步内终止。
-
从上面求解最优解X的算法可知，在以下区间约束集合：
$$ I = \left\{ {i:{x_i} = {l_i}且{g_i}\left( \mathit{\boldsymbol{X}} \right) > 0} \right\}\mathop \cup \left\{ {i:{x_i} = {u_i}且{g_i}\left( \mathit{\boldsymbol{X}} \right) < 0} \right\} $$ 可以找到$I$的积极集$S = I - \left( {L\left( {\mathit{\boldsymbol{\bar X}}} \right)\mathop \cup U\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right)$，从而得到等式约束的平差模型：
$$ \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}}\\ {{\mathit{\boldsymbol{B}}_S}\mathit{\boldsymbol{X}} = {\mathit{\boldsymbol{d}}_S}} \end{array}} \right. $$ (16) 式中，约束方程系数矩阵BS为s×n阶矩阵，$s$为集合$S$的元素个数，${{\mathit{\boldsymbol{B}}_S}}$由单位向量${\mathit{\boldsymbol{b}}_i}$组成（${\mathit{\boldsymbol{b}}_i}$除了第$i$个分量为1外，其余分量为0），$i \in S$；dS为第i分量等于li或ui的向量。设${{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PL}}$为模型（16）的无约束最小二乘估计，其单位权方差的无偏估计为：
$$ \sigma _S^2 = \frac{{{{\left( {\mathit{\boldsymbol{L}} - A{{\hat X}_S}} \right)}^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{{\hat X}_S}} \right)}}{{m - n - s}} $$ (17) 式中，$s$为式（16）中${{\mathit{\boldsymbol{B}}_S}}$的行数；$m$为观测向量的维数；$n$为参向量的维数。由带有等式约束的平差算法[22]可知，式（16）的参数解${{{\mathit{\boldsymbol{\hat X}}}_S}}$可以表示为：
$$ {{\mathit{\boldsymbol{\hat X}}}_S} = {{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\mathit{\boldsymbol{P}}_S}\left( {{\mathit{\boldsymbol{B}}_S}{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\mathit{\boldsymbol{d}}_S}} \right) = {{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{\hat \phi}}} $$ (18) $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\hat X}_{{\rm{LS}}}}}} = {\rm{E}}\left[ {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\rm{E}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}} \right)} \right]{{\left[ {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\rm{E}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}} \right)} \right]}^{\rm{T}}} = }\\ {{\sigma ^2}{{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{ - 1}}} \end{array} $$ (19) 其中，
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_S} = {{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{ - 1}}\mathit{\boldsymbol{B}}_S^{\rm{T}}\mathit{\boldsymbol{Q}}_S^{ - 1}}\\ {{\mathit{\boldsymbol{Q}}_S} = {\mathit{\boldsymbol{B}}_S}{{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{ - 1}}\mathit{\boldsymbol{B}}_S^{\rm{T}}}\\ {{\mathit{\boldsymbol{\hat \phi}}} = {\mathit{\boldsymbol{B}}_S}{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\mathit{\boldsymbol{d}}_S}} \end{array} $$ 因为$S$无法事先确定，式（18）只是一个事后的形式解，只能在作精度分析时使用。由于${\rm{E}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{\rm{S}}}} \right) = \mathit{\boldsymbol{X}}$，令${\mathit{\boldsymbol{\phi}}} = {\mathit{\boldsymbol{B}}_S}\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{d}}_S}$，利用${\mathit{\boldsymbol{\hat \phi}}} - {\mathit{\boldsymbol{\phi}}} = {\mathit{\boldsymbol{B}}_S}\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - \mathit{\boldsymbol{X}}} \right)$可以计算出：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\hat X}_S}}} = {\rm{E}}\left[ {\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - \mathit{\boldsymbol{X}}} \right) - {\mathit{\boldsymbol{P}}_S}\left( {{\mathit{\boldsymbol{\hat \phi}}} - {\mathit{\boldsymbol{\phi}}} } \right)} \right]{{\left[ {\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - \mathit{\boldsymbol{X}}} \right) - {\mathit{\boldsymbol{P}}_S}\left( {{\mathit{\boldsymbol{\hat \phi}}} - {\mathit{\boldsymbol{\phi}}}} \right)} \right]}^{\rm{T}}} = }\\ {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}} - {\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{B}}_S}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}} - {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}}\mathit{\boldsymbol{B}}_S^{\rm{T}}\mathit{\boldsymbol{P}}_S^{\rm{T}} + {\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{B}}_S}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}}\mathit{\boldsymbol{B}}_S^{\rm{T}}\mathit{\boldsymbol{P}}_S^{\rm{T}}} \end{array} $$ (20) 又因为
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{B}}_S}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}} = {\sigma ^2}{\mathit{\boldsymbol{Q}}_S}\mathit{\boldsymbol{P}}_S^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}}\mathit{\boldsymbol{B}}_S^{\rm{T}} = {\sigma ^2}{\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{Q}}_S}}\\ {{\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{B}}_S}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}}\mathit{\boldsymbol{B}}_S^{\rm{T}}\mathit{\boldsymbol{P}}_S^{\rm{T}} = {\sigma ^2}{\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{Q}}_S}\mathit{\boldsymbol{P}}_S^{\rm{T}}} \end{array} $$ 故
$$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_S}}} = {\sigma ^2}\left[ {{{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}})}^{ - 1}} - {\mathit{\boldsymbol{P}}_S}{\mathit{\boldsymbol{Q}}_S}\mathit{\boldsymbol{P}}_S^{\rm{T}}} \right] $$ (21) 式（21）可以用来对具有线性区间约束的平差模型的参数估计的精度进行评估计算。
-
如图 1所示的测边网，P1、P2为已知点，其坐标分别为(48 580.285, 600 500.496)和(47 943.013, 66 225.845)，单位为m。假设算例1中的点P3、P4、P5、P6的真实坐标如表 1所示，边长观测值如表 2所示，由前期的观测结果可知P3、P4、P5、P6的近似坐标(见表 1)以及它们相应的点位精度。先验约束区间D2是通过点位精度进行推算的，D1、D3是特别设定的约束区间，其中，D1是把区间约束放大了6倍，是为了说明当约束区间不准确时对参数估计的影响；D3是在D2的基础上人为增加的一个约束，目的是检测算法能否真正找到这个有效约束。
图 1 测边三角网
Figure 1. Distance-Measuring Triangulation Network
表 1 真实坐标与近似坐标/m
Table 1. True Coordinates and Approximate Coordinates/m
坐标 真实坐标 近似坐标 P3 P4 P5 P6 P3 P4 P5 P6 x 53 743.151 48 681.398 43 767.234 40 843.239 53 743.674 48 680.496 43 768.794 40 840.905 y 61 003.810 55 018.270 57 968.590 64 867.876 61 006.568 55 018.806 57 966.087 64 870.541 表 2 算例1边长观测值
Table 2. Edge Length Observation Values of Example 1
边号 观测边长/m 1 5 760.706 2 7 804.611 3 5 187.314 4 7 838.890 5 5 483.157 6 5 731.827 7 5 438.409 8 7 493.319 9 8 884.539 10 7 228.509 设$F\left( \mathit{\boldsymbol{X}} \right) = \left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}} }\right\|_\mathit{\boldsymbol{P}}^2$，$m = {\left( {\mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{X}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat X}} - \mathit{\boldsymbol{X}}} \right)$。相对于近似坐标的改正数构成的未知向量为：
$$ \mathit{\boldsymbol{X}} = {[{x_3}, {\rm{}}{y_3}, {\rm{}}{x_4}, {\rm{}}{y_4}, {\rm{}}{x_5}, {\rm{}}{y_5}, {\rm{}}{x_6}, {\rm{}}{y_6}]^{\rm{T}}} $$ 通过计算可以得到相应的系数矩阵和观测向量如下：
$$ {\mathit{\boldsymbol{A}}_1} = \left[ {\begin{array}{*{20}{c}} {0.7434}&{ - 0.6689}&0&0&0&0&0&0\\ {0.9952}&{0.0975}&0&0&0&0&0&0\\ {0.6457}&{0.7636}&{ - 0.6457}&{ - 0.7636}&0&0&0&0\\ 0&0&{0.0183}&{ - 0.9998}&0&0&0&0\\ 0&0&{0.8575}&{ - 0.5145}&{ - 0.8575}&{0.5145}&0&0\\ 0&0&0&0&{ - 0.8848}&{ - 0.4660}&0&0\\ 0&0&0&0&{0.3904}&{ - 0.9206}&{ - 0.3904}&{0.9206}\\ 0&0&0&0&0&0&{ - 0.9708}&{0.4917}\\ 0&0&0&0&0&0&{0.9823}&{0.1874}\\ \end{array}} \right] $$ $$ {\mathit{\boldsymbol{L}}_1} = {\left[ {1.501{\rm{}}3, - 0.815{\rm{}}7, - 2.606{\rm{}}1, 0.551{\rm{}}1, 3.714{\rm{}}0, 0.240{\rm{}}7, - 6.282{\rm{}}3, - 3.390{\rm{}}8, - 1.760{\rm{}}3} \right]^{\rm{T}}} $$ 方案1：设${\mathit{\boldsymbol{d}}_1} = {\left[ {18, 18, 18, 18, 18, 18, 18, 18} \right]^{\rm{T}}}$，约束区间为：
$$ {D_1} = \left\{ {\mathit{\boldsymbol{X}}: - {\mathit{\boldsymbol{d}}_1} \le \mathit{\boldsymbol{X}} \le {\mathit{\boldsymbol{d}}_1}} \right\} $$ 得到的区间约束参数估计为${{{\mathit{\boldsymbol{\hat X}}}_1}}$。
方案2：设${\mathit{\boldsymbol{d}}_2} = {\left[ {3, 3, 3, 3, 3, 3, 3, 3} \right]^{\rm{T}}}$, 约束区间为：
$$ {D_2} = \left\{ {\mathit{\boldsymbol{X}}: - {\mathit{\boldsymbol{d}}_2} \le \mathit{\boldsymbol{X}} \le {\mathit{\boldsymbol{d}}_2}} \right\} $$ 得到的区间约束参数估计为${{{\mathit{\boldsymbol{\hat X}}}_2}}$。
方案3：设${\mathit{\boldsymbol{l}}_3} = [ - 3, - 3, - 3, - 3, - 3, - 3, - 3, $${\left. { - 2.6650} \right]^{\rm{T}}}$, ${\mathit{\boldsymbol{u}}_3} = {\left[ {3, 3, 3, 3, 3, 3, 3, 3} \right]^{\rm{T}}}$, 约束区间为：
$$ {D_3} = \left\{ {\mathit{\boldsymbol{X}}:{\mathit{\boldsymbol{l}}_3} \le \mathit{\boldsymbol{X}} \le {\mathit{\boldsymbol{u}}_3}} \right\} $$ 得到的区间约束参数估计为${{\mathit{\boldsymbol{\hat X}}}_3}$。
算法的外部迭代终止条件是得到的积极集${I_k} = {I_{k - 1}}$, 内部迭代算法的终止条件设置为：$\left| {{\mathit{\boldsymbol{X}}_k} - {\mathit{\boldsymbol{X}}_{k - 1}}} \right| < 0.000{\rm{}}01$。
算法结果分析与说明如下：
1）方案1和方案2的约束区间都比较大，算法中没有得到有效约束，这时得到的解算结果与最小二乘估计的结果一致，换句话说，就是最优解在区间D1和D2之中，因此有${{\mathit{\boldsymbol{\hat X}}}_1} = {{\mathit{\boldsymbol{\hat X}}}_2} = {{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}$。
2）方案3中，约束区间D3的一个端点与真值相同（${y_8} = - 2.665{\rm{}}0$），人为地生成一个有效约束，与方案1、2不同的是，这个有效约束不是算法自动寻找的，但计算结果与有效约束非常一致，这进一步说明了算法的有效性。
3）从方案3的结果可以看出，此时的有效约束为${y_8} = - 2.665{\rm{}}0$，相对于式（16），有${\mathit{\boldsymbol{B}}_S} = \left[ {0, 0, 0, 0, 0, 0, 0, 1} \right]$，${d_S} = - 2.665{\rm{}}0$，由式（18）有：
$$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_S} = {{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA)}}}^{ - 1}}\mathit{\boldsymbol{B}}_S^{\rm{T}}\left( {{\mathit{\boldsymbol{B}}_S}{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}} - {\mathit{\boldsymbol{d}}_S}} \right) = }\\ {[ - 0.522{\rm{}}3, - 2.801{\rm{}}2, 0.883{\rm{}}1, - 0.554{\rm{}}6, }\\ { - 1.598{\rm{}}7, 2.502{\rm{}}4, 2.336{\rm{}}9{]^{\rm{T}}}} \end{array} $$ 与论文算法的结果一致，由式（21）可以对${{\mathit{\boldsymbol{\hat X}}}_3}$进行精度评估。
4）算法的收敛速度非常快，方案1和方案2都只进行了2次外部迭代和8次内部迭代（k=1时）；方案3进行了3次外部迭代和7次内部迭代（k=2时），其中k表示外部迭代的序号。
算例1说明对于列满秩的系数矩阵A，法矩阵N=ATPA为正定矩阵，二次规划问题有唯一的最优解，从表 3的计算结果可以看出，虽然区间约束不同，本文算法却与最小二乘算法结果一致，它们都非常接近真值，说明观测信息充分时，不需要补充先验信息，只要最小二乘解的最优解在可行域（约束区间）$D$上，由唯一性可知参数解与最小二乘是一致的（先验信息没有起到作用）。
表 3 系数矩阵非病态时算法结果比较
Table 3. Results Comparison when Normal Equation Matrix is non Ill-posed
比较项 真值 最小二乘算法 方案1 方案2 方案3 ${\mathit{\boldsymbol{\hat X}}}$ -0.523 0 -0.516 0 -0.516 0 -0.516 0 -0.522 3 -2.758 0 -2.786 0 -2.786 0 -2.786 0 -2.801 2 0.902 0 0.923 3 0.923 3 0.923 3 0.883 1 -0.536 0 -0.560 6 -0.560 6 -0.560 6 -0.554 6 -1.560 0 -1.577 4 -1.577 4 -1.577 4 -1.598 7 2.503 0 2.456 2 2.456 2 2.456 2 2.502 4 2.334 0 2.327 0 2.327 0 2.327 0 2.336 9 -2.665 0 -2.728 6 -2.728 6 -2.728 6 -2.665 0 $F\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.008 5 0.004 0 0.004 0 0.004 0 0.005 1 m 0 0.008 5 0.008 5 0.008 5 0.004 1 -
修改算例1中已知点P2的坐标为(48 570.013, 60 555.845)，让它非常靠近P1点，导致算法1中的系数矩阵病态，用来分析病态问题中算法的性能。此时相应的边长观测也会发生变化，见表 4，相应的平差模型的系数矩阵和观测向量如下：
$$ {\mathit{\boldsymbol{A}}_2} = \left[ {\begin{array}{*{20}{c}} {0.995 8}&{0.092 0}&0&0&0&0&0&0\\ {0.995 2}&{0.097 5}&0&0&0&0&0&0\\ {0.645 7}&{0.763 6}&{-0.645 7}&{-0.763 6}&0&0&0&0\\ 0&0&{0.018 3}&{-0.999 8}&0&0&0&0\\ 0&0&{0.857 5}&{-0.514 5}&{-0.857 5}&{0.514 5}&0&0\\ 0&0&0&0&{ -0.884 8}&{-0.466 0}&0&0\\ 0&0&0&0&{0.390 4}&{-0.920 6}&{-0.390 4}&{0.920 6}\\ 0&0&0&0&0&0&{-0.870 8}&{0.491 7}\\ 0&0&0&0&0&0&{-0.871 0}&{0.491 3}\\ \end{array}} \right] $$ $$ {\mathit{\boldsymbol{L}}_2} = {\left[ { - 0.775{\rm{}}1, - 0.739{\rm{}}6, - 2.629{\rm{}}5, 0.555{\rm{}}6, 3.643{\rm{}}0, 0.214{\rm{}}7, - 6.285{\rm{}}7, - 3.326{\rm{}}6, - 3.324{\rm{}}8} \right]^{\rm{T}}} $$ 表 4 算例2边长观测值
Table 4. Edge Length Observation Values of Example 2
边号 观测边长/m 1 45.075 2 5 222.056 3 5 187.391 4 7 838.867 5 5 483.162 6 5 731.756 7 5 438.383 8 7 493.316 9 8 884.603 10 8 839.687 采用与算例1相同的3种方案，算法结果分析与说明如下：
1）算例2中，虽然$\mathit{\boldsymbol{N}} = {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}}$是一个正定矩阵(此例中$\mathit{\boldsymbol{P}}$为单位矩阵)，但有一特征值接近0，$\mathit{\boldsymbol{N}}$的条件数为1.3533×106，属于病态法方程，这时最小二乘解偏离真值较大。
2）方案1中的可行区间D1由于设置过大, 没有形成有效约束，从表 5的结果中可以看出，由方案1得到的参数估计解${{\mathit{\boldsymbol{\hat X}}}_1}$与最小二乘解${{{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}}$结果一致。
表 5 法矩阵病态时算法的结果比较
Table 5. Results Comparison when Normal Equation Matrix is Ill-posed
比较项 真值 最小二乘方法 截断奇异值法 岭估计法
(L曲线法)方案1 方案2 方案3 ${\mathit{\boldsymbol{\hat X}}}$ -0.523 0 -1.347 3 -0.535 0 -0.552 7 -1.347 3 -0.510 5 -0.489 3 -2.758 0 6.162 5 -2.372 9 -2.263 9 6.162 5 -2.630 3 -2.853 2 0.902 0 10.443 9 1.358 5 1.382 6 10.443 9 1.084 4 0.847 1 -0.536 0 -0.364 5 -0.530 7 -0.497 5 -0.364 5 -0.535 6 -0.540 0 -1.560 0 2.869 2 -1.331 3 -1.328 4 2.869 2 -1.458 1 -1.567 8 2.503 0 -5.908 4 2.067 1 2.080 8 -5.908 4 2.307 9 2.516 2 2.334 0 -5.331 9 1.907 1 1.935 5 -5.331 9 2.125 7 2.314 7 -2.665 0 -16.214 1 -3.387 3 -3.328 9 -16.214 1 -3.000 0 -2.665 0 $F\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.004 3 1.686 1×10-5 0.001 2 0.004 4 1.686 1×10-5 0.001 2 0.001 3 m 0 504.044 1 1.303 2 1.308 9 504.044 1 0.253 8 0.013 8 3）方案2中，约束区间D2与约束区间D1相比已经有了较大的缩小，从计算结果可以看出，已经产生了一个有效约束（此时积极集I={8}为非空集合），有效约束显著改善了模型的病态性，而最小二乘解的估计${{{\mathit{\boldsymbol{\hat X}}}_{{{LS}}}}}$不在约束区间内，与先验信息明显不符，严重失真，这说明当平差模型出现严重病态时，加入区间约束可以显著改善其病态性。
4）方案3中，约束区间D3的一个端点与真值相同（y6=-2.6650），人为地生成了一个有效约束，与方案2不同的是，这个有效约束不是算法自动寻找的，但计算结果与有效约束非常一致，这进一步说明了算法的有效性。
5）算法的收敛速度非常快，方案1运算进行的外部迭代数k=2, 且仅在k=1进行了内部迭代，次数q=10；方案2、方案3运算进行的外部迭代数k=3, 且仅在k=2时进行了内部迭代，次数q=7。
6）同算例1一样，可以用式（16）计算出${{\mathit{\boldsymbol{\hat X}}}_1}$、${{\mathit{\boldsymbol{\hat X}}}_2}$、${{\mathit{\boldsymbol{\hat X}}}_3}$，并用式（21）对它们进行精度评估。
此算例说明，利用参数的先验约束信息可以改善平差模型的病态性。因此，充分利用参数的先验信息对参数进行约束，可以弥补观测不充分所导致的病态问题，保证解的稳定性。
-
大地测量中存在一些用区间、集合以及椭球形式描述的约束信息，这些约束信息的统计性质不清楚，只有一个模糊变化的范围，如参数的可行空间、噪声变化的范围等。充分利用它们可以对部分参数进行约束，从而保证参数解的唯一性和稳定性。本文利用区间约束来描述参数的先验信息，提供了一种基于共轭梯度积极集的平差新方法，其算法简单。通过实例说明，利用参数的先验约束信息可以改善平差模型的病态性。在测量数据处理中，充分利用参数的先验信息建立平差模型，可以弥补观测不充分导致的病态问题，提高参数估计的效率。
An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints
-
摘要: 主要研究参数带有区间约束的平差算法，通过把平差问题转化成一个带有区间约束的二次规划问题，利用积极集对二次规划问题进行划分与重组，结合无约束共轭梯度优化算法，给出了带有区间约束的平差算法，并同时给出了参数解的精度评估。由于投影梯度法可以迅速改变积极约束集的构成，新的算法比经典的积极集法效率更高，可以降低模型的不适定性，保持参数先验信息中的统计、几何或物理意义，适合于求解大规模的带有区间约束的平差问题。Abstract: This paper mainly studies the adjustment model with interval constraints, in which the adjustment problem is transformed into a quadratic programming problem with interval constraints. A new adjustment algorithm with interval constraints is presented, and the accuracy of the parameter solution is evaluated, in which the active set algorithm and unconstrained conjugate gradient optimization algorithm are used to partition and reconstruct the quadratic problems. Because the projection gradient method can rapidly change the composition of the active constraint set, the new algorithm is more efficient than the classical active set method, which can reduce the uncertainty of the model, and maintain the statistical, geometric or physical meaning of priori information, and is suitable for solving large-scale adjustment problems with interval constraints.
-
Key words:
- active set /
- conjugate gradients method /
- interval constraints /
- adjustment model /
- ill-posed problem
-
表 1 真实坐标与近似坐标/m
Table 1. True Coordinates and Approximate Coordinates/m
坐标 真实坐标 近似坐标 P3 P4 P5 P6 P3 P4 P5 P6 x 53 743.151 48 681.398 43 767.234 40 843.239 53 743.674 48 680.496 43 768.794 40 840.905 y 61 003.810 55 018.270 57 968.590 64 867.876 61 006.568 55 018.806 57 966.087 64 870.541 
下载: 导出CSV
表 2 算例1边长观测值
Table 2. Edge Length Observation Values of Example 1
边号 观测边长/m 1 5 760.706 2 7 804.611 3 5 187.314 4 7 838.890 5 5 483.157 6 5 731.827 7 5 438.409 8 7 493.319 9 8 884.539 10 7 228.509
下载: 导出CSV
表 3 系数矩阵非病态时算法结果比较
Table 3. Results Comparison when Normal Equation Matrix is non Ill-posed
比较项 真值 最小二乘算法 方案1 方案2 方案3 ${\mathit{\boldsymbol{\hat X}}}$ -0.523 0 -0.516 0 -0.516 0 -0.516 0 -0.522 3 -2.758 0 -2.786 0 -2.786 0 -2.786 0 -2.801 2 0.902 0 0.923 3 0.923 3 0.923 3 0.883 1 -0.536 0 -0.560 6 -0.560 6 -0.560 6 -0.554 6 -1.560 0 -1.577 4 -1.577 4 -1.577 4 -1.598 7 2.503 0 2.456 2 2.456 2 2.456 2 2.502 4 2.334 0 2.327 0 2.327 0 2.327 0 2.336 9 -2.665 0 -2.728 6 -2.728 6 -2.728 6 -2.665 0 $F\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.008 5 0.004 0 0.004 0 0.004 0 0.005 1 m 0 0.008 5 0.008 5 0.008 5 0.004 1
下载: 导出CSV
表 4 算例2边长观测值
Table 4. Edge Length Observation Values of Example 2
边号 观测边长/m 1 45.075 2 5 222.056 3 5 187.391 4 7 838.867 5 5 483.162 6 5 731.756 7 5 438.383 8 7 493.316 9 8 884.603 10 8 839.687
下载: 导出CSV
表 5 法矩阵病态时算法的结果比较
Table 5. Results Comparison when Normal Equation Matrix is Ill-posed
比较项 真值 最小二乘方法 截断奇异值法 岭估计法
(L曲线法)方案1 方案2 方案3 ${\mathit{\boldsymbol{\hat X}}}$ -0.523 0 -1.347 3 -0.535 0 -0.552 7 -1.347 3 -0.510 5 -0.489 3 -2.758 0 6.162 5 -2.372 9 -2.263 9 6.162 5 -2.630 3 -2.853 2 0.902 0 10.443 9 1.358 5 1.382 6 10.443 9 1.084 4 0.847 1 -0.536 0 -0.364 5 -0.530 7 -0.497 5 -0.364 5 -0.535 6 -0.540 0 -1.560 0 2.869 2 -1.331 3 -1.328 4 2.869 2 -1.458 1 -1.567 8 2.503 0 -5.908 4 2.067 1 2.080 8 -5.908 4 2.307 9 2.516 2 2.334 0 -5.331 9 1.907 1 1.935 5 -5.331 9 2.125 7 2.314 7 -2.665 0 -16.214 1 -3.387 3 -3.328 9 -16.214 1 -3.000 0 -2.665 0 $F\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.004 3 1.686 1×10-5 0.001 2 0.004 4 1.686 1×10-5 0.001 2 0.001 3 m 0 504.044 1 1.303 2 1.308 9 504.044 1 0.253 8 0.013 8
下载: 导出CSV
-
[1] 朱建军, 丁晓利, 陈永奇.集成地质、力学信息和监测数据的滑坡动态模型[J].测绘学报, 2003, 32(3):261-266 doi: 10.3321/j.issn:1001-1595.2003.03.015 Zhu Jianjun, Ding Xiaoli, Chen Yongqi. Dynamic Landsliding Model with Integration of Monitoring Information and Mechanic Information[J]. Acta Geodaetica et Cartographica Sinica, 2003, 32(3):261-266 doi: 10.3321/j.issn:1001-1595.2003.03.015 [2] 张勤, 黄观文, 王利, 等.附有系统参数和附加约束条件的GPS城市沉降监测网数据处理方法研究[J].武汉大学学报·信息科学版, 2009, 34(3):269-272 http://ch.whu.edu.cn/CN/abstract/abstract1199.shtml Zhang Qin, Huang Guanwen, Wang Li, et al. Datum Design Study of GPS Height Monitoring Network with Systematic Parameters and Constraints[J]. Geomatics and Information Science of Wuhan University, 2009, 34(3):269-272 http://ch.whu.edu.cn/CN/abstract/abstract1199.shtml [3] 赵少荣, 陶本藻, 于正林.论变形测量数据的反演[J].测绘学报, 1992, 21(3):161-172 doi: 10.3321/j.issn:1001-1595.1992.03.001 Zhao Shaorong, Tao Benzao, Yu Zhenglin. On the Inversion of Deformation Survey Data[J]. Acta Geodaetica et Cartographica Sinica, 1992, 21(3):161-172 doi: 10.3321/j.issn:1001-1595.1992.03.001 [4] Mohammad B. Deformation Monitoring Using Different Least Squares Adjustment Methods:A Simulated Study[J]. KSCE Journal of Civil Engineering, 2016, 20(2):855-862 doi: 10.1007/s12205-015-0454-5 [5] 王乐洋, 许才军, 汪建军.附有病态约束矩阵的等式约束反演问题研究[J].测绘学报, 2009, 38(5):397-401 doi: 10.3321/j.issn:1001-1595.2009.05.004 Wang Leyang, Xu Caijun, Wang Jianjun. Research on Equality Constraint Inversion with Ill-posed Constraint Matrix[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5):397-401 doi: 10.3321/j.issn:1001-1595.2009.05.004 [6] 谢建, 朱建军.等式约束病态模型的正则化解及其统计性质[J].武汉大学学报·信息科学版, 2013, 38(12):1440-1444 http://ch.whu.edu.cn/CN/abstract/abstract2826.shtml Xie Jian, Zhu Jianjun. A Regularized Solution and Statistical Properties of Ill-Posed Problem with Equality Constraints[J]. Geomatics and Information Science of Wuhan University, 2013, 38(12):1440-1444 http://ch.whu.edu.cn/CN/abstract/abstract2826.shtml [7] 谢建, 朱建军.等式约束对病态问题的影响及约束正则化方法[J].武汉大学学报·信息科学版, 2015, 40(10):1344-1348 http://ch.whu.edu.cn/CN/abstract/abstract3345.shtml Xie Jian, Zhu Jianjun. Influence of Equality Constraints on Ill-conditioned Problems and Constrained Regularization Method[J]. Geomatics and Information Science of Wuhan University, 2015, 40(10):1344-1348 http://ch.whu.edu.cn/CN/abstract/abstract3345.shtml [8] 朱建军, 谢建.附不等式约束平差的一种简单迭代算法[J].测绘学报, 2011, 40(2):209-212 http://d.old.wanfangdata.com.cn/Periodical/chxb201102013 Zhu Jianjun, Xie Jian. A Simple Iterative Algorithm for Inequality Constrained Adjustment[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(2):209-212 http://d.old.wanfangdata.com.cn/Periodical/chxb201102013 [9] 宋迎春, 左廷英, 朱建军.带有线性不等式约束平差模型的算法研究[J].测绘学报, 2008, 37(4):433-437 doi: 10.3321/j.issn:1001-1595.2008.04.006 Song Yingchun, Zuo Tingying, Zhu Jianjun. Research on Algorithm of Adjustment Model with Linear Inequality Constrained Parameters[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(4):433-437 doi: 10.3321/j.issn:1001-1595.2008.04.006 [10] 宋迎春, 谢雪梅, 陈晓林.不确定性平差模型的平差准则与解算方法[J].测绘学报, 2015, 44(2):135-141 http://d.old.wanfangdata.com.cn/Periodical/chxb2015020005 Song Yingchun, Xie Xuemei, Chen Xiaolin. Adjustment Criterion and Algorithm in Adjustment Model with Uncertainty[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(2):135-141 http://d.old.wanfangdata.com.cn/Periodical/chxb2015020005 [11] 曾安敏, 杨元喜, 欧阳桂崇.附加约束条件的序贯平差[J].武汉大学学报·信息科学版, 2008, 33(2):183-186 http://ch.whu.edu.cn/CN/abstract/abstract1551.shtml Zeng Anmin, Yang Yuanxi, Ouyang Guichong. Sequential Adjustment with Constraints Among Parameters[J]. Geomatics and Information Science of Wuhan University, 2008, 33(2):183-186 http://ch.whu.edu.cn/CN/abstract/abstract1551.shtml [12] 刘根友, 欧吉坤.具有坐标函数约束的动态定位算法[J].武汉大学学报·信息科学版, 2004, 29(5):389-393 http://ch.whu.edu.cn/CN/abstract/abstract4601.shtml Liu Genyou, Ou Jikun. Kinematic Positioning Algorithm with Coordinate Function Constraint[J]. Geomatics and Information Science of Wuhan University, 2004, 29(5):389-393 http://ch.whu.edu.cn/CN/abstract/abstract4601.shtml [13] Peng J, Zhang H, Shong S, et al. An Aggregate Constraint Method for Inequality-Constrained Least Square Problem[J].Journal of Geodesy, 2006, 79(12):705-713 doi: 10.1007/s00190-006-0026-z [14] Song Yingchun, Zhu Jianjun, Li Zhiwei. The Least-Squares Estimation of Adjustment Model Constrained by Some Non-Negative Parameters[J]. Survey Review, 2010, 42(315):62-71 doi: 10.1179/003962610X12572516251367 [15] 周斌, 高立, 戴彧虹.求解大规模带边界约束二次规划问题的单调投影梯度法[J].中国科学A辑:数学, 2006, 36(5):556-570 http://d.old.wanfangdata.com.cn/Periodical/zgkx-ca200605005 Zhou Bin, Gao Li, Dai Yuhong. A Monotonic Gradient Projection Method in Large-Scale Quadratic Programming with Boundary Constraints[J]. Science in China Series A:Mathematics, 2006, 36(5):556-570 http://d.old.wanfangdata.com.cn/Periodical/zgkx-ca200605005 [16] Gao Yuelin, Xu Chengxian, Yang Chuansheng. A Global Optimality Method for Solving the Nonconvex Quadratic Programming Problem with Additional Box Constraints[J]. Journal of Engineering Mathematics, 2002, 19(1):99-102 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gcsxxb200201013 [17] 付文龙, 杜廷松, 翟军臣.基于D.C.分解的一类箱型约束的非凸二次规划的新型分支定界算法[J].数学研究, 2013, 46(3):311-318 doi: 10.3969/j.issn.1006-6837.2013.03.013 Fu Wenlong, Du Tingsong, Zhai Junchen. A New Branch and Bound Algorithm Based on D.C. Decomposition About Nonconvex Quadratic Programming with Box Constrained[J]. Journal of Mathematical Study, 2013, 46(3):311-318 doi: 10.3969/j.issn.1006-6837.2013.03.013 [18] 朱克强, 贺力群.大规模简单界约束的凸二次规划新算法[J].北方交通大学学报, 1998, 22(3):98-103 http://d.old.wanfangdata.com.cn/Periodical/bfjtdxxb199803022 Zhu Keqiang, He Liqun. New Algorithm of Large Scale Convex Quadratic Programs with Simple Bound Constraints[J]. Journal of Northern Jiaotong University, 1998, 22(3):98-103 http://d.old.wanfangdata.com.cn/Periodical/bfjtdxxb199803022 [19] Hestenes M R, Stiefel E. Methods of Conjugate Gradients for Solving Linear Systems[J]. Journal of Research of the National Bureau of Standards, 1952, 49(6):409-436 doi: 10.6028/jres.049.044 [20] Murty B S N, Husain A. Orthogonality Correction in the Conjugate-Gradient Method[J]. J Comput & Appl Math, 1983, 9:299-304 http://www.sciencedirect.com/science/article/pii/0377042783900018 [21] 肖运海.求解大规模优化问题的几种方法[D].长沙: 湖南大学, 2007 http://d.wanfangdata.com.cn/Thesis/Y1209515 Xiao Yunhai. Research on the Iteration Method for Several Kinds of Consistent and Inconsistent Constrained Matrix Equation[D]. Changsha: Hunan University, 2007 http://d.wanfangdata.com.cn/Thesis/Y1209515 [22] Rao C R. Linear Statistical Inference and Its Applications[M]. New York:John Wiley and Son, 1973 [23] 谢雪梅, 宋迎春, 肖兆兵.参数带有区间约束的平差模型迭代算法[J].测绘学报, 2018, 47(8):1141-1147 http://d.old.wanfangdata.com.cn/Periodical/chxb201808015 Xie Xuemei, Song Yingchun, Xiao Zhaobing. Parameter Estimate Algorithm in Adjustment Model with Interval Constraint[J]. Acta Geodaetica et Cartographica Sinica, 2018, 47(8):1141-1147 http://d.old.wanfangdata.com.cn/Periodical/chxb201808015 -
点击查看大图
图(1) / 表(5)
计量
- 文章访问数: 1203
- HTML全文浏览量: 102
- PDF下载量: 93
- 被引次数: 0
