梁新茹，高長(zhǎng)生，荊武興

(哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱 150001)

0 引 言

近些年來,高超聲速飛行器因其自身具有機(jī)動(dòng)能力強(qiáng)、機(jī)動(dòng)范圍廣、能夠長(zhǎng)時(shí)飛行等特點(diǎn),引起了越來越廣泛的關(guān)注[1-3]。與傳統(tǒng)飛行器相比,高超聲速飛行器具有飛行空域復(fù)雜、飛行范圍廣等特性[4-5],導(dǎo)致在高超聲速飛行器跟蹤過程中探測(cè)噪聲獲取的統(tǒng)計(jì)信息不準(zhǔn)確,同時(shí)探測(cè)噪聲容易產(chǎn)生異常值,使得實(shí)際測(cè)量噪聲呈現(xiàn)出非平穩(wěn)非高斯特性。

針對(duì)實(shí)際噪聲呈現(xiàn)非平穩(wěn)非高斯這一問題,文獻(xiàn)[6]指出閃爍噪聲作為非高斯噪聲中的一種特殊分布,通常是由異常值引起的,其分布具有重尾特性[7]。但是在傳統(tǒng)的濾波框架內(nèi),通常假設(shè)噪聲分布為高斯分布[8-9]。針對(duì)探測(cè)過程中出現(xiàn)的異常值問題,先驗(yàn)的高斯噪聲建模將會(huì)損失建模的精度,使得建模成本和建模誤差大大增加。針對(duì)重尾分布建模問題,目前主要采用學(xué)生t (Student’s t, ST)分布及其魯棒算法進(jìn)行描述。文獻(xiàn)[10]指出多元ST分布能夠?yàn)榫哂挟惓Ｖ档臄?shù)據(jù)提供更為完善的描述,ST分布可以把高斯分布當(dāng)作一種特殊情況包含在內(nèi)。文獻(xiàn)[11]提供一種基于ST分布近似噪聲分布的后驗(yàn)概率密度函數(shù)(Probability density function, PDF)的濾波算法。文獻(xiàn)[12]設(shè)計(jì)了一種ST分布與粒子濾波相結(jié)合的策略。但是對(duì)于非平穩(wěn)重尾噪聲分布來說,單一的ST分布難以準(zhǔn)確地描述。文獻(xiàn)[6,13]提出一種混合概率可自適應(yīng)變化的高斯-學(xué)生t混合(Gaussian-student’s t mixture, GSTM)分布以適應(yīng)噪聲的非平穩(wěn)變化。但是上述文獻(xiàn)均沒有考慮噪聲先驗(yàn)信息未知的情況,且需要噪聲的先驗(yàn)信息,才能保證算法的成立。

高超聲速飛行器飛行空域廣,空域環(huán)境復(fù)雜,導(dǎo)致噪聲的先驗(yàn)統(tǒng)計(jì)信息復(fù)雜,噪聲可能是時(shí)變、多峰的,無法給出準(zhǔn)確的先驗(yàn)統(tǒng)計(jì)信息。針對(duì)這一問題,文獻(xiàn)[14-16]提出利用逆Wishart分布對(duì)噪聲的協(xié)方差進(jìn)行建模,以解決噪聲統(tǒng)計(jì)特性不精確的問題。文獻(xiàn)[17]在協(xié)方差建模的基礎(chǔ)上,進(jìn)一步提出將噪聲的均值和協(xié)方差聯(lián)合分布近似為獨(dú)立的高斯-逆Wishart分布。文獻(xiàn)[18]將逆Wishart分布和ST建模結(jié)合,文獻(xiàn)[19]將高斯-逆Wishart分布應(yīng)用在ST建模中,完成對(duì)噪聲的建模。同時(shí),為了進(jìn)一步提高濾波算法對(duì)于時(shí)變且先驗(yàn)信息未知噪聲的魯棒性,需要解決ST分布中自由度采用固定常值這一問題,文獻(xiàn)[20]不再采用固定的自由度參數(shù),而是用伽馬函數(shù)來刻畫自由度參數(shù),提高了算法對(duì)于測(cè)量異常值的魯棒性。

本文針對(duì)復(fù)雜空域飛行的高超聲速飛行器探測(cè)過程出現(xiàn)的噪聲異常值干擾,以及異常值的隨機(jī)不確定性導(dǎo)致量測(cè)噪聲的非平穩(wěn)非高斯與統(tǒng)計(jì)信息未知一系列問題,在GSTM分布的基礎(chǔ)上提出了一種改進(jìn)魯棒高斯-學(xué)生t混合分布濾波(Robust Gaussian-student’s t mixture distribution filtering, RGSTMF)算法。首先,引入了GSTM分布,使得算法能夠自適應(yīng)地調(diào)整模型概率來適應(yīng)高超聲速飛行器探測(cè)噪聲的非平穩(wěn)特性。其次,結(jié)合高斯-逆Wishart分布,利用伽馬函數(shù)對(duì)GSTM分布進(jìn)行改進(jìn),改進(jìn)后的GSTM分布可以利用高斯-逆Wishart分布代替GSTM中的零均值假設(shè)來描述實(shí)際應(yīng)用中的均值未知且時(shí)變的噪聲。伽馬函數(shù)代替了原有GSTM分布中固定的自由度以解決噪聲非高斯程度不確定的問題。因此,改進(jìn)后的GSTM分布可以實(shí)現(xiàn)在測(cè)量噪聲先驗(yàn)統(tǒng)計(jì)信息未知的情況下,對(duì)噪聲進(jìn)行準(zhǔn)確建模。最后,將改進(jìn)后的GSTM分布與變分貝葉斯算法相結(jié)合,推導(dǎo)出RGSTMF濾波算法中參數(shù)后驗(yàn)更新的顯式表達(dá)式,從而構(gòu)成完整的濾波算法。

1 高超聲速飛行器助推滑翔段建模

本文主要對(duì)高超聲速飛行器助推滑翔段的跟蹤問題進(jìn)行研究,以HTV-2飛行器作為典型代表。參考文獻(xiàn)[21-22],在彈道坐標(biāo)系下對(duì)高超聲速飛行器滑翔段建模,得到公式如下:

(1)

(2)

式中:關(guān)于升力系數(shù)CL和阻力系數(shù)CD的計(jì)算,本文引用文獻(xiàn)[23]中的數(shù)據(jù),將CL和CD擬合成馬赫數(shù)Ma和攻角αT的函數(shù),如下:

CL=-0.177 3+0.046 4αT+0.085 6e0.032 1Ma+

(3)

(4)

本文的大氣密度采用指數(shù)模型

ρT(h)=ρ0e-h/H

(5)

式中:ρ0取值為ρ0=1.46×10-3kg/m3;H=6 970 m。

地球自轉(zhuǎn)相關(guān)的高階項(xiàng)如下:

(6)

式中:ωe表示地球自轉(zhuǎn)角速率的大小。

2 RGSTMF算法

2.1 GSTM分布

在高超聲速飛行器實(shí)際跟蹤過程中,由于目標(biāo)特性和飛行環(huán)境的復(fù)雜性,跟蹤量測(cè)值極易出現(xiàn)異常值,這將導(dǎo)致噪聲分布的非高斯閃爍特性。因此引入GSTM分布對(duì)噪聲進(jìn)行建模,GSTM分布公式如下:

pN(x)=N(x;μ,Σ)

(7)

pST(x)=ST(x;μ,Σ,v)

(8)

pGSTM(x)=δpN(x)+(1-δ)pST(x)

(9)

式中:N(·)表示高斯分布;μ表示分布的均值;Σ表示分布的協(xié)方差;ST(·)表示ST分布;v表示ST分布的自由度;δ表示混合概率。

參考文獻(xiàn)[6],混合概率δ滿足

p(δ)=B(δ;e,1-e)

(10)

式中:B(·)表示Beta分布;e表示Beta的先驗(yàn)形狀參數(shù)。

從式(9)中可以看出當(dāng)混合概率δ=1時(shí),GSTM分布服從高斯分布,當(dāng)混合概率δ=0時(shí),GSTM分布表現(xiàn)為ST分布,因此高斯分布和ST分布是GSTM分布的一種特殊形式。GSTM分布在應(yīng)對(duì)非平穩(wěn)噪聲方面具有一定的優(yōu)勢(shì),具有更加廣泛的應(yīng)用。

引入隱變量ξ,式(9)可以被改寫為

pGSTM(x|ξ)=δN(x;μ,Σ)+(1-δ)N(x;μ,Σ/ξ)

(11)

式中:隱變量ξ滿足

(12)

因?yàn)樵肼暯y(tǒng)計(jì)特性未知且時(shí)變,固定自由度參數(shù)v的使用會(huì)導(dǎo)致噪聲建模的不準(zhǔn)確性,因此,將自由度參數(shù)v建模成為伽馬函數(shù)[18]。

引用伯努利隨機(jī)變量y,式(11)可以被改寫為

pGSTM(x|ξ,y)=[N(x;μ,Σ)]y+

[N(x;μ,Σ/ξ)](1-y)

(13)

式中:伯努利隨機(jī)變量y滿足

p(y|δ)=δy(1-δ)y

(14)

因此,式(9)可以寫作

p(ξ|v)p(δ)p(v)dξdys.t.y∈{0, 1}

(15)

式中:關(guān)于隱變量ξ的積分為定積分,關(guān)于伯努利隨機(jī)變量y的積分為不定積分,使得最終的積分結(jié)果為y的函數(shù),用于積分外側(cè)的求和運(yùn)算。

2.2 問題描述

構(gòu)建高超聲速飛行器非線性跟蹤系統(tǒng):

xk=f(xk-1)+ωk-1

(16)

zk=h(xk)+υk

(17)

式中:xk∈Rn表示k時(shí)刻的狀態(tài)向量;f(·)表示狀態(tài)函數(shù);h(·)表示量測(cè)函數(shù);zk∈Rm表示k時(shí)刻的量測(cè)向量;wk-1表示系統(tǒng)噪聲,滿足高斯分布;υk表示量測(cè)噪聲。

對(duì)于高超聲速飛行器來說,由于飛行器飛行空域復(fù)雜,目標(biāo)特性復(fù)雜,在跟蹤過程中,容易出現(xiàn)異常值,這會(huì)導(dǎo)致目標(biāo)的量測(cè)噪聲出現(xiàn)非高斯特性。同時(shí)由于異常值的非平穩(wěn)特性,導(dǎo)致目標(biāo)量測(cè)噪聲存在統(tǒng)計(jì)特性未知的情況。因此,在實(shí)際應(yīng)用過程中,量測(cè)噪聲υk的表達(dá)形式和統(tǒng)計(jì)信息沒有先驗(yàn)的假設(shè),量測(cè)噪聲υk可能是非高斯分布,甚至可能是多峰的、傾斜的分布。本文的目的是在沒有先驗(yàn)統(tǒng)計(jì)特性信息的條件下,對(duì)量測(cè)噪聲進(jìn)行建模,在此基礎(chǔ)上推導(dǎo)完成濾波算法。

2.3 先驗(yàn)分布的選擇

根據(jù)貝葉斯公式,為了求解濾波的后驗(yàn)狀態(tài)估計(jì),需要對(duì)一步預(yù)測(cè)和似然函數(shù)進(jìn)行建模。

由于過程噪聲的高斯特性,將一步預(yù)測(cè)建模為高斯分布:

p(xk|z1:k-1)=N(xk;xk/k-1,Pk/k-1)

(18)

式中:xk/k-1表示k時(shí)刻狀態(tài)估計(jì)的一步預(yù)測(cè)均值;Pk/k-1表示k時(shí)刻狀態(tài)估計(jì)的一步預(yù)測(cè)協(xié)方差,同時(shí),下文中出現(xiàn)的下角標(biāo)(·)k/k-1均表示該參數(shù)(·)的一步預(yù)測(cè)值,下角標(biāo)(·)1:k-1表示該參數(shù)(·)的1～k-1時(shí)刻的值。

由于測(cè)量過程中異常值的影響,量測(cè)噪聲呈現(xiàn)出非高斯特性,因此利用GSTM分布來刻畫似然函數(shù)。同時(shí)考慮到在目標(biāo)測(cè)量時(shí),飛行器的飛行空域和自身輻射特性的復(fù)雜性會(huì)導(dǎo)致量測(cè)信息統(tǒng)計(jì)特性未知,因此將量測(cè)噪聲建模成均值未知且時(shí)變的GSTM分布,根據(jù)量測(cè)噪聲得到似然函數(shù)建模公式如下:

[N(xk;h(xk)+Ak,Βk/ξk)](1-yk)·

p(δk)p(vk)dξkdyks.t.yk∈{0,1}

(19)

式中:與式(15)相同,隱變量ξk的積分為定積分,關(guān)于伯努利隨機(jī)變量yk的積分為不定積分,使得最終的積分結(jié)果為yk的函數(shù),用于積分外側(cè)的求和運(yùn)算。Ak和Βk分別表示量測(cè)噪聲均值和協(xié)方差,參考文獻(xiàn)[17],將均值和協(xié)方差的聯(lián)合函數(shù)建模成為高斯-逆Wishart分布,公式如下:

p(Ak|Bk)=N(Ak;ak/k-1,bk/k-1Bk)

(20)

p(Bk|z1:k-1)=IW(Bk;uk/k-1,Uk/k-1)

(21)

式中:IW(·)表示逆Wishart分布,公式如下

IW(B;λ,Ψ)=

(22)

根據(jù)2.1小節(jié)中的建模思路,將ST分布中的自由度參數(shù)vk表示為伽馬分布。

p(vk)=G(vk;ck/k-1,dk/k-1)

(23)

混合概率δk滿足p(δk)=B(δk;e0,1-e0),其中,e0∈[0,1]表示PDFp(δk)的先驗(yàn)形狀參數(shù)。

根據(jù)上文的建模,得到相對(duì)應(yīng)的概率圖模型,如圖1所示。

圖1 概率圖模型Fig.1 Graphical model

在Kalman濾波框架內(nèi),已知一步預(yù)測(cè)的建模,利用容積濾波器(Cubature Kalman filter, CKF)濾波方法得到時(shí)間更新方程:

(24)

(25)

式中:Qk-1表示過程噪聲協(xié)方差;2m表示粒子數(shù);ζi和wi的計(jì)算公式如下

Sk-1/k-1=chol(Pk-1/k-1)

(26)

(27)

(28)

式中:chol(·)表示Cholesky分解,[1]i表示點(diǎn)集[1]的第i組列向量,[1]=[Im,-Im],Im為m維的單位矩陣。

在飛行器濾波跟蹤算法的實(shí)際應(yīng)用過程中,p(Bk|Bk-1),p(Ak|Ak-1),p(vk|vk-1)無法獲取,本文參考文獻(xiàn)[24],引入遺忘因子ρ來描述系統(tǒng)的動(dòng)態(tài)不確定性,ρ的取值范圍為[0.9,1)。

uk/k-1=ρ(uk-1/k-1-md-1)+md+1

(29)

Uk/k-1=ρUk-1/k-1

(30)

ak/k-1=ak-1/k-1,bk/k-1=ρbk-1/k-1

(31)

ck/k-1=ρck-1/k-1,dk/k-1=ρdk-1/k-1

(32)

式中:md表示矩陣維度。

2.4 后驗(yàn)分布的變分逼近

根據(jù)2.3小節(jié)中對(duì)于一步預(yù)測(cè)和似然函數(shù)的建模,為了得到濾波算法的量測(cè)更新,需要對(duì)聯(lián)合后驗(yàn)PDFp(Θk|z1:k)進(jìn)行推導(dǎo),其中Θk包含一步預(yù)測(cè)和似然函數(shù)中的未知參數(shù),定義Θk為變分變量,Θk={xk,Ak,Βk,yk,δk,ξk,vk}。

2.4.1變分貝葉斯

由于聯(lián)合后驗(yàn)PDF中各參數(shù)的耦合,對(duì)于后驗(yàn)參數(shù)的估計(jì)難以得到解析解,因此,引入變分貝葉斯(Variational Bayesian, VB)算法[25],利用自由因子分布近似聯(lián)合分布,聯(lián)合后驗(yàn)PDFp(xk,Ak,Βk,yk,δk,ξk,vk|z1:k)可以寫作:

p(Θk|z1:k)≈q(xk)q(Ak)q(Βk)q(yk)q(δk)·

q(ξk)q(vk)

(33)

式中:q(·)表示p(·)的近似后驗(yàn)概率分布。

為了得到后驗(yàn)PDF,需要最小化q(xk),q(Ak),q(Βk),q(yk),q(δk),q(ξk),q(vk)和聯(lián)合后驗(yàn)PDFp(Θk|z1:k)之間的Kullback-Leibler散度(Kullback-Leibler divergence, KLD)[26]。

{q(xk),q(Ak),q(Βk),q(yk),q(δk),q(ξk),q(vk)}=

argmin KLD(q(xk),q(Ak),q(Βk),q(yk),

q(δk),q(ξk),q(vk)||p(Θk|z1:k))

(34)

式中:KLD用來衡量?jī)蓚€(gè)概率分布之間的差異程度[27],計(jì)算公式為:

(35)

根據(jù)式(35)可以看出,當(dāng)兩個(gè)概率分布之間的差異性越小,KLD越小,因此,可以通過最小化KLD,求取概率分布p(x)的近似分布q(x)。

根據(jù)VB計(jì)算方法,式(34)的最優(yōu)解滿足下面的公式[24]:

lnq(θ)=EΘ(-θ)[lnp(Θk,z1:k)]+cθ

(36)

式中:θ∈Θ并且滿足Θ(-θ)∪θ=Θ;cθ表示獨(dú)立于θ的常值。

由于Θ中的各個(gè)參數(shù)之間相互耦合,無法直接得到其解析解,因此需要定點(diǎn)迭代求解式(36)。利用q(j)(Θ(-θ))的結(jié)果,近似求解后驗(yàn)PDFq(θ)在j+1次迭代時(shí)的PDFq(j+1)(θ),直到收斂到式(36)的局部最優(yōu)值[24]。下面在2.4.2小節(jié)中推導(dǎo)各個(gè)參數(shù)的迭代更新過程。

2.4.2變分迭代過程

已知聯(lián)合分布p(Θk,z1:k)可以被分解為

p(Θk,z1:k)=p(z1:k-1)N(xk;xk/k-1,Pk/k-1)·

[N(xk;h(xk)+Ak,Βk)]yk[N(xk;h(xk)+

Ak,Βk/ξk)](1-yk)N(Ak;ak/k-1,bk/k-1Bk)·

IW(Bk;uk/k-1,Uk/k-1)δkyk(1-δk)ykB(δk;e0,

(37)

1) 令θ={Ak,Bk}

在第j+1次循環(huán)時(shí),{Ak,Bk}的聯(lián)合后驗(yàn)PDF更新為高斯-逆Wishart分布:

(38)

式中:下角標(biāo)(·)k/k均表示該參數(shù)(·)的估計(jì)值。

(39)

(40)

式中:E(j)[x]表示參數(shù)x在j次迭代后的數(shù)學(xué)期望。

(41)

(42)

(43)

(44)

式中:

(45)

Sk/k-1=chol(Pk/k-1)

(46)

(47)

(48)

2) 令θ={yk}

在第j+1次循環(huán)時(shí),yk的后驗(yàn)PDF更新為Bernoulli分布:

(49)

(50)

(51)

根據(jù)式(49)～(51)和Bernoulli分布的性質(zhì),可以求解得到:

(52)

3) 令θ={δk}

在第j+1次循環(huán)時(shí),δk的后驗(yàn)PDF更新為Beta分布

(53)

式中:

(54)

(55)

根據(jù)式(53)～(55)和Beta分布的性質(zhì),可以求解得到:

(56)

(57)

4) 令θ=(ξk)

(58)

式中:

(59)

(60)

根據(jù)式(58)～(60)和伽馬分布的性質(zhì),可以求解得到:

(61)

(62)

5) 令θ={vk}

在第j+1次循環(huán)時(shí),vk后驗(yàn)PDF更新為

(63)

式中:

(64)

(65)

6) 令θ={xk}

在第j+1次循環(huán)時(shí),xk的后驗(yàn)PDF更新為

(66)

式中:

(67)

(68)

(69)

(70)

式中:

(71)

(72)

3 仿真校驗(yàn)

3.1 高超聲速飛行器HTV-2的建模仿真

根據(jù)小節(jié)1中對(duì)高超聲速飛行器HTV-2助推滑翔段的動(dòng)力學(xué)描述,以及考慮飛行器在實(shí)際飛行過程中的飛行走廊約束,完成對(duì)HTV-2的彈道規(guī)劃[28-29]。

表1 HTV-2高超聲速飛行器彈道仿真參數(shù)Table 1 Trajectory simulation parameters of hypersonic vehicle HTV-2

仿真結(jié)果圖彈道如圖2所示:

圖2 彈道形態(tài)Fig.2 Ballistic configuration

3.2 仿真場(chǎng)景設(shè)置

采用地基雷達(dá)觀測(cè),觀測(cè)方程如下:

(73)

式中:Aak,Eek,Rrk分別為雷達(dá)探測(cè)的高低角方位角和斜距;υk為量測(cè)噪聲。

(74)

式中:w.p.表示“with probability”;p1,p2,p3表示異常值概率。R表達(dá)式如下:

(75)

式中:σA=σE=0.03°,σR=20 m。

(76)

根據(jù)第1小節(jié)對(duì)高超聲速飛行器的動(dòng)力學(xué)建模分析,在高超聲速飛行器的滑翔段,飛行器主要受重力和氣動(dòng)力的控制,因此可以從動(dòng)力學(xué)的角度出發(fā),設(shè)計(jì)合理的模型參數(shù)對(duì)動(dòng)力學(xué)進(jìn)行建模,利用擴(kuò)維后的動(dòng)力學(xué)模型進(jìn)行狀態(tài)估計(jì)[8,28,30-31]。

本文參考文獻(xiàn)[28],采用彈道系數(shù)模型作為跟蹤模型。定義彈道系數(shù)為

(77)

(78)

為了比較不同算法對(duì)目標(biāo)的狀態(tài)估計(jì)精度,給出狀態(tài)的均方根誤差(RMSE)和平均均方根誤差(ARMSE)的定義如下:

RMSEk,X=

(79)

ARMSEX=

(80)

為了對(duì)比算法的估計(jì)性能,本文對(duì)比了5種濾波算法:傳統(tǒng)的CKF算法,量測(cè)噪聲采用高斯建模的自適應(yīng)變分貝葉斯容積卡爾曼濾波(Adaptive variational Bayesian cubature Kalman filtering, AVBCKF-R)算法[32],采用ST建模的魯棒變分貝葉斯學(xué)生t濾波(Robust variational Bayesian Student’s t filt-ering, RVBST)算法[33],未包含測(cè)量統(tǒng)計(jì)信息建模的GSTMF算法[6],以及本文中提出的改進(jìn)RGSTMF算法。同時(shí),為了更好地比較GSTM分布的建模優(yōu)勢(shì),將AVBCKF-R算法和RVBST算法中的噪聲均值按照本文提出的方法進(jìn)行建模。

3.2.1仿真場(chǎng)景1

在仿真場(chǎng)景1中,設(shè)置噪聲異常值變化幅度較小且變化平緩,設(shè)置異常值概率分別為:p1=0.01,p2=0.02,p3=0.04。在場(chǎng)景1中設(shè)置異常值所占的概率較小,噪聲重尾分布較輕。仿真結(jié)果見圖3～圖4和表2。

表2 不同濾波算法的ARMSEsTable 2 The ARMSEs of different filtering algorithms

圖3 位置變量的RMSEFig.3 The RMSEs of position

圖4 速度變量的RMSEFig.4 The RMSEs of velocity

從仿真圖和表格可以看出,由于缺少對(duì)噪聲的建模,CKF算法精度表現(xiàn)最差,本文提出的RGSTMF算法精度最優(yōu)。同時(shí)可以看出,與噪聲高斯建模的AVBCKF-R算法和ST建模的RVBST算法相比,本文提出的RGSTMF算法精度更高,這說明GSTM分布與高斯分布和ST分布相比,能更好地?cái)M合非平穩(wěn)噪聲,這與文獻(xiàn)[5]中的結(jié)論一致。此外,與文獻(xiàn)[5]中的GSTM算法相比,RGSTMF算法包含對(duì)噪聲統(tǒng)計(jì)特性的描述,因此在面對(duì)統(tǒng)計(jì)信息未知的量測(cè)噪聲時(shí),具有更好的魯棒性和更高的計(jì)算精度。

3.2.2仿真場(chǎng)景2

在仿真場(chǎng)景2中,設(shè)置噪聲變化幅度較大,設(shè)置異常值概率分別為:p1=0.1,p2=0.2,p3=0.4。在場(chǎng)景2中設(shè)置異常值所占的概率較大,噪聲重尾程度較重。仿真結(jié)果見圖5～圖6和表3。

表3 不同濾波算法的ARMSETable 3 The ARMSEs of different filtering algorithms

圖5 位置變量的RMSEFig.5 The RMSEs of position

圖6 速度變量的RMSEFig.6 The RMSEs of velocity

從仿真圖可以看出,在噪聲異常值進(jìn)一步增加的情況下,CKF算法已經(jīng)出現(xiàn)了較大的誤差,GSTM分布的誤差也進(jìn)一步加大,而文中所提出的RGSTMF算法仍然能保持較高的精度。

綜合仿真1和仿真2,可以看出對(duì)比現(xiàn)有的濾波算法,本文所提出的RGSTMF算法濾波精度最高。與噪聲高斯建模的AVBCKF-R算法和ST建模的RVBST算法相比,所提出的RGSTMF算法無論是在面對(duì)異常值概率較低、變化較為緩慢的噪聲分布,還是異常值占比較高、噪聲突變的噪聲分布,都能夠表現(xiàn)出較為良好的跟蹤精度,這說明GSTM分布對(duì)非平穩(wěn)噪聲具有更好的擬合精度。同時(shí),可以看出,采用對(duì)噪聲均值建模的方法可以進(jìn)一步提升RGSTMF算法的魯棒性,使得其與CKF和GSTMF相比,能在噪聲統(tǒng)計(jì)信息未知的情況下保持較好的跟蹤精度。

3.3 遺忘因子的影響

本小節(jié)探討影響算法的可調(diào)節(jié)參數(shù)遺忘因子ρ對(duì)算法精度的影響,仿真結(jié)果見圖7和圖8。

圖7 場(chǎng)景1中不同的遺忘因子對(duì)RMSE的影響Fig.7 RMSE performances of different forgetting factors in Scenario 1

圖8 場(chǎng)景2中不同的遺忘因子對(duì)RMSE的影響Fig.8 RMSE performances of different forgetting factors in Scenario 2

從仿真圖可以看出,當(dāng)遺忘因子選取為ρ=0.9,0.92,0.94,0.96,0.98時(shí),算法的精度沒有出現(xiàn)較大的差異,但是當(dāng)算法采用ρ=1時(shí),濾波誤差開始增大,并且隨著噪聲異常值概率的增加,誤差越明顯。這是因?yàn)楫?dāng)ρ=1,算法p(Bk|Bk-1),p(Ak|Ak-1),p(vk|vk-1)中的參數(shù)采用固定的建模參數(shù),難以應(yīng)對(duì)時(shí)變的噪聲。

4 結(jié) 論

針對(duì)高超聲速飛行器飛行過程中出現(xiàn)的統(tǒng)計(jì)特性未知的非平穩(wěn)非高斯探測(cè)噪聲,本文提出一種改進(jìn)RGSTMF算法。利用高斯-逆Wishart分布構(gòu)建時(shí)變且未知的均值模型用于代替原始GSTM中的零均值假設(shè),同時(shí)引入伽馬函數(shù)分布來刻畫不確定的非高斯程度,改進(jìn)的RGSTM分布使得濾波算法對(duì)非平穩(wěn)噪聲的建模精度得到進(jìn)一步的提高;將變分貝葉斯算法與RGSTM分布相結(jié)合,解析推導(dǎo)出了后驗(yàn)概率更新的顯式表達(dá)式。仿真結(jié)果驗(yàn)證了本文提出方案的可行性,與其他先進(jìn)算法相比,本文中所提出的RGSTMF算法能夠在先驗(yàn)信息未知的情況下,更好地刻畫噪聲時(shí)變且未知的均值和非高斯程度,在面對(duì)統(tǒng)計(jì)信息未知的非平穩(wěn)非高斯程度不確定的噪聲時(shí),表現(xiàn)出良好的魯棒性和跟蹤精度,并且隨著非高斯程度的進(jìn)一步增加,算法仍能保持良好的跟蹤精度。