文章快速检索     高级检索
  大地测量与地球动力学  2022, Vol. 42 Issue (1): 1-4  DOI: 10.14075/j.jgg.2022.01.001

引用本文  

王乐洋, 罗鑫磊. 非线性函数协方差传播的Stein Monte Carlo方法[J]. 大地测量与地球动力学, 2022, 42(1): 1-4.
WANG Leyang, LUO Xinlei. Stein Monte Carlo Method for Nonlinear Function Covariance Propagation[J]. Journal of Geodesy and Geodynamics, 2022, 42(1): 1-4.

项目来源

国家自然科学基金(41874001, 42174011)。

Foundation support

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

第一作者简介

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

About the first author

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

文章历史

收稿日期:2021-04-16
非线性函数协方差传播的Stein Monte Carlo方法
王乐洋1,2     罗鑫磊1     
1. 东华理工大学测绘工程学院,南昌市广兰大道418号,330013;
2. 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室,南昌市广兰大道418号, 330013
摘要:基于泰勒级数展开的近似函数法在求解非线性函数的中误差时需要进行复杂的导数计算,已有的Monte Carlo法虽然可以避免导数运算,但在模拟次数的选择上不具有客观性,且无法直接控制模拟结果。因此,将Stein两阶段法融入非线性函数的协方差传播理论中,并与Monte Carlo方法结合,设计了一套非线性函数协方差传播的Stein Monte Carlo算法流程。将该方法用于二维多项式函数和GNSS基线向量的协方差传播计算中,实验结果验证了其有效性,为非线性模型协方差传播的计算提供了一种新思路。
关键词非线性函数协方差传播Stein两阶段法Monte Carlo方法

传统的非线性模型协方差传播方法主要是基于泰勒级数展开的近似函数法。但是当解决多维复杂的非线性函数协方差传播问题时,近似函数法获取非线性函数的Jacobi矩阵和Hessian矩阵较为困难,其应用具有局限性。为避免导数运算,近年来有学者提出了近似概率分布法[1-3],其中Monte Carlo(MC)法是一种较为常见的方法。MC法无需了解非线性函数的构造,根据先验信息随机生成大量样本点来近似非线性函数输出量的概率分布,可以避免导数运算,其主要用于求解参数估值的置信区间和方差-协方差矩阵[4-9]。然而,由于无法确定模拟次数与模拟结果精度之间的关系,现有MC法的模拟次数通常会预先给定,增加了不必要的计算成本。

为进一步发展非线性函数协方差传播理论方法,保证模拟结果精度,同时减少不必要的模拟次数,本文将Stein两阶段方法[10]与MC法结合进行非线性函数的协方差传播计算,并通过2个算例来验证本文方法的可行性和有效性。

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

两阶段抽样方法结合MC法可以用于确定给定模拟结果精度时MC法的模拟次数。下面对SMC方法进行简要介绍。

选择初始MC模拟数n1,通过式(1)计算样本方差:

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

通过式(2)确定样本数量n2,计算样本数量为n1+n2时的样本均值:

${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)

式中,tn1-1(1-α/2)表示当自由度为n1-1时,对于给定的显著性水平αt分布的上分位点值。本文中设定α=0.05,$ \left\lfloor p \right\rfloor $为不大于p的最大整数,δ为给定的数值容差。当n2>0时,继续进行MC模拟,则最终的估计值为x=x(n1+n2)。目前,数值容差δ的确定方法有2种:1)数值容差的一般形式为δ=0.5×10l,其中l一般为非正整数,具体取值可根据数值的有效数字确定[11],本文取2位有效数字。为了得到准确的结果,本文中$ \delta =\frac{1}{5}{\rm{ }} \times 0.5 \times {10^l}$。2)Stein证明了该估计值位于[x-δ, x+δ] 区间内的概率为1-α,因此可以根据需要直接指定δ的值。

本文参考文献[12]的思想,将上述SMC方法建立在批处理模式的基础上。首先确定每一批次的MC模拟次数M(本文中M=1 000),经过M次MC模拟后可以得到M个值,将其视为样本值并计算这批样本值的均值和方差,记为第一批次的结果。将上述过程执行n次,可以得到n个均值和方差的结果。将批次n视为样本的数量,由中心极限定理可知,当批次足够大时,均值和方差都可以视为服从正态分布的变量,所以均值和方差的期望即为待求量。

根据以上内容,将非线性函数协方差传播的SMC方法的算法步骤总结如下:

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)

式中,$ {{\mathit{\boldsymbol{\bar \beta }}}^{\left( h \right)}}$为第h批次中M组参数估值的均值,$ \mathit{\boldsymbol{\hat \beta }}_{\left( i \right)}^{\left( h \right)}$为第h批次中第i个参数估值,u2($ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( h \right)}}$) 为第h批次中M组参数估值的方差;diag(·)为将矩阵的对角线元素组成向量。

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)

式中,max(·)表示取所有元素的最大值。若n2≤0,那么n1批次所有样本的均值和方差将作为最终结果;若n2>0,则继续执行n2次步骤2),最终结果将由n1+n2批次所有样本计算得到。

2 算例与分析

本节共设计2个算例说明本文方法的有效性,分别为二维多项式函数(算例1)和GNSS基线向量协方差传播(算例2)。其中,算例1的非线性模型强度弱且观测值不相关,算例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)

式中,y1y2是经非线性函数转换后的值;x1x2x3视为服从正态分布的观测值,其随机模型为:

$\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个算例均将107次MC模拟结果的近似值视为真值,每批次模拟次数为1 000,107次MC模拟的批次数视为10 000。分别将107次MC法均值和δ=0.001时SMC法的相应结果列于表 1,将一阶泰勒近似法(FTT)、二阶泰勒近似法(STT)、MC法和SMC方法得到的标准差结果列于表 2

表 1 2种方法求得的坐标估值结果 Tab. 1 The results of coordinate estimation calculated by two methods

表 2 4种方法求得的标准差结果 Tab. 2 The results of standard deviation calculated by four methods
2.2 算例2

在工程测量中,GNSS基线向量网与地面控制网的联合平差可以在二维高斯平面坐标系中完成,因此,需要将GNSS基线向量的随机模型传递到高斯平面坐标系中。随机模型的传递过程如下:首先由GNSS网中的固定点坐标和两点之间的基线向量得到某观测点的空间直角坐标;然后通过空间坐标和大地坐标之间的转换及高斯投影正算得到该点的高斯平面坐标,并通过协方差传播律得到GNSS平面基线向量的随机模型。因此,若想得到平面基线向量的方差-协方差矩阵,首先需要获得高斯平面坐标的方差-协方差矩阵。

空间直角坐标(X, Y, Z)与大地坐标(B, L)之间具有如下关系[13]

$\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)

式中,ab分别表示WGS 84椭球的长半轴和短半轴。

大地坐标(B, L)与高斯平面坐标(x, y) 之间具有如下非线性关系[13]

$\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)

式中,l表示经差;其中相关参数的计算公式如下[13]

$\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)

设某GNSS网中的2个测站SP,其中,P测站坐标和PS之间的基线向量的方差-协方差矩阵信息(数据源于文献[14])如下:

$\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)

由式(9)和式(11)可知,空间直角坐标与高斯平面坐标之间具有复杂的非线性关系,若通过近似函数法求解平面基线向量的方差-协方差矩阵,计算过程复杂且难以实现,因此本算例仅对比SMC法与MC法的结果以验证本文方法的有效性。表 3给出了2种方法求得的S点坐标估值及其标准差,其中,SMC法的δ=0.000 01。

表 3 2种方法求得的坐标估值及标准差结果 Tab. 3 The results of coordinate estimation and standard deviation obtained by two methods
2.3 算例分析

表 1表 3可见,SMC法和MC法坐标估值的计算结果基本一致。然而,本文方法的重点在于获取待定参数的后验精度信息,下面对其作进一步分析。

从计算结果的有效性来看,对于非线性强度较弱的二维多项式函数模型,表 2给出了4种方法求得的标准差结果,对比4种方法发现,由于只考虑了一阶精度,FTT法的精度最差;若顾及二阶精度,STT法精度较FTT法有所提升。由于MC法顾及了高阶项的精度信息,因此得到的精度结果最优。SMC法和MC法的结果仅在小数点后第4位有微小差别,因此可以认为2种方法的计算结果基本一致,验证了SMC方法的可行性。对于非线性强度较强、函数形式复杂且观测值具有相关性的GNSS基线向量模型,表 3给出了2种方法求得的标准差结果,可以看到,SMC法与MC法得到的标准差结果一致,同样可以验证SMC方法的有效性。由于SMC法包含了高阶项的精度信息,因此在实际应用中,通过SMC法可以将GNSS基线向量网的随机模型更精确地传递到高斯平面坐标系中,从而得到更为准确的联合平差结果。

从计算效率来看,MC法的模拟次数是预先指定的,无法确定达到指定计算精度时的模拟次数,这无疑会增加不必要的计算成本。SMC方法则可以很好地解决这一问题。从上述2个算例的计算结果可以看到,通过预先指定数值容差,SMC方法能够在保证计算结果正确的前提下降低计算成本,避免了MC法模拟次数盲目选择的缺陷,大大提高了计算效率。

3 结语

为了进一步丰富和发展非线性函数协方差传播理论,本文提出SMC方法,并给出SMC法用于非线性函数协方差传播的算法步骤。二维多项式函数和GNSS基线向量协方差传播的算例结果表明,将SMC法用于非线性函数协方差传播是有效的,其可以同时顾及计算效率和计算精度,获得合理的精度结果。此外,SMC法原理简单,易于程序实现,适用性强,在处理实际多维复杂的非线性函数协方差传播问题时,具有一定的应用价值。

参考文献
[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