基于軌跡跟蹤預測的反HTV2 中制導研究*

孫雪琪，張銳

（1.北京電子工程總體研究所，北京 100854；2.中國航天科工集團第二研究院 研究生院，北京 100854）

0 引言

臨近空間高超聲速滑翔飛行器（hypersonic technology vehicle 2）HTV2 能夠在臨近空間實現(xiàn)高速高機動無動力滑翔，實施遠距離快速打擊。面對這種高速、高機動性目標，攔截彈通常采用復合制導攔截策略。攔截彈需要對目標機動進行穩(wěn)定跟蹤預測，以便攔截彈以最佳的相對幾何關系進入末制導［1］。目前國內(nèi)外對于跟蹤臨近空間高超聲速滑翔飛行器的研究有很多，主要分為基于動力學方程和運動學方程2 類。文獻［2］基于目標動力學方程，引入了2 個氣動力參數(shù)，并使用擴展卡爾曼濾波對該氣動參數(shù)進行辨識濾波。文獻［3］對臨近空間高速高機動目標的跟蹤算法展開了深入研究，主要針對快速準確的航跡起始和穩(wěn)定的目標航跡維持來進行展開。文獻［4］介紹了幾種機動目標跟蹤模型和濾波算法，并針對影響跟蹤濾波性能的因素進行了分析，引出了自適應跟蹤濾波算法研究的必要性。文獻［5］提出了一種自適應網(wǎng)格交互多模型算法，該方法能處理自適應時變模型集合，隨時調(diào)整當前時刻使用的模型數(shù)量，并通過數(shù)值仿真結(jié)果驗證了所提算法有效性。

攔截彈使用目標跟蹤濾波器可以提取目標機動信息［6］，之后利用這些信息對目標軌跡進行預測，獲取彈目交會區(qū)域，進而計算出攔截彈所要滿足的中制導約束［7-8］。目前，針對高超聲速滑翔飛行器軌跡預測的主要方法也分為基于動力學和運動學模型兩大類。文獻［6］基于目標運動學模型，使用最小二乘法對目標彈道角和速度進行多項式擬合，以此預測目標運動狀態(tài)。文獻［9］基于目標動力學模型提出一種多層遞階預測軌跡方法，該方法借鑒多層遞階預測理論對預測模型進行隨機補償，將軌跡預測問題分解成氣動參數(shù)和模型誤差的混合預測以及在此基礎上對目標軌跡的預測。文獻［10］則針對臨近空間高超聲速滑翔飛行器機動能力強、軌跡靈活多變的特點，提出了一種基于典型控制規(guī)律的高超聲速滑翔飛行器軌跡預測方法。

攔截彈通過對目標進行跟蹤和預測，獲取了目標的運動狀態(tài)，之后將其用于制導律中，可以補償目標機動帶來的誤差。目前常用的制導設計方法有滑模制導律、最優(yōu)制導律以及改進的比例制導律。例如，文獻［6］用擴展卡爾曼濾波估計目標信息，將預測命中時刻目標運動軌跡作為虛擬目標點，采用針對虛擬目標的中制導攔截策略攔截目標。文獻［11］針對臨近空間防御作戰(zhàn)問題，設計了考慮零控攔截的中制導段最優(yōu)彈道修正算法。文獻［12］以導彈逆軌攔截高速運動目標為背景，運用間接高斯偽譜法設計帶攻擊角度約束的最優(yōu)中制導律。文獻［13］基于導彈和目標相對運動方程，設計了視線角約束自適應滑模中制導律。

本文針對HTV2 目標，提出了一種基于軌跡跟蹤預測的反HTV2 中制導算法。使用自適應網(wǎng)格交互多模型（adaptive grid interacting multiple model，AGIMM）濾波器對目標進行跟蹤，之后對預測交班點處目標彈道角和速度進行預測，并將其作為中制導約束，通過最優(yōu)中制導獲取滿足多約束條件的最優(yōu)彈道。

1 HTV2 運動特性

要對HTV2 目標進行攔截，需要先對其運動特性進行分析，選取合適的攔截段和攔截策略。HTV2 的運動過程包括主動段、再入段、無動力滑翔段和末制導段。其飛行軌跡靈活多變，橫縱向機動性能強。典型的縱向運動有平衡滑翔和跳躍滑翔，橫向運動有擺動式和轉(zhuǎn)彎式機動［1］。平衡滑翔是幾乎沒有起伏的平穩(wěn)滑翔，此時飛行器在大氣層內(nèi)進行無動力滑翔過程中所受的升力縱平面分量、重力以及離心力達到平衡。只有當選取合適的攻角，使式（1）成立時，飛行器才會進行平衡滑翔運動［2-3］。

式中：L為飛行器所受升力；σ為飛行器的傾側(cè)角；v為速度；m為目標質(zhì)量；r為目標地心距。

當平衡滑翔條件不滿足時，運動軌跡會變?yōu)樘S滑翔。跳躍滑翔的控制量攻角α經(jīng)常采用常值、分段線性或線性等簡單函數(shù)來描述，其跳躍幅度呈衰減振蕩形式［3］。圖1 為攻角取常值時，目標的彈道圖。從圖中可以看出，此時縱向為跳躍滑翔機動，側(cè)向為轉(zhuǎn)彎機動。

圖1 HTV2 常值攻角下的跳躍滑翔彈道Fig.1 HTV2 jump glide trajectory under constant angle of attack

目標側(cè)向機動主要靠傾側(cè)角控制，傾側(cè)角可以改變升力在縱向和側(cè)向上的分量，進而改變飛行器的側(cè)向運動姿態(tài)，實現(xiàn)側(cè)向機動突防。典型的側(cè)向機動一般為轉(zhuǎn)彎機動和擺動式機動。側(cè)向最大機動過載為2。圖2 為目標進行側(cè)向擺動式機動的彈道。

圖2 HTV2 側(cè)向擺動機動彈道圖Fig.2 HTV2 lateral swing maneuver trajectory diagram

根據(jù)上述分析可知，HTV2 的機動形式多樣，機動時機不定，導致目標彈道多變。因此，攔截彈要想成功攔截目標，需要對目標機動進行穩(wěn)定跟蹤和預測，以便根據(jù)目標機動做出及時調(diào)整。

2 目標跟蹤濾波器

2.1 機動模型建立

假設目標做勻轉(zhuǎn)速運動，角速率Ω和速度向量v相互垂直，即=0，因此有

式中：ω為轉(zhuǎn)彎速率。

根據(jù)二階馬爾科夫過程，勻轉(zhuǎn)速模型（constant turn，CT）轉(zhuǎn)彎模型可建模為［4］

式中：W為高斯白噪聲。

采樣時間為T，其離散形式變?yōu)?/p>

當ω＞0 時，運動方程描述向左轉(zhuǎn)彎；當ω＜0時，運動方程描述向右轉(zhuǎn)彎；當ω=0 時，模型近似常值速度運動［5-6］。

建立量測方程如下：

式中：Z(k)=ym，ym為目標沿y軸方向的位置量測值；H=(1，0，0)；V(k)為量測噪聲，設定為零均值白噪聲。

2.2 跟蹤問題可觀性分析

對于如下線性系統(tǒng)：

若滿足式（11）條件，則該線性系統(tǒng)是完全可觀測的。

故該系統(tǒng)滿足條件（11），是完全可觀測的。

2.3 AGIMM 跟蹤濾波器

目前，AGIMM 濾波器是跟蹤高超聲速滑翔飛行器最有效的算法之一。相對于傳統(tǒng)的模型參數(shù)固定的交互式多模型（interacting multiple model，IMM）算法，AGIMM 算法在跟蹤過程中能夠自適應調(diào)整轉(zhuǎn)彎角速率，使得模型集能夠隨目標機動性的變化實時地調(diào)整動態(tài)，從而實現(xiàn)對目標的精確跟蹤［6］。

AGIMM 濾波器由3 種勻轉(zhuǎn)彎模型組成，分別為向左轉(zhuǎn)彎模型、向右轉(zhuǎn)彎模型和中心轉(zhuǎn)彎模型。其中，中心轉(zhuǎn)彎模型的轉(zhuǎn)彎速率是上一幀3 種模型概率加權和［6］，即

式中：ωL(k)，ωC(k)，ωR(k)為3 種模型的轉(zhuǎn)彎角速率；μL(k)，μC(k)，μR(k)為3 種模型的模型概率。

模型調(diào)整，分為以下3 種情況。

（1）當μC(k)=max{μL(k)，μC(k)，μR(k)}成 立時，中心模型為最優(yōu)模型。此時，向左轉(zhuǎn)彎速率調(diào)整見式（15），向右轉(zhuǎn)彎速率調(diào)整見式（16）。

式 中 ：λL(k)=max {ωC(k) -ωL(k)，δω}；λR(k)=max {ωR(k) -ωC(k)，δω}，δω為最小網(wǎng)格間距；t1為無效模型的閾值。

（2）當μL(k)=max{μL(k)，μC(k)，μR(k)}成 立時，向左轉(zhuǎn)彎模型為最優(yōu)模型。此時，向左轉(zhuǎn)彎速率調(diào)整見式（16），向右轉(zhuǎn)彎速率調(diào)整見式（17）。

式中：t2為無效模型的閾值。

（3）當μR(k)=max{μL(k)，μC(k)，μR(k)}成 立時，向右轉(zhuǎn)彎模型為最優(yōu)模型。此時，向左轉(zhuǎn)彎速率調(diào)整見式（18），向右轉(zhuǎn)彎速率調(diào)整見式（19）。

通過上述3 種情況對模型進行調(diào)整之后，按IMM 算法對目標進行濾波跟蹤。AGIMM 算法對初始模型集具有較強的依賴性，需要根據(jù)目標的運動特性設置合適的初始參數(shù)，提高濾波快速性。

3 目標軌跡預測

目前，針對高超聲速滑翔飛行器軌跡預測的方法主要分為基于動力學和運動學模型兩大類?；趧恿W模型的方法通常先通過跟蹤估計出氣動參數(shù)、升阻比或控制參數(shù)的變化趨勢，然后研究其變化規(guī)律或者借助函數(shù)擬合工具給出未來時刻的變化趨勢，進而代入動力學模型進行軌跡預測。但是該方法需要通過前期情報搜集和飛行器反設計，對飛行器的本體參數(shù)精確已知，并且需要根據(jù)飛行器發(fā)射位置、光電特性等方面的信息對飛行器進行有效識別［7］。基于運動學模型的方法，首先通過濾波估計目標的運動參數(shù)，為預報提供準確的初值。之后，用解析法計算或數(shù)值方法迭代獲取目標運動軌跡。該方法需要目標的持續(xù)跟蹤信息，但不需要高超聲速滑翔飛行器的先驗信息［8-10］。

本文基于目標運動學模型對目標軌跡進行預測，無需獲取高超聲速滑翔飛行器的先驗信息。通過第2 節(jié)所給AGIMM 濾波器對目標進行跟蹤，可以獲取目標的速度信息，之后使用式（20）進行計算，獲取目標的彈道角信息。

式中：θT為目標彈道傾角；ψvT為目標彈道偏角；vT為目標速度；vTx，vTy，vTz為目標在發(fā)射慣性坐標系下的分速度。

由于目標運動在短時間內(nèi)具有延展性，據(jù)此可以建立目標運動模型，根據(jù)目標跟蹤結(jié)果，使用最小二乘法來計算模型系數(shù)［7］。一般情況下，可以將目標彈道角和速度使用多項式進行描述，如式（21）所示。其中，模型系數(shù)會隨著量測信息的變化而不斷更新。

實際上，當目標進行縱向跳躍滑翔時，彈道傾角會呈現(xiàn)出跳躍衰減的變化趨勢；當目標進行橫向擺動式機動時，彈道偏角會呈現(xiàn)出跳躍變化的趨勢。此時，可以建立正弦衰減模型來描述目標的運動狀態(tài)。如式（22）所示。

計算得到模型系數(shù)后，將預測交班時間代入式（21）或（22），即可獲取交班時刻的預測彈道角。其中，交班時間的預測可通過式（23）計算得出。

4 反HTV2 中制導算法

由第1 節(jié)的分析可知，HTV2 的典型彈道在縱向平面是跳躍滑翔式的。當?shù)孛胬走_探測到HTV2后，使用雷達對其進行跟蹤，當判斷出目標從高點下降后，發(fā)射攔截彈對目標進行攔截。初制導段結(jié)束后，攔截彈進入中制導段，中制導段包括有動力和無動力段2 個階段，其中有動力段主要是繼續(xù)增加導彈速率，并結(jié)合彈體的姿態(tài)控制，將導彈的速度快速指向需要的方向。在姿態(tài)控制階段結(jié)束后，攔截彈進入中制導的無動力階段。

在縱向平面內(nèi)，攔截彈在中末交班時刻需要滿足側(cè)窗探測視線角約束、視線角速率約束［12］，如式（24）～（26）所示。

式中：tf為中末交班時刻；qε為彈目視線傾角；qεd為交班時刻的期望視線傾角；qεmin為側(cè)窗可探測的最小視線傾角；qεmax為側(cè)窗可探測的最大視線傾角；為視線傾角角速率。

此外，攔截彈在中末交班時刻也需要對零控脫靶量進行約束，使其在末制導修偏能力范圍之內(nèi)，保證目標可以被捕獲，約束形式如式（27）所示。當交班時刻零控脫靶量為0 時，認為末制導階段不加控制也能實現(xiàn)對目標的直接碰撞。

式中：ZEM為中末交班時刻的零控脫靶量；ZEMmax為末制導所能修偏的最大脫靶量。

根據(jù)文獻［6］中的研究，假設目標不機動，目標和攔截彈的速度比，那么在中末交班過程中攔截彈和目標存在零控攔截條件，且該條件由式（28）唯一確定：

θ(tf)=θd=qεd+arcsin(psin(θT(tf) -qεd))，（28）式中：θ(tf)為交班點處導彈的彈道傾角；θT(tf)為交班點處的目標彈道傾角；p=為交班點處的彈目速度比，可通過第2，3 節(jié)軌跡跟蹤和預測方法獲取。

在側(cè)向平面內(nèi)，攔截彈同樣需要滿足交班時刻彈目視線角和視線角速率約束。此外，考慮到末制導階段采用側(cè)窗探測的方式，攔截彈在中制導段內(nèi)，需要根據(jù)目標機動情況，進行彈道規(guī)劃，使目標能始終在攔截彈的同一側(cè)，以便為交班時刻導引頭捕獲目標創(chuàng)造良好的條件。側(cè)向平面約束條件，如式（29）～（31）所示。

式中：qβ為彈目 視線偏 角；qβd為 交班時刻的期望視線偏角；qβmin為側(cè)窗可探測的最小視線偏角；qβmax為側(cè)窗可探測的最大視線偏角為視線偏角角速率。

在縱向平面內(nèi)，使用最優(yōu)中制導律可以實現(xiàn)上述多約束條件。首先，建立彈目縱向相對運動方程：

式中：ay為導彈縱向加速度；v為導彈速度，aTy為目標縱向加速度。

對式（32）第2 項求導，則有

設定f=Cδ，其中δ=aTy，由于目標加速度大小是有界的，可推知擾動量f有界。

攔截彈在中末交班時刻期望的狀態(tài)量x(tf)=(0，0，θd)T，選取性能指標函數(shù)如下：

式中：G=，?為一正數(shù)，tj為剩余交班時間；F，Q和G分別表示末端狀態(tài)量、狀態(tài)量和控制量的權重系數(shù)。

在側(cè)向平面內(nèi)，使用最優(yōu)中制導律對視線角和視線角速率進行約束，即可實現(xiàn)中制導策略。令狀態(tài)量x=(qβ-qβd，)T，x(tf)=(0，0)T。建立如下狀態(tài)方程：

性能指標函數(shù)和參數(shù)選取與縱向平面相同。對于該最優(yōu)控制問題，常用高斯偽譜法來求解最優(yōu)制導指令［13-16］。圖3 為基于軌跡跟蹤預測的攔截彈中制導算法流程圖。

圖3 基于軌跡跟蹤預測的攔截彈中制導算法Fig.3 Midcourse guidance algorithm of interceptor missile based on trajectory tracking prediction

5 仿真校驗

假設HTV2 采用恒攻角控制，進行縱向跳躍滑翔。攔截彈縱向和側(cè)向均采用上述中制導策略，末制導使用比例導引。仿真開始時刻，目標位置為xT=820 km，yT=40 km，zT=3 km，目標速率為vT=4 500 m/s，目標彈道角為θT=8.78°，ψvT=0°，目標采用恒攻角控制，αT=10°。導彈位置為x=0，y=0，z=0，速度為vx=0，vy=0，vz=0，導引頭探測距離為Rj=80 km，攔截彈交班時刻期望視線角為qεd=5°，qβd=0°，側(cè)窗探測范圍為qε∈[2°，30°]，qβ∈[-4°，4°]，零控脫靶量約束參數(shù)為ZEMmax=1.2 km。目標在經(jīng)過65 s 的飛行后，已處于高度下降飛行狀態(tài)，進入地面雷達的探測范圍，此時發(fā)射攔截導彈，隨后，將此發(fā)射時刻記為0 時刻。

5.1 目標跟蹤濾波器性能仿真校驗

針對上述跳躍滑翔HTV2 目標，分別使用AGIMM 濾波器、基于CT 模型的IMM 濾波器與基于Singer 模型的卡爾曼濾波器（Kalman filter，KF）對目標進行跟蹤估計。其中，IMM 濾波器選用和AGIMM濾波器相同的模型集，但是模型集的參數(shù)是固定不變的。圖4為目標y向加速度估計對比圖，圖5為目標y向位置跟蹤對比圖，圖6為目標y向速度估計對比圖。

圖4 y 向加速度濾波估計Fig.4 y-direction acceleration estimation by filtering

圖5 y 向位置濾波跟蹤Fig.5 y-direction position estimation by filtering

圖6 y 向速度濾波跟蹤Fig.6 y-direction velocity estimation by filtering

仿真結(jié)果采用均方根誤差（root mean square error，RMSE）進行評估，其計算公式為

經(jīng)過100 次蒙特卡羅仿真，計算出各濾波算法的均方根平均誤差如表1 所示。

表1 各濾波算法的RMSE 比較Table 1 RMSE comparison of filtering algorithms

根據(jù)仿真圖以及表1 可知，相比另外2 種濾波器而言，使用AGIMM 濾波器跟蹤該HTV2 目標，所獲取的目標信息更接近真實值，具有更高的濾波精度。

5.2 目標軌跡預測仿真校驗

當攔截彈進入中制導主動段后，攔截彈開始對目標的彈道角進行預測，之后計算出零控攔截交班所期望的彈道角，將其應用于制導律當中。本文使用最小二乘法對目標進行預測，預測結(jié)果與目標真實運動信息的對比圖如圖7～9 所示。圖7 為目標彈道傾角預測值和真實值的對比，圖8 為目標彈道偏角預測值和真實值的對比，圖9 為目標速度對比圖。

圖7 目標彈道傾角預測值與真實值對比圖Fig.7 Comparison diagram of predicted and real target trajectory

圖8 目標彈道偏角預測值與真實值對比圖Fig.8 Comparison diagram of predicted and real target trajectory inclination angle

圖9 目標速度預測值與真實值對比圖Fig.9 Comparison diagram of predicted and real target trajectory deviation angle

從圖中可以看出，彈道角和速度的預測值與真實值之間具有相同的變化走向，兩者之間的誤差如表2 所示。

表2 目標預測值與真實值誤差平均值Table 2 Average value of the error between the predicted value and the true value

5.3 反HTV2 中制導算法仿真校驗

假設攔截彈交班時刻期望視線角為qεd=5°，qβd=0°，側(cè)窗探測范圍為qε∈[2°，30°]，qβ∈[-4°，4°]，零控脫靶量約束參數(shù)為ZEMmax=1.2 km。選取高斯節(jié)點數(shù)為N=15，性能指標函數(shù)權重系數(shù)選取Q=15I，F(xiàn)=3I，ζ=20，其中I為單位矩陣。

針對該HTV2 目標，所得攔截彈彈道仿真圖如圖10 所示。該彈道中末交班時間為79.61 s，交班高度為32.97 km，交班時視線角速率為0.002 4（°）/s，視線角為qεd=0.56°，qβd=-1.69°，零控脫靶量為ZEM=36.20 m。可知交班時刻視線角速率較小，接近于0，視線角在側(cè)窗約束范圍內(nèi)，且零控脫靶量也滿足約束。攔截彈在91.28 s 成功攔截到目標，最終脫靶量為0.381 0 m，小于0.5 m，滿足直接碰撞攔截的要求。

圖10 攔截彈道仿真圖Fig.10 Intercept trajectory simulation diagram

對上述算例做100 次蒙特卡羅仿真，結(jié)果如圖11 所示，統(tǒng)計結(jié)果如表3 所示。

圖11 蒙特卡羅仿真脫靶量Fig.11 Monte Carlo results of miss distance

表3 100 次蒙特卡羅打靶統(tǒng)計結(jié)果Table 3 Monte Carlo shooting statistics

當目標以不同常值攻角跳躍滑翔時，攔截彈攔截彈道簇如圖12 所示。此時攔截彈在中制導階段目標彈道傾角變化如圖13 所示。彈道偏角變化如圖14 所示。

圖12 目標攻角變化時對應攔截彈道簇Fig.12 Intercept ballistic cluster when the target angle of attack changes

圖13 目標攻角變化對應攔截彈彈道傾角變化圖Fig.13 Interceptor inclination angle diagram when the target angle of attack changes

圖14 目標攻角變化對應攔截彈彈道偏角變化圖Fig.14 Interceptor deflection angle diagram when the target angle of attack changes

通過第1 節(jié)的分析可以得知，HTV2 的機動方式除了在縱向進行跳躍滑翔以外，在側(cè)向也會采用擺動式機動進行突防。當目標分別從高度60 km 開始進行側(cè)向機動，目標傾側(cè)角初始指令正、負皆有，則攔截彈道簇如圖15 所示，這些彈道的脫靶量均小于0.5 m，可實現(xiàn)直接碰撞毀傷目標。此時攔截彈對應的彈道傾角變化如圖16 所示，彈道偏角變化如圖17 所示。

圖15 不同側(cè)向機動時機對應攔截彈道簇Fig.15 Different lateral maneuver timing lead to interceptor trajectory cluster

圖16 側(cè)向機動時對應攔截彈道傾角變化Fig.16 lateral maneuver lead to the change of the interceptor inclination angle

圖17 側(cè)向機動時對應攔截彈道偏角變化Fig.17 Lateral maneuver lead to the change of the interceptor deflection angle

通過上述仿真驗證可知，本文提出的基于軌跡跟蹤預測的反HTV2 中制導算法可以滿足中末交班時的視線角、視線角速率和零控攔截彈道角約束。當目標進行機動時，該制導律可以根據(jù)目標變化調(diào)整彈道角約束，從而滿足零控脫靶量要求。

6 結(jié)束語

本文針對HTV2 目標，提出了一種基于軌跡跟蹤預測的反HTV2 中制導算法。該算法采用AGIMM 濾波器在中制導階段對HTV2 目標進行跟蹤。將濾波結(jié)果作為量測值用于軌跡預測中，使用基于遺忘因子的最小二乘法對目標軌跡進行預測，獲取中末交班時刻的預測彈道角和速度并用于計算攔截彈中制導約束條件。本文在考慮了視線角、視線角速率、零控攔截等約束條件的基礎上，使用最優(yōu)制導律進行制導。為了驗證上述中制導算法的性能，進行了仿真驗證。針對給定的目標初始條件，首先進行了目標軌跡跟蹤和目標軌跡預測的仿真，之后使用所設計的中制導算法進行攔截彈道仿真。通過仿真結(jié)果可知，本文所提基于軌跡跟蹤預測的反HTV2 中制導算法可以滿足中末交班的視線角、視線角速率和零控攔截彈道角約束，且能夠成功攔截不同機動形式下的HTV2 目標。