彈道再入目標(biāo)軌跡跟蹤中非線性濾波算法研究＊

田亞菲 莫 驊 于 飛

（1．蘭州大學(xué)信息學(xué)院 蘭州 730000）（2．總參謀部通信訓(xùn)練基地 張家口 075100）

1 引言

彈道再入目標(biāo)的軌跡跟蹤問題一直以來受到人們的普遍關(guān)注，它是彈道導(dǎo)彈防御系統(tǒng)的核心技術(shù)，直接影響對(duì)來襲導(dǎo)彈攔截的成功率。針對(duì)跟蹤目標(biāo)的動(dòng)力學(xué)模型和雷達(dá)的觀測(cè)模型具有典型的非線性特性，國內(nèi)外相關(guān)文獻(xiàn)提出許多非線性濾波算法，在算法選擇上存在明顯不同的觀點(diǎn)［1～3］。常見的算法中，大致可分為兩大類［4］：一類是基于點(diǎn)的濾波算法，主要包括擴(kuò)展卡爾曼濾波器（EKF）和無跡卡爾曼濾波器（UKF）及對(duì)應(yīng)的改進(jìn)算法；另一類是基于密度的濾波算法，包括貝葉斯粒子濾波器（PF）及對(duì)應(yīng)的改進(jìn)算法。

本文依據(jù)彈道再入目標(biāo)動(dòng)力學(xué)模型和雷達(dá)觀測(cè)模型特點(diǎn)，從算法基礎(chǔ)理論出發(fā)對(duì)比分析EKF和UKF算法的異同點(diǎn)，通過仿真實(shí)驗(yàn)全面比較了EKF、UKF和PF三種基本算法的性能，得出到在高斯噪聲條件下，EKF要優(yōu)于UKF和PF，在閃爍噪聲條件下，PF要優(yōu)于EKF和UKF的結(jié)論。

2 基于點(diǎn)的濾波算法理論

2．1 線性遞歸估計(jì)算法的基本思想

本質(zhì)上講，EKF和UKF在狀態(tài)更新階段是一個(gè)基于最小均方誤差（MMSE）條件下的線性遞歸估計(jì)算法。一個(gè)非線性動(dòng)態(tài)系統(tǒng)的狀態(tài)方程和觀測(cè)方程可用下面兩式來描述。

使用線性回歸估計(jì)方法可以從當(dāng)前觀測(cè)值z(mì)k中估計(jì)出k時(shí)刻的狀態(tài)變量xk，其中經(jīng)典的算法是線性最小均方誤差估計(jì)法（LMMSE）。觀測(cè)值序列的概率密度函數(shù)未知的平穩(wěn)過程在前二階矩已知的條件下，可以通過遞歸迭代的算法完成估計(jì)，此時(shí)基于觀測(cè)值的估計(jì)是觀測(cè)值的線性變換。若符合高斯約束條件時(shí)，LMMSE算法是基于MMSE約束條件下的最優(yōu)估計(jì)。

用zk表示從初始時(shí)刻到k時(shí)刻所得到的一組觀測(cè)數(shù)據(jù)，LMMSE 估計(jì)中E＊xk｜zk具有如下解析形式［5］：

如果以上式中狀態(tài)一步預(yù)測(cè)、一步預(yù)測(cè)誤差協(xié)方差，觀測(cè)值預(yù)測(cè)、觀測(cè)預(yù)測(cè)估計(jì)協(xié)方差Sk，狀態(tài)估計(jì)值與觀測(cè)估計(jì)值互協(xié)方差都能通過給定的k-1時(shí)刻的狀態(tài)估計(jì)值和誤差協(xié)方差Pk-1｜k-1（即初值）求出，此LMMSE算法就是真正意義上的遞歸運(yùn)算。若系統(tǒng)的動(dòng)力學(xué)方程和觀測(cè)方程能滿足線性和高斯條件，以上計(jì)算公式就是著名的kalman濾波算法。

對(duì)于具有觀測(cè)方程（2）的非線性系統(tǒng)

和都依賴于k-1時(shí)刻以前的觀測(cè)值z(mì)k-1和相應(yīng)時(shí)刻LMMSE 估計(jì)值。為實(shí)現(xiàn)真正的遞歸線性估計(jì)，要求k時(shí)刻的估計(jì)值僅從和Pk-1｜k-1就能求得，即：

以上兩式中，Pred f（·），，Pk-1｜k-1表示｛，Pk-1｜k-1｝通過非線性函數(shù)f（·）的傳遞作用來近似得到一步狀態(tài)預(yù)測(cè)值和對(duì)應(yīng)的協(xié)方差Pk｜k-1。Pred h（·），，Pk｜k-1表示 通過狀態(tài)預(yù)測(cè)值來近似得到觀測(cè)預(yù)測(cè)值和相應(yīng)的誤差矩陣。近似運(yùn)算由于舍棄了部分信息，從而導(dǎo)致信息傳遞的不完整，使非線性濾波器的性能發(fā)生退化。通常采用兩種方法來克服近似運(yùn)算帶來的弊端：一是盡可能地使用最優(yōu)的近似，這樣就能避免過多的舍棄信息，保證預(yù)測(cè)的準(zhǔn)確性，例如采用高階的泰勒級(jí)數(shù)展開方法（TSE）；另外一種方法是用一種適當(dāng)?shù)姆植紒肀平蔷€性函數(shù)而不進(jìn)行舍棄，這樣的預(yù)測(cè)中包含了所有非線性系統(tǒng)的統(tǒng)計(jì)信息，例如無跡變換方法（UT）［6］。

無論觀測(cè)方程是否滿足線性條件，只要能實(shí)現(xiàn)非線性系統(tǒng)的近似線性化，得到預(yù)測(cè)值的一階矩和二階矩信息，LMMSE估計(jì)器就可以非常理想地直接應(yīng)用在非線性系統(tǒng)的狀態(tài)更新上，這就是EKF和UKF本質(zhì)上相同之處。兩者的區(qū)別僅在于在式（7）～式（8）中對(duì)非線性系統(tǒng)的近似處理上采用不同的策略。

2．2 EKF的線性化策略

EKF的近似策略是對(duì)非線性函數(shù)的Taylor展開式進(jìn)行一階或二階的線性化截?cái)?，忽略高階項(xiàng)。假設(shè)過程噪聲和量測(cè)噪聲是0均值且互不相關(guān)的的高斯噪聲，在一階EKF中，非線性系統(tǒng)的狀態(tài)方程（1）和觀測(cè)方程（2）通過以下的步驟完成線性化：

F和H分別是函數(shù)f（·）和h（·）對(duì)狀態(tài)變量x求偏導(dǎo)的雅克比矩陣，即和W和V分別是函數(shù)f（·）和h（·）分別對(duì)過程噪聲w和量測(cè)噪聲v求偏導(dǎo)數(shù)的雅克比矩陣，即W［i，j］＝根據(jù)以上線性化變換，預(yù)測(cè)及狀態(tài)更新過程描述如下

式（11）～式（15）與LMMSE 的遞歸線性運(yùn)算式共同構(gòu)成了一階EKF算法公式。

EKF的估計(jì)精度嚴(yán)重依賴于系統(tǒng)的非線性程度。為提高近似精度，可采用二階EKF。雖然二階EKF可以比一階EKF獲得更好的濾波性能，但是計(jì)算量較大，同時(shí)要求非線性函數(shù)的海森矩陣在數(shù)學(xué)上存在，而在計(jì)算海森矩陣中會(huì)產(chǎn)生計(jì)算誤差。因此，針對(duì)彈道再入目標(biāo)觀測(cè)方程的強(qiáng)非線性特性，在EKF線性化近似之前，可以采用坐標(biāo)系轉(zhuǎn)換的方式來盡可能地減弱其非線性特性［7］。對(duì)于極坐標(biāo)下所描述的觀測(cè)方程，可以通過以下兩式轉(zhuǎn)化成直角坐標(biāo)系的方程［8］。

轉(zhuǎn)化后的觀測(cè)變量成為偽線性的形式，減弱了雷達(dá)觀測(cè)變量在初始極坐標(biāo)系下的非線性程度。

2．3 UKF的線性化策略

在基于近似非線性函數(shù)的概率分布比近似其本身更容易［9～10］的思想指導(dǎo)下，Julier等提出了基于unscented transformation的Kalman濾波算法，即UKF。x為nx維的隨機(jī)變量，y＝g（x）是一個(gè)非線性變換，假設(shè)x的均值和方差分別表示為和Cx，算法描述如下：

第一步：計(jì)算sigma點(diǎn)的值χi及其權(quán)重Wi：

其中κ是縮放比例系數(shù)，是矩陣（nx＋κ），Cx的方根的第i行（列）。

第二步：通過非線性函數(shù)傳播sigma點(diǎn)：

第三步：計(jì)算y的均值和方差，以及x和y的協(xié)方差：

非線性系統(tǒng)的狀態(tài)方程（1）和觀測(cè)方程（2）通過UT 完成線性化變換：

和分別表示在k-1 時(shí)刻預(yù)測(cè)值和更新值，和是在k-1時(shí)刻過程噪聲和觀測(cè)噪聲。狀態(tài)更新階段同EKF一樣使用LMMSE估計(jì)算法。

盡管采用有限個(gè)采樣點(diǎn)簡(jiǎn)化了近似過程，但是計(jì)算采樣點(diǎn)的方根運(yùn)算增加了計(jì)算量，且根式中的矩陣一旦失去正定勢(shì)必導(dǎo)致濾波器發(fā)散。

3 基于密度的PF算法理論

Gordon和Salmond在1993年提出了一種基于序列重要抽樣濾波思想的Bootstrap非線性濾波方法［11］，從此奠定了粒子濾波（PF）算法的基礎(chǔ)。該算法的基本原理是對(duì)于一個(gè)平穩(wěn)的時(shí)變動(dòng)態(tài)系統(tǒng)，設(shè)k-1時(shí)刻系統(tǒng)的后驗(yàn)概率密度為p（xk-1｜zk-1），根據(jù)一定的原則選擇n個(gè)隨機(jī)樣本點(diǎn)（粒子），k時(shí)刻獲得觀測(cè)信息后，經(jīng)過預(yù)測(cè)更新和狀態(tài)更新，n個(gè)粒子的后驗(yàn)概率密度近似為p（xk｜zk），隨著粒子數(shù)目的增加，粒子的概率密度函數(shù)就逐漸逼近狀態(tài)變量的概率密度函數(shù)，此時(shí)粒子濾波器的估計(jì)就達(dá)到了最優(yōu)貝葉斯估計(jì)的效果［12］。

基本的粒子濾波算法流程如下：

第一步初始化：取k＝0，按p（x0）抽取N個(gè)樣本點(diǎn)，i＝1，2，…，N；

第四步歸一化權(quán)重：

第五步重采樣：根據(jù)各自歸一化權(quán)重的大小復(fù)制或者舍棄樣本，得到N個(gè)近似服從p分布的樣本。

進(jìn)博會(huì)食品及農(nóng)產(chǎn)品展區(qū)6萬平方米，分為乳制品、蔬果農(nóng)產(chǎn)品、肉制品和水產(chǎn)品、休閑食品甜食調(diào)味品、酒類和飲料等五大專業(yè)展區(qū)。徜徉其間，從參觀者身上你會(huì)感受到中國對(duì)美好生活的追求和向往是那么強(qiáng)烈！中國當(dāng)自強(qiáng)！中國農(nóng)業(yè)當(dāng)自強(qiáng)！

第七步k＝k＋1，重復(fù)第二至六步進(jìn)行迭代運(yùn)算直到結(jié)束。

PF算法通過尋找一組在狀態(tài)空間傳播的隨機(jī)樣本來對(duì)概率密度函數(shù)完成近似，以樣本均值代替積分運(yùn)算，從而獲得狀態(tài)變量最小方差分布。它與EKF、UKF不同，對(duì)狀態(tài)變量沒有高斯約束，因此近年來成為了非線性非高斯系統(tǒng)狀態(tài)估計(jì)的“最優(yōu)”估計(jì)器。但是PF算法濾波精度與粒子數(shù)量、重要性函數(shù)、重采樣方法的選取有很大關(guān)系［13］。近年來國內(nèi)外出現(xiàn)了大量的改進(jìn)算法，但計(jì)算量巨大、實(shí)時(shí)性差的缺點(diǎn)始終未能克服，在一定程度上限制了其在彈道再入目標(biāo)軌跡跟蹤問題上的工程應(yīng)用。

4 仿真研究

4．1 彈道再入目標(biāo)的數(shù)學(xué)模型及仿真環(huán)境描述

在不影響模型精確度的前提下，假定地球表面為平面，目標(biāo)為質(zhì)點(diǎn)且不考慮自旋。相關(guān)文獻(xiàn)對(duì)彈道再入目標(biāo)的數(shù)學(xué)模型進(jìn)行了較為系統(tǒng)的分析［7，14］。本文用如下兩式描述動(dòng)力學(xué)方程和觀測(cè)方程。

xk＝［xk，˙xk，yk，˙yk］′是狀態(tài)變量，各元素分別表示目標(biāo)在直角坐標(biāo)系下x軸和y軸方向上的位置和速度分量。zk＝［rk，θk］′是觀測(cè)變量，各元素分別表示雷達(dá)對(duì)目標(biāo)的測(cè)量距離和角度。wk是過程噪聲，代表觀測(cè)噪聲，各元素分別是距離和角度的觀測(cè)誤差。假設(shè)過程噪聲和觀測(cè)噪聲互不相關(guān)，且服從不依賴于初值x0的0均值高斯分布，它們的協(xié)方差矩陣通過以下三個(gè)矩陣式描述，其中Q是過程噪聲協(xié)方差矩陣，q是過程噪聲的方差，R是觀測(cè)噪聲協(xié)方差矩陣。

式（18）中非線性函數(shù)Ψ（·）定義為

式（20）右端后兩項(xiàng)描述的是空氣阻力和重力對(duì)目標(biāo)狀態(tài)變量的影響作用。g是重力加速度，f（xk，β）是空氣動(dòng)力阻力函數(shù)，β為彈道系數(shù)，一般的，該系數(shù)是時(shí)間變量t的函數(shù)，即β（t）。文獻(xiàn)［15］討論了參數(shù)β（t）的建模方式對(duì)非線性濾波算法的性能影響。為簡(jiǎn)化模型，本文采用中的彈道系數(shù)為確知參數(shù)。式（20）中的狀態(tài)轉(zhuǎn)移矩陣F和驅(qū)動(dòng)矩陣G分別表示為

空氣動(dòng)力阻力函數(shù)f（xk，β）定義為式（21）。

ρ（·）是隨海拔高度呈指數(shù)衰落的空氣密度函數(shù)。觀測(cè)方程中的非線性測(cè)量函數(shù)h（·）表示為

仿真中雷達(dá)掃描周期T＝2s，仿真跟蹤步數(shù)60步，目標(biāo)彈道系數(shù)β＝40000kg／ms2，蒙特卡洛仿真次數(shù)為1000次，目標(biāo)初始狀態(tài)［332000m，2290cos（190°）m／s，98000m，2290sin（190°）m／s］，高斯噪聲環(huán)境下觀測(cè)噪聲的標(biāo)準(zhǔn)差分別是σr＝100m 和σε＝0．05rad，過程噪聲方差q＝2。

在高斯噪聲條件和閃爍噪聲條件下，分別使用EKF、UKF及PF進(jìn)行濾波，通過比較不同方向上位置、速度的估計(jì)精度和運(yùn)算時(shí)間來對(duì)比三種算法的基本性能。EKF采用坐標(biāo)變換的方法減弱觀測(cè)方程的非線性特性，使用一階近似，PF 中粒子數(shù)為500，采用殘差重采樣方法。

4．2 仿真結(jié)果及分析

圖1 三種算法在X 軸上位置的RMSE

圖2 三種算法在X 軸上速度的RMSE

圖3 三種算法在Y 軸上位置的RMSE

圖4 三種算法在Y 軸上速度的RMSE

圖1比較了高斯條件下三種濾波算法在X和Y軸方向上位置與速度的均方根誤差（RMSE），由于各濾波器的初值已給定，從初始時(shí)刻開始三種濾波器就能得到良好的跟蹤效果。圖1中，EKF 和PF對(duì)X軸方向上位置變量的跟蹤精度幾乎一致，要略優(yōu)于UKF。從理論上講，如果適當(dāng)增加PF的粒子數(shù)量和嘗試不同的重采樣算法，會(huì)提高其跟蹤精度，但是會(huì)帶來運(yùn)算量的劇增。在仿真的第60步，三種濾波器的結(jié)果趨于穩(wěn)定狀態(tài)，分別收斂于60m、70m 和120m 的誤差水平。圖2表明，EKF和PF對(duì)X方向上速度估計(jì)的性能要優(yōu)于UKF，誤差僅是后者的1／3左右。圖3中盡管UKF在第40步左右出現(xiàn)190m 的誤差，但是三種濾波器在第50步后，誤差穩(wěn)定于160m，可以認(rèn)為三種濾波器的性能相當(dāng)。圖4可以得到和圖2基本相同的結(jié)論。

通過對(duì)仿真結(jié)果的分析可看出PF 由于在粒子數(shù)量和重采樣方法選擇策略兩個(gè)條件限制下，并未能達(dá)到理想的跟蹤精度，但仍體現(xiàn)出良好的性能，基本與EKF 的性能相當(dāng)。相同條件的仿真中EKF則要優(yōu)于UKF，這并非與文獻(xiàn)［10］的論述相悖，只是由于在彈道再入目標(biāo)軌跡跟蹤問題上，若彈道系數(shù)預(yù)先確知，對(duì)觀測(cè)方程采用適當(dāng)?shù)膫尉€性變換，削弱其非線性程度，EKF 是可以達(dá)到了較為精確的濾波效果的。理論上講，UKF 能夠得到比EKF更高的近似精度，而并非在任何情形下都要優(yōu)于EKF。UT 作為一種線性化的技術(shù)，只是為了獲得濾波器狀態(tài)更新時(shí)必要的信息，而并不能完全取代其它的線性化技術(shù)。

表1給出了MATLAB 中測(cè)量的三種濾波器相對(duì)于EKF在一次蒙特卡洛仿真中的運(yùn)算時(shí)間。

表1 三種算法的相對(duì)運(yùn)行時(shí)間

三種濾波算法中，UKF的運(yùn)算時(shí)間為EKF 的1．85 倍，運(yùn)行時(shí)間與EKF 相當(dāng)，估計(jì)精度不及EKF。PF盡管在運(yùn)行精度上和EKF相當(dāng)，運(yùn)行時(shí)間卻是其91．05倍。

表2比較了閃爍噪聲環(huán)境下三種濾波算法的X 和Y 軸方向上位置和速度的均方根誤差的均值（ARMSE）。閃爍噪聲采用高斯噪聲和具有“厚尾”特性的拉普拉斯分布噪聲加權(quán)和來實(shí)現(xiàn)。在嚴(yán)重閃爍噪聲干擾條件下，仿真程序中閃爍噪聲權(quán)重Weight＝0．9。

表2中，閃爍噪聲條件下EKF 出現(xiàn)了極大偏差，UKF發(fā)散，濾波器無法正常工作，因此基于點(diǎn)的濾波器在此條件下失去了對(duì)目標(biāo)的跟蹤性能。PF的ARMSE盡管相對(duì)于高斯環(huán)境下略有升高，但仍然收斂。

表2 三種算法在閃爍噪聲環(huán)境下的ARMSE

5 結(jié)語

本文基于彈道再入目標(biāo)的軌跡跟蹤問題，比較了EKF、UKF 和PF 的性能。通過分析LMMSE濾波器的基本原理，可知所有基于點(diǎn)的濾波算法在狀態(tài)更新階段的本質(zhì)是相同的，僅在線性化策略上有所區(qū)別。盡管在高斯噪聲條件下非線性濾波問題中，UKF要比EKF 具有更高階的信息，其性能優(yōu)于EKF，但在彈道再入目標(biāo)軌跡跟蹤問題中，如果目標(biāo)的彈道系數(shù)確知，EKF 具有更高的估計(jì)精度和更短的運(yùn)算時(shí)間，能夠接近優(yōu)化后的PF 算法的性能。因此高斯條件下在濾波器的選擇上本文更傾向于EKF 算法。閃爍噪聲條件下，基于高斯假設(shè)的基于點(diǎn)的濾波算法失去了理論基礎(chǔ)，濾波器性能出現(xiàn)偏差甚至是發(fā)散，此時(shí)基于密度的PF 算法就成為了“最優(yōu)”的估計(jì)器，同時(shí)對(duì)彈道再入目標(biāo)滯空時(shí)間短暫與PF 算法運(yùn)算時(shí)間長的矛盾應(yīng)該予以重視。

應(yīng)當(dāng)指出的是，初始狀態(tài)對(duì)非線性濾波器的穩(wěn)定性非常關(guān)鍵，彈道再入目標(biāo)的初始狀態(tài)的確定方法對(duì)濾波算法的性能會(huì)產(chǎn)生很大影響，這個(gè)問題有待進(jìn)一步研究。當(dāng)彈道系數(shù)作為狀態(tài)變量的一個(gè)待估參量時(shí)，其建模方法對(duì)非線性濾波器的性能會(huì)產(chǎn)生影響。下一步將繼續(xù)研究這些因素對(duì)彈道目標(biāo)軌跡跟蹤問題中非線性濾波器選擇策略的影響。

［1］Farina A，Ristic B，Benvenuti D．Tracking a Ballistic Target：Comparison of Several Nonlinear Filters［J］．IEEE Transactions on Aerospace and Electronic Systems，2002，38（3）：854-867．

［2］江曉東，謝京穩(wěn)，郭軍海．三種非線性濾波在再入彈道估計(jì)中的分析研究［J］．戰(zhàn)術(shù)導(dǎo)彈技術(shù)，2012（1）：43-49．

［3］劉鑫，郭云飛，薛安克．幾種濾波器跟蹤性能的比較［J］．計(jì)算機(jī)與數(shù)字工程，2013（3）：391-393．

［4］Zhao Z L，Chen H M，Chen G S，et al．Comparison of Several Ballistic Target Tracking Filters［M］．New York：IEEE，2006：2197-2202．

［5］Zhao Z L，Li X R，Jilkov V P．Best Linear Unbiased Filtering with Nonlinear Measurements forTarget Tracking［J］．Aerospace and Electronic Systems，IEEE Transactions on，2004，40（4）：1324-1336．

［6］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking：Approximation Techniques for Nonlinear Filtering［C］／／Pefense and Security．Intemational Society for Optics and Photonics，2004：537-550．

［7］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking—PartⅢ：Measurement Models［C］／／Internation Symposium on Optical Science and Technology．International Society for Optics and Photonics，2001：423-446．

［8］Yaakov B S，Li X R．Estimation with Applications to Tracking and Navigation：Theory Algorithms and Software［M］．New York：John Wiley ＆Sons，Inc．，2001．

［9］Julier S J，Uhlmann J K．Unscented Filtering and Non-linear Estimation［J］．Proceedings of the IEEE，2004，92（3）：401-422．

［10］Julier S J，Uhlmann J K，Whyte H F．A New Method for the Nonlinear Transformation of Means and Covariances inFilters and Estimators［J］．Automatic Control，IEEE Transactions on，2000，45（3）：477-482．

［11］Gordon N J，Salmond D J，Smith A F M．Novel Approach to Nonlinear／non-Gaussian BayesianState Estimation［J］．Iee Proceedings F Radar and Signal Processing，1993，140（2）：107-113．

［12］Doucet A，F(xiàn)reitas N，Gordon N．Sequential Monte Carlo Methods in Practice［M］．New York：Springer-Verlag，2001．

［13］馮馳，王萌，汲清波．粒子濾波器重采樣算法的分析與比較［J］．系統(tǒng)仿真學(xué)報(bào)，2009（4）：1101-1105．

［14］Li X R，Jilkov V P．A Survey of Maneuvering Target Tracking-partⅡ：Ballistic Target Models［C］／／Proceedings of 2001SPIE Conference on Signal and Data Processing of Small Targets，2001，4473：559-581．

［15］Ristic B，F(xiàn)arina A，Benvenuti D．Tracking a Ballistic Object on Reentry：Performance Bounds and Comparison of NonlinearFilters［C］／／Radar，Sonar and Navigatio，IEE Proceedings，2003，150（2）：65-70．