-
摘要: 不等式约束部分变量含误差(partial errors-in-variables, PEIV)模型目前主要采用线性化方法和非线性规划类算法, 前者计算效率较低, 后者基于最优化理论, 计算复杂, 未能与经典平差理论建立联系, 难以在测量实际中推广。在整体最小二乘准则下, 根据最优解的Kuhn-Tucker条件, 将不等式约束整体最小二乘解的计算转化为二次规划问题, 并提出改进的Jacobian迭代法求解二次规划。所提方法不需要对观测方程线性化, 与经典最小二乘法具有相同的形式, 易于编程实现。数值实例表明, 所提方法形式简洁, 具有良好的计算效率, 是经典最小二乘平差理论的有益拓展。
-
关键词:
- 不等式约束 /
- 部分变量含误差模型 /
- 整体最小二乘 /
- Jacobian迭代法 /
- 经典最小二乘
Abstract:Objectives The inequality constrained partial errors-in-variables(ICPEIV) model is mainly solved by linear approximation method and nonlinear programming algorithms. The linear approximation method is computationally inefficient and the nonlinear programming algorithms are complicated because they are based on optimization theory. The nonlinear programming algorithms are impracticable to apply in surveying fields because the connections between nonlinear programming methods and classical adjustment have not been established.Methods Under the total least squares(TLS) criterion, the inequality constrained TLS problem is transformed into the quadratic programming according to the Kuhn-Tucker condition. Then an improved Jacobian iteration approach is proposed to solve the quadratic programming.Results The proposed method does not require the linearization process and has the same form with the classical least squares which is easy to code.Conclusions The numerical examples show that the proposed method is efficient in computation and concise in form and it is a beneficial extension of classical least squares theory. -
在大地测量坐标转换、曲线曲面拟合、摄影测量后方交会、大地测量反演等领域[1-6]，平差系统的系数矩阵和观测向量均含有随机误差，这类模型称为变量含误差（errors-in-variables, EIV）模型。令随机误差的平方和最小，求EIV模型最优解的算法称为整体最小二乘（total least squares, TLS）方法。当随机误差都服从独立同分布时，称为等权TLS。当误差相关、精度不相等时，称为加权TLS(weighted TLS, WTLS）。Fang[7]推导了WTLS问题获得最优解的条件，并提出了3类典型算法。Mahboub[8]给出了构建系数阵权阵的几点规则和特殊结构下WTLS的迭代方案，并运用到线性回归和二维仿射变换。Shen等[9]将WTLS问题转换为非线性加权最小二乘问题，利用Gauss-Newton法得到的解在形式上与经典最小二乘解（least squares, LS）相同，并采用蒙特卡罗模拟得到了单位权方差的近似无偏估计。Amiri-Simkooei等[10]也给出了与经典LS解形式相同的WTLS解。
当参数间存在可靠的先验信息时，将其纳入EIV模型能提高参数估计的精度和可靠性。Schaffrin等[11]将含线性约束的TLS问题转化为特征值问题，得到了方差的近似值，并扩展到同时含有线性和二次约束的情形。Zhang等[12]提出了解不等式约束WTLS(inequality constrained WTLS, ICWTLS）问题的穷举搜索法，其计算量随约束的增加迅速增大。Fang[13]采用起作用集方法和序列二次规划（sequential quadratic programming, SQP）算法求解ICWTLS问题。
Xu等[14]将EIV模型扩展成部分变量含误差（partial errors-in-variables, PEIV）模型并推导了其算法和近似方差。Shi等[15]给出了PEIV模型的经典LS形式的解，并分析了算法计算量与系数阵中含误差元素个数间的关系。Zeng等[16]提出了不等式约束PEIV(inequality constrained PEIV, ICPEIV）模型的线性近似方法，将其转化为二次规划子问题，其计算效率较低。ICPEIV模型在WTLS准则下可归结为约束非线性规划问题，采用最优化理论的乘子法、惩罚函数法、积极约束法等，计算过程较复杂，与经典平差理论相去甚远[17]。解不等式约束最小二乘（inequality constrained LS, ICLS）问题的Jacobian迭代算法与解等式约束LS问题的算法相似，ICLS问题可纳入经典平差模型的范畴。
为了简化ICPEIV模型的计算，得到与经典平差形式一致的算法，本文提出了解ICPEIV模型的经典LS方法。首先，根据ICPEIV模型获得最优解的一阶必要条件，将模型参数和系数阵元素的估计进行分离。模型参数的估计归结为解ICLS问题，系数阵元素的估计根据两类参数的关系得到。然后，提出扩展的Lagrange方法，将ICLS问题归结为经典的等式约束平差问题，不需要对观测方程线性化，也无需利用规划类算法，与经典LS方法形式一致。
1 ICPEIV模型及WTLS解
1.1 ICPEIV模型的最优性条件
ICPEIV模型的函数模型可以表示为[18]:
$$ \mathit{\boldsymbol{y}} = \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) + {\mathit{\boldsymbol{e}}_y} $$ (1a) $$ \mathit{\boldsymbol{a}} = \mathit{\boldsymbol{\bar a}} + {\mathit{\boldsymbol{e}}_a} $$ (1b) $$ \mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0 $$ (1c) 式中，y和ey分别表示n维观测向量及其误差向量；β表示m维模型参数；In是n维单位矩阵；⊗是Kronecker积符号，其定义为A⊗B=[aij×B], aij是矩阵A的第i行第j列元素；a是已知t维列向量；a和ea分别是a的真值和误差向量；h是已知的nm维常数向量；B是nm×t维常数矩阵；a、B的定义见文献[15]; G是s×m维约束矩阵；d是s维常数向量。令随机误差e=(eaT eyT)T，其均值为0，且ea和ey互不相关，它们的方差分别为：
$$ D\left( {{\mathit{\boldsymbol{e}}_y}} \right) = \sigma _0^2{\mathit{\boldsymbol{W}}^{ - 1}} $$ (1d) $$ D\left(\boldsymbol{e}_{a}\right)=\sigma_{0}^{2} \boldsymbol{\omega}^{-1} $$ (1e) 式中，σ02为单位权方差；W为n维观测值权矩阵；ω为t维随机误差权矩阵。
ICPEIV模型的整体最小二乘解等价于如下约束非线性规划问题的最优解：
$$ \left\{ {\begin{array}{*{20}{l}} {\min :\mathit{\Phi }(\mathit{\boldsymbol{\bar a}},\mathit{\boldsymbol{\beta }}) = {{(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}})}^{\rm{T}}}\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}) + }\\ {\quad {{\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)}^{\rm{T}}} \times }\\ {\quad \mathit{\boldsymbol{W}}\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\quad \mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0} \end{array}} \right. $$ (2) 由最优化理论得到Lagrange目标函数为：
$$ \begin{array}{l} \psi (\mathit{\boldsymbol{\bar a}},\mathit{\boldsymbol{\beta }}) = {\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right)^{\rm{T}}} \times \\ \quad \mathit{\boldsymbol{W}}\left( {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\bar a}}) - \mathit{\boldsymbol{y}}} \right) + \\ \quad {(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}})^{\rm{T}}}\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}) - {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}(\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} - \mathit{\boldsymbol{\eta }}) \end{array} $$ (3) 式中，λ是s维Lagrange乘子向量；η是非负基变量。将$ \psi(\bar{{\boldsymbol{a}}}, {\boldsymbol{\beta}}) $分别对$\bar{{\boldsymbol{a}}} $、β、λ求偏导，令其在最优解$ \hat{\bar{{\boldsymbol{a}}}} \text { 和 } \hat{{\boldsymbol{\beta}}} $处的值为0，可得：
$$ \begin{array}{c} \frac{{\partial \psi (\mathit{\boldsymbol{\bar a}},\mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\bar a}}}} = 2\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{\hat {\bar a}}} - \mathit{\boldsymbol{a}}) + 2{\mathit{\boldsymbol{B}}^{\rm{T}}}\left( {\mathit{\boldsymbol{\hat \beta }} \otimes {\mathit{\boldsymbol{I}}_n}} \right) \times \\ \mathit{\boldsymbol{W}}\left( {\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\hat {\bar a}}}) - \mathit{\boldsymbol{y}}} \right) = 0 \end{array} $$ (4a) $$ \begin{array}{l} \frac{{\partial \psi (\mathit{\boldsymbol{\bar a}},\mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\beta }}}} = 2\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{h}}_1^{\rm{T}} + {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_1^{\rm{T}}}\\ {\mathit{\boldsymbol{h}}_2^{\rm{T}} + {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_2^{\rm{T}}}\\ \vdots \\ {\mathit{\boldsymbol{h}}_m^{\rm{T}} + {{\mathit{\boldsymbol{\hat {\bar a}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_m^{\rm{T}}} \end{array}} \right]\mathit{\boldsymbol{W}} \times \\ \left\{ {\left( {{{\mathit{\boldsymbol{\hat \beta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)(\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B\hat {\bar a}}}) - \mathit{\boldsymbol{y}}} \right\} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = 0 \end{array} $$ (4b) $$ \frac{{\partial \psi (\mathit{\boldsymbol{\bar a}},\mathit{\boldsymbol{\beta }})}}{{\partial \mathit{\boldsymbol{\lambda }}}} = \mathit{\boldsymbol{G\hat \beta }} - \mathit{\boldsymbol{d}} + \mathit{\boldsymbol{\eta }} = 0 $$ (4c) 式中，$ \boldsymbol{h}^{\mathrm{T}}=\left[\begin{array}{llll} \boldsymbol{h}_{1}^{\mathrm{T}} & \boldsymbol{h}_{2}^{\mathrm{T}} & \cdots & \boldsymbol{h}_{m}^{\mathrm{T}} \end{array}\right], {\boldsymbol{B}}^{\mathrm{T}}=\left[\begin{array}{llll} {\boldsymbol{B}}_{1}^{\mathrm{T}} & {\boldsymbol{B}}_{2}^{\mathrm{T}} & \cdots & {\boldsymbol{B}}_{m}^{\mathrm{T}} \end{array}\right] $, hi是n维列向量，Bi是n×t维矩阵(i=1, 2…m)。由式（4a）、（4b）、（4c）可得ICPEIV模型获得最优解的Kuhn-Tucker(KT）条件为：
$$ \begin{array}{l} \left( {\mathit{\boldsymbol{\omega }} + \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}}} \right)\mathit{\boldsymbol{\hat {\bar a}}} - \mathit{\boldsymbol{\omega a}} + \\ \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{W}}\left( {\sum\limits_{i = 1}^m {{\mathit{\boldsymbol{h}}_i}} {{\mathit{\boldsymbol{\hat \beta }}}_i}} \right) - \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{Wy}} = 0 \end{array} $$ (5a) $$ \begin{array}{*{20}{l}} {\left( {{\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{h}}} + {\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{B}}} + {\mathit{\boldsymbol{N}}_{\mathit{\boldsymbol{Bh}}}} + {\mathit{\boldsymbol{N}}_{\mathit{\boldsymbol{hB}}}}} \right)\mathit{\boldsymbol{\hat \beta }} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{h}}} - {\mathit{\boldsymbol{u}}_\mathit{\boldsymbol{B}}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = 0} \end{array} $$ (5b) $$ \mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\hat \beta }} - {d_i} \ge 0,i = 1,2 \cdots s $$ (5c) $$ {\lambda _i}\left( {\mathit{\boldsymbol{g}}_i^{\rm{T}}\mathit{\boldsymbol{\hat \beta }} - {d_i}} \right) = 0,i = 1,2 \cdots s $$ (5d) $$ {\lambda _i} \ge 0,i = 1,2 \cdots s $$ (5e) 式中，$ \boldsymbol{S}_{{\boldsymbol{\beta}}}=\sum\limits_{i=1}^{m} {\boldsymbol{B}}_{i} \hat{{\boldsymbol{\beta}}}_{i}=\left(\hat{{\boldsymbol{\beta}}}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{n}\right) {\boldsymbol{B}}, {\boldsymbol{N_{h}}}, {\boldsymbol{N_{B}}} 、{\boldsymbol{N_{B h}}}、\boldsymbol{N}_{{\boldsymbol{h B}}}, \boldsymbol{u}_{{\boldsymbol{h}}}, \boldsymbol{u}_{{\boldsymbol{B}}}$的含义见文献[15]。令GT= $ \left[\begin{array}{llll} {\boldsymbol{g}}_{1} & {\boldsymbol{g}}_{2} & \cdots & {\boldsymbol{g}} \end{array}\right], {\boldsymbol{g}}_{i}^{\mathrm{T}}=\left[\begin{array}{llll} g_{i 1} & g_{i 2} & \cdots & g_{i m} \end{array}\right] \text { 是 } {\boldsymbol{G}} $的第i行，常数项d=[d1 d2 … ds]T。若$ (\hat{\bar{{\boldsymbol{a}}}}, \hat{{\boldsymbol{\beta}}}) $是式（2）的局部最优解，则满足KT必要条件，但式（5a)~(5e）并未直接提供ICPEIV模型的最优化算法。本文通过已有的线性近似方法和SQP方法，提出可直接利用经典LS原理的算法。
1.2 线性近似方法求ICPEIV模型的WTLS解
给定模型参数和系数阵中随机元素的近似值$ \boldsymbol{\beta}_{a}^{0}=\left[\begin{array}{ll} {\boldsymbol{\beta}}_{0}^{\mathrm{T}} & \bar{{\boldsymbol{a}}}_{0}^{\mathrm{T}} \end{array}\right]^{\mathrm{T}}, \text { 令 } \Delta {\boldsymbol{\beta}}={\boldsymbol{\beta}}_{0}-{\boldsymbol{\beta}}, \Delta {\boldsymbol{a}}=\bar{{\boldsymbol{a}}}_{0}-\bar{{\boldsymbol{a}}}, \Delta {\boldsymbol{d}}={\boldsymbol{d}}- {\boldsymbol{G \beta}}^{0}$。忽略二阶及高阶项，则式（1a)~(1e）可以线性化为[16]:
$$ \left\{ {\begin{array}{*{20}{l}} {\Delta \mathit{\boldsymbol{y}} = {\mathit{\boldsymbol{A}}_0}\Delta \mathit{\boldsymbol{\beta }} + \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\mathit{\boldsymbol{B}}\Delta \mathit{\boldsymbol{a}} + {\mathit{\boldsymbol{e}}_y}}\\ {{{\mathit{\boldsymbol{\bar a}}}_0} - \mathit{\boldsymbol{a}} = \Delta \mathit{\boldsymbol{a}} + {\mathit{\boldsymbol{e}}_a}}\\ { - \mathit{\boldsymbol{G}}\Delta \mathit{\boldsymbol{\beta }} \ge \Delta \mathit{\boldsymbol{d}}} \end{array}} \right. $$ (6) 式中，$ \Delta {\boldsymbol{y}}={\boldsymbol{A}}_{0} {\boldsymbol{\beta}}_{0}- {\boldsymbol{y}}, {\boldsymbol{A}}_{0}=\operatorname{ivec}\left({\boldsymbol{h}}+{\boldsymbol{B}} \bar{{\boldsymbol{a}}}_{0}\right) $, ivec（）是将n×m维列向量转化为n×m维矩阵的操作算子。令$ \Delta {\boldsymbol{\beta}}_{a}=\left[\begin{array}{c} \Delta {\boldsymbol{\beta}} \\ \Delta {\boldsymbol{ a}} \end{array}\right], \Delta {\boldsymbol{y}}_{a}=\left[\begin{array}{c} \Delta {\boldsymbol{y}} \\ \bar{{\boldsymbol{a}}}_{0}- {\boldsymbol{a}} \end{array}\right], {\boldsymbol{G}}_{a}=-\left[\begin{array}{ll} {\boldsymbol{G}} & 0 \end{array}\right] $, $ \Delta {\boldsymbol{d}}_{a}={\boldsymbol{d}}- {\boldsymbol{G}}_{a} {\boldsymbol{\beta}}_{a}^{0}, {\boldsymbol{A}}_{a}=\left[\begin{array}{cc} {\boldsymbol{A}}_{0} & \left({\boldsymbol{\beta}}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{n}\right) {\boldsymbol{B}} \\ 0 & {\boldsymbol{I}}_{t} \end{array}\right] $，则式（6）可以写成：
$$ \left\{ {\begin{array}{*{20}{l}} {\Delta {\mathit{\boldsymbol{y}}_a} = {\mathit{\boldsymbol{A}}_a}\Delta {\mathit{\boldsymbol{\beta }}_a} + \mathit{\boldsymbol{e}}}\\ {{\mathit{\boldsymbol{G}}_a}\Delta {\mathit{\boldsymbol{\beta }}_a} \ge \Delta {\mathit{\boldsymbol{d}}_a}} \end{array}} \right. $$ (7) 式（7）是经典的ICLS模型，可以用线性互补算法、积极约束法等求解。线性化方法包含内外迭代，外迭代是由近似值进行线性化的过程，内迭代是组成ICLS模型、采用二次规划（quadratic programming, QP）算法得到新的近似值的过程。当近似值距离最优解较远时，每代入一次新的近似值，就要计算一个n+t维二次规划问题，在大规模问题中，计算量会显著增加。
1.3 SQP算法求ICPEIV模型的WTLS解
令$ \boldsymbol{u}^{\mathrm{T}}=\left[\begin{array}{ll} \bar{{\boldsymbol{a}}}^{\mathrm{T}} & \boldsymbol{\beta}^{\mathrm{T}} \end{array}\right], {\boldsymbol{G}}_{\vartheta}=\left[\begin{array}{ll} 0 & {\boldsymbol{G}} \end{array}\right] $，则式（2）可以表示成约束非线性规划问题：
$$ \left\{ {\begin{array}{*{20}{l}} {\min f(\mathit{\boldsymbol{u}})}\\ {{\rm{s}}{\rm{.t}}{\rm{. }}{\mathit{\boldsymbol{G}}_\vartheta }\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{d}} \ge 0} \end{array}} \right. $$ (8) 给定初始值u(j)，计算函数f(u)在u(j)处的梯度g(u(j))和Hessian矩阵H (u(j))。将目标函数和约束条件在u(j)处用泰勒级数展开，保留二次项，则式（8）转化为QP问题[13]:
$$ \left\{\begin{array}{l} \min f\left(\boldsymbol{u}^{(j)}+\mathrm{d} \boldsymbol{u}\right)=\frac{1}{2} \mathrm{~d} \boldsymbol{u}^{\mathrm{T}}\boldsymbol{ H}\left(\boldsymbol{u}^{(j)}\right) \mathrm{d} \boldsymbol{u}+ \\ \ \ \ \ \ \ \ \ \ \ \ \ \boldsymbol{g}^{\mathrm{T}}\left(\boldsymbol{u}^{(j)}\right) \mathrm{d} \boldsymbol{u}+f\left(\boldsymbol{u}^{(j)}\right) \\ \text { s.t. } \boldsymbol{G}_{\vartheta} \boldsymbol{u}^{(j)}-\boldsymbol{d}+\boldsymbol{G}_{\vartheta} \mathrm{d} \boldsymbol{u} \geqslant 0 \end{array}\right. $$ (9) 采用积极约束法求解，得到du的最优值。令新的近似值u(j+1)=u(j)+du(j=1, 2…N），形成新的QP问题。计算步骤为：
1) 令$ \boldsymbol{u}^{(0)}=\left[\begin{array}{ll} \bar{{\boldsymbol{a}}}^{(0) \mathrm{T}} & \boldsymbol{\beta}^{(0) \mathrm{T}} \end{array}\right]^{\mathrm{T}}, \text { 计算 } {\boldsymbol{g}}^{\mathrm{T}}\left(\boldsymbol{u}^{(0)}\right) \text { 和 } $ H(u(0))，在可行集内选取初值du0，利用积极约束法求解式（9）的QP问题；
2）若‖du‖ ≤ ε（ε是任意给定的小的正数），则停止迭代；若‖du‖ > ε，令u(j) = u(j - 1) + du，转步骤1）。
g(u(j))和H (u(j))的计算见文献[18]。SQP算法是一种牛顿型算法，当目标函数和约束方程在最优解的邻域内二次可微、且约束条件线性无关时，若一阶KT条件和二阶充分条件成立，则算法具有二阶收敛速率。SQP方法需要求解目标函数对参数的微分和Hessian矩阵，计算较为复杂。针对线性化和SQP算法的不足，本文提出一种无需线性化、与经典LS平差模型形式一致的迭代算法。
2 ICPEIV模型的经典LS算法
2.1 ICPEIV模型的变量分离算法
由式（5a）可知，若给定参数β的一个估值$ \hat {\boldsymbol{\beta}} $，可得到相应的估计$ \hat {\bar {\boldsymbol{a}}} $为[14]:
$$ \begin{array}{l} \mathit{\boldsymbol{\hat {\bar a}}} = {\left( {\mathit{\boldsymbol{\omega }} + \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}}} \right)^{ - 1}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\mathit{\boldsymbol{\omega a}} - \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{h}}_\mathit{\boldsymbol{\beta }}} + \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}\mathit{\boldsymbol{Wy}}} \right) \end{array} $$ (10) 式中，$ \boldsymbol{h}_{{\boldsymbol{\beta}}}=\sum\limits_{i=1}^{m} \boldsymbol{h}_{i} \hat{\boldsymbol{\beta}}_{i} $。当$ \hat {\bar {\boldsymbol{a}}} $的元素个数t小于观测值的个数n时，采用式（10）计算$ \hat {\bar {\boldsymbol{a}}} $值，只需求t维矩阵的逆，计算量较小。若t > n，用式（11）计算的$ \hat {\bar {\boldsymbol{a}}} $值与式（10）等价，却只需求n维矩阵的逆，此时计算量小[15]。
$$ \mathit{\boldsymbol{\hat {\bar a}}} = \mathit{\boldsymbol{a}} + {\mathit{\boldsymbol{\omega }}^{ - 1}}\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\beta }}^{\rm{T}}{\mathit{\boldsymbol{E}}^{ - 1}}(\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A\hat \beta }}) $$ (11) 式中，$ {\boldsymbol{E}}=\boldsymbol{W}^{-1}+\boldsymbol{S}_{\beta} \boldsymbol{\omega}^{-1} \boldsymbol{S}_{\beta}^{\mathrm{T}} \text { 令 } \bar{{\boldsymbol{A}}}=\operatorname{ivec}({\boldsymbol{h}}+{\boldsymbol{B}} \bar{{\boldsymbol{a}}}) , $$ \text { 则 有 }\left(\hat{{\boldsymbol{\beta}}}^{\mathrm{T}} \otimes {\boldsymbol{I}}_{n}\right)({\boldsymbol{h}}+ {\boldsymbol{B}} \bar{{\boldsymbol{a}}})=\bar{{\boldsymbol{A}}} \hat{{\boldsymbol{\beta}}} $。因此式（2）等价于求解约束最优化问题：
$$ \left\{ {\begin{array}{*{20}{l}} {\min \mathit{\Phi }(\mathit{\boldsymbol{\beta }}) = {{(\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{\bar a}})}^{\rm{T}}}\mathit{\boldsymbol{\omega }}(\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{\bar a}}) + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}^{\rm{T}}}\mathit{\boldsymbol{W}}(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0} \end{array}} \right. $$ (12) 当给定a值时，式（12）等价于QP问题：
$$ \left\{ {\begin{array}{*{20}{l}} {\min \mathit{\Phi }(\mathit{\boldsymbol{\beta }}) = {{(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}^{\rm{T}}}\mathit{\boldsymbol{W}}(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} \ge 0} \end{array}} \right. $$ (13) 综上所述，可以分别估计β和a的值。给定$ \hat {\bar {\boldsymbol{a}}} $的近似值，组成QP问题（式（13）），求解得到新的$ \hat {\boldsymbol{\beta}} $以后，利用式（10）或式（11）求解新的$ \hat {\bar {\boldsymbol{a}}} $近似值。求解ICPEIV模型转换为求解ICLS模型和代入公式这两个步骤的循环过程。因此，将求解IC-PEIV模型简化成标准LS形式的关键是将ICLS问题（式（13））的求解转换为标准LS形式。
2.2 ICLS平差模型的改进Jacobian迭代法
式（13）是典型的QP问题，根据最优化理论可以得到其获得最优解的KT条件为[17]:
$$ \mathit{\boldsymbol{N\beta }} = {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} + {{\mathit{\boldsymbol{\bar A}}}^{\rm{T}}}\mathit{\boldsymbol{Wy}} $$ (14a) $$ {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}(\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}}) = 0 $$ (14b) $$ \mathit{\boldsymbol{\lambda }} \ge 0 $$ (14c) 式中，$ \boldsymbol{N}=\bar{{\boldsymbol{A}}}^{\mathrm{T}} \boldsymbol{W} \bar{{\boldsymbol{A}}}_{\circ} \text { 由 } \hat{\bar{{\boldsymbol{A}}}^{\mathrm{T}}} \boldsymbol{W} {\boldsymbol{y}}=\boldsymbol{u}_{{\boldsymbol{h}}}+\boldsymbol{u}_{{\boldsymbol{B}}} \text { 和 } \hat{\bar{{\boldsymbol{A}}}^{\mathrm{T}}} \boldsymbol{W} \hat{\bar{{\boldsymbol{A}}}}= \boldsymbol{N}_{h}+\boldsymbol{N}_{{\boldsymbol{B}}}+\boldsymbol{N}_{{\boldsymbol{Bh}} }+\boldsymbol{N}_{{\boldsymbol{h B}}}$，容易得到式（14a)~(14e）与式（5b)~(5e）等价。由式（14a）可得：
$$ \mathit{\boldsymbol{\beta }} = {\mathit{\boldsymbol{N}}^{ - 1}}\left( {{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} + {{\mathit{\boldsymbol{\bar A}}}^{\rm{T}}}\mathit{\boldsymbol{Wy}}} \right) $$ (15) 将式（15）分别代入式（1c）和式（14b），并令$ {\boldsymbol{D}}={\boldsymbol{G N}}^{-1} {\boldsymbol{G}}^{\mathrm{T}}, {\boldsymbol{l}}= {\boldsymbol{d}}- {\boldsymbol{G \hat{\beta}}}_{\mathrm{LS}}, \text { 且 } \hat{{\boldsymbol{\beta}}}_{\mathrm{LS}}= {\boldsymbol{N}}^{-1} \overline{{\boldsymbol{A}}^{\mathrm{T}}} {\boldsymbol{W y}} $为无约束LS解，可得：
$$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{D\lambda }} - \mathit{\boldsymbol{l}} \ge 0}\\ {{\mathit{\boldsymbol{\lambda }}^{\rm{T}}}(\mathit{\boldsymbol{D\lambda }} - \mathit{\boldsymbol{l}}) = 0} \end{array}} \right. $$ (16) 设Dij为D的第i行第j列元素，li为l的第i个分量，λ的第k次迭代值为λ(k)=(λ1(k) λ2(k) … λs(k))T，用Jacobian迭代法求解下列方程：
$$ \mathit{\boldsymbol{D\lambda }} - \mathit{\boldsymbol{l}} = 0 $$ (17) 由此得到第k+1次迭代值为$ {\boldsymbol{\lambda}}^{*(k+1)}= \left(\lambda_{1}^{*(k+1)} \lambda_{2}^{*(k+1)} \cdots \lambda_{s}^{*(k+1)}\right)^{\mathrm{T}}{ }_{\circ} \quad \text { 若 } {\boldsymbol{\lambda}}^{*(k+1)}>0, \text { 令 } {\boldsymbol{\lambda}}^{(k+1)}= {\boldsymbol{\lambda}}^{*(k+1)}$, 易知λ(k + 1)是式（16）的可行解。若λ*(k + 1)的第i个分量λi*(k + 1) < 0，对应的第i个方程满足：
$$ \begin{gathered} D_{i, 1} \lambda_{1}^{*(k+1)}+D_{i, 2} \lambda_{2}^{*(k+1)}+\cdots D_{i, i-1} \lambda_{i-1}^{*(k+1)}+ \\ D_{i, i+1} \lambda_{i+1}^{*(k+1)}+\cdots D_{i, s} \lambda_{s}^{*(k+1)}-l_{i}=-D_{i, i} \lambda_{i}^{*(k+1)}>0 \end{gathered} $$ (18) 假设方阵D为正定矩阵，则主对角元素Di, i > 0，式（18）成立。此时，设λi(k+1)=0，其他分量不变，即λj(k+1)=max(λj*(k+1), 0)(j=1, 2…s)，满足式（16）。重复上述迭代过程，直到最后两次迭代值相等。
经典的等式约束间接平差模型为：
$$ \left\{ {\begin{array}{*{20}{l}} {\min \mathit{\Phi }(\mathit{\boldsymbol{\beta }}) = {{(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}^{\rm{T}}}\mathit{\boldsymbol{W}}(\mathit{\boldsymbol{\bar A\beta }} - \mathit{\boldsymbol{y}})}\\ {{\rm{ s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{G\beta }} - \mathit{\boldsymbol{d}} = 0} \end{array}} \right. $$ (19) 式（19）中等式约束对应的Lagrange乘子应满足的关系与式（17）在形式上完全一致，二者的差别在于，式（19）的Lagrange乘子由式（17）直接求解得到，对λ的取值没有任何要求，最终解λ严格满足式（17）。对于ICLS问题（式（13）），需要迭代求解式（17），获得非负解。由上述分析可知，当λi(k+1) > 0时，对应的第i个方程严格满足：
$$ D_{i, 1} \lambda_{1}^{(k+1)}+D_{i, 2} \lambda_{2}^{(k+1)}+\cdots+D_{i, s} \lambda_{s}^{(k+1)}-l_{i}=0 $$ (20) 此时λi对应的不等式约束为有效约束，等同于等式约束。当λi(k+1)=0时，对应的第i个方程满足：
$$ D_{i, 1} \lambda_{1}^{(k+1)}+D_{i, 2} \lambda_{2}^{(k+1)}+\cdots+D_{i, s} \lambda_{s}^{(k+1)}-l_{i}>0 $$ (21) 此时λi对应的不等式约束为无效约束，对参数估计的结果没有影响。由式（20）和式（21）易知，式（16）恒成立。由上述改进的Jacobian迭代法得到式（16）的非负解，代入式（15），得到式（13）的ICLS解。因此，解ICLS问题的改进Jacobian迭代法是对等式约束LS问题的扩展。求解式（16）的Jacobian迭代法的计算步骤如下：
1）令初始解λ(0)=[0 0…0]T;
2）计算式（17）的第k+1(k=0, 1…N, N为迭代次数）次Jacobian迭代解为：
$$ {\mathit{\boldsymbol{\lambda }}^{*(k + 1)}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{D}}} \right){\mathit{\boldsymbol{\lambda }}^{(k)}} + {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{l}} $$ (22) 式中，M是由D的主对角元素组成的对角阵，则式（22）的分量形式为：
$$ \lambda_{i}^{*(k+1)}=\frac{l_{i}-\sum\limits_{j=1}^{i-1} D_{i j} \lambda_{j}^{(k)}-\sum\limits_{j=i+1}^{s} D_{i j} \lambda_{j}^{(k)}}{D_{i i}} $$ (23) 取$ \lambda_{i}^{(k+1)}=\max \left(0, \lambda_{i}^{*(k+1)}\right) $;
3）若λ(k+1)≠λ(k)，则转步骤2）；否则，停止迭代。令λ=λ(k+1)，根据式（15）求得ICLS解。
由于式（17）是求解不等式约束对应的Lagrange乘子，法方程的维数等于不等式约束的个数，当约束个数较少时，方程的计算量较小，因此，它更适应于不等式约束个数较少的情况。
2.3 ICPEIV模型的经典最小二乘迭代算法
ICPEIV模型的计算步骤如下：
1）给定已知量a、y、h、B、ω、W;
2）令$ \bar{{\boldsymbol{a}}}^{(0)}={\boldsymbol{a}}, \bar{{\boldsymbol{A}}}^{(0)}=\operatorname{ivec}\left({\boldsymbol{h}}+{\boldsymbol{B}} \bar{{\boldsymbol{a}}}^{(0)}\right) $组成QP问题（式（13）），采用改进的Jacobian迭代法求解，得到估值β(0);
3）根据t和n的大小以及使得计算量较小的原则，分别采用式（10）或式（11）更新a(k)
4）如果a或β的误差小于预先确定的容许值，即$ \left\|\hat{\bar{{\boldsymbol{a}}}}^{(k)}-\hat{\bar{{\boldsymbol{a}}}}^{(k-1)}\right\| \leqslant \varepsilon_{1} \text { 或 }\left\|\hat{{\boldsymbol{\beta}}}^{(k)}-\hat{{\boldsymbol{\beta}}}^{(k-1)}\right\| \leqslant \varepsilon_{2}\left(\varepsilon_{1}, \varepsilon_{2}\right. $为足够小的正数）时，停止迭代，否则，转步骤2）。
本文提出的经典LS迭代法也包含内外迭代，外迭代次数是指需要计算的QP问题的次数，内迭代次数是指采用Jacobian迭代法求解每个QP问题的平均迭代次数。而在线性化方法和SQP算法中，外迭代次数分别是指需要计算的QP问题（式（7）和式（9））的次数，内迭代次数指采用积极约束法求解相应QP问题的平均迭代次数，即计算相应等式约束LS问题的个数。
本文方法求解QP问题的参数仅为模型参数，而线性化和SQP方法求解QP问题的参数包括模型参数和系数阵元素。前者的QP问题的维数小于后者，特别是在系数阵元素中随机量较多的情况下。在线性化和SQP方法的内迭代中，求解QP问题一般采用计算效率较好的积极约束法，内迭代一次是指计算一次等式约束LS问题，其复杂度远大于本文方法中的一次内迭代，因为本文方法的一次Jacobian迭代仅仅是公式代入过程。由于本文方法的内外迭代的复杂度都小于线性近似和SQP方法，在比较计算效率时，应同时比较内外迭代次数和计算时间。
3 算例分析
3.1 数值算例1
算例1取自文献[13]，对应的约束EIV模型为：
$$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{y}} = \left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_A}} \right)\mathit{\boldsymbol{\beta }} + {\mathit{\boldsymbol{e}}_y}}\\ {\mathit{\boldsymbol{G\beta }} \ge \mathit{\boldsymbol{d}}} \end{array}} \right. $$ (24) 系数矩阵、右端观测项、一般约束矩阵、约束常数项及区间约束条件见文献[13]。设系数阵元素和右端向量的随机误差都服从独立同分布，即Qe=I25为25阶单位矩阵。不考虑系数阵误差，求解相应的ICLS解作为线性化方法的初值，即$ \hat{\boldsymbol{\beta}} _{{\rm{ICLS}}}= $[-0.100 0 -0.100 0 0.215 2 0.350 2]T。将式（24）转换成式（1a)~(1e），分别采用线性化方法和SQP方法（解式（7）和式（9）均采用积极约束法）及本文经典LS方法求解，迭代终止标准均取ε=1.0×10-8。算例1的参数估值、内外迭代次数、计算时间等模拟实验结果见表 1，其中，$\hat{{\boldsymbol{\beta}}}_{\mathrm{LA}}、\hat{{\boldsymbol{\beta}}}_{\mathrm{SQP}} \text { 和 } \hat{{\boldsymbol{\beta}}}_{\mathrm{CLS}} $分别表示采用线性化方法、SQP算法和经典LS方法的ICWTLS解。
表 1 算例1的模拟实验结果Table 1. Simulation Results of Example 1ICWTLS解 $ \hat {\boldsymbol{\beta}}_1 $ $ \hat {\boldsymbol{\beta}}_2 $ $ \hat {\boldsymbol{\beta}}_3 $ $ \hat {\boldsymbol{\beta}}_4 $ 外迭代次数 平均内迭代次数 计算时间/s $ \hat {\boldsymbol{\beta}}_{{\rm{LA}}} $ −0.100 000 −0.100 000 0.168 581 0.399 765 23 2 0.038 $ \hat {\boldsymbol{\beta}}_{{\rm{SQP}}} $ −0.100 000 −0.100 000 0.168 579 0.399 766 5 2 0.020 $ \hat {\boldsymbol{\beta}}_{{\rm{CLS}}} $ −0.100 000 −0.100 000 0.168 547 0.399 777 22 132 0.025 从表 1中可以看出，本文方法得到的结果与文献[13]中基于式（24）分别采用积极约束法和SQP方法得到的结果相比，它们在小数点后第4位完全一致。将3种方法得到的估值（$ \hat {\bar {\boldsymbol{a}}} , \hat {\boldsymbol{\beta}} $）代入式（2），可以得到Φ($ \hat {\bar {\boldsymbol{a}}} $LA，$ \hat {\boldsymbol{\beta}} $LA)=Φ($ \hat {\bar {\boldsymbol{a}}} $SQP，$ \hat {\boldsymbol{\beta}} $SQP)=0.139 732，Φ($ \hat {\bar {\boldsymbol{a}}} $CLS，$ \hat {\boldsymbol{\beta}} $CLS)=0.139 737。由此可见，在相同的迭代终止标准下，3种方法的计算精度相当。线性化方法和本文方法的外迭代次数相近，本文方法的内迭代次数远大于线性化方法，SQP方法具有二阶收敛性，外迭代次数最少。
对于外迭代，线性化方法和SQP方法需要求解24维的QP问题，而本文方法是4维，外迭代的复杂度明显减小。线性化方法和SQP方法每次内迭代都要解24维的等式约束问题，而本文方法每次内迭代只需要解11维线性方程组。
从计算时间上看，本文方法小于线性近似方法，具有更高的计算效率。由于不等式约束的个数较多，本文方法迭代次数远大于SQP方法，计算时间与SQP方法相比缺乏优势。本文方法的总迭代次数远大于其他两种方法，但计算时间和SQP法相近，进一步验证了本文方法内外迭代的复杂度明显较低。
3.2 数值算例2
由§2.2分析可知，本文方法的计算量随约束个数的减少而下降。算例2采用与算例1相同的观测数据，仅保留3个一般不等式约束，去掉区间约束，计算不考虑系数阵误差的ICLS解为$ \hat {\boldsymbol{\beta}} $ICLS=[0.129 9 -0.575 7 0.425 1 0.243 8]T。将其作为线性近似方法的初始近似值，并采用与算例1相同的计算方案，结果见表 2。
表 2 算例2的模拟实验结果Table 2. Simulation Results of Example 2ICWTLS解 $ \hat {\boldsymbol{\beta}}_1 $ $ \hat {\boldsymbol{\beta}}_2 $ $ \hat {\boldsymbol{\beta}}_3 $ $ \hat {\boldsymbol{\beta}}_4 $ 外迭代次数 平均内迭代次数 计算时间/s $ \hat {\boldsymbol{\beta}}_{{\rm{LA}}} $ 0.127 524 −0.576 759 0.426 986 0.243 459 4 2 0.018 $ \hat {\boldsymbol{\beta}}_{{\rm{SQP}}} $ −0.100 000 −0.100 000 0.168 579 0.399 766 5 2 0.020 $ \hat {\boldsymbol{\beta}}_{{\rm{CLS}}} $ 0.127 524 −0.576 759 0.426 986 0.243 459 14 27 0.009 从表 2可知，当约束个数从算例1的11个变为3个以后，线性化方法和SQP方法需求解24维的QP问题，而本文方法只需求解4维的QP问题，且采用改进的Jacobian迭代法只需求解3维线性方程组。本文方法的内外迭代次数及计算时间都明显减少，验证了本文方法的计算效率随约束个数的减少而减少。因此，在约束个数较少的测绘应用中具有更好的适用性。3种方法在相同的迭代终止标准下得到的结果一致，但本文方法的计算效率比线性化方法和SQP方法更好。
4 结语
目前，求解ICPEIV模型的线性化方法计算效率偏低，转换成QP问题后需采用规划类算法求解。若采用非线性规划算法，如SQP方法，要求目标函数的梯度和Hessian矩阵，形式复杂，难以推广应用。针对上述问题，本文提出了求解ICPEIV模型的经典LS方法，将模型参数和系数阵元素分别估计，模型参数的估计归结为解QP问题，系数阵元素的估计则利用其与模型参数的关系计算。针对QP问题，以KT条件为基础，设计了改进的Jacobian迭代法将QP问题转化为线性方程的求解，它是求解等式约束LS模型的Lagrange乘子法的拓展。本文方法无需对观测方程线性化，也无需掌握复杂的规划类算法，形式简洁，与经典LS一致，易于推广应用。数值实验表明，本文方法具有良好的计算效率和广泛的适用性，能得到与线性化方法和SQP方法精度相当的WTLS解。特别是不等式约束的个数较少时，计算效率明显提高。本文方法将ICPEIV模型及相应的WTLS理论纳入经典LS理论的框架，是对经典LS理论的有益拓展和补充。
-
表 1 算例1的模拟实验结果
Table 1 Simulation Results of Example 1
ICWTLS解 $ \hat {\boldsymbol{\beta}}_1 $ $ \hat {\boldsymbol{\beta}}_2 $ $ \hat {\boldsymbol{\beta}}_3 $ $ \hat {\boldsymbol{\beta}}_4 $ 外迭代次数 平均内迭代次数 计算时间/s $ \hat {\boldsymbol{\beta}}_{{\rm{LA}}} $ −0.100 000 −0.100 000 0.168 581 0.399 765 23 2 0.038 $ \hat {\boldsymbol{\beta}}_{{\rm{SQP}}} $ −0.100 000 −0.100 000 0.168 579 0.399 766 5 2 0.020 $ \hat {\boldsymbol{\beta}}_{{\rm{CLS}}} $ −0.100 000 −0.100 000 0.168 547 0.399 777 22 132 0.025 
下载: 导出CSV
表 2 算例2的模拟实验结果
Table 2 Simulation Results of Example 2
ICWTLS解 $ \hat {\boldsymbol{\beta}}_1 $ $ \hat {\boldsymbol{\beta}}_2 $ $ \hat {\boldsymbol{\beta}}_3 $ $ \hat {\boldsymbol{\beta}}_4 $ 外迭代次数 平均内迭代次数 计算时间/s $ \hat {\boldsymbol{\beta}}_{{\rm{LA}}} $ 0.127 524 −0.576 759 0.426 986 0.243 459 4 2 0.018 $ \hat {\boldsymbol{\beta}}_{{\rm{SQP}}} $ −0.100 000 −0.100 000 0.168 579 0.399 766 5 2 0.020 $ \hat {\boldsymbol{\beta}}_{{\rm{CLS}}} $ 0.127 524 −0.576 759 0.426 986 0.243 459 14 27 0.009
下载: 导出CSV
-
[1] Neitzel F. Generalization of Total Least-Squares on Example of Unweighted and Weighted 2D Similarity Transformation[J]. Journal of Geodesy, 2010, 84 (12): 751-762 doi: 10.1007/s00190-010-0408-0
[2] Akyilmaz O. Total Least Squares Solution of Coordinate Transformation[J]. Survey Review, 2007, 39 (303): 68-80 doi: 10.1179/003962607X165005
[3] Zhou Y, Fang X. A Mixed Weighted Least Squares and Weighted Total Least Squares Adjustment Method and Its Geodetic Applications[J]. Survey Review, 2016, 48 (351): 421-429 doi: 10.1179/1752270615Y.0000000040
[4] Tong X H, Jin Y M, Zhang S L, et al. Bias-Corrected Weighted Total Least-Squares Adjustment of Condition Equations[J]. Journal of American Surveying Engineering, 2014, 141 (2): 04014013 http://d.wanfangdata.com.cn/periodical/26f7db7fb994669a82ad317964b29beb
[5] 陈义, 陆珏, 郑波. 总体最小二乘方法在空间后方交会中的应用[J]. 武汉大学学报·信息科学版, 2008, 33 (12): 1 271-1 274 http://ch.whu.edu.cn/article/id/1776 Chen Yi, Lu Jue, Zheng Bo. Application of Total Least Squares to Space Resection[J]. Geomatics and Information Science of Wuhan University, 2008, 33(12): 1 271-1 274 http://ch.whu.edu.cn/article/id/1776
[6] Xu C J, Wang L Y, Wen Y M, et al. Strain Rates in the Sichuan-Yunnan Region Based upon the Total Least Squares Heterogeneous Strain Model from GPS Data[J]. Terrestrial Atmospheric and Oceanic Sciences, 2011, 22 (2): 133-147 doi: 10.3319/TAO.2010.07.26.02(TibXS)
[7] Fang X. Weighted Total Least Squares: Necessary and Sufficient Conditions, Fixed and Random Parameters[J]. Journal of Geodesy, 2013, 87(8): 733-749 doi: 10.1007/s00190-013-0643-2
[8] Mahboub V. On Weighted Total Least-Squares for Geodetic Transformations[J]. Journal of Geodesy, 2012, 86 (5): 359-367 doi: 10.1007/s00190-011-0524-5
[9] Shen Y Z, Li B F, Chen Y. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85 (4): 229-238 doi: 10.1007/s00190-010-0431-1
[10] Amiri-Simkooei A, Jazaeri S. Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]. Journal of Geodetic Sciences, 2012, 2 (2): 113-124 doi: 10.2478/v10156-011-0036-5
[11] Schaffrin B, Felus Y A. An Algorithmic Approach to the Total Least-Squares Problem with Linear and Quadratic Constraints[J]. Studia Geophysica et Geodaetica, 2009, 53 (1): 1-16 doi: 10.1007/s11200-009-0001-2
[12] Zhang S L, Tong X H, Zhang K. A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J]. Journal of Geodesy, 2013, 87 (1): 23-28 doi: 10.1007/s00190-012-0575-2
[13] Fang X. On Non-combinatorial Weighted Total Least Squares with Inequality Constraints[J]. Journal of Geodesy, 2014, 88 (8): 805-816 doi: 10.1007/s00190-014-0723-y
[14] Xu P L, Liu J N, Shi C. Total Least Squares Adjustment in Partial Errors-in-Variables Models: Algorithm and Statistical Analysis[J]. Journal of Geodesy, 2012, 86 (8): 661-675 doi: 10.1007/s00190-012-0552-9
[15] Shi Y, Xu P L, Liu J N, et al. Alternative Formulae for Parameter Estimation in Partial Errors-inVariables Models[J]. Journal of Geodesy, 2015, 89 (1): 13-16 doi: 10.1007/s00190-014-0756-2
[16] Zeng W X, Liu J N, Yao Y B. On Partial Errors-inVariables Models with Inequality Constraints of Parameters and Variables[J]. Journal of Geodesy, 2015, 89 (2): 111-119 doi: 10.1007/s00190-014-0775-z
[17] 陈宝林. 最优化理论与算法[M]. 北京: 清华大学出版社, 2005 Chen Baolin. Theory and Algorithm in Optimization[M]. Beijing: Tsinghua University Press, 2005
[18] 谢建, 龙四春, 周璀. 不等式约束PEIV模型的最优性条件及SQP算法[J]. 武汉大学学报·信息科学版, 2020, 45 (7): 1 002-1 007 doi: 10.13203/j.whugis20180297 Xie Jian, Long Sichun, Zhou Cui. Optimality Conditions of Inequality Constrained Partial EIV Model and the SQP Algorithm[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 1 002-1 007 doi: 10.13203/j.whugis20180297
-
期刊类型引用(1)
1. Jianjun ZHU，Leyang WANG，Jun HU，Bofeng LI，Haiqiang FU，Yibin YAO. Recent Advances in the Geodesy Data Processing. Journal of Geodesy and Geoinformation Science. 2023(03): 33-45 . 
必应学术
其他类型引用(2)
计量
- 文章访问数: 1514
- HTML全文浏览量: 474
- PDF下载量: 96
- 被引次数: 3
