一種基于廣義延拓逼近的動(dòng)態(tài)濾波方法初步研究*

劉 成，李 芳，高為廣，王 威

(1.北京跟蹤與通信技術(shù)研究所，北京 100094；2.中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100101)

1 概 述

廣義延拓逼近方法自1986年首次提出以來(lái)[1]展示了強(qiáng)大的生命力，應(yīng)用橫跨天文、力學(xué)、電磁學(xué)、機(jī)械、電子、通信和導(dǎo)航等多個(gè)領(lǐng)域，不斷證明了其科學(xué)性和有效性。起初，廣義延拓逼近方法的提出是為了解決大型天文望遠(yuǎn)鏡設(shè)計(jì)中的實(shí)際問(wèn)題，其后很快在一些光學(xué)與射電望遠(yuǎn)鏡、雷達(dá)天線、彈性力學(xué)以及電磁場(chǎng)理論等領(lǐng)域得到實(shí)踐，并被用于構(gòu)造有限單元內(nèi)的插值函數(shù)[2]。1991年，文[3]首次將廣義延拓逼近方法引入太陽(yáng)物理學(xué)，證明了延拓邊界元空間的逼近性質(zhì)和邊界元近似解的收斂性，并于1993年出版[4]。2004年，文[5]重新對(duì)方法進(jìn)行了系統(tǒng)歸納和整理。隨后，廣義延拓逼近方法開始受到科研工作者更加廣泛的關(guān)注，并被推廣至電子[6]、電力[7-8]、機(jī)械[9-12]、特種設(shè)備[13]、氣象觀測(cè)[14]等工程領(lǐng)域以及導(dǎo)彈制導(dǎo)等軍事領(lǐng)域[15-16]。

同時(shí)，廣義延拓逼近方法開始在導(dǎo)航與定位領(lǐng)域廣泛應(yīng)用。文[17]將其用于全球定位系統(tǒng)信號(hào)遮擋情況下，對(duì)偽距測(cè)量值進(jìn)行外推，以維持用戶基本定位能力，結(jié)果表明，可在15 s內(nèi)維持優(yōu)于5 m的偽距外推精度。文[18]對(duì)接收機(jī)鐘漂進(jìn)行建模以增強(qiáng)定位，結(jié)果表明，除維持遮擋情況下的定位外，還可提高正常信號(hào)情況下約24%的定位精度。文[19-24]先后開展了廣義延拓逼近方法在導(dǎo)航衛(wèi)星廣播星歷和精密星歷插值中的應(yīng)用，均表明精度優(yōu)于一般插值方法。文[24]的研究表明，在選取合適參數(shù)時(shí)，廣義延拓可達(dá)到毫米級(jí)的軌道插值精度，較拉格朗日(Lagrange)插值法高出56.6%～83.5%，較切比雪夫插值法高出20.9%～73.5%，且避免了傳統(tǒng)插值方法的高階振蕩現(xiàn)象。文[25-27]研究了廣義延拓逼近方法在衛(wèi)星精密鐘差中的插值應(yīng)用，獲得的絕大多數(shù)數(shù)據(jù)優(yōu)于0.03 ns、全部數(shù)據(jù)優(yōu)于0.1 ns的插值精度。文[28-29]利用廣義延拓進(jìn)行高程擬合，結(jié)果表明，精度優(yōu)于傳統(tǒng)二次多項(xiàng)式曲面法和多面函數(shù)法。文[30-33]將其用于導(dǎo)航信號(hào)的頻率及載波跟蹤，證明了廣義延拓方法具有更優(yōu)的信號(hào)跟蹤性能。文[34-35]將其用于對(duì)載波測(cè)量值周跳的檢測(cè)和修復(fù)，不同的研究結(jié)果均表明，廣義延拓方法能夠準(zhǔn)確探測(cè)和修復(fù)1周以上的周跳。文[36-37]將其用于接收機(jī)碼鑒相器中的碼相位測(cè)量，結(jié)果表明廣義延拓方法具有更高的精度和穩(wěn)定性。文[38-41]探討了廣義延拓方法在衛(wèi)星導(dǎo)航以及室內(nèi)導(dǎo)航算法中的優(yōu)化應(yīng)用，進(jìn)一步拓展了其應(yīng)用場(chǎng)景。

上述諸多研究和應(yīng)用均獲得了理想效果，證明了廣義延拓方法的科學(xué)性和正確性。然而，無(wú)論是軌道、鐘差、高程還是偽距、載波、碼相關(guān)峰，其特點(diǎn)都是一維、靜態(tài)數(shù)據(jù)下的插值或外推計(jì)算，局限于廣義延拓方法的一般形式。廣義延拓方法能利用插值多項(xiàng)式實(shí)現(xiàn)對(duì)狀態(tài)量的外推預(yù)測(cè)，然而，將其與濾波方法相結(jié)合，應(yīng)用于動(dòng)態(tài)數(shù)據(jù)處理領(lǐng)域，需要解決其測(cè)量修正問(wèn)題并構(gòu)建相應(yīng)的誤差模型，目前鮮有相關(guān)研究性工作。本文在上述研究工作的基礎(chǔ)上對(duì)傳統(tǒng)廣義延拓方法進(jìn)行了拓展，提出將其推廣至動(dòng)態(tài)濾波數(shù)據(jù)處理領(lǐng)域，研究和構(gòu)造了基于廣義延拓插值方程的濾波模型與框架，并推導(dǎo)給出了求解方法與流程。所提出的廣義延拓逼近濾波方法既具備多維狀態(tài)量的數(shù)據(jù)處理能力，也能滿足狀態(tài)量動(dòng)態(tài)更新的要求，豐富了廣義延拓方法的內(nèi)涵。分別利用仿真數(shù)據(jù)和北斗衛(wèi)星導(dǎo)航系統(tǒng)(BeiDouSatellite Navigation System, BDS)實(shí)測(cè)數(shù)據(jù)對(duì)方法進(jìn)行了初步驗(yàn)證，結(jié)果表明，在高斯噪聲、非高斯噪聲以及實(shí)際導(dǎo)航應(yīng)用環(huán)境下均獲得了比卡爾曼濾波更好的性能，具有較大的發(fā)展?jié)摿蛻?yīng)用價(jià)值。

2 廣義延拓逼近數(shù)學(xué)方法

廣義延拓逼近方法的數(shù)學(xué)本質(zhì)是基于 “編結(jié)” 思想將數(shù)據(jù)集劃分為多個(gè)相互連結(jié)的單元延拓域，并分段進(jìn)行兼顧插值與擬合的數(shù)值處理，如圖1，記某一單元延拓域Δe內(nèi)的廣義延拓二次插值函數(shù)U(t)為

圖1 延拓域及其對(duì)應(yīng)函數(shù)

U(t)=a0+a1t+a2t2,t∈[t1,t2]，

(1)

其中，ti為數(shù)據(jù)對(duì)應(yīng)的時(shí)刻；a0,a1,a2為廣義延拓多項(xiàng)式系數(shù)，可由延拓域中的數(shù)學(xué)模型確定：

(2)

也就是說(shuō)，一方面要求廣義延拓插值多項(xiàng)式在延拓域邊界點(diǎn)上滿足插值條件，保持邊界點(diǎn)的數(shù)據(jù)特征和約束作用，并使得各段數(shù)據(jù)之間的變化具有更好的協(xié)調(diào)性。另一方面，利用分段插值區(qū)域的周圍結(jié)點(diǎn)(包括內(nèi)點(diǎn))實(shí)現(xiàn)分段區(qū)域內(nèi)外部數(shù)據(jù)的殘差最小擬合，充分結(jié)合了插值和擬合方法的優(yōu)點(diǎn)。

廣義延拓逼近方法能根據(jù)實(shí)際情況，靈活構(gòu)建具體的數(shù)學(xué)模型。例如，當(dāng)對(duì)某個(gè)邊界域節(jié)點(diǎn)的數(shù)值精度具有較高置信度時(shí)，可以將其增加作為插值約束點(diǎn)，以獲得更準(zhǔn)確的擬合結(jié)果；反之，可以放棄某個(gè)插值節(jié)點(diǎn)，將其退化為擬合點(diǎn)，降低對(duì)整體結(jié)果的影響，也可以采用三次或更高次的廣義延拓多項(xiàng)式，以獲得更高的數(shù)值處理精度。

3 基于廣義延拓逼近的動(dòng)態(tài)濾波方法

3.1 基本原理

基于廣義延拓逼近的動(dòng)態(tài)濾波方法的基本思想是利用廣義延拓法構(gòu)造一個(gè)具有一定先驗(yàn)數(shù)據(jù)長(zhǎng)度、可靈活設(shè)置和優(yōu)選插值約束點(diǎn)的移動(dòng)窗口，對(duì)系統(tǒng)狀態(tài)量進(jìn)行高精度、具有抗差能力的逼近和預(yù)測(cè)，以取代傳統(tǒng)卡爾曼濾波方法中的狀態(tài)量時(shí)間更新(狀態(tài)轉(zhuǎn)移)過(guò)程，如圖2。

與一般卡爾曼濾波方法相比，廣義延拓逼近濾波方法使用了數(shù)據(jù)窗口，能夠起到類似事后平滑的效果?？柭鼮V波的數(shù)學(xué)本質(zhì)是基于馬爾科夫鏈的一步狀態(tài)轉(zhuǎn)移，利用最新歷元n的狀態(tài)量預(yù)測(cè)下一歷元n+ 1的狀態(tài)量。它的優(yōu)點(diǎn)在于大幅減小了計(jì)算量，但同時(shí)也導(dǎo)致過(guò)于依賴當(dāng)前歷元狀態(tài)量的精度；換言之，即使根據(jù)誤差協(xié)方差矩陣P判斷出當(dāng)前歷元n的狀態(tài)量最優(yōu)估計(jì)值具有很大誤差，也只能繼續(xù)利用該狀態(tài)量估計(jì)值進(jìn)行下一步預(yù)測(cè)。相比之下，廣義延拓逼近濾波利用一個(gè)數(shù)據(jù)窗口進(jìn)行預(yù)測(cè)，具有更好的平滑性和穩(wěn)定性。

與一般的開窗濾波方法相比，廣義延拓逼近濾波方法可以在數(shù)據(jù)窗口內(nèi)靈活選擇插值約束點(diǎn)，從而實(shí)現(xiàn)更好的抗差能力。例如，當(dāng)當(dāng)前歷元n的定位誤差較大時(shí)，可將其退化作為擬合點(diǎn)，以弱化其作用，并從窗口信息中另行選擇更加精確的狀態(tài)量估計(jì)值作為插值約束點(diǎn)(如圖2)，從而在符合載體實(shí)際運(yùn)動(dòng)軌跡的同時(shí)，充分發(fā)揮可靠先驗(yàn)信息的約束作用，規(guī)避粗差影響，提高狀態(tài)轉(zhuǎn)移的精度和可靠性。

圖2 廣義延拓動(dòng)態(tài)濾波原理示意圖

3.2 插值約束方程的構(gòu)造與狀態(tài)量的時(shí)間轉(zhuǎn)移

3.2.1 基本形式

若已知一段窗口長(zhǎng)度為L(zhǎng)的先驗(yàn)狀態(tài)量最優(yōu)估計(jì)值x^i(i=n-L+ 1,n-L+ 2, …,n)，通過(guò)將窗口范圍內(nèi)的某一歷元j(j∈[n-L+ 1,n])的狀態(tài)量最優(yōu)估計(jì)值x^j作為插值點(diǎn)，其他歷元的狀態(tài)量最優(yōu)估計(jì)值作為擬合點(diǎn)，可建立一個(gè)具有插值約束條件的廣義延拓逼近方程：

(3)

(4)

(5)

(6)

3.2.2 擴(kuò)展形式

上節(jié)給出了廣義延拓逼近方程的一般形式。實(shí)際上，方程也可采用更高階的插值多項(xiàng)式形式：

(7)

或者可以進(jìn)一步增加和采用某另一歷元k(k∈[n-L+1,n]且k≠j)的狀態(tài)量最優(yōu)估計(jì)值，從而構(gòu)成具有多個(gè)插值約束條件的方程形式：

(8)

以上擴(kuò)展形式可根據(jù)具體需求靈活構(gòu)造。

3.3 完整動(dòng)態(tài)濾波框架

利用廣義延拓逼近方法完成狀態(tài)量的插值平滑更新后，結(jié)合測(cè)量更新處理步驟，即可建立完整的優(yōu)化動(dòng)態(tài)濾波框架模型，如圖3。

圖3 廣義延拓動(dòng)態(tài)濾波模型

(9)

進(jìn)行更新，其中，Q為系統(tǒng)轉(zhuǎn)移過(guò)程噪聲；A為系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣。當(dāng)廣義延拓插值方程采用二次多項(xiàng)式時(shí)，對(duì)應(yīng)狀態(tài)量是勻速變化模型，有

(10)

其中，Δt為兩時(shí)刻之間的時(shí)間間隔。當(dāng)廣義延拓插值方程采用三次多項(xiàng)式時(shí)，對(duì)應(yīng)狀態(tài)量是勻加速變化模型，有

(11)

在卡爾曼濾波測(cè)量更新過(guò)程中，分別計(jì)算n+ 1歷元的權(quán)增益矩陣Kn+1、狀態(tài)量最優(yōu)估計(jì)值X^n+1及其協(xié)方差矩陣Pn+1：

(12)

其中，C為測(cè)量關(guān)系矩陣；R為測(cè)量噪聲協(xié)方差矩陣。

4 求解方法與步驟

對(duì)廣義延拓動(dòng)態(tài)濾波逼近方程的求解，可采用線性化展開的形式。具體地，將(3)式中的插值約束方程改寫為

(13)

將(13)式代入最小二乘公式可得

(14)

(15)

(16)

(17)

(18)

(19)

5 仿真計(jì)算

5.1 信號(hào)去噪處理仿真

分別對(duì)正弦信號(hào)在白噪聲和有色噪聲情況下的濾波性能進(jìn)行仿真。仿真時(shí)間T=0～30 s，步長(zhǎng)Δt=0.1 s，重復(fù)計(jì)算1 000次。利用蒙特卡洛模擬對(duì)精度結(jié)果進(jìn)行統(tǒng)計(jì)，通過(guò)與廣泛應(yīng)用的卡爾曼濾波器進(jìn)行比較，對(duì)廣義延拓逼近濾波的性能進(jìn)行分析和評(píng)估。

5.1.1 白噪聲濾波性能

首先，加入滿足r～N(0, 0.1)分布的高斯白噪聲進(jìn)行數(shù)值仿真，考察和分析濾波方法的基本特性。兩種濾波器的初始狀態(tài)量為X0=[0 0]，初始協(xié)方差為P0=diag([0.1 2])，初始過(guò)程噪聲Q0=[0.0 0.000 5; 0.000 5 0.01]，初始測(cè)量噪聲為R0=0.01，廣義延拓逼近濾波的窗口長(zhǎng)度L=20，插值約束點(diǎn)為當(dāng)前歷元觀測(cè)點(diǎn)。兩種濾波方法的效果如圖4，誤差概率統(tǒng)計(jì)如圖5。其中，卡爾曼濾波的均方根誤差約0.083，廣義延拓逼近濾波的均方根誤差約0.069，提升約17%。

圖4 白噪聲下的濾波仿真結(jié)果

圖5 白噪聲下的濾波誤差概率統(tǒng)計(jì)

5.1.2 有色噪聲濾波性能

在上述高斯白噪聲基礎(chǔ)上，加入有色噪聲e(k)=0.2(x(k)+0.9x(k-1))，濾波器誤差參數(shù)和窗口數(shù)據(jù)長(zhǎng)度均保持不變。廣義延拓逼近方法和卡爾曼濾波方法均基于噪聲高斯分布的理想假設(shè)，但在實(shí)際應(yīng)用中，理想高斯分布的情況很難真實(shí)存在，因此，考察兩種方法在有色噪聲影響下的性能具有實(shí)際意義。兩種濾波方法的效果如圖6，誤差概率統(tǒng)計(jì)如圖7。其中，卡爾曼濾波的均方根誤差約0.113，廣義延拓逼近濾波的均方根誤差約0.102，提升約10%。

圖6 疊加有色噪聲后的濾波仿真結(jié)果

圖7 疊加有色噪聲后的濾波誤差概率統(tǒng)計(jì)

5.2 衛(wèi)星導(dǎo)航實(shí)測(cè)數(shù)據(jù)處理驗(yàn)證

2020年1月，采用自主研制的全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System, GNSS)多模導(dǎo)航接收機(jī)在北京地區(qū)測(cè)量得到的北斗衛(wèi)星導(dǎo)航系統(tǒng)觀測(cè)數(shù)據(jù)(包括偽距、多普勒、載波和廣播星歷)對(duì)方法進(jìn)行實(shí)測(cè)驗(yàn)證，并與卡爾曼濾波器以及經(jīng)典希格(Sage)開窗濾波器處理結(jié)果進(jìn)行比較。幾種方法的定位輸出結(jié)果如圖8，其中，廣義延拓逼近濾波的窗口長(zhǎng)度取L=15。

從圖8可以看出，廣義延拓濾波方法有能力獲得比經(jīng)典卡爾曼濾波和希格開窗濾波更優(yōu)的性能，具有應(yīng)用潛力。以實(shí)時(shí)動(dòng)態(tài)(Real-time Kinemaic, RTK)載波定位解算結(jié)果作為基準(zhǔn)，統(tǒng)計(jì)3種方法的東西向平面誤差，結(jié)果見表1。

表1 濾波定位誤差結(jié)果對(duì)比

圖8 北斗導(dǎo)航濾波定位結(jié)果對(duì)比

6 結(jié) 論

本文在廣義延拓逼近理論的基礎(chǔ)上，提出了一種動(dòng)態(tài)濾波更新方法。研究和構(gòu)造了基于廣義延拓插值方程的濾波模型與框架，并推導(dǎo)給出了其求解方法與流程。方法既具備多維狀態(tài)量的處理能力，也能夠滿足高低階狀態(tài)量動(dòng)態(tài)更新的要求。在一段移動(dòng)平滑窗口內(nèi)構(gòu)造和求解廣義延拓插值多項(xiàng)式，實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)量的時(shí)間轉(zhuǎn)移，具有更好的平滑效果，拓展了經(jīng)典廣義延拓?cái)?shù)學(xué)方法的內(nèi)涵。與一般開窗濾波方法相比，廣義延拓逼近濾波方法數(shù)據(jù)窗口內(nèi)插值點(diǎn)的位置及數(shù)量均可靈活選取，從而可通過(guò)鎖定最新觀測(cè)點(diǎn)或最準(zhǔn)確的先驗(yàn)觀測(cè)點(diǎn)，實(shí)現(xiàn)更高的精度和抗差能力。分別利用正弦信號(hào)數(shù)值仿真和北斗衛(wèi)星導(dǎo)航實(shí)測(cè)數(shù)據(jù)對(duì)方法進(jìn)行了驗(yàn)證，結(jié)果表明，在高斯噪聲和非高斯噪聲仿真下均獲得了比卡爾曼濾波更優(yōu)的性能，驗(yàn)證了將廣義延拓逼近理論與濾波方法相結(jié)合的可行性。另外需要指出的是，目前初步構(gòu)造廣義延拓濾波模型仍是基于卡爾曼濾波誤差假設(shè)條件的近似，相信隨著后續(xù)理論和模型的完善，方法將具備更優(yōu)的性能與更大的應(yīng)用價(jià)值。