機動發(fā)射的彈道導彈飛行諸元的快速計算

韋文書，荊武興，高長生

（哈爾濱工業(yè)大學航天工程系，150001 哈爾濱）

式中:φ為地表一點的地心緯度;B為該點的地理緯度;a、b分別為地球的長半軸和短半軸.

O、T兩點地心角和球面方位角能夠通過球面三角關系得到:

機動發(fā)射的彈道導彈飛行諸元的快速計算

韋文書，荊武興，高長生

（哈爾濱工業(yè)大學航天工程系，150001 哈爾濱）

基于解析預報與數(shù)值尋優(yōu)相結合的方法，本文研究機動發(fā)射的彈道導彈飛行諸元快速設計問題.首先，利用線性回歸方法建立了導彈諸參數(shù)與發(fā)射信息之間的多項式擬合關系式，利用該式可以得到飛行方案諸參數(shù)迭代初值;然后，以彈道落點縱向偏差、橫向偏差及彈道頂點高度偏差平方和極小為指標函數(shù)，采用Levenberg-Marquardt方法搜索諸參數(shù)精確解;最后，以民兵III彈道導彈為例，對所提方法進行了仿真驗證，研究結果表明，該方法計算量少，適合已定參數(shù)的導彈大范圍機動發(fā)射的快速計算，且在落點精度要求不高的情況下，可以用解析表達式的解作為彈道飛行方案諸元的值，對目標進行打擊.

解析預報;機動發(fā)射;Levenberg-Marquardt方法;多級彈道導彈;線性回歸方法

導彈飛行方案設計的目的是找到滿足所有系統(tǒng)要求的最適合的彈道［1］.對于機動發(fā)射的彈道導彈來說，針對不同的發(fā)射點和目標點等攻擊信息能夠快速選取飛行方案諸參數(shù)，將大大的提高導彈的作戰(zhàn)效能.所以，機動發(fā)射彈道導彈的飛行方案快速設計是1個重要的、亟待解決的內(nèi)容.目前，飛行方案設計包括直接法、間接法和智能優(yōu)化算法等［2－8］，然而這些方法通常需要大量的計算［9］，故對于機動發(fā)射的彈道導彈來說不能夠迅速計算其彈道諸參數(shù).

描述飛行器空間彈道的方法有3類:解析法、數(shù)值積分和函數(shù)逼近法.解析法計算速度快，實時性好，但解析解求取困難且精度受限;數(shù)值積分法精度高，缺點是計算復雜，實時性好;函數(shù)逼近法介于解析法和數(shù)值積分法之間.目前，國內(nèi)外有很多學者致力于函數(shù)逼近法的研究［10－13］.如陸平［14］在研究多級火箭上升段制導方案過程中，將飛行高度擬合成1個多項式，從而設計了二級火箭的反饋控制律.

王明海［15］等將落地射程偏差和方位角偏差的二元迭代分解成兩個一元迭代，從而提供了彈道導彈彈道設計的一種快速迭代方法，該方法對多變量、有約束的彈道設計有一定缺陷;黃岳［16］等以發(fā)射方位角為射擊諸元，對球面方位角迭代得到發(fā)射方位角，該方法計算時間較長，不適合機動發(fā)射的彈道導彈彈道設計;劉剛、何麟書等［17］基于微積分、變分原理提出一種多級火箭飛行方案設計的工程方法，該方法考慮工程約束較少.

本文研究機動平臺發(fā)射的三級固體彈道導彈全彈道飛行方案設計，不同的發(fā)射點緯度、射向、射程對應不同的飛行方案彈道，彈道諸元的快速選取對導彈的效能有重大影響.為了減少大量計算，本文將解析預報與優(yōu)化方法相結合，能夠快速給出最優(yōu)的期望方案.首先，對彈道進行分段，根據(jù)每段的飛行特性，設定了相應的約束條件，設計了飛行方案;其次，通過開環(huán)仿真得到了飛行方案諸元與發(fā)射點、目標點等信息的數(shù)據(jù)庫，利用多項式擬合他們之間的解析關系式，從而，對一定打擊范圍內(nèi)任意的發(fā)射點、目標點及期望彈道頂點高度，能夠迅速給出飛行方案諸元的近似值;然后，以解析表達式得到的近似值為初值，縱向、橫向、及高度偏差的平方和極小為目標函數(shù)，采用L-M法進行迭代計算，得到諸元的精確值;最后，對文中方法進行了仿真，結果表明:采用解析預報能快速給出諸元概略值，以此為初值，采用數(shù)值方法只需少量計算即可得到精確值;當落點精度要求不高時，可以用解析預報的值直接作為諸元的真值來打擊目標.

1 數(shù)學模型

彈道導彈運動模型基于六自由度彈道仿真模型建立.考慮地球自轉、標準橢球模型，不考慮級間分離干擾誤差，給定發(fā)射方位角與主動段飛行程序角后，導彈即可按預定彈道飛行.

發(fā)射坐標系［18］中導彈的數(shù)學模型為

其中:R為導彈相對發(fā)射系的位置;V是導彈相對發(fā)射系的速度;T為推力;m為t時刻質量;D為氣動阻力;L為氣動總升力，T為推力幅值;Isg為比沖;ac為科氏加速度;ae為牽連加速度.

引力加速度的表達式是

其中:r為地心到飛行器的矢量;ωe為地球自轉角速度矢量，則

式中:r為地心到飛行器的距離;μ為地球引力常數(shù);Re為地球赤道平均半徑;J2為帶諧系數(shù);φ為地心緯度.

科氏加速度的表達式為

牽連加速度的表達式為

彈道導彈一般是軸對稱的，定義其體坐標系為obxbybzb.由氣動力的定義可以得到，氣動阻力和氣動總升力為

式中CD、CL分別為阻力系數(shù)和總升力系數(shù)，它們均為攻角α和馬赫數(shù)Ma的函數(shù).

發(fā)動機的推力矢量沿著obxb軸方向，利用發(fā)射坐標系到體坐標系的轉換矩陣A，可得到推力在發(fā)射坐標系中的投影，其中

式中:s代表 sin;c代表 cos;γ、φ、ψ分別為滾轉角、偏航角、俯仰角.

2 工程描述

對于多級固體彈道導彈，其主動段飛行方案有3種模式:連續(xù)推力模式、連續(xù)推力+滑行(小推力末端修正)、間隔推力模式.本文選取連續(xù)推力模式.

飛行程序通常指的是主動段飛行時程序角的變化規(guī)律.其選擇原則為:

1)垂直飛行段，要確定轉彎時刻，其選取比較靈活，一般不能小于由發(fā)射條件所確定的規(guī)定時間;

2)跨音速段的攻角應接近零值，這樣可減小氣動特性劇烈變化而造成的影響;

3)導彈飛過跨音速，進入高動壓區(qū)，應考慮殼體溫度和載荷的限制，為了使導彈較輕，又不至于造成折斷的危險，需要限制攻角;

4)程序要簡單易行，引入瞄準段，即此段內(nèi)程序角為常值，這樣，導彈就可以在規(guī)定的各種射程內(nèi)使用;

5)級間分離時，應盡量保證可靠分離和有利于下一級起控.因此，程序選擇時，適當提高第一級分離高度，并在第一級發(fā)動機關機前后的瞬間采用零攻角飛行程序;

6)飛行程序要考慮頭部再入大氣層的限制，如動壓、再入點彈道傾角等.

綜合以上原則，三級導彈的主動段飛行方案為

其中:t1為轉彎開始時刻;t2為重力轉彎開始時刻;t3為一級關機時刻;t4為二級關機時刻;t5為三級關機時刻;φ1f為一級關機時刻俯仰角;φ2f為二級關機時刻俯仰角;θ為彈道傾角;k1為二級俯仰角變化率;k2為三級俯仰角變化率;

式中tm為攻角絕對值等于極值的時刻.

3 計算方法

3.1 諸元解析預報

發(fā)射方位角是彈道導彈的重要射擊諸元，通常根據(jù)發(fā)射點和目標點的位置信息，以球面方位角為初值，進行迭代搜索發(fā)射方位角.本文用球面方位角確定發(fā)射坐標系，選取偏航角來修正彈道，克服地球旋轉帶來的橫向偏差.

對于參數(shù)已定的多級導彈，采用文中的數(shù)學模型和飛行方案，選取不同的發(fā)射點地理緯度B0、球面方位角A0、轉彎段最大負攻角絕對值αm、二級等俯仰角速率k1，三級等俯仰角速率k2，常值偏航角ψ等物理量，進行開環(huán)仿真計算，得到對應的彈道頂點高度H與射程角β.挑選滿足約束條件的彈道并存儲相應的諸元參數(shù)于數(shù)據(jù)庫中，利用線性回歸技術得到攻擊信息與諸元之間的關系式.機動發(fā)射情況下，能夠利用解析關系式迅速獲得滿足約束的諸元近似值.通過研究諸元與H、β的關系，發(fā)現(xiàn)用二次多項式擬合各參數(shù)之間的關系能夠達到較理想結果.

本文以民兵三導彈為例，其具體參數(shù)見表1.

表1 民兵3導彈參數(shù)

以式(1)及如下的式(2)～式(4)為約束條件，按照下列方式建立數(shù)據(jù)庫:球面方位角、緯度間隔選取10°;轉彎段負攻角極值間隔選取1°;二、三級俯仰角變化率間隔選取0.005.對數(shù)據(jù)進行處理可以得到:

分析以上各式，不難發(fā)現(xiàn):諸元與發(fā)射點、目標點的位置沒有直接關系，只與發(fā)射點地理緯度、發(fā)射點和目標點之間的射程、發(fā)射方位角有關.因此，以上解析表達式能夠使民兵三導彈大范圍內(nèi)機動快速發(fā)射，或者快速打擊射程范圍內(nèi)的任意目標.對于其他機動發(fā)射的彈道導彈，可以采用本方法，建立數(shù)據(jù)庫，然后利用線性回歸技術處理數(shù)據(jù)庫數(shù)據(jù)，得到諸元與攻擊信息之間的表達式，在變點發(fā)射時，能夠快速獲得諸元概略值.

3.2 諸元數(shù)值尋優(yōu)

當發(fā)射點為O(λo，Bo)、目標點為T(λT，BT)，發(fā)射點、目標點的地心緯度可用下式計算:

式中:φ為地表一點的地心緯度;B為該點的地理緯度;a、b分別為地球的長半軸和短半軸.

O、T兩點地心角和球面方位角能夠通過球面三角關系得到:

因為遠程彈道導彈的射程不超過20 000 km，所以β在I、II象限內(nèi)取值，A0在4個象限內(nèi)取值，故β、A0的計算公式為

式中:sgn為符號函數(shù);as代表 arcsin;ac代表arccos.

設定期望的彈道高度，利用式(6)～(12)可得一組初值，然后選取參數(shù)X=［αmk1k2ψ］T為變量，以

的值等于零為目標函數(shù).最終以數(shù)值方法對參數(shù)尋優(yōu).式中 G(X)= ［ΔL Δz ΔH］T.ΔL為縱向偏差;Δz為橫向偏差，即r在oz軸投影;ΔH彈道高度與期望高度偏差.

由解非線性最小二乘問題的L-M方法可得到［19］:

其中:i=0，1，…;I為4階單位陣;λi為加權系數(shù)，且0 ＜ λi＜1.

全局收斂的L-M方法計算步驟如下:

步驟1 取 ρ，σ ∈ (0，1)和 λi＞ 0，X0∈R4×n，置 i=0;

步驟2 若G(Xi)=0，停止計算;

步驟3 求解方程組(A+λiI)di=－B得到di;

步驟4 由Armijo搜索求步長，令τi是滿足下面不等式的最小非負整數(shù)τ:

令 τI= ρτ，置 Xi+1=Xi+ εidi+1，i=i+1.

步驟5 目標函數(shù)的最優(yōu)值為0，可以按式λi= ‖G(Xi)‖1+σ，σ ∈［0，1］更新 λk的值，轉步驟2.

4 仿真分析

4.1 仿真條件

以民兵三導彈為例，選擇三級固體彈道導彈，其起飛重量為34 500，彈長為18.26 m，其他重要參數(shù)如表1所示.

4.2 仿真結果及分析

設置四組仿真場景，其中發(fā)射點O和目標點T，發(fā)射點和目標點高程Ho、HT，期望彈道頂點高度H.利用多項式逼近的結果、以多項式結果為初值，式(10)為性能指標，利用如表2所示的攻擊參數(shù)進行L-M迭代，結果如表3所示.

表2 攻擊參數(shù)列表

表3 不同特征分類結果比較

由表2可知，第一組和第二組為相同發(fā)射點攻擊不同目標，第二組和第三組是不同發(fā)射點攻擊相同目標，第三組和第四組是相同發(fā)射點、相同目標點、不同飛行高度，第一組和第三組為不同發(fā)射點、不同目標點.

對比表3的結果，文中方法能夠解決導彈機動發(fā)射情況下彈道諸元快速計算問題.如果在不追求精度的情況下，可以采用本文的擬合函數(shù)進行計算，得到的解析表達式能夠使落點達到精度5 km以內(nèi);若希望得到精確的彈道參數(shù)，則可以擬合公式的解為初值，利用L-M方法進行4～6次迭代計算即可.

以第一組場景為例，給出本文方法得到的仿真曲線，如圖1～圖2所示.由圖1可知，導彈能夠準確擊中目標;由圖2可知，彈道頂點高度與期望高度一致.

圖1 縱向射面內(nèi)曲線

圖2 速度和高度曲線

5 結論

本文基于解析預報與參數(shù)優(yōu)化，對三級固體彈道導彈飛行方案的諸元進行了快速設計.利用多項式函數(shù)擬合了諸元與彈道的方位角、緯度、彈道頂點高度、射程等參數(shù)的函數(shù)關系;對于任意的發(fā)射點、目標點和期望的頂點高度，利用這個函數(shù)關系可以迅速得到諸元的近似值;以彈道諸元為變量，以縱向偏差、橫向偏差、及頂點高度偏差的平方和極小為指標，用L-M方法可以得到諸元的精確值，文中將近似值、L-M方法的解進行了對比分析.本文的結論如下:

1)用多項式對彈道參數(shù)和導彈攻擊參數(shù)進行擬合，能夠得到近似解，在落點精度要求不高的情況下，可以用此解直接作為彈道飛行方案諸元的值，從而實現(xiàn)機動平臺發(fā)射的三級彈道導彈飛行方案參數(shù)快速選取.

2)多項式擬合的方法，即函數(shù)逼近的方法，大量應用于彈道落點預測問題，本文以逼近函數(shù)的解為初值，利用L-M方法迭代尋優(yōu)，節(jié)省大量計算，適合快速計算.

3)本文雖然是針對三級固體彈道導彈進行研究，但本文的方法具有通用性，適合多級固體或液體彈道導彈飛行方案設計.

［1］ROH W R，KIM Y.Trajectory optimization for multistage launch vehicle using time FEM versus collocation［R］.Reston VA:AIAA，1999:AIAA－99－4072.

［2］LU P，PAN B F.Highly constrained optimal launch ascent guidance［J］.Journal of Guidance，Control，and Dynamics.2010，33(2):404 －414.

［3］DUKEMAN G A.Atmospheric ascent guidance for rocket-powered launch vehicles［C］//AIAA Guidance，Navigation，and Control Conference and Exhibit.Monterey:［s.n.］，2002:AIAA 2002 －4559.

［4］LEUNG M S K，CALISE A J.Hybrid approach to nearoptimal launch vehicle guidance［J］.Journal of Guidance，Control，and Dyanmics，1994，17(5):881 －888.

［5］GATH P F，CALISE A J.Optimization of launch vehicle ascent trajectories with path constraints and coast arcs［J］.Journal of Guidance，Control and Dynamics，2001，24(2):296 －304.

［6］CHRISTOPHER L D，WILLIAM W H，ANIL V R.Direct trajectory optimization using a variable low-order adaptive pseudospectral method［J］.Journal of Spacecraft and Rockets，2011，48(3):433 －444.

［7］LISA S，MARC M.General adaptive guidance using nonlinear programming constraint-solving methods［J］.Journal of Guidance，Control，and Dynamics，1993，16(3):517－522.

［8］PARK S Y.Launch vehicle trajectories with a dynamic pressure constraint［J］.Journal of Spacecraft and Rockets，1998，35(6):765 －773.

［9］JOHN T B.Survey of numerical methods for trajectory optimization［J］.Journal of Guidance，Control，and Dynamics，1998，21(2):193 －203.

［10］OLIVER M，MARKUS M.GPS based prediction of the instantaneous impact point for sounding rockets［J］.Aerospace Science and Technology，2002，(6):283 －294.

［11］LI Y C，KIRUBARAJAN T.Trajectory and launch point estimation for ballistic missiles from boost phase LOS measurements［C］//1999 IEEE Aerospace Conference，Aspen，CO，Proceedings.Piscataway:IEEE，1999:425－442.

［12］BENSON A，THORVALDSEN T，RAO Y.Direct trajectory optimization and costate estimation via an orthtogonal collocation method［J］.Journal of Guidance，Control and Dynamics，2006，29(6):1435 －1440.

［13］ROSS I M，F(xiàn)AHROO F.Legendre pseudo-spectral approximations of optimal control problems［C］//Lecture Notes in Control and Information Sciences.Berlin:［s.n.］，2003，295:C101 － C115.

［14］LU P.Trajectory optimization and guidance for an advanced launch system［C］//30th Aerospace Sciences Meeting ＆ Exhibit.Reno，NV:［s.n.］，1992:1 －6.

［15］王明海，李邦杰.彈道導彈彈道設計的一種快速迭代方法［J］.飛行力學，2007，25(2):75 －77.

［16］黃岳，肖凡，張濤.彈道導彈發(fā)射方位角快速計算新方法［J］.四川兵工學報，2009，30(6):51 －52.

［17］劉剛，MURAD Y，何麟書.一種多級火箭飛行方案設計的工程方法［J］.航空學報，2005，26(1):50－53.

［18］賈沛然，陳克俊，何力.遠程火箭彈道學［M］.長沙:國防科技大學出版社，1993:7－8.

［19］馬昌鳳.最優(yōu)化方法及其 Matlab程序設計［M］.北京:科學出版社，2010:87 －96.

A rapid method for flight program design of the ballistic missile launched on mobile platform

WEI Wen-shu，JING Wu-xing，GAO Chang-sheng

(Dept.of Aerospace Engineering，Harbin Institute of Technology，150001 Harbin，China)

Combining analytical prediction and optimization of parameters，a rapid method for basic firing data of flight program of the ballistic missile launched on mobile platform is studied.First，the analytical expressions between the parameters about flight program and attack information are obtained by using polynomial function approximation，by which the initial value of the firing data can be calculated.Second，considering the parameters of flight program as variables and the least square of the placement＇s aberration as the performance index，the Levenberg-Marquardt(L-M)Method is used to obtain the optimal values.Finally，as an example，the proposed method is provided to attacking targets with the LMG-30G.The simulation results show that，because of the less computational effort，the proposed method is suitable for the rapid design of the Big Wide Intratheater Mobility for the missile.When the precision of the placement is not highly demanded，the solution of analytical expressions can be used directly as the value of flight parameters to attack the target.

analytical prediction;launched on mobile platform;Levenberg-Marquardt method;the multi-stage ballistic missile;polynomial function approximation

V412

A

0367－6234(2012)11－0007－06

2011－12－05.

國家自然科學基金資助項目(11172077).

韋文書(1987—)，男，博士研究生;

荊武興(1965—)，男，教授，博士生導師.

韋文書，weiwenshu11@hotmail.com.

(編輯 張 宏)