基于粒子群優化算法的高空風彈道修正技術

程光輝，荊武興，許元男，葉家銘

（1. 哈爾濱工業大學航天學院，哈爾濱，150000；2. 中國運載火箭技術研究院，北京，100076）

0 引 言

導彈在稠密大氣層中，由于風干擾影響會產生附加氣流攻角。大動壓區的風干擾會導致彈體法向的氣動載荷增大，引發姿態不穩定[1]。在設計標準軌跡時，為了改善飛行條件，需要考慮由風引起的氣流攻角對氣動載荷的影響。

余夢倫針對CZ-2E火箭，提出了火箭高空風彈道修正的方法，介紹了高空風修正的原理和方法以及相關的計算模型，達到了改善火箭飛行環境的目的[2]。宋征宇使用基于彈道修正的被動控制技術，以動壓q和總攻角η的乘積q η×表示法向風載，通過火箭型號飛行的實際數據，表明了風修正技術能夠明顯提高飛行的可靠性[3]。Guanghui Cheng改進了高空風彈道修正技術，針對可重復使用火箭大氣層內返回段，將統計風場的風速均值和風向均值引入到標準軌跡設計中，并將法向風載q η×當作一個優化指標，基于粒子群優化算法設計標準軌跡，提高了可重復使用火箭的飛行可靠性[4]。在太陽帆的軌跡規劃研究領域，傳統的根據太陽風平均特性進行電動太陽帆軌跡修正的方法與高空風彈道修正技術原理相同，Caruso基于太陽風動壓的實時測量值改進了傳統方法，部分抵消了太陽風不確定性對太陽帆的影響[5]。

Eberhart通過模仿鳥群覓食行為，建立了初始的粒子群優化算法[6]。粒子群優化算法由于具有易于收斂和方便操作的特點，已廣泛應用于工程實踐中[7]，用以解決多變量優化問題。然而，初始的粒子群優化算法易陷入局部極小值陷阱。學者們從3個方面改善粒子群優化算法的搜索能力：a）設計新的拓撲結構；b）引入新的策略或技術；c）與其他算法組合使用[8]。李晶基于動態權重改進的粒子群優化算法設計了慣性穩定平臺自適應分數階參數尋優方法，改善了系統的動態性能和抗干擾能力[9]。將本文的軌跡規劃問題轉化為一個多變量優化問題，把較弱的風場當作已知信息，利用改進的粒子群優化算法，可以方便地對該優化問題進行求解。

為降低法向風載，改善導彈的飛行條件，本文分析了水平風場模型對彈體運動的影響，給出飛行時序，提出了一種降低法向風載的高空風彈道修正技術，并基于改進的粒子群優化算法，規劃了多種工況下的標準軌跡，提高了飛行穩定性。

1 動力學模型與風場模型

1.1 動力學模型

針對導彈的軌跡規劃問題，建立質心平動動力學方程：

式中r為質心到參考坐標系原點的位置矢量；v為速度矢量；P為發動機的推力矢量；Q為氣動力矢量；m為彈體質量；g為考慮J2項攝動的重力加速度矢量；Isp為發動機比沖。

本文將式（1）中的動力學模型投影到地面發射慣性系下進行仿真分析。地面發射慣性系的原點是發射點O，Ox指向初始瞄準方向，Oy垂直向上，Oxyz構成右手直角坐標系。在導彈發射后，該參考坐標系在慣性空間保持不變。

1.2 水平風的影響

有風時，空氣流對彈體有附加速度。地面發射坐標系與發射慣性系的定義相近，不同之處在于地面發射坐標系相對于地面靜止不動，因此地面發射坐標系是一個隨地球轉動的動坐標系。將風速矢量在地面發射坐標系下表示為

式中W為按高度插值的風速大小；AW為與射向角定義方式相同的風向；A0為發射點處的射向角。

迎風速度矢量為

式中V為導彈在發射系下的速度矢量。

彈體坐標系的原點為彈體質心Ob，Obxb沿彈體外殼的對稱軸指向彈頭，Obyb在彈體的主對稱平面內且垂直于Obxb，Obxbybzb構成右手直角坐標系。地面發射坐標系到彈體系的轉換矩陣如下：

式中γ為滾轉角，在設計標準軌跡時，設 0γ=；ψ為偏航角；φ為俯仰角。

M1,M2,M3分別為繞x,y,z軸轉動的基變換矩陣。

將式（3）中的迎風速度矢量從發射系投影到彈體系上，可得：

速度坐標系的原點為彈體質心Ob，Obxv指向導彈的速度方向，Obyv在彈體主對稱面內且垂直于Obxv，Obxv y v zv構成右手直角坐標系。考慮風的影響后，速度坐標系到彈體坐標系的轉換矩陣如下：

式中αW為風影響下的攻角；βW為風影響下的側滑角。 將速度矢量由速度系投影到彈體系上，可得：

聯立式（6）和式（8），可得風影響下的側滑角和攻角為

總攻角定義為

法向風載定義為

式中ρ為大氣密度，其值隨高度變化。

1.3 風場模型

1.3.1 基于統計數據的風場模型

相對打擊目標處，發射場附近的高空風數據獲取難度較低，根據過往的風場觀測數據，可以給出某地在某季節按高度節點符合正態分布的風場數據。在射前進行蒙特卡洛打靶仿真飛行實驗時，一般采用如表1所示的符合正態分布的風場數據。

表1 符合正態分布的風場數據 Tab.1 Data of the Wind Field According with Normal Distribution

1.3.2 考慮風切變的風場模型

在極限拉偏仿真時，一般采用考慮風切變的風場模型。平穩風相對于時間和高度變化緩慢。切變風的速度大小隨高度急劇增大，然后急劇減小，一般以三角波的形式，在一定高度范圍內加入切變風。風剖面如圖1所示。

圖1 風剖面 Fig.1 Wind Profile Figure

考慮風切變的水平風場數據如表2所示，平穩風風速稍大于表1中的風速均值。

表2 考慮風切變的風場數據 Tab.2 Data of Wind Field with Shear Wind

2 彈道修正技術

在設計標準軌跡時引入預置風場，選擇合適的指令攻角變化律，基于粒子群優化算法，將風載最大值通過罰函數法轉化為優化指標，可降低法向風載最大值，改善飛行器的抗干擾性能。

2.1 飛行時序

參考文獻[10]，制定表3所示的飛行時序。由于本文關注法向風載對飛行軌跡的影響，采用了以指令攻角和側滑角為主的飛行時序設計思路。

表3 飛行時序 Tab.3 Flight Profile

在一級垂直上升段，指令俯仰角為90°，指令偏航角為0°。

在一級攻角轉彎段，無風干擾的情況下，指令攻角的變化規律為

式中t為飛行時間；αm1為一級攻角轉彎段指令攻角最大值；tm1為指令攻角達到最大值的時刻。

在一級重力轉彎段，無風干擾的情況下，指令攻角與指令側滑角均為0。在設計標準軌跡時，為避免氣動載荷對彈體運動產生較大的影響，本階段需包括跨音速段和動壓最大值點。

在二級攻角轉彎段，無風干擾的情況下，指令攻角變化律如式（13）所示。

式中αm2為二級攻角轉彎段指令攻角最大值；tm31為指令攻角達到最大值的時刻；tm32為指令攻角由最大值開始減小的時刻，即在tm31～tm32時段內，指令攻角保持最大值αm2。

2.2 高空風彈道修正技術

傳統的高空風彈道修正技術要求在跨音速段和最大動壓段的氣流攻角為零，即彈體縱軸與迎風速度矢量共線。為避免指令攻角因為引入風場而發生突變，將風速與一個隨時間變化的“梯形系數”相乘，使得風速由0開始增加至預置風速，經過一段時間，再緩慢減小至0。設在規劃標準軌跡時引入的風速大小為Wsm，各飛行階段實際使用的預置風場的風速大小為Ws，則有：

隨時間變化的“梯形系數”KW（t）見表4。一級重力轉彎段包含跨音速點和最大動壓點，是風干擾影響最大的階段，系數KW（t）設為1。

表4 預置風速的梯形系數 Tab.4 Trapezoidal Coefficients for the Preset Wind Speed

一般基于工程經驗或氣象統計資料給定預置風場的風速與風向數據。本文中預置風場的風速有2種來源：a）統計風場平均風速，見表1，設 smW為0.25倍的風速平均值；b）平穩風，見表2，設 smW為0.25倍的平穩風風速。預置風場的風向有2種來源：a）不隨高度變化的恒定風向，范圍為[0,360)° °；b）統計風場中隨高度變化的風向平均值。對于不同區域的發射需求，需要根據當地氣象統計資料，確定多條用于設計標準軌跡的預置風場數據，在氣象部門給出當前風場預測結果后，選取最為接近的預置風場對應的諸元進行裝訂。

若引入與彈道平面存在夾角的風場，會使得軌跡在偏航方向產生偏差，則需在二級飛行階段引入非零的指令側滑角，對飛行軌跡在偏航方向的偏差進行修正。二級飛行階段指令側滑角的變化規律與式（13）中指令攻角的變化規律相近，為

基于式（9）和無風干擾時指令攻角和指令側滑角的變化規律，風影響下的攻角和側滑角增量為

將式（16）中風攻角和側滑角增量與不考慮風干擾的指令攻角和指令側滑角相結合，使得標準軌跡中指令攻角和指令側滑角按照式（17）進行修正。式（17）為本文提出的“引入風攻角增量進行彈道修正”的方法。

式中Kalf為攻角修正系數。

綜上所述，在使用高空風彈道修正技術設計標準軌跡時，需要優化的變量有：

“引入風攻角增量進行彈道修正”的方法比較依賴于當地氣象統計資料。若將式（11）定義的法向風載引入到優化指標中，在無風工況下，亦能優化出法向風載較小的飛行軌跡，即是“引入法向風載最大值作為優化指標”的方法。

以飛行終點即頭體分離處的狀態為約束，將式（11）定義的法向風載引入到優化指標中，基于罰函數法設計如下優化指標：

式中 下角標f，E分別表示實際和期望的終端狀態；h為飛行高度，km；θ為速度傾角；σ速度偏角，rad；Kh，Kθ，Kσ，KL分別為對應的罰函數系數，Kh=1，Kθ= 57.30，Kσ= 57.30，KL=1 × 10-3。KL越大，降低法向風載最大值的效果越明顯。

式（19）中，除法向風載LW外，其他優化指標可根據實際工程需要進行選取。

基于式（18）和式（19），本文所定義的軌跡優化問題為

優化指標F的值越小，對應的解越優。

3 粒子群優化算法

粒子群優化算法通過粒子群內部信息交換的方式解決優化問題。在式（18）構成的N維解空間中，設有M個粒子。第i個粒子在第t次迭代中的位置和速度矢量為

式中d為區間[1,N]內的正整數。

優化變量的位置和速度邊界如下：

式中

第i個粒子的歷史最優解記作pi，所有粒子的最優解記作pg。第i個粒子的速度矢量vi和位置矢量xi更新公式如下

式中r1和r2為[0,1]內的均勻分布隨機數。

慣性權重w較小時，粒子群優化算法的局部搜索能力較強，w較大時，粒子群優化算法的全局搜索能力較強。慣性權重隨迭代次數線性下降的變化規律[11]如下：

式中 慣性權重最大值wmax= 0.9；慣性權重最小值wmin=0.2；tmax為最大迭代次數。

時變的加速因子c1和c2有利于改善粒子群優化算法的搜索效率，加速因子隨迭代次數線性變化的規律[12]如下

式中 時變加速因子的最大值為cmax=2.5，最小值為cmin=0.2。

為賦予粒子群優化算法跳出局部極小值的能力，當某個粒子靠近當前最優粒子且速度極小時，重新初始化該粒子的速度矢量。設置“重新初始化”規則[4]為

式中ε為一個較小的正數，本文取為 1 × 10-5。

式（27）中歸一化后的距離Dxi和速度Dvi為

即當D xi≤ε且D vi≤ε時，使用均勻分布隨機數重新初始化第i個粒子的速度矢量。

粒子群優化算法的流程如圖2所示。

圖2 粒子群優化算法流程 Fig.2 The Flow Chart of the Particle Swarm Optimization Algorithm

4 仿真分析

本文瞄準頭體分離點終端約束，設計多種工況下的標準軌跡，驗證本文提出的高空風彈道修正技術的有效性。

4.1 仿真參數

發射系原點的大地緯度和大地經度均為0°，發射方位角為270°，海平面高度為0。采用冷發射模式，初始高度為50 m，初始速度為15 m/s。制導周期為 0.1 s。

頭體分離點約束中，高度為147.3556 km，地面發射慣性系下的速度傾角為43.1525°，速度偏角為0°。

飛行參數見表5。

表5 飛行參數 Tab.5 Flight parameters

式（18）中優化變量的范圍見表6。為限制轉彎時的法向過載值，基于表3中的飛行時序，選取了合適的時間變量范圍。

表6 優化變量的范圍 Tab.6 Ranges of Optimized Variables

4.2 仿真結果與分析

為驗證本文提出的改進彈道修正技術在多種工況下的有效性，本節選取了5種仿真場景：a）工況1，預置風場的風速Wsm為0；b）工況2，預置風場取自表，即風速Wsm為0.25倍的風速均值，風向為風向平均值；c）工況3，預置風場取自表，即風速Wsm為0.25倍的風速均值，風向為風向平均值偏轉90 °；d）工況4，預置風場取自表，即風速Wsm為0.25倍的平穩風風速，風向為270 °；e）預置風場取自表，即風速Wsm為0.25倍的平穩風風速，風向為360 deg。

將式（17）中“引入風攻角增量進行彈道修正”作為變量1，將式（19）中“引入法向風載最大值作為優化指標”視作變量2。變量1與變量2共同組成了本文提出的改進高空風彈道修正法。將引入2種變量的操作進行組合，得到如表7所示的標準軌跡分組。

表7 標準軌跡分組 Tab.7 Groups of the Standard Trajectories

工況1（無風工況）的各組法向風載隨時間變化曲線見圖3。圖中法向風載的變化趨勢與總攻角的變化趨勢相近。

圖3 法向風載隨時間變化曲線（工況1） Fig.3 Normal Wind Load Curve Concerning Time (Scenario 1)

工況2的各組法向風載隨時間變化曲線見圖4。由于引入了預置風場，在不使用彈道修正技術的D組，一級重力轉彎段法向風載最大值為1482.6 Pa。在使用了改進彈道修正技術的A組，法向風載最大值減小為981.1 Pa。

圖4 法向風載隨時間變化曲線（工況2） Fig.4 Normal Wind Load Curve Concerning Time (Scenario 2)

工況3的各組法向風載隨時間變化曲線見圖5。

圖5 法向風載隨時間變化曲線（工況3） Fig.5 Normal Wind Load Curve Concerning Time (Scenario 3)

工況4的各組法向風載隨時間變化曲線見圖6。考慮風切變的風場模型中，平穩風的風速稍大于統計風場的平穩風速，因此圖6中一級重力轉彎段的法向風載值大于圖4和圖5。

圖6 法向風載隨時間變化曲線（工況4） Fig.6 Normal Wind Load Curve Concerning Time (Scenario 4)

工況5的各組法向風載隨時間變化曲線見圖7。

圖7 法向風載隨時間變化曲線（工況5） Fig.7 Normal Wind Load Curve Concerning Time (Scenario 5)

各工況下的法向風載最大值見表8。采用“引入風攻角增量進行彈道修正”和“引入法向風載最大值作為優化指標”的A組相對于D組，法向風載最大值分別減小了43.91%、33.83%、55.77%、43.98%和35.10%。僅使用“引入法向風載最大值作為優化指標”的B組相對于D組，法向風載最大值分別減小了31.87%、1.83%、50.98%、16.21%和28.07%。

表8 法向風載最大值 Tab.8 Maximum Values of Normal Wind Load

5 結 論

本文為了降低法向風載，提高導彈的飛行穩定性，提出了一種改進的高空風彈道修正法。以動壓與總攻角的乘積表示法向風載，在規劃標準軌跡時引入預置風場，計算風攻角和風側滑角增量，并引入攻角修正系數，修正了指令攻角和指令側滑角，將法向風載最大值當作優化指標，基于改進粒子群優化算法，規劃標準軌跡。仿真結果表明，相對于不使用高空風彈道修正法的案例，改進的彈道修正法使得法向風載最大值減小了33.83%～55.77%，能夠提高飛行穩定性。