臨近空間高超聲速滑翔彈雙通道跟蹤算法

喻晨龍，譚賢四，曲智國，王 紅，謝 非

(1.空軍預(yù)警學(xué)院防空預(yù)警裝備系，武漢 430019； 2.長江勘測規(guī)劃設(shè)計研究院，武漢 430019)

0 引 言

當(dāng)前可能用于作戰(zhàn)的臨近空間高超聲速武器主要有兩類：①高超聲速助推滑翔導(dǎo)彈，典型有HTV-2和AHW；②高超聲速巡航導(dǎo)彈，典型有X-43A、X-51A、“鋯石”、“布拉莫斯-2”[1-3]。面對日益臨近的潛在威脅，如何更好的跟蹤這些目標(biāo)成為亟待解決的現(xiàn)實問題。

關(guān)于臨近空間高超聲速滑翔彈的跟蹤研究主要包括跟蹤模型和濾波算法。當(dāng)跟蹤目標(biāo)的運動狀態(tài)是非線性時，需要用到EKF、UKF和CKF等非線性濾波算法[4-6]，當(dāng)跟蹤目標(biāo)處于強機動狀態(tài)時，常常采用交互多模型(IMM)濾波算法[7-8]。然而，這些都取決于對目標(biāo)運動狀態(tài)的數(shù)學(xué)描述，也就是目標(biāo)的狀態(tài)模型，好的狀態(tài)模型不僅可以提高跟蹤精度，還能夠簡化跟蹤過程提高運算速度。

跟蹤的狀態(tài)模型的建立有兩種方法:①從飛行受力的角度考慮，結(jié)合推力模型、氣動力模型、大氣模型等建模得到目標(biāo)運動學(xué)方程[9-10]。文獻[11]對再入飛行器的氣動力進行了深入分析，提出把氣動力投影到速度坐標(biāo)系實現(xiàn)解耦，然后利用非線性Kalman濾波器進行狀態(tài)估計。②在彈道分析的基礎(chǔ)上，對目標(biāo)的某些運動參數(shù)變化規(guī)律合理假設(shè)，對隨機過程噪聲加以描述[12-13]。文獻[14]指出臨近空間高超聲速滑翔目標(biāo)的位置軌跡曲線呈現(xiàn)明顯的周期特性,用正弦自相關(guān)隨機過程來描述目標(biāo)的加速度變化規(guī)律建立了Sine Wave模型。

關(guān)于目標(biāo)的機動模型，Li和Jilkov等[15-18]的一系列文章對此進行了詳細論述，他們指出跟蹤的狀態(tài)模型可以歸為兩大類:①坐標(biāo)解耦的Markov模型，這包括CA、CV、CS、Singer、Jerk等;②坐標(biāo)耦合的空間模型，這包括2D水平運動模型和3D運動模型，也就是CT模型。當(dāng)跟蹤模型的狀態(tài)量參數(shù)屬于同一個位置域或角度域時，模型是線性的，然而，空間上的CT模型通常轉(zhuǎn)向率未知，需要進行狀態(tài)擴維，使得狀態(tài)模型非線性，大大降低了跟蹤精度。關(guān)于目標(biāo)的三維跟蹤方法，張翔宇等在文獻[19]中提出將量測數(shù)據(jù)轉(zhuǎn)換到經(jīng)緯高坐標(biāo)系，在經(jīng)緯平面采用基于點跡歸并和凝聚的線性Kalman濾波器，在高度方向采用Singer+CA+CT組成的IMM濾波器并行跟蹤，然后合并數(shù)據(jù)輸出三維跟蹤濾波結(jié)果。

在本文中，我們通過角度域的參數(shù)描述目標(biāo)的運動狀態(tài)，建立狀態(tài)模型，它們具有比位置參數(shù)更明顯的變化規(guī)律，能更準(zhǔn)確地近似目標(biāo)的真實運動。在角度域中，由于代表縱向機動的航跡傾角存在明顯周期性，代表橫向機動的航向角呈近似線性，提出在縱向通道和橫向通道分別建立角度域參數(shù)的Markov模型。雖然航跡傾角和航向角的規(guī)律顯著，但是在雷達觀測坐標(biāo)系下是不可見的。在縱向觀測平面，傾角變化率與航跡傾角變化率近似，可以用傾角替代航跡傾角建立Markov模型，進而估算出角度變化率，把它輸入到隨后的空間CT模型中，使得其狀態(tài)模型線性。在橫向觀測平面，方位角是可觀測角，基于方位角建立Markov模型，進而估算出每一時刻目標(biāo)的橫向偏移。通過量測轉(zhuǎn)換，雙通道并行跟蹤和輸出合并三個步驟，實現(xiàn)了目標(biāo)位置狀態(tài)的三維跟蹤。最后設(shè)置了一個仿真場景，分析了系統(tǒng)參數(shù)的取值范圍和算法的魯棒性，與現(xiàn)有的臨近空間高超聲速目標(biāo)跟蹤算法相比，該算法具有較高的跟蹤精度。

1 目標(biāo)運動特性

1.1 參考坐標(biāo)系和坐標(biāo)平面

坐標(biāo)平面如圖2所示，圖2(a)、(b)分別表示縱向觀測平面(Longitudinal observation plane, LOP)和橫向觀測平面(Horizontal observation plane,HOP)。LOP與ENU坐標(biāo)系的OSySzS平面重合，HOP與ENU坐標(biāo)系的OSxSyS平面重合。Tra表示目標(biāo)彈道，(x,y,z)表示目標(biāo)位置，(r0,θ,?)代表測量得到的距離、俯仰角、方位角，γ表示航跡傾角，β表示航程角，φ表示傾角。P′表示質(zhì)心的投影，R為地球半徑，r為地心距。

圖1 ECEF和ENU坐標(biāo)系Fig.1 ECEF and ENU Coordinate system

圖2 觀測平面Fig.2 Observation planes

1.2 運動軌跡特性分析

臨近空間高超聲速滑翔彈飛行全程主要受到地球重力，氣動力和離心力的作用，當(dāng)這三個力在縱向上的分量的合力為0時，構(gòu)成平衡滑翔條件，軌跡為無震蕩的平坦光滑曲線，運動參數(shù)(速度、高度、航跡傾角)之間存在特定的解析關(guān)系，當(dāng)合力不為0時，不滿足平衡滑翔條件，軌跡為跳躍震蕩曲線，運動參數(shù)之間不存在近似的解析關(guān)系，只能依據(jù)目標(biāo)的三自由度彈道方程，通過數(shù)值仿真分析運動參數(shù)的變化規(guī)律[20]。

目標(biāo)位置域的參數(shù)包括高度和速度，目標(biāo)角度域的參數(shù)包括航跡傾角和航向角。高度衰減震蕩，速度成階梯式下降，航跡傾角為小量且周期性變化，航向角臺階式增大。我們還可以得到如下結(jié)論：

1)軌跡可以分解到縱向和橫向平面。控制量攻角決定著飛行器縱向上的機動，控制量傾側(cè)角決定著飛行器橫向上的機動，在縱向上有機動強度很大的跳躍運動，在橫向上則是機動強度較弱的偏移運動。機動強度越大跟蹤模型的數(shù)學(xué)描述越復(fù)雜，跟蹤的關(guān)鍵在縱向上。

2)角度域的參數(shù)具有較為確定的顯著規(guī)律。航跡傾角反映著目標(biāo)的縱向機動，始終保持著周期性變化，航向角反映著目標(biāo)的橫向機動，隨著傾側(cè)角的增大偏移越來越大。高度和速度不能直接反映目標(biāo)橫縱向機動運動。

3)當(dāng)傾側(cè)角σ=0時，目標(biāo)軌跡橫向上不發(fā)生偏移，航向角的變化率維持恒零，僅在縱向上有跳躍機動。當(dāng)傾側(cè)角σ不為零時，目標(biāo)軌跡不僅在縱向上有跳躍機動，在橫向上也有偏轉(zhuǎn)機動。目標(biāo)處于跳躍滑翔狀態(tài)時，在縱向上的跳躍機動是始終存在的，在橫向上的偏移轉(zhuǎn)彎則受到傾側(cè)角控制，對應(yīng)到角度域的參數(shù)上來說，航跡傾角的變化率是始終存在的，航向角的變化率受傾側(cè)角控制。

因此，在跟蹤高超聲速目標(biāo)時我們可以從角度域參數(shù)著手，把目標(biāo)軌跡投影到縱向觀測平面和橫向觀測平面，基于航跡傾角和航向角分別建立縱向通道和橫向通道的跟蹤模型，然后并行跟蹤濾波。但是，航跡傾角和航向角在ENU坐標(biāo)系下是不可見的，不能直接用于建模，要在LOP和HOP平面內(nèi)尋找替代量，參與模型的建立。

如圖2(a)，在LOP中有如下關(guān)系：

(1)

對式(1)的第1個式子求導(dǎo)：

(2)

從式(2)可知航跡傾角變化率、傾角變化率、航程角變化率之間存在加法關(guān)系，需要明確這三個角中的哪一個可以作為替代量，用來建立縱向通道的跟蹤模型。

圖3 縱向觀測平面的角度變化率Fig.3 Change rate of angles in longitudinal observation plane

(3)

(4)

因此，在HOP平面，用方位角代替航向角，把方位角變化率描述為正弦自相關(guān)隨機過程，建立橫向通道的機動模型。

2 跟蹤模型與算法

在上文中，我們找到了角度域中的傾角和方位角參數(shù)，代替原來的航跡傾角和航向角，分別在縱向通道和橫向通道建立跟蹤模型。在建立角度域模型時，均把角度變化率描述為零均值正弦自相關(guān)函數(shù)，構(gòu)造了角速度正弦波(Angular velocity sine wave, AVSW)模型。在縱向通道，為了得到目標(biāo)的位置參數(shù)，在AVSW后端嫁接了一個CT模型，把角度變化率輸入到CT模型中。在橫向通道，直接輸出方位角估計。系統(tǒng)通過量測轉(zhuǎn)換，雙通道并行濾波和輸出合并，實現(xiàn)了目標(biāo)的三維跟蹤，其結(jié)構(gòu)如圖4所示。

圖4(a)描述了系統(tǒng)的總體結(jié)構(gòu)，在跟蹤過程中，首先把量測數(shù)據(jù)由球坐標(biāo)轉(zhuǎn)換到直角坐標(biāo)，然后分類輸入到縱向通道和橫向通道，最后通過一個數(shù)據(jù)合并模塊后，輸出目標(biāo)的三維狀態(tài)估計。圖4(b)描述了縱向通道的AVSW-CT模型結(jié)構(gòu)，輸入是經(jīng)過量測轉(zhuǎn)換后的LOP平面數(shù)據(jù)，該模型包括一個角度子濾波器和一個位置子濾波器，其中角度子濾波器的狀態(tài)模型為傾角變化率的Markov模型，位置子濾波器的狀態(tài)模型為CT結(jié)構(gòu)的空間模型。AVSW模型作用于角度域，CT模型作用于位置域。圖4(c)描述了橫向通道的AVSW模型結(jié)構(gòu)，輸入為目標(biāo)的方位角觀測值，該模型包含一個角度濾波器，其狀態(tài)模型為方位角變化率的Markov模型。圖4(d)描述了雙通道的合并輸出方法，縱向通道輸出目標(biāo)縱向的位置和速度，橫向通道輸出橫向的方位角，先得到目標(biāo)的橫向偏移，然后得到目標(biāo)的三維位置估計。其詳細的論述見下文。

2.1 縱向通道AVSW-CT模型

不論控制量如何變化，縱向上的跳躍機動始終存在，AVSW-CT模型用來描述目標(biāo)在縱向上的機動運動，它由一個角度子濾波器和一個位置子濾波器組成。

1)角度子濾波器的狀態(tài)模型

圖4 系統(tǒng)架構(gòu)Fig.4 System frame

角度濾波器的狀態(tài)模型也就是AVSW模型，其推到過程如下。假設(shè)目標(biāo)的傾角變化率為一零均值的正弦隨機信號，其自相關(guān)函數(shù)為：

(5)

這是一有色噪聲信號，其功率譜密度為：

(6)

因此傾角變化率的二階Markov模型為：

(7)

(8)

那么，傾角變化率的離散時間狀態(tài)方程為：

Xω(k+1)=AωXω(k)+νωc(k)

(9)

其中：

(10)

過程噪聲νωc具有的協(xié)方差為：

(11)

其中：

(12)

2)角度子濾波器的量測模型

縱向通道的跟蹤在LOP平面進行，首先需要把以球坐標(biāo)形式表示的ENU坐標(biāo)系下的觀測數(shù)據(jù)投影到LOP平面。這里采用無偏轉(zhuǎn)換[14]，實現(xiàn)球-直坐標(biāo)轉(zhuǎn)換。

然而，得到的位置屬性的量測數(shù)據(jù)不能直接輸入到角度子濾波器，需要根據(jù)式(1)的第二個式子計算，把它轉(zhuǎn)換為角度格式。

那么，角度子濾波器的等效量測模型為：

Zω,k=HωXω,k+κω,k

(13)

角度子濾波器輸出角度變化率，輸入到位置子濾波器，使得其狀態(tài)模型線性。

3)位置子濾波器的狀態(tài)模型

位置子濾波器的狀態(tài)模型也就是CT模型，是LOP中角度關(guān)系在位置維度的映射。在圖2(a)中有如下關(guān)系：

(14)

位置的連續(xù)時間狀態(tài)方程為：

(15)

Xa(k+1)=AaXa(k)+νac(k)

(16)

式中：

(17)

過程噪聲νac具有的協(xié)方差為：

(18)

其中：

(19)

4)位置子濾波器的量測模型

位置子濾波器的等效量測模型為：

Za,k=Ha,kXa,k+κa,k

(20)

位置子濾波器的輸出為目標(biāo)的位置和速度。

2.2 橫向通道AVSW模型

當(dāng)傾側(cè)角為零時，目標(biāo)僅在縱向平面跳躍，當(dāng)傾側(cè)角不為零時，目標(biāo)除了在縱向上跳躍，還在橫向上偏轉(zhuǎn)，AVSW模型用來描述目標(biāo)在橫向上的機動運動，它包含一個角度濾波器。

1)角度濾波器的狀態(tài)模型

角度濾波器的狀態(tài)模型和上文所述的AVSW模型相同。假設(shè)方位角變化率為一零均值的正弦隨機信號，其自相關(guān)函數(shù)為：

(21)

(22)

X?(k+1)=A?X?(k)+ν?c(k)

(23)

2)角度濾波器的量測模型

橫向通道的輸入為目標(biāo)的方位角觀測序列，等效量測模型表示為：

Z?,k=H?X?,k+κ?,k

(24)

角度濾波器的輸出為方位角估計值。

2.3 輸出合并

(25)

因此k時刻目標(biāo)在x軸上的橫向的位置估計為：

xk=Xa(1,k)tan(X?(1,k))

(26)

2.4 跟蹤算法流程

綜上可見，整個系統(tǒng)均可采用線性的kalman濾波器跟蹤，跟蹤算法的步驟如下：

輸入：ENU坐標(biāo)系下的觀測序列。

輸出：ENU坐標(biāo)系下的狀態(tài)估計。

步驟1：量測數(shù)據(jù)預(yù)處理

步驟1.1：將方位角原始測量數(shù)據(jù)輸入到橫向通道，直接獲得角度格式的量測序列{z?,k}。

步驟1.2：將觀測值從球坐標(biāo)系轉(zhuǎn)換到直角坐標(biāo)系，提取y軸和z軸數(shù)據(jù)，組成位置格式的量測序列{za,k}。根據(jù)式(1)的第2個式子，將位置格式的量測序列{za,k}轉(zhuǎn)換成角度格式的量測序列{zω,k}。

步驟2：線性濾波器并行跟蹤

步驟2.1：橫向通道跟蹤

假設(shè)k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為X?,k-1/k-1和P?,k-1/k-1。線性Kalman濾波：

(27)

(28)

(29)

步驟2.2：縱向通道跟蹤

步驟2.2.1：角度子濾波器跟蹤：

假設(shè)在k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為Xω,k-1/k-1和Pω,k-1/k-1。線性Kalman濾波：

(30)

(31)

(32)

步驟2.2.2：驅(qū)動位置子濾波器跟蹤：

提取角度子濾波器輸出狀態(tài)的第二個分量ω(k-1)=Xω(2,k-1)。利用角度變化率分別生成狀態(tài)轉(zhuǎn)移矩陣Aa,k-1和過程噪聲協(xié)方差矩陣Qa,k-1。

假設(shè)在k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為Xa,k-1/k-1和Pa,k-1/k-1。線性Kalman濾波：

(33)

(34)

(35)

步驟3：合并輸出。

根據(jù)式(26)生成k時刻在x軸上的位置分量，合并狀態(tài)量，輸出目標(biāo)的位置和速度估計。

3 仿真校驗

3.1 仿真場景假設(shè)

參考美國官方公布的HTV-2的基本試驗情況，本文假設(shè)如下仿真場景：升力式滑翔目標(biāo)由火箭助推至80 km高空后釋放，分離時目標(biāo)的位置為E110°N0°，再入攻角為11.6°，傾側(cè)角為20°，速度為20 Ma，航跡傾角和航向角為0°。假設(shè)某型雷達部署在E112°N10°，其參數(shù)如表1所示。

表1 雷達參數(shù)Table 1 Radar parameters

由于地球曲率影響，雷達在ENU坐標(biāo)系下只觀測到了目標(biāo)的部分軌跡，如圖5所示。

圖5 ENU坐標(biāo)系下的目標(biāo)軌跡Fig.5 Trajectory in ENU coordinate system

從圖5可看出，目標(biāo)在(-975 km, -225 km)處穿越地平面進入雷達視線范圍，然后在縱向上跳躍運動，同時在橫向上發(fā)生偏移，最后在(770 km，-160 km)處進入地平面以下，離開雷達視線范圍。整個過程中目標(biāo)的飛行速度很快，觀測的時間很短，不到400 s。

3.2 模型參數(shù)設(shè)置

首先設(shè)置縱向通道的模型參數(shù)α0。假設(shè)AVSW-CT模型中的角度子濾波器的角度分布的誤差標(biāo)準(zhǔn)差為10-8，位置子濾波器的位置分布的誤差標(biāo)準(zhǔn)差為103，參數(shù)α0分別為10-1、10-2、 10-3、 10-4、 10-5，進行100次蒙特卡洛仿真，得到結(jié)果如圖6和圖7所示。

圖6 縱向通道角度子濾波器輸出誤差(α0)Fig.6 Output error of longitudinal channel angle sub-filter (α0)

圖7 縱向通道位置子濾波器輸出誤差(α0)Fig.7 Output error of longitudinal channel position sub-filter (α0)

圖6為角度子濾波器的輸出誤差，圖7為位置子濾波器的輸出誤差。從圖6和圖7中對比可知，AVSW-CT模型具有較好的跟蹤精度，當(dāng)α0≤10-2時，角速度和位置的估計更精確。

然后設(shè)置橫向通道的模型參數(shù)α1。假設(shè)AVSW模型中的角度濾波器的過程噪聲誤差標(biāo)準(zhǔn)差為10-4，參數(shù)α1分別為10-1、10-2、 10-3、 10-4、 10-5，進行100次蒙特卡洛仿真，得到結(jié)果如圖8所示。

圖8 橫向通道角度濾波器輸出誤差(α1)Fig.8 Output error of horizontal channel angle filter (α1)

圖8為角度濾波器的輸出誤差。從圖8中對比可知，AVSW模型具有較好的跟蹤精度，當(dāng)α1≤10-2時，方位角估計更精確。

3.3 仿真試驗結(jié)果與分析

1)魯棒性分析

通過仿真試驗分析模型過程噪聲對跟蹤精度的影響，首先是縱向通道的角度子濾波器。假設(shè)AVSW-CT模型中的參數(shù)α0為10-3，角度子濾波器的過程噪聲標(biāo)準(zhǔn)差σω分別為10-6、 10-7、 10-8、 10-9、 10-10，進行100次蒙特卡洛仿真，得到結(jié)果如圖9所示。

圖9 縱向通道角度子濾波器輸出誤差(σω)Fig.9 Output error of longitudinal channel angle sub-filter (σω)

圖9為角度子濾波器的輸出誤差。從圖中的對比可知，當(dāng)σω≤10-8時，跟蹤精度更高。

然后分析縱向通道的位置子濾波器。假設(shè)AVSW-CT模型中的參數(shù)α0為10-3，角度子濾波器的過程噪聲標(biāo)準(zhǔn)差σω取10-8，位置子濾波器的過程噪聲標(biāo)準(zhǔn)差σa分別為10、 102、 103、 104、 105進行100次蒙特卡洛仿真，得到結(jié)果如圖10所示。

圖10 縱向通道位置子濾波器輸出誤差(σa)Fig.10 Output error of longitudinal channel position sub-filter (σa)

圖10為位置子濾波器的輸出誤差。從圖中的對比可知，噪聲對跟蹤精度的影響較明顯，當(dāng)σa∈[1001000]時，跟蹤精度最佳。

最后分析橫向通道的角度濾波器。假設(shè)AVSW模型中的參數(shù)α1為10-3，角度濾波器的過程噪聲標(biāo)準(zhǔn)差σ?分別為10-1, 10-2, 10-3, 10-4, 10-5，進行100次蒙特卡洛仿真，得到結(jié)果如圖11所示。

圖11 橫向通道角度濾波器輸出誤差(σ?)Fig.11 Output error of horizontal channel angle filter (σ?)

圖11為角度濾波器的輸出誤差。從圖中的對比可知，噪聲對跟蹤精度的影響較明顯，當(dāng)σ?∈[10-410-3]時，跟蹤精度最佳。

2)跟蹤精度對比

通過與其他算法比較觀察文中算法的適用性。假設(shè)算法1中采用Singer模型，σa=0.1，α=1/20。算法2中采用文獻[14]的Sine Wave模型，σa=0.1，α=1/18。算法3為本文算法，σω=10-8,σa=103,σ?=10-4,α0=α1=10-3。算法4采用IMM(Singer + CA +CV)。進行100次蒙特卡洛仿真，得到結(jié)果如圖12所示。

圖12 系統(tǒng)輸出誤差Fig.12 Output error of system

圖12為系統(tǒng)的輸出誤差，可看出在文中所設(shè)參數(shù)場景下，文中的雙通道并行跟蹤濾波算法具有更好的跟蹤效果。

4 結(jié) 論

本文研究了臨近空間高超聲速滑翔彈的運動特性，分析了運動參數(shù)的變化規(guī)律，從角度域參數(shù)入手分別建立了目標(biāo)橫向和縱向的機動模型，提出了一種雙通道并行跟蹤算法，實現(xiàn)了三維跟蹤。仿真對比表明，本文所提的算法具有較強的穩(wěn)定性，具有較高的跟蹤精度。