留言板

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

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

区间约束平差模型的共轭梯度积极集算法

谢雪梅 宋迎春 夏玉国

谢雪梅, 宋迎春, 夏玉国. 区间约束平差模型的共轭梯度积极集算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
引用本文: 谢雪梅, 宋迎春, 夏玉国. 区间约束平差模型的共轭梯度积极集算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
XIE Xuemei, SONG Yingchun, XIA Yuguo. An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints[J]. Geomatics and Information Science of Wuhan University, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
Citation: XIE Xuemei, SONG Yingchun, XIA Yuguo. An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints[J]. Geomatics and Information Science of Wuhan University, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325

区间约束平差模型的共轭梯度积极集算法

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

国家自然科学基金 41574006

国家自然科学基金 41674009

国家自然科学基金 41674012

详细信息
    作者简介:

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

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

An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints

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, lecturer, specializes in theories and methods of surveying data processing. E-mail: xiexuemei_2003@126.com

    Corresponding author: SONG Yingchun, PhD, professor. E-mail: csusyc@qq.com
图(1) / 表(5)
计量
  • 文章访问数:  1094
  • HTML全文浏览量:  80
  • PDF下载量:  92
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-07-01
  • 刊出日期:  2019-09-05

区间约束平差模型的共轭梯度积极集算法

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

    国家自然科学基金 41574006

    国家自然科学基金 41674009

    国家自然科学基金 41674012

    作者简介:

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

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

摘要: 主要研究参数带有区间约束的平差算法,通过把平差问题转化成一个带有区间约束的二次规划问题,利用积极集对二次规划问题进行划分与重组,结合无约束共轭梯度优化算法,给出了带有区间约束的平差算法,并同时给出了参数解的精度评估。由于投影梯度法可以迅速改变积极约束集的构成,新的算法比经典的积极集法效率更高,可以降低模型的不适定性,保持参数先验信息中的统计、几何或物理意义,适合于求解大规模的带有区间约束的平差问题。

English Abstract

谢雪梅, 宋迎春, 夏玉国. 区间约束平差模型的共轭梯度积极集算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
引用本文: 谢雪梅, 宋迎春, 夏玉国. 区间约束平差模型的共轭梯度积极集算法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
XIE Xuemei, SONG Yingchun, XIA Yuguo. An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints[J]. Geomatics and Information Science of Wuhan University, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
Citation: XIE Xuemei, SONG Yingchun, XIA Yuguo. An Active Set Algorithm of Conjugate Gradients for Adjustment Model with Interval Constraints[J]. Geomatics and Information Science of Wuhan University, 2019, 44(9): 1274-1281. doi: 10.13203/j.whugis20170325
  • 大地测量数据处理中,参数本身附加的先验约束信息和物理约束信息常常被用来解决因观测信息不充分导致的病态问题[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维未知参数;Lm维观测向量;Am×n维系数矩阵(mn);em维随机误差向量;$\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=ATPAc=-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,初始点为X0F0=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)

      对应指标集IJ重新划分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)

      式中,XIX中对应下标$i \in I$的分量xi构成的向量,xi=lixi=uiXI的值保持固定,而让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)

      由于NJJN的主子阵,由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-1g(Xk)=NXk +cIk-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]$,NJJs×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)

      式中,约束方程系数矩阵BSs×n阶矩阵,$s$为集合$S$的元素个数,${{\mathit{\boldsymbol{B}}_S}}$由单位向量${\mathit{\boldsymbol{b}}_i}$组成(${\mathit{\boldsymbol{b}}_i}$除了第$i$个分量为1外,其余分量为0),$i \in S$;dS为第i分量等于liui的向量。设${{\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所示的测边网,P1P2为已知点,其坐标分别为(48 580.285, 600 500.496)和(47 943.013, 66 225.845),单位为m。假设算例1中的点P3P4P5P6的真实坐标如表 1所示,边长观测值如表 2所示,由前期的观测结果可知P3P4P5P6的近似坐标(见表 1)以及它们相应的点位精度。先验约束区间D2是通过点位精度进行推算的,D1D3是特别设定的约束区间,其中,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的约束区间都比较大,算法中没有得到有效约束,这时得到的解算结果与最小二乘估计的结果一致,换句话说,就是最优解在区间D1D2之中,因此有${{\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)对它们进行精度评估。

      此算例说明,利用参数的先验约束信息可以改善平差模型的病态性。因此,充分利用参数的先验信息对参数进行约束,可以弥补观测不充分所导致的病态问题,保证解的稳定性。

    • 大地测量中存在一些用区间、集合以及椭球形式描述的约束信息,这些约束信息的统计性质不清楚,只有一个模糊变化的范围,如参数的可行空间、噪声变化的范围等。充分利用它们可以对部分参数进行约束,从而保证参数解的唯一性和稳定性。本文利用区间约束来描述参数的先验信息,提供了一种基于共轭梯度积极集的平差新方法,其算法简单。通过实例说明,利用参数的先验约束信息可以改善平差模型的病态性。在测量数据处理中,充分利用参数的先验信息建立平差模型,可以弥补观测不充分导致的病态问题,提高参数估计的效率。

参考文献 (23)

目录

    /

    返回文章
    返回