非线性函数协方差传播的Stein Monte Carlo方法

王乐洋 罗鑫磊

1 东华理工大学测绘工程学院，南昌市广兰大道418号，330013 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)计算样本方差：

(1)

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

(2)

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

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

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

(3)

(4)

(5a)

(5b)

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

(6)

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

2 算例与分析

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

2.1 算例1

使用二维多项式函数模型，定义的二维多项式函数如下：

(7)

式中，y1、y2是经非线性函数转换后的值；x1、x2、x3视为服从正态分布的观测值，其随机模型为：

(8)

为给出精度评定结果的参照基准，2个算例均将107次MC模拟结果的近似值视为真值，每批次模拟次数为1 000，107次MC模拟的批次数视为10 000。分别将107次MC法均值和δ=0.001时SMC法的相应结果列于表1，将一阶泰勒近似法(FTT)、二阶泰勒近似法(STT)、MC法和SMC方法得到的标准差结果列于表2。

表1 2种方法求得的坐标估值结果

表2 4种方法求得的标准差结果

2.2 算例2

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

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

(9)

(10)

式中，a和b分别表示WGS 84椭球的长半轴和短半轴。

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

(11)

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

(12)

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

(13)

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

表3 2种方法求得的坐标估值及标准差结果

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法原理简单，易于程序实现，适用性强，在处理实际多维复杂的非线性函数协方差传播问题时，具有一定的应用价值。