留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

附不等式约束平差模型的一种快速搜索算法

谢雪梅 宋迎春 肖兆兵

谢雪梅, 宋迎春, 肖兆兵. 附不等式约束平差模型的一种快速搜索算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
引用本文: 谢雪梅, 宋迎春, 肖兆兵. 附不等式约束平差模型的一种快速搜索算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
XIE Xuemei, SONG Yingchun, XIAO Zhaobing. A Fast Search Algorithm in Adjustment Model with Inequality Constraint[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
Citation: XIE Xuemei, SONG Yingchun, XIAO Zhaobing. A Fast Search Algorithm in Adjustment Model with Inequality Constraint[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435

附不等式约束平差模型的一种快速搜索算法

doi: 10.13203/j.whugis20160435
基金项目: 

国家自然科学基金 41574006

国家自然科学基金 41674009

国家自然科学基金 41674012

详细信息
    作者简介:

    谢雪梅, 博士生, 讲师, 主要从事测量数据处理的理论与方法研究。xiexuemei_2003@126.com

    通讯作者: 宋迎春, 博士, 教授。csusyc@qq.com
  • 中图分类号: P207

A Fast Search Algorithm in Adjustment Model with Inequality Constraint

Funds: 

The National Natural Science Foundation of China 41574006

The National Natural Science Foundation of China 41674009

The National Natural Science Foundation of China 41674012

More Information
    Author Bio:

    XIE Xuemei, PhD candidate, lecturer, specializes in surveying data processing. E-mail:xiexuemei_2003@126.com

    Corresponding author: SONG Yingchun, PhD, professor. E-mail:csusyc@qq.com
图(1) / 表(3)
计量
  • 文章访问数:  984
  • HTML全文浏览量:  48
  • PDF下载量:  159
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-04-23
  • 刊出日期:  2018-09-05

附不等式约束平差模型的一种快速搜索算法

doi: 10.13203/j.whugis20160435
    基金项目:

    国家自然科学基金 41574006

    国家自然科学基金 41674009

    国家自然科学基金 41674012

    作者简介:

    谢雪梅, 博士生, 讲师, 主要从事测量数据处理的理论与方法研究。xiexuemei_2003@126.com

    通讯作者: 宋迎春, 博士, 教授。csusyc@qq.com
  • 中图分类号: P207

摘要: 大地测量中常存在一些先验不等式约束信息,充分利用它们可以保证参数解的唯一性和稳定性。然而,现有的不等式约束平差算法主要是基于优化理论,算法通常比较复杂,需要选取有效约束或建立罚函数。在最小二乘平差准则基础上,把不等式约束看成是一个可行域,借助Fisher函数在可行域中快速搜索使误差平方和达到最小的最优解,推导出了可行解为最优解的充分必要条件。建立了基于Wolfe-Powell算法的非精确快速搜索算法,从而减小了搜索算法的计算量,得到了一种新的不等式约束平差计算方法。该算法的平差准则与最小二乘平差准则一致,不需要矩阵求逆运算,可适用于维数较大的平差问题解算。

English Abstract

谢雪梅, 宋迎春, 肖兆兵. 附不等式约束平差模型的一种快速搜索算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
引用本文: 谢雪梅, 宋迎春, 肖兆兵. 附不等式约束平差模型的一种快速搜索算法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
XIE Xuemei, SONG Yingchun, XIAO Zhaobing. A Fast Search Algorithm in Adjustment Model with Inequality Constraint[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
Citation: XIE Xuemei, SONG Yingchun, XIAO Zhaobing. A Fast Search Algorithm in Adjustment Model with Inequality Constraint[J]. Geomatics and Information Science of Wuhan University, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
  • 大地测量实际问题中,除了观测信息,还存在着一些先验约束信息[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)

      式中,Am×n(mn)维设计矩阵;Ln维观测向量;en维随机误差向量;X=[x1  x2xn]Tn维参数向量;Gs×n维矩阵;ds维常向量。模型(1)被称为参变量带有不等式约束的平差模型。模型(1)中,考虑约束条件进行最小二乘平差,即解‖L-AX2的最小值问题,可以等价地求解下面的二次规划问题:

      $$ \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=ATPAc=-(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,μ(Xh(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))≤GXd,若ρ(X)=0,有Gr(X)=-ρ(X)=0,G(X+r(X))=Gd,上述两种情形都说明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(Xr(X)=||Mg(X)||>0。

      2) 情形2:μ(X)>0。因X是可行解,对任何的1≤jn,有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≤jn,有μ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}$LTPLLTPL是一个常值,这说明对所有的可行解有:

      $$ 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))≤GXkd,若ρ(Xk)=0,有Gr(Xk)=-ρ(Xk)=0,G(Xk+αr(Xk))=GXkd。上述两种情形都说明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已经通过其他定位方式得到了先验信息,它在圆ABC的内部,3个球的半径分别为RA=50, RB=50, RC=10。

      图  1  已知点和未知点的位置

      Figure 1.  Location of Known and Unknown Points

      通过线性化观测方程和约束方程,可以得到带有不等式约束的平差模型:L=AX+e, GXd,其中,X为近似坐标的改正数,LA分别为线性化后的观测向量和观测矩阵,Gd分别为距离约束线性化后的系数矩阵和常数向量,相应的数据见表 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,因此,它不是最优目标函数值。在满足约束条件GXd时,罚函数法、有效约束法以及规划类方法的目标函数值都为3.729 4,3个解都满足约束条件。因为本文算法是一种搜索算法,不是寻找满足K-T条件的解,因此,得到的目标值要比其他方法得到的目标值更小。

      表 1  观测方程及约束方程的相关数据

      Table 1.  Data About Observation Equation and Constraint Equation

      ALGd
      0.607 6-0.434 8-0.664 7-1.113 7-0.949 90.308 90.048 1-2.649 6
      0.607 6-0.434 8-0.664 6-1.058 00.931 30.363 70.021 410.346 7
      0.976 0-0.095 2-0.195 80.409 1-0.680 90.111 4-8.810 86.189 2
      -0.625 6-0.759 9-0.176 90.585 3
      -0.515 7-0.309 0-0.799 1-3.106 5
      -0.056 7-0.810 30.583 24.447 8
      -0.500 6-0.775 20.385 32.867 6
      0.667 3-0.391 70.633 44.506 1
      -0.821 8-0.557 4-0.118 30.250 2

      表 2  法方程非病态下几种算法结果的比较

      Table 2.  Results Comparison of Several Algorithms Under the Normal Equation is not Ill

      参数真值最小二乘法罚函数法有效约束法规划类方法本文方法
      2.029 51.030 31.983 61.983 61.983 51.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 04.213 94.225 04.225 04.225 04.225 2
      $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$3.839 00.104 23.729 43.729 43.729 43.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 00.124 30.797 40.797 40.797 30.796 7
      1.451 12.065 41.358 51.358 51.358 71.359 3

      表 3是取表 1中系数矩阵A的前3行和相应的观测向量L的前3行构成的观测方程,此时,约束方程仍然是表 1中的GXd。由于有意地修改了A的前两行,使之高度相关,导致观测方程病态,这时的最小二乘解异常,虽然增加了约束信息,但罚函数法、有效约束方法、规划类方法要求矩阵A是非病态的。当出现病态时,在确定有效约束时,算法将会出现异常,得到的有效约束并不是真正的有效约束,本来应该是约束不等式GXd中第1个约束是有效约束,却变成第3个式子是有效约束了。最优值虽然相同,但是本文算法不要求A非病态,解更接近于真值。

      表 3  法方程病态情形下几种算法结果的比较

      Table 3.  Results Comparison of Several Algorithms Under the Normal Equation is Ill

      参数真值最小二乘法罚函数法有效约束法规划类方法本文方法
      2.029 533.986 61.488 3-31.107 41.488 41.165 0
      ${\mathit{\boldsymbol{\hat X}}}$-3.018 2-801.458 1-13.902 8776.016 4-13.904 4-6.029 6
      4.378 0557.000 012.089 2-534.457 112.090 36.643 5
      $f\left( {\mathit{\boldsymbol{\hat X}}} \right)$2.046 11.90×10-150.001 50.001 50.001 50.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 1290.788 46.189 26.189 26.189 23.344 3
    • 大地测量平差模型中的参数通常存在一些不确定的附加信息或先验信息(参数内在的相关性、精度、几何或物理信息),充分利用它们可以对部分参数进行约束,从而保证参数解的唯一性和稳定性。这些附加信息虽然没有实际观测值的绝对精度或可靠性,但可以降低模型的不适定性,选择合适的解空间,保持参数先验值和验后估值的统计、几何或物理意义,有效提高参数估计的效率。本文针对参数带有不等式约束信息,提供了一种新的平差方法,其算法实际仍然是一种规划类的方法,寻找使K-T条件成立的解,因此平差解的精度评估与规划类方法、有效约束方法的精度评估方法相同。本文算法中没有用到N=APAT的逆矩阵,因此A的病态性对算法影响较小。本文算法简单,可以很好地应用于测量数据处理,提供参数估计的准确率。

参考文献 (17)

目录

    /

    返回文章
    返回