2. 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室,南昌市广兰大道418号, 330013
传统的非线性模型协方差传播方法主要是基于泰勒级数展开的近似函数法。但是当解决多维复杂的非线性函数协方差传播问题时,近似函数法获取非线性函数的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,
本文参考文献[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{\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)分别计算均值的方差
$\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) |
式中,y1、y2是经非线性函数转换后的值;x1、x2、x3视为服从正态分布的观测值,其随机模型为:
$\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。
在工程测量中,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) |
式中,a和b分别表示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个测站S和P,其中,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。
由表 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) |
2. Key Laboratory of Mine Environmental Monitoring and Improving around Poyang Lake, MNR, 418 Guanglan Road, Nanchang 330013, China