﻿ 非线性函数协方差传播的Stein Monte Carlo方法
 文章快速检索 高级检索
 大地测量与地球动力学  2022, Vol. 42 Issue (1): 1-4  DOI: 10.14075/j.jgg.2022.01.001

### 引用本文

WANG Leyang, LUO Xinlei. Stein Monte Carlo Method for Nonlinear Function Covariance Propagation[J]. Journal of Geodesy and Geodynamics, 2022, 42(1): 1-4.

### Foundation support

National Natural Science Foundation of China, No. 41874001, 42174011.

### 第一作者简介

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

### 文章历史

1. 东华理工大学测绘工程学院，南昌市广兰大道418号，330013;
2. 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室，南昌市广兰大道418号, 330013

1 Stein Monte Carlo(SMC)方法与非线性函数协方差传播

 ${s^2}=\frac{1}{{{n_1} - 1}}\sum\limits_{i=1}^{{n_1}} {{{({x_i} - \bar x)}^2}}$ (1)

 ${n_2}=\max \left( {\left\lfloor {\frac{{{s^2}({t_{{n_1} - 1}}{{\left( {1 - \alpha /2} \right)}^2})}}{{{\delta ^2}}}} \right\rfloor - {n_1} + 1, 0} \right)$ (2)

1) 确定初始阶段MC模拟的批次数n1(本文根据文献[12]方法，取n1=10)，每批次MC模拟次数为M，数值容差为δ

2) h=1时，根据观测量l的先验均值l和先验方差-协方差矩阵Dy随机生成M组观测值li(i=1, …, M)，将一组观测值代入非线性方程β=φ(l)，得到一组参数估值$\mathit{\boldsymbol{\hat \beta }}$；重复计算M次，然后使用式(3)和式(4)计算M组参数估值的均值${\mathit{\boldsymbol{\bar \beta }}^{(1)}}$和方差u2(${\mathit{\boldsymbol{\hat \beta }}^{(1)}}$)：

 ${{\mathit{\boldsymbol{\bar \beta }}}^{(h)}} = \frac{1}{M}\sum\limits_{i = 1}^M {\mathit{\boldsymbol{\hat \beta }}_{(i)}^{(h)}}$ (3)
 $\begin{array}{*{20}{c}} {{u^2}({{\mathit{\boldsymbol{\hat \beta }}}^{(h)}}) = }\\ {{\rm{diag}}\left( {\frac{1}{{M - 1}}\sum\limits_{i = 1}^M {((\mathit{\boldsymbol{\hat \beta }}_{(i)}^{(h)} - \mathit{\boldsymbol{\bar \beta }}{^{(h)}}){{(\mathit{\boldsymbol{\hat \beta }}_{(i)}^{(h)} - \mathit{\boldsymbol{\bar \beta }}{^{(h)}})}^{\rm{T}}})} } \right)} \end{array}$ (4)

3) 将步骤2)执行n1次，得到n1组均值和方差，通过式(5a)和式(5b)分别计算均值的方差$\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\bar\beta }}^2$和方差的方差$\mathit{\boldsymbol{S}}_{{\mathit{\boldsymbol{u}}^2}({\mathit{\boldsymbol{\hat \beta }}})}^2$

 $\begin{array}{l} \mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\bar\beta }}^2 = {\mathop{\rm diag}\nolimits} \left( {\frac{1}{{M - 1}}\sum\limits_{h = 1}^{{n_1}} {\left( {\left( {{{\mathit{\boldsymbol{\overline \beta }} }^{(h)}} - \frac{1}{{{n_1}}}\sum\limits_{h = 1}^{{n_1}} {{{\mathit{\boldsymbol{\overline \beta }} }^{(h)}}} } \right)} \right.} } \right. \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {{{\left( {{{\mathit{\boldsymbol{\overline \beta }} }^{(h)}} - \frac{1}{{{n_1}}}\sum\limits_{h = 1}^{{n_1}} {{{\mathit{\boldsymbol{\overline \beta }} }^{(h)}}} } \right)}^{\rm{T}}}} \right)} \right) \end{array}$ (5a)
 $\begin{array}{l} \;\;\;\mathit{\boldsymbol{S}}_{{u^2}(\mathit{\boldsymbol{\hat\beta }})}^2 = {\mathop{\rm diag}\nolimits} \left( {\frac{1}{{M - 1}}\sum\limits_{h = 1}^{{n_1}} {\left( {\left( {{u^2}\left( {{{\mathit{\boldsymbol{\widehat \beta }}}^{(h)}}} \right) - } \right.} \right.} } \right.\\ \left. {\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{{n_1}}}\sum\limits_{h = 1}^{{n_1}} {{u^2}} \left( {{{\mathit{\boldsymbol{\widehat \beta }}}^{(h)}}} \right)} \right)\left( {{u^2}\left( {{{\mathit{\boldsymbol{\widehat \beta }}}^{(h)}}} \right) - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left. {{{\left. {\frac{1}{{{n_1}}}\sum\limits_{h = 1}^{{n_1}} {{u^2}} \left( {{{\mathit{\boldsymbol{\widehat \beta }}}^{(h)}}} \right)} \right)}^{\rm{T}}}} \right)} \right) \end{array}$ (5b)

4) 由式(6)确定第二阶段MC模拟的批次数n2:

 $\begin{array}{l} {n_2} = \max \left( {\left\lfloor {\frac{{\mathit{\boldsymbol{S}}_\mathit{\boldsymbol{\bar\beta }}^2\left( {{t_{{n_1} - 1}}{{\left( {1 - \alpha /2} \right)}^2}} \right)}}{{{\delta ^2}}}} \right\rfloor - {n_1} + 1, } \right.\\ \left. {\;\;\;\;\left\lfloor {\frac{{\mathit{\boldsymbol{S}}_{{\mathit{\boldsymbol{u}}^2}\left( {\mathit{\boldsymbol{\hat \beta }}} \right)}^2({t_{{n_1} - 1}}{{\left( {1 - \alpha /2} \right)}^2})}}{{{\delta ^2}}}} \right\rfloor - {n_1} + 1, 0} \right) \end{array}$ (6)

2 算例与分析

2.1 算例1

 $\left\{ \begin{array}{l} {y_1} = 0.6 - 0.28{x_1}{x_2} + 0.25{x_1}{x_3} + 0.36{x_1}{x_2^2} + \\ \;\;\;\;\;0.12{x_2}x_3^2 + 0.49x_2^2 - 0.17x_1^3 + 0.03x_2^3\\ {y_2} = 0.43 + 0.2{x_1}{x_2} - 0.4{x_2}\sqrt {{x_3}} - \\ \;\;\;\;\;\;0.04x_1^2{x_2} + 0.15{x_1}x_2^2 \end{array} \right.$ (7)

 $\left[ {\begin{array}{*{20}{c}} {{x_1}}\\ {{x_2}}\\ {{x_3}} \end{array}} \right]\sim N\left( {\left[ {\begin{array}{*{20}{c}} {1.0}\\ {5.0}\\ {2.6} \end{array}} \right], \left[ {\begin{array}{*{20}{c}} {0.02}&0&0\\ 0&{0.05}&0\\ 0&0&{0.04} \end{array}} \right]} \right)$ (8)

2.2 算例2

 $\left\{ \begin{array}{l} L = {\rm{arcsin}}\left( {\frac{Y}{{\sqrt {{X^2} + {Y^2}} }}} \right)\\ {\rm{tan}}B = \frac{{Z + N{e^2}{\rm{cos}}B}}{{\sqrt {{X^2} + {Y^2}} }} \end{array} \right.$ (9)
 $\left\{ \begin{array}{l} e = \frac{{\sqrt {{a^2} - {b^2}} }}{a}\\ N = \frac{a}{{\sqrt {1 - {e^2}{{\sin }^2}B} }} \end{array} \right.$ (10)

 $\left\{ \begin{array}{l} x = {x_0} + \frac{N}{{2\rho '{'^2}}}\sin B\cos Bl'{'^2} + \frac{N}{{24\rho '{'^4}}}\sin B{\cos ^3}B\cdot\\ \;\;\;\;(5 - {t^2} + 9{\eta ^2} + 4{\eta ^4})l'{'^4} + \frac{N}{{720\rho '{'^6}}}\sin B{\cos ^5}B\cdot\\ \;\;\;\;(61 - 58{t^2} + {t^4})l'{'^6}\\ y = \frac{N}{{\rho ''}}\cos Bl'' + \frac{N}{{6\rho '{'^3}}}{\cos ^3}B(1 - {t^2} + {\eta ^2})l'{'^3} + \\ \;\;\;\;\;\frac{N}{{120\rho '{'^5}}}{\cos ^5}B(5 - 18{t^2} + {t^4} + 14{\eta ^2}-58{\eta ^2}{t^2})l'{'^5} \end{array} \right.$ (11)

 $\left\{ \begin{array}{l} t = \tan B\\ {\eta ^2} = e{'^2}{\cos ^2}B\\ e{'^2} = \frac{{\sqrt {{a^2} - {b^2}} }}{b}\\ {x_0} = {a_0}B - \frac{{{a_2}}}{2}\sin 2B + \frac{{{a_4}}}{4}\sin 4B - \\ \;\;\;\;\frac{{{a_6}}}{6}\sin 6B + \frac{{{a_8}}}{8}\sin 8B \end{array} \right.$ (12)

 $\left\{ \begin{array}{l} {\left[ {X, Y, Z} \right]^{\rm{T}}} = {\left[ { - 1\;500\;000, {\rm{ }}800\;000, 2\;300\;000} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{D}}_{\Delta X\Delta Y\Delta Z}} = \\ \left[ {\begin{array}{*{20}{c}} {7.239\;7 \times {{10}^{ - 4}}}&{7.280\;0 \times {{10}^{ - 6}}}&{7.520\;0 \times {{10}^{ - 6}}}\\ {7.280\;0 \times {{10}^{ - 6}}}&{6.762\;0 \times {{10}^{ - 4}}}&{7.290\;0 \times {{10}^{ - 6}}}\\ {7.520\;0 \times {{10}^{ - 6}}}&{7.290\;0 \times {{10}^{ - 6}}}&{7.310\;0 \times {{10}^{ - 4}}} \end{array}} \right] \end{array} \right.$ (13)

2.3 算例分析

3 结语

 [1] Wang L Y, Zhao Y W. Unscented Transformation with Scaled Symmetric Sampling Strategy for Precision Estimation of Total Least Squares[J]. Studia Geophysica et Geodaetica, 2017, 61(3): 1-27 (0) [2] Wang L Y, Zhao Y W. Scaled Unscented Transformation of Nonlinear Error Propagation: Accuracy, Sensitivity, and Applications[J]. Journal of Surveying Engineering, 2018, 144(1) (0) [3] Wang L Y, Zou C Y. Accuracy Analysis and Applications of the Sterling Interpolation Method for Nonlinear Function Error Propagation[J]. Measurement, 2019, 146: 55-64 DOI:10.1016/j.measurement.2019.06.017 (0) [4] 李朝奎, 黄力民, 曾卓乔, 等. 基于Monte-Carlo方法的强非线性函数方差估计[J]. 中国有色金属学报, 2000, 10(4): 613-617 (Li Chaokui, Huang Limin, Zeng Zhuoqiao, et al. Variance Estimation of Intensive Nonlinear Functions Based on Monte-Carlo Method[J]. The Chinese Journal of Nonferrous Metals, 2000, 10(4): 613-617 DOI:10.3321/j.issn:1004-0609.2000.04.035) (0) [5] Shen Y Z, Li B F, Chen Y. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238 DOI:10.1007/s00190-010-0431-1 (0) [6] Malengo A, Pennecchi F. A Weighted Total Least-Squares Algorithm for any Fitting Model with Correlated Variables[J]. Metrologia, 2013, 50(6): 654-662 DOI:10.1088/0026-1394/50/6/654 (0) [7] 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) (0) [8] Niemeier W, Tengen D. Uncertainty Assessment in Geodetic Network Adjustment by Combining GUM and Monte-Carlo-Simulations[J]. Journal of Applied Geodesy, 2017, 11(2): 67-76 (0) [9] Alkhatib H, Schuh W D. Integration of the Monte Carlo Covariance Estimation Strategy into Tailored Solution Procedures for Large-Scale Least Squares Problems[J]. Journal of Geodesy, 2007, 81(1): 53-66 (0) [10] Stein C. A Two-Sample Test for a Linear Hypothesis whose Power is Independent of the Variance[J]. The Annals of Mathematical Statistics, 1945, 16(3): 243-258 DOI:10.1214/aoms/1177731088 (0) [11] 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 (0) [12] Wübbeler G, Harris P M, Cox M G, et al. A Two-Stage Procedure for Determining the Number of Trials in the Application of a Monte Carlo Method for Uncertainty Evaluation[J]. Metrologia, 2010, 47(3): 317-324 DOI:10.1088/0026-1394/47/3/023 (0) [13] 孔祥元, 郭际明, 刘宗泉. 大地测量学基础[M]. 武汉: 武汉大学出版社, 2010 (Kong Xiangyuan, Guo Jiming, Liu Zongquan. Foundation of Geodesy[M]. Wuhan: Wuhan University Press, 2010) (0) [14] Ghilani C D. Adjustment Computations: Spatial Data Analysis[M]. New Jersey: John Wiley and Sons, 2017 (0)
Stein Monte Carlo Method for Nonlinear Function Covariance Propagation
WANG Leyang1,2     LUO Xinlei1
1. Faculty of Geomatics, East China University of Technology, 418 Guanglan Road, Nanchang 330013, China;
2. Key Laboratory of Mine Environmental Monitoring and Improving around Poyang Lake, MNR, 418 Guanglan Road, Nanchang 330013, China
Abstract: The standard deviations of the observation function can be calculated by the law of covariance propagation. However, the approximate function method based on Taylor series expansion requires complex derivative operations when solving the standard deviation of nonlinear function. The Monte Carlo method can avoid the derivative operation, but it is not objective in the selection of the number of simulations. Moreover, the Monte Carlo method cannot directly control the simulation results. To overcome these disadvantages, we introduce the Stein two-stage method into the covariance propagation theory of nonlinear functions. Combined with the Monte Carlo method, we design the Stein Monte Carlo (SMC) algorithm flow of nonlinear function covariance propagation. We apply the SMC method in the two-dimensional polynomial function and covariance propagation of GNSS baseline vector. Results verify the effectiveness of the SMC method. This method provides a new idea for the covariance propagation of nonlinear models.
Key words: nonlinear function; covariance propagation; Stein two-stage scheme; Monte Carlo method