非線性函數(shù)協(xié)方差傳播的Stein Monte Carlo方法

王樂(lè)洋 羅鑫磊

1 東華理工大學(xué)測(cè)繪工程學(xué)院，南昌市廣蘭大道418號(hào)，330013 2 自然資源部環(huán)鄱陽(yáng)湖區(qū)域礦山環(huán)境監(jiān)測(cè)與治理重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號(hào),330013

傳統(tǒng)的非線性模型協(xié)方差傳播方法主要是基于泰勒級(jí)數(shù)展開的近似函數(shù)法。但是當(dāng)解決多維復(fù)雜的非線性函數(shù)協(xié)方差傳播問(wèn)題時(shí)，近似函數(shù)法獲取非線性函數(shù)的Jacobi矩陣和Hessian矩陣較為困難，其應(yīng)用具有局限性。為避免導(dǎo)數(shù)運(yùn)算，近年來(lái)有學(xué)者提出了近似概率分布法[1-3]，其中Monte Carlo(MC)法是一種較為常見(jiàn)的方法。MC法無(wú)需了解非線性函數(shù)的構(gòu)造，根據(jù)先驗(yàn)信息隨機(jī)生成大量樣本點(diǎn)來(lái)近似非線性函數(shù)輸出量的概率分布，可以避免導(dǎo)數(shù)運(yùn)算，其主要用于求解參數(shù)估值的置信區(qū)間和方差-協(xié)方差矩陣[4-9]。然而，由于無(wú)法確定模擬次數(shù)與模擬結(jié)果精度之間的關(guān)系，現(xiàn)有MC法的模擬次數(shù)通常會(huì)預(yù)先給定，增加了不必要的計(jì)算成本。

為進(jìn)一步發(fā)展非線性函數(shù)協(xié)方差傳播理論方法，保證模擬結(jié)果精度，同時(shí)減少不必要的模擬次數(shù)，本文將Stein兩階段方法[10]與MC法結(jié)合進(jìn)行非線性函數(shù)的協(xié)方差傳播計(jì)算，并通過(guò)2個(gè)算例來(lái)驗(yàn)證本文方法的可行性和有效性。

1 Stein Monte Carlo(SMC)方法與非線性函數(shù)協(xié)方差傳播

兩階段抽樣方法結(jié)合MC法可以用于確定給定模擬結(jié)果精度時(shí)MC法的模擬次數(shù)。下面對(duì)SMC方法進(jìn)行簡(jiǎn)要介紹。

選擇初始MC模擬數(shù)n1，通過(guò)式(1)計(jì)算樣本方差：

(1)

通過(guò)式(2)確定樣本數(shù)量n2，計(jì)算樣本數(shù)量為n1+n2時(shí)的樣本均值：

(2)

本文參考文獻(xiàn)[12]的思想，將上述SMC方法建立在批處理模式的基礎(chǔ)上。首先確定每一批次的MC模擬次數(shù)M(本文中M=1 000)，經(jīng)過(guò)M次MC模擬后可以得到M個(gè)值，將其視為樣本值并計(jì)算這批樣本值的均值和方差，記為第一批次的結(jié)果。將上述過(guò)程執(zhí)行n次，可以得到n個(gè)均值和方差的結(jié)果。將批次n視為樣本的數(shù)量，由中心極限定理可知，當(dāng)批次足夠大時(shí)，均值和方差都可以視為服從正態(tài)分布的變量，所以均值和方差的期望即為待求量。

根據(jù)以上內(nèi)容，將非線性函數(shù)協(xié)方差傳播的SMC方法的算法步驟總結(jié)如下：

1)確定初始階段MC模擬的批次數(shù)n1(本文根據(jù)文獻(xiàn)[12]方法，取n1=10)，每批次MC模擬次數(shù)為M，數(shù)值容差為δ。

(3)

(4)

(5a)

(5b)

4)由式(6)確定第二階段MC模擬的批次數(shù)n2:

(6)

式中，max(·)表示取所有元素的最大值。若n2≤0，那么n1批次所有樣本的均值和方差將作為最終結(jié)果；若n2>0，則繼續(xù)執(zhí)行n2次步驟2)，最終結(jié)果將由n1+n2批次所有樣本計(jì)算得到。

2 算例與分析

本節(jié)共設(shè)計(jì)2個(gè)算例說(shuō)明本文方法的有效性，分別為二維多項(xiàng)式函數(shù)(算例1)和GNSS基線向量協(xié)方差傳播(算例2)。其中，算例1的非線性模型強(qiáng)度弱且觀測(cè)值不相關(guān)，算例2的非線性模型強(qiáng)度強(qiáng)且觀測(cè)值具有相關(guān)性。

2.1 算例1

使用二維多項(xiàng)式函數(shù)模型，定義的二維多項(xiàng)式函數(shù)如下：

(7)

式中，y1、y2是經(jīng)非線性函數(shù)轉(zhuǎn)換后的值；x1、x2、x3視為服從正態(tài)分布的觀測(cè)值，其隨機(jī)模型為：

(8)

為給出精度評(píng)定結(jié)果的參照基準(zhǔn)，2個(gè)算例均將107次MC模擬結(jié)果的近似值視為真值，每批次模擬次數(shù)為1 000，107次MC模擬的批次數(shù)視為10 000。分別將107次MC法均值和δ=0.001時(shí)SMC法的相應(yīng)結(jié)果列于表1，將一階泰勒近似法(FTT)、二階泰勒近似法(STT)、MC法和SMC方法得到的標(biāo)準(zhǔn)差結(jié)果列于表2。

表1 2種方法求得的坐標(biāo)估值結(jié)果

表2 4種方法求得的標(biāo)準(zhǔn)差結(jié)果

2.2 算例2

在工程測(cè)量中，GNSS基線向量網(wǎng)與地面控制網(wǎng)的聯(lián)合平差可以在二維高斯平面坐標(biāo)系中完成，因此，需要將GNSS基線向量的隨機(jī)模型傳遞到高斯平面坐標(biāo)系中。隨機(jī)模型的傳遞過(guò)程如下：首先由GNSS網(wǎng)中的固定點(diǎn)坐標(biāo)和兩點(diǎn)之間的基線向量得到某觀測(cè)點(diǎn)的空間直角坐標(biāo)；然后通過(guò)空間坐標(biāo)和大地坐標(biāo)之間的轉(zhuǎn)換及高斯投影正算得到該點(diǎn)的高斯平面坐標(biāo)，并通過(guò)協(xié)方差傳播律得到GNSS平面基線向量的隨機(jī)模型。因此，若想得到平面基線向量的方差-協(xié)方差矩陣，首先需要獲得高斯平面坐標(biāo)的方差-協(xié)方差矩陣。

空間直角坐標(biāo)(X,Y,Z)與大地坐標(biāo)(B,L)之間具有如下關(guān)系[13]：

(9)

(10)

式中，a和b分別表示W(wǎng)GS 84橢球的長(zhǎng)半軸和短半軸。

大地坐標(biāo)(B,L)與高斯平面坐標(biāo)(x,y)之間具有如下非線性關(guān)系[13]：

(11)

式中，l表示經(jīng)差；其中相關(guān)參數(shù)的計(jì)算公式如下[13]：

(12)

設(shè)某GNSS網(wǎng)中的2個(gè)測(cè)站S和P，其中，P測(cè)站坐標(biāo)和PS之間的基線向量的方差-協(xié)方差矩陣信息(數(shù)據(jù)源于文獻(xiàn)[14])如下：

(13)

由式(9)和式(11)可知，空間直角坐標(biāo)與高斯平面坐標(biāo)之間具有復(fù)雜的非線性關(guān)系，若通過(guò)近似函數(shù)法求解平面基線向量的方差-協(xié)方差矩陣，計(jì)算過(guò)程復(fù)雜且難以實(shí)現(xiàn)，因此本算例僅對(duì)比SMC法與MC法的結(jié)果以驗(yàn)證本文方法的有效性。表3給出了2種方法求得的S點(diǎn)坐標(biāo)估值及其標(biāo)準(zhǔn)差，其中，SMC法的δ=0.000 01。

表3 2種方法求得的坐標(biāo)估值及標(biāo)準(zhǔn)差結(jié)果

2.3 算例分析

由表1和表3可見(jiàn)，SMC法和MC法坐標(biāo)估值的計(jì)算結(jié)果基本一致。然而，本文方法的重點(diǎn)在于獲取待定參數(shù)的后驗(yàn)精度信息，下面對(duì)其作進(jìn)一步分析。

從計(jì)算結(jié)果的有效性來(lái)看，對(duì)于非線性強(qiáng)度較弱的二維多項(xiàng)式函數(shù)模型，表2給出了4種方法求得的標(biāo)準(zhǔn)差結(jié)果，對(duì)比4種方法發(fā)現(xiàn)，由于只考慮了一階精度，F(xiàn)TT法的精度最差；若顧及二階精度，STT法精度較FTT法有所提升。由于MC法顧及了高階項(xiàng)的精度信息，因此得到的精度結(jié)果最優(yōu)。SMC法和MC法的結(jié)果僅在小數(shù)點(diǎn)后第4位有微小差別，因此可以認(rèn)為2種方法的計(jì)算結(jié)果基本一致，驗(yàn)證了SMC方法的可行性。對(duì)于非線性強(qiáng)度較強(qiáng)、函數(shù)形式復(fù)雜且觀測(cè)值具有相關(guān)性的GNSS基線向量模型，表3給出了2種方法求得的標(biāo)準(zhǔn)差結(jié)果，可以看到，SMC法與MC法得到的標(biāo)準(zhǔn)差結(jié)果一致，同樣可以驗(yàn)證SMC方法的有效性。由于SMC法包含了高階項(xiàng)的精度信息，因此在實(shí)際應(yīng)用中，通過(guò)SMC法可以將GNSS基線向量網(wǎng)的隨機(jī)模型更精確地傳遞到高斯平面坐標(biāo)系中，從而得到更為準(zhǔn)確的聯(lián)合平差結(jié)果。

從計(jì)算效率來(lái)看，MC法的模擬次數(shù)是預(yù)先指定的，無(wú)法確定達(dá)到指定計(jì)算精度時(shí)的模擬次數(shù)，這無(wú)疑會(huì)增加不必要的計(jì)算成本。SMC方法則可以很好地解決這一問(wèn)題。從上述2個(gè)算例的計(jì)算結(jié)果可以看到，通過(guò)預(yù)先指定數(shù)值容差，SMC方法能夠在保證計(jì)算結(jié)果正確的前提下降低計(jì)算成本，避免了MC法模擬次數(shù)盲目選擇的缺陷，大大提高了計(jì)算效率。

3 結(jié) 語(yǔ)

為了進(jìn)一步豐富和發(fā)展非線性函數(shù)協(xié)方差傳播理論，本文提出SMC方法，并給出SMC法用于非線性函數(shù)協(xié)方差傳播的算法步驟。二維多項(xiàng)式函數(shù)和GNSS基線向量協(xié)方差傳播的算例結(jié)果表明，將SMC法用于非線性函數(shù)協(xié)方差傳播是有效的，其可以同時(shí)顧及計(jì)算效率和計(jì)算精度，獲得合理的精度結(jié)果。此外，SMC法原理簡(jiǎn)單，易于程序?qū)崿F(xiàn)，適用性強(qiáng)，在處理實(shí)際多維復(fù)雜的非線性函數(shù)協(xié)方差傳播問(wèn)題時(shí)，具有一定的應(yīng)用價(jià)值。