-
大地测量实际问题中，除了观测信息，还存在着一些先验约束信息[1-6]，利用它们可以降低模型的不适定性，保持参数先验值和验后估值的统计、几何或物理意义，它们在平差模型中的重要体现就是约束条件[7-10]。不等式约束比等式约束更能准确地描述先验信息，针对先验信息建立科学的不等式约束模型，研究其解算方法，提高参数解的精度，越来越受到国内外学者的关注。解带有不等式约束平差模型的关键是对于不等式约束的处理方法。文献[10-11]采用了凝聚函数法，文献[12]基于有效约束算法，都是采用把不等式约束转化为有效的等式约束进行平差解算的方法。这类方法可以得到一个显式解，并能对解进行精度分析，但此类方法对于选取有效的等式约束过程比较复杂。文献[9]利用凸二次规划理论给出了不等式约束平差的规划类算法，尝试利用优化理论解决不等式约束平差问题，这些方法不同于常规的平差方法，需要掌握一定的优化理论知识，有一定的复杂性。文献[8]基于最优化计算理论中罚函数方法提出了另一种处理不等式约束的平差算法，然而，罚函数方法也不是常规的平差方法。本文针对带有不等式约束的平差模型，在最小二乘平差准则基础上，把不等式约束看成是一个可行域，利用Fisher函数在可行域中快速搜索使误差平方和达最小的参数估计算法。该算法的平差准则与最小二乘平差准则一致，容易被测量工作者接受，可以对维数较大的平差问题进行解算。
-
当参数向量的不确定性用一个区间来描述时，可以建立平差模型为：
$$ \left\{ \begin{array}{l} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}\\ \mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{d}} \end{array} \right.,\;\;\;权阵\;\mathit{\boldsymbol{P}} $$ (1) 式中，A为m×n(m≥n)维设计矩阵；L为n维观测向量；e为n维随机误差向量；X=[x1 x2…xn]T为n维参数向量；G为s×n维矩阵；d为s维常向量。模型(1)被称为参变量带有不等式约束的平差模型。模型(1)中，考虑约束条件进行最小二乘平差，即解‖L-AX‖2的最小值问题，可以等价地求解下面的二次规划问题：
$$ \begin{array}{*{20}{c}} {\min \mathit{\boldsymbol{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)}\\ {{\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{d}}} \end{array} $$ (2) 由于目标函数F(X)=(L-AX)TP(L-AX)=XT(ATPA)X-2LTPAX+LTPL中，LTPL为一个常量，若设N=ATPA，c=-(LTPA)T，则二次规划问题(2)可以转化为如下的二次规划问题：
$$ \begin{array}{*{20}{c}} {\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}}.\;\mathit{\boldsymbol{GX}} \le \mathit{\boldsymbol{d}}} \end{array} $$ (3) 二次规划问题(3)的K-T条件为：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{NX}} + \mathit{\boldsymbol{c}} + {\mathit{\boldsymbol{G}}^{\rm{T}}}\lambda = 0,\;\;\mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{d}} \le 0,}\\ {\lambda \ge 0,\;\lambda \left( {\mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{d}}} \right) = 0} \end{array} $$ (4) 令
$$ \left\{ \begin{array}{l} \mathit{\boldsymbol{B}} = {\left( {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}} \right)^{ - 1}} - \mathit{\boldsymbol{G}}\\ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{G}} = \mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{B}}\\ g\left( \mathit{\boldsymbol{X}} \right) = - \nabla f\left( \mathit{\boldsymbol{X}} \right) = - \left( {\mathit{\boldsymbol{NX}} + \mathit{\boldsymbol{c}}} \right)\\ h\left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{GX}} - \mathit{\boldsymbol{d}}\\ \begin{array}{*{20}{c}} {\mu \left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{B}}g\left( \mathit{\boldsymbol{X}} \right) = - \mathit{\boldsymbol{B}}\left( {\mathit{\boldsymbol{NX}} + \mathit{\boldsymbol{c}}} \right) = }\\ {{{\left( {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\lambda = \lambda } \end{array} \end{array} \right. $$ 式中，▽f(X)为f(X)的梯度。则二次规划问题(3)的K-T条件可以改写为：
$$ \left\{ \begin{array}{l} g\left( \mathit{\boldsymbol{X}} \right) - {\mathit{\boldsymbol{G}}^{\rm{T}}}\lambda = 0,h\left( \mathit{\boldsymbol{X}} \right) \le 0\\ \mu \left( \mathit{\boldsymbol{X}} \right) \ge 0,\mu \left( \mathit{\boldsymbol{X}} \right)h\left( \mathit{\boldsymbol{X}} \right) = 0 \end{array} \right. $$ (5) -
设h(X)的第j个分量为hj(X)，μ(X)的第j个分量为μj(X)，定义Fisher函数：
$$ \begin{array}{*{20}{c}} {{\rho _j}\left( \mathit{\boldsymbol{X}} \right) = \sqrt {\mu _j^2\left( \mathit{\boldsymbol{X}} \right) + h_j^2\left( \mathit{\boldsymbol{X}} \right)} - {\mu _j}\left( \mathit{\boldsymbol{X}} \right) + {h_j}\left( \mathit{\boldsymbol{X}} \right),}\\ {j = 1,2 \cdots n} \end{array} $$ (6) ρ(X)是由ρj(X)组成的列向量，函数ρ(X)有如下两个重要性质：
性质1:ρ(X)≥0。
性质2:ρ(X)=0的充要条件是h(X)≤0，μ(X)≥0，μ(X)h(X)=0。
下面分析模型(3)最优解的一些判别方法。
1) 若可行解X是式(3)的最优解，因为二次规划问题(3)是一个凸二次规划，则K-T条件(5)成立，即存在λ≥0，使得式(5)成立。即：
$$ \left\{ \begin{array}{l} g\left( {\mathit{\boldsymbol{\bar X}}} \right) - {\mathit{\boldsymbol{G}}^{\rm{T}}}\lambda = 0,h\left( {\mathit{\boldsymbol{\bar X}}} \right) \le 0\\ \mu \left( {\mathit{\boldsymbol{\bar X}}} \right) \ge 0,\mu \left( {\mathit{\boldsymbol{\bar X}}} \right)h\left( {\mathit{\boldsymbol{\bar X}}} \right) = 0 \end{array} \right. $$ (7) 再由ρ(X)的性质2，有ρ(X)=0，另外，由M的定义及式(6)有：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right)g\left( {\mathit{\boldsymbol{\bar X}}} \right) = }\\ {g\left( {\mathit{\boldsymbol{\bar X}}} \right) - {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{\lambda }} = 0} \end{array} $$ (8) 2) 若可行解X使得ρ(X)=0，Mg(X)=0，由ρ(X)=0，可知h(X)≤0，μ(X)≥0，μ(X)·h(X)=0，由Mg(X)=0，可知g(X)-GTλ=0，这说明可行解X使得式(6)成立，因此X是式(3)的最优解。
综合上面分析可得出，可行解X是式(3)的最优解的充分必要条件是：
$$ \rho \left( {\mathit{\boldsymbol{\bar X}}} \right) = 0,Mg\left( {\mathit{\boldsymbol{\bar X}}} \right) = 0 $$ (9) 定义函数：
$$ r\left( {\mathit{\boldsymbol{\bar X}}} \right) = \mathit{\boldsymbol{M}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - {\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {\mathit{\boldsymbol{\bar X}}} \right) $$ (10) 显然，若X是式(3)的最优解，式(8)成立，此时r(X)=0。
反之，若r(X)=0，有：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{G}}r\left( {\mathit{\boldsymbol{\bar X}}} \right) = \mathit{\boldsymbol{GM}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {\mathit{\boldsymbol{\bar X}}} \right) = \mathit{\boldsymbol{G}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - }\\ {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{B}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - \mathit{\boldsymbol{G}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {\mathit{\boldsymbol{\bar X}}} \right) = - \rho \left( {\mathit{\boldsymbol{\bar X}}} \right) = 0} \end{array} $$ (11) 由式(10)及ρ(X)=0，可以得到Mg(X)=0，即式(9)成立。由前面的分析知道X是式(3)的最优解，因此，又得到了一个结论：可行解X是式(3)的最优解的充分必要条件是r(X)=0。
-
若X是式(3)的可行解，但不是最优解，则式(9)不会成立，若ρ(X)≠0，以及由ρ(X)性质1知ρ(X)>0，根据式(11)的推导，有Gr(X)=-ρ(X) < 0，此时，G(X+r(X))≤GX≤d，若ρ(X)=0，有Gr(X)=-ρ(X)=0，G(X+r(X))=G≤d，上述两种情形都说明X+r(X)满足不等式约束条件(3)，它是式(3)的可行解X。利用泰勒展开式，再考虑到N的非定性，有：
$$ \begin{array}{*{20}{c}} {f\left( {\mathit{\boldsymbol{\bar X}} + r\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right) = f\left( {\mathit{\boldsymbol{\bar X}}} \right) + \nabla f\left( {\mathit{\boldsymbol{\bar X}}} \right)r\left( {\mathit{\boldsymbol{\bar X}}} \right) + }\\ {\frac{1}{2}{\nabla ^2}f\left( {\mathit{\boldsymbol{\bar X}}} \right){r^2}\left( {\mathit{\boldsymbol{\bar X}}} \right) \le f\left( {\mathit{\boldsymbol{\bar X}}} \right) + \nabla f\left( {\mathit{\boldsymbol{\bar X}}} \right)r\left( {\mathit{\boldsymbol{\bar X}}} \right)} \end{array} $$ (12) 式中，▽f(X)为f(X)的梯度。因为g(X)=-▽f(X)，先讨论gT(X)r(X)的正负性，因为：
$$ \begin{array}{*{20}{c}} {{g^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right)r\left( {\mathit{\boldsymbol{\bar X}}} \right) = {g^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right)\left( {\mathit{\boldsymbol{M}}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - {\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {\mathit{\boldsymbol{\bar X}}} \right)} \right) = }\\ {{g^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right){\mathit{\boldsymbol{M}}^2}g\left( {\mathit{\boldsymbol{\bar X}}} \right) - {g^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right){\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {\mathit{\boldsymbol{\bar X}}} \right) = }\\ {\left\| {\mathit{\boldsymbol{M}}g\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right\| - {\mu ^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right)\rho \left( {\mathit{\boldsymbol{\bar X}}} \right)} \end{array} $$ (13) 下面根据μ(X)，分3种情况进行讨论。
1) 情形1：μ(X)=0。此时有μ(X)h(X)=0，而由约束条件式(3)可知h(X)≤0，因此，根据ρ(X)的性质2有ρ(X)=0，又因X不是最优解，所以由前面的分析知Mg(X)≠0，故gT(X)·r(X)=||Mg(X)||>0。
2) 情形2：μ(X)>0。因X是可行解，对任何的1≤j≤n，有hj(X)≤0，若hj(X) < 0，有-2 μj(X)hj(X)>0，因为：
$$ \begin{array}{*{20}{c}} {{\rho _j}\left( {\mathit{\boldsymbol{\bar X}}} \right) = \sqrt {\mu _j^2\left( {\mathit{\boldsymbol{\bar X}}} \right) + h_j^2\left( {\mathit{\boldsymbol{\bar X}}} \right)} - {\mu _j}\left( {\mathit{\boldsymbol{\bar X}}} \right) + {h_j}\left( {\mathit{\boldsymbol{\bar X}}} \right) = \sqrt {\mu _j^2\left( {\mathit{\boldsymbol{\bar X}}} \right) + h_j^2\left( {\mathit{\boldsymbol{\bar X}}} \right)} - }\\ {\sqrt {\mu _j^2\left( {\mathit{\boldsymbol{\bar X}}} \right) - 2{\mu _j}\left( {\mathit{\boldsymbol{\bar X}}} \right){h_j}\left( {\mathit{\boldsymbol{\bar X}}} \right) + h_j^2\left( {\mathit{\boldsymbol{\bar X}}} \right)} < 0} \end{array} $$ 与ρ(X)的性质1矛盾，所以，hj(X)=0，ρj(X)= $\sqrt {\mu _j^2\left( {\mathit{\boldsymbol{\bar X}}} \right)} $-μj(X)=0，此即ρ(X)=0，同情形1的分析有gT(X)r(X)=||Mg(X)||>0。
3) 情形3：μ(X) < 0。对任何的1≤j≤n，有μj(X) < 0，由ρ(X)的性质1和性质2，立即可知ρj(X)>0，则有：
$$ \begin{array}{*{20}{c}} {g_j^{\rm{T}}\left( {\mathit{\boldsymbol{\bar X}}} \right){r_j}\left( {\mathit{\boldsymbol{\bar X}}} \right) = \left\| {\mathit{\boldsymbol{M}}g_j^{\rm{T}}\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right\| - \mu _j^{\rm{T}}\left( {\mathit{\boldsymbol{\bar X}}} \right)\rho \left( {\mathit{\boldsymbol{\bar X}}} \right) > }\\ {\left\| {\mathit{\boldsymbol{M}}g_j^{\rm{T}}\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right\| > 0} \end{array} $$ 综合上面3种情形，可知：
$$ \nabla f\left( {\mathit{\boldsymbol{\bar X}}} \right)r\left( {\mathit{\boldsymbol{\bar X}}} \right) = - {g^{\rm{T}}}\left( {\mathit{\boldsymbol{\bar X}}} \right)r\left( {\mathit{\boldsymbol{\bar X}}} \right) < 0 $$ (14) 结合式(12)有：
$$ f\left( {\mathit{\boldsymbol{\bar X}} + r\left( {\mathit{\boldsymbol{\bar X}}} \right)} \right) < f\left( {\mathit{\boldsymbol{\bar X}}} \right) $$ (15) 由X+r(X)的可行性和式(15)可知，r(X)是往目标函数减小方向的增量。因此，可以构建下面的算法求二次规划问题(3)的最优解。
设初始点为X0，迭代终止误差为ε，k=0。
1) 确定投影矩阵：
$$ \mathit{\boldsymbol{M}} = \mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{G}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{G}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{G}} $$ 2) 计算ρ(Xk)，其中，
$$ \begin{array}{*{20}{c}} {{\rho _j}\left( {{\mathit{\boldsymbol{X}}_k}} \right) = \sqrt {\mu _j^2\left( {{\mathit{\boldsymbol{X}}_k}} \right) + h_j^2\left( {{\mathit{\boldsymbol{X}}_k}} \right)} - {\mu _j}\left( {{\mathit{\boldsymbol{X}}_k}} \right) + }\\ {{h_j}\left( {{\mathit{\boldsymbol{X}}_k}} \right),j = 1,2 \cdots n} \end{array} $$ 3) 计算
$$ g\left( {{\mathit{\boldsymbol{X}}_k}} \right) = - \left( {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_k} + c} \right) $$ $$ r\left( {{\mathit{\boldsymbol{X}}_k}} \right) = \mathit{\boldsymbol{M}}g\left( {{\mathit{\boldsymbol{X}}_k}} \right) - {\mathit{\boldsymbol{B}}^{\rm{T}}}\rho \left( {{\mathit{\boldsymbol{X}}_k}} \right) $$ 假若||r(Xk)|| < ε，Xk为K-T点，则算法终止。
4) 令Xk+1=Xk+r(Xk)，k=k+1，转步骤2)。
-
由前面的分析可知，算法中每一步得到的Xk都是可行解。若在有限步以后有r(Xk)=0，则Xk就是问题(3)的最优解；否则，就会产生一个可行解序列{Xk}，由前面的分析有：
$$ f\left( {{\mathit{\boldsymbol{X}}_0}} \right) < f\left( {{\mathit{\boldsymbol{X}}_1}} \right) < \cdots < f\left( {{\mathit{\boldsymbol{X}}_k}} \right) < \cdots $$ (16) 因为问题(2)的目标函数F(X)=(L-AX)TP·(L-AX)≥0，而式(3)的目标函数f(X)=$\frac{1}{2}$·F(X)-$\frac{1}{2}$LTPL，LTPL是一个常值，这说明对所有的可行解有：
$$ f\left( \mathit{\boldsymbol{X}} \right) \ge - \frac{1}{2}{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $$ (17) 式(16)和式(17)说明序列{f(Xk)是一个单调递降有下界的序列，其极限存在，因此，算法是收敛的。若在算法的步骤4)中，令Xk+1=Xk+αr(Xk)，若Xk是式(3)的可行解，但不是最优解，则式(9)不会成立，若ρ(Xk)≠0，以及由ρ(Xk)性质1知ρ(Xk)>0，根据式(11)的推导，有Gr(Xk)=-ρ(Xk) < 0。此时，G(Xk+αr(Xk))≤GXk≤d，若ρ(Xk)=0，有Gr(Xk)=-ρ(Xk)=0，G(Xk+αr(Xk))=GXk≤d。上述两种情形都说明Xk+αr(Xk)满足不等式约束条件式(3)，它是式(3)的可行解。为了保证步长α的增量中不跨过最小值点，要求步长α取到一元函数φ(α)=f(Xk+αr(Xk))的极小值处，由φ′(α)=0，可知:
$$ \begin{array}{*{20}{c}} {\alpha = - \frac{{{{\left( {\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_k} + \mathit{\boldsymbol{c}}} \right)}^{\rm{T}}}r\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}{{{r^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_k}r\left( {{\mathit{\boldsymbol{X}}_k}} \right)}} = }\\ {\frac{{{g^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)r\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}{{{r^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{X}}_k}r\left( {{\mathit{\boldsymbol{X}}_k}} \right)}}} \end{array} $$ (18) 如果每一步都取到最小值点的地方，是一种精确的搜索方法，这显然会增加计算量。在实际计算中为了减小计算量，可以采用非精确的搜索方法，只要求步长α使得f(Xk+αr(Xk))比f(Xk)有一定的下降，计算较容易实现。常用的非精确算法有Armijo-Goldstein型线性搜索[13]、Wolfe-Powell型线性搜索[14-16]。若存在σ1和σ2满足0 < σ1 < 0.5，σ1 < σ2 < 1，使得：
$$ \begin{array}{*{20}{c}} {f\left( {{\mathit{\boldsymbol{X}}_k} + {\alpha _k}r\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right) \le f\left( {{\mathit{\boldsymbol{X}}_k}} \right) + }\\ {{\sigma _1}{\alpha _k}\nabla {f^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)r\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \end{array} $$ (19) $$ \nabla {f^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k} + {\alpha _k}r\left( {{\mathit{\boldsymbol{X}}_k}} \right)} \right) \ge {\sigma _2}\nabla {f^{\rm{T}}}\left( {{\mathit{\boldsymbol{X}}_k}} \right)r\left( {{\mathit{\boldsymbol{X}}_k}} \right) $$ (20) 对于非精确的搜索方法，可以采用Wolfe-Powell型线性搜索算法[17]。
-
为了验证算法的有效性，本文建立了一个GPS定位的算例，有3个已知点和一个未知点，如图 1所示, 已知点A(-2 191 820，5 182 440，2 993 206)、B(-2 191 860，5 182 440，2 993 206)、C(-2 191 840，5 182 480，2 993 206)，未知点P近似坐标为(-2 191 840，5 182 480，2 993 206)。P已经通过其他定位方式得到了先验信息，它在圆A、B、C的内部，3个球的半径分别为RA=50, RB=50, RC=10。
图 1 已知点和未知点的位置
Figure 1. Location of Known and Unknown Points
通过线性化观测方程和约束方程，可以得到带有不等式约束的平差模型：L=AX+e, GX≤d，其中，X为近似坐标的改正数，L和A分别为线性化后的观测向量和观测矩阵，G和d分别为距离约束线性化后的系数矩阵和常数向量，相应的数据见表 1。为了更好地分析算法的效果，特意将矩阵A的第一行数据与第二行数据修改为几乎完全相同。实例的解算采用带有步长的搜索法，其中，设置初始可行解为X0=[1 1 1]T，迭代终止误差为ε=0.000 1，目标函数值$F\left( {\mathit{\boldsymbol{\hat X}}} \right){\rm{ = }}{\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}}} \right)^{\rm{T}}}\mathit{\boldsymbol{P}}\left( {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}}} \right)$。表 2给出了系数矩阵非病态情形下5种不同算法的结果比较，在观测信息充分的情形下，系数矩阵非病态，满足罚函数法、有效约束法以及规划类方法的条件，所以解都没有出现异常。从表 2可以看出，最小二乘目标函数值为0.104 2, 它的可行域没有限制，真值所对应的目标函数值为3.839 0，因此，它不是最优目标函数值。在满足约束条件GX≤d时，罚函数法、有效约束法以及规划类方法的目标函数值都为3.729 4，3个解都满足约束条件。因为本文算法是一种搜索算法，不是寻找满足K-T条件的解，因此，得到的目标值要比其他方法得到的目标值更小。
表 1 观测方程及约束方程的相关数据
Table 1. Data About Observation Equation and Constraint Equation
A L G d 0.607 6 -0.434 8 -0.664 7 -1.113 7 -0.949 9 0.308 9 0.048 1 -2.649 6 0.607 6 -0.434 8 -0.664 6 -1.058 0 0.931 3 0.363 7 0.021 4 10.346 7 0.976 0 -0.095 2 -0.195 8 0.409 1 -0.680 9 0.111 4 -8.810 8 6.189 2 -0.625 6 -0.759 9 -0.176 9 0.585 3 -0.515 7 -0.309 0 -0.799 1 -3.106 5 -0.056 7 -0.810 3 0.583 2 4.447 8 -0.500 6 -0.775 2 0.385 3 2.867 6 0.667 3 -0.391 7 0.633 4 4.506 1 -0.821 8 -0.557 4 -0.118 3 0.250 2 表 2 法方程非病态下几种算法结果的比较
Table 2. Results Comparison of Several Algorithms Under the Normal Equation is not Ill
参数 真值 最小二乘法 罚函数法 有效约束法 规划类方法 本文方法 2.029 5 1.030 3 1.983 6 1.983 6 1.983 5 1.982 7 ${\mathit{\boldsymbol{\hat X}}}$ -3.018 2 -2.544 5 -3.135 5 -3.135 5 -3.135 5 -3.135 2 4.378 0 4.213 9 4.225 0 4.225 0 4.225 0 4.225 2 $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 3.839 0 0.104 2 3.729 4 3.729 4 3.729 4 3.722 9 -2.649 6 -1.562 0 -2.649 6 -2.649 6 -2.649 6 -2.649 6 $\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.886 0 0.124 3 0.797 4 0.797 4 0.797 3 0.796 7 1.451 1 2.065 4 1.358 5 1.358 5 1.358 7 1.359 3 表 3是取表 1中系数矩阵A的前3行和相应的观测向量L的前3行构成的观测方程，此时，约束方程仍然是表 1中的GX≤d。由于有意地修改了A的前两行，使之高度相关，导致观测方程病态，这时的最小二乘解异常，虽然增加了约束信息，但罚函数法、有效约束方法、规划类方法要求矩阵A是非病态的。当出现病态时，在确定有效约束时，算法将会出现异常，得到的有效约束并不是真正的有效约束，本来应该是约束不等式GX≤d中第1个约束是有效约束，却变成第3个式子是有效约束了。最优值虽然相同，但是本文算法不要求A非病态，解更接近于真值。
表 3 法方程病态情形下几种算法结果的比较
Table 3. Results Comparison of Several Algorithms Under the Normal Equation is Ill
参数 真值 最小二乘法 罚函数法 有效约束法 规划类方法 本文方法 2.029 5 33.986 6 1.488 3 -31.107 4 1.488 4 1.165 0 ${\mathit{\boldsymbol{\hat X}}}$ -3.018 2 -801.458 1 -13.902 8 776.016 4 -13.904 4 -6.029 6 4.378 0 557.000 0 12.089 2 -534.457 1 12.090 3 6.643 5 $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 2.046 1 1.90×10-15 0.001 5 0.001 5 0.001 5 0.001 5 -2.649 6 -253.062 5 -5.126 9 -5.126 9 -5.127 4 -2.649 6 $\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.886 0 -247.911 8 -3.411 9 -3.411 6 -3.412 4 -0.965 8 1.451 1 290.788 4 6.189 2 6.189 2 6.189 2 3.344 3 -
大地测量平差模型中的参数通常存在一些不确定的附加信息或先验信息(参数内在的相关性、精度、几何或物理信息)，充分利用它们可以对部分参数进行约束，从而保证参数解的唯一性和稳定性。这些附加信息虽然没有实际观测值的绝对精度或可靠性，但可以降低模型的不适定性，选择合适的解空间，保持参数先验值和验后估值的统计、几何或物理意义，有效提高参数估计的效率。本文针对参数带有不等式约束信息，提供了一种新的平差方法，其算法实际仍然是一种规划类的方法，寻找使K-T条件成立的解，因此平差解的精度评估与规划类方法、有效约束方法的精度评估方法相同。本文算法中没有用到N=APAT的逆矩阵，因此A的病态性对算法影响较小。本文算法简单，可以很好地应用于测量数据处理，提供参数估计的准确率。
-
摘要: 大地测量中常存在一些先验不等式约束信息，充分利用它们可以保证参数解的唯一性和稳定性。然而，现有的不等式约束平差算法主要是基于优化理论，算法通常比较复杂，需要选取有效约束或建立罚函数。在最小二乘平差准则基础上，把不等式约束看成是一个可行域，借助Fisher函数在可行域中快速搜索使误差平方和达到最小的最优解，推导出了可行解为最优解的充分必要条件。建立了基于Wolfe-Powell算法的非精确快速搜索算法，从而减小了搜索算法的计算量，得到了一种新的不等式约束平差计算方法。该算法的平差准则与最小二乘平差准则一致，不需要矩阵求逆运算，可适用于维数较大的平差问题解算。Abstract: There usually exists some prior information with inequality constraint in the survey of adjustment model. The uniqueness and stability of the solution can be guaranteed by making full use of it. However, the existing adjustment algorithms with inequality constrain, which are mainly based on optimization theory, are usually complex. They need to select the effective constraint or establish penalty function. This paper mainly studies the adjustment model with inequality constraint, in which the inequality constraint is considered as a feasible region on the basis of the least squares adjustment rule and the Fisher function is used to search the optimal solution that minimizes the sum of squared errors, and sufficient necessary conditions for the optimal feasible solution are derived. A non-precise fast search based on Wolfe-Powell algorithm is given in the feasible region, which reduces the computational complexity, a new adjustment algorithm with inequality constraint is presented. The given algorithm, in which the adjustment criterion is consistent with that of the least squares adjustment criteria, does not require matrix inversion operation, and can solve some of the large dimension adjustment problem with inequality constrain.
-
Key words:
- prior information /
- inequality constraint /
- adjustment model /
- quadratic programming /
- Fisher function
-
表 1 观测方程及约束方程的相关数据
Table 1. Data About Observation Equation and Constraint Equation
A L G d 0.607 6 -0.434 8 -0.664 7 -1.113 7 -0.949 9 0.308 9 0.048 1 -2.649 6 0.607 6 -0.434 8 -0.664 6 -1.058 0 0.931 3 0.363 7 0.021 4 10.346 7 0.976 0 -0.095 2 -0.195 8 0.409 1 -0.680 9 0.111 4 -8.810 8 6.189 2 -0.625 6 -0.759 9 -0.176 9 0.585 3 -0.515 7 -0.309 0 -0.799 1 -3.106 5 -0.056 7 -0.810 3 0.583 2 4.447 8 -0.500 6 -0.775 2 0.385 3 2.867 6 0.667 3 -0.391 7 0.633 4 4.506 1 -0.821 8 -0.557 4 -0.118 3 0.250 2 
下载: 导出CSV
表 2 法方程非病态下几种算法结果的比较
Table 2. Results Comparison of Several Algorithms Under the Normal Equation is not Ill
参数 真值 最小二乘法 罚函数法 有效约束法 规划类方法 本文方法 2.029 5 1.030 3 1.983 6 1.983 6 1.983 5 1.982 7 ${\mathit{\boldsymbol{\hat X}}}$ -3.018 2 -2.544 5 -3.135 5 -3.135 5 -3.135 5 -3.135 2 4.378 0 4.213 9 4.225 0 4.225 0 4.225 0 4.225 2 $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 3.839 0 0.104 2 3.729 4 3.729 4 3.729 4 3.722 9 -2.649 6 -1.562 0 -2.649 6 -2.649 6 -2.649 6 -2.649 6 $\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.886 0 0.124 3 0.797 4 0.797 4 0.797 3 0.796 7 1.451 1 2.065 4 1.358 5 1.358 5 1.358 7 1.359 3
下载: 导出CSV
表 3 法方程病态情形下几种算法结果的比较
Table 3. Results Comparison of Several Algorithms Under the Normal Equation is Ill
参数 真值 最小二乘法 罚函数法 有效约束法 规划类方法 本文方法 2.029 5 33.986 6 1.488 3 -31.107 4 1.488 4 1.165 0 ${\mathit{\boldsymbol{\hat X}}}$ -3.018 2 -801.458 1 -13.902 8 776.016 4 -13.904 4 -6.029 6 4.378 0 557.000 0 12.089 2 -534.457 1 12.090 3 6.643 5 $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 2.046 1 1.90×10-15 0.001 5 0.001 5 0.001 5 0.001 5 -2.649 6 -253.062 5 -5.126 9 -5.126 9 -5.127 4 -2.649 6 $\mathit{\boldsymbol{G}}\left( {\mathit{\boldsymbol{\hat X}}} \right)$ 0.886 0 -247.911 8 -3.411 9 -3.411 6 -3.412 4 -0.965 8 1.451 1 290.788 4 6.189 2 6.189 2 6.189 2 3.344 3
下载: 导出CSV
-
[1] 舒梦珵, 王彦飞.储层重力密度反演后验约束正则化方法[J].地球物理学报, 2015, 58(6):2079-2086 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=DQWX201506022&dbname=CJFD&dbcode=CJFQ Shu Mengcheng, Wang Yanfei. The Posterior Constrained Regularization Method for Reservoir Density Inversion[J]. Chinese Journal Geophysics, 2015, 58(6):2079-2086 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=DQWX201506022&dbname=CJFD&dbcode=CJFQ [2] Wang Y F, Ma S Q, Yang H, et.al. On the Effective Inversion by Imposing a Priori Information for Retrieval of Land Surface Parameters[J]. Science in China D, 2009, 52(4):540-549 doi: 10.1007/s11430-009-0036-9 [3] 朱建军, 丁晓利, 陈永奇.集成地质、力学信息和监测数据的滑坡动态模型[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 [4] 左廷英, 宋迎春, 朱建军.带有先验约束信息边坡变形监测滤波算法[J].湖南大学学报(自然科学版), 2011, 38(2):18-22 http://d.old.wanfangdata.com.cn/Periodical/hndxxb201102004 Zuo Tingying, Song Yingchun, Zhu Jianjun. Filter Algorithm for Slope Monitoring with Prior Constrained Information[J]. Journal of Hunan University (Natural Sciences), 2011, 38(2):18-22 http://d.old.wanfangdata.com.cn/Periodical/hndxxb201102004 [5] 张勤, 黄观文, 王利, 等.附有系统参数和附加约束条件的GPS城市沉降监测网数据处理方法研究[J].武汉大学学报·信息科学版, 2009, 34(3):269-272 http://ch.whu.edu.cn/CN/Y2009/V34/I3/269 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/Y2009/V34/I3/269 [6] Li Bofeng, Shen Yunzhong, Feng Yanming. Fast GNSS Ambiguity Resolution as an Ill-Posed Problem[J]. Journal of Geodesy, 2010, 84(11):683-698 doi: 10.1007/s00190-010-0403-5 [7] 王乐洋, 许才军, 汪建军.附有病态约束矩阵的等式约束反演问题研究[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 [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].测绘学报, 2007, 36(1):50-55 http://d.old.wanfangdata.com.cn/Periodical/chxb200701009 Peng Junhuan, Zhang Yali, Zhang Hongping, et a1.The Solution of Inequality-Constrained Least Squares Problem and Its Statistical Properties[J]. Acta Geodaetica et Cartographica Sinica, 2007, 36(1):50-55 http://d.old.wanfangdata.com.cn/Periodical/chxb200701009 [11] Peng J. An Aggregate Constraint Method for Inequality-Constrained Least Squares Problem[J].Journal of Geodesy, 2006, 79(12):705-713 doi: 10.1007/s00190-006-0026-z [12] 冯光财, 朱建军, 陈正阳, 等.基于有效约束的附不等式约束平差的一种新算法[J].测绘学报, 2007, 36(2):120-123 http://d.old.wanfangdata.com.cn/Periodical/chxb200702001 Feng Guangcai, Zhu Jianjun, Chen Zhengyang, et a1.A New Approach to Inequality Constrained Least-Squares Adjustment[J]. Acta Geodaetica et Cartogra phica Sinica, 2007, 36(2):120-123 http://d.old.wanfangdata.com.cn/Periodical/chxb200702001 [13] Goldstein A A, Price J F. An Effective Algorithm for Minimization[J]. Numerische Mathematic, 1967, 10(2):184-189 http://d.old.wanfangdata.com.cn/OAPaper/oai_doaj-articles_124b6c8b96541f8dfa12eaba2310f8dc [14] Wolfe P. Convergence Conditions for Ascent Methods[J]. SIAM Review, 1969, 11(2):226-235 doi: 10.1137/1011036 [15] Wolfe P. Convergence Conditions for Ascent Methods Ⅱ[J]. SIAM Review, 1971, 13(2):185-188 doi: 10.1137/1013035 [16] Powell M J D. On the Convergence of the Variable Metric Algorithm[J]. IMA Journal of Applied Mathematics, 1971, 7(1):21-36 doi: 10.1093/imamat/7.1.21 [17] 袁亚湘.非线性优化计算方法[M].北京:科学出版社, 2008 Yuan Yaxiang. Nonlinear Optimization Computatio-nal Methods[M]. Beijing:The Science Publishing Company, 2008 -
点击查看大图
图(1) / 表(3)
计量
- 文章访问数: 1268
- HTML全文浏览量: 91
- PDF下载量: 162
- 被引次数: 0
