李春月，廖育榮，倪淑燕，陳 帥

0 引 言

彈道導(dǎo)彈具有射程遠(yuǎn)、速度快、精度高、突防能力強(qiáng)、殺傷威力大、效費比高等優(yōu)點，已成為現(xiàn)代戰(zhàn)爭非常重要的進(jìn)攻性武器。為了遏制日益加劇的彈道導(dǎo)彈威脅，世界各國都在大力發(fā)展導(dǎo)彈防御系統(tǒng)。雷達(dá)是導(dǎo)彈防御系統(tǒng)中的核心探測器，對彈道導(dǎo)彈的實時跟蹤精度直接影響到導(dǎo)彈防御系統(tǒng)對導(dǎo)彈攔截的成功率。因此，高精度彈道導(dǎo)彈實時跟蹤算法是目前研究的熱點問題。

彈道目標(biāo)實時跟蹤問題本質(zhì)上為高維非線性系統(tǒng)的最優(yōu)狀態(tài)估計問題。卡爾曼濾波是目前最優(yōu)狀態(tài)估計中應(yīng)用最為廣泛的一種算法，經(jīng)典的非線性卡爾曼濾波算法主要有擴(kuò)展卡爾曼濾波（Extended KalmanFilter，EKF）算法[1]和無跡卡爾曼濾波（Unscented Kalman Filter，UKF）算法[2,3]。其中，EKF算法的核心思想是將非線性的系統(tǒng)狀態(tài)函數(shù)做線性化處理，具體方法是利用泰勒級數(shù)展開法對非線性狀態(tài)方程進(jìn)行一階展開，然后再用離散卡爾曼濾波進(jìn)行最優(yōu)狀態(tài)估計，這種算法對非線性較強(qiáng)的系統(tǒng)估計精度不夠高；UKF算法的出發(fā)點是基于“對概率分布進(jìn)行近似要比對非線性函數(shù)近似容易很多”的思想，采用Sigma點的分布近似表示非線性函數(shù)的分布，有效提高了估計精度，但是對于高維非線性系統(tǒng)（維數(shù)n≥4），UKF算法中的自由調(diào)節(jié)參數(shù)κ＜ 0，使得中心采樣點的權(quán)值κ＜ 0，從而使 UKF 算法在濾波過程中可能會出現(xiàn)協(xié)方差非正定情況，導(dǎo)致最終收斂結(jié)果不穩(wěn)定甚至發(fā)散，并且隨著系統(tǒng)維數(shù)的增加，采樣點與中心點距離不斷增大，導(dǎo)致非局部效應(yīng)，嚴(yán)重影響 UKF 算法的估計精度。2009年，Arasaratnam等人[4,5]采用Spherical -Radial規(guī)則近似非線性函數(shù)傳遞的后驗均值和協(xié)方差的方法，依據(jù)高斯濾波框架提出容積卡爾曼濾波（Cubature Kalman Filter，CKF）算法，相比于UKF算法，CKF算法有嚴(yán)格的數(shù)學(xué)推導(dǎo)過程，降低了運算量，提升了估計精度，得到廣泛應(yīng)用[6～9]。文獻(xiàn)[10]對UKF算法和CKF算法的估計精度和數(shù)值穩(wěn)定性做了比較，指出當(dāng)系統(tǒng)狀態(tài)維數(shù)小于3維時，UKF算法估計精度高于CKF算法，當(dāng)系統(tǒng)狀態(tài)維數(shù)等于3維時，二者估計精度相當(dāng)，當(dāng)系統(tǒng)狀態(tài)維數(shù)大于3維時，無論估計精度和數(shù)值穩(wěn)定性，CKF算法都高于 UKF算法；受到Spherical -Radial規(guī)則的啟發(fā)，文獻(xiàn)[11]提出了更多容積點的CKF算法，以提高估計精度，但是也相應(yīng)增加了計算量；文獻(xiàn)[12]為了進(jìn)一步提高 CKF 算法的估計精度，提出迭代 CKF 算法，并將其應(yīng)用于再入彈道目標(biāo)狀態(tài)估計，其效果優(yōu)于 CKF 算法。

針對如何進(jìn)一步提高對再入彈道目標(biāo)的實時跟蹤精度，本文提出將另一種基于Spherical Simplex-Radial準(zhǔn)則的SSRCKF算法用于再入彈道目標(biāo)實時跟蹤中。該算法將Spherical Simplex準(zhǔn)則引入CKF算法，以一種新的容積點選擇方法來計算球面積分，進(jìn)而近似計算高斯域下的貝葉斯積分，有效改善了CKF算法的估計精度，將其用于再入彈道目標(biāo)跟蹤中，提高彈道目標(biāo)的實時跟蹤精度，最后通過仿真分析，驗證了算法的有效性。

1 彈道目標(biāo)跟蹤數(shù)學(xué)模型

1.1 再入段動力學(xué)模型

對再入彈道目標(biāo)進(jìn)行動力學(xué)建模，首先，對再入彈道目標(biāo)進(jìn)行受力分析。設(shè)彈道目標(biāo)一般保持零度攻角再入，受到的空氣動力表現(xiàn)為大氣阻力，大氣阻力加速度方向與再入速度方向相反。定義以彈道測量雷達(dá)為原點的北天東（North-Up-East）坐標(biāo)系，O-xyz，Ox指向正北方向，Oy垂直于地球表面指向正上方，Oz指向正東方向，由于再入彈道導(dǎo)彈速度快，再入時間短，忽略地球自轉(zhuǎn)角速度的影響。對再入目標(biāo)建模如圖 1所示，假設(shè)地球為標(biāo)準(zhǔn)球體，EO為地球質(zhì)心，ag和af分別為引力加速度和大氣阻力加速度，彈道目標(biāo)的位置、速度分量和彈道系數(shù)構(gòu)成七維狀態(tài)向量，給出再入彈道目標(biāo)的系統(tǒng)狀態(tài)方程。

圖1 雷達(dá)觀測再入彈道目標(biāo)示意Fig.1 Radar Observation Reentry Ballistic Target Signal

式 中 [ω1( t ) … ω7(t )]Τ為 系 統(tǒng) 狀 態(tài) 噪 聲 ， 滿 足ω～0(,)Q的統(tǒng)計特性 ；μ為地球重力常數(shù)，μ= G M = 3 .986× 1 014m3/s2；為 地 球 半 徑 ，Re= 6 378137m ； R(t)為目標(biāo)到地心距離，()Vt為目標(biāo)再入速度，ρ(t)為大氣密度模型，ρ(t ) = ρ0exp ( ?(R (t) ? Re)H0)，H0為大氣密度標(biāo)高，H0=6700km，ρ0為海平面的大氣密度，ρ0= 1 .225kg/m3；β( t)為彈道系數(shù)模型，有[13]：

式中DC為氣動阻力系數(shù)；m為彈道目標(biāo)質(zhì)量；A為彈道目標(biāo)有效截面積。為了保證彈道系數(shù)的非負(fù)性，防止濾波器發(fā)散采用式（2）的指數(shù)模型；()tγ的變化率可以用零均值高斯白噪聲表示，即：

由此可以得到再入彈道目標(biāo)的7維狀態(tài)方程，狀態(tài)量，為了保證濾波器的精度，采用四階龍格-庫塔方法對狀態(tài)方程進(jìn)行離散，最終得到離散非線性系統(tǒng)狀態(tài)方程：

1.2 量測模型

以觀測雷達(dá)為坐標(biāo)系原點，在北天東坐標(biāo)系下建立再入彈道目標(biāo)量測模型。如圖1所示，雷達(dá)對彈道目標(biāo)的距離為 r，俯仰角為η，方位角為ε，有z = [r,ε,η]Τ。由觀測值相對彈道目標(biāo)位置的幾何關(guān)系建立非線性量測方程：

式中 vk為雷達(dá)測量噪聲向量，，滿足v～0(,)R的統(tǒng)計特性，且kv，kω和kx互不相關(guān)。由此得到系統(tǒng)的量測方程：

2 球面單徑容積卡爾曼濾波跟蹤算法

2.1 Spherical Simplex-Radial準(zhǔn)則

高斯域下的貝葉斯濾波可以一般化為下面的積分公式：

的求解一般可以采用一系列函數(shù)通過點集加權(quán)求和的近似方法，即：

式中 N為總點數(shù)； xi， wj分別為正交點集和權(quán)重。對式（7）采用Spherical Radial變換： x =ry且 yΤy = 1；xΤx =r2，則式（7）可重寫為

式中 σ（·）為球面區(qū)域= { y ∈ RnyΤy = 1}的面積微元。通過式（9）可以看出，高斯域下的貝葉斯濾波公式通過Spherical-Radial準(zhǔn)則變換，被分解為一個球面積分和一個徑向積分。這兩個積分同樣無法直接求得，此時再利用式（8），球面積分用包含sN個點的球面積分準(zhǔn)則計算[4]：

徑向積分可由包含rN個點的高斯求積準(zhǔn)則求得：

下面以這兩個積分為基礎(chǔ)介紹 Spherical Simplex準(zhǔn)則和Radial準(zhǔn)則。

2.1.1 Spherical Simplex準(zhǔn)則

首先計算球面積分，球面積分 ()Sr的一種高效率的計算方法是：選取一系列包含 n個正則單行頂點的點集 aj=[aj,1, aj,2, … ,aj,n]Τ, j = 1 ,2,… ,n +1作為容積點[14,15]（n為狀態(tài)維數(shù)），采用中心對稱的容積準(zhǔn)則來近似球面面積。其中容積點的每個元素,jia 選取規(guī)則為

利用ja構(gòu)造Spherical Simplex準(zhǔn)則形式如下：

式中為單位球的表面積，，且

計算徑向積分S(r)rn?1e?r2dr ，與式（13）相對應(yīng)，通過與式（11）進(jìn)行匹配可以得到 Nr=1的Radial準(zhǔn)則[12]：

2.1.2 Radial準(zhǔn)則

由 Γ ( n + 1 ) = n Γ ( n)對第 2個等式進(jìn)行化簡，解出，權(quán)重 wr,1=Γ( n/2)/2。

非線性函數(shù)與高斯概率密度的乘積的積分是高斯域 下 貝 葉 斯 濾 波 器 的 核 心 ， 當(dāng)(x ) = e?xΤx，w2(x ) = N (x ; 0, I)時，通過式（7）可得：

由式（13）、式（14）和式（20）可以求得：

把式（13）、式（16）代入式（18），即可得到 Spherical Simplex-Radial準(zhǔn)則，用 ()Qf表示：

權(quán)重為 N (x;μ,P)時，μ，P分別為x的均值與協(xié)方差，對式（19）進(jìn)行線性變換，得到一般形式Spherical Simplex-Radial 準(zhǔn)則：

2.2 SSRCKF算法

基于Spherical Simplex-Radial 準(zhǔn)則的容積卡爾曼濾波采用一系列等權(quán)值的容積點對高斯域下的貝葉斯濾波進(jìn)行非線性逼近。其中容積點的選取規(guī)則為

式中 2m n= ，（n為狀態(tài)空間維數(shù)）；為權(quán)值；,j = 1,2,… ,n +1，由式（12）得到；的第i列。

算法具體實現(xiàn)步驟如下：

a）步驟1：濾波器初始化。

SSRCKF濾波器初始化如下：

循環(huán)，完成以下步驟。

b）步驟2：時間更新。

首先對進(jìn)行 Cholesky分解，有公式，取下三角陣，計算容積點：

狀態(tài)方程傳遞容積點：

估計k時刻狀態(tài)預(yù)測值：

估計k時刻的先驗估計誤差協(xié)方差：

c）步驟3：量測更新。

對先驗估計誤差協(xié)方差k?P進(jìn)行Cholesky分解：

計算容積點：

量測方程傳遞容積點：

估計k時刻量測預(yù)測值：

估計k時刻的量測誤差協(xié)方差：

估計k時刻的一步預(yù)測互相關(guān)誤差協(xié)方差：

計算k時刻濾波增益：

計算k時刻的狀態(tài)估計值：

計算k時刻的后驗估計誤差協(xié)方差：

通過完成上述3個步驟，完成了由k-1時刻到k時刻的后驗估計值與后驗估計誤差協(xié)方差的更新。

3 仿真實驗與分析

在Matlab（R2010b）環(huán)境下對雷達(dá)觀測載入彈道目標(biāo)系統(tǒng)進(jìn)行建模仿真，對提出的基于 Spherical Simplex-Radial 準(zhǔn)則的容積卡爾曼濾波算法進(jìn)行驗證。首先基于狀態(tài)模型，在雷達(dá)站坐標(biāo)系下生成帶有狀態(tài)噪聲的載入彈道目標(biāo)的真實軌跡，用四階龍格-庫塔方法對狀態(tài)方程進(jìn)行離散時，離散步長h=1；然后根據(jù)量測模型生成帶有量測噪聲的測量信息用于濾波。再入目標(biāo)初始狀態(tài)以及彈道系數(shù)常值為

雷達(dá)每秒對再入彈道目標(biāo)進(jìn)行一次測量，狀態(tài)噪聲協(xié)方差Q與量測噪聲協(xié)方差R分別為

通過仿真，得到再入彈道目標(biāo)的運動軌跡，如圖2所示。

圖2 再入彈道目標(biāo)的運動軌跡Fig.2 The Trajectory of the Ballistic Trajectory

給定濾波初值如下：

分別以UKF，CKF和SSRCKF算法對再入彈道目標(biāo)進(jìn)行實時跟蹤，設(shè)蒙特卡洛打靶次數(shù)為500。根據(jù)速度和位置均方根誤差對比 3種算法的性能。位置均方根誤差定義為

式中 N為蒙特卡洛打靶次數(shù)；，分別為第 n次打靶時 k時刻再入彈道目標(biāo)的真實值與濾波值，速度均方根誤差與式（37）方法相同。

仿真結(jié)果如圖3、圖4所示。圖3為彈道目標(biāo)速度均方根誤差，圖4為彈道目標(biāo)位置均方根誤差。從圖3、圖4中，可以看出，SSRCKF算法相比于UKF算法與CKF算法有更高的再入彈道目標(biāo)實時跟蹤精度。

圖3 彈道目標(biāo)速度均方根誤差Fig.3 Velocity Root-mean-square Error of the Ballistic Target

圖4 彈道目標(biāo)位置均方根誤差Fig.4 Root Mean Square Error of Trajectory Target

統(tǒng)計 200～250 s，UKF、CKF、SSRCKF 3 種算法對再入彈道目標(biāo)實時跟蹤的速度與位置的均方根誤差，求出其平均值，如表1所示。由表1可以看出：相同條件下對再入彈道目標(biāo)進(jìn)行實時跟蹤，SSRCKF算法在定位精度上比CKF算法提高了約4.5 m，比UKF算法提高了5 m；在定速精度上比CKF算法提高了約0.6 m/s，比UKF算法提高了0.7 m/s，充分證明了算法的有效性。

表1 200～250 s速度與位置均方根誤差平均值Tab. 1 Mean Square Root Mean Square Error of 200～250 s

4 結(jié) 論

本文首先對再入彈道目標(biāo)進(jìn)行動力學(xué)建模，然后以推導(dǎo)的基于Spherical Simplex-Radial準(zhǔn)則的容積卡爾曼濾波算法對再入彈道目標(biāo)進(jìn)行實時跟蹤，該算法在沒有明顯提高計算量的前提下，有效提高了對再入彈道目標(biāo)的實時跟蹤精度；最后通過仿真驗證，證明所提出的方法較UKF與經(jīng)典的CKF算法，有更好的性能。