陳林秀, 郝明瑞, 趙佳佳
(復(fù)雜系統(tǒng)控制與智能協(xié)同技術(shù)重點實驗室, 北京 100074)
基于只測角體制的飛行器以其導(dǎo)引頭作用距離遠、不主動發(fā)射電磁波、隱蔽性好等優(yōu)點,已成為打擊雷達等輻射源目標的主戰(zhàn)武器。為應(yīng)對目標雷達突然關(guān)機,飛行器制導(dǎo)系統(tǒng)設(shè)計的一個關(guān)鍵問題就是如何通過導(dǎo)引頭探測到的信息對目標進行被動定位,即基于角度觀測信號的被動定位濾波算法設(shè)計。
現(xiàn)有的目標定位濾波算法大多采用卡爾曼濾波(KF)算法,基本KF算法是一種線性最小方差估計。但是在被動定位系統(tǒng)中,不論在直角坐標系中還是極坐標系中,目標定位系統(tǒng)都是一個強非線性系統(tǒng)。近年來出現(xiàn)了眾多關(guān)于非線性濾波的理論和方法,其中Arasaratnam等[1]于2009年提出的容積卡爾曼濾波(CKF)算法在濾波精度、算法穩(wěn)定性和計算量等方面表現(xiàn)較為突出。研究成果表明該算法在處理高維濾波問題時,可獲得比擴展卡爾曼濾波(EKF)、中心差分卡爾曼濾波(CDKF)、無跡卡爾曼濾波(UKF)等非線性估計算法更高的濾波精度、穩(wěn)定性和相對低的計算復(fù)雜度[2-5]。然而,CKF算法和基本線性濾波方法一樣,均假定過程噪聲及觀測噪聲為高斯白噪聲,實際上在被動定位應(yīng)用中,過程噪聲來源于運動學建模,其建模準確度較高;而觀測噪聲受電磁波傳播空間、目標所處環(huán)境、電磁干擾以及天線罩等因素的影響,往往是含有時變有色噪聲、跳變噪聲等多種特性的混合噪聲,在此情況下,若仍將其當作白噪聲對待則可能會導(dǎo)致濾波精度下降甚至濾波器發(fā)散[6-10]。文獻[11]推導(dǎo)了有色噪聲條件下的EKF算法,該算法結(jié)構(gòu)簡單,但存在復(fù)雜的雅克比矩陣計算問題,且濾波精度不高,受參數(shù)影響較大。文獻[12]利用歷元噪聲的相關(guān)性特征構(gòu)建多步相關(guān)的噪聲協(xié)方差陣,通過線性變換得到了改進的狀態(tài)協(xié)方差和增益矩陣。該算法能減弱有色噪聲對系統(tǒng)的擾動,但是過程噪聲與觀測噪聲之間的協(xié)方差計算相當復(fù)雜,當考慮相關(guān)步數(shù)增多時算法計算量急劇增大。目前,被動定位背景下的濾波估計算法大多基于EKF或UKF濾波算法,且對傳感器觀測噪聲的適應(yīng)性還有待提升。
本文以基本CKF算法為基礎(chǔ),在建立的被動定位濾波模型基礎(chǔ)上,探索能夠適應(yīng)混合噪聲的CKF算法,以期進一步提升算法對混合觀測噪聲的適應(yīng)性和魯棒性。
選取制導(dǎo)坐標系Oxyz下飛行器- 目標相對位置矢量、相對速度矢量、目標加速度矢量在其3個軸上的投影分量(Δx、Δy、Δz)、(Δvx、Δvy、Δvz)、(atx、aty、atz)作為狀態(tài)向量,表示為
X=[ΔxΔyΔzΔvxΔvyΔvzatxatyatz]T,
式中:X表示系統(tǒng)的狀態(tài)向量;下標t表示目標。則系統(tǒng)狀態(tài)方程表示為
(1)
式中:A為狀態(tài)的轉(zhuǎn)移矩陣;G為飛行器加速度分配矩陣;U表示飛行器加速度矢量在制導(dǎo)坐標系3個軸上的分量,U=[axayaz]T;Γ為系統(tǒng)過程噪聲分配矩陣;W表示過程噪聲,W=[wxwywz]T為零均值高斯分布白噪聲向量,且各分量間彼此相互獨立,wx、wy、wz分別為制導(dǎo)坐標系3個軸上的過程噪聲。
考慮艦船目標的弱機動性,本文中目標模型采用Singer加速度模型[13],用α表示目標機動時間常數(shù)的倒數(shù),即目標機動頻率。則相關(guān)矩陣表示為
(2)
G=[03×3-I3×303×3]T,
(3)
Γ=[03×303×3I3×3]T.
(4)
離散化后的狀態(tài)方程為
xk=Fk|k-1xk-1+Bk-1uk-1+Γk|k-1Wk-1,
(5)
式中:xk表示系統(tǒng)在離散時間k的狀態(tài)向量,k表示離散的時間量;Fk|k-1表示離散時間k-1到k的狀態(tài)轉(zhuǎn)移矩陣,
(6)
T為濾波周期即仿真步長;Bk-1表示離散時間k-1的飛行器加速度分配矩陣,
(7)
uk-1表示離散時間k-1的飛行器加速度矢量在制導(dǎo)坐標系3個軸上的分量;Γk|k-1表示離散時間k-1到k的系統(tǒng)過程噪聲分配矩陣,
(8)
被動傳感器可獲得目標視線高低角qf和視線方位角qh,如圖1所示。
圖1 被動傳感器所測信息Fig.1 Information measured by passive sensor
qf、qh的理論計算公式為
(9)
式中:hf(x)為計算視線高低角理論值的非線性函數(shù);hh(x)為計算視線方位角理論值的非線性函數(shù)。
考慮傳感器的觀測誤差,可建立如下觀測方程:
(10)
基本CKF算法是依據(jù)貝葉斯濾波器的遞推過程,通過假設(shè)各概率密度函數(shù)均服從高斯分布獲得高斯域的貝葉斯濾波框架,在高斯濾波框架上應(yīng)用三度容積規(guī)則近似得到的,文獻[1]給出了其具體推導(dǎo)過程。此處,結(jié)合第1節(jié)建立的被動濾波定位模型,給出被動跟蹤問題中的基本CKF過程。
考慮如(5)式和(10)式建立的非線性隨機系統(tǒng),寫成一般形式如(11)式所示:
(11)
在此條件下,CKF算法的計算流程如下:
2) 時間更新:
k|k-1=Fk|k-1k-1+Bk-1uk-1,
(12)
(13)
3) 觀測更新:
①計算Cubature點:
(14)
(15)
②傳播Cubature點:
Zi,k|k-1=h(Xi,k|k-1),
(16)
(17)
互協(xié)方差陣Pxz,k|k-1和觀測預(yù)測協(xié)方差陣Pzz,k|k-1為
(18)
(19)
于是,可以得到觀測矩陣為
Hk=((Pk|k-1)-1Pxz,k|k-1)T,
(20)
濾波增益為
Kk=Pxz,k|k-1(Pzz,k|k-1)-1,
(21)
以及后驗狀態(tài)估計值及其協(xié)方差陣為
k=k|k-1+Kk(zk-k|k-1),
(22)
由以上濾波過程可知,獲得Pxz,k|k-1之后,可結(jié)合Pk|k-1計算得到觀測矩陣Hk,于是非線性系統(tǒng)中的濾波問題可轉(zhuǎn)化至線性系統(tǒng)中進行分析和處理。
2.2.1 建立有色噪聲模型
傳感器觀測噪聲中,有色隨機噪聲為主要分布噪聲,因此本節(jié)假設(shè)觀測噪聲vk為有色隨機噪聲特征,在基本CKF算法的基礎(chǔ)上推導(dǎo)得到有色噪聲條件下的容積卡爾曼濾波(CKF-CMN)算法。傳感器觀測噪聲模型中對有色噪聲的描述為:正態(tài)分布的白噪聲通過時間常數(shù)為τ的1階濾波器后的輸出(標準差為σ)表示為[σ,τ]。根據(jù)歐拉法求解微分方程可以獲得有色噪聲序列,同時可以將有色噪聲寫成如下形式:
vk=(1-T/τ)·vk-1+(T/τ)·rk,
(23)
式中:rk表示1階濾波器輸入的正態(tài)分布白噪聲,rk的方差為Rr.將(23)式寫成(24)式所示的1階馬爾可夫形式:
vk=φvk-1+ξk,
(24)
式中:φ=1-T/τ;ξk=(T/τ)·rk,ξk的方差Rk=(T/τ)2·Rr.
2.2.2 有色噪聲的白化
考慮與(11)式等價的線性系統(tǒng):
(25)
式中:wk-1與(11)式中的Γk|k-1Wk-1等價,vk=φvk-1+ξk.為了使CKF算法在有色觀測噪聲條件下仍然保持優(yōu)良的濾波性能,首先對有色噪聲進行白化,構(gòu)造新的觀測量zk-φzk-1,經(jīng)等價替換及整理得
zk-φzk-1=(HkFk|k-1-φHk-1)xk-1+
HkBk-1uk-1+Hkwk-1+ξk.
(26)
令
(27)
(28)
(29)
則觀測方程可以改寫為
(30)
(31)
(32)
(33)
通過以上分析可知,將有色噪聲白化后,系統(tǒng)的觀測噪聲模型便符合基本KF器的觀測噪聲為白噪聲這一假設(shè),然而與基本KF假設(shè)不同的是,此時觀測噪聲和系統(tǒng)噪聲之間存在一定的相關(guān)性。下面考慮采用待定系數(shù)去相關(guān)法對上述相關(guān)性進行消除。
2.2.3 基于待定系數(shù)法的噪聲去相關(guān)方法
完成有色觀測噪聲白色化后,由(30)式可將(25)式等效為(34)式所示的一般形式:
(34)
(34)式中的狀態(tài)方程改寫為
(35)
(36)
(37)
選擇
(38)
(39)
(40)
綜上所述,可得到CKF-CMN算法計算過程如下:
3)根據(jù)(23)式~(40)式,構(gòu)造新的觀測方程和狀態(tài)方程,按以下步驟對濾波過程進行更新:
濾波增益:
(41)
狀態(tài)估計及其協(xié)方差:
(42)
針對混合觀測噪聲中的跳變噪聲以及有色噪聲可能時變給CKF-CMN算法帶來的影響,在CKF-CMN算法基礎(chǔ)上設(shè)計自適應(yīng)濾波算法,實時對混合噪聲的噪聲協(xié)方差進行估計,并對觀測噪聲中的有害信息進行剔除。
2.3.1 對噪聲協(xié)方差的估計
(43)
(44)
(45)
2.3.2 有害觀測信息的剔除
H0:觀測信息有效;
H1:觀測信息無效,剔除。
采用檢驗統(tǒng)計量Dk
(46)
圖2 飛行器與目標運動軌跡Fig.2 Moving trajectories of aircraft and target
本文仿真時間為100 s,仿真步長為0.01 s,蒙特卡洛仿真次數(shù)N為50次。假設(shè)被動傳感器測角噪聲為時變有色噪聲和跳變噪聲的混合噪聲。對于有色噪聲[σ,τ](參考2.2節(jié),σ、τ分別表示其標準差和時間常數(shù))。某次仿真中其變化特性如圖3所示,圖中Tc表示噪聲方差發(fā)生變化的時刻,在Tc=50 s處,有色噪聲特性由[0.5°,0.02 s]變?yōu)閇1.0°,0.02 s],跳變噪聲幅值為1.5°,周期為5 s,某次仿真中其變化特性如圖4所示,二者的混合噪聲如圖5所示。
圖3 時變有色噪聲Fig.3 Time-varying colored noise
圖4 跳變噪聲Fig.4 Jumping noise
圖5 被動傳感器觀測噪聲隨時間變化Fig.5 Passive sensor observation noise vs. time
在觀測噪聲為混合噪聲條件下采用基本CKF、CKF-CMN以及ACKF算法在同一條件下進行對比仿真,結(jié)果如圖6~圖10所示,圖中縱軸為目標位置分量的均方根誤差(RMSE),即x軸、y軸和z軸方向目標位置真實值和估計值的均方根誤差(RMSEx、RMSEy、RMSEz),RMSEx、RMSEy、RMSEz根據(jù)(47)式計算得到:
圖6 x軸方向被動定位估計誤差Fig.6 Estimated errors of passive positioning in x direction
圖10 z軸方向被動定位估計誤差(放大圖)Fig.10 Estimated errors of passive positioning in z direction (enlarged view)
(47)
圖7 x軸方向被動定位估計誤差(放大圖)Fig.7 Estimated errors of passive positioning in x direction (enlarged view)
圖8 y軸方向被動定位估計誤差Fig.8 Estimated errors of passive positioning in y direction
圖9 z軸方向被動定位估計誤差Fig.9 Estimated errors of passive positioning in z direction
取圖中第1~100 s的數(shù)據(jù)進行計算,得到x軸、y軸、z軸3個方向RMSE的平均值如表1所示。
由圖6~圖10可知:當被動定位系統(tǒng)的觀測數(shù)據(jù)中存在時變有色噪聲及跳變噪聲時,基本KF算法的精度和適用性受到明顯制約,難以實現(xiàn)高精度的目標位置估計。在觀測量包含混合噪聲時,3種被動定位算法的精度依次為CKF算法 表1 3種濾波定位算法精度分析 本文從混合觀測噪聲條件下的被動定位精度和濾波穩(wěn)定性提升需求出發(fā),在分析CKF算法的基礎(chǔ)上,通過重新構(gòu)建觀測方程和狀態(tài)方程提高所建立的濾波器對混合噪聲的適應(yīng)性。得出主要結(jié)論如下: 1)針對混合噪聲中的有色噪聲使得基本CKF算法定位精度下降的問題,設(shè)計了CKF-CMN算法。通過觀測重構(gòu)、待定系數(shù)去相關(guān)的思想,構(gòu)建新的濾波狀態(tài)方程和觀測方程,將有色噪聲轉(zhuǎn)變?yōu)榘自肼暎沟脼V波估計精度得到提升。 2)針對有色噪聲的時變性以及跳變噪聲使得濾波性能受損的問題,在CKF-CMN算法的基礎(chǔ)上設(shè)計了ACKF算法,采用實時噪聲方差估計和有害信息剔除法使得被動定位精度得到進一步提升。 3)50次蒙特卡洛仿真實驗結(jié)果表明,ACKF算法的被動濾波定位精度得到提升且魯棒性更好,基本CKF、CKF-CMN、ACKF 3種算法的定位精度依次為CKF算法4 結(jié)論