留言板

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

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

不等式约束PEIV模型的经典最小二乘方法

谢建 龙四春 周璀

谢建, 龙四春, 周璀. 不等式约束PEIV模型的经典最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
引用本文: 谢建, 龙四春, 周璀. 不等式约束PEIV模型的经典最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
XIE Jian, LONG Si-chun, ZHOU Cui. Classical Least Squares Method for Inequality Constrained PEIV Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
Citation: XIE Jian, LONG Si-chun, ZHOU Cui. Classical Least Squares Method for Inequality Constrained PEIV Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196

不等式约束PEIV模型的经典最小二乘方法

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

国家自然科学基金 41704007

国家自然科学基金 41877283

国家自然科学基金 42074016

详细信息
    作者简介:

    谢建, 博士, 主要从事测量平差及数据处理研究。hsiejian841006@163.com

  • 中图分类号: P207

Classical Least Squares Method for Inequality Constrained PEIV Model

Funds: 

The National Natural Science Foundation of China 41704007

The National Natural Science Foundation of China 41877283

The National Natural Science Foundation of China 42074016

More Information
    Author Bio:

    XIE Jian, PhD, specializes in surveying adjustment and data processing. E-mail: hsiejian841006@163.com

  • 摘要: 不等式约束部分变量含误差(partial errors-in-variables, PEIV)模型目前主要采用线性化方法和非线性规划类算法, 前者计算效率较低, 后者基于最优化理论, 计算复杂, 未能与经典平差理论建立联系, 难以在测量实际中推广。在整体最小二乘准则下, 根据最优解的Kuhn-Tucker条件, 将不等式约束整体最小二乘解的计算转化为二次规划问题, 并提出改进的Jacobian迭代法求解二次规划。所提方法不需要对观测方程线性化, 与经典最小二乘法具有相同的形式, 易于编程实现。数值实例表明, 所提方法形式简洁, 具有良好的计算效率, 是经典最小二乘平差理论的有益拓展。
  • 表  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] 曾文宪, 刘泽邦, 方兴, 李玉兵.  通用EIV平差模型的线性化估计算法 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1284-1290. doi: 10.13203/j.whugis20200243
    [2] 陶叶青, 王坚, 刘超.  附不等式约束的地表沉降时间序列自回归EIV模型 . 武汉大学学报 ● 信息科学版, 2020, 45(9): 1455-1460. doi: 10.13203/j.whugis20180268
    [3] 谢建, 龙四春, 周璀.  不等式约束PEIV模型的最优性条件及SQP算法 . 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
    [4] 刘超, 王彬, 赵兴旺, 余学祥.  三维坐标转换的高斯-赫尔默特模型及其抗差解法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
    [5] 谢雪梅, 宋迎春, 肖兆兵.  附不等式约束平差模型的一种快速搜索算法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1349-1354. doi: 10.13203/j.whugis20160435
    [6] 王乐洋, 李海燕, 陈晓勇.  拟牛顿修正法解算不等式约束加权总体最小二乘问题 . 武汉大学学报 ● 信息科学版, 2018, 43(1): 127-132. doi: 10.13203/j.whugis20150333
    [7] 谢建, 龙四春, 李黎, 李博超.  不等式约束加权整体最小二乘的凝聚函数法 . 武汉大学学报 ● 信息科学版, 2018, 43(10): 1526-1530. doi: 10.13203/j.whugis20160507
    [8] 姚宜斌, 黄书华, 张良, 胡羽丰, 李国平.  求解三维坐标转换参数的整体最小二乘新方法 ! . 武汉大学学报 ● 信息科学版, 2015, 40(7): 853-857. doi: 10.13203/j.whugis20130302
    [9] 魏二虎, 殷志祥, 李广文, 李智强.  虚拟观测值法在三维坐标转换中的应用研究 . 武汉大学学报 ● 信息科学版, 2014, 39(2): 152-156. doi: 10.13203/j.whugis20120648
    [10] 姚宜斌, 黄书华, 孔建, 何军泉.  空间直线拟合的整体最小二乘算法 . 武汉大学学报 ● 信息科学版, 2014, 39(5): 571-574. doi: 10.13203/j.whugis20120104
    [11] 姚宜斌, 孔建.  顾及设计矩阵随机误差的最小二乘组合新解法 . 武汉大学学报 ● 信息科学版, 2014, 39(9): 1028-1032. doi: 10.13203/j.whugis20130030
    [12] 姚宜斌, 黄书华, 陈家君.  求解自回归模型参数的整体最小二乘新方法 . 武汉大学学报 ● 信息科学版, 2014, 39(12): 1463-1466.
    [13] 陈传法, 蔡乾广.  DEM快速构建的最小二乘配置法 . 武汉大学学报 ● 信息科学版, 2013, 38(1): 86-89.
    [14] 葛旭明, 伍吉仓.  三维基准转换的约束加权混合整体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2012, 37(2): 178-182.
    [15] 鲁铁定, 周世健.  总体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1351-1354.
    [16] 邱卫宁, 齐公玉, 田丰瑞.  整体最小二乘求解线性模型的改进算法 . 武汉大学学报 ● 信息科学版, 2010, 35(6): 708-710.
    [17] 孔建, 姚宜斌, 吴寒.  整体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(6): 711-714.
    [18] 宋迎春, 朱建军, 罗德仁, 曾联斌.  附非负约束平差模型的最小二乘估计 . 武汉大学学报 ● 信息科学版, 2008, 33(9): 907-909.
    [19] 鲁铁定, 陶本藻, 周世健.  基于整体最小二乘法的线性回归建模和解法 . 武汉大学学报 ● 信息科学版, 2008, 33(5): 504-507.
    [20] 张勤, 陶本藻.  基于同伦法的非线性最小二乘平差统一模型 . 武汉大学学报 ● 信息科学版, 2004, 29(8): 708-710.
  • 加载中
计量
  • 文章访问数:  788
  • HTML全文浏览量:  301
  • PDF下载量:  60
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-12-02
  • 刊出日期:  2021-09-18

不等式约束PEIV模型的经典最小二乘方法

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

    国家自然科学基金 41704007

    国家自然科学基金 41877283

    国家自然科学基金 42074016

    作者简介:

    谢建, 博士, 主要从事测量平差及数据处理研究。hsiejian841006@163.com

  • 中图分类号: P207

摘要: 不等式约束部分变量含误差(partial errors-in-variables, PEIV)模型目前主要采用线性化方法和非线性规划类算法, 前者计算效率较低, 后者基于最优化理论, 计算复杂, 未能与经典平差理论建立联系, 难以在测量实际中推广。在整体最小二乘准则下, 根据最优解的Kuhn-Tucker条件, 将不等式约束整体最小二乘解的计算转化为二次规划问题, 并提出改进的Jacobian迭代法求解二次规划。所提方法不需要对观测方程线性化, 与经典最小二乘法具有相同的形式, 易于编程实现。数值实例表明, 所提方法形式简洁, 具有良好的计算效率, 是经典最小二乘平差理论的有益拓展。

English Abstract

谢建, 龙四春, 周璀. 不等式约束PEIV模型的经典最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
引用本文: 谢建, 龙四春, 周璀. 不等式约束PEIV模型的经典最小二乘方法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
XIE Jian, LONG Si-chun, ZHOU Cui. Classical Least Squares Method for Inequality Constrained PEIV Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
Citation: XIE Jian, LONG Si-chun, ZHOU Cui. Classical Least Squares Method for Inequality Constrained PEIV Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(9): 1291-1297. doi: 10.13203/j.whugis20190196
  • 在大地测量坐标转换、曲线曲面拟合、摄影测量后方交会、大地测量反演等领域[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方法形式一致。

    • 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)

      式中,yey分别表示n维观测向量及其误差向量;β表示m维模型参数;Inn维单位矩阵;⊗是Kronecker积符号,其定义为AB=[aij×B], aij是矩阵A的第i行第j列元素;a是已知t维列向量;aea分别是a的真值和误差向量;h是已知的nm维常数向量;Bnm×t维常数矩阵;aB的定义见文献[15]; Gs×m维约束矩阵;ds维常数向量。令随机误差e=(eaT  eyT)T,其均值为0,且eaey互不相关,它们的方差分别为:

      $$ 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为单位权方差;Wn维观测值权矩阵;ω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] $, hin维列向量,Bin×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原理的算法。

    • 给定模型参数和系数阵中随机元素的近似值$ \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维二次规划问题,在大规模问题中,计算量会显著增加。

    • 令$ \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平差模型形式一致的迭代算法。

    • 由式(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形式。

    • 式(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)

      DijD的第i行第j列元素,lil的第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乘子,法方程的维数等于不等式约束的个数,当约束个数较少时,方程的计算量较小,因此,它更适应于不等式约束个数较少的情况。

    • ICPEIV模型的计算步骤如下:

      1)给定已知量ayhBω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)根据tn的大小以及使得计算量较小的原则,分别采用式(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方法,在比较计算效率时,应同时比较内外迭代次数和计算时间。

    • 算例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 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

      表 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法相近,进一步验证了本文方法内外迭代的复杂度明显较低。

    • 由§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 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

      表 2可知,当约束个数从算例1的11个变为3个以后,线性化方法和SQP方法需求解24维的QP问题,而本文方法只需求解4维的QP问题,且采用改进的Jacobian迭代法只需求解3维线性方程组。本文方法的内外迭代次数及计算时间都明显减少,验证了本文方法的计算效率随约束个数的减少而减少。因此,在约束个数较少的测绘应用中具有更好的适用性。3种方法在相同的迭代终止标准下得到的结果一致,但本文方法的计算效率比线性化方法和SQP方法更好。

    • 目前,求解ICPEIV模型的线性化方法计算效率偏低,转换成QP问题后需采用规划类算法求解。若采用非线性规划算法,如SQP方法,要求目标函数的梯度和Hessian矩阵,形式复杂,难以推广应用。针对上述问题,本文提出了求解ICPEIV模型的经典LS方法,将模型参数和系数阵元素分别估计,模型参数的估计归结为解QP问题,系数阵元素的估计则利用其与模型参数的关系计算。针对QP问题,以KT条件为基础,设计了改进的Jacobian迭代法将QP问题转化为线性方程的求解,它是求解等式约束LS模型的Lagrange乘子法的拓展。本文方法无需对观测方程线性化,也无需掌握复杂的规划类算法,形式简洁,与经典LS一致,易于推广应用。数值实验表明,本文方法具有良好的计算效率和广泛的适用性,能得到与线性化方法和SQP方法精度相当的WTLS解。特别是不等式约束的个数较少时,计算效率明显提高。本文方法将ICPEIV模型及相应的WTLS理论纳入经典LS理论的框架,是对经典LS理论的有益拓展和补充。

参考文献 (18)

目录

    /

    返回文章
    返回