遠程攔截飛行諸元快速裝訂算法

楊繼龍，周 軍，呼衛軍

(西北工業大學 精確制導與控制研究所，西安 710072)

遠程攔截飛行諸元快速裝訂算法

楊繼龍，周 軍，呼衛軍

(西北工業大學 精確制導與控制研究所，西安 710072)

傳統諸元解算僅強調終端打擊精度而忽略目標的高速運動趨勢，易引起攔截末段過載飽和，影響攔截效能。為解決此問題，針對多級推進的地基遠程攔截彈結合解析計算和數據擬合算法，設計一種基于預測攔截點(PIP點)的在線諸元計算方法。首先，該方法通過引入攔截彈到達PIP點時彈目速度交會角約束，基于橢圓彈道理論實現對關機點參數解析求解，進而利用目標距PIP點飛行時間約束設計收斂策略，快速確定發射時刻，求解出發射方位角初值。然后，結合離線建立的近似解析多項式，完成飛行諸元初值求解。最后，利用牛頓法僅少數步迭代即可確定出攔截彈發射諸元精確值。六自由度仿真表明，該方法不僅具有很好的在線能力，攔截精度高，而且適用于大范圍機動發射。

預測命中點；遠程攔截；關機點參數；交會角；機動發射

0 引言

隨著導彈技術的發展，攻防對抗將成為未來戰爭的主要形態。目前，多國正積極構建以大氣層外遠程攔截彈為中堅的多層次彈道導彈防御系統，如美國GBI、俄羅斯S400等。諸元的有效解算是此類系統重要戰斗保障。因此，針對遠程中段攔截彈設計更快速、高效的諸元解算方法顯得十分必要。目前，因國外導彈技術封鎖，少見國外相關諸元資料，且由于諸元解算具有強針對性，國內諸元計算算法設計主要還集中在應用層面上，采用積分、射表擬合以及自由段解析等方法。積分法[1-3]使用簡化模型，考慮約束較少，需要進一步改進；射表擬合法[4]則需要綜合運用均勻實驗設計、逐步回歸分析等理論，離線計算量大，尤其當諸元涉及參數較多時，數據存儲量巨大；自由段解析法形式復雜，精度較差，一般不作單獨使用。近年來，提出一種擬合+迭代[5-7]的諸元計算思路，當選定合適擬合函數式，求解合理的初值，通過有限步迭代即可獲得精確的諸元參數，一定程度上綜合了前3種方法的利弊。但以上的多種計算方法主要針對彈道導彈設計，只強調對地面固定目標的打擊精度，并不適用于攔截高速運動目標的遠程動能攔截彈。

本文研究機動發射的三級固體助推地基中段遠程攔截彈對彈道導彈攔截諸元計算問題，考慮地球自轉，為其設計一條通過發射點與PIP點，且滿足攔截剩余飛行時間，以及彈目速度交會角約束的標準彈道，并確定相應的發射時間、發射方位角以及飛行程序角等諸元參數：首先，根據設定的目標PIP點位置與速度，確定攔截彈彈道約束角度，利用橢圓軌道方程，并結合通過離線開環飛行仿真數據庫所擬合出飛行程序參數與諸元參數之間近似解析式，在大范圍機動發射情況下，實現對目標外推彈道上滿足發射時間約束和飛行參數取值范圍約束的PIP點的選取，并快速求解該點所對應的攔截彈發射時刻與諸元計算初值。然后，利用有限次數值積分迭代，即可確定飛行參數精確值。

本文方法主要基于解析計算結合數據擬合，具備諸元快速在線解算能力，且通過引入反映目標運動趨勢的交會角約束，在順軌與逆軌攔截高速目標情況下，該約束角度對平滑末段攔截彈指令過載、提升攔截效能具有重要實際意義。

1 攔截彈關機點參數求解

1.1 問題描述

遠程攔截彈飛行主要分為助推段與自由飛行段。諸元解算階段，本文主要考慮連續助推方式推送攔截彈頭滿足命中橢圓彈道參數要求，保證攔截彈頭無控情況下自由滑行，精確與目標碰撞。

彈道示意圖見圖1。

設OkP之間橢圓弧段為命中彈道，則Ok點滿足命中彈道的參數即為所求的關機點參數，但因Ok點位置隨飛行程序參數變化較大，不易確定。考慮以攔截彈助推段關機點平均地心距ravg作為關機點地心距rk初值，關機點平均射程角θavg為關機點射程角Ok初值，θt作為初始射程角，且滿足：

(1)

此時，利用PIP點位置、速度矢所確定的終端交會角約束，求解命中彈道參數Vk、γk，作為攔截彈關機點參數，且上述ravg與θavg可根據2.2節擬合解析式數據儲備中計算作為初值，為精確求解，后續解算需進一步更新。

1.2 命中彈道求解

在實際作戰中，若已知目標PIP點速度矢與位置矢，可確定命中時刻攔截彈所需當地彈道傾角約束γpm，考慮典型順軌與逆軌場景，表述如下：Vt、Vm分別為命中時刻目標與攔截彈速度矢，P為PIP點，γpt、γpm分別為攔截點攔截彈與目標當地彈道傾角，xPy設為攔截射面，xPz設為攔截點所在切平面，VtL為目標速度Vt在攔截面內投影。由圖2分析，選取：

(1) 廣義逆軌攔截，約束角度取γpm=-γpt；

(2) 廣義順軌攔截，角度約束取γpm=γpt。

可知，在攔截面內，彈目視線角速率在PIP點收斂為零，且目標相對攔截面法向速度分量最小。因此，該約束角度選取方式對減小攔截末段彈目視線角速率以及平滑末段攔截彈過載指令具有實際意義。

確定攔截彈命中時刻約束角度γpm，假設攔截彈在理想平方反比引力場作用下，攔截彈關機點的彈道傾角γk可解析求解，推導如下：

根據橢圓軌道方程相關公式[8]：

(2)

(3)

h=rkVkcosγk

(4)

P=rkv0cos2γk

(5)

(6)

其中，μ為引力常數；e為軌道偏心率；v0為能量參數；h為動量矩；P為半通徑；f為真近點角；r為地心距；已知P與h為常值[9]。

當到達PIP點時彈道傾角為γpm，設此刻攔截彈的真近點角為 ，則由式(6)：

(7)

根據真近點角與γpm之間的關系：

(8)

可得

(9)

代入式(2)～式(9)中，平方得

(10)

將式(4)代入式(10)，取K1=rp/rk得

(11)

回代式(5)～式(11)，得

(12)

令K2=cos2γpm(2K1-2K12)，則變為

(13)

對兩點邊值問題的求解[10-11]，在γk已知的情況下，求解需要速度：

(14)

代入式(3)，并取K3=K1(1-cosθt)，即可得

(15)

代入式(15)～ 式(13) ，并交叉相乘，并設如下變量：

結果整理為

Atan2γk+Btanγk+C=0

(16)

把終端角度約束問題轉化成一元二次方程。根據式(14)，關機點的速度Vk由初始彈道傾角γk決定，若γk選取不當，所求的攔截彈速度達到逃逸速度，不可達PIP點，由v0<2限定攔截彈沿橢圓彈道飛行[12]，得

(17)

得出Vk上限，可看出Vkmax為逃逸速度。將式(17)代入式(14)得γk，取值為γk∈(γkmin,γkmax)。

(18)

由上式推導可看出，γk的求解取決于關機點地心距與PIP點地心距的比值K1、兩地心距間的夾角θt以及終端當地彈道傾角約束γpm，即

(19)

以赤道半徑a0=6 378 140 m為參考高度，利用數值計算，定量分析γk分別與K1、θt、γpm之間的關系，可看出在一般攔截范圍內，總有

Δ=B2-4AC>0

(20)

PIP高度為1 000 km時，Δ取值見圖3。

但上述推導公式中，角度約束以cosγpm形式出現，故γk有兩解，帶入原方程驗證，確定可行解。

1.3 發射時刻確定

當選定PIP點，以目標發射時間為基準零值，設目標當前時刻tCur，自由飛行到PIP點時間為tp，則確定tLeft=tp-tCur為剩余飛行時間，設攔截器發射時刻tLau，則需要滿足約束條件：

(21)

即若設攔截彈發射時刻相對目標發射時間為tILau，則tLau=tILau-tCur。不考慮地球自轉時,根據PIP點狀態與發射點位置，由式(14)與式(16)解析求解關機點參數，繼而利用自由段飛行時間[10]，直接求解發射時刻tLau。

(22)

在考慮地球自轉情況下，選擇不同發射時刻，攔截彈發射點與PIP相對位置存在變化[13]， 需要對tLau修正，直到|tLaun-tLaun-1|≤ε。

1.4 關機點參數迭代求解

綜述以上內容，本節主要基于發射點與目標外推彈道，引入攔截終端角度約束，結合橢圓軌道方程，考慮地球自轉情況下，在目標彈道上選取滿足攔截彈發射時刻tLau>0的PIP點，確定攔截彈在rk=ravg情況下對應的關機點參數γk、Vk，具體解算流程如圖4所示。

1.5 關機點參數求解分析

在上述仿真解算流程中，針對rk=ravg時關機點參數求解時，涉及tLau迭代計算，該迭代快速性與穩健性將直接影響到諸元計算的可用性以及作戰中的攔截效能，本節將作進一步分析。

1.5.1 關機點參數求解收斂性分析

(1)tLau迭代。主要補償地球自轉引起的相對位置變化，在求解過程中，僅tLau為自變量。選擇tLau初值后，根據式 (16) 式與式(22) ，解析計算tf，結合 (21)式約束關系，確定新tLau，循環進行實現收斂。

(2)針對tLau<0情況，在一般攔截場景中，目標飛行時間單調遞增。當PIP點依飛行時序選取，目標數千公里飛行過程中，完全可選出tLau。

1.5.2 關機點參數求解快速性分析

(1)根據下文表6中針對tLau仿真結果可知，tLau不僅大范圍選取，且僅2～3步可收斂到0.5 s誤差范圍內，不僅快速，而且驗證上述tLau收斂特性。

(2)在迭代解算過程中，直接解析求解，并不涉及數值積分。

(3)當內環解算tLau<0時，將基于當前PIP點直接外推|tLau|以上時間，僅一次更新，即可滿足tLau≥0。在后文中，默認選定PIP點均滿足tLau≥0。

2 諸元近似解析式擬合

2.1 方位角初值選取

諸元計算中，為尋求精確度較高的方位角初值，需充分考慮遠程攔截特點。目標慣性系下，某時刻狀態可經算法外推得到。當攔截tLau確定，PIP點相對于攔截彈發射點位置隨之確定。設發射點經度為λ0，地心緯度為φ0，PIP在導彈發射時刻對應經度為λ1、φ1，令Δλ=λ1-λ0，據球面三角形得

(23)

事實上，因地球自轉，攔截彈主動段彈道存在指向射面外的初速度引起彈道變形。假設發射系中的攔截彈主動段彈道為一平面曲線[14]，在主動段關機點k作切平面，如圖5所示。因發射點與主動段關機點較接近，它們切平面近似重合。V為關機點發射系射向速度，Vt為關機點發射慣性系射向速度。A0為式(23)方位角求解初值。

由下式推導出主動段彈道變形角為ΔA0：

(24)

Ve=r0cosφ0we

(25)

(26)

其中，r0為發射點的地心距，Ve為向東初始速度，即修正后的發射方位角為

(27)

修正結果有足夠精度作為諸元精確迭代計算初值。2.2 諸元近似擬合解析式

遠程攔截所選PIP點在大氣層外，不考慮再入段，可采用更簡單的迭代計算。本文分析遠程攔截彈典型飛行程序：一級垂直發射，攻角轉彎，重力轉彎，分離段定軸飛行；二三級常值角速率轉彎。攻角轉彎段轉彎遵循以下規律：

(28)

設定發射參數與飛行程序參數范圍及間隔，在開環仿真中，對不滿足飛行過程中過載、動壓等約束條件的彈道予以剔除，考慮實際應用中作戰范圍 僅限于北半球取值。通過離線開環仿真，獲取對應的關機點參數數據庫，利用二次多項式擬合：

根據以上表達式,結合上述參數取值所保存的數據庫，離線對參數ai、bi、ci線性回歸解算，在選定發射點與初始方位角情況下，可根據關機點約束，快速求解飛行程序參數初始值。

2.3 諸元解算流程

(30)

利用牛頓法有限次迭代[1]，可得諸元精確解。具體諸元流程如圖6所示。

3 仿真算例

本文將參考地基中段攔截彈GBI性能參數建立六自由度助推段模型，設置典型的攔截場景，驗證本文設計的諸元計算方法有效性與快速性，攔截彈具體參數如表 1所示。

表1 攔截彈參數

以目標發射時間為參考基準，設tCur=180 s時目標以表3狀態開始自由飛行。

考慮J2項引力攝動外推彈道，獲取實時PIP點與對應tLeft。設置3組發射點分別模擬順軌、逆軌以及攔截面與目標彈道平面呈一夾角的攔截情況，見表4。

根據圖6仿真解算流程，先大步長迭代，選定PIP點，因不同攔截場景中存在多個可用PIP點，本節僅列寫目標最早可攔截點(即首次出現tLau>0時對應PIP點)的諸元初值(對應a)及利用牛頓迭代確定的諸元精確值(對應b)，并給出迭代次數，說明諸元初值可用性及求解精確值的快速性。諸元初值與精確值見表5。

表2 攔截彈初始條件

表3 目標初始狀態(地心慣性系)

表4 發射點設置

為本文方法完備性，列出3組攔截場景中tLau收斂特性，如表6所示。給定任一發射時刻初值，tLau僅2～3步收斂到0.5 s誤差范圍內。

通過六自由度彈道仿真，繪制典型攔截場景中三維仿真攔截圖，給出攔截時刻彈目交會約束角度誤差以及PIP點位置，見圖7和圖8。

表5 諸元初值與精確值

表6 tLau迭代結果

為驗證諸元算法的適應性，僅以上述逆軌攔截場景為例，進一步外推目標彈道，選取PIP點，即在不同的tLau時，驗證該裝訂算法。

由仿真可知，當首次出現tLau>0后，依目標飛行時序選取多組PIP點，依然可實現諸元解算。在實際攻防對抗中，通過本文算法對目標彈道逐點解算，確定可攔截弧段，為組織多次攔截提供參考。

4 結論

(1)綜合考慮目標PIP點位置與速度，首次引入攔截時刻彈目速度交會角度約束，基于橢圓軌道方程，實現關機點參數解析求解。

(2)通過文中表格及三維圖可看出，該方法不僅適用于機動發射條件下不同場景的諸元裝訂任務，而且諸元初值具有較高精度，僅2～3次迭代，實現脫靶量收斂到50 m內，快速且高效。

(3)尤其在逆軌攔截場景中，通過引入彈目速度交會約束角度，使得攔截末段彈目近 相對飛行，有助于末端導引頭對目標捕獲，提高攔截效能，具有工程應用價值。

[1] 王海麗,陳磊,胡小平.彈道導彈基本諸元的快速裝訂算法研究[J].國防科技大學學報,1999,21(2): 5-8.

[2] 雍恩米,唐國金.基于攝動制導的彈道導彈發射諸元的仿真算法[J].系統仿真學報,2005,17(5): 1048-1051.

[3] 陳春,劉玉文,楊啟仁.彈道積分決定射擊諸元的建模及精度分析[J].彈道學報,2001,3(1): 81-85.

[4] 多招平.彈道式導彈射擊諸元計算方法研究[D].中國航天科工集團第四總體設計部,2003.

[5] 潘剛.遠程彈道導彈快速諸元計算方法研究[D].國防科技大學研究生院,2002.

[6] 郝佳新,丁振國,陳克俊.彈道導彈基本諸元的快速裝訂算法研究[J].計算機仿真,2008,25(4): 25-29.

[7] 韋文書,荊武興,高長生.機動發射的彈道導彈飛行諸元的快速計算[J].哈爾濱工業大學學報,2012,44(11): 7-12.[8] LU P,PAN B F.Highly constrained optimal launch ascent guidance[J].Journal of Guidance,Control And Dynamics,2010,33(2):404-414.

[9] 賈沛然,陳克俊,何力.遠程火箭彈道學[M].長沙: 國防科技大學出版社,1993.

[10] Zarchan P.Tactical and strategic missile guidance[M].Progress in Astronautics and Aeronautics,AIAA,Washington,DC,2002.

[11] Steven L Nelson,Zarchan P.Alternative approach to solution of Lambert's problem[J].Journal of Guidance,Control and Dynamics,1994,15(4):1003-1009.

[12] Karthikeyan K,Ashok Joshi.Solution of Lambert's theorem with additional terminal constraints[R].AIAA 2014-2543.

[13] 荊武興,李羅剛,高長生．反導攔截飛行方案及時間窗口快速搜索算法[J].系統工程與電子技術,2013,35(6): 1256-1261.

[14] 黃岳,肖凡,張濤.彈道導彈發射方位角快速計算新方法[J].四川兵工學報,2009,30(6): 51-52.

[15] 劉君,陳克俊,謝愈,等.助推滑翔飛行器發射諸元計算方法研究[J].彈箭與制導學報,2012,32(6): 33-36.

[16] 賈沛然，陳克俊，何力.遠程火箭彈道學[M].長沙: 國防科技大學出版社,1993.

[17] 張濤,馬瑞萍,陳陽,等.彈道導彈諸元檢驗方法研究[J].航天控制,2011,29(1): 76-79.

(編輯：崔賢彬)

Effective algorithm for calculating flight firing data for long range interception

YANG Ji-long,ZHOU Jun,HU Wei-jun

(Institute of Precision Guidance and Control,Northwestern Polytechnical University,Xi'an 710072,China)

The traditional methods of firing data which only focus on final hitting accuracy and ignore the target motion trend with high speed are easy to lead to overloading saturation during the terminal phase of intercept,even decrease intercepting efficiency.To solve this problem,in view of ground-based long range interceptor propelled by solid fuel engine with multi-stages,this paper combines analytical calculation and data fitting to design an online firing data calculation algorithm based on prediction intercept point (PIP) .First of all,the algorithm analytically solves the shutdown point parameters based on the ellipse trajectory theory by introducing the speed intersection angle between interceptor and target when interceptor arrives at PIP,designs convergence strategy to quickly determine launch time by using constraint condition of flight time from target to PIP，and calculate the azimuth initial value.Then,initial firing data solution was completed by combining with the established offline approximate analytic polynomial.Finally,the exact firing data value was determined through a few iterations by using Newton method.6DOF simulation show that the algorithm not only has a good online ability to intercept high precision,but applies to a wide range of mobile launchers.

PIP；long range interception；shutdown point parameters；intersection angle；mobile launch

2015-11-10；

2015-11-27。

楊繼龍(1990—)，男，碩士生，主要從事飛行器制導、控制與仿真技術研究。E-mail:365776138@qq.com

V448

A

1006-2793(2017)01-0128-07

10.7673/j.issn.1006-2793.2017.01.023