在大地测量数据处理中，参数通常会有一些附加信息或者先验信息，这些信息能够在一定程度上给出参数解的可行范围，充分利用先验信息对参数建立某种约束，能够有效降低大地测量问题中的不适定性，保证解的稳定性和唯一性。一些学者利用先验信息建立了附等式约束和附不等式约束的平差模型^[1-3]，如朱建军等^[2]将罚函数与测量平差中的零权和无限权思想结合，提出一种简单的附不等式约束的平差迭代方法；王乐洋等^[3]提出附不等式约束的大地测量反演理论，有效解决了反演解的非唯一性问题。但这些前期得到的关于参数解的附加信息或者先验信息有时需要以区间、集合或模糊数形式表达^[4]，如果仍用不等式来表达这种区间约束，就会大大增加约束方程的个数，从而增加不等式约束解算的难度。对传统的平差模型附加这种区间约束，可以得到一种参数带区间约束的平差模型^[5]，这种平差模型可以转换成边界约束的二次规划问题。在数学领域，有一些学者提出解决这种边界约束二次规划的算法，如内点算法^[6]、分解方法^[7]、子空间截断牛顿法^[8]，但这些算法只注重解的最优性，忽视了平差的另一任务——精度评定。本文尝试将最优化理论中的子空间截断牛顿法引入到测量平差中，用病态问题对这种算法进行验证，并给出精度评估方法，扩展了现有的误差理论与测量平差方法。

1 参数带区间约束的平差模型与截断牛顿法

对传统的平差模型附加区间约束，可以建立参数带区间约束的平差模型：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} ,\mathit{\boldsymbol{ \boldsymbol{\varDelta} }} \sim N\left( {0,{\sigma ^2}{\mathit{\boldsymbol{P}}^{ - 1}}} \right)\\ {l_i} \le {x_i} \le {u_i},i = 1,2, \cdots ,n \end{array} \right. $

(1)

其中，B是m×n(m≥n)维系数矩阵，L是m维观测向量，Δ是m维随机误差向量，l={l₁, …,l_n}是参数X的下界，u ={u₁, …,u_n}是参数X的上界。对模型(1)进行最小二乘平差，即求解如下二次规划问题：

$ \left\{ \begin{array}{l} \min F\left( \mathit{\boldsymbol{X}} \right) = {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{PV}} = {\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{BX}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{BX}}} \right)\\ {\rm{s}}.\;{\rm{t}}.\;{l_i} \le {x_i} \le {u_i},i = 1,2, \cdots ,n \end{array} \right. $

(2)

在目标函数F(X)=(L -BX)^TP(L -BX)=X^T(B^TPB)X-2L^TPBX+L^TPL中，L^TPL是一个常量，因此二次规划问题(2)可以转化为如下边界约束的二次规划问题：

$ \left\{ \begin{array}{l} \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}}.\;{l_i} \le {x_i} \le {u_i},i = 1,2, \cdots ,n \end{array} \right. $

(3)

其中，N =B^TPB，c =B^TPL。g(X)=NX -c是f(X)在X处的梯度，H =N是f(X)在X处的Hesse矩阵^[8]。对于不带区间约束的二次规划minf(X)=$\frac{1}{2}$X^TNX -c^TX，修正牛顿法^[9]的一般迭代公式为：

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

(4)

式中，d_k是n维搜索方向，α_k是通过一维搜索满足f(X_k+α_kd_k)=min_α>0f(X_k+αd_k)得到的步长。

针对无约束优化的凸二次规划，牛顿法虽然具有局部二次收敛性，但在实际运用中存在不少问题，如Hesse矩阵求逆、选取初始值困难，修正牛顿法部分地克服了这些不足并具有非常强的整体收敛性。但牛顿方向主要具有局部性质，当迭代点远离最优解时求解牛顿方程的计算量非常大。为了减少计算量，Dembo等^[9]在修正牛顿法的基础上提出无约束优化的截断牛顿法(TNCG)，这种方法是利用共轭梯度法，按一种截断机制求出一个截断牛顿方向p。令g_p(X)和H_p分别是对应截断牛顿方向p的梯度和Hesse矩阵。截断牛顿法求截断牛顿方向p的具体步骤^[8-9]如下：

1) 令t=0，p₀=0，r₀=-g_p(X)，d₀=r₀以及δ₀=r₀^Tr₀；

2) 令q_t=H_pd_t，如果d_t^Tq_t≤εδ_t，停止迭代得到截断牛顿方向p：如果t=0，p =d₀，否则p =p_t，其中ε是充分小的常数，一般取ε=10^-6；

3) 令α_t=r_t^Tr_t/d_t^Tq_t，p_t+1=p_t+α_td_t以及r_t+1=r_t-α_tq_t，如果‖r_t+1‖/‖g_p(X)‖≤η，停止迭代并得到p =p_t+1；

4) 令β_t=r_t+1^Tr_t+1/r_t^Tr_t，d_t+1=r_t+1+β_td_t，δ_t+1=r_t+1^Tr_t+1+β_t²δ_t以及t=t+1，并转到步骤2)。

上述算法是截断牛顿法的内部迭代过程，目的是求出截断牛顿方向p，并将p作为搜索方向，其余外部迭代步骤与修正牛顿法相同。文献[15]证明了无约束的截断牛顿法的收敛性在线性收敛和二次收敛性之间，这种收敛性依赖于强制序列{η_k}的选取，一般具有超线性收敛性或更快的收敛性。强制序列{η_k}一般由用户给定，也可以由算法外部迭代自动产生：η_k=min{1/k, ‖g(X)‖},k=1, 2, …。但当参数受到区间约束时，搜索方向就不能完全由截断牛顿方向决定。针对这个问题，子空间截断牛顿法将搜索方向分成4个部分，截断牛顿方向只是搜索方向的一部分。

2 基于积极集的子空间截断牛顿法

定义有效指标集I ={i:l_i≤x_i≤l_i+b}∪{i:u_i-b≤x_i≤u_i}，非有效指标集J ={i:l_i+b≤x_i≤u_i-b}，将指标在I中的参数称为有效参数，将指标在J中的参数称为非有效参数，其中b是很小的常数，可取b=10^-6，再将指标集I分成3个部分：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{I}}_1} = \left\{ {i:{x_i} = {l_i}} \right\} \cup \left\{ {i:{x_i} = {u_i}} \right\} \cap \left\{ {i:\left( {{l_i} + {u_i} - 2{x_i}} \right){g_i}\left( x \right) \ge 0} \right\}\\ {\mathit{\boldsymbol{I}}_2} = \left\{ {i:{l_i} \le {x_i} \le {l_i} + b} \right\} \cup \left\{ {i:{u_i} - b \le {x_i} \le {u_i}} \right\} \cap \left\{ {i:\left( {{l_i} + {u_i} - 2{x_i}} \right){g_i}\left( x \right) < 0} \right\}\\ {\mathit{\boldsymbol{I}}_3} = \left\{ {i:{l_i} < {x_i} \le {l_i} + b} \right\} \cup \left\{ {i:{u_i} - b \le {x_i} < {u_i}} \right\} \cap \left\{ {i:\left( {{l_i} + {u_i} - 2{x_i}} \right){g_i}\left( x \right) \ge 0} \right\} \end{array} \right. $

(5)

其中，指标在I₁中的有效变量是最速下降方向指向可行域外侧的，指标在I₂中的有效变量是最速下降方向指向可行域内部的，指标在I₃中的有效变量是最速下降方向指向可行域边界的^[8]。

子空间截断牛顿法的关键在于求解搜索方向，其搜索方向根据有效集和非有效集的划分被分为4个部分。

第一步，求解非有效变量的搜索方向。这一部分的搜索方向由截断牛顿方向决定，类似于无约束优化，得到的截断牛顿方向p需要通过转化得到非有效参数的搜索方向。具体方法是令p左乘n×n_J维Z矩阵，即d =Zp，此时对应有效指标集的搜索方向为0，n_J代表指标集J中元素的个数，z_sj是Z矩阵中第s行第j列元素，由式(6)得到：

$ {z_{sj}} = \left\{ \begin{array}{l} 1,\;\;如果\;s = {J_j}\\ 0,\;\;其他 \end{array} \right. $

(6)

其中，J_j代表指标集J中第j个元素。

第二步，求解有效变量的搜索方向。对应I₁、I₂、I₃这3个有效指标集的搜索方向，由式(7)得到：

$ {\mathit{\boldsymbol{d}}_i} = \left\{ \begin{array}{l} 0,i \in {\mathit{\boldsymbol{I}}_1}\\ - {g_i}\left( \mathit{\boldsymbol{X}} \right),i \in {\mathit{\boldsymbol{I}}_2}\\ - {\lambda _i}{g_i}\left( \mathit{\boldsymbol{X}} \right),i \in {\mathit{\boldsymbol{I}}_3} \end{array} \right. $

(7)

其中，${\lambda _i} = \left\{ \begin{array}{l} \left( {{x_i}-{l_i}} \right)/{g_i}\left( \mathit{\boldsymbol{X}} \right), {\mathit{l}_i} < {x_i} \le {l_i} + b\\ \;\;\;\;\;\;\;\;\;\;\;{\rm{且}}{\mathit{x}_i}-{g_i}\left( \mathit{\boldsymbol{X}} \right) \le {l_i}\\ \left( {{x_i}-{u_i}} \right)/{g_i}\left( \mathit{\boldsymbol{X}} \right), {\mathit{u}_i} - b < {x_i} \le {u_i}\\ \;\;\;\;\;\;\;\;\;\;\;\;{\rm{且}}{\mathit{x}_i} - {g_i}\left( \mathit{\boldsymbol{X}} \right) \ge {u_i}\\ 1, {\rm{其他}} \end{array} \right.$

子空间截断牛顿法的步长α通过投影搜索法得到，能够保证参数在约束区间范围内，即满足f(P_Ω[X +αd])≤ f(P_Ω[X])+μαg(X)^Td，其中μ∈(0, 0.5)是给定的常数，P_Ω是Ω ={X :l_i≤x_i≤u_i,i=1, 2, …,n}的投影算子，投影公式如(8)，满足要求的步长α通过Armijio非精确搜索技巧确定：

$ {\rm{mid}}\left( {\mathit{\boldsymbol{X}},\mathit{\boldsymbol{l}},\mathit{\boldsymbol{u}}} \right) = \left\{ \begin{array}{l} {l_i},{x_i} < {l_i}\\ {x_i},{x_i} \in \left[ {{l_i},{u_i}} \right]\\ {u_i},{x_i} > {u_i} \end{array} \right. $

(8)

子空间截断牛顿算法步骤如下：

1) 选取初始点X₀，使得l <X₀ <u，计算g(X₀)和H，令k=0；

2) 划分指标集I_1k、I_2k、I_3k、J_k；

3) 构建Z_k，g_p(X_k)=Z_k^Tg(X_k)，H_pk=Z_kH，η_k=min{1/k, ‖g(X)‖}，用TNCG算法解算得到截断牛顿方向p_k，令d_k=Z_kp_k，继续计算式(7)得到完整的搜索方向d_k；

4) 选取参数μ∈(0, 0.5)且令α=1，求最小非负整数w(w=0, 1, …)，从w=0开始取值，直到满足：

$ \begin{array}{*{20}{c}} {f\left( {{\mathit{\boldsymbol{P}}_\mathit{\Omega }}\left[ {{\mathit{\boldsymbol{X}}_k} + {2^{ - w}}\bar \alpha {\mathit{\boldsymbol{d}}_k}} \right]} \right) \le f\left( {{\mathit{\boldsymbol{P}}_\mathit{\Omega }}\left[ {{\mathit{\boldsymbol{X}}_k}} \right]} \right) + }\\ {\mu {2^{ - w}}\bar \alpha g{{\left( {{\mathit{\boldsymbol{P}}_\mathit{\Omega }}\left[ {{\mathit{\boldsymbol{X}}_k}} \right]} \right)}^{\rm{T}}}{\mathit{\boldsymbol{d}}_k}} \end{array} $

令α_k=2^-wα，X_k+1=P_Ω[X_k+α_kd_k]，计算g(X_k+1)；

5) 如果满足‖X_k+1-X_k‖<$\epsilon$，取$\epsilon$=10^-10，停止迭代，最终结果为X_k+1；否则令k=k+1，转步骤2)。

当搜索方向d_k≠0时，这始终是一个下降方向，即f(X_k)是严格单调递减的，而且集合I₂、I₃在大部分迭代中是空的。梁昔明等^[8]证明算法得到的序列{X_k}的任一聚点均是问题(3)的KT点，即证明了子空间截断牛顿法具有整体收敛性。

针对边界约束的二次规划问题(3)，子空间截断牛顿法利用有效集的思想将参数划分成有效参数和非有效参数进行迭代，最终参数得到有效约束，从而参数带区间约束的平差模型(1)可以转化为带等式约束的平差模型：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }},\mathit{\boldsymbol{ \boldsymbol{\varDelta} }} \sim N\left( {0,{\sigma ^2}{\mathit{\boldsymbol{P}}^{ - 1}}} \right)\\ \mathit{\boldsymbol{CX}} = {\mathit{\boldsymbol{L}}_x} \end{array} \right. $

(9)

式中，若有效参数个数为m_x，则矩阵C是m_x×n维约束矩阵，其中对应有效参数x_i的元素为1，其他元素为0；L_x是m_x维约束向量，元素是对应的l_i或u_i。

平差模型(9)的参数估值与参数协因数矩阵^[10]分别为：

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\hat X}}}_C} = \left( {\mathit{\boldsymbol{N}}_b^{ - 1} - \mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{N}}_c^{ - 1}\mathit{\boldsymbol{CN}}_b^{ - 1}} \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PL}} + }\\ {\mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{N}}_c^{ - 1}{\mathit{\boldsymbol{L}}_x}} \end{array} $

(10)

$ {\mathit{\boldsymbol{Q}}_X} = \mathit{\boldsymbol{N}}_b^{ - 1} - \mathit{\boldsymbol{N}}_b^{ - 1}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{N}}_c^{ - 1}\mathit{\boldsymbol{CN}}_b^{ - 1} $

(11)

其中，N_b=B^TPB,N_c=CN_b^-1C^T。平差模型(9)与参数带区间约束的平差模型(1)是等价的，因此参数协因数矩阵(11)也为平差模型(1)的参数协因数矩阵。

3 算例分析

算例中测边网数据引自文献[11]。该测边网得到的系数矩阵是病态矩阵，P₁点到P₁₀点是已知点，其坐标以及与未知点之间的观测边长如表 1所示；K₁、K₂、K₃为未知点，其真实坐标与近似坐标如表 2所示，各观测量为等精度观测。根据前期观测得到的点位精度为σ=±0.01 m，约束区间D₁是由点位精度推算的，为了便于算法分析，将约束区间D₁放大2倍、3倍得到D₂和D₃。3个未知点之间的观测边长分别为d_K1K2=88.340 2 m，d_K1K3=73.355 1 m。相对于未知点近似坐标的改正数参数估值为$\widehat X$ ={$\widehat x$₁, $\widehat y$₁, $\widehat z$₁, $\widehat x$₂, $\widehat y$₂, $\widehat z$₂, $\widehat x$₃, $\widehat y$₃, $\widehat z$₃}^T，测边网如图 1所示。

方案1：令θ ={0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01}^T，l₁=-3θ，u₁=3θ，约束区间为D₁={X :l₁≤X ≤u₁}；方案2：l₂=-6θ，u₂=6θ，约束区间为D₂={X:l₂≤X ≤u₂}；方案3：l₃=-9θ，u₃=9θ，约束区间为D₃={X:l₃≤X ≤u₃}。

表 3是求解测边网未知点坐标改正数的算法结果。从中可以看出：

1) 经典最小二乘法的结果偏离真值较大，m=32.367 5，严重失真。

2) 岭估计(L曲线)^[12]的参数估值与真值的残差平方和m=0.013 4，对于由前期点位精度推算的约束区间D₁，本文算法的m=0.004 4 < 0.013 4，这说明本文算法参数估值比岭估计(L曲线)的结果好；从另一方面来看，若D₁是一个可靠的先验信息，最小二乘解和岭估计解并不是D₁的可行解。

3) 对于放大2倍和3倍的约束区间，本文算法虽然没有比岭估计(L曲线)的结果好，但m也远小于经典最小二乘法，当先验信息有效时，本文算法能够改善病态性。

4) 由约束区间D₁下的算法结果可得：

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} 0&1&0&0&0&0&0&0&0\\ 0&0&1&0&0&0&0&0&0\\ 0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&0&1 \end{array}} \right]\\ {\mathit{\boldsymbol{L}}_x} = \left[ {\begin{array}{*{20}{c}} { - 0.0300}\\ {0.0300}\\ { - 0.0300}\\ {0.0300} \end{array}} \right] \end{array} \right. $

(12)

根据式(10)，解得：$\mathit{\boldsymbol{\widehat X}}$_C=[-0.020 5, -0.030 0, 0.030 0, 0.003 6, 0.005 3, -0.002 4, -0.030 0, -0.007 1, 0.030 0]^T, 与本文算法的结果完全一致，则由式(11)可以解得平差算法的参数协因数阵Q_X。

5) 测边网的算法结果进一步说明，本文算法在充分利用先验信息的基础上，能够有效改善平差模型的病态性，并能给出精度评定。

6) 随着先验约束区间的增大，算法的精度会逐渐降低，对病态性的改善程度也会逐渐减弱。

4 结语

大地测量数据处理过程中，利用可靠的先验信息对参数附加一定的约束，能够降低平差模型的不适定性。针对参数带有区间约束的平差模型，本文基于积极集的思想，提供了一种关于区间约束解算的新算法。子空间截断牛顿法可以迅速改变积极约束集的构成，得到区间约束的积极集(也称为有效约束)，它不需要化为不等式约束而增加算法的复杂性，可以直接对区间约束进行迭代解算。在精度评估上，可以利用积极集构成的等式约束，通过带有等式约束的平差理论，对参数的精度进行评估，解决了不等式约束平差算法的精度评估问题，通过最优估计理论拓展了现有的误差理论与测量平差方法。