留言板

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

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

非线性平差精度评定的自适应蒙特卡罗法

王乐洋 赵英文

王乐洋, 赵英文. 非线性平差精度评定的自适应蒙特卡罗法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
引用本文: 王乐洋, 赵英文. 非线性平差精度评定的自适应蒙特卡罗法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
WANG Leyang, ZHAO Yingwen. Adaptive Monte Carlo Method for Precision Estimation of Nonlinear Adjustment[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
Citation: WANG Leyang, ZHAO Yingwen. Adaptive Monte Carlo Method for Precision Estimation of Nonlinear Adjustment[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064

非线性平差精度评定的自适应蒙特卡罗法

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

国家自然科学基金 41874001

国家自然科学基金 41664001

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

国家重点研发计划 2016YFB0501405

详细信息
    作者简介:

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

  • 中图分类号: P207

Adaptive Monte Carlo Method for Precision Estimation of Nonlinear Adjustment

Funds: 

The National Natural Science Foundation of China 41874001

The National Natural Science Foundation of China 41664001

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, associate professor, specializes in geodetic inversion and geodetic data processing.E-mail::wleyang@163.com

  • 摘要: 针对现有的非线性平差精度评定理论中,蒙特卡罗法模拟次数的选择不具有客观性,无法对结果进行直接控制,以及没有同时考虑到平差参数估值、随机量改正数和单位权方差估值的有偏性等问题,把自适应蒙特卡罗法融入到非线性平差精度评定理论中。通过基于自适应蒙特卡罗法的估值偏差计算和参数估值协方差阵计算,设计了非线性平差精度评定一套理论完整的算法流程。基于对偶变量的思想,提出了参数估值偏差计算的对偶自适应蒙特卡罗法。直线拟合模型和椭圆拟合模型两个算例结果表明,非线性平差精度评定的自适应蒙特卡罗法能获得稳定且合理的精度评定结果,具有更强的适用性;对偶自适应蒙特卡罗法计算估值偏差的收敛速度更快,效率更高。
  • 图  1  自适应蒙特卡罗法计算流程图

    Figure  1.  The Flowchart of Adaptive Monte Carlo Method

    图  2  直线拟合模型参数估值均值及方差曲线图

    Figure  2.  The Curves of Means and Variances of Parameters Estimates of Straight Line Fitting Model

    图  3  椭圆拟合模型参数估值均值及方差曲线图

    Figure  3.  The Curves of Means and Variances of Parameters Estimations of Ellipse Fitting Model

    表  1  MC、AMC、AAMC 3种方法精度对比

    Table  1.   The Comparison for Accuracy of MC, AMC and AAMC

    均值 标准差 MC AMC AAMC
    均值未知 均值的标准差 $\sqrt {\sigma _z^2/M'} $ $s_{\bar z}^{\left( h \right)}$ $s_{{{\bar z}^ \pm }}^{\left( h \right)}$
    标准差的标准差 $\sqrt {\frac{{\sigma _z^2}}{{2\left( {M'-1} \right)}} + \frac{{\sigma _z^2}}{{8{{\left( {M'-1} \right)}^2}}}} $ su(z)(h) -
    均值已知 标准差的标准差 $\sqrt {\frac{{\sigma _z^2}}{{2M'}} + \frac{{\sigma _z^2}}{{8{{\left( {M'} \right)}^2}}}} $ su(z)(h) -
    下载: 导出CSV

    表  2  直线拟合模型参数偏差计算结果

    Table  2.   The Biases of Parameters Estimates of Straight Line Fitting Model

    参数 项目 ξ1 ξ2
    / 参数真值 2 1
    0.1I20(h=6) 参数估值 2.002 4 0.986 5
    MC 0.002 4(±1×10-5) -0.013 5(±7×10-5)
    AMC 0.002 2±2(3)×10-4 -0.013 8±13(20)×10-4
    AAMC 0.002 5±2(3)×10-5 -0.013 6±1(2)×10-4
    0.5I20(h=52) 参数估值 2.012 5 0.931 2
    MC 0.012 5(±2×10-5) -0.068 8(±2×10-4)
    AMC 0.012 7±2.5(2.5)×10-4 -0.069 1±1.5(1.5)×10-3
    AAMC 0.012 5±4(5)×10-5 -0.068 7±2(3)×10-4
    I20(h=328) 参数估值 2.025 7 0.858 4
    MC 0.025 7(±4×10-5) -0.141 6(±2×10-4)
    AMC 0.025 9±1.5(1.4)×10-4 -0.141 7±9.2(8.9)×10-4
    AAMC 0.025 7±4.4(4.1)×10-5 -0.141 3±2.5(2.4)×10-4
    下载: 导出CSV

    表  3  直线拟合模型参数标准差及协方差计算结果

    Table  3.   The Standard Deviations and Covariances of Parameters Estimates of Straight Line Fitting Model

    项目 方法 ${\hat \sigma _{{{\hat \xi }_1}}}$ ${\hat \sigma _{{{\hat \xi }_2}}}$ ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$
    0.1I20(h=6) FTT 0.078 2 0.484 9 -0.033 6
    MC 0.078 3(±8×10-6) 0.485 4(±5×10-5) -0.033 7
    AMC 0.078 7±4(23)×10-5 0.487 2±23(14)×10-4 -0.034 0
    0.5I20(h=52) FTT 0.178 0 1.102 0 -0.174 3
    MC 0.179 3(±2×10-5) 1.1074(±1×10-4) -0.176 7
    AMC 0.179 3±6(18)×10-5 1.107 7±25(11)×10-4 -0.1768
    I20(h=328) FTT 0.258 3 1.594 1 -0.367 0
    MC 0.261 6(±3×10-5) 1.608 6(±2×10-4) -0.376 4
    AMC 0.261 6±7(10)×10-5 1.608 2±25(6)×10-4 -0.376 3
    下载: 导出CSV

    表  4  直线拟合模型估值偏差计算结果

    Table  4.   The Biases of Estimates of Straight Line Fitting Model

    项目 估值 MC AMC(h=340) AAMC(h=340)
    ξ1 0.658 0 0.005 8(±2×10-6) 0.005 8±6(7)×10-5 0.005 8±2.0(2.5)×10-5
    ξ2 0.551 2 -0.013 1(±5×10-6) -0.013 0±1.8(2.0)×10-4 -0.013 1±5(6)×10-5
    σ02 1.538 6 -0.010 8(±1×10-5) -0.010 5±5.0(0.5)×10-4 -
    bv - 0.004 7 0.004 6 -
    下载: 导出CSV

    表  5  直线拟合模型参数估值标准差和协方差计算结果

    Table  5.   The Standard Deviations and Covariances of Parameters Estimates of Straight Line Fitting Model

    标准差或协方差 FTT MC AMC(h=64)
    ${\hat \sigma _{{{\hat \xi }_1}}}$ 0.119 5 0.124 9(±1×10-5) 0.124 9±3(11)×10-5
    ${\hat \sigma _{{{\hat \xi }_2}}}$ 0.349 1 0.360 3(±4×10-5) 0.360 1±2(3)×10-4
    ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$ -0.032 2 -0.035 2 -0.035 2
    下载: 导出CSV

    表  6  椭圆拟合模型估值偏差计算结果

    Table  6.   The Biases of Estimations of Ellipse Fitting Model

    参数 估值 MC AMC(h=297) AAMC(h=297)
    ξ1 -0.059 8 -0.010 8(±8×10-6) -0.010 9±3.0(3.2)×10-4 -0.010 6±6.3(6.5)×10-5
    ξ2 -0.194 2 0.001 4(±7×10-6) 0.001 3±2.9(3.0)×10-4 0.001 4±4.9(5.1)×10-5
    ξ3 13.108 7 0.092 5(±9×10-6) 0.092 3±3.8(4.0)×10-4 0.092 5±9.5(9.6)×10-5
    ξ4 11.513 1 0.048 8(±9×10-6) 0.049 2±3.6(3.6)×10-4 0.048 8±7.1(6.8)×10-5
    σ02 1.046 4 -0.009 5(±9×10-6) -0.009 5±3.5(0.4)×10-4 -
    bv - 0.022 6 0.0226 -
    下载: 导出CSV

    表  7  椭圆拟合模型参数估值标准差和协方差计算结果

    Table  7.   The Standard Deviations and Covariances of Parameters Estimations of Ellipse Fitting Model

    参数 FTT MC AMC(h=196)
    ${\hat \sigma _{{{\hat \xi }_1}}}$ 0.534 9 0.546 8(±5×10-5) 0.547 4±3.0(2.8)×10-4
    ${\hat \sigma _{{{\hat \xi }_2}}}$ 0.515 4 0.520 3(±5×10-5) 0.520 4±2.7(2.6)×10-4
    ${\hat \sigma _{{{\hat \xi }_3}}}$ 0.632 7 0.682 3(±7×10-5) 0.682 4±5.0(3.4)×10-4
    ${\hat \sigma _{{{\hat \xi }_4}}}$ 0.606 1 0.624 6(±6×10-5) 0.624 7±3.9(3.2)×10-4
    ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$ -0.003 3 -0.004 6 -0.004 7
    ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_3}}}$ -0.081 6 -0.095 1 -0.095 0
    ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_4}}}$ 0.017 2 0.019 1 0.018 6
    ${\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_3}}}$ 0.022 9 0.027 9 0.028 8
    ${\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_4}}}$ -0.148 3 -0.151 2 -0.151 1
    ${\hat \sigma _{{{\hat \xi }_3}{{\hat \xi }_4}}}$ -0.108 7 -0.108 4 -0.109 1
    下载: 导出CSV
  • [1] 陶本藻.形变反演模型的非线性平差[J].武汉大学学报·信息科学版, 2001, 26(6):504-508 http://ch.whu.edu.cn/CN/abstract/abstract5222.shtml

    Tao Benzao. Non-linear Adjustment of Deformation Inversion Model[J]. Geomatics and Information Science of Wuhan University, 2001, 26(6):504-508 http://ch.whu.edu.cn/CN/abstract/abstract5222.shtml
    [2] 李朝奎.非线性模型空间测量数据处理理论及其应用[D].长沙: 中南大学, 2001

    Li Chaokui. Theory and Application of Data Processing in Space of Nonlinear Models[D]. Changsha: Central South University, 2001
    [3] 刘国林.非线性最小二乘与测量平差[M].北京:测绘出版社, 2002

    Liu Guolin. Nonlinear Least Squares and Surveying Adjustment[M]. Beijing:Surveying and Mapping Press, 2002
    [4] 王新洲.非线性模型参数估计理论与应用[M].武汉:武汉大学出版社, 2002

    Wang Xinzhou. The Theory and Application of Parameter Estimation of Nonlinear Model[M]. Wuhan:Wuhan University Press, 2002
    [5] 张松林.非线性半参数模型最小二乘估计理论及应用研究[D].武汉: 武汉大学, 2003 http://cdmd.cnki.com.cn/Article/CDMD-10486-2006031323.htm

    Zhang Songlin. The Theoretical and Application Research on Nonlinear Semiparametric Model[D]. Wuhan: Wuhan University, 2003 http://cdmd.cnki.com.cn/Article/CDMD-10486-2006031323.htm
    [6] 唐利民.非线性最小二乘的不适定性及算法研究[D].长沙: 中南大学, 2011

    Tang Limin. Research on the Ill-posed and Solving Methods of Nonlinear Least Squares Problem[D]. Changsha: Central South University, 2011
    [7] 薛树强, 杨元喜, 党亚民.测距定位方程非线性平差的封闭牛顿迭代公式[J].测绘学报, 2014, 43(8):771-777 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QKC20142014092200007705

    Xue Shuqiang, Yang Yuanxi, Dang Yamin. A Closed-form of Newton Iterative Formula for Nonlinear Adjustment of Distance Equations[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(8):771-777 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QKC20142014092200007705
    [8] Xue S, Yang Y, Dang Y. A Closed-form of Newton Method for Solving Over-Determined Pseudo-Distance Equations[J]. Journal of Geodesy, 2014, 88(5):441-448 doi:  10.1007/s00190-014-0695-y
    [9] Box M J. Bias in Nonlinear Estimation[J]. Journal of the Royal Statistical Society, Series B (Metho-dological), 1971, 33(2):171-201 doi:  10.1111/rssb.1971.33.issue-2
    [10] Teunissen P J G, Knickmeyer E H. Nonlinearity and Least Squares[J]. CISM Journal ASCGC, 1988, 42(4):321-330 http://d.old.wanfangdata.com.cn/Periodical/xtgcydzjs201803008
    [11] Teunissen P J G. First and Second Moments of Non-Linear Least-Squares Estimators[J]. Bulletin Géodésique, 1989, 63(3):253-262 doi:  10.1007/BF02520475
    [12] 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
    [13] Schaffrin B, Snow K. Total Least-Squares Regularization of Tykhonov Type and an Ancient Racetrack in Corinth[J]. Linear Algebra and Its Applications, 2010, 432(8):2061-2076 doi:  10.1016/j.laa.2009.09.014
    [14] Fang X. Weighted Total Least Squares Solutions for Applications in Geodesy[D]. Hanover: Leibniz University of Hanover, 2011
    [15] Shen Y, Li B, 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
    [16] Xu P, Liu J, 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
    [17] 姚宜斌, 孔建.顾及设计矩阵随机误差的最小二乘组合新解法[J].武汉大学学报·信息科学版, 2014, 39(9):1028-1032 http://ch.whu.edu.cn/CN/abstract/abstract3065.shtml

    Yao Yibin, Kong Jian. A New Combined LS Method Considering Random Errors of Design Matrix[J]. Geomatics and Information Science of Wuhan University, 2014, 39(9):1028-1032 http://ch.whu.edu.cn/CN/abstract/abstract3065.shtml
    [18] Tong X, Jin Y, Zhang S, et al. Bias-Corrected Weighted Total Least-Squares Adjustment of Condition Equations[J]. Journal of Surveying Engineering, 2014, 141(2):04014013 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=26f7db7fb994669a82ad317964b29beb
    [19] Zeng W, Liu J, Yao Y. On Partial Errors-in-Variables Models with Inequality Constraints of Parameters and Variables[J]. Journal of Geodesy, 2015, 89(2):111-119 doi:  10.1007/s00190-014-0775-z
    [20] Fang X. Weighted 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
    [21] Schaffrin B. Adjusting the Errors-in-Variables Model: Linearized Least-Squares vs. Nonlinear Total Least-Squares[C]//Ⅷ Hotine-Marussi Symposium on Mathematical Geodesy. Switzerland: Springer International Publishing, 2015: 301-307
    [22] Fang X, Wu Y. On the Errors-in-Variables Model with Equality and Inequality Constraints for Selected Numerical Examples[J]. Acta Geodaetica et Geophysica, 2016, 51(3):515-525 doi:  10.1007/s40328-015-0141-5
    [23] 王乐洋, 余航, 陈晓勇. Partial EIV模型的解法[J].测绘学报, 2016, 45(1):22-29 http://d.old.wanfangdata.com.cn/Periodical/chxb201601004

    Wang Leyang, Yu Hang, Chen Xiaoyong. An Algorithm for Partial EIV Model[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(1):22-29 http://d.old.wanfangdata.com.cn/Periodical/chxb201601004
    [24] 王乐洋, 赵英文, 陈晓勇, 等.多元总体最小二乘问题的牛顿解法[J].测绘学报, 2016, 45(4):411-417 http://d.old.wanfangdata.com.cn/Periodical/chxb201604005

    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 http://d.old.wanfangdata.com.cn/Periodical/chxb201604005
    [25] 曾文宪, 方兴, 刘经南, 等.通用EIV平差模型及其加权整体最小二乘估计[J].测绘学报, 2016, 45(8):890-894 http://d.old.wanfangdata.com.cn/Periodical/chxb201608002

    Zeng Wenxian, Fang Xing, Liu Jingnan, et al. Weighted Total Least Squares of Universal EIV Adjustment Model[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(8):890-894 http://d.old.wanfangdata.com.cn/Periodical/chxb201608002
    [26] 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
    [27] Amiri-Simkooei A R, Zangeneh-Nejad F, Asgari J. On the Covariance Matrix of Weighted Total Least-Squares Estimates[J]. Journal of Surveying Engineering, 2016, 142(3):04015014 doi:  10.1061/(ASCE)SU.1943-5428.0000153
    [28] JCGM 101: 2008. Evaluation of Measurement Data-Supplement 1 to the "Guide to the Expression of Uncertainty in Measurement"-Propagation of Distributions Using a Monte Carlo Method[S]. Sèvres: JCGM, 2008
    [29] Sneeuy N, Krumm F, Roth M. Adjustment Theory[M]. Stuttgart:University of Stuttgart, 2015
    [30] Ratkowsky D A. Nonlinear Regression Modeling[M]. New York:Marcel Dekker, 1983
    [31] Robert C, Casella G. Monte Carlo Statistical Methods[M]. Berlin:Springer Science & Business Media, 2013
    [32] 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
  • [1] 杨鹏飞, 赵庆志, 苏静, 姚宜斌.  黄土高原地区PWV影响因素分析及精度评定 . 武汉大学学报 ● 信息科学版, 2022, 47(9): 1470-1478. doi: 10.13203/j.whugis20200390
    [2] 王乐洋, 邹传义.  乘性误差模型参数估计及精度评定的Sterling插值方法 . 武汉大学学报 ● 信息科学版, 2022, 47(2): 219-225. doi: 10.13203/j.whugis20200052
    [3] 王乐洋, 邹传义.  PEIV模型参数估计理论及其应用研究进展 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1273-1283, 1297. doi: 10.13203/j.whugis20200312
    [4] 曾文宪, 刘泽邦, 方兴, 李玉兵.  通用EIV平差模型的线性化估计算法 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1284-1290. doi: 10.13203/j.whugis20200243
    [5] 聂建亮, 刘晓云, 田婕, 李秀明, 赵大江, 黄功文, 张海平.  利用自适应水准网动态平差建立山东垂直运动速度模型 . 武汉大学学报 ● 信息科学版, 2020, 45(4): 620-625. doi: 10.13203/j.whugis20180296
    [6] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [7] 谢建, 龙四春, 李黎, 李博超.  不等式约束加权整体最小二乘的凝聚函数法 . 武汉大学学报 ● 信息科学版, 2018, 43(10): 1526-1530. doi: 10.13203/j.whugis20160507
    [8] 刘超, 王彬, 赵兴旺, 余学祥.  三维坐标转换的高斯-赫尔默特模型及其抗差解法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1320-1327. doi: 10.13203/j.whugis20160348
    [9] 王君刚, 陈俊平, 王解先, 章洁君, 宋雷.  对流层经验改正模型在中国区域的精度评估 . 武汉大学学报 ● 信息科学版, 2016, 41(12): 1656-1663. doi: 10.13203/j.whugis20140696
    [10] 薛树强, 杨元喜.  抗差高斯-雅克比组合平差法 . 武汉大学学报 ● 信息科学版, 2015, 40(7): 932-937. doi: 10.13203/j.whugis20130247
    [11] 姚宜斌, 刘强, 江国焰, 张良.  华北地区应变率及其精度评定 . 武汉大学学报 ● 信息科学版, 2015, 40(10): 1317-1323. doi: 10.13203/j.whugis20140001
    [12] 姚宜斌, 孔建.  顾及设计矩阵随机误差的最小二乘组合新解法 . 武汉大学学报 ● 信息科学版, 2014, 39(9): 1028-1032. doi: 10.13203/j.whugis20130030
    [13] 游为, 范东明, 傅淑娟.  同伦函数与填充函数相结合的非线性最小二乘平差模型 . 武汉大学学报 ● 信息科学版, 2010, 35(2): 185-188.
    [14] 吴江飞, 黄珹.  非线性自适应抗差滤波定轨算法 . 武汉大学学报 ● 信息科学版, 2008, 33(2): 187-190.
    [15] 张勤, 陶本藻.  基于同伦法的非线性最小二乘平差统一模型 . 武汉大学学报 ● 信息科学版, 2004, 29(8): 708-710.
    [16] 陶本藻.  形变反演模型的非线性平差 . 武汉大学学报 ● 信息科学版, 2001, 26(6): 504-508.
    [17] 王新洲.  非线性模型平差中单位权方差的估计 . 武汉大学学报 ● 信息科学版, 2000, 25(4): 358-361.
    [18] 袁修孝.  GPS辅助光束法平差的理论精度 . 武汉大学学报 ● 信息科学版, 1998, 23(4): 394-398.
    [19] 於宗俦, 于正林.  分组平差与精度评定 . 武汉大学学报 ● 信息科学版, 1988, 13(3): 32-46.
    [20] 孔祥元.  用力学原理进行工测三维测边网平差和精度评定 . 武汉大学学报 ● 信息科学版, 1986, 11(3): 126-135.
  • 加载中
图(3) / 表(7)
计量
  • 文章访问数:  1620
  • HTML全文浏览量:  103
  • PDF下载量:  252
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-08-16
  • 刊出日期:  2019-02-05

非线性平差精度评定的自适应蒙特卡罗法

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

    国家自然科学基金 41874001

    国家自然科学基金 41664001

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

    国家重点研发计划 2016YFB0501405

    作者简介:

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

  • 中图分类号: P207

摘要: 针对现有的非线性平差精度评定理论中,蒙特卡罗法模拟次数的选择不具有客观性,无法对结果进行直接控制,以及没有同时考虑到平差参数估值、随机量改正数和单位权方差估值的有偏性等问题,把自适应蒙特卡罗法融入到非线性平差精度评定理论中。通过基于自适应蒙特卡罗法的估值偏差计算和参数估值协方差阵计算,设计了非线性平差精度评定一套理论完整的算法流程。基于对偶变量的思想,提出了参数估值偏差计算的对偶自适应蒙特卡罗法。直线拟合模型和椭圆拟合模型两个算例结果表明,非线性平差精度评定的自适应蒙特卡罗法能获得稳定且合理的精度评定结果,具有更强的适用性;对偶自适应蒙特卡罗法计算估值偏差的收敛速度更快,效率更高。

English Abstract

王乐洋, 赵英文. 非线性平差精度评定的自适应蒙特卡罗法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
引用本文: 王乐洋, 赵英文. 非线性平差精度评定的自适应蒙特卡罗法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
WANG Leyang, ZHAO Yingwen. Adaptive Monte Carlo Method for Precision Estimation of Nonlinear Adjustment[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
Citation: WANG Leyang, ZHAO Yingwen. Adaptive Monte Carlo Method for Precision Estimation of Nonlinear Adjustment[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 206-213, 220. doi: 10.13203/j.whugis20170064
  • 测量领域中常常遇到非线性模型,非线性平差一直是国内外研究的热点问题。相比于参数估计理论,非线性平差精度评定理论有待进一步发展。以非线性最小二乘法为基础,文献[1-8]从不同的应用背景研究了非线性模型参数估计的平差算法;文献[9]首先给出了非线性最小二乘参数估值的偏差公式;文献[10-11]给出了非线性最小二乘估值前两阶矩的计算公式;文献[3-4]分别讨论了顾及二次项的参数估值协方差阵计算和非线性最小二乘的统计性质与精度评定。近10年来,顾及函数模型系数矩阵误差的总体最小二乘法[12]得到了广泛的研究。以EIV(errors-in-variables)模型[13]为基础,总体最小二乘法本质上为一种非线性平差方法。文献[13-26]给出了不同的总体最小二乘参数估计方法。在总体最小二乘精度评定方面,参数估值一阶近似协方差阵可基于EIV模型[15, 17, 27]、多元EIV模型[24]、等式或不等式约束的EIV模型[20, 22]、高斯-赫尔默特(Gauss-Helmert,GH)模型[13, 25]推导得出。文献[16]基于partial EIV模型进行了总体最小二乘估值偏差以及一阶参数估值协方差阵;文献[19]给出不等式约束下的一阶近似协方差阵;文献[18]推导了基于条件方程的总体最小二乘平差估值,的偏差及协方差阵计算公式。依据总体最小二乘估值文献[15]使用蒙特卡罗(Monte Carlo,MC)法计算单位权方差估值的偏差,并利用偏差改正的单位权方差估值进一步计算参数估值的均方误差矩阵;文献[27]针对模拟数据使用蒙特卡罗法计算总体最小二乘参数估值协方差阵。总体来看,非线性平差参数估值的协方差阵公式主要局限于一阶精度,并且蒙特卡罗法在非线性平差精度评定中的研究较少;此外尚没有文献在提高精度评定结果的同时综合考虑所有非线性平差估值有偏性。

    为了进一步发展非线性平差精度评定理论,获取稳定的蒙特卡罗模拟结果且避免事先选择模拟次数,本文引入自适应蒙特卡罗(adaptive Monte Carlo,AMC)法[28]进行非线性平差精度评定。先给出3自适应蒙特卡罗法进行精度评定的流程,并提出计算参数估值偏差的对偶自适应蒙特卡罗(antithetic adaptive Monte Carlo, AAMC)法。然后通过算例验证了本文方法进行非线性平差精度评定的可行性和有效性。

    • 作为一般的非线性平差模型,高斯-赫尔默特模型定义为隐函数[29]

      $$ f\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{v}},\mathit{\boldsymbol{\xi }}} \right) = 0 $$ (1)

      其随机模型可假设为:

      $$ \mathit{\boldsymbol{v}} \sim N\left( {0,\sigma _0^2\mathit{\boldsymbol{Q}}} \right) $$ (2)

      式中,ln×1的观测向量;v为误差向量;ξm×1的未知参数向量;f(·)为非线性映射;σ02为未知单位权方差,P=Q-1QP为观测向量协因数阵和权阵。

      线性化后的模型(1)的参数估值为[13-14]

      $$ \mathit{\boldsymbol{\hat \xi }} = - {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{BQ}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{BQ}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{w}} + {{\mathit{\boldsymbol{\hat \xi }}}^0} $$ (3)

      式中,$\mathit{\boldsymbol{A}} = {\left. {\frac{{\partial f}}{{\partial \mathit{\boldsymbol{\xi }}}}} \right|_{\mathit{\boldsymbol{l}}-{\mathit{\boldsymbol{v}}^0};{\mathit{\boldsymbol{\xi }}^0}}}$;$\mathit{\boldsymbol{B}} = {\left. {\frac{{\partial f}}{{\partial \mathit{\boldsymbol{v}}}}} \right|_{\mathit{\boldsymbol{l}}-{\mathit{\boldsymbol{v}}^0}, {\mathit{\boldsymbol{\xi }}^0}}}$;w=f(lv0, ξ0)-Bv0,上标“0”表示近似值。

      式(3)的迭代过程见文献[13-14]。参数估值的一阶近似协方差阵(first-order truncated Taylor ekpansion, FTT)为[13]:

      $$ {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat \xi }}}} = \hat \sigma _0^2{\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{BQ}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} \right)}^{ - 1}}\mathit{\boldsymbol{A}}} \right)^{ - 1}} $$ (4)

      式中,AB在$\mathit{\boldsymbol{\hat \xi }}$和$\mathit{\boldsymbol{\hat v}}$处取值,$\hat \sigma _0^2 = {\mathit{\boldsymbol{\hat v}}^{\text{T}}}\mathit{\boldsymbol{P\hat v/}}\left( {m-n} \right)$为单位权方差估值。

    • 通过分批次增加蒙特卡罗的模拟次数, 直到输出结果满足数值容差,达到统计意义上稳定性的方法为自适应蒙特卡罗法。其计算流程如图 1所示。

      图  1  自适应蒙特卡罗法计算流程图

      Figure 1.  The Flowchart of Adaptive Monte Carlo Method

      图 1中,x=g(l)为已知的非线性函数,x为输出量(标量),ll的均值,Dll的协方差阵,g(·)为非线性映射,h表示批次,M为每批次的模拟次数,δ为给定的数值容差,x(h)u(x(h))分别为第h批次中x的均值与标准差,sx(h)su(x)(h)为均值的标准差以及标准差的标准差,Ex(h)Dx(h)分别为h×M个输出量的均值和方差,max(sx(h), su(x)(h))表示sx(h)su(x)(h)中的最大值。

    • 考虑到非线性平差估计结果的有偏性,本文进一步研究非线性平差精度评定的自适应蒙特卡罗法。

      令式(3)的迭代计算过程表示为非线性映射g(·),参数估值、改正数和单位权方差估值向量化为$\mathit{\boldsymbol{\hat x = }}{\left[{{{\mathit{\boldsymbol{\hat \xi }}}^{\text{T}}}\;\;\;{{\mathit{\boldsymbol{\hat v}}}^{\text{T}}}\;\;\;\;\hat \sigma _0^2} \right]^{\text{T}}}$,则估值向量和观测值向量之间的关系为$\mathit{\boldsymbol{\hat x}} = g\left( \mathit{\boldsymbol{l}} \right)$。在计算估值偏差时,观测向量l的均值和协方差阵用其平差值$\mathit{\boldsymbol{\hat l}} = \mathit{\boldsymbol{l}}-\mathit{\boldsymbol{\hat v}}$和后验协方差阵${\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{l}}} = \hat \sigma _0^2\mathit{\boldsymbol{Q}}$近似表示。除每批次所生成的估值向量不同外,第h批次均值和均值标准差(h>1)的计算具有规律性[28],即为:

      $$ {{\mathit{\boldsymbol{\hat x}}}^{\left( h \right)}} = \frac{1}{M}\sum\limits_{i = 1}^M {\mathit{\boldsymbol{x}}_{\left( i \right)}^{\left( h \right)}} $$ (5)
      $$ s_{{{\mathit{\boldsymbol{\bar x}}}_k}}^{\left( h \right)} = \sqrt {\frac{1}{{h\left( {h - 1} \right)}}\sum\limits_{i = 1}^h {{{\left( {\mathit{\boldsymbol{\bar x}}_k^{\left( i \right)} - \frac{1}{h}\sum\limits_{j = 1}^h {\mathit{\boldsymbol{\bar x}}_k^{\left( j \right)}} } \right)}^2}} } $$ (6)

      式中,k=1, 2…m+n+1;x(h)表示h批次的估值向量均值;$s_{{{\mathit{\boldsymbol{\bar x}}}_k}}^{\left( h \right)}$表示前h个均值第k个分量的标准差;x(i)(h)表示第h批次的第i个估值向量随机数;xk(h)表示第h批次的估值向量第k个分量的均值。为减少每批次估值向量的存储,本文给出每批次均值与最终均值的递推公式为:

      $$ \mathit{\boldsymbol{E}}_\mathit{\boldsymbol{x}}^{\left( h \right)} = \frac{{h - 1}}{h}\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{x}}^{\left( {h - 1} \right)} + \frac{1}{h}{{\mathit{\boldsymbol{\bar x}}}^{\left( h \right)}} $$ (7)

      根据式(7),参数估值、改正数和单位权方差估值的均值分别为${\mathit{\boldsymbol{E}}_{\hat \xi }}$、${\mathit{\boldsymbol{E}}_{\hat v}}$和${\mathit{\boldsymbol{E}}_{\hat \sigma _0^2}}$,则其相应的偏差和参数估值偏差百分比分别为:

      $$ {\mathit{\boldsymbol{b}}_{\hat \xi }} = {\mathit{\boldsymbol{E}}_{\hat \xi }} - \mathit{\boldsymbol{\hat \xi }} $$ (8)
      $$ {\mathit{\boldsymbol{b}}_{\hat v}} = {\mathit{\boldsymbol{E}}_{\hat v}} $$ (9)
      $$ {b_{\hat \sigma _0^2}} = {\mathit{\boldsymbol{E}}_{\hat \sigma _0^2}} - \hat \sigma _0^2 $$ (10)
      $$ {\theta _{{b_{{{\hat \xi }_i}}}}} = \left( {{\mathit{\boldsymbol{b}}_{{{\mathit{\boldsymbol{\hat \xi }}}_i}}}/{{\hat \xi }_i}} \right) \times 100\% \left( {i = 1,2 \cdots m} \right) $$ (11)

      式中,${\mathit{\boldsymbol{b}}_{\hat \xi }}$、${\mathit{\boldsymbol{b}}_{\hat v}}$、${b_{\hat \sigma _0^2}}$和${\theta _{{b_{\hat \xi i}}}}$分别表示参数估值偏差、改正数偏差、单位权方差估值偏差和参数估值偏差百分比。

      根据偏差改正后的估值向量$\mathit{\boldsymbol{\hat x}} = \mathit{\boldsymbol{\hat x}}- {\left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{b}}_{\hat \xi }^{\text{T}}}&{\mathit{\boldsymbol{b}}_{\tilde v}^{\text{T}}}&{{b_{\tilde \sigma _0^2}}} \end{array}} \right]^{\text{T}}} = {\left[{\begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde \xi }}}^{\text{T}}}}&{{{\mathit{\boldsymbol{\tilde v}}}^{\text{T}}}}&{\tilde \sigma _0^2} \end{array}} \right]^{\text{T}}}$,重新构建观测向量l的均值$\mathit{\boldsymbol{\bar l'}} = \mathit{\boldsymbol{l}}-\mathit{\boldsymbol{\tilde v}}$和协方差阵${\mathit{\boldsymbol{D'}}_l} = \tilde \sigma _0^2\mathit{\boldsymbol{Q}}$。第h批次参数估值的协方差阵、标准差及标准差的标准差为:

      $$ \mathit{\boldsymbol{U}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( h \right)}}} \right) = \frac{1}{M}\sum\limits_{i = 1}^M {\left( {\mathit{\boldsymbol{\hat \xi }}_{\left( i \right)}^{\left( h \right)} - \mathit{\boldsymbol{\tilde \xi }}} \right){{\left( {\mathit{\boldsymbol{\hat \xi }}_{\left( i \right)}^{\left( h \right)} - \mathit{\boldsymbol{\tilde \xi }}} \right)}^{\rm{T}}}} $$ (12)
      $$ \mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{\hat \xi }}_{i }^{\left( h \right)}} \right) = \sqrt {\mathit{\boldsymbol{U}}{{\left( {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( h \right)}}} \right)}_{ii}}} ,i = 1,2 \cdots m $$ (13)
      $$ \mathit{\boldsymbol{s}}_{\mathit{\boldsymbol{u}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}_k}} \right)}^{\left( h \right)} = \sqrt {\frac{1}{{h\left( {h - 1} \right)}}\sum\limits_{i = 1}^h {\left( {\mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{\hat \xi }}_k^{\left( i \right)}} \right)}\right. } - \frac{1}{h}\sum\limits_{j = 1}^h {\mathit{\boldsymbol{u}}{{\left( {\mathit{\boldsymbol{\hat \xi }}_k^{\left( j \right)}} \right)}^2}} } $$ (14)

      式中,k=1, 2…m;$\mathit{\boldsymbol{U}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( h \right)}}} \right)$、$\mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{\hat \xi }}_i^{\left( h \right)}} \right)$和$s_{\mathit{\boldsymbol{u}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}_k}} \right)}^{\left( h \right)}$分别表示第h批次参数估值的协方差阵,参数估值第k个分量的标准差及前h个标准差的标准差。

      如均值递推公式、协方差阵递推公式为:

      $$ \mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat \xi }}}^{\left( h \right)} = \frac{{h - 1}}{h}\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{\hat \xi }}}^{\left( {h - 1} \right)} + \frac{1}{h}\mathit{\boldsymbol{U}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( h \right)}}} \right) $$ (15)

      根据以上内容,非线性平差精度评定的自适应蒙特卡罗法总结为算法1。其具体步骤为:

      1) h=1,据lDl,随机生成Ml(i),据式(3)、式(5)计算x(1)

      2) h=h+1,重复计算x(h)sx(h)

      3) 当2max(sx1(h), sx2(h)sxm+n+1(h)) < δ时,终止循环,据式(7)-(11)计算${\mathit{\boldsymbol{E}}_x} = \mathit{\boldsymbol{E}}_x^{\left( h \right)}$、${\mathit{\boldsymbol{b}}_{\hat \xi }}$、${\mathit{\boldsymbol{b}}_{\hat v}}$、${b_{\hat \sigma _0^2}}$和${\theta _{{b_{{{\hat \xi }_i}}}}}$;

      4) h=1,据l′和Dl,随机生成M′组l(i),据式(3)、式(12)计算$\mathit{\boldsymbol{U}}\left( {{{\mathit{\boldsymbol{\hat \xi }}}^{\left( 1 \right)}}} \right)$;

      5) h=h+1,重复4)计算$\mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{\hat \xi }}_i^{\left( h \right)}} \right)$、$s_{u\left( {{{\hat \xi }_i}} \right)}^{\left( h \right)}$;

      6) 当$2\max \left( {s_{u\left( {{{\hat \xi }_1}} \right)}^{\left( h \right)}, s_{u\left( {{{\hat \xi }_2}} \right)}^{\left( h \right)} \cdots s_{u\left( {{{\hat \xi }_m}} \right)}^{\left( h \right)}} \right) < \delta $时,终止循环,计算${\mathit{\boldsymbol{D}}_{\hat \xi }} = \mathit{\boldsymbol{D}}_{\hat \xi }^{\left( h \right)}$。

      蒙特卡罗法和自适应蒙特卡罗法模拟次数的选择不同。蒙特卡罗法模拟次数M′事先给定(如1×106),不分批次计算。算法1每批次蒙特卡罗模拟次数M可为100/(1-p)(p为包含概率,如取0.95)和1×104之间的最大值[28],本文取1×104图 1中的判断条件可表示为max(sx(h), su(x)(h)) < δ/2。因为sx(h)su(x)(h)为均值的标准以及与标准差的标准差,δ/2可认为是均值和标准差的最低精度,所以自适应蒙特卡罗法能同时给出输出结果及大致精度。本文中δ取值为2倍的均值或标准差的精度ε,即δ=2ε

    • 根据偏差的非对称性度量方法[28],在自适应蒙特卡罗法中设计根据一组随机误差生成两组呈完全负相关的观测向量

      $$ {\mathit{\boldsymbol{l}}^ - } = \mathit{\boldsymbol{\bar l}} - {\rm{d}}\mathit{\boldsymbol{l}},{\mathit{\boldsymbol{l}}^ + } = \mathit{\boldsymbol{\bar l}} + {\rm{d}}\mathit{\boldsymbol{l}} $$ (16)

      式中,dl表示随机误差向量;ll+即为呈负相关的观测向量。

      假设ll+对应的参数估值为${\mathit{\boldsymbol{\hat \xi }}^-}$和${\mathit{\boldsymbol{\hat \xi }}^ + }$,如果其每个分量之间的相关系数corr($\mathit{\boldsymbol{\hat \xi }}_i^-, \mathit{\boldsymbol{\hat \xi }}_i^ + $)(i=1,2…m)都小于0,则均值(${\mathit{\boldsymbol{\hat \xi }}^-} + {\mathit{\boldsymbol{\hat \xi }}^ + }$)/2每个分量的方差要小于两个参数估值相互独立情况下的方差,其中corr(a, b)表示随机量ab之间的相关系数。鉴于蒙特卡罗模拟方差缩减技术中的对偶变量法[31]原理,本文把采用式(16)生成随机样本点的自适应蒙特卡罗法称为对偶自适应蒙特卡罗法。

      假设根据式(16)生成M/2对负相关的观测量ll+,则相应的参数估值矩阵为:

      $$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varXi} }}}^ - } = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat \xi }}_{\left( 1 \right)}^ - }&{\mathit{\boldsymbol{\hat \xi }}_{\left( 2 \right)}^ - }& \cdots &{\mathit{\boldsymbol{\hat \xi }}_{\left( {M/2} \right)}^ - } \end{array}} \right]^{\rm{T}}} $$ (17)
      $$ {{\mathit{\boldsymbol{ \boldsymbol{\hat \varXi} }}}^ + } = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat \xi }}_{\left( 1 \right)}^ + }&{\mathit{\boldsymbol{\hat \xi }}_{\left( 2 \right)}^ + }& \cdots &{\mathit{\boldsymbol{\hat \xi }}_{\left( {M/2} \right)}^ + } \end{array}} \right]^{\rm{T}}} $$ (18)

      以及参数估值分量之间的相关系数为:

      $$ \begin{array}{*{20}{c}} {{\rho _j} = }\\ {\frac{{\sum\limits_{i = 1}^{M/2} {\left( {\mathit{\hat \Xi }_{ij}^ - - \frac{1}{{M/2}}\sum\limits_{i = 1}^{M/2} {\mathit{\hat \Xi }_{ij}^ - } } \right)} \left( {\mathit{\hat \Xi }_{ij}^ + - \frac{1}{{M/2}}\sum\limits_{i = 1}^{M/2} {\mathit{\hat \Xi }_{ij}^ + } } \right)}}{{\sqrt {\sum\limits_{i = 1}^{M/2} {{{\left( {\mathit{\hat \Xi }_{ij}^ - - \frac{1}{{M/2}}\sum\limits_{i = 1}^{M/2} {\mathit{\hat \Xi }_{ij}^ - } } \right)}^2}} \sum\limits_{i = 1}^{M/2} {{{\left( {\mathit{\hat \Xi }_{ij}^ + - \frac{1}{{M/2}}\sum\limits_{i = 1}^{M/2} {\mathit{\hat \Xi }_{ij}^ + } } \right)}^2}} } }}} \end{array} $$ (19)

      式中,j=1,2…m

      对偶自适应蒙特卡罗法对应的计算公式为:

      $$ {\left( {{{\mathit{\boldsymbol{\bar \xi }}}^ \pm }} \right)^{\left( h \right)}} = \frac{1}{M}\sum\limits_{i = 1}^{M/2} {\left( {{{\left( {\mathit{\boldsymbol{\hat \xi }}_{\left( i \right)}^ - } \right)}^{\left( h \right)}} + {{\left( {\mathit{\boldsymbol{\hat \xi }}_{\left( i \right)}^ + } \right)}^{\left( h \right)}}} \right)} $$ (20)
      $$ s_{\mathit{\boldsymbol{\bar \xi }}_k^ \pm }^{\left( h \right)} = \sqrt {\frac{1}{{h\left( {h - 1} \right)}}\sum\limits_{i = 1}^h {{{\left( {{{\left( {\mathit{\boldsymbol{\bar \xi }}_k^ \pm } \right)}^{\left( i \right)}} - \frac{1}{h}\sum\limits_{j = 1}^h {{{\left( {\mathit{\boldsymbol{\bar \xi }}_k^ \pm } \right)}^{\left( j \right)}}} } \right)}^2}} } $$ (21)
      $$ \mathit{\boldsymbol{E}}_{{{\mathit{\boldsymbol{\hat \xi }}}^ \pm }}^{\left( h \right)} = \frac{{h - 1}}{h}\mathit{\boldsymbol{E}}_{{{\mathit{\boldsymbol{\hat \xi }}}^ \pm }}^{\left( {h - 1} \right)} + \frac{1}{h}{\left( {{{\mathit{\boldsymbol{\bar \xi }}}^ \pm }} \right)^{\left( h \right)}} $$ (22)

      式中,k=1, 2…m;(${\mathit{\boldsymbol{\bar \xi }}^ \pm }$)(h)表示由第h批次${\mathit{\boldsymbol{\hat \xi }}^-}$和${\mathit{\boldsymbol{\hat \xi }}^ + }$计算的参数估值均值;$s_{\bar \xi _k^ \pm }^{\left( h \right)}$表示前h个批次所有参数估值均值(${\mathit{\boldsymbol{\bar \xi }}^ \pm }$)(h)的标准差;$\mathit{\boldsymbol{E}}_{{{\hat \xi }^ \pm }}^{\left( h \right)}$表示由h×(M/2)个${\mathit{\boldsymbol{\hat \xi }}^-}$和h×(M/2)个${\mathit{\boldsymbol{\hat \xi }}^ + }$计算的参数估值均值。

      参数估值偏差及其偏差百分比公式为:

      $$ {\mathit{\boldsymbol{b}}_{{{\mathit{\boldsymbol{\hat \xi }}}^ \pm }}} = {\mathit{\boldsymbol{E}}_{{{\mathit{\boldsymbol{\hat \xi }}}^ \pm }}} - \mathit{\boldsymbol{\hat \xi }} $$ (23)
      $$ {\theta _{{b_{\hat \xi _i^ \pm }}}} = \left( {{\mathit{\boldsymbol{b}}_{\hat \xi _i^ \pm }}/\mathit{\boldsymbol{\hat \xi }}_i^ \pm } \right) \times 100\% \left( {i = 1,2 \cdots m} \right) $$ (24)

      式中,${\mathit{\boldsymbol{b}}_{{{\hat \xi }^ \pm }}}$表示参数估值偏差;${\theta _{{b_{\hat \xi _i^ \pm }}}}$表示参数估值各个分量的偏差百分比。

      对偶自适应蒙特卡罗法计算参数估值偏差的具体流程总结为算法2。具体步骤为:

      1) 据lDl,随机生成M/2组dl(i),构建l(i)l(i)+,据式(3)、式(19)计算ρi

      2) 当ρi < 0时,h=1,计算(ξ±)(1)

      3) h=h+1,重复2)计算(ξ±)(h)、$s_{\hat \xi _i^ \pm }^{\left( h \right)}$;

      4) 当$2\max \left( {s_{\hat \xi _1^ \pm }^{\left( h \right)}, s_{\hat \xi _2^ \pm }^{\left( h \right)} \cdots s_{\hat \xi _m^ \pm }^{\left( h \right)}} \right) < \delta $时,终止循环,据式(22)-(24)计算${\mathit{\boldsymbol{E}}_{{{\hat \xi }^ \pm }}} = \mathit{\boldsymbol{E}}_{{{\hat \xi }^ \pm }}^{\left( h \right)}$、${\mathit{\boldsymbol{b}}_{{{\hat \xi }^ \pm }}}$和${\theta _{{b_{\hat \xi \bar \xi _i^ \pm }}}}$。

      算法1至算法2在进行偏差计算时,如果单位权方差估值偏差过大,则需要利用偏差改正后的单位权方差估值重新构建观测向量协方差阵,迭代计算偏差。

    • 本文把MC、AMC和AAMC输出结果的标准差公式列于表 1中。其中,输出量统一用标量z表示,zσz2表示获得的均值和方差,sz(h)su(z)(h)su(z)(h)对应于AMC中均值的标准差和标准差的标准差,sz±(h)对应于AAMC中均值的标准差,其中su(z)(h)su(z)(h)的差别在于计算方差的自由度分别为M-1和M。令M′=hM可获得由误差传播定律确定的AMC均值与标准差各自的标准差,则AAMC的均值标准差为$\sqrt {\left( {1 + \rho } \right)\sigma _z^2/\left( {hM} \right)} $。

      表 1  MC、AMC、AAMC 3种方法精度对比

      Table 1.  The Comparison for Accuracy of MC, AMC and AAMC

      均值 标准差 MC AMC AAMC
      均值未知 均值的标准差 $\sqrt {\sigma _z^2/M'} $ $s_{\bar z}^{\left( h \right)}$ $s_{{{\bar z}^ \pm }}^{\left( h \right)}$
      标准差的标准差 $\sqrt {\frac{{\sigma _z^2}}{{2\left( {M'-1} \right)}} + \frac{{\sigma _z^2}}{{8{{\left( {M'-1} \right)}^2}}}} $ su(z)(h) -
      均值已知 标准差的标准差 $\sqrt {\frac{{\sigma _z^2}}{{2M'}} + \frac{{\sigma _z^2}}{{8{{\left( {M'} \right)}^2}}}} $ su(z)(h) -

      对于现有的精度评定方法,相比于一阶近似协方差阵,本文AMC既能获得稳定的数值结果,又考虑到非线性平差估值的有偏性。对于参数估值的偏差,AAMC生成随机误差的次数为AMC的一半,且ρi < 0使得AAMC结果精度更高或收敛更快。

    • 本文采用的直线拟合模型定义为:

      $$ {\mathit{\boldsymbol{l}}_y} - {\mathit{\boldsymbol{v}}_y} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{l}}_x}}&{ - {\mathit{\boldsymbol{v}}_x}}&\mathit{\boldsymbol{e}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\xi _1}}\\ {{\xi _2}} \end{array}} \right] $$ (25)

      式中,lyn/2×1的纵坐标观测向量;vyly的误差向量;lxn/2×1的横坐标观测向量,vxlx的误差向量;en/2×1的元素全为1的向量;ξ1ξ2为模型参数。

      由于横坐标也含有观测误差,式(25)可作为EIV或Partial EIV模型求取其总体最小二乘解,但基于GH模型的算法具有一般性,且有利于平差结果的统计分析[32],首先假设模拟的坐标服从均值为

      $$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\bar l}}}_x} = {\left[ {\begin{array}{*{20}{c}} 1&2&3&4&5&6&7&8&9&{10} \end{array}} \right]^{\rm{T}}}\\ {{\mathit{\boldsymbol{\bar l}}}_y} = {\left[ {\begin{array}{*{20}{c}} 3&5&7&9&{11}&{13}&{15}&{17}&{19}&{21} \end{array}} \right]^{\rm{T}}} \end{array} \right. $$ (26)

      且协方差阵为σ02I20的正态分布,σ02取值分别为0.1、0.5和1。

      为了验证自适应蒙特卡罗法的有效性,本文把5×107次的蒙特卡罗模拟结果作为参照标准。观测向量的均值、单位权方差估值和参数估值分别用其真值代替。MC、AMC和AAMC结果和相应精度以及式(4)计算的一阶近似协方差阵(5×107次结果的均值)分别列于表 2表 3中。AMC的δ为0.005,AAMC使用与其相同的模拟次数。

      表 2  直线拟合模型参数偏差计算结果

      Table 2.  The Biases of Parameters Estimates of Straight Line Fitting Model

      参数 项目 ξ1 ξ2
      / 参数真值 2 1
      0.1I20(h=6) 参数估值 2.002 4 0.986 5
      MC 0.002 4(±1×10-5) -0.013 5(±7×10-5)
      AMC 0.002 2±2(3)×10-4 -0.013 8±13(20)×10-4
      AAMC 0.002 5±2(3)×10-5 -0.013 6±1(2)×10-4
      0.5I20(h=52) 参数估值 2.012 5 0.931 2
      MC 0.012 5(±2×10-5) -0.068 8(±2×10-4)
      AMC 0.012 7±2.5(2.5)×10-4 -0.069 1±1.5(1.5)×10-3
      AAMC 0.012 5±4(5)×10-5 -0.068 7±2(3)×10-4
      I20(h=328) 参数估值 2.025 7 0.858 4
      MC 0.025 7(±4×10-5) -0.141 6(±2×10-4)
      AMC 0.025 9±1.5(1.4)×10-4 -0.141 7±9.2(8.9)×10-4
      AAMC 0.025 7±4.4(4.1)×10-5 -0.141 3±2.5(2.4)×10-4

      表 3  直线拟合模型参数标准差及协方差计算结果

      Table 3.  The Standard Deviations and Covariances of Parameters Estimates of Straight Line Fitting Model

      项目 方法 ${\hat \sigma _{{{\hat \xi }_1}}}$ ${\hat \sigma _{{{\hat \xi }_2}}}$ ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$
      0.1I20(h=6) FTT 0.078 2 0.484 9 -0.033 6
      MC 0.078 3(±8×10-6) 0.485 4(±5×10-5) -0.033 7
      AMC 0.078 7±4(23)×10-5 0.487 2±23(14)×10-4 -0.034 0
      0.5I20(h=52) FTT 0.178 0 1.102 0 -0.174 3
      MC 0.179 3(±2×10-5) 1.1074(±1×10-4) -0.176 7
      AMC 0.179 3±6(18)×10-5 1.107 7±25(11)×10-4 -0.1768
      I20(h=328) FTT 0.258 3 1.594 1 -0.367 0
      MC 0.261 6(±3×10-5) 1.608 6(±2×10-4) -0.376 4
      AMC 0.261 6±7(10)×10-5 1.608 2±25(6)×10-4 -0.376 3

      本文采用文献[30]中的直线拟合算例,其相应的坐标观测向量和权阵为

      $$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{l}}_x} = {\left[ {\begin{array}{*{20}{c}} { - 1}&0&1&2&3&4&5 \end{array}} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{l}}_y} = {\left[ {\begin{array}{*{20}{c}} {1.3}&{0.8}&{0.9}&{1.2}&{2.0}&{3.5}&{4.1} \end{array}} \right]^{\rm{T}}} \end{array} \right. $$ (27)
      $$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{P}}_y} = {\rm{diag}}\left( {{{\left[ {\begin{array}{*{20}{c}} 2&8&7&5&{10}&8&6 \end{array}} \right]}^{\rm{T}}}} \right)\\ {\mathit{\boldsymbol{P}}_x} = {\rm{diag}}\left( {{{\left[ {\begin{array}{*{20}{c}} 3&9&8&4&5&7&{10} \end{array}} \right]}^{\rm{T}}}} \right) \end{array} \right. $$ (28)

      式中,diag(·)表示把向量转化成对角矩阵。计算结果分别列于表 4表 5中。AMC偏差计算的δ为0.001,协方差阵计算的δ为0.000 5。MC的参数估值均值和方差(浅灰色线)、AMC每批次参数估值均值和方差(灰色线)和AAMC每批次参数估值的均值(黑色线)的相关关系如图 2所示。

      表 4  直线拟合模型估值偏差计算结果

      Table 4.  The Biases of Estimates of Straight Line Fitting Model

      项目 估值 MC AMC(h=340) AAMC(h=340)
      ξ1 0.658 0 0.005 8(±2×10-6) 0.005 8±6(7)×10-5 0.005 8±2.0(2.5)×10-5
      ξ2 0.551 2 -0.013 1(±5×10-6) -0.013 0±1.8(2.0)×10-4 -0.013 1±5(6)×10-5
      σ02 1.538 6 -0.010 8(±1×10-5) -0.010 5±5.0(0.5)×10-4 -
      bv - 0.004 7 0.004 6 -

      表 5  直线拟合模型参数估值标准差和协方差计算结果

      Table 5.  The Standard Deviations and Covariances of Parameters Estimates of Straight Line Fitting Model

      标准差或协方差 FTT MC AMC(h=64)
      ${\hat \sigma _{{{\hat \xi }_1}}}$ 0.119 5 0.124 9(±1×10-5) 0.124 9±3(11)×10-5
      ${\hat \sigma _{{{\hat \xi }_2}}}$ 0.349 1 0.360 3(±4×10-5) 0.360 1±2(3)×10-4
      ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$ -0.032 2 -0.035 2 -0.035 2

      图  2  直线拟合模型参数估值均值及方差曲线图

      Figure 2.  The Curves of Means and Variances of Parameters Estimates of Straight Line Fitting Model

    • 本文采用的椭圆拟合模型定义为:

      $$ {\left( {\left( {{l_{{x_i}}} - {v_{{x_i}}} - {\xi _1}} \right)/{\xi _3}} \right)^2} + {\left( {\left( {{l_{{y_i}}} - {v_{{y_i}}} - {\xi _2}} \right)/{\xi _4}} \right)^2} = 1 $$ (29)

      式中,i=1, 2…n/2;lyivyilxivxi的定义同式(25);ξ1ξ2ξ3ξ4为模型参数。

      本文采用文献[29]缩小10倍后的等精度观测向量:

      $$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{l}}_x} = {\left[ {\begin{array}{*{20}{c}} 0&5&9&{12}&{13}&{ - 13}&{ - 10}&{ - 5}&0 \end{array}} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{l}}_y} = {\left[ {\begin{array}{*{20}{c}} {12}&{11}&8&0&{ - 5}&{ - 5}&6&{10}&{ - 11} \end{array}} \right]^{\rm{T}}} \end{array} \right. $$ (30)

      并使用与直线拟合模型同样的处理方式,偏差和协方差阵结果分别列于表 6表 7中。AMC的δ为0.001。参数估值均值和方差曲线图如图 3所示。

      表 6  椭圆拟合模型估值偏差计算结果

      Table 6.  The Biases of Estimations of Ellipse Fitting Model

      参数 估值 MC AMC(h=297) AAMC(h=297)
      ξ1 -0.059 8 -0.010 8(±8×10-6) -0.010 9±3.0(3.2)×10-4 -0.010 6±6.3(6.5)×10-5
      ξ2 -0.194 2 0.001 4(±7×10-6) 0.001 3±2.9(3.0)×10-4 0.001 4±4.9(5.1)×10-5
      ξ3 13.108 7 0.092 5(±9×10-6) 0.092 3±3.8(4.0)×10-4 0.092 5±9.5(9.6)×10-5
      ξ4 11.513 1 0.048 8(±9×10-6) 0.049 2±3.6(3.6)×10-4 0.048 8±7.1(6.8)×10-5
      σ02 1.046 4 -0.009 5(±9×10-6) -0.009 5±3.5(0.4)×10-4 -
      bv - 0.022 6 0.0226 -

      表 7  椭圆拟合模型参数估值标准差和协方差计算结果

      Table 7.  The Standard Deviations and Covariances of Parameters Estimations of Ellipse Fitting Model

      参数 FTT MC AMC(h=196)
      ${\hat \sigma _{{{\hat \xi }_1}}}$ 0.534 9 0.546 8(±5×10-5) 0.547 4±3.0(2.8)×10-4
      ${\hat \sigma _{{{\hat \xi }_2}}}$ 0.515 4 0.520 3(±5×10-5) 0.520 4±2.7(2.6)×10-4
      ${\hat \sigma _{{{\hat \xi }_3}}}$ 0.632 7 0.682 3(±7×10-5) 0.682 4±5.0(3.4)×10-4
      ${\hat \sigma _{{{\hat \xi }_4}}}$ 0.606 1 0.624 6(±6×10-5) 0.624 7±3.9(3.2)×10-4
      ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}$ -0.003 3 -0.004 6 -0.004 7
      ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_3}}}$ -0.081 6 -0.095 1 -0.095 0
      ${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_4}}}$ 0.017 2 0.019 1 0.018 6
      ${\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_3}}}$ 0.022 9 0.027 9 0.028 8
      ${\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_4}}}$ -0.148 3 -0.151 2 -0.151 1
      ${\hat \sigma _{{{\hat \xi }_3}{{\hat \xi }_4}}}$ -0.108 7 -0.108 4 -0.109 1

      图  3  椭圆拟合模型参数估值均值及方差曲线图

      Figure 3.  The Curves of Means and Variances of Parameters Estimations of Ellipse Fitting Model

    • 表 7中,${\hat \sigma _{{{\hat \xi }_1}}}、{\hat \sigma _{{{\hat \xi }_2}}}、{\hat \sigma _{{{\hat \xi }_3}}}$和${\hat \sigma _{{{\hat \xi }_4}}}$代表参数估值标准差,${\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_2}}}、{\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_3}}}、{\hat \sigma _{{{\hat \xi }_1}{{\hat \xi }_4}}}、{\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_3}}}、{\hat \sigma _{{{\hat \xi }_2}{{\hat \xi }_4}}}、{\hat \sigma _{{{\hat \xi }_3}{{\hat \xi }_4}}}$代表参数估值间的协方差,$\left\| {{\mathit{\boldsymbol{b}}_{\hat v}}} \right\| = {\left( {b_{{v_{{x_1}}}}^2 + \cdots b_{{v_{{x_{n/2}}}}}^2 + b_{{v_{{y_1}}}}^2 + \cdots + b_{{v_{{y_{n/2}}}}}^2} \right)^{1/2}}$代表改正数偏差二范数。从表 2表 3可以看出,不同σ02下AMC和AAMC都能收敛且获得与MC极为接近的偏差或协方差阵结果,这表明了本文AMC和AAMC的可行性。

      表 2表 4表 6表明3类蒙特卡罗法偏差结果相一致,表 3表 5表 7表明,MC和AMC协方差阵结果基本上小数点后前3位相同,但是在绝对值上都要大于一阶近似协方差阵的结果,这与文献[15]的结果相一致。表 2-7中,3类蒙特卡罗法结果括号内的标准差根据表 1第3列公式确定,其中AMC和AAMC结果括号外的标准差为由式(6)、式(14)和式(21)获得的最终结果。两类自适应蒙特卡罗法结果的两种标准差结果基本一致,表明式(6)、式(14)和式(21)可反映结果的精度。

      图 2图 3表明,随着模拟次数增加,偏差和方差曲线由最初的波动较大逐渐接近并重合于5×107次蒙特卡罗结果对应的曲线。AMC能够在一定数值容差下获得统计意义上稳定的偏差和标准差及其相应的精度,其也可以根据数值容差事先判断结果的大致精度。而预先选定模拟次数的MC对此往往难以实现。

      表 2表 4表 6表明,AAMC的偏差精度较AMC高出一个数量级,这是两类模型参数的相关系数都接近于-1导致的。图 2图 3表明,AAMC的均值变化曲线波动性更小,能够很快收敛。因此,AAMC在偏差计算方面较有优势。

      表 3表 5表 7表明,参数估值的一阶近似协方差阵绝对值结果要小于MC和AMC的结果,其往往会高估非线性平差参数估值的精度。表 2表 4表 6中模型参数估值偏差百分比都有大于或接近于1%的分量,最大为18.1%,最小为0.1%。根据文献[30]给出的1%的经验法则,本文两个模型都具有非线性性态,有必要考虑非线性平差估值的有偏性。

    • 本文在阐述非线性平差原理以及自适应蒙特卡罗原理的基础上,给出了自适应蒙特卡罗法进行非线性平差精度评定的算法流程,提出了快速计算参数估值偏差的对偶自适应蒙特卡罗法。直线拟合模型和椭圆拟合模型的算例结果表明:

      1) 本文的自适应蒙特卡罗法能够顾及估值的有偏性,获得统计意义上稳定及合理的精度评定结果;

      2) 对偶自适应蒙特卡罗法较自适应蒙特卡罗法的偏差计算结果精度更高,收敛更快;当参数估值的相关系数小于零时,对偶自适应蒙特卡罗法可用于快速计算参数估值偏差;

      3) 当参数估值偏差大于或接近于1%时,非线性平差估值的有偏性可以在精度评定中加以考虑。

      有关数值容差更为合理的选取方法以及自适应蒙特卡罗法的进一步应用是接下来要研究的内容。

参考文献 (32)

目录

    /

    返回文章
    返回