留言板

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

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

总体最小二乘回归预测模型的方差分量估计

王乐洋 孙坚强

王乐洋, 孙坚强. 总体最小二乘回归预测模型的方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
引用本文: 王乐洋, 孙坚强. 总体最小二乘回归预测模型的方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
WANG Leyang, SUN Jianqiang. Variance Components Estimation for Total Least-Squares Regression Prediction Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
Citation: WANG Leyang, SUN Jianqiang. Variance Components Estimation for Total Least-Squares Regression Prediction Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450

总体最小二乘回归预测模型的方差分量估计

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

国家自然科学基金 41874001

国家自然科学基金 41664001

江西省杰出青年人才资助计划 20162BCB23050

国家重点研发计划 2016YFB0501405

详细信息
    作者简介:

    王乐洋,博士,教授,主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P207

Variance Components Estimation for Total Least-Squares Regression Prediction Model

Funds: 

The National Natural Science Foundation of China 41874001

The National Natural Science Foundation of China 41664001

the Support Program for Outstanding Youth Talents in Jiangxi Province 20162BCB23050

the National Key Research and Development Program of China 2016YFB0501405

More Information
    Author Bio:

    WANG Leyang, PhD, professor, specializes in the geodetic inversion and geodetic data processing. E-mail: wleyang@163.com

  • 摘要: 回归预测模型是对传统回归模型的进一步扩展,不仅涉及回归模型的固定参数估计,而且将模型预测纳入平差的部分内容,更加符合实际解算需求。针对在回归模型预测中经常出现待预测非公共点(自变量)含有观测误差和随机模型不准确的问题,基于EIV(errors-in-variables)模型提出了一种同时顾及所有变量观测误差的整体解法。同时,将方差-协方差分量估计方法引入回归预测模型解算中,以修正随机模型与待预测非公共点的先验协因数阵,并推导了相关计算公式和迭代算法。算例试验表明,该方法能够有效估计各类观测数据的方差分量,为获取更合理的参数估计与更高的模型预测精度提供了可行手段。另外,通过设计多种对比方案可知,该方法的预测效果较好,尤其是针对观测数据与系数矩阵中随机元素之间存在一定相关性的情况。
  • 图  1  方案5求得的方差分量收敛图

    Figure  1.  Convergence of Variance Components Derived from Scheme 5

    图  2  3次试验求得的方差分量收敛图

    Figure  2.  Convergences of Variance Components Derived from Three Tests

    图  3  本文方法求得的方差分量变化图

    Figure  3.  Convergence of the Variance Components Performed from the Presented Scheme

    表  1  模拟的观测数据

    Table  1.   Simulated Observation Data

    点号 观测数据
    yi xi xp(i)
    1 3.593 800 0.980 963 10.974 960
    2 5.411 958 1.984 376 11.923 399
    3 7.299 314 2.989 785 12.890 481
    4 9.025 751 3.997 224 13.897 391
    5 11.404 835 5.013 869 14.922 087
    6 13.593 908 6.029 296 15.944 041
    7 15.475 807 7.050 611 16.951 830
    8 17.226 161 8.065 007 17.943 819
    9 19.691 151 9.052 208 18.943 993
    10 21.647 393 10.017 662 19.960 737
    下载: 导出CSV

    表  2  5种方案解算的结果

    Table  2.   Results of Five Solution Schemes

    估值 方案1 方案2 方案3 方案4 方案5
    ${\beta _1}$ 2.005 3 2.003 2 2.002 6 2.002 6 2.002 6
    ${\beta _2}$ 1.371 8 1.435 2 1.452 2 1.452 2 1.452 2
    ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.128 3 0.064 9 0.047 8 0.047 8 0.047 8
    ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.587 2 0.491 9 0.465 5 0.377 5 0.377 5
    下载: 导出CSV

    表  3  μ1=10、μ2=10时5种方案解算的结果

    Table  3.   Results of Five Solution Schemes When μ1=10、μ2=10

    估值 方案1 方案2 方案3 方案4 方案5
    ${\beta _1}$ 2.005 3 2.005 0 2.002 6 2.002 6 2.002 6
    ${\beta _2}$ 1.371 8 1.381 2 1.452 2 1.452 2 1.452 2
    ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.128 3 0.118 9 0.047 8 0.047 8 0.047 8
    ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.587 2 0.572 3 0.465 5 0.377 5 0.377 5
    下载: 导出CSV

    表  4  模拟的观测数据、权值及其对应的相关系数

    Table  4.   Simulated Observations, Weights and the Corresponding Correlation Coefficients

    点号 观测数据 权值 相关系数
    yi xi xp(i) Pyi Pxi Pxp(i) ρxiyi ρxixp(i)
    1 11.288 5 4.995 0 10.005 1 1.0 1000.0 1000.0 -0.165 956 -0.314 723
    2 12.288 7 5.396 5 10.498 6 1.8 1000.0 1000.0 0.440 649 0.405 792
    3 13.150 3 5.803 3 10.991 8 4.0 500.0 500.0 -0.599 771 -0.373 013
    4 13.878 5 6.194 7 11.499 6 8.0 800.0 800.0 -0.395 335 -0.413 375
    5 14.748 5 6.579 9 11.985 0 20.0 200.0 200.0 -0.706 488 -0.132 359
    6 15.504 9 6.977 0 12.480 3 20.0 80.0 80.0 -0.715 323 0.402 459
    7 16.343 0 7.401 0 13.016 3 70.0 60.0 60.0 -0.627 480 -0.221 501
    8 17.102 2 7.837 1 13.509 1 70.0 20.0 20.0 -0.308 879 0.046 881
    9 17.884 9 8.142 8 13.784 5 100.0 1.8 1.8 -0.206 465 0.457 506
    10 18.699 1 8.521 8 14.318 9 500.0 1.0 1.0 0.077 633 0.464 888
    下载: 导出CSV

    表  5  4种方案解算的结果

    Table  5.   Results of Four Schemes

    估值 方案1 方案2 方案3 方案4 真值
    ${\beta _1}$ 1.997 9 1.997 9 2.001 0 2.001 3 2
    ${\beta _2}$ 1.545 8 1.545 8 1.526 2 1.523 9 1.5
    ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.045 9 0.045 9 0.026 2 0.023 9
    ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.545 4 0.484 0 0.462 9 0.461 0
    下载: 导出CSV
  • [1] Adcok R. A Problem in Least Squares[J]. Analyst, 1878, 5: 53-54 doi:  10.2307/2635758
    [2] Pearson K. On Lines and Planes of Closest Fit to Systems of Points in Space[J]. Philosophical Magazine and Journal of Science, 2010, 2(11): 559-572
    [3] York D. Least-Squares Fitting of a Straight Line[J]. Canadian Journal of Physics, 1996, 44(5): 1 079-1 086
    [4] Golub G H, van Loan C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6): 883-893 doi:  10.1137/0717073
    [5] Fang Xing, Wang Jin, Li Bofeng, et al. On Total Least Squares for Quadratic Form Estimation[J]. Studia Geophysica et Geodaetica, 2015, 59(3): 366-379 doi:  10.1007/s11200-014-0267-x
    [6] Shen Yunzhong, Li Bofeng, Chen Yi. An Iterative Solution of Weight Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238 doi:  10.1007/s00190-010-0431-1
    [7] Fang Xing. Weighted Total Least Squares Solutions for Applications in Geodesy[D]. Germany: Leibniz Universität Hannover, 2011
    [8] Fang Xing. Weight Total Least-Squares with Constraints: A Universal Formula for Geodetic Symmetrical Transformations[J]. Journal of Geodesy, 2015, 89(5): 459-469 doi:  10.1007/s00190-015-0790-8
    [9] Jazaeri S, Amiri-Simkooei A R, Sharfi M A. Iterative Algorithm for Weight Total Least Squares Adjustment[J]. Survey Review, 2014, 46(334): 19-27 doi:  10.1179/1752270613Y.0000000052
    [10] Amiri-Simkooei A R, Jazaeri S. Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]. Jouranal of Geodetic Science, 2012, 2(2): 113-124 doi:  10.2478/v10156-011-0036-5
    [11] 汪奇生, 杨德宏, 杨建文.基于总体最小二乘的线性回归迭代算法[J].大地测量与地球动力学, 2013, 33(6): 112-114, 120 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201306026.htm

    Wang Qisheng, Yang Dehong, Yang Jianwen. An Iterative Algorithm of Linear Regression Based on Total Least Squares[J]. Journal of Geodesy and Geodynamics, 2013, 33(6): 112-114, 120 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201306026.htm
    [12] 孔建, 姚宜斌, 吴寒.整体最小二乘的迭代解法[J].武汉大学学报∙信息科学版, 2010, 35(6): 711-714 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201006024.htm

    Kong Jian, Yao Yibin, Wu Han. Iterative Method for Total Least-Squares[J]. Geomatics and Information Science of Wuhan University, 2010, 35(6):711-714 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201006024.htm
    [13] Schaffrin B, Wieser A. On Weighted Total Least-Squares Adjustment for Linear Regression[J]. Journal of Geodesy, 2008, 82(7):415-421 doi:  10.1007/s00190-007-0190-9
    [14] 王乐洋, 赵英文, 陈晓勇, 等.多元总体最小二乘问题的牛顿解法[J].测绘学报, 2016, 45(4): 411-417 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201604006.htm

    Wang Leyang, Zhao Yingwen, Chen Xiaoyong, et al. A Newton Algorithm for Multivariate Total Least Squares Problems[J].Acta Geodaetica et Cartographica Sinica, 2016, 45(4): 411-417 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201604006.htm
    [15] 许光煜. Partial EIV模型的总体最小二乘方法及应用研究[D].南昌: 东华理工大学, 2016

    Xu Guangyu. The Total Least Squares Method and Its Application of Partial Errors-in-Variables Model[D]. Nanchang: East China University of Technology, 2016
    [16] 陶武勇, 鲁铁定, 许光煜, 等.稳健总体最小二乘Helmert方差分量估计[J].大地测量与地球动力学, 2017, 37(11): 1 193-1 197 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201711019.htm

    Tao Wuyong, Lu Tieding, Xu Guangyu, et al. Helmert Variance Component Estimation for Robust Total Least Squares[J]. Journal of Geodesy and Geodynamics, 2017, 37(11): 1 193-1 197 https://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201711019.htm
    [17] 王乐洋, 许光煜.加权总体最小二乘平差随机模型的验后估计[J].武汉大学学报∙信息科学版, 2016, 41(2): 255-261 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201602018.htm

    Wang Leyang, Xu Guangyu. Application of Posteriori Estimation of a Stochastic Model on the Weighted Total Least Squares Problem[J]. Geomatics and Information Science of Wuhan University, 2016, 41(2): 255-261 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201602018.htm
    [18] 杨元喜, 张晓东.基于严密Helmert方差分量估计的动态Kalman滤波[J].同济大学学报(自然科学版), 2009, 37(9): 1 241-1 245 https://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ200909021.htm

    Yang Yuanxi, Zhang Xiaodong. Variance Component Estimation of Helmert Type-Based Dynamic Kalman Filtering[J]. Journal of Tongji University (Natural Science), 2009, 37(9): 1 241-1 245 https://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ200909021.htm
    [19] 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
    [20] 乐科军.基于Helmert方差分量估计的半参数回归模型若干算法研究[D].长沙: 中南大学, 2009

    Le Kejun. A Study on the Algorithms of Semiparametric Regression Model Based on Helmert Variance Component Estimation[D]. Changsha: Central South University, 2009
    [21] 赵俊, 郭建锋.方差分量估计的通用公式[J].武汉大学学报·信息科学版, 2013, 38(5): 580-583 http://ch.whu.edu.cn/article/id/2639

    Zhao Jun, Guo Jianfeng. Auniversal Formula of Variance Component Estimation[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 580-583 http://ch.whu.edu.cn/article/id/2639
    [22] Xu Peiliang, 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
    [23] Amiri-Simkooei A R. Least-Squares Variance Component Estimation: Theory and GPS Applications[D]. Delft: Delft University of Technology, 2007
    [24] 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
    [25] Amiri-Simkooei A R. Non-negative Least-Squares Variance Component Estimation with Application to GPS Time Series[J]. Journal of Geodesy, 2016, 90(5): 451-466 doi:  10.1007/s00190-016-0886-9
    [26] 王乐洋, 温贵森.相关观测PEIV模型的最小二乘方差分量估计[J].测绘地理信息, 2018, 43(1): 8-14 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG201801002.htm

    Wang Leyang, Wen Guisen. Least Squares Variance Components Estimation of PEIV Model with Correlated Observations[J]. Journal of Geomatics, 2018, 43(1): 8-14 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG201801002.htm
    [27] 王乐洋, 温贵森. Partial EIV模型的非负最小二乘方差分量估计[J].测绘学报, 2017, 46(7): 857-865 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201707009.htm

    Wang Leyang, Wen Guisen. Non-negative Least Squares Variance Component Estimation of Partial EIV Model[J]. Acta Geodetica et Cartographica Sinica, 2017, 46(7): 857-865 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201707009.htm
    [28] 陶叶青, 蔡安宁, 杨娟.方差分量估计的抗差总体最小二乘算法[J].测绘工程, 2018, 27(3): 1-5 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201803001.htm

    Tao Yeqing, Cai Anning, Yang Juan. Solution for Robust Total Least Squares Based on Variance Component Estimation[J]. Engineering of Surveying and Mapping, 2018, 27(3): 1-5 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201803001.htm
    [29] 刘志平, 张书毕.方差-协方差分量估计的概括平差因子法[J].武汉大学学报∙信息科学版, 2013, 38(8): 925-929 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201308011.htm

    Liu Zhiping, Zhang Shubi. Variance-Covariance Component Estimation Method Based on Generalization Adjustment Factor[J]. Geomatics and Information Science of Wuhan University, 2013, 38(8): 925-929 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201308011.htm
    [30] 王乐洋, 余航.附有相对权比的加权总体最小二乘联合平差方法[J].武汉大学学报∙信息科学版, 2019, 44(8): 1 233-1 240 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201908018.htm

    Wang Leyang, Yu Hang. Weighted Total Least Squares Method for Joint Adjustment Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1 233-1 240 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201908018.htm
    [31] 王苗苗, 李博峰.无缝线性回归与预测模型[J].测绘学报, 2016, 45(12): 1 396-1 405 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201612003.htm

    Wang Miaomiao, Li Bofeng. Seamless Linear Regression and Prediction Model[J]. Acta Geodetica et Cartographica Sinica, 2016, 45(12): 1 396-1 405 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201612003.htm
    [32] 李博峰, 沈云中, 李微晓.无缝三维基准转换模型[J].中国科学:地球科学, 2012, 42(7): 1 047-1 054 https://www.cnki.com.cn/Article/CJFDTOTAL-JDXK201207012.htm

    Li Bofeng, Shen Yunzhong, Li Weixiao. The Seamless Model for Three-Dimensional Datum Transformation[J]. Science China: Earth Sciences, 2012, 42(7): 1 047-1 054 https://www.cnki.com.cn/Article/CJFDTOTAL-JDXK201207012.htm
    [33] 李博峰.无缝仿射基准转换模型的方差分量估计[J].测绘学报, 2016, 45(1): 30-35 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201601006.htm

    Li Bofeng. Variance Component Estimation in the Seamless Affine Transformation Model[J]. Acta Geodetica et Cartographica Sinica, 2016, 45(1): 30-35 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201601006.htm
    [34] 李博峰, 沈云中, 楼立志.基于等效残差的方差-协方差分量估计[J].测绘学报, 2010, 39(4): 349-354, 363 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201004006.htm

    Li Bofeng, Shen Yunzhong, Lou Lizhi. Variance-Covariance Component Estimation Based on the Equivalent Residuals[J]. Acta Geodetica et Cartographica Sinica, 2010, 39(4): 349-354, 363 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201004006.htm
  • [1] 吕志鹏, 隋立芬.  基于变量投影法的自回归模型方差分量估计 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
    [2] 李林林, 周拥军, 周瑜.  二维直角建筑物规则化的加权总体最小二乘平差方法 . 武汉大学学报 ● 信息科学版, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
    [3] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [4] 刘春阳, 王坚, 王彬, 刘超, 刘纪平.  基于中位参数法相关观测的抗差加权整体最小二乘算法 . 武汉大学学报 ● 信息科学版, 2019, 44(3): 378-384. doi: 10.13203/j.whugis20160516
    [5] 苏勇, 范东明, 蒲星钢, 游为, 肖东升, 于冰.  联合GOCE卫星数据和GRACE法方程确定SWJTU-GOGR01S全球重力场模型 . 武汉大学学报 ● 信息科学版, 2018, 43(3): 457-463. doi: 10.13203/j.whugis20150100
    [6] 温扬茂, 冯怡婷.  地震破裂模型约束的中国阿里地震三维形变场 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1369-1375. doi: 10.13203/j.whugis20160450
    [7] 王乐洋, 李海燕, 陈晓勇.  拟牛顿修正法解算不等式约束加权总体最小二乘问题 . 武汉大学学报 ● 信息科学版, 2018, 43(1): 127-132. doi: 10.13203/j.whugis20150333
    [8] 王乐洋, 余航.  火山Mogi模型反演的总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
    [9] 王乐洋, 许光煜, 陈晓勇.  附有相对权比的PEIV模型总体最小二乘平差 . 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
    [10] 王乐洋, 余航.  总体最小二乘联合平差 . 武汉大学学报 ● 信息科学版, 2016, 41(12): 1683-1689. doi: 10.13203/j.whugis20140670
    [11] 王乐洋, 吴飞, 吴良才.  GPS高程转换的总体最小二乘拟合推估模型 . 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
    [12] 汪奇生, 杨德宏, 杨腾飞.  EIV模型参数估计的新方法 . 武汉大学学报 ● 信息科学版, 2016, 41(3): 356-360. doi: 10.13203/j.whugis20140182
    [13] 姚宜斌, 黄书华, 张良, 胡羽丰, 李国平.  求解三维坐标转换参数的整体最小二乘新方法 ! . 武汉大学学报 ● 信息科学版, 2015, 40(7): 853-857. doi: 10.13203/j.whugis20130302
    [14] 王乐洋, 许才军.  总体最小二乘研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(7): 850-856.
    [15] 刘经南, 曾文宪, 徐培亮.  整体最小二乘估计的研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(5): 505-512.
    [16] 鲁铁定, 周世健.  总体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1351-1354.
    [17] 童小华, 刘大杰, 龚健雅.  GIS中宗地面积的平差模型与方案探讨 . 武汉大学学报 ● 信息科学版, 2003, 28(2): 182-187.
    [18] 余学祥, 吕伟才.  基于标准化残差的相关观测抗差估计模型 . 武汉大学学报 ● 信息科学版, 1999, 24(1): 75-78.
    [19] 黄声享, 王金岭.  多余观测分量与可靠性度量指标研究 . 武汉大学学报 ● 信息科学版, 1997, 22(2): 114-118.
    [20] 许才军.  卫星网与地面网联合平差中随机模型的估计 . 武汉大学学报 ● 信息科学版, 1989, 14(3): 35-45.
  • 加载中
图(3) / 表(5)
计量
  • 文章访问数:  34
  • HTML全文浏览量:  4
  • PDF下载量:  17
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-11-10
  • 刊出日期:  2021-02-05

总体最小二乘回归预测模型的方差分量估计

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

    国家自然科学基金 41874001

    国家自然科学基金 41664001

    江西省杰出青年人才资助计划 20162BCB23050

    国家重点研发计划 2016YFB0501405

    作者简介:

    王乐洋,博士,教授,主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P207

摘要: 回归预测模型是对传统回归模型的进一步扩展,不仅涉及回归模型的固定参数估计,而且将模型预测纳入平差的部分内容,更加符合实际解算需求。针对在回归模型预测中经常出现待预测非公共点(自变量)含有观测误差和随机模型不准确的问题,基于EIV(errors-in-variables)模型提出了一种同时顾及所有变量观测误差的整体解法。同时,将方差-协方差分量估计方法引入回归预测模型解算中,以修正随机模型与待预测非公共点的先验协因数阵,并推导了相关计算公式和迭代算法。算例试验表明,该方法能够有效估计各类观测数据的方差分量,为获取更合理的参数估计与更高的模型预测精度提供了可行手段。另外,通过设计多种对比方案可知,该方法的预测效果较好,尤其是针对观测数据与系数矩阵中随机元素之间存在一定相关性的情况。

English Abstract

王乐洋, 孙坚强. 总体最小二乘回归预测模型的方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
引用本文: 王乐洋, 孙坚强. 总体最小二乘回归预测模型的方差分量估计[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
WANG Leyang, SUN Jianqiang. Variance Components Estimation for Total Least-Squares Regression Prediction Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
Citation: WANG Leyang, SUN Jianqiang. Variance Components Estimation for Total Least-Squares Regression Prediction Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
  • 在测量学领域中,数据处理是必不可少的环节。虽然数据处理的问题种类很多,但最终都可归结为对模型中的参数进行估计。最小二乘估计理论(least-squares,LS)具有良好的统计性质且使用简单,被广泛应用于各个领域。但在实际参数估计中,由于系数矩阵也含有观测值或者模型线性化等原因,使得系数矩阵为非常数阵(含有随机误差),此时使用LS估计方法显然不合理。近年来,同时顾及系数矩阵误差和观测值误差的总体最小二乘(total least-squares,TLS)方法[1-4]得到了迅速发展。TLS理论依据的函数模型一般为EIV(errors-in-variables)模型[5-10]。文献[11-12]研究了基于EIV模型推导的独立等精度下的TLS解法。文献[6]和文献[13]推导了观测数据不相关情形下的加权TLS解法。文献[14]基于多元变量EIV模型推导了多元加权TLS算法,解决了系数矩阵中随机元素与观测矩阵存在相关性的问题。文献[15]不限定观测数据之间的相关性,推导了不等精度情况下的3种加权TLS方法。

    随机模型的不准确会使平差结果产生较大影响,方差分量估计法通过二次定权来改善平差解算结果,从而达到修正随机模型的目的,该方法也称验后估计法。方差分量估计法主要包括赫尔默特法(Helmert)[16-21]、最小范数二次无偏估计(minimum norm quadratic unbiased estimation,MINQUE)[22]、最优不变二次无偏估计(best invariant quadratic unbiased estimation,BIQUE)、最小二乘方差分量估计(least-squares variance component estimation,LSVCE)[23-28]和方差-协方差分量估计[29-30]等。文献[16]基于EIV模型推导了稳健TLS赫尔默特方差分量估计方法。文献[22]以EIV模型为基础,推导了总体最小二乘MINQUE的相关公式,并从参数的可估性、模型的稳定性和偏差分析等方面进行了讨论。文献[26]顾及了观测量之间的相关信息,基于PEIV(partial EIV)模型推导了相关观测PEIV模型最小二乘方差估计器的迭代算法,并通过公式推导证明了所提出的算法与现有方法的等价性。文献[27]将非负约束法与最小二乘问题相结合,给出了PEIV模型非负LSVCE方法,有效避免了估计过程中可能出现的负方差现象。文献[28]中,为了缓解粗差造成参数估计结果严重有偏,等价权函数被用于TLS抗差估计算法的实现当中,进而利用LSVCE方法实现了先验随机模型的修正,进一步提升了抗差效果。文献[29]借助平差因子扩宽了方差分量估计的理论研究范围,提出了一种基于概括平差因子的方差分量估计新方法,解决了方差-协方差分量估计理论存在的问题。文献[30]在判别函数最小化的理论基础上,通过确定多源大地测量数据在联合平差问题中的权重占比,成功推导了附有相对权比的WTLS(weighted TLS)联合平差方法。以上研究在进行模型预测时均没有考虑待预测自变量(例如在基准变换中为非公共点)的观测误差。文献[31]基于EIV模型构建了TLS的无缝线性回归与预测模型,将待预测非公共点的观测误差纳入平差模型中联合求解。文献[32]拓展了新颖的无缝三维基准转换模型,除了求解常规的转换参数外,还充分利用该模型的特点给出了非公共点坐标观测误差的解决方案。文献[33]以仿射基准变换为例,考虑所有点位误差,对其进行联合处理,并推导了方差分量估计公式。

    通过以上分析可知,已有方法未考虑到随机模型不准确的情况,也未考虑到观测数据之间的相关性。为此,本文不限定观测数据与系数矩阵中随机元素之间的相关性,提出了一种能同时顾及所有变量误差的整体解法,并结合方差-协方差分量估计方法修正随机模型,提高平差解算结果,从而提高预测效果。通过设置不同类型的协因数阵进行算例试验,验证了本文方法的可行性及有效性。

    • 基于EIV模型的TLS回归预测模型可表示如下:

      $$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} = \left( {\mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{\beta }}}\\ {{\mathit{\boldsymbol{y}}_\mathit{\boldsymbol{p}}} = \left( {{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}} + {\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}}}} \right)\mathit{\boldsymbol{\beta }}} \end{array}} \right. $$ (1)

      式中,$\mathit{\boldsymbol{y}} \in {{\bf{R}}^{n \times 1}}$为模型中的观测向量;${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} \in {{\bf{R}}^{n \times 1}}$为对应于观测随机变量的误差项;$\mathit{\boldsymbol{A}} \in {{\bf{R}}^{n \times m}}$为受随机误差污染的系数矩阵;${\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}} \in {{\bf{R}}^{n \times m}}$为与系数矩阵A相对应的误差矩阵;${\mathit{\boldsymbol{y}}_\mathit{\boldsymbol{p}}} \in {{\bf{R}}^{n \times 1}}$为模型预测时的待估观测向量;${\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}$为模型预测中的系数矩阵;${\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}}} \in {{\bf{R}}^{k \times m}}$为对应于矩阵${\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}$的误差矩阵;$\mathit{\boldsymbol{\beta }} \in {{\bf{R}}^{m \times 1}}$为待估的固定参数。

      随机模型如下:

      $$ \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}} \end{array}} \right] \backsim N\left( {\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}} \end{array}} \right], \sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right]} \right) $$ (2)

      式中,${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}} \in {{\bf{R}}^{t \times 1}}$为通过抽取矩阵A中随机元素并以此排列组成的误差项;${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}} \in {{\bf{R}}^{r \times 1}}$为对应于矩阵${\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}$中随机元素所形成的误差项;${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}$为实际观测向量y的协因数阵;${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}$和${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}$为实际观测向量y和系数矩阵A中随机元素的协因数阵;${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}、{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}$为系数矩阵中随机元素的协因数阵;${\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}、{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}$为系数矩阵中随机元素的互协因数阵;$\sigma _0^2$为单位权方差。

    • 本文使用Gauss-Newton法求解非线性模型(1)的参数。设第j次迭代后参数的估值为${\left( j \right)}$,则$\mathit{\boldsymbol{\beta }}$可表示为:

      $$ \mathit{\boldsymbol{\beta }} = {\mathit{\boldsymbol{\beta }}_{\left( j \right)}} + {\rm{\Delta }}\mathit{\boldsymbol{\beta }} $$ (3)

      结合式(1)和式(3),使用${{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}$的第j次迭代值${{\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_{\left( j \right)}}}}}$取代${{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}$,以克服计算过程中出现的迭代不收敛或不正确等现象,即${{\mathit{\boldsymbol{A}}_{\left( j \right)}} = \mathit{\boldsymbol{A}} + {\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_{\left( j \right)}}}}}$,由此可得:

      $$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}} + {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}} + {\mathit{\boldsymbol{A}}_{\left( j \right)}}{\rm{\Delta }}\mathit{\boldsymbol{\beta }} = }\\ {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}} + \left( {\mathit{\boldsymbol{\beta }}_{\left( j \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right){\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{A}}} + {\mathit{\boldsymbol{A}}_{\left( j \right)}}{\rm{\Delta }}\mathit{\boldsymbol{\beta }}} \end{array} $$ (4)

      令${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{A}}} = {\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}, {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}}} = {\mathit{\boldsymbol{B}}_2}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}, \mathit{\boldsymbol{V}} = {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}$,其中${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{A}}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right), {\mathit{\boldsymbol{e}}_{{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}}} = {\rm{vec}}\left( {{\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_\mathit{\boldsymbol{p}}}}}} \right)$,${\rm{vec}}\left( \cdot \right)$为拉直运算符,${{\mathit{\boldsymbol{B}}_1}}、{{\mathit{\boldsymbol{B}}_2}}$是由常数构成的固定矩阵,为$mn \times t$维矩阵,则$\mathit{\boldsymbol{V}}$可表示如下:

      $$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}} + \left( {\mathit{\boldsymbol{\beta }}_{\left( j \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right){\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}} + {\mathit{\boldsymbol{A}}_{\left( j \right)}}{\rm{\Delta }}\mathit{\boldsymbol{\beta }} - \mathit{\boldsymbol{y}} $$ (5)

      式中,$ \otimes $为克罗内克积。式(5)可视为最小二乘配置的规范表达。将每次迭代更新的参数增量${{\rm{\Delta }}\mathit{\boldsymbol{\beta }}}$视为待估的固定参数(也称倾向参数),${{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}$和${{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}$视为随机信号参数,将期望值为${\bf{0}}$的${{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}$和${{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}$作为虚拟观测值,则得到下式:

      $$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{V}} = {\mathit{\boldsymbol{A}}_{\left( j \right)}}{\rm{\Delta }}\mathit{\boldsymbol{\beta }} + \left( {\mathit{\boldsymbol{\beta }}_{\left( j \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right){\mathit{\boldsymbol{B}}_1}{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}} - \left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}}} \right)}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}} = {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}} - {{\bf{0}}_{t \times 1}}}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}} = {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}} - {{\bf{0}}_{r \times 1}}} \end{array}} \right. $$ (6)

      式中,${y - A{\beta _{\left( j \right)}}}$可视为观测值与近似值之间的差值;$\mathit{\boldsymbol{V}}$为观测向量的残差向量。令$\mathit{\boldsymbol{l}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}}$,则${\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{l}}} = {\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}$。由此,TLS问题可转换为如下的最小二乘问题:

      $$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{V}}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}}}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{\left( j \right)}}}&{\left( {\mathit{\boldsymbol{\beta }}_{\left( j \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right){\mathit{\boldsymbol{B}}_1}{\rm{}}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{I}}_t}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{I}}_r}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\rm{\Delta }}\mathit{\boldsymbol{\beta }}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {y - A{\beta _{\left( j \right)}}}\\ \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}} \end{array}} \right] $$ (7)

      令$\mathit{\boldsymbol{\bar V}} = \left[{\begin{array}{*{20}{c}} \mathit{\boldsymbol{V}}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}}}\\ {{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}}} \end{array}} \right]$,${\mathit{\boldsymbol{G}}_{\left( j \right)}} = \left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{\left( j \right)}}}&{\left( {\beta _{\left( j \right)}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_n}} \right){\mathit{\boldsymbol{B}}_1}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{I}}_t}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{I}}_r}} \end{array}} \right]$,${\mathit{\boldsymbol{L}}_{\left( j \right)}} = \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\beta }}_{\left( j \right)}}}\\ \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}} \end{array}} \right]$,${{\mathit{\boldsymbol{\hat X}}}_j} = \left[{\begin{array}{*{20}{c}} {{\rm{\Delta }}\mathit{\boldsymbol{\beta }}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}} \end{array}} \right]$,求解式与间接平差类似,即:

      $$ {{\mathit{\boldsymbol{\hat X}}}_j} = {\left( {\mathit{\boldsymbol{G}}_{\left( j \right)}^{\rm{T}}{\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{G}}_{\left( j \right)}}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_{\left( j \right)}^{\rm{T}}{\mathit{\boldsymbol{G}}^{ - 1}}{\mathit{\boldsymbol{L}}_{\left( j \right)}} $$ (8)

      当待定参数满足收敛条件时,迭代终止,则与系数矩阵AAp对应的误差向量为eaep,将ep重构成为EAp,由式(5)可以获得最终模型预测的观测向量值。

      j+1次迭代计算结果满足收敛阈值条件,则参数的后验协方差阵为:

      $$ {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat X\hat X}}}} = \sigma _0^2{\left( {\mathit{\boldsymbol{G}}_{\left( j \right)}^{\rm{T}}{\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{G}}_{\left( j \right)}}} \right)^{ - 1}} = \sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat \beta \hat \beta }}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat \beta }}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{\hat \beta }}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}}}}\\ {{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}\mathit{\boldsymbol{\hat \beta }}}}}&{{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}}}}&{{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}}}}\\ {{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}\mathit{\boldsymbol{\hat \beta }}}}}&{{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{a}}}}}}&{{\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}{{\mathit{\boldsymbol{\hat e}}}_\mathit{\boldsymbol{p}}}}}} \end{array}} \right] $$ (9)
    • 先验随机模型的准确与否影响最终的平差结果,但由于实际数据采集中无法获取较为精确的随机模型,因此为了提高预测效果,通常使用方差分量估计方法对随机模型进行修正。§1.1中的TLS回归预测模型及求解算法主要涉及了最小二乘配置思想,将所有随机误差项进行归类,从而实现最小二乘估计问题的转化。借助方差-协方差分量估计思想[34],通过不断修正随机模型,以获取更为合理的估计结果。平差中观测变量的协方差矩阵为:

      $$ \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{yy}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\sigma _1^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&{\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right] $$ (10)

      使用如下表示:

      $$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{Q}}_1} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right], {\mathit{\boldsymbol{Q}}_2} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{Q_{aa}}}&{{Q_{ap}}}\\ \mathit{\boldsymbol{0}}&{{Q_{pa}}}&{{Q_{pp}}} \end{array}} \right], {\mathit{\boldsymbol{Q}}_3} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&{{Q_{ya}}}&\mathit{\boldsymbol{0}}\\ {{Q_{ay}}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right]\\ {\mathit{\boldsymbol{P}}_1} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{yy}}}}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right], {\mathit{\boldsymbol{P}}_2} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right], {\mathit{\boldsymbol{P}}_3} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {{\mathit{\boldsymbol{P}}_{\mathit{\boldsymbol{ay}}}}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}} \end{array}} \right] \end{array} \right. $$ (11)

      式中,Q1Q2Q3分别为观测值之间、系数矩阵随机元素与模型预测系数矩阵随机元素以及观测值与系数矩阵随机元素之间的协因数阵或互协因数阵;P1P2P3分别为对应的权阵。令$\mathit{\boldsymbol{\theta }} = \left[{\hat \sigma _1^2\ \ \hat \sigma _2^2\ \ \hat \sigma _3^2} \right]$,$\mathit{\boldsymbol{R}} = {\left( {\mathit{\boldsymbol{G}}_{\left( j \right)}^{\rm{T}}{\mathit{\boldsymbol{Q}}^{ - 1}}{\mathit{\boldsymbol{G}}_{\left( j \right)}}} \right)^{ - 1}}\mathit{\boldsymbol{G}}_{\left( j \right)}^{\rm{T}}{\mathit{\boldsymbol{G}}^{ - 1}}{\mathit{\boldsymbol{L}}_{\left( j \right)}} - \mathit{\boldsymbol{E}}$。采用方差-协方差分量估计理论可得:

      $$ \sum\limits_{i = 1}^3 E \left( {{{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{\bar V}}} \right) = \sum\limits_{i = 1}^3 {\sum\limits_{j = 1}^3 {{\rm{tr}}} } \left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_i}{\mathit{\boldsymbol{\theta }}_j}} \right) $$ (12)

      即有:

      $$ \left[ {\begin{array}{*{20}{c}} {{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_1}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_2}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_3}} \right)}\\ {{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_1}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_2}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_3}} \right)}\\ {{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_3}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_1}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_3}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_2}} \right)}&{{\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_3}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_3}} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\hat \sigma _1^2}\\ {\hat \sigma _2^2}\\ {\hat \sigma _3^2} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{RL}}}\\ {{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{RL}}}\\ {{\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_3}\mathit{\boldsymbol{RL}}} \end{array}} \right] $$ (13)

      式中,${\hat \sigma _1^2}$、${\hat \sigma _2^2}$、${\hat \sigma _3^2}$分别表示${{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}}$、${{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{a}}}}$以及互协因数阵所对应的方差分量估值;${{\rm{tr}}\left( \cdot \right)}$表示矩阵求迹运算。

      由式(13)可知,各类观测向量所对应的方差分量估值为:

      $$ \mathit{\boldsymbol{\hat \theta }} = {\mathit{\boldsymbol{S}}^{ - 1}}\mathit{\boldsymbol{W}} $$ (14)

      式中,${\mathit{\boldsymbol{S}}_{ij}} = {\rm{tr}}\left( {{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{Q}}_j}} \right)$;${\mathit{\boldsymbol{W}}_i} = {\mathit{\boldsymbol{L}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_i}\mathit{\boldsymbol{RL}}$。

      TLS回归预测模型的方差-协方差分量估计步骤如下:

      1)数据采集与准备,即yAQ

      2)通过LS准则求解参数初值${\mathit{\boldsymbol{\beta }}_0}, {\mathit{\boldsymbol{E}}_{{\mathit{\boldsymbol{A}}_0}}} = \mathit{\boldsymbol{0}}$;

      3)经过首次TLS平差,得到$\mathit{\boldsymbol{V}}、{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}}、{\mathit{\boldsymbol{V}}_{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}}}$,迭代阈值取1×10-10

      4)按式(14)计算相应的方差分量因子$\hat \sigma _1^2、\hat \sigma _2^2、\hat \sigma _3^2$,更新D,从而更新QiPi

      $$ \left[ {\begin{array}{*{20}{c}} {\sigma _1^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&{\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {\sigma _3^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}}&{\sigma _2^2{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right] $$ (15)

      5)重复步骤3)、4),当满足收敛条件$\sigma _1^2 = \sigma _2^2 = \sigma _3^2$且收敛阈值为1×10-10时,终止循环。

      本文充分借助了最小二乘配置思想,实现了复杂TLS问题向LS问题的转化。在TLS回归预测模型的场景下,同时考虑了模型固定参数与同模型预测相关的非公共点随机误差的处理问题,由于可能存在不准确的先验随机模型,进而拓展了TLS回归预测模型方差-协方差分量估计的迭代算法。与文献[28]相比,由于顾及了模型预测时非公共点可能存在的误差,因此其理论更为严密,预测效果更佳。相较于文献[33],本文进一步考虑了系数矩阵中随机元素与观测向量相关的情形;相比于文献[31]和文献[32],本文为TLS回归预测模型的求解方法提供了一种新的思路,其结合了方差-协方差分量估计思想,在一定程度上提升了模型的预测精度。另外,由于此方差分量估计方法简化了公式推导过程,在基准变换和高程拟合等方面应用较为方便。

    • 为了体现本文方法的优势,采用两种不同类型的协因数阵进行试验。

      类型1:

      $$ \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}} \end{array}} \right] \backsim N\left( {\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}} \end{array}} \right], \sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ya}}}}}&\mathit{\boldsymbol{0}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ay}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right]} \right) $$

      类型2:

      $$ \left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{a}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{p}}}} \end{array}} \right] \backsim N\left( {\left[ {\begin{array}{*{20}{l}} \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}} \end{array}} \right], \sigma _0^2\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yy}}}}}&\mathit{\boldsymbol{0}}&\mathit{\boldsymbol{0}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{ap}}}}}\\ \mathit{\boldsymbol{0}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{pp}}}}} \end{array}} \right]} \right) $$
    • 模拟一条直线yi = β1xi + β2β1β2的真值分别为2、1.5。xi的取值范围为1~10,间隔为1;xp(i)的取值范围为11~20,间隔为1,根据上述直线方程式分别计算相应的yi值。本文暂不涉及上述两类数据之间的相关性。直线拟合回归预测模型为:

      $$ \left\{ \begin{array}{l} \left[ {\begin{array}{*{20}{l}} {{y_1}}\\ {{y_2}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots }\\ {{y_n}} \end{array}} \right] - \left[ {\begin{array}{*{20}{l}} {{e_{{y_1}}}}\\ {{e_{{y_2}}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots }\\ {{e_{{y_n}}}} \end{array}} \right] = \left( {\left[ {\begin{array}{*{20}{c}} {{x_1}}&1\\ {{x_2}}&1\\ \vdots & \vdots \\ {{x_n}}&1 \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_{{x_1}}}}&0\\ {{e_{{x_2}}}}&0\\ \vdots & \vdots \\ {{e_{{x_n}}}}&0 \end{array}} \right]} \right)\left[ {\begin{array}{*{20}{l}} {{\beta _1}}\\ {{\beta _2}} \end{array}} \right]\\ \left[ {\begin{array}{*{20}{l}} {{y_{p\left( 1 \right)}}}\\ {{y_{p\left( 2 \right)}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \vdots }\\ {{y_{p\left( n \right)}}} \end{array}} \right] = \left( {\left[ {\begin{array}{*{20}{c}} {{x_{p\left( 1 \right)}}}&1\\ {{x_{p\left( 2 \right)}}}&1\\ \vdots & \vdots \\ {{x_{p\left( n \right)}}}&1 \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{e_{p\left( 1 \right)}}}&0\\ {{e_{p\left( 2 \right)}}}&0\\ \vdots & \vdots \\ {{e_{p\left( n \right)}}}&0 \end{array}} \right]} \right)\left[ {\begin{array}{*{20}{l}} {{\beta _1}}\\ {{\beta _2}} \end{array}} \right] \end{array} \right. $$ (16)

      式中,$\mathit{\boldsymbol{\beta }} = {\left[{\begin{array}{*{20}{c}} {{\beta _1}}&{{\beta _2}} \end{array}} \right]^{\rm{T}}}$表示未知的待求参数。模拟的因变量的方差阵为:

      $$ {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{yy}}}} = {\left( {0.1 \times 2} \right)^2}{\mathit{\boldsymbol{I}}_n} $$ (17)

      考虑到自变量之间存在相关性,本文由两点之间的距离决定自变量的相关系数,主要依据下式:

      $$ {\rho _{ij}} = \frac{1}{{1 + {{\left( {{S_{ij}}/{S_0}} \right)}^2}}} $$ (18)

      式中,${{S_{ij}}}$代表两个坐标点之间的欧氏距离;${{S_0}}$取值为10。

      根据式(18)可确定自变量的相关系数,从而得到相关的方差阵,即:

      $$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{aa}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ap}}}}}\\ {{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{pa}}}}}&{{\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{ap}}}}} \end{array}} \right] = {\left( {0.1 \times 1} \right)^2}\left[ {\begin{array}{*{20}{c}} {{\rho _{1, 1}}}&{{\rho _{1, 2}}}& \cdots &{{\rho _{1, 2n}}}\\ {{\rho _{2, 1}}}&{{\rho _{2, 2}}}& \cdots &{{\rho _{2, 2n}}}\\ \vdots & \vdots &{}& \vdots \\ {{\rho _{2n, 1}}}&{{\rho _{2n, 2}}}& \cdots &{{\rho _{2n, 2n}}} \end{array}} \right] $$ (19)

      进而可以得到协因数阵Q

      依据本文方法,给定固定常数阵B1为:

      $$ {\mathit{\boldsymbol{B}}_1} = \left[ {{\mathit{\boldsymbol{I}}_{10}};\mathop 0\limits_{10 \times 10} } \right] $$ (20)

      通过MATLAB中的mvnrnd函数设置因变量和自变量真值,分别添加期望为0、方差阵如式(17)和式(19)给定的随机误差,由此得到一套模拟的观测数据(见表 1)。

      表 1  模拟的观测数据

      Table 1.  Simulated Observation Data

      点号 观测数据
      yi xi xp(i)
      1 3.593 800 0.980 963 10.974 960
      2 5.411 958 1.984 376 11.923 399
      3 7.299 314 2.989 785 12.890 481
      4 9.025 751 3.997 224 13.897 391
      5 11.404 835 5.013 869 14.922 087
      6 13.593 908 6.029 296 15.944 041
      7 15.475 807 7.050 611 16.951 830
      8 17.226 161 8.065 007 17.943 819
      9 19.691 151 9.052 208 18.943 993
      10 21.647 393 10.017 662 19.960 737

      为了说明本文方法的优势,采用如下5种方案进行对比分析:

      方案1:常规最小二乘方法(LS)。

      方案2:未顾及相关性的总体最小二乘方法(TLS)。

      方案3:未顾及非公共点随机误差的总体最小二乘Helmert方差分量估计(TLS-VCE)。

      方案4:采用文献[33]中的算法(L-TLS-VCE)。

      方案5:采用本文方法(P-TLS-VCE)。

      根据模拟的数据,分别采用上述5种方案进行算法实现,由此获得参数估计值、参数估计值${\mathit{\boldsymbol{\hat \beta }}}$与真值(${\mathit{\boldsymbol{\hat \beta }}_{{\rm{true }}}}$)之间的差值二范数、被预测观测向量${\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{ }}}}$与真值(${\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}\_{\rm{true }}}}$)之间的差值二范数(见表 2)。方案5的方差因子迭代变化如图 1所示。

      表 2  5种方案解算的结果

      Table 2.  Results of Five Solution Schemes

      估值 方案1 方案2 方案3 方案4 方案5
      ${\beta _1}$ 2.005 3 2.003 2 2.002 6 2.002 6 2.002 6
      ${\beta _2}$ 1.371 8 1.435 2 1.452 2 1.452 2 1.452 2
      ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.128 3 0.064 9 0.047 8 0.047 8 0.047 8
      ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.587 2 0.491 9 0.465 5 0.377 5 0.377 5

      图  1  方案5求得的方差分量收敛图

      Figure 1.  Convergence of Variance Components Derived from Scheme 5

      对比表 2中的二范数结果可以发现,由于方案2是同时顾及系数矩阵误差和观测值误差的平差方法,因此求得的参数估值结果要优于只考虑观测值误差的方案1的结果。同时,对比5种方案的结果可以得出,方案3~5通过验后估计各类观测量的验前方差,以此修正协因数阵,因此求得的参数估值结果比方案1和方案2更优。此外,对比3种方差分量估计方法的结果(方案3、4、5)可以看出,3种方法都可以得到较好的平差结果,但方案4和方案5不仅实现了先验随机模型的修正,而且对待预测非公共点坐标的协因数阵亦进行了修正,因此获得的预测观测向量的二范数相对较小。方案4和方案5求得的参数估值和预测因变量的二范数结果均一致,这是由于当Qay=0时,本文方法退化成L-TLS-VCE法,因此两种方法求得的二范数结果均一致。

      考虑到上述试验存在随机性,为进一步分析本文方法在不同先验方差下均能取得较好的结果,并证明同文献[33]存在等价的效果,分别设计9次试验,且任一试验均分配因变量和自变量不同的先验方差,即:

      $$ \left\{ {\begin{array}{*{20}{l}} {\sigma _1^2 = {\mu _1} \times 0.01}\\ {\sigma _2^2 = {\mu _2} \times 0.04} \end{array}} \right. $$ (21)

      式中,μ1的取值范围为1~10;μ2为与μ1相对应的逆序取值,即10~1。

      分别取μ1=1, 5.5, 10,μ2=10, 5.5, 1,3次试验求得的方差分量迭代变化情况见图 2表 3展示了当μ1=10、μ2=10时5种方案解算的结果。

      图  2  3次试验求得的方差分量收敛图

      Figure 2.  Convergences of Variance Components Derived from Three Tests

      表 3  当μ1=10、μ2=10时5种方案解算的结果

      Table 3.  Results of Five Solution Schemes When μ1=10、μ2=10

      估值 方案1 方案2 方案3 方案4 方案5
      ${\beta _1}$ 2.005 3 2.005 0 2.002 6 2.002 6 2.002 6
      ${\beta _2}$ 1.371 8 1.381 2 1.452 2 1.452 2 1.452 2
      ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.128 3 0.118 9 0.047 8 0.047 8 0.047 8
      ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.587 2 0.572 3 0.465 5 0.377 5 0.377 5

      根据图 2表 3中的计算结果可得,即便μ1μ2的取值发生改变,本文方法均能取得较为合理的预测效果,且求得的预测值的二范数几乎不随μ1μ2的取值改变而变化;同时文献[33]的测试结果也佐证了所提方法的有效性和合理性。

    • 假设真值x的范围为5~8.6,间隔为0.4;真值y的范围为11.5~18.7,间隔为0.8;真值xp的范围为10~14.5,间隔为0.5。同前文一样,假设上述两组数据是通过不同手段获取得到的,且考虑两者之间的相关性。直线拟合模型如式(16)所示,β1β2的真值分别为2和1.5。两组数据的相关系数由式(18)生成,并通过公式${\sigma _{{x_i}{y_i}}} = {\rho _{{x_i}{y_i}}}{\sigma _{{x_i}}}{\sigma _{{y_i}}}$和${\sigma _{{x_i}{x_{p\left( i \right)}}}} = {\rho _{{x_i}{x_{p\left( i \right)}}}}{\sigma _{{x_i}}}{\sigma _{{x_{p\left( i \right)}}}}$计算此时的协因数阵。

      令$\sigma _1^2 = 0.04$,$\sigma _2^2 = 0.02$,$\sigma _3^2 = 0.08$,由此可以获得与之对应的模拟的观测数据、权值及其相关系数(见表 4)。设计4种方案如下:

      表 4  模拟的观测数据、权值及其对应的相关系数

      Table 4.  Simulated Observations, Weights and the Corresponding Correlation Coefficients

      点号 观测数据 权值 相关系数
      yi xi xp(i) Pyi Pxi Pxp(i) ρxiyi ρxixp(i)
      1 11.288 5 4.995 0 10.005 1 1.0 1000.0 1000.0 -0.165 956 -0.314 723
      2 12.288 7 5.396 5 10.498 6 1.8 1000.0 1000.0 0.440 649 0.405 792
      3 13.150 3 5.803 3 10.991 8 4.0 500.0 500.0 -0.599 771 -0.373 013
      4 13.878 5 6.194 7 11.499 6 8.0 800.0 800.0 -0.395 335 -0.413 375
      5 14.748 5 6.579 9 11.985 0 20.0 200.0 200.0 -0.706 488 -0.132 359
      6 15.504 9 6.977 0 12.480 3 20.0 80.0 80.0 -0.715 323 0.402 459
      7 16.343 0 7.401 0 13.016 3 70.0 60.0 60.0 -0.627 480 -0.221 501
      8 17.102 2 7.837 1 13.509 1 70.0 20.0 20.0 -0.308 879 0.046 881
      9 17.884 9 8.142 8 13.784 5 100.0 1.8 1.8 -0.206 465 0.457 506
      10 18.699 1 8.521 8 14.318 9 500.0 1.0 1.0 0.077 633 0.464 888

      方案1:附有相关观测下的TLS方法(PEIV-TLS)。

      方案2:采用文献[31]中的算法(G-TLS)。

      方案3:未顾及相关性的本文方法(P-TLS-VCE)。

      方案4:顾及相关性的本文方法。

      根据上述观测数据,分别采用以上4种设计方案进行算法实现,其计算结果见表 5。本文方法获得的方差分量迭代变化如图 3所示。结合图 3表 5可以看出,本文提出的回归预测模型的方差分量估计方法由于考虑了系数矩阵与观测值相关,因此求得的参数估计值与被预测观测向量之间的差值二范数结果相对较小。

      表 5  4种方案解算的结果

      Table 5.  Results of Four Schemes

      估值 方案1 方案2 方案3 方案4 真值
      ${\beta _1}$ 1.997 9 1.997 9 2.001 0 2.001 3 2
      ${\beta _2}$ 1.545 8 1.545 8 1.526 2 1.523 9 1.5
      ${\left\| {\mathit{\boldsymbol{\hat \beta }} - {\mathit{\boldsymbol{\beta }}_{{\rm{true }}}}} \right\|}$ 0.045 9 0.045 9 0.026 2 0.023 9
      ${\left\| {{{\mathit{\boldsymbol{\hat y}}}_\mathit{\boldsymbol{p}}} - {\mathit{\boldsymbol{y}}_{\mathit{\boldsymbol{p}}{\rm{\_true }}}}} \right\|}$ 0.545 4 0.484 0 0.462 9 0.461 0

      图  3  本文方法求得的方差分量变化图

      Figure 3.  Convergence of the Variance Components Performed from the Presented Scheme

    • 1)对于类型1,由图 1表 2的计算结果和上述分析可得,方案4和方案5由于顾及了待预测非公共点的观测误差,因此求得的被预测观测向量更接近真值。由于大地测量观测技术手段多样,此时应考虑各类观测值的精度问题,若直接采用TLS法进行计算可能并不恰当,因此应进一步采用验后估计方法,以实现随机模型的调整,此时待预测非公共点的数值能够得到更好的修正,进而获得的预测观测向量将更接近真值。

      2)从图 2表 3中的结果可知,当先验方差不合理时,未考虑方差分量估计的TLS法求得的二范数结果均被影响,但本文方法均能得到较合理的结果,表明本文方法能够实现观测数据不相关情形下各类观测值之间正确的权比求解,进一步测试了所提方法的可行性与有效性。

      3)当顾及观测值与系数矩阵之间的相关信息(即Qya=0),且对于不合适的先验随机模型以及模型预测时未考虑待预测非公共点的观测误差,本文利用最小二乘配置原理,结合方差-协方差估计理论实现了算法的迭代求解过程。从类型2中的计算结果可得,本文方法较其他方法能够得到较准确的参数估值结果,且被预测的观测向量的精度更佳,这可能归因于方差分量估计方法可以同时实现随机模型和待预测非公共点的协因数阵修正,因此获得的预测效果与精度得到了进一步改善。

    • 较为合理的参数估计结果以及各类观测误差的有效处理是提升回归模型预测效果的关键。针对传统的模型预测中未考虑待预测非公共点的观测误差和合适随机模型的获取问题,本文实现了模型参数估计与待预测非公共点误差改正数的求解问题,借助方差-协方差估计思想缓解了随机模型不准确的问题。另外,增加了相关观测问题的处理,在获取更好的模型预测结果的基础上,其理论更为严密。所采用的方差分量估计方法同时实现了随机模型与待预测非公共点的协因数阵的修正,参数估计结果与待预测非公共点误差改正数精度更高。本研究对于拓展与完善总体最小二乘理论与方法具有一定的理论与实用价值。

参考文献 (34)

目录

    /

    返回文章
    返回