基于粒子群優(yōu)化算法的高空風(fēng)彈道修正技術(shù)

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

（1. 哈爾濱工業(yè)大學(xué)航天學(xué)院，哈爾濱，150000；2. 中國運載火箭技術(shù)研究院，北京，100076）

0 引 言

導(dǎo)彈在稠密大氣層中，由于風(fēng)干擾影響會產(chǎn)生附加氣流攻角。大動壓區(qū)的風(fēng)干擾會導(dǎo)致彈體法向的氣動載荷增大，引發(fā)姿態(tài)不穩(wěn)定[1]。在設(shè)計標(biāo)準(zhǔn)軌跡時，為了改善飛行條件，需要考慮由風(fēng)引起的氣流攻角對氣動載荷的影響。

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

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

為降低法向風(fēng)載，改善導(dǎo)彈的飛行條件，本文分析了水平風(fēng)場模型對彈體運動的影響，給出飛行時序，提出了一種降低法向風(fēng)載的高空風(fēng)彈道修正技術(shù)，并基于改進(jìn)的粒子群優(yōu)化算法，規(guī)劃了多種工況下的標(biāo)準(zhǔn)軌跡，提高了飛行穩(wěn)定性。

1 動力學(xué)模型與風(fēng)場模型

1.1 動力學(xué)模型

針對導(dǎo)彈的軌跡規(guī)劃問題，建立質(zhì)心平動動力學(xué)方程：

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

本文將式（1）中的動力學(xué)模型投影到地面發(fā)射慣性系下進(jìn)行仿真分析。地面發(fā)射慣性系的原點是發(fā)射點O，Ox指向初始瞄準(zhǔn)方向，Oy垂直向上，Oxyz構(gòu)成右手直角坐標(biāo)系。在導(dǎo)彈發(fā)射后，該參考坐標(biāo)系在慣性空間保持不變。

1.2 水平風(fēng)的影響

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

式中W為按高度插值的風(fēng)速大??；AW為與射向角定義方式相同的風(fēng)向；A0為發(fā)射點處的射向角。

迎風(fēng)速度矢量為

式中V為導(dǎo)彈在發(fā)射系下的速度矢量。

彈體坐標(biāo)系的原點為彈體質(zhì)心Ob，Obxb沿彈體外殼的對稱軸指向彈頭，Obyb在彈體的主對稱平面內(nèi)且垂直于Obxb，Obxbybzb構(gòu)成右手直角坐標(biāo)系。地面發(fā)射坐標(biāo)系到彈體系的轉(zhuǎn)換矩陣如下：

式中γ為滾轉(zhuǎn)角，在設(shè)計標(biāo)準(zhǔn)軌跡時，設(shè) 0γ=；ψ為偏航角；φ為俯仰角。

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

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

速度坐標(biāo)系的原點為彈體質(zhì)心Ob，Obxv指向?qū)椀乃俣确较?，Obyv在彈體主對稱面內(nèi)且垂直于Obxv，Obxv y v zv構(gòu)成右手直角坐標(biāo)系?？紤]風(fēng)的影響后，速度坐標(biāo)系到彈體坐標(biāo)系的轉(zhuǎn)換矩陣如下：

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

聯(lián)立式（6）和式（8），可得風(fēng)影響下的側(cè)滑角和攻角為

總攻角定義為

法向風(fēng)載定義為

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

1.3 風(fēng)場模型

1.3.1 基于統(tǒng)計數(shù)據(jù)的風(fēng)場模型

相對打擊目標(biāo)處，發(fā)射場附近的高空風(fēng)數(shù)據(jù)獲取難度較低，根據(jù)過往的風(fēng)場觀測數(shù)據(jù)，可以給出某地在某季節(jié)按高度節(jié)點符合正態(tài)分布的風(fēng)場數(shù)據(jù)。在射前進(jìn)行蒙特卡洛打靶仿真飛行實驗時，一般采用如表1所示的符合正態(tài)分布的風(fēng)場數(shù)據(jù)。

表1 符合正態(tài)分布的風(fēng)場數(shù)據(jù) Tab.1 Data of the Wind Field According with Normal Distribution

1.3.2 考慮風(fēng)切變的風(fēng)場模型

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

圖1 風(fēng)剖面 Fig.1 Wind Profile Figure

考慮風(fēng)切變的水平風(fēng)場數(shù)據(jù)如表2所示，平穩(wěn)風(fēng)風(fēng)速稍大于表1中的風(fēng)速均值。

表2 考慮風(fēng)切變的風(fēng)場數(shù)據(jù) Tab.2 Data of Wind Field with Shear Wind

2 彈道修正技術(shù)

在設(shè)計標(biāo)準(zhǔn)軌跡時引入預(yù)置風(fēng)場，選擇合適的指令攻角變化律，基于粒子群優(yōu)化算法，將風(fēng)載最大值通過罰函數(shù)法轉(zhuǎn)化為優(yōu)化指標(biāo)，可降低法向風(fēng)載最大值，改善飛行器的抗干擾性能。

2.1 飛行時序

參考文獻(xiàn)[10]，制定表3所示的飛行時序。由于本文關(guān)注法向風(fēng)載對飛行軌跡的影響，采用了以指令攻角和側(cè)滑角為主的飛行時序設(shè)計思路。

表3 飛行時序 Tab.3 Flight Profile

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

在一級攻角轉(zhuǎn)彎段，無風(fēng)干擾的情況下，指令攻角的變化規(guī)律為

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

在一級重力轉(zhuǎn)彎段，無風(fēng)干擾的情況下，指令攻角與指令側(cè)滑角均為0。在設(shè)計標(biāo)準(zhǔn)軌跡時，為避免氣動載荷對彈體運動產(chǎn)生較大的影響，本階段需包括跨音速段和動壓最大值點。

在二級攻角轉(zhuǎn)彎段，無風(fēng)干擾的情況下，指令攻角變化律如式（13）所示。

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

2.2 高空風(fēng)彈道修正技術(shù)

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

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

表4 預(yù)置風(fēng)速的梯形系數(shù) Tab.4 Trapezoidal Coefficients for the Preset Wind Speed

一般基于工程經(jīng)驗或氣象統(tǒng)計資料給定預(yù)置風(fēng)場的風(fēng)速與風(fēng)向數(shù)據(jù)。本文中預(yù)置風(fēng)場的風(fēng)速有2種來源：a）統(tǒng)計風(fēng)場平均風(fēng)速，見表1，設(shè) smW為0.25倍的風(fēng)速平均值；b）平穩(wěn)風(fēng)，見表2，設(shè) smW為0.25倍的平穩(wěn)風(fēng)風(fēng)速。預(yù)置風(fēng)場的風(fēng)向有2種來源：a）不隨高度變化的恒定風(fēng)向，范圍為[0,360)° °；b）統(tǒng)計風(fēng)場中隨高度變化的風(fēng)向平均值。對于不同區(qū)域的發(fā)射需求，需要根據(jù)當(dāng)?shù)貧庀蠼y(tǒng)計資料，確定多條用于設(shè)計標(biāo)準(zhǔn)軌跡的預(yù)置風(fēng)場數(shù)據(jù)，在氣象部門給出當(dāng)前風(fēng)場預(yù)測結(jié)果后，選取最為接近的預(yù)置風(fēng)場對應(yīng)的諸元進(jìn)行裝訂。

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

基于式（9）和無風(fēng)干擾時指令攻角和指令側(cè)滑角的變化規(guī)律，風(fēng)影響下的攻角和側(cè)滑角增量為

將式（16）中風(fēng)攻角和側(cè)滑角增量與不考慮風(fēng)干擾的指令攻角和指令側(cè)滑角相結(jié)合，使得標(biāo)準(zhǔn)軌跡中指令攻角和指令側(cè)滑角按照式（17）進(jìn)行修正。式（17）為本文提出的“引入風(fēng)攻角增量進(jìn)行彈道修正”的方法。

式中Kalf為攻角修正系數(shù)。

綜上所述，在使用高空風(fēng)彈道修正技術(shù)設(shè)計標(biāo)準(zhǔn)軌跡時，需要優(yōu)化的變量有：

“引入風(fēng)攻角增量進(jìn)行彈道修正”的方法比較依賴于當(dāng)?shù)貧庀蠼y(tǒng)計資料。若將式（11）定義的法向風(fēng)載引入到優(yōu)化指標(biāo)中，在無風(fēng)工況下，亦能優(yōu)化出法向風(fēng)載較小的飛行軌跡，即是“引入法向風(fēng)載最大值作為優(yōu)化指標(biāo)”的方法。

以飛行終點即頭體分離處的狀態(tài)為約束，將式（11）定義的法向風(fēng)載引入到優(yōu)化指標(biāo)中，基于罰函數(shù)法設(shè)計如下優(yōu)化指標(biāo)：

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

式（19）中，除法向風(fēng)載LW外，其他優(yōu)化指標(biāo)可根據(jù)實際工程需要進(jìn)行選取。

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

優(yōu)化指標(biāo)F的值越小，對應(yīng)的解越優(yōu)。

3 粒子群優(yōu)化算法

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

式中d為區(qū)間[1,N]內(nèi)的正整數(shù)。

優(yōu)化變量的位置和速度邊界如下：

式中

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

式中r1和r2為[0,1]內(nèi)的均勻分布隨機數(shù)。

慣性權(quán)重w較小時，粒子群優(yōu)化算法的局部搜索能力較強，w較大時，粒子群優(yōu)化算法的全局搜索能力較強。慣性權(quán)重隨迭代次數(shù)線性下降的變化規(guī)律[11]如下：

式中 慣性權(quán)重最大值wmax= 0.9；慣性權(quán)重最小值wmin=0.2；tmax為最大迭代次數(shù)。

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

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

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

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

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

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

粒子群優(yōu)化算法的流程如圖2所示。

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

4 仿真分析

本文瞄準(zhǔn)頭體分離點終端約束，設(shè)計多種工況下的標(biāo)準(zhǔn)軌跡，驗證本文提出的高空風(fēng)彈道修正技術(shù)的有效性。

4.1 仿真參數(shù)

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

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

飛行參數(shù)見表5。

表5 飛行參數(shù) Tab.5 Flight parameters

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

表6 優(yōu)化變量的范圍 Tab.6 Ranges of Optimized Variables

4.2 仿真結(jié)果與分析

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

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

表7 標(biāo)準(zhǔn)軌跡分組 Tab.7 Groups of the Standard Trajectories

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

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

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

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

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

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

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

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

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

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

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

表8 法向風(fēng)載最大值 Tab.8 Maximum Values of Normal Wind Load

5 結(jié) 論

本文為了降低法向風(fēng)載，提高導(dǎo)彈的飛行穩(wěn)定性，提出了一種改進(jìn)的高空風(fēng)彈道修正法。以動壓與總攻角的乘積表示法向風(fēng)載，在規(guī)劃標(biāo)準(zhǔn)軌跡時引入預(yù)置風(fēng)場，計算風(fēng)攻角和風(fēng)側(cè)滑角增量，并引入攻角修正系數(shù)，修正了指令攻角和指令側(cè)滑角，將法向風(fēng)載最大值當(dāng)作優(yōu)化指標(biāo)，基于改進(jìn)粒子群優(yōu)化算法，規(guī)劃標(biāo)準(zhǔn)軌跡。仿真結(jié)果表明，相對于不使用高空風(fēng)彈道修正法的案例，改進(jìn)的彈道修正法使得法向風(fēng)載最大值減小了33.83%～55.77%，能夠提高飛行穩(wěn)定性。