郭瑋林，鮮 勇，張大巧，凌王輝

（火箭軍工程大學 七系，西安 710025）

由于高超聲速飛行器助推段彈道設計約束條件多，動力學過程復雜，終端入軌條件苛刻，就目前研究情況而言，實現助推段彈道優(yōu)化設計本身就是一個難點問題[1-3]，若要實現助推段彈道的快速高精度計算滿足高超聲速飛行器機動發(fā)射要求更是一個巨大的挑戰(zhàn)。

文獻[4]針對高超聲速飛行器助推段彈道設計方面進行了研究，但彈道計算速度和精度有待提高。文獻[5]聯合間接法和直接法可在3 min內計算生成高超聲速飛行器助推段最優(yōu)軌跡，但彈道計算耗時較長未達到快速計算的要求。文獻[6]主要是運用牛頓迭代法實現對助推段彈道的快速計算，但牛頓法對初值敏感且易陷入局部最優(yōu)。而L-M算法作為牛頓法的改進，融合了牛頓迭代局部快速收斂和梯度下降法全局迭代迅速的優(yōu)點，通過引入阻尼因子，能夠有效避免雅可比矩陣奇異時迭代不收斂問題[7]。

BP神經網絡是一種多層前饋網絡，網絡結構簡單實時性好，只要有足夠的隱層和隱節(jié)點，就可以逼近任意的非線性映射關系[8]?；贐P神經網絡建立發(fā)射點及終端入軌點狀態(tài)量與彈道參數的映射關系，可以實現彈道諸元參數快速預測，為后續(xù)L-M方法數值尋優(yōu)計算提供彈道參數初值。

因此本文采用基于BP神經網絡和L-M方法的聯合算法實現助推段彈道快速計算。首先綜合考慮各項約束條件設計了助推段飛行程序和彈道優(yōu)化模型；然后利用 BP神經網絡預測計算諸元參數初值，并且建立了基于BP神經網絡和L-M方法的聯合數值計算模型；最后通過仿真實驗驗證了聯合算法能夠有效地實現助推段彈道的快速高精度計算。

1 助推段彈道設計模型

1.1 助推段飛行程序設計

以三級固體火箭作為高超聲速飛行器助推段的運載器，假定地球是以角速度ω自轉的橢球，且飛行器運動姿態(tài)能夠達到瞬時平衡。

高超聲速飛行器的助推段關機點高度相比彈道導彈較低，所以一級飛行程序采取二次攻角轉彎的設計方式，即飛行器分別在跨音速前和跨音速后進行一次攻角轉彎，跨音速時零攻角飛行?？v向俯仰程序設計如下：

式中：φcx表示飛行程序角；0～t1為垂直上升段；t1～t2為跨音速前攻角轉彎段；t2～t3為跨音速飛行段，攻角為零；t3～t4為跨音速后攻角轉彎段；t4～tk1為一二級分離前的等程序瞄準段；a1、a2、a～1、a～2為攻角設計參數；θ為彈道傾角，ωz為地球自轉角速度分量。

二級和三級飛行程序采用分段線性化方法進行設計，具體形式如公式(3)，其中：為二級等程序飛行段；為二級轉彎段，k1為該段的程序轉彎角速率；為級間分離段；為三級轉彎段，為該段的程序轉彎角速率；為三級等程序飛行段，tk3為三級關機點。

1.2 助推段彈道優(yōu)化模型

1）助推段約束條件

由于受到彈體結構、執(zhí)行機構性能的限制，必須限制飛行器轉彎角速率、過載和動壓，以滿足執(zhí)行機構性能約束和穩(wěn)定飛行的要求。為減小被敵方紅外偵測系統識別的可能性，飛行器低彈道飛行時，高度須低于 90 km；同時，為順利實現入軌飛行，助推段關機點還必須滿足終端高度、速度和當地彈道傾角約束。具體約束形式如下所示：

2）優(yōu)化變量的確定

所以最終高超聲速飛行器助推段彈道參數確定為3個優(yōu)化變量，具體如下：

式中，Tcz為一級發(fā)動機垂直上升段工作時間。

3）目標函數

為滿足高超聲速飛行器終端入軌條件，采取罰函數方式將終端約束轉換為目標函數，即：

式中，a、b、c為權系數。

2 BP神經網絡參數初值計算模型

以某發(fā)射區(qū)域為例進行仿真分析（假定地球為橢球，所以發(fā)射點經度改變對終端參數無影響）。為實現高超聲速飛行器助推段彈道快速計算，采用 BP神經網絡算法，通過訓練樣本，建立發(fā)射點緯度、高程、射向及關機點高度、速度、彈道傾角與助推段彈道諸元參數的函數關系式，即建立以下映射關系：

2.1 BP神經網絡結構

根據BP神經網絡基本原理，若確定了網絡層數、各層節(jié)點數、傳遞函數、初始權系數和閥值、學習方法等，BP網絡結構也就得以確定。

由式(7)可知，神經網絡輸入層節(jié)點為6個，輸出層節(jié)點為3個，所以高超聲速飛行器彈道模型為六輸入三輸出的神經網絡結構。

根據Kosmogorov定理，一個具有合適參數和結構的3層BP神經網絡可以完成任意的n維到m維的映射，所以可以選取3層BP神經網絡結構進行仿真計算。隱層傳遞函數選用logsig函數，輸出層傳遞函數選用purelin線性函數，學習方法選用L-M方法。

采用網絡結構增長型方法進行仿真計算，確定隱層節(jié)點數量。通過測試誤差的比較，選取隱層節(jié)點數為12個。最終確定的BP神經網絡結構如圖1所示。

圖1 BP神經網絡結構Fig.1 BP neural network structure

2.2 BP神經網絡樣本生成及訓練

隨機選取該發(fā)射區(qū)域發(fā)射點緯度、高程和射向信息，利用優(yōu)化算法求解滿足飛行約束和入軌終端約束的高超聲速飛行器助推段彈道，從而得出5000組訓練樣本。由于粒子群算法具有收斂速度快、易實現的特點[9-11]，所以采用粒子群智能算法對彈道諸元參數進行優(yōu)化求解，生成訓練樣本。

利用圖1神經網絡結構對樣本進行訓練，得到BP神經網絡輸入層至隱層的權值W1和閾值b1以及隱層至輸出層的權值W2和閾值b2。因此可以得到發(fā)射點位置射向、關機點高度、速度、彈道傾角與彈道諸元參數的函數關系式。根據 BP神經網絡結構和各層傳遞函數可得：

由式(8)和式(9)推導BP神經網絡函數解析式的最終形式如下：

由式(10)可以快速預測計算發(fā)射區(qū)域任意發(fā)射點位相應助推段彈道諸元參數。

2.3 BP神經網絡仿真計算與誤差分析

訓練神經網絡的根本目的是確保訓練好的網絡模型對非訓練樣本具有好的泛化能力。為測試以上神經網絡模型性能，再次隨機選取發(fā)射區(qū)域500組測試樣本，測試樣本發(fā)射點位為高超聲速飛行器助推段終端要求的標準高度、標準速度和標準彈道傾角為將代入BP神經網絡解析式(10)即可計算得到彈道諸元參數，將該諸元參數代入高超聲速飛行器標準彈道得到 BP神經網絡預測的助推段終端高度、速度和彈道傾角與標準終端高度、速度和彈道傾角對比得到 BP神經網絡終端高度計算偏差、速度計算偏差和彈道傾角計算偏差，如圖2～4所示。基于BP神經網絡方法的樣本仿真統計結果如表1所示。

圖2 BP神經網絡終端高度計算偏差Fig.2 Calculated deviation of terminal height based on BP neural network

圖3 BP神經網絡終端速度計算偏差Fig.3 Calculated deviation of terminal velocity based on BP neural network

圖4 BP神經網絡終端彈道傾角計算偏差Fig.4 Calculated deviation of terminal inclination angle based on BP neural network

表1 基于BP神經網絡的樣本仿真結果Tab.1 Simulation results of samples based on BP neural network

由圖2至圖4和表1的仿真結果可知，基于BP神經網絡求解的高超聲速飛行器助推段彈道具有一定的精度，但入軌精度并不高，特別是終端高度最大入軌偏差達到386.92 m，高度平均偏差也有121.56 m，所以僅僅基于 BP神經網絡計算飛行器助推段彈道精度是不夠的。

BP神經網絡解算的助推段彈道終端參數偏差較大，不能滿足高精度入軌要求，但可以利用 BP神經網絡解析式快速預測計算諸元參數，將其作為后續(xù)數值尋優(yōu)計算的彈道參數初值，加快數值計算的收斂速度，以求得滿足終端參數精度要求的助推段彈道。

3 BP神經網絡與L-M算法聯合計算模型建立

牛頓法具有收斂速度快、精度高等優(yōu)點，但對初值敏感且易陷入局部最優(yōu)，而L-M算法作為牛頓法的改進，通過引入阻尼因子λ，能夠有效避免雅可比矩陣奇異或病態(tài)時迭代不收斂情況發(fā)生。因此可以利用BP神經網絡快速預測彈道諸元參數初值，然后采用L-M算法進一步數值尋優(yōu)計算，以解決機動發(fā)射條件下高超聲速飛行器彈道快速計算的問題。

下面研究建立基于BP神經網絡與L-M算法的聯合數值尋優(yōu)計算模型。

設助推段終端高度偏差Δhf、速度偏差Δvf和彈道傾角偏差ΔfΘ是諸元參數Tcz、k1、k2的三維偏差函數f(x)，即

式中，x為狀態(tài)變量，將f(x)在x處進行一階泰勒展開可得如下公式：

式中，x0為彈道諸元參數初值，Δx為迭代步長，Q為助推段終端參數偏差對諸元參數的偏導數矩陣，可采用彈道求差法計算求解。Q矩陣形式如下：

其中，彈道諸元參數初值x0的計算采用BP神經網絡解析式(10)求解，即得：

助推段終端偏差函數的最小二乘模型為

將式(12)代入式(15)得到聯合算法數值計算模型如下：

由文獻[12]可知式(16)的解可表示為

式中，λ為L-M算法的阻尼因子，λI為阻尼項。

在迭代計算過程中：當助推段終端參數偏差較大時，可以選取較大的λ，使得L-M算法具有梯度下降法下降量大、迭代迅速的特點；當助推段終端參數偏差較小時，可以選取較小的λ，使得L-M算法具有牛頓法二階收斂的特點，從而快速收斂到最優(yōu)解。

Yamashita和Fukushima在文獻[13]證明了若阻尼因子時，那么在一個比非奇異性要弱的局部誤差有界條件下，L-M算法將保持二次收斂。因此針對阻尼因子λ的選取，可采用如下自適應策略：

式中，1δ、2δ為終端參數偏差量范數的閥值，且0＜δ2＜δ1。

基于聯合算法計算模型，根據式(17)進行多次迭代數值尋優(yōu)計算，最終解得滿足助推段終端入軌精度要求的彈道諸元參數。

基于BP神經網絡與L-M方法的聯合算法計算流程如圖5所示，聯合算法具體計算步驟如下：

1）將機動發(fā)射點緯度、高程和射向信息與終端入軌點高度、速度和彈道傾角參數輸入到訓練好的BP神經網絡解析式中，計算得到彈道諸元初始值

2）再將諸元參數初值代入L-M算法數值計算模型中進行迭代尋優(yōu)計算，生成新的彈道諸元參數

4）若助推段終端參數滿足飛行器終端入軌精度要求或迭代達到最大次數，則退出計算，輸出結果；否則，返回步驟2）。

圖5 聯合算法計算流程圖Fig.5 Flow chart of the joint algorithm

4 仿真實驗結果與分析

隨機選取發(fā)射區(qū)域 500組發(fā)射點位作為仿真樣本，采用基于BP神經網絡和L-M方法的聯合算法對高超聲速飛行器助推段彈道進行仿真計算得到樣本助推段終端參數偏差值?；诼摵纤惴ǖ慕K端高度、速度和彈道傾角計算偏差分別如圖6～8所示，助推段彈道計算耗時Ts仿真如圖 9所示，基于聯合算法的樣本仿真統計結果如表2所示。本次仿真實驗計算在普通微機上實現。

圖6 聯合算法終端高度計算偏差Fig.6 Calculated deviation of terminal height based on the joint algorithm

圖7 聯合算法終端速度計算偏差Fig.7 Calculated deviation of terminal velocity based on the joint algorithm

圖8 聯合算法終端彈道傾角計算偏差Fig.8 Calculated deviation of terminal inclination angle based on the joint algorithm

圖9 彈道計算耗時仿真圖Fig.9 Simulation of trajectory calculation time

由圖6至圖9和表2仿真結果可知，經仿真計算生成的500條助推段彈道均滿足終端約束要求，并且終端高度偏差在2 m以內，終端速度偏差在0.1 m/s以內，終端彈道傾角偏差在0.01°以內，保證高超聲速飛行器以較高的精度入軌。同時助推段彈道計算最大耗時為 2.84 s，平均耗時0.8211 s，彈道計算效率高、速度快。所以基于BP神經網絡和L-M方法的聯合算法能夠較好地滿足機動發(fā)射條件下高超聲速飛行器助推段彈道快速高精度計算的要求。

表2 基于聯合算法的樣本仿真結果Tab.2 Simulation results of samples based on joint algorithm

5 結 論

本文研究了機動發(fā)射條件下的高超聲速飛行器助推段彈道快速計算問題?；?BP神經網絡建立了發(fā)射點及終端入軌點狀態(tài)量與助推段彈道諸元參數的函數解析式，并采用 BP神經網絡方法預測彈道參數初值，最終推導出基于BP神經網絡和L-M算法的聯合數值尋優(yōu)計算模型。由仿真結果可知，采用聯合算法計算得到的助推段彈道終端高度、速度和彈道傾角偏差分別在2 m、0.1 m/s和0.01°以內，彈道計算最大耗時不超過3.0 s。結果表明基于BP神經網絡和L-M方法的聯合算法對于機動發(fā)射高超聲速飛行器助推段彈道快速高精度計算具有較強的魯棒性和適應性。本文通過引入神經網絡算法，利用BP神經網絡為L-M算法數值尋優(yōu)計算提供參數初值，提高了數值算法的智能化水平，較好地解決了算法計算速度的問題。

參考文獻（References）：

[1]Wu G, Liu L, Wang Y, et al. Ascent trajectory optimization of hypersonic vehicle based on improved Particle Swarm algorithm[C]//IEEE Chinese Automation Congress. 2016: 115-120.

[2]Dalle D J, Torrez S M, Driscoll J F, et al. Minimum-fuel ascent of a hypersonic vehicle using surrogate optimization[J]. Journal of Aircraft, 2014, 51(6): 1973-1986.

[3]Feng L, Liu L, Wang Y. Trajectories optimization of hypersonic vehicle based on a hybrid optimization algorithm of PSO and SQP[C]//IEEE Control and Decision Conference. 2015: 4518-4522.

[4]任京濤. 助推滑翔導彈上升段多終端約束彈道設計及制導方法研究[D]. 哈爾濱: 哈爾濱工業(yè)大學碩士論文, 2013.Ren J T. Research on reconfigurable control system design of heavy launch vehicle[D]. Harbin: Harbin Institute of Technology Master Degree Dissertation,2013.

[5]崔乃剛, 傅瑜, 盧寶剛. 助推-滑翔飛行器助推段最優(yōu)軌跡快速生成方法[J]. 彈道學報, 2013,25(2): 33-38.Cui N G, Fu Y, Lu B G. A fast generation method of optional trajectory of boost phase for boost-glide vehicle [J]. Journal of ballistics, 2013, 25(2): 33-38.

[6]劉開封, 陳穎, 何念念, 等. 基于正交實驗與牛頓方法的滑翔導彈助推彈道設計[J]. 導彈與航天運載技術, 2014, 334(4): 50-54.Liu K F, Chen Y, He N N, et al. Boost trajectory design of boost-glide missile based on orthogonal experimental design and Newton's iterative method[J]. Missiles and Space Vehicles, 2014, 334(4): 50-54.

[7]Ma C F, Jiang L H. Some research on Levenberg-Marquardt method for the nonlinear equations[J].Applied Mathematics and Computation, 2007, 184:1032-1040.

[8]鮮勇, 李少朋, 李邦杰. 基于BP神經網絡的固體導彈耗盡關機姿態(tài)調制方法的研究[J]. 兵工學報,2015, 36 (4): 668-673.Xian Y, Li S P, Li B J. An approach to attitude angle modulation of solid missiles under the condition of depleted shutdown based on BP neural network[J]. Acta Armamentarii, 2015, 36(4): 668-673.

[9]Ran M P, Wang Q. Spacecraft rendezvous trajectory optimization method based on EPSO[J]. Journal of Astronautics, 2013, 34(9): 1195-1201.

[10]Pontani M, Conway B A. Optimal finite-thrust rendezvous trajectories found via particle swarm algorithm[J]. Journal of Spacecraft and Rockets,2014, 94(1): 434-445.

[11]Zhang W, Huang X, Gao X Z, et al. Multi-objective longitudinal trajectory optimization for hypersonic reentry glide vehicle based on PSO algorithm[C]//IEEE Control Conference. 2015: 2350-2356.

[12]崔乃剛, 張亮, 韋常柱, 等. 可重復使用運載器大姿態(tài)機動自抗擾控制[J]. 中國慣性技術學報,2017, 25(3): 387-394.Cui N G, Zhang L, Wei C Z, et al. Active disturbance rejection control for reusable launch vehicle with large attitude maneuver[J]. Journal of Chinese Inertial Technology, 2017, 25(3): 387-394.

[13]Yamashita N, Fukushima M. On the rate of convergence of the Levenberg-Marquardt method[J]. Computing (Supplement 15), 2001: 239-249.