留言板

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

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

基于变量投影法的自回归模型方差分量估计

吕志鹏 隋立芬

吕志鹏, 隋立芬. 基于变量投影法的自回归模型方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
引用本文: 吕志鹏, 隋立芬. 基于变量投影法的自回归模型方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
LÜ Zhipeng, SUI Lifen. Variance Component Estimation of Autoregressive Model Based on Variable Projection Method[J]. Geomatics and Information Science of Wuhan University, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
Citation: LÜ Zhipeng, SUI Lifen. Variance Component Estimation of Autoregressive Model Based on Variable Projection Method[J]. Geomatics and Information Science of Wuhan University, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352

基于变量投影法的自回归模型方差分量估计

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

国家自然科学基金 41674016

国家自然科学基金 41274016

详细信息
    作者简介:

    吕志鹏, 博士生, 主要从事空间大地测量数据处理研究。lvzhipeng1989@qq.com

  • 中图分类号: P207

Variance Component Estimation of Autoregressive Model Based on Variable Projection Method

Funds: 

The National Natural Science Foundation of China 41674016

The National Natural Science Foundation of China 41274016

More Information
    Author Bio:

    LÜ Zhipeng, PhD candidate, specializes in theories and methods of spatial geodetic data processing. E-mail:lvzhipeng1989@qq.com

  • 摘要: 在自回归(autoregressive,AR)模型中,系数矩阵与观测向量中的随机误差同源。针对AR模型平差时观测权阵分配不合理、随机模型不准确的情况,采用变量投影法提取系数矩阵和观测向量构成的增广矩阵中的随机量,将变量误差(errors-in-variables,EIV)模型转化为非线性高斯-赫尔默特(Gauss-Helmert,GH)模型,利用非线性最小二乘理论得到一种结构总体最小二乘(structural total least squares,STLS)算法,并与最小二乘方差分量估计(least squares variance component estimation,LS-VCE)相结合推导出STLS问题的一种方差分量估计算法,将其应用到AR模型的方差分量估计。通过实测算例对算法有效性进行了验证,取得了与已有方法一致的结果。该算法观测权阵的构造十分简洁,同时也可用于协方差分量的估计。
  • 图  1  方差分量估值变化图

    Figure  1.  Changes of Estimated Variance Factors

    表  1  沉降观测数据/mm

    Table  1.   Settlement Observation Data/mm

    序号 高程 序号 高程 序号 高程 序号 高程
    1 26.33 10 25.46 19 28.09 28 26.57
    2 26.27 11 26.12 20 26.78 29 28.36
    3 26.43 12 27.28 21 28.66 30 27.94
    4 25.56 13 26.67 22 26.75 31 26.81
    5 26.82 14 27.95 23 27.24 32 28.50
    6 26.56 15 26.74 24 28.02 33 27.68
    7 25.93 16 27.53 25 26.81 34 26.57
    8 26.43 17 25.31 26 28.50 35 28.36
    9 26.52 18 26.90 27 27.68 36 27.94
    下载: 导出CSV

    表  2  不同算法得到的参数估值及其精度

    Table  2.   Estimates of Model Parameters and Their Deviations from Different Methods

    项目 参数 LS估计 PEIV模型解法 Amiri-Simkooei改进算法 本文算法
    估值 ${{\hat \xi }_1}$ 0.635 1 1.179 1 1.075 5 1.075 5
    ${{\hat \xi }_2}$ 0.327 8 0.041 9 0.005 5 0.005 5
    ${{\hat \xi }_3}$ 0.041 1 -0.214 4 -0.075 2 -0.075 2
    标准差 ${{\hat \sigma }_{{\xi _1}}}$ 0.150 7 0.112 5 0.034 7 0.034 7
    ${{\hat \sigma }_{{\xi _2}}}$ 0.133 1 0.097 4 0.032 7 0.032 7
    ${{\hat \sigma }_{{\xi _3}}}$ 0.145 5 0.095 3 0.031 5 0.031 5
    下载: 导出CSV

    表  3  不同算法得到的方差分量估值及其精度

    Table  3.   Estimates of Variance Factors and Their Deviations from Different Methods

    项目 参数 Amiri-Simkooei改进算法 本文算法
    估值 $\hat \sigma _f^2$ 0.178 1 0.178 1
    $\hat \sigma _s^2$ 1.015 6 1.015 6
    $\hat \sigma _t^2$ 0.081 3 0.081 3
    标准差 ${{\hat \sigma }_{\sigma _f^2}}$ 0.083 5 0.083 5
    ${{\hat \sigma }_{\sigma _s^2}}$ 0.423 2 0.423 2
    ${{\hat \sigma }_{\sigma _t^2}}$ 0.039 9 0.039 9
    下载: 导出CSV

    表  4  不同算法得到的单位权方差估值及其精度

    Table  4.   Estimates of Unit Weight Variance and Their Deviations from Different Methods

    项目 参数 LS估计 PEIV模型解法
    估值 ${\hat \sigma _0^2}$ 0.804 7 0.643 4
    标准差 ${{\hat \sigma }_{\hat \sigma _0^2}}$ 0.167 2 0.106 9
    下载: 导出CSV
  • [1] Van Loan S, Van Dewalle J. The Total Least-Squares Problem:Computational Aspects and Analysis[M]. Philadelphia:Society for Industrial and Applied Mathematics, 1991
    [2] 刘经南, 曾文宪, 徐培亮.整体最小二乘估计的研究进展[J].武汉大学学报·信息科学版, 2013, 38(5):505-512 http://ch.whu.edu.cn/CN/abstract/abstract2641.shtml

    Liu Jingnan, Zeng Wenxian, Xu Peiliang. Overview of Total Least Squares Methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5):505-512 http://ch.whu.edu.cn/CN/abstract/abstract2641.shtml
    [3] 王乐洋, 许才军.总体最小二乘研究进展[J].武汉大学学报·信息科学版, 2013, 38(7):850-856 http://ch.whu.edu.cn/CN/abstract/abstract2703.shtml

    Wang Leyang, Xu Caijun. Progress in Total Least Squares[J]. Geomatics and Information Science of Wuhan University, 2013, 38(7):850-856 http://ch.whu.edu.cn/CN/abstract/abstract2703.shtml
    [4] Fang Xing. A Structured and Constrained Total Least-Squares Solution with Cross-Covariances[J]. Studia Geophysica et Geodaetica, 2014, 58(1):1-16 doi:  10.1007/s11200-012-0671-z
    [5] Rosen J B, Park H, Glick J. Total Least Norm Formulation and Solution for Structured Problems[J]. SIAM Journal on Matrix Analysis & Applications, 1996, 17(1):110-126 http://www.researchgate.net/publication/239061482_Total_Least_Norm_Formulation_and_Solution_for_Structured_Problems
    [6] Van Huffel S, Park H, Rosen J B. Formulation and Solution of Structured Total Least Norm Problems for Parameter Estimation[J]. IEEE Transactions on Signal Processing, 1996, 44(10):2464-2474 doi:  10.1109/78.539031
    [7] Markovsky I, Huffel S V. On Weighted Structured Total Least Squares[C]. International Conference on Large-Scale Scientific Computing, Sozopol, Bulgaria, 2005
    [8] Xu Peiliang, Liu Jiangnan, Shi Chuang. 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
    [9] 王乐洋, 余航, 陈晓勇. Partial EIV模型的解法[J].测绘学报, 2016, 45(1):22-29 doi:  10.11947/j.AGCS.2016.20140560

    Wang Leyang, Yu Hang, Chen Xiaoyong. An Algorithm for Partial EIV Model[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1):22-29 doi:  10.11947/j.AGCS.2016.20140560
    [10] 王乐洋, 许光煜, 温贵森.一种相关观测的Partial EIV模型求解方法[J].测绘学报, 2017, 46(8):44-53 http://d.old.wanfangdata.com.cn/Periodical/chxb201708007

    Wang Leyang, Xu Guangyu, Wen Guisen. A Method for Partial EIV Model with Correlated Observations[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(8):44-53 http://d.old.wanfangdata.com.cn/Periodical/chxb201708007
    [11] Fang Xing. 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
    [12] Amiri-Simkooei A R, Jazaeri S. Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]. Journal of Geodetic Science, 2012, 2(2):113-124 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.2478/v10156-011-0036-5
    [13] 姚宜斌, 黄书华, 陈家君.求解自回归模型参数的整体最小二乘新方法[J].武汉大学学报·信息科学版, 2014, 39(12):1463-1466 http://ch.whu.edu.cn/CN/abstract/abstract3142.shtml

    Yao Yibin, Huang Shuhua, Chen Jiajun. A New Method of TLS to Solving the Autoregressive Model Parameter[J]. Geomatics and Information Science of Wuhan University, 2014, 39(12):1463-1466 http://ch.whu.edu.cn/CN/abstract/abstract3142.shtml
    [14] 姚宜斌, 熊朝晖, 张豹, 等.顾及设计矩阵误差的AR模型新解法[J].测绘学报, 2017, 46(11):1795-1801 doi:  10.11947/j.AGCS.2017.20170004

    Yao Yibin, Xiong Zhaohui, Zhang Bao, et al.A New Method to Solving AR Model Parameters Considering Random Errors of Design Matrix[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(11):1795-1801 doi:  10.11947/j.AGCS.2017.20170004
    [15] Teunissen P J G, Amiri-Simkooei A R. Least-Squares Variance Component Estimation[J]. Journal of Geodesy, 2008, 82(2):65-82 doi:  10.1007/s00190-007-0157-x
    [16] Amiri-Simkooei A R. Least-Squares Variance Component Estimation: Theory and GPS Applications[D]. Delft: Delft University of Technology, 2007
    [17] Amiri-Simkooei A R. Application of Least Squares Variance Component Estimation to Errors-In-Variables Models[J]. Journal of Geodesy, 2013, 87(10-12):935-944 doi:  10.1007/s00190-013-0658-8
    [18] Mahboub V. Variance Component Estimation in Errors-In-Variables Models and a Rigorous Total Least-Squares Approach[J]. Studia Geophysica et Geodaetica, 2014, 58(1):17-40 doi:  10.1007/s11200-013-1150-x
    [19] Xu Peiling, Liu Jingnan. Variance Components in Errors-In-Variables Models:Estimability, Stability and Bias Analysis[J]. Journal of Geodesy, 2014, 88(8):719-734 doi:  10.1007/s00190-014-0717-9
    [20] Wang Leyang, Xu Guangyu. Variance Component Estimation for Partial Errors-In-Variables Models[J]. Studia Geophysica et Geodaetica, 2016, 60(1):35-55 doi:  10.1007/s11200-014-0975-2
    [21] 王乐洋, 许光煜, 陈晓勇.附有相对权比的PEIV模型总体最小二乘平差[J].武汉大学学报·信息科学版, 2017, 42(6):857-863 http://ch.whu.edu.cn/CN/abstract/abstract5764.shtml

    Wang Leyang, Xu Guangyu, Chen Xiaoyong. Total Least Squares Adjustment of Partial Errors-In-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6):857-863 http://ch.whu.edu.cn/CN/abstract/abstract5764.shtml
    [22] 王乐洋, 温贵森.相关观测PEIV模型的最小二乘方差分量估计[J].测绘地理信息, 2018, 43(1):8-14 http://d.old.wanfangdata.com.cn/Periodical/chxxygc201801002

    Wang Leyang, Wen Guisen. Least Squares Variance Components Estimation of PEIV Model with Correlated Observations[J]. Journal of Geomatics, 2018, 43(1):8-14 http://d.old.wanfangdata.com.cn/Periodical/chxxygc201801002
    [23] 王乐洋, 温贵森. Partial EIV模型的非负最小二乘方差分量估计[J].测绘学报, 2017, 46(7):857-865 http://d.old.wanfangdata.com.cn/Periodical/chxb201707008

    Wang Leyang, Wen Guisen. Non-negative Least Squares Variance Estimation of Partial EIV Model[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(7):857-865 http://d.old.wanfangdata.com.cn/Periodical/chxb201707008
    [24] Wang Leyang, Zhao Yingwen. Unscented Transformation with Scaled Symmetric Sampling Strategy for Precision Estimation of Total Least Squares[J]. Studia Geophysica et Geodaetica, 2017, 61(3):385-411 doi:  10.1007/s11200-016-1113-0
    [25] Wang Leyang, Zhao Yingwen. The Scaled Unscented Transformation for Nonlinear Error Propagation Accuracy, Sensitivity and Applications[J]. Journal of Surveying Engineering, 2018, 144(1):04017022 doi:  10.1061/(ASCE)SU.1943-5428.0000243
    [26] Wang Leyang, Zhao Yingwen. Second Order Approximating Function Method for Precision Estimation of Total Least Squares[J]. Journal of Surveying Engineering, 2019, 145(1):04018011 doi:  10.1061/(ASCE)SU.1943-5428.0000266
    [27] 王新洲, 陶本藻, 邱卫宁, 等.高等测量平差[M].北京:测绘出版社, 2006

    Wang Xinzhou, Tao Benzao, Qiu Weining, et al. Advanced Surveying Adjustment[M]. Beijing:Surveying and Mapping Press, 2006
    [28] Teunissen P J G. Nonlinear Least-Squares[J]. Manuscripta Geodaetica, 1990, 15(3):137-150 http://d.old.wanfangdata.com.cn/Periodical/gpxygpfx201411036
  • [1] 王乐洋, 谷旺旺, 赵雄, 许光煜, 高华.  利用偏差改正的方差分量估计方法确定联合反演相对权比 . 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
    [2] 王乐洋, 邹传义.  PEIV模型参数估计理论及其应用研究进展 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1273-1283, 1297. doi: 10.13203/j.whugis20200312
    [3] 汪友军, 胡俊, 刘计洪, 孙倩.  融合InSAR和GNSS的三维形变监测:利用方差分量估计的改进SISTEM方法 . 武汉大学学报 ● 信息科学版, 2021, 46(10): 1598-1608. doi: 10.13203/j.whugis20210113
    [4] 吕志鹏, 隋立芬.  基于变量投影的结构总体最小二乘算法 . 武汉大学学报 ● 信息科学版, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
    [5] 王乐洋, 孙坚强.  总体最小二乘回归预测模型的方差分量估计 . 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
    [6] 许才军, 邓长勇, 周力璇.  利用方差分量估计的地震同震滑动分布反演 . 武汉大学学报 ● 信息科学版, 2016, 41(1): 37-44. doi: 10.13203/j.whugis20150500
    [7] 孙文, 吴晓平, 王庆宾, 刘晓刚, 朱志大.  基于方差分量估计的正则化配置法及其在多源重力数据融合中的应用 . 武汉大学学报 ● 信息科学版, 2016, 41(8): 1087-1092. doi: 10.13203/j.whugis20140159
    [8] 郭春喜, 聂建亮, 王斌, 蒋光伟.  利用自适应最小二乘配置的GPS水准与重力似大地水准面的拟合 . 武汉大学学报 ● 信息科学版, 2013, 38(1): 44-47.
    [9] 刘志平, 张书毕.  方差-协方差分量估计的概括平差因子法 . 武汉大学学报 ● 信息科学版, 2013, 38(8): 925-929.
    [10] 徐天河, 居向明.  基于方差分量估计的CHAMP重力场恢复方法 . 武汉大学学报 ● 信息科学版, 2007, 32(3): 242-245.
    [11] 王新洲, 游扬声, 刘星, 史文中.  GIS叠置图层方差分量的极大似然估计 . 武汉大学学报 ● 信息科学版, 2005, 30(10): 892-896.
    [12] 杨元喜, 徐天河.  基于移动开窗法协方差估计和方差分量估计的自适应滤波 . 武汉大学学报 ● 信息科学版, 2003, 28(6): 714-718.
    [13] 王仲锋.  导线网方差分量估计的综合研究 . 武汉大学学报 ● 信息科学版, 2001, 26(2): 112-117.
    [14] 盛乐山.  三维网平差的方差分量估计 . 武汉大学学报 ● 信息科学版, 1990, 15(1): 91-98.
    [15] 许才军.  卫星网与地面网联合平差中随机模型的估计 . 武汉大学学报 ● 信息科学版, 1989, 14(3): 35-45.
    [16] 吴晓清.  权因子法方差分量估计 . 武汉大学学报 ● 信息科学版, 1989, 14(4): 8-15.
    [17] 张万鹏.  方差分量估计理论及其在边角网平差中的应用 . 武汉大学学报 ● 信息科学版, 1988, 13(1): 9-21.
    [18] 冯延明.  水准测量往返测较差方差之估计 . 武汉大学学报 ● 信息科学版, 1987, 12(4): 20-27.
    [19] 张宝庭.  方差-协方差-分量估计 . 武汉大学学报 ● 信息科学版, 1982, 7(1): 42-57.
    [20] 雷雨, 赵丹宁.  附有时变参数的周期项模型在极移预报中的应用 . 武汉大学学报 ● 信息科学版, 0, 0(0): 0-0. doi: 10.13203/j.whugis20200007
  • 加载中
图(1) / 表(4)
计量
  • 文章访问数:  870
  • HTML全文浏览量:  166
  • PDF下载量:  62
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-28
  • 刊出日期:  2020-02-05

基于变量投影法的自回归模型方差分量估计

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

    国家自然科学基金 41674016

    国家自然科学基金 41274016

    作者简介:

    吕志鹏, 博士生, 主要从事空间大地测量数据处理研究。lvzhipeng1989@qq.com

  • 中图分类号: P207

摘要: 在自回归(autoregressive,AR)模型中,系数矩阵与观测向量中的随机误差同源。针对AR模型平差时观测权阵分配不合理、随机模型不准确的情况,采用变量投影法提取系数矩阵和观测向量构成的增广矩阵中的随机量,将变量误差(errors-in-variables,EIV)模型转化为非线性高斯-赫尔默特(Gauss-Helmert,GH)模型,利用非线性最小二乘理论得到一种结构总体最小二乘(structural total least squares,STLS)算法,并与最小二乘方差分量估计(least squares variance component estimation,LS-VCE)相结合推导出STLS问题的一种方差分量估计算法,将其应用到AR模型的方差分量估计。通过实测算例对算法有效性进行了验证,取得了与已有方法一致的结果。该算法观测权阵的构造十分简洁,同时也可用于协方差分量的估计。

English Abstract

吕志鹏, 隋立芬. 基于变量投影法的自回归模型方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
引用本文: 吕志鹏, 隋立芬. 基于变量投影法的自回归模型方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
LÜ Zhipeng, SUI Lifen. Variance Component Estimation of Autoregressive Model Based on Variable Projection Method[J]. Geomatics and Information Science of Wuhan University, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
Citation: LÜ Zhipeng, SUI Lifen. Variance Component Estimation of Autoregressive Model Based on Variable Projection Method[J]. Geomatics and Information Science of Wuhan University, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
  • 变量误差(errors-in-variables,EIV)模型及其总体最小二乘(total least squares,TLS)估计研究已逾百年,取得了一系列重要成果,并在科学和工程中得到了广泛应用[1]。测绘领域对TLS估计的应用集中于线性回归、基准变换和摄影测量空间后方交会等问题[2-3],在这些应用中,EIV模型的系数矩阵存在结构特征,即系数矩阵存在固定量和随机量以及不同位置的随机量之间存在相关性。对于系数矩阵存在固定列的情况,可以采用基于矩阵QR分解的混合LS-TLS(least squares-total least squares)估计求解[1]。而更一般的情况,可以采用各种结构总体最小二乘(structural total least squares, STLS)估计[4-7]求解。Xu等[8]提出了Partial EIV (PEIV)模型,通过构造虚拟观测值推导出一种STLS问题的算法,王乐洋等[9-10]对上述算法进行了改进,提高了算法的计算效率,且适用于相关观测PEIV模型的参数估计。也可以通过构造秩亏的观测权阵,利用加权总体最小二乘(weighted total least squares,WTLS)估计[11-12]求解STLS问题。

    自回归(autoregressive, AR)模型是一类重要的时间序列模型,它在科学研究和工程实践中有广泛的应用,如形变分析与预报、导航定位、极移预测、陀螺漂移分析等[13-14]。AR模型的求解实际上是一个STLS问题。文献[13]推导了一种非线性的AR模型TLS迭代算法,但难以采用协方差传播定律进行精度评定;文献[14]通过引入虚拟观测值,将TLS问题转化为经典的最小二乘问题,采用间接平差方法进行参数求解。AR模型观测值可能来源于不同类型的观测或平差,甚至来自不同期观测,因而具有不同的精度。已有算法[10, 13-14]认为,系数矩阵和观测向量中的随机量拥有相同的验前单位权方差并不合适,需要确定合理的观测权矩阵以确保参数的最优估计,因此需要进行方差分量估计,对随机模型进行修正。

    对于高斯-马尔可夫(Gauss-Markov,GM)模型,文献[15]提出了最小二乘方差分量估计(least squares variance component estimation, LS-VCE),文献[16]将LS-VCE应用到GPS的数据处理中,说明了这种方法的有效性。对于EIV模型,文献[17]给出了EIV模型的LS-VCE算法;文献[18]基于LS-VCE原理提出了一种严密的TLS估计算法,实现了参数和方差分量同步估计,但这种算法并不适用于包含固定列的系数矩阵,因而作者进一步给出了两种替代算法;文献[19]推导了EIV模型的MINQUE(minimum norm quadratic unbiased estimator)算法;文献[20]研究了PEIV模型的Helmert方差分量估计算法;文献[21]以PEIV模型为基础,推导了附有相对权比的TLS平差算法,并给出了两种相对权比的确定方法。上述各种算法并未考虑系数矩阵和观测向量的相关性,故不适用于协方差分量的估计,不能直接应用于AR模型的方差分量估计。文献[22]进一步分析了相关观测PEIV模型的LS-VCE算法;文献[23]基于非负最小二乘理论给出了PEIV模型的非负LS-VCE算法。文献[20-23]均是将系数矩阵和观测向量中的随机量作为两类数据在PEIV模型下进行方差分量估计,需要进行简单的多类观测数据扩展,否则也不能直接应用于AR模型的方差分量估计。

    为此,本文提出了一种STLS问题的LS-VCE算法,并将其应用于AR模型的方差分量估计。该方法通过变量投影法[7]将系数矩阵和观测向量构成的增广矩阵展开成仿射矩阵形式,将EIV模型在参数近似值和增广矩阵中随机量近似值处展开成非线性GH(Gauss-Helmert)模型,由此将TLS问题转化为LS问题,并采用附有系统参数的条件平差法进行迭代求解。同时根据LS-VCE原理对STLS问题的随机模型进行分析,计算随机量的方差和协方差分量估值,在迭代过程中更新观测权阵,达到修正参数估值的效果,并通过实测算例验证了本文算法的有效性。

    • EIV模型考虑了系数矩阵和观测向量误差,其函数模型如下:

      $$ \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} - \left( {A - {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right) \cdot \mathit{\boldsymbol{x}} = 0 $$ (1)

      式中,$\mathit{\boldsymbol{y}} \in {{\bf{R}}^{m \times 1}}$表示观测向量;${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}$表示观测向量的误差向量;$\mathit{\boldsymbol{A}} \in {{\bf{R}}^{m \times n}}$表示系数矩阵;${{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}$表示系数矩阵的误差矩阵;$\mathit{\boldsymbol{x}} \in {{\bf{R}}^{n \times 1}}$表示参数向量。对于STLS问题,其系数矩阵和观测向量构成的增广矩阵$\mathit{\boldsymbol{C}} = \left[ {\mathit{\boldsymbol{A}}\;\;\mathit{\boldsymbol{y}}} \right]$具有某种结构,其误差矩阵${\rm{\Delta }}\mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \end{array}} \right]$应具有相同结构。假设C仅存在某种仿射结构,即不考虑随机量之间存在非线性关系的情况,对于随机量之间存在非线性关系时,可以先将其进行线性化然后再求解。将C中独立随机量组成的向量设为$\mathit{\boldsymbol{p}} \in {{\bf{R}}^{{n_p}}}$,${{n_p}}$为独立随机量个数。故C的仿射函数形式如下[7]

      $$ C = \mathit{\boldsymbol{S}}\left( \mathit{\boldsymbol{p}} \right) = {\mathit{\boldsymbol{S}}_0} + \sum\limits_{j = 1}^{{n_p}} {{\mathit{\boldsymbol{S}}_j}} {p_j}\;\;\;\;\;\left( {{\mathit{\boldsymbol{S}}_0}, {\mathit{\boldsymbol{S}}_j} \in {{\rm{R}}^{m \times \left( {n + 1} \right)}}} \right) $$ (2)

      式中,${\mathit{\boldsymbol{S}}_0}$为固定量仿射结构矩阵,${\mathit{\boldsymbol{S}}_0}$中的非零元素为C中的固定量;${\mathit{\boldsymbol{S}}_j} = \left( {j = 1, 2 \cdots {n_p}} \right)$为随机量${p_j} = \left( {j = 1, 2 \cdots {n_p}} \right)$对应的仿射结构矩阵,它描述了${p_j}$在C中的分布。故可将STLS问题的数学模型变形为:

      $$ \left\{ {\begin{array}{*{20}{l}} {f\left( {{\rm{\Delta }}\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{p}} - {\rm{\Delta }}\mathit{\boldsymbol{p}}} \right){\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = 0}\\ {{\rm{\Delta }}\mathit{\boldsymbol{p}} \sim \left( {0, {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}} \right)} \end{array}} \right. $$ (3)

      式中,${\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{x}}\\ { - 1} \end{array}} \right]$为x的增广向量;${{\rm{\Delta }}\mathit{\boldsymbol{p}}}$为p的误差向量;QPP的协方差阵。对式(3)中的函数模型进一步变形得:

      $$ \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{p}} - {\rm{\Delta }}\mathit{\boldsymbol{p}}} \right){\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{x}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} = 0 $$ (4)

      其中,

      $$ \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{x}} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}_1}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}}&{{\mathit{\boldsymbol{S}}_2}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}}& \ldots &{{\mathit{\boldsymbol{S}}_{{n_\mathit{\boldsymbol{p}}}}}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}} \end{array}} \right] = \left( {\mathit{\boldsymbol{x}}_{{\rm{ext}}}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{\psi }} $$ (5)

      式中,$\mathit{\boldsymbol{\psi }} = \left[ {\begin{array}{*{20}{c}} {{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_1}} \right)}&{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_2}} \right)}& \ldots &{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_{{n_\mathit{\boldsymbol{p}}}}}} \right)}\end{array}} \right]$;$ \otimes $为克罗内克乘积算子;Im表示m × m的单位矩阵;vec(·)为矩阵向量化算子。

    • EIV模型是一种非线性模型,故将式(4)在参数近似值${{{\hat x}^{i - 1}}}$和增广矩阵C中随机量误差近似值${\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^{i - 1}}$处进行一阶泰勒展开,得:

      $$ f\left( {{\rm{\Delta }}\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{x}}} \right) \approx \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} - {{\hat A}^{i - 1}}\mathit{\boldsymbol{x}} + \hat y_c^{i - 1} = 0 $$ (6)

      式中,${{\hat A}^{i - 1}} = \mathit{\boldsymbol{A}} - \hat E_\mathit{\boldsymbol{A}}^{i - 1}$ ($\hat E_\mathit{\boldsymbol{A}}^{i - 1}$的构造方法见式(2));$\hat y_c^{i - 1} = \mathit{\boldsymbol{y}} - \hat E_\mathit{\boldsymbol{A}}^{i - 1}{{\hat x}^{i - 1}}$。式(6)可视为经典平差中附加系统参数的条件平差模型,可采用经典平差原理进行求解。令${\rm{\Delta }}\hat y_c^{i - 1} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}}$,故${\rm{\Delta }}\hat y_c^{i - 1} \sim \left( {0, \;\mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right)} \right)$,将式(6)变为经典平差中的参数平差模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {{\rm{\Delta }}\hat y_c^{i - 1} = {{\hat A}^{i - 1}}\mathit{\boldsymbol{x}} - \hat y_c^{i - 1}}\\ {{\rm{\Delta }}\hat y_c^{i - 1} \sim \left( {0, \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right)} \right)} \end{array}} \right. $$ (7)

      根据参数平差法得:

      $$ {{\hat x}^i} = {[{{\hat A}^{i - 1}}^{\rm{T}}{({\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}})^{ - 1}}{{\hat A}^{i - 1}}]^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}{({\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}})^{ - 1}}\hat y_c^{i - 1} $$ (8)

      其中,

      $$ {\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}} = {\mathit{\boldsymbol{Q}}_{{\rm{\Delta }}\hat y_c^{i - 1}}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right) $$ (9)

      同时,可得$\hat y_c^{i - 1}$的估值$\hat y_c^{i}$及其误差的估值$\;{\mathit{\boldsymbol{V}}_{\hat y_c^i}}$为:

      $$ \left\{ \begin{array}{l} \hat y_c^i = {\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}}\hat y_c^{i - 1} = {{\hat A}^{i - 1}}{{\hat x}^i} = \left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}}_A^{i - 1}} \right){{\hat x}^i}\\ {\mathit{\boldsymbol{V}}_{\hat y_c^i}} = \mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}^ \bot \hat y_c^{i - 1} = \hat y_c^{i - 1} - \hat y_c^i = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{{\rm A}}}{{\hat x}^i} = \hat e_\mathit{\boldsymbol{y}}^{i - 1} - \mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}^{i - 1}{{\hat x}^i} \end{array} \right. $$ (10)

      式中,${\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}} = {{\hat A}^{i - 1}}{({{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}{{\hat A}^{i - 1}})^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}$;$\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}^ \bot = {\mathit{\boldsymbol{I}}_m} - {{\hat A}^{i - 1}}{({{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}{{\hat A}^{i - 1}})^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}$,均为正交投影矩阵。

      已知x估值${{\hat x}^i}$后,则式(6)变为经典平差中的条件平差模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{G}}\left( {{{\hat x}^i}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} + \left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\hat x}^i}} \right) = 0}\\ {{\rm{\Delta }}\mathit{\boldsymbol{p}} \sim \left( {0, {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}} \right)} \end{array}} \right. $$ (11)

      根据条件平差法得:

      $$ {\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^i} = - {\mathit{\boldsymbol{Q}}_p}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^i}} \right)\mathit{\boldsymbol{Q}}_{\hat y_c^i}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\hat x}^i}} \right) $$ (12)

      式中,${\mathit{\boldsymbol{Q}}_{\hat y_c^i}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^i}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^i}} \right)$。已知${\rm{\Delta }}\mathit{\boldsymbol{p}}$的估值${\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^i}$后,由式(2)得:

      $${\rm{\Delta }}{\hat C^i} = \left[ {\begin{array}{*{20}{c}} {\hat E_\mathit{\boldsymbol{A}}^i}&{\hat e_\mathit{\boldsymbol{y}}^i} \end{array}} \right] = \mathit{\boldsymbol{S}}\left( {{\rm{\Delta }}{{\hat p}^i}} \right) = \sum\limits_{j = 1}^{{n_\mathit{\boldsymbol{p}}}} {{\mathit{\boldsymbol{S}}_j}} {\rm{\Delta }}\hat p_j^i$$ (13)

      则${{\hat A}^i} = \mathit{\boldsymbol{A}} - \hat E_\mathit{\boldsymbol{A}}^i$。

      根据最小二乘原理,上述STLS估计的单位权方差可以表示为:

      $${{\hat \sigma }^2} = {\rm{\Delta }}{{\hat p}^i}^{\rm{T}}\mathit{\boldsymbol{P}}{\rm{\Delta }}{{\hat p}^i}/r$$ (14)

      式中,P为权矩阵;r为自由度,即r = m - n

      对于参数估值的精度,由参数平差模型式(7)和参数估值迭代表达式(8),根据协因数传播定律,可得:

      $${\mathit{\boldsymbol{Q}}_{\hat x\hat x}} = {[{{\hat A}^i}^{\rm{T}}{(\mathit{\boldsymbol{Q}}_{\hat y}^i)^{ - 1}}{{\hat A}^i}]^{ - 1}}$$ (15)

      上述推导过程仅考虑到了一阶误差的影响,对于更高阶误差所引起的参数估计以及精度评定偏差参见文献[24-26]。

    • LS-VCE由Teunissen等提出[15-16],该方法采用最小二乘准则进行方差分量估计,在线性GM模型假设下,可以得到方差分量的线性最优无偏估计。各种方差分量估计的目的均是对随机模型进行修正,区别在于所采用的估计原则以及对统计信息的依赖程度等方面的差异,并没有哪一种方差分量估计方法的估计特性较其他方法绝对占优。LS-VCE的优点在于,将方差分量估计问题纳入到最小二乘平差理论体系中,有利于方差分量估计的精度评定以及结合最小二乘理论进行其他方面的拓展研究。

      式(3)中的随机模型可以分解为更一般的分量形式:

      $${\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sum\limits_{k = 1}^d {{\sigma _k}{\mathit{\boldsymbol{Q}}_k}} $$ (16)

      式中,${{\sigma _k}}(k=1, 2…d)$表示随机量p的方差或协方差分量;d为方差分量的个数;Q0为随机模型中的已知部分;Qk$(k=1, 2…d)$为${{\sigma _k}}$对应的协因数矩阵且为线性无关的对称矩阵,并要求式(16)中的线性组合Qp为正定矩阵。

      在实际平差中,系数矩阵和观测向量中的随机量可能具有不同的方差分量,如果考虑到随机量之间的相关性,也可能具有不同的协方差分量,此时采用唯一的单位权方差进行平差并不合理。Teunissen提出的LS-VCE适用于GM模型,并可以进一步推广到EIV模型。将STLS问题转换为LS问题,在式(8)和式(12)的迭代计算过程中,Qp为不断修正的协方差阵,可以通过方差分量估计确定未知的方差和协方差分量,以修正下一轮迭代中所给定的权,从而实现优化参数估值的目的。假设STLS迭代结束后得到的参数估值和随机量误差估值分别为${\hat x}$和${\rm{\Delta }}\hat p$,将它们代入式(7)得:

      $$ \left\{ {\begin{array}{*{20}{l}} {{\rm{\Delta }}{{\hat y}_c} = \hat A\mathit{\boldsymbol{x}} - {{\hat y}_c}}\\ {{\rm{\Delta }}{{\hat y}_c} \sim \left( {0, \sigma _0^2\mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right)} \right)} \end{array}} \right. $$ (17)

      式中,$\hat A = \mathit{\boldsymbol{A}} - {{\hat E}_A}$;${\rm{\Delta }}{{\hat y}_c} = \mathit{\boldsymbol{G}}\left( {\hat x} \right){\rm{\Delta }}\mathit{\boldsymbol{p}}$;${{\hat y}_c} = \mathit{\boldsymbol{y}} - {{\hat E}_A}\hat x$。

      式(17)可视为GM模型,其中${{\hat y}_c}$、${\rm{\Delta }}{{\hat y}_c}$、${\hat A}$分别为GM模型的观测向量、观测向量的误差向量和系数矩阵,因而可以直接将LS-VCE应用于式(17)中的平差问题。

      将式(16)代入式(9)使${{\hat y}_c}$的协方差阵变形为:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{_{{{\hat y}_c}}}} = {\mathit{\boldsymbol{Q}}_{{\rm{\Delta }}{{\hat y}_c}}} = \mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) = }\\ {\mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_0}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) + \mathit{\boldsymbol{G}}\left( {\hat x} \right)\left( {\sum\limits_{k = 1}^d {{\sigma _k}{\mathit{\boldsymbol{Q}}_k}} } \right){\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) = }\\ {\mathit{\boldsymbol{Q}}_0^\mathit{\boldsymbol{G}} + \sum\limits_{k = 1}^d {{\sigma _k}} \mathit{\boldsymbol{Q}}_k^\mathit{\boldsymbol{G}}} \end{array} $$ (18)

      设矩阵${{\hat B}^{\rm{T}}}$为${\hat A}$的零空间,即${{\hat B}^{\rm{T}}}\hat A = 0$,${{\hat A}^{\rm{T}}}\hat B = 0$,同时引入向量t。由文献[15-16]可知,将${{\hat B}^{\rm{T}}}$左乘式(17)中的函数模型可得${{\hat B}^{\rm{T}}}$、t和${{\hat y}_c}$的关系如下:

      $$ \mathit{\boldsymbol{t}} = {{\hat B}^{\rm{T}}}{{\hat y}_c}, \;{\rm{E}}\left( \mathit{\boldsymbol{t}} \right) = 0 $$ (19)

      式中,E(·)表示数学期望。由协因数传播定律得:

      $$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{t}}} = {\rm{E}}\left( {\mathit{\boldsymbol{t}}{\mathit{\boldsymbol{t}}^{\rm{T}}}} \right) = {{\hat B}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{_{{{\hat y}_c}}}}\hat B $$ (20)

      式(20)即为方差分量估计的函数模型。为了消除式(20)中的重复观测方程,将vh(·)算子(由对称矩阵所对应的上三角矩阵或者下三角矩阵按列排列构成列向量)应用于上式,得:

      $$ {\rm{E}}\left( {{\mathit{\boldsymbol{y}}_{{\rm{vh}}}}} \right) = {\mathit{\boldsymbol{A}}_{{\rm{vh}}}}\mathit{\boldsymbol{\sigma }} $$ (21)

      式中,yvhAvh为对式(20)进行vh(·)操作后的向量和矩阵;$\mathit{\boldsymbol{\sigma }}$为未知的方差和协方差分量组成的向量。由最小二乘准则并经过一系列推导和简化[15-16]得到方差和协方差分量估值及其协因数阵如下:

      $$ \left\{ {\begin{array}{*{20}{l}} {\hat \sigma = {\mathit{\boldsymbol{N}}^{ - 1}}\mathit{\boldsymbol{l}}}\\ {{\mathit{\boldsymbol{Q}}_{\hat \sigma }} = {\mathit{\boldsymbol{N}}^{ - 1}}} \end{array}} \right. $$ (22)

      其中,Nl分别为计算方差分量的法方程的系数矩阵和常数向量,其元素分别为:

      $$ \left\{ \begin{array}{l} {n_{{j_1}{j_2}}} = \frac{1}{2}{\rm{tr}}\left( {\mathit{\boldsymbol{Q}}_{{j_1}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot \mathit{\boldsymbol{Q}}_{{j_2}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot } \right)\\ {l_{{j_1}}} = \frac{1}{2}\mathit{\boldsymbol{V}}_{{{\hat y}_c}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{Q}}_{{j_1}}^{\bf{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}{\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} - \frac{1}{2}{\rm{tr}}\left( {\mathit{\boldsymbol{Q}}_0^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot \mathit{\boldsymbol{Q}}_{{j_1}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat y}}}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot } \right) \end{array} \right. $$ (23)

      式中,$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}$;$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\;\;\;{\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}\hat x$;${j_1}, {j_2} = 1, 2 \ldots d$。由以上推导可知,STLS问题的最小二乘方差分量估计法是由内外两个迭代循环过程构成,其中内循环用于计算方差和协方差分量的估值,外循环用于计算参数和随机误差的估值。具体的迭代步骤为:

      1) 数据准备:AyQpψ,内循环收敛条件ε = 1 × 10-10,外循环收敛条件τ = 1 × 10-8

      2) 给定方差和协方差分量的初值$\mathit{\boldsymbol{\hat \sigma }} = {{\mathit{\boldsymbol{\hat \sigma }}}^0} = \left[ {\begin{array}{*{20}{c}} {\sigma _1^0}&{\sigma _2^0}& \ldots &{\sigma _t^0} \end{array}} \right]$。

      3) 给定参数初值$\hat x = {{\hat x}^0} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{y}}$。

      4) 设置外循环迭代次数s1 = 0,并进入外循环。

      5) 用式(16)更新${\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sum\limits_{k = 1}^d {\hat \sigma _k^{{s_1}}} {\mathit{\boldsymbol{Q}}_k}$。

      6) 设置内循环迭代次数s2 = 0,并进入内循环:

      (a) 用式(5)更新$\mathit{\boldsymbol{G}}\left( {{{\hat x}^{{s_2}}}} \right)$,并用式(9)更新${\mathit{\boldsymbol{Q}}_{\hat y_c^{{S_2}}}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{{s_2}}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{{s_2}}}} \right)$。

      (b) 用式(12)计算得${\rm{\Delta }}{{\hat p}^{{s_2}}}$。

      (c) 用式(13)计算得$\hat E_\mathit{\boldsymbol{A}}^{{s_2}}$,则有${{\hat A}^{{s_2}}} = A - \hat E_\mathit{\boldsymbol{A}}^{{s_2}}$。

      (d) 用式(8)计算得${{\hat x}^{{s_2} + 1}}$。

      (e) 重复内循环,直到$\;\left\| {{{\hat x}^{{s_2} + 1}} - {{\hat x}^{{s_2}}}} \right\| \le \varepsilon $。

      7) 用式(18)更新${\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}}$。

      8) 计算$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}$,${\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}\hat x$。

      9) 用式(23)计算得Nl

      10) 用式(22)计算得${{\hat \sigma }^{{s_1} + 1}} = {\mathit{\boldsymbol{N}}^{ - 1}}\mathit{\boldsymbol{l}}$,${\mathit{\boldsymbol{Q}}_{{{\hat \sigma }^{{S_1} + 1}}}} = {\mathit{\boldsymbol{N}}^{ - 1}}$。

      11) 重复外循环,直到$\left\| {{{\mathit{\boldsymbol{\hat \sigma }}}^{{s_1} + 1}} - {{\mathit{\boldsymbol{\hat \sigma }}}^{{s_1}}}} \right\| \le \tau $。

    • 现以某一建筑物某个时段的36期沉降观测[27]为例,高程观测值如表 1所示。文献[27]通过对参数的显著性进行检验,确定采用3阶AR模型对沉降数据进行分析和预测。为了便于比较,本文同样采用3阶AR模型,其对应的EIV模型如下:

      $$ \left\{ \begin{array}{l} H_4^{\rm{*}} = H_1^{\rm{*}}{\xi _1} + H_2^{\rm{*}}{\xi _2} + H_3^{\rm{*}}{\xi _3}\\ H_5^{\rm{*}} = H_2^{\rm{*}}{\xi _1} + H_3^{\rm{*}}{\xi _2} + H_4^{\rm{*}}{\xi _3}\\ H_6^{\rm{*}} = H_3^{\rm{*}}{\xi _1} + H_4^{\rm{*}}{\xi _2} + H_5^{\rm{*}}{\xi _3}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ H_{36}^{\rm{*}} = H_{33}^{\rm{*}}{\xi _1} + H_{34}^{\rm{*}}{\xi _2} + H_{35}^{\rm{*}}{\xi _3} \end{array} \right. $$ (24)

      表 1  沉降观测数据/mm

      Table 1.  Settlement Observation Data/mm

      序号 高程 序号 高程 序号 高程 序号 高程
      1 26.33 10 25.46 19 28.09 28 26.57
      2 26.27 11 26.12 20 26.78 29 28.36
      3 26.43 12 27.28 21 28.66 30 27.94
      4 25.56 13 26.67 22 26.75 31 26.81
      5 26.82 14 27.95 23 27.24 32 28.50
      6 26.56 15 26.74 24 28.02 33 27.68
      7 25.93 16 27.53 25 26.81 34 26.57
      8 26.43 17 25.31 26 28.50 35 28.36
      9 26.52 18 26.90 27 27.68 36 27.94

      式中,${\xi _1}$、${\xi _2}$、${\xi _3}$表示AR模型的模型参数;$H_i^{\rm{*}}\left( {i = 1, 2 \ldots 36} \right)$表示36期沉降观测的高程;*表示随机量。

      利用本文算法进行解算时,还需要给出ψ。首先由系数矩阵和观测向量构成的增广矩阵C可表示为

      $$ C = \left[ {\begin{array}{*{20}{c}} {H_1^{\rm{*}}}&{H_2^{\rm{*}}}&{H_3^{\rm{*}}}&{H_4^{\rm{*}}}\\ {H_2^{\rm{*}}}&{H_3^{\rm{*}}}&{H_4^{\rm{*}}}&{H_5^{\rm{*}}}\\ {H_3^{\rm{*}}}&{H_4^{\rm{*}}}&{H_5^{\rm{*}}}&{H_6^{\rm{*}}}\\ \vdots & \vdots & \vdots & \vdots \\ {H_{33}^{\rm{*}}}&{H_{34}^{\rm{*}}}&{H_{35}^{\rm{*}}}&{H_{36}^{\rm{*}}} \end{array}} \right] $$ (25)

      然后根据式(2)构造C的仿射结构矩阵S1~S30。不难发现,依次保留C各主对角线元素并置为1,即为S1~S30,进而可得$\mathit{\boldsymbol{\psi }} = \left[ {\begin{array}{*{20}{c}} {{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_1}} \right)}&{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_2}} \right)}& \ldots &{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_{30}}} \right)} \end{array}} \right]$。

      假设上述36组沉降数据来源于3个阶段的观测,其中第1阶段观测了12组数据,即${\mathit{\boldsymbol{H}}_f} = {\left[ {\begin{array}{*{20}{c}} {{H_1}}&{{H_2}}& \ldots &{{H_{12}}} \end{array}} \right]^{\rm{T}}}$;第2阶段也观测了12组数据,即${\mathit{\boldsymbol{H}}_s} = {\left[ {\begin{array}{*{20}{c}} {{H_{13}}}&{{H_{14}}}& \ldots &{{H_{24}}} \end{array}} \right]^{\rm{T}}}$;第3阶段同样观测了12组数据,即${\mathit{\boldsymbol{H}}_t} = {\left[ {\begin{array}{*{20}{c}} {{H_{25}}}&{{H_{26}}}& \ldots &{{H_{36}}} \end{array}} \right]^{\rm{T}}}$。参照式(16),假定函数模型式(24)对应的随机模型如下:

      $$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sigma _f^2{\mathit{\boldsymbol{Q}}_f} + \sigma _s^2{\mathit{\boldsymbol{Q}}_s} + \sigma _t^2{\mathit{\boldsymbol{Q}}_t} $$ (26)

      式中,$\sigma _f^2$、$\sigma _s^2$、$\sigma _t^2$分别为第1阶段、第2阶段、第3阶段沉降观测的方差,其对应的协因数矩阵分别为${\mathit{\boldsymbol{Q}}_f}$、${\mathit{\boldsymbol{Q}}_s}$、${\mathit{\boldsymbol{Q}}_t}$;${\mathit{\boldsymbol{Q}}_0}$表示随机模型的已知部分,其中,

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_0} = {0_{36}}, \;{\mathit{\boldsymbol{Q}}_f} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}} \end{array}} \right], }\\ {{\mathit{\boldsymbol{Q}}_s} = \left[ {\begin{array}{*{20}{c}} {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{\mathit{\boldsymbol{I}}_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}} \end{array}} \right], {\mathit{\boldsymbol{Q}}_t} = \left[ {\begin{array}{*{20}{c}} {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{\mathit{\boldsymbol{I}}_{12}}} \end{array}} \right]} \end{array} $$ (27)

      式中,I12表示12×12的单位矩阵;012表示12×12的零矩阵;036表示36×36的零矩阵。

      文献[17-19]的算法并不适用于AR模型的方差分量估计,在此对文献[17]中的算法进行改进,称为Amiri-Simkooei改进算法,即考虑系数矩阵和观测向量之间的相关性,建立132×132的观测权阵Q,同时进一步修改文献[17]中${\mathit{\boldsymbol{Q}}_{\tilde y}}$的表达式为${\mathit{\boldsymbol{Q}}_{\tilde y}} = {{\hat B}^{\rm{T}}}\mathit{\boldsymbol{Q}}\hat B$(其中,$\hat B = \left[ {{{\hat x}^T} \otimes {\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{I}}_m}} \right]$])[17],由于AR模型的系数矩阵和观测向量误差同源,观测权阵Q将变成非对角秩亏矩阵。

      假设方差分量的迭代初值分别为$\sigma _f^2\left( 0 \right) = 1$,$\sigma _s^2\left( 0 \right) = 1$,$\sigma _t^2\left( 0 \right) = 1$,分别采用最小二乘估计(LS估计)、PEIV模型解法[10]、Amiri-Simkooei改进算法[17]、本文算法进行求解,参数估计的迭代收敛阈值设为ε = 1 × 10-10,方差分量估计的迭代收敛阈值设为τ = 1 × 10-8

      表 2列出了各种算法得到的参数估值及其精度(标准差),表 3列出了Amiri-Simkooei改进算法和本文算法得到的方差分量估值及其精度,表 4列出了LS估计、PEIV模型解法得到的单位权方差估值及其精度,图 1给出了本文算法的方差分量收敛曲线。

      表 2  不同算法得到的参数估值及其精度

      Table 2.  Estimates of Model Parameters and Their Deviations from Different Methods

      项目 参数 LS估计 PEIV模型解法 Amiri-Simkooei改进算法 本文算法
      估值 ${{\hat \xi }_1}$ 0.635 1 1.179 1 1.075 5 1.075 5
      ${{\hat \xi }_2}$ 0.327 8 0.041 9 0.005 5 0.005 5
      ${{\hat \xi }_3}$ 0.041 1 -0.214 4 -0.075 2 -0.075 2
      标准差 ${{\hat \sigma }_{{\xi _1}}}$ 0.150 7 0.112 5 0.034 7 0.034 7
      ${{\hat \sigma }_{{\xi _2}}}$ 0.133 1 0.097 4 0.032 7 0.032 7
      ${{\hat \sigma }_{{\xi _3}}}$ 0.145 5 0.095 3 0.031 5 0.031 5

      表 3  不同算法得到的方差分量估值及其精度

      Table 3.  Estimates of Variance Factors and Their Deviations from Different Methods

      项目 参数 Amiri-Simkooei改进算法 本文算法
      估值 $\hat \sigma _f^2$ 0.178 1 0.178 1
      $\hat \sigma _s^2$ 1.015 6 1.015 6
      $\hat \sigma _t^2$ 0.081 3 0.081 3
      标准差 ${{\hat \sigma }_{\sigma _f^2}}$ 0.083 5 0.083 5
      ${{\hat \sigma }_{\sigma _s^2}}$ 0.423 2 0.423 2
      ${{\hat \sigma }_{\sigma _t^2}}$ 0.039 9 0.039 9

      表 4  不同算法得到的单位权方差估值及其精度

      Table 4.  Estimates of Unit Weight Variance and Their Deviations from Different Methods

      项目 参数 LS估计 PEIV模型解法
      估值 ${\hat \sigma _0^2}$ 0.804 7 0.643 4
      标准差 ${{\hat \sigma }_{\hat \sigma _0^2}}$ 0.167 2 0.106 9

      图  1  方差分量估值变化图

      Figure 1.  Changes of Estimated Variance Factors

      表 2可知,LS估计和PEIV模型解法所得的参数估值差异明显,在假定36期观测等精度的前提下,PEIV模型解法给出的模型估值更加可靠。在引入3个方差分量σf2σs2σt2后,Amiri-Simkooei改进算法和本文算法所得的参数估值产生了进一步的改变,这说明观测权阵对参数估计具有显著影响,合理的观测权阵是参数估计可靠性的前提。同时,本文算法和Amiri-Simkooei改进算法的计算结果一致,说明了本文算法的有效性以及WTLS估计[12]和本文给出的STLS估计具有等价性,但本文算法显著减少了待估参数个数,也无需重新构造秩亏的观测权阵,因而在一定程度上可以提高算法的计算效率。通过表 3表 4可知,本例引入的3个方差分量σf2σs2σt2估值存在显著差异,并且与PEIV模型解法得到的单位权方差估值也存在差异,因而可以得出结论:文献[27]对沉降数据进行处理时做出的等精度观测假设并不合适,由于对沉降观测没有进一步的了解,因此无法给出合理的随机模型对表 1中的沉降数据进行可靠的分析。

      需要注意的是,本文给出的AR模型方差分量估计方法最终归结为一个非凸优化问题,所以并不能保证得到其全局最优解,其全局收敛性与模型的非线性程度、观测方程的不吻合程度等相关,并且在算法设计上采用了高斯法,因而收敛性也与初值相关[28]

    • AR模型是一类重要的时间序列模型,它在测量数据处理中有广泛应用,AR模型的求解本质上是一个STLS问题。针对AR模型平差时观测权阵分配不合理、随机模型不准确的情况,本文首先基于变量投影法给出了一种STLS估计算法,然后将其与LS-VCE相结合,推导出了STLS问题的一种方差分量估计算法。该算法可以直接应用于AR模型的方差分量估值,并通过实测算例验证了它的有效性。结论如下:

      1)本文给出的STLS估计算法可以认为是文献[7]的一种替代算法,该算法在一定程度上简化了推导过程,提高了计算效率,并给出了参数的一阶近似精度。

      2)本文将STLS问题转化为附加参数的条件平差问题进行求解,即将TLS问题纳入到最小二乘平差理论体系中,因而能够直接和现有的方差分量估计算法相结合,如最小二乘方差分量估计。本文算法适用于多类数据的方差分量估计,并且也可以进行协方差分量的估计。

      3)本文只讨论了EIV模型系数矩阵状态良好、观测量只存在随机误差时的方差分量估计问题,EIV模型下方差分量的稳健估计算法和正则化估计算法是下一步的研究重点。

参考文献 (28)

目录

    /

    返回文章
    返回