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

0 引 言

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

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

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

1 彈道目標跟蹤數(shù)學模型

1.1 再入段動力學模型

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

圖1 雷達觀測再入彈道目標示意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)為目標到地心距離，()Vt為目標再入速度，ρ(t)為大氣密度模型，ρ(t ) = ρ0exp ( ?(R (t) ? Re)H0)，H0為大氣密度標高，H0=6700km，ρ0為海平面的大氣密度，ρ0= 1 .225kg/m3；β( t)為彈道系數(shù)模型，有[13]：

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

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

1.2 量測模型

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

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

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

2.1 Spherical Simplex-Radial準則

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

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

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

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

徑向積分可由包含rN個點的高斯求積準則求得：

下面以這兩個積分為基礎介紹 Spherical Simplex準則和Radial準則。

2.1.1 Spherical Simplex準則

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

利用ja構造Spherical Simplex準則形式如下：

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

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

2.1.2 Radial準則

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

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

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

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

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

2.2 SSRCKF算法

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

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

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

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

SSRCKF濾波器初始化如下：

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

b）步驟2：時間更新。

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

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

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

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

c）步驟3：量測更新。

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

計算容積點：

量測方程傳遞容積點：

估計k時刻量測預測值：

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

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

計算k時刻濾波增益：

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

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

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

3 仿真實驗與分析

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

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

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

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

給定濾波初值如下：

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

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

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

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

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

統(tǒng)計 200～250 s，UKF、CKF、SSRCKF 3 種算法對再入彈道目標實時跟蹤的速度與位置的均方根誤差，求出其平均值，如表1所示。由表1可以看出：相同條件下對再入彈道目標進行實時跟蹤，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 結 論

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