喻晨龍,譚賢四,曲智國,王 紅,謝 非
(1.空軍預警學院防空預警裝備系,武漢 430019; 2.長江勘測規(guī)劃設計研究院,武漢 430019)
當前可能用于作戰(zhàn)的臨近空間高超聲速武器主要有兩類:①高超聲速助推滑翔導彈,典型有HTV-2和AHW;②高超聲速巡航導彈,典型有X-43A、X-51A、“鋯石”、“布拉莫斯-2”[1-3]。面對日益臨近的潛在威脅,如何更好的跟蹤這些目標成為亟待解決的現(xiàn)實問題。
關于臨近空間高超聲速滑翔彈的跟蹤研究主要包括跟蹤模型和濾波算法。當跟蹤目標的運動狀態(tài)是非線性時,需要用到EKF、UKF和CKF等非線性濾波算法[4-6],當跟蹤目標處于強機動狀態(tài)時,常常采用交互多模型(IMM)濾波算法[7-8]。然而,這些都取決于對目標運動狀態(tài)的數(shù)學描述,也就是目標的狀態(tài)模型,好的狀態(tài)模型不僅可以提高跟蹤精度,還能夠簡化跟蹤過程提高運算速度。
跟蹤的狀態(tài)模型的建立有兩種方法:①從飛行受力的角度考慮,結合推力模型、氣動力模型、大氣模型等建模得到目標運動學方程[9-10]。文獻[11]對再入飛行器的氣動力進行了深入分析,提出把氣動力投影到速度坐標系實現(xiàn)解耦,然后利用非線性Kalman濾波器進行狀態(tài)估計。②在彈道分析的基礎上,對目標的某些運動參數(shù)變化規(guī)律合理假設,對隨機過程噪聲加以描述[12-13]。文獻[14]指出臨近空間高超聲速滑翔目標的位置軌跡曲線呈現(xiàn)明顯的周期特性,用正弦自相關隨機過程來描述目標的加速度變化規(guī)律建立了Sine Wave模型。
關于目標的機動模型,Li和Jilkov等[15-18]的一系列文章對此進行了詳細論述,他們指出跟蹤的狀態(tài)模型可以歸為兩大類:①坐標解耦的Markov模型,這包括CA、CV、CS、Singer、Jerk等;②坐標耦合的空間模型,這包括2D水平運動模型和3D運動模型,也就是CT模型。當跟蹤模型的狀態(tài)量參數(shù)屬于同一個位置域或角度域時,模型是線性的,然而,空間上的CT模型通常轉向率未知,需要進行狀態(tài)擴維,使得狀態(tài)模型非線性,大大降低了跟蹤精度。關于目標的三維跟蹤方法,張翔宇等在文獻[19]中提出將量測數(shù)據轉換到經緯高坐標系,在經緯平面采用基于點跡歸并和凝聚的線性Kalman濾波器,在高度方向采用Singer+CA+CT組成的IMM濾波器并行跟蹤,然后合并數(shù)據輸出三維跟蹤濾波結果。
在本文中,我們通過角度域的參數(shù)描述目標的運動狀態(tài),建立狀態(tài)模型,它們具有比位置參數(shù)更明顯的變化規(guī)律,能更準確地近似目標的真實運動。在角度域中,由于代表縱向機動的航跡傾角存在明顯周期性,代表橫向機動的航向角呈近似線性,提出在縱向通道和橫向通道分別建立角度域參數(shù)的Markov模型。雖然航跡傾角和航向角的規(guī)律顯著,但是在雷達觀測坐標系下是不可見的。在縱向觀測平面,傾角變化率與航跡傾角變化率近似,可以用傾角替代航跡傾角建立Markov模型,進而估算出角度變化率,把它輸入到隨后的空間CT模型中,使得其狀態(tài)模型線性。在橫向觀測平面,方位角是可觀測角,基于方位角建立Markov模型,進而估算出每一時刻目標的橫向偏移。通過量測轉換,雙通道并行跟蹤和輸出合并三個步驟,實現(xiàn)了目標位置狀態(tài)的三維跟蹤。最后設置了一個仿真場景,分析了系統(tǒng)參數(shù)的取值范圍和算法的魯棒性,與現(xiàn)有的臨近空間高超聲速目標跟蹤算法相比,該算法具有較高的跟蹤精度。
坐標平面如圖2所示,圖2(a)、(b)分別表示縱向觀測平面(Longitudinal observation plane, LOP)和橫向觀測平面(Horizontal observation plane,HOP)。LOP與ENU坐標系的OSySzS平面重合,HOP與ENU坐標系的OSxSyS平面重合。Tra表示目標彈道,(x,y,z)表示目標位置,(r0,θ,?)代表測量得到的距離、俯仰角、方位角,γ表示航跡傾角,β表示航程角,φ表示傾角。P′表示質心的投影,R為地球半徑,r為地心距。
圖1 ECEF和ENU坐標系Fig.1 ECEF and ENU Coordinate system
圖2 觀測平面Fig.2 Observation planes
臨近空間高超聲速滑翔彈飛行全程主要受到地球重力,氣動力和離心力的作用,當這三個力在縱向上的分量的合力為0時,構成平衡滑翔條件,軌跡為無震蕩的平坦光滑曲線,運動參數(shù)(速度、高度、航跡傾角)之間存在特定的解析關系,當合力不為0時,不滿足平衡滑翔條件,軌跡為跳躍震蕩曲線,運動參數(shù)之間不存在近似的解析關系,只能依據目標的三自由度彈道方程,通過數(shù)值仿真分析運動參數(shù)的變化規(guī)律[20]。
目標位置域的參數(shù)包括高度和速度,目標角度域的參數(shù)包括航跡傾角和航向角。高度衰減震蕩,速度成階梯式下降,航跡傾角為小量且周期性變化,航向角臺階式增大。我們還可以得到如下結論:
1)軌跡可以分解到縱向和橫向平面??刂屏抗ソ菦Q定著飛行器縱向上的機動,控制量傾側角決定著飛行器橫向上的機動,在縱向上有機動強度很大的跳躍運動,在橫向上則是機動強度較弱的偏移運動。機動強度越大跟蹤模型的數(shù)學描述越復雜,跟蹤的關鍵在縱向上。
2)角度域的參數(shù)具有較為確定的顯著規(guī)律。航跡傾角反映著目標的縱向機動,始終保持著周期性變化,航向角反映著目標的橫向機動,隨著傾側角的增大偏移越來越大。高度和速度不能直接反映目標橫縱向機動運動。
3)當傾側角σ=0時,目標軌跡橫向上不發(fā)生偏移,航向角的變化率維持恒零,僅在縱向上有跳躍機動。當傾側角σ不為零時,目標軌跡不僅在縱向上有跳躍機動,在橫向上也有偏轉機動。目標處于跳躍滑翔狀態(tài)時,在縱向上的跳躍機動是始終存在的,在橫向上的偏移轉彎則受到傾側角控制,對應到角度域的參數(shù)上來說,航跡傾角的變化率是始終存在的,航向角的變化率受傾側角控制。
因此,在跟蹤高超聲速目標時我們可以從角度域參數(shù)著手,把目標軌跡投影到縱向觀測平面和橫向觀測平面,基于航跡傾角和航向角分別建立縱向通道和橫向通道的跟蹤模型,然后并行跟蹤濾波。但是,航跡傾角和航向角在ENU坐標系下是不可見的,不能直接用于建模,要在LOP和HOP平面內尋找替代量,參與模型的建立。
如圖2(a),在LOP中有如下關系:
(1)
對式(1)的第1個式子求導:
(2)
從式(2)可知航跡傾角變化率、傾角變化率、航程角變化率之間存在加法關系,需要明確這三個角中的哪一個可以作為替代量,用來建立縱向通道的跟蹤模型。
圖3 縱向觀測平面的角度變化率Fig.3 Change rate of angles in longitudinal observation plane
(3)
(4)
因此,在HOP平面,用方位角代替航向角,把方位角變化率描述為正弦自相關隨機過程,建立橫向通道的機動模型。
在上文中,我們找到了角度域中的傾角和方位角參數(shù),代替原來的航跡傾角和航向角,分別在縱向通道和橫向通道建立跟蹤模型。在建立角度域模型時,均把角度變化率描述為零均值正弦自相關函數(shù),構造了角速度正弦波(Angular velocity sine wave, AVSW)模型。在縱向通道,為了得到目標的位置參數(shù),在AVSW后端嫁接了一個CT模型,把角度變化率輸入到CT模型中。在橫向通道,直接輸出方位角估計。系統(tǒng)通過量測轉換,雙通道并行濾波和輸出合并,實現(xiàn)了目標的三維跟蹤,其結構如圖4所示。
圖4(a)描述了系統(tǒng)的總體結構,在跟蹤過程中,首先把量測數(shù)據由球坐標轉換到直角坐標,然后分類輸入到縱向通道和橫向通道,最后通過一個數(shù)據合并模塊后,輸出目標的三維狀態(tài)估計。圖4(b)描述了縱向通道的AVSW-CT模型結構,輸入是經過量測轉換后的LOP平面數(shù)據,該模型包括一個角度子濾波器和一個位置子濾波器,其中角度子濾波器的狀態(tài)模型為傾角變化率的Markov模型,位置子濾波器的狀態(tài)模型為CT結構的空間模型。AVSW模型作用于角度域,CT模型作用于位置域。圖4(c)描述了橫向通道的AVSW模型結構,輸入為目標的方位角觀測值,該模型包含一個角度濾波器,其狀態(tài)模型為方位角變化率的Markov模型。圖4(d)描述了雙通道的合并輸出方法,縱向通道輸出目標縱向的位置和速度,橫向通道輸出橫向的方位角,先得到目標的橫向偏移,然后得到目標的三維位置估計。其詳細的論述見下文。
不論控制量如何變化,縱向上的跳躍機動始終存在,AVSW-CT模型用來描述目標在縱向上的機動運動,它由一個角度子濾波器和一個位置子濾波器組成。
1)角度子濾波器的狀態(tài)模型
圖4 系統(tǒng)架構Fig.4 System frame
角度濾波器的狀態(tài)模型也就是AVSW模型,其推到過程如下。假設目標的傾角變化率為一零均值的正弦隨機信號,其自相關函數(shù)為:
(5)
這是一有色噪聲信號,其功率譜密度為:
(6)
因此傾角變化率的二階Markov模型為:
(7)
(8)
那么,傾角變化率的離散時間狀態(tài)方程為:
Xω(k+1)=AωXω(k)+νωc(k)
(9)
其中:
(10)
過程噪聲νωc具有的協(xié)方差為:
(11)
其中:
(12)
2)角度子濾波器的量測模型
縱向通道的跟蹤在LOP平面進行,首先需要把以球坐標形式表示的ENU坐標系下的觀測數(shù)據投影到LOP平面。這里采用無偏轉換[14],實現(xiàn)球-直坐標轉換。
然而,得到的位置屬性的量測數(shù)據不能直接輸入到角度子濾波器,需要根據式(1)的第二個式子計算,把它轉換為角度格式。
那么,角度子濾波器的等效量測模型為:
Zω,k=HωXω,k+κω,k
(13)
角度子濾波器輸出角度變化率,輸入到位置子濾波器,使得其狀態(tài)模型線性。
3)位置子濾波器的狀態(tài)模型
位置子濾波器的狀態(tài)模型也就是CT模型,是LOP中角度關系在位置維度的映射。在圖2(a)中有如下關系:
(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)
位置子濾波器的輸出為目標的位置和速度。
當傾側角為零時,目標僅在縱向平面跳躍,當傾側角不為零時,目標除了在縱向上跳躍,還在橫向上偏轉,AVSW模型用來描述目標在橫向上的機動運動,它包含一個角度濾波器。
1)角度濾波器的狀態(tài)模型
角度濾波器的狀態(tài)模型和上文所述的AVSW模型相同。假設方位角變化率為一零均值的正弦隨機信號,其自相關函數(shù)為:
(21)
(22)
X?(k+1)=A?X?(k)+ν?c(k)
(23)
2)角度濾波器的量測模型
橫向通道的輸入為目標的方位角觀測序列,等效量測模型表示為:
Z?,k=H?X?,k+κ?,k
(24)
角度濾波器的輸出為方位角估計值。
(25)
因此k時刻目標在x軸上的橫向的位置估計為:
xk=Xa(1,k)tan(X?(1,k))
(26)
綜上可見,整個系統(tǒng)均可采用線性的kalman濾波器跟蹤,跟蹤算法的步驟如下:
輸入:ENU坐標系下的觀測序列。
輸出:ENU坐標系下的狀態(tài)估計。
步驟1:量測數(shù)據預處理
步驟1.1:將方位角原始測量數(shù)據輸入到橫向通道,直接獲得角度格式的量測序列{z?,k}。
步驟1.2:將觀測值從球坐標系轉換到直角坐標系,提取y軸和z軸數(shù)據,組成位置格式的量測序列{za,k}。根據式(1)的第2個式子,將位置格式的量測序列{za,k}轉換成角度格式的量測序列{zω,k}。
步驟2:線性濾波器并行跟蹤
步驟2.1:橫向通道跟蹤
假設k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為X?,k-1/k-1和P?,k-1/k-1。線性Kalman濾波:
(27)
(28)
(29)
步驟2.2:縱向通道跟蹤
步驟2.2.1:角度子濾波器跟蹤:
假設在k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為Xω,k-1/k-1和Pω,k-1/k-1。線性Kalman濾波:
(30)
(31)
(32)
步驟2.2.2:驅動位置子濾波器跟蹤:
提取角度子濾波器輸出狀態(tài)的第二個分量ω(k-1)=Xω(2,k-1)。利用角度變化率分別生成狀態(tài)轉移矩陣Aa,k-1和過程噪聲協(xié)方差矩陣Qa,k-1。
假設在k-1時刻的狀態(tài)向量和協(xié)方差矩陣分別為Xa,k-1/k-1和Pa,k-1/k-1。線性Kalman濾波:
(33)
(34)
(35)
步驟3:合并輸出。
根據式(26)生成k時刻在x軸上的位置分量,合并狀態(tài)量,輸出目標的位置和速度估計。
參考美國官方公布的HTV-2的基本試驗情況,本文假設如下仿真場景:升力式滑翔目標由火箭助推至80 km高空后釋放,分離時目標的位置為E110°N0°,再入攻角為11.6°,傾側角為20°,速度為20 Ma,航跡傾角和航向角為0°。假設某型雷達部署在E112°N10°,其參數(shù)如表1所示。
表1 雷達參數(shù)Table 1 Radar parameters
由于地球曲率影響,雷達在ENU坐標系下只觀測到了目標的部分軌跡,如圖5所示。
圖5 ENU坐標系下的目標軌跡Fig.5 Trajectory in ENU coordinate system
從圖5可看出,目標在(-975 km, -225 km)處穿越地平面進入雷達視線范圍,然后在縱向上跳躍運動,同時在橫向上發(fā)生偏移,最后在(770 km,-160 km)處進入地平面以下,離開雷達視線范圍。整個過程中目標的飛行速度很快,觀測的時間很短,不到400 s。
首先設置縱向通道的模型參數(shù)α0。假設AVSW-CT模型中的角度子濾波器的角度分布的誤差標準差為10-8,位置子濾波器的位置分布的誤差標準差為103,參數(shù)α0分別為10-1、10-2、 10-3、 10-4、 10-5,進行100次蒙特卡洛仿真,得到結果如圖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模型具有較好的跟蹤精度,當α0≤10-2時,角速度和位置的估計更精確。
然后設置橫向通道的模型參數(shù)α1。假設AVSW模型中的角度濾波器的過程噪聲誤差標準差為10-4,參數(shù)α1分別為10-1、10-2、 10-3、 10-4、 10-5,進行100次蒙特卡洛仿真,得到結果如圖8所示。
圖8 橫向通道角度濾波器輸出誤差(α1)Fig.8 Output error of horizontal channel angle filter (α1)
圖8為角度濾波器的輸出誤差。從圖8中對比可知,AVSW模型具有較好的跟蹤精度,當α1≤10-2時,方位角估計更精確。
1)魯棒性分析
通過仿真試驗分析模型過程噪聲對跟蹤精度的影響,首先是縱向通道的角度子濾波器。假設AVSW-CT模型中的參數(shù)α0為10-3,角度子濾波器的過程噪聲標準差σω分別為10-6、 10-7、 10-8、 10-9、 10-10,進行100次蒙特卡洛仿真,得到結果如圖9所示。
圖9 縱向通道角度子濾波器輸出誤差(σω)Fig.9 Output error of longitudinal channel angle sub-filter (σω)
圖9為角度子濾波器的輸出誤差。從圖中的對比可知,當σω≤10-8時,跟蹤精度更高。
然后分析縱向通道的位置子濾波器。假設AVSW-CT模型中的參數(shù)α0為10-3,角度子濾波器的過程噪聲標準差σω取10-8,位置子濾波器的過程噪聲標準差σa分別為10、 102、 103、 104、 105進行100次蒙特卡洛仿真,得到結果如圖10所示。
圖10 縱向通道位置子濾波器輸出誤差(σa)Fig.10 Output error of longitudinal channel position sub-filter (σa)
圖10為位置子濾波器的輸出誤差。從圖中的對比可知,噪聲對跟蹤精度的影響較明顯,當σa∈[1001000]時,跟蹤精度最佳。
最后分析橫向通道的角度濾波器。假設AVSW模型中的參數(shù)α1為10-3,角度濾波器的過程噪聲標準差σ?分別為10-1, 10-2, 10-3, 10-4, 10-5,進行100次蒙特卡洛仿真,得到結果如圖11所示。
圖11 橫向通道角度濾波器輸出誤差(σ?)Fig.11 Output error of horizontal channel angle filter (σ?)
圖11為角度濾波器的輸出誤差。從圖中的對比可知,噪聲對跟蹤精度的影響較明顯,當σ?∈[10-410-3]時,跟蹤精度最佳。
2)跟蹤精度對比
通過與其他算法比較觀察文中算法的適用性。假設算法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次蒙特卡洛仿真,得到結果如圖12所示。
圖12 系統(tǒng)輸出誤差Fig.12 Output error of system
圖12為系統(tǒng)的輸出誤差,可看出在文中所設參數(shù)場景下,文中的雙通道并行跟蹤濾波算法具有更好的跟蹤效果。
本文研究了臨近空間高超聲速滑翔彈的運動特性,分析了運動參數(shù)的變化規(guī)律,從角度域參數(shù)入手分別建立了目標橫向和縱向的機動模型,提出了一種雙通道并行跟蹤算法,實現(xiàn)了三維跟蹤。仿真對比表明,本文所提的算法具有較強的穩(wěn)定性,具有較高的跟蹤精度。