張意國，趙長見，趙俊鋒，賈生偉

（中國運(yùn)載火箭技術(shù)研究院，北京，100076）

0 引 言

對于大部分運(yùn)載火箭和彈道導(dǎo)彈而言，助推器工作結(jié)束后需要進(jìn)行分離并被拋掉以降低結(jié)構(gòu)重量。對于一般飛行器而言，在子級助推器分離時刻飛行速度不足以達(dá)到第一宇宙速度，使得分離后的助推器在地球引力的作用下將落入地球表面。助推器再入過程由于沒有主動控制力的作用，在重力及非線性氣動力的作用下將出現(xiàn)復(fù)雜的姿態(tài)運(yùn)動。對于助推器再入姿態(tài)運(yùn)動模態(tài)（即無控狀態(tài)下姿態(tài)周期性運(yùn)動規(guī)律）[1]進(jìn)行預(yù)示，有助于對助推器回收系統(tǒng)進(jìn)行設(shè)計(jì)，同時也有助于發(fā)射場安控工作的組織實(shí)施。王景國等[2]人采用了傳統(tǒng)的基于3DOF模型的落點(diǎn)預(yù)測方法對火箭整流罩的殘骸落點(diǎn)進(jìn)行估計(jì)，其中通過大量的試驗(yàn)數(shù)據(jù)進(jìn)行聚類分析結(jié)合工程經(jīng)驗(yàn)對阻力系數(shù)和參考面積這兩個關(guān)鍵彈道參數(shù)進(jìn)行了計(jì)算。陳彬[3]等人對助推器再入運(yùn)動的單通道運(yùn)動特性進(jìn)行了分析，分析了助推器再入過程姿態(tài)運(yùn)動的穩(wěn)定性。

為全面準(zhǔn)確地描述助推器再入自由運(yùn)動過程的姿態(tài)運(yùn)動特性，有必要開展其三通道運(yùn)動模態(tài)的研究。對助推器三通道運(yùn)動模態(tài)的研究建立在助推器再入三通道耦合非線性模型的基礎(chǔ)上，相較于單通道線性化模型[3]，三通道運(yùn)動學(xué)模型能夠提供更高的模型精度，同時也存在不止一組平衡點(diǎn)，對助推器再入動力學(xué)特性的描述可以更準(zhǔn)確。但也因?yàn)榉蔷€性模型強(qiáng)耦合以及強(qiáng)非線性的特點(diǎn)，使得基于小擾動理論的線性系統(tǒng)理論對一些系統(tǒng)特性無法展開準(zhǔn)確地分析。

本文采用了基于演化搜索算法的改進(jìn)郭濤算法，解決高階、強(qiáng)耦合非線性動力學(xué)模型的多平衡點(diǎn)求解問題[4～7]，并在此基礎(chǔ)上對助推器再入過程平衡點(diǎn)的穩(wěn)定性及演變規(guī)律進(jìn)行了分析。

1 助推器再入非線性模型

在再入彈頭和自旋穩(wěn)定飛行器動力學(xué)的研究中，常使用基于經(jīng)典歐拉角描述的運(yùn)動學(xué)模型，以便于更好地描述大姿態(tài)角運(yùn)動下的運(yùn)動規(guī)律[8,9]。對于再入助推器而言，其存在大范圍的姿態(tài)角運(yùn)動，且姿態(tài)角運(yùn)動存在高動態(tài)特性，適合于使用基于經(jīng)典歐拉角及其角速度描述的姿態(tài)動力學(xué)模型來進(jìn)行動力學(xué)分析工作。

經(jīng)典歐拉角的定義如圖1所示，包括進(jìn)動角、章動角及自旋角，是描述助推器體坐標(biāo)系與準(zhǔn)速度坐標(biāo)系之間轉(zhuǎn)換關(guān)系的一組歐拉角，相應(yīng)坐標(biāo)系描述如下：

a）助推器體坐標(biāo)系（Oxbybzb）。

坐標(biāo)系原點(diǎn)O位于質(zhì)心，O xb軸沿助推器縱軸向前為正，Oyb在對稱平面內(nèi)指向上為正，Ozb軸與其余兩軸滿足右手定則。

b）準(zhǔn)速度坐標(biāo)系（Ox5y5z5）。

坐標(biāo)系原點(diǎn)O位于質(zhì)心，O x5軸沿速度方向向前為正，O y5位于包含Ox5的豎直平面內(nèi)向上為正，O z5軸與其余兩軸滿足右手定則。

由準(zhǔn)速度系轉(zhuǎn)到助推器體坐標(biāo)系為1-3-1轉(zhuǎn)序：首先Ox5y5z5坐標(biāo)系沿Ox5旋轉(zhuǎn)ζ（進(jìn)動角）至Ox5y' z'坐標(biāo)系；然后繞Oz'軸旋轉(zhuǎn)δ（章動角）至Oxby'' z'坐標(biāo)系；最后繞Oxb軸旋轉(zhuǎn)λ（自轉(zhuǎn)角）至Oxbybzb坐標(biāo)系。

圖1 經(jīng)典歐拉角定義Fig.1 Classical Euler Angle Definitionζ—進(jìn)動角；δ—章動角；λ—自轉(zhuǎn)角

助推器再入飛行過程沒有控制力，只在重力和氣動力的作用下自由運(yùn)動，而空氣動力一般在彈體系下描述，是與攻角α及側(cè)滑角β相關(guān)的變量，根據(jù)坐標(biāo)轉(zhuǎn)換關(guān)系可得到α，β以ζ，δ，λ表示的表達(dá)式：

考慮到助推器姿態(tài)運(yùn)動相對質(zhì)心運(yùn)動為快時變運(yùn)動，在分析姿態(tài)運(yùn)動時可以認(rèn)為助推器速度的大小及方向保持不變，準(zhǔn)速度系可視為慣性坐標(biāo)系，則再入助推器的角速度（xω，yω，zω）可以用δ、ζ及λ的角速度之和表示：

式中xω，yω，zω為彈體系下歐拉角的3個分量；ζ˙為進(jìn)動角速度；λ˙為自旋角速度；δ˙為進(jìn)動角速度。

對式（2）進(jìn)行求導(dǎo)并進(jìn)行轉(zhuǎn)換，可以得到以歐拉角表示的助推器非線性動力學(xué)模型：

式中xω˙，yω˙，zω˙為彈體系下歐拉角速度bω˙的3個分量；δ˙˙為進(jìn)動角加速度；ζ˙˙為章動角加速度；λ˙˙為自旋角加速度。

可根據(jù)歐拉方程求得，矢量形式為

式中 I為彈體轉(zhuǎn)動慣量矢量；ω為彈體角速度矢量；M為合外力矩矢量。

2 非線性模型平衡點(diǎn)求解

對于助推器再入飛行無控姿態(tài)動力學(xué)系統(tǒng)，為了得到系統(tǒng)在固定參數(shù)下的平衡點(diǎn)，需要將模型中的所有狀態(tài)變量關(guān)于時間的導(dǎo)數(shù)設(shè)為0。以上標(biāo)“P”表示平衡點(diǎn)處的各狀態(tài)變量，則得到平衡點(diǎn)處的運(yùn)動方程：

式中 Ip為彈體平衡點(diǎn)轉(zhuǎn)動慣量矢量； ωp為彈體平衡點(diǎn)角速度矢量； Mp為平衡點(diǎn)合外力矩矢量。

式（6）即為平衡狀態(tài)的運(yùn)動方程，展開可得到平衡方程的標(biāo)量形式為

3 郭濤算法

針對高階強(qiáng)非線性方程組，在理論上不存在解析解，對于該問題的求解一般采用數(shù)值方法。比較常用的數(shù)值方法包括牛頓法和遺傳算法等，但這些算法對初值敏感，求解結(jié)果具有局部性，無法實(shí)現(xiàn)全局多解的求解，甚至在一些非線性較強(qiáng)的方程組求解過程中會出現(xiàn)直接發(fā)散的現(xiàn)象。所以有必要研究對于該類方程組的有效求解方法。

武漢大學(xué)的郭濤提出了一種演化搜索算法，其在求解非線性優(yōu)化問題及多解問題時有獨(dú)特的優(yōu)勢，能夠保證解的全局性，同時收斂速度較快，是一種適合解決本文問題的方法，其特點(diǎn)為：

a）提出了一種基于全局的群體搜索策略，理論上能夠保證解的全局性。

b）提出了種群迭代優(yōu)化策略，采用遺傳算法中父代雜交與隨機(jī)子空間搜索并行策略，保證了在種群優(yōu)化過程中兼顧搜索的遍歷性。

c）提出了末等淘汰策略，種群迭代中每次只淘汰最劣等個體，能夠較好地控制種群的整體性。

郭濤算法同時在理論上保證了全局性及遍歷性，這樣就保證了算法能夠依概率一收斂到全局最優(yōu)解。

該方法的主要求解策略及流程如圖2所示。

圖2 郭濤算法的求解過程Fig.2 Guo Tao Algorithm Solution Process

a）步驟一：問題的描述。

郭濤算法求解的第1步，將常見的非線性多維方程組多解問題抽象為有約束下的最優(yōu)化問題。

設(shè)非線性方程組為

假定該方程組有解，其求解問題可以轉(zhuǎn)化為

其中，S為搜索空間， S ?Rn，當(dāng)fitness(x):S→R稱為目標(biāo)函數(shù)，fitness(x)= 0 時，x為方程組的精確解，通常S是一個n維長方體，即由所確定；Φ為可行點(diǎn)集或可行區(qū)域。

b）步驟二：確定編碼方式。

確定問題的數(shù)學(xué)模型后，需要制定研究對象的編碼方式，對于復(fù)雜的非線性問題，采用浮點(diǎn)數(shù)編碼能夠提高算法求解的準(zhǔn)確程度。

c）步驟三：種群初始化。

d）步驟四：個體評價。

本文的優(yōu)化問題可歸納為連續(xù)優(yōu)化問題，因?yàn)榧s束及解空間在實(shí)數(shù)空間上是處處連續(xù)的，因此針對該問題的求解過程能夠保證較好的局域性及精確程度。

用邏輯函數(shù) b etter(x1, x2)表示測試點(diǎn) x1優(yōu)于測試點(diǎn)：

定義：

用定義的邏輯運(yùn)算函數(shù)，求解種群中的最優(yōu)個體xbest和 最 差 值xworst， 使 得 ： ?x∈P，better(xbest, x) ,better(x, xworst)；如果 b etter(xworst,xbest)，轉(zhuǎn)到步驟八。

e）步驟五：子代優(yōu)化（父體雜交）。

以相同概率選取上述兩種方法中的一種決定搜索空間V，并在V中隨機(jī)選取一點(diǎn)x。

f）步驟六：劣汰策略。

如果better(x, xworst)，則xworst= x 。

g）步驟七：轉(zhuǎn)到步驟四。

h）步驟八：輸出 xbest。

i）步驟九：結(jié)束。

4 助推器再入三通道動力學(xué)特性分析

4.1 平衡點(diǎn)求解

求解平衡狀態(tài)非線性方程組時，氣動力系數(shù)通過攻角α和側(cè)滑角β在（-180°，+180°）范圍內(nèi)插值得到，系統(tǒng)具有定義域范圍大和參數(shù)強(qiáng)非線性的特點(diǎn)（各特征馬赫數(shù)下法向力系數(shù)NC和俯仰力矩系數(shù)MZC的變化規(guī)律見圖 3、圖 4）。傳統(tǒng)解非線性方程組的方法如半解析法和牛頓法只能解決小定義域和參數(shù)定常的問題，不適用于本文的模型。郭濤算法采用變空間搜索及劣勢淘汰策略，對方程組解的搜索具有全局性，不受定義域大小的限制。同時由于郭濤算法基于演化搜索的思想，只要求模型參數(shù)可求，不要求其定?；蚩晌?，使得郭濤算法具有更強(qiáng)的適應(yīng)性，適用于本文的參數(shù)強(qiáng)非線性模型。

圖3 NC 隨攻角變化Fig.3 CN’s Relationship to the Change in Angle of Attack

圖4 MZC 隨攻角變化Fig.4 CMZ’s Relationship to the Change in Angle of Attack

根據(jù)非線性微分方程組的特性和助推器氣動軸對稱的特點(diǎn)，可知若{ζ ˙ ,δ, λ}為方程組的一個解，那么{ζ ˙ ,? δ,π+λ}，{? ζ ˙, δ, ?λ}，{?ζ˙,? δ, π ?λ}也是方程組的解，且在這 4個解處，系統(tǒng)具有相似的姿態(tài)運(yùn)動規(guī)律，可以認(rèn)為它們代表1個解。因此在下文中僅給出一種解的形式作為相同動態(tài)特性解的代表，且如{0 ,π,0} 等顯而易見的解也不再記錄。最后通過郭濤算法仿真計(jì)算得到助推器在不同特征馬赫數(shù)下對應(yīng)的非線性姿態(tài)運(yùn)動的平衡點(diǎn)，見表1。

表1 郭濤算法解不同特征Ma下對應(yīng)非線性模型平衡點(diǎn)Tab.1 Guo Tao Algorithm’s Solutions Depend on Diferrent Ma

續(xù)表1

由表1可知，利用郭濤算法求解出每個特征馬赫數(shù)下助推器姿態(tài)運(yùn)動非線性方程組存在4組解，每組解對應(yīng)的適應(yīng)值函數(shù)均接近于零，表明算法精度較高，結(jié)果可信。

通過經(jīng)典歐拉角的定義可知：ζ˙表征助推器的旋轉(zhuǎn)速度；δ表征助推器縱軸與速度軸夾角，等價于總攻角；λ表征助推器縱平面與總攻角平面夾角。通過表 1對不同特征馬赫數(shù)下的平衡點(diǎn)分析可以得到平衡點(diǎn)隨馬赫數(shù)變化的規(guī)律，見圖5、圖6。由圖5以及圖6的分析結(jié)果知，不同于線性系統(tǒng)，助推器非線性系統(tǒng)存在多組平衡點(diǎn)。不考慮對稱分布的相似平衡點(diǎn)，在再入馬赫數(shù) M a∈ (1 ,9)的情況下分布著 4組平衡點(diǎn)：其中平衡點(diǎn)1對應(yīng)δ=0°；平衡點(diǎn)2對應(yīng)δ≈90°，隨著馬赫數(shù)的增加章動角變化不大，而相應(yīng)的進(jìn)動角速度隨馬赫數(shù)增加而增加；平衡點(diǎn)3、4對應(yīng)δ隨馬赫數(shù)的增加而增加，由δ≈110°向δ≈160°變化。

圖5 平衡點(diǎn)對應(yīng)章動角隨馬赫數(shù)變化規(guī)律Fig.5 Pδ’s Relationship to the Change in Angle of Attack

圖6 平衡點(diǎn)對應(yīng)進(jìn)動角速度隨馬赫數(shù)變化規(guī)律Fig.6 ζ ˙p ’s Relationship to the Change in Angle of Attack

4.2 平衡點(diǎn)穩(wěn)定性分析

助推器再入姿態(tài)動力學(xué)是強(qiáng)耦合高維非線性系統(tǒng)，很難直接根據(jù)李雅普諾夫穩(wěn)定性的定義來判斷其穩(wěn)定性。對于此類非線性系統(tǒng)，李雅普諾夫第一法（或稱間接法）在穩(wěn)定性的分析中具有著廣泛的應(yīng)用。

李雅普諾夫第一法（間接法）：

若x=0是非線性系統(tǒng) x˙ =f( x)的一個平衡點(diǎn)，其中f →Rn是連續(xù)可微的，設(shè)雅克比矩陣為

如果D的所有特征值都滿足 R eλi< 0，則原點(diǎn)是漸近穩(wěn)定的；

如果D至少有一個特征值滿足 R eλi> 0，則原點(diǎn)是不穩(wěn)定的；

如果有某一解的特征根實(shí)部在零點(diǎn)而其余特征根的實(shí)部為負(fù)值，此時則無法判斷該平衡點(diǎn)的穩(wěn)定性。

運(yùn)用Lyapunov的第1種方法求解本文再入非線性系統(tǒng)平衡點(diǎn)的穩(wěn)定性，可分析得到平衡點(diǎn)的穩(wěn)定性隨馬赫數(shù)變化規(guī)律，可得到再入助推器姿態(tài)動力學(xué)平衡點(diǎn)穩(wěn)定性的變化規(guī)律，見圖7。

圖7 δ解曲線的穩(wěn)定性Fig.7 Stability of the δ’s Solution

由圖7可知，在助推器再入飛行絕大部分馬赫數(shù)域內(nèi)，非線性系統(tǒng)具有一組穩(wěn)定平衡點(diǎn)和3組不穩(wěn)定平衡點(diǎn)，這意味著在合適的初始條件下系統(tǒng)會被鎖定在這一平衡狀態(tài)下。而在Ma∈（0.76,0.84）這一馬赫數(shù)域內(nèi)，系統(tǒng)的2組不穩(wěn)定解曲線發(fā)生了性態(tài)變化，演變成穩(wěn)定平衡解曲線，此時系統(tǒng)存在3組穩(wěn)定平衡點(diǎn)，意味著在合適的初始條件下系統(tǒng)可能被鎖定在其中任意一種平衡狀態(tài)下。此處以Ma=0.8的工況為例，將不同初始狀態(tài)對應(yīng)的相軌跡在δ-δ˙平面上進(jìn)行投影，所得結(jié)果如圖8所示。

a）相軌跡在δ-δ˙平面上進(jìn)行投影

續(xù)圖8

圖8 中不同曲線代表的是不同初始條件對應(yīng)的相圖在δ-δ˙平面的投影，從圖 8中可以看到在仿真工況下的吸引子包括2個漸進(jìn)穩(wěn)定的平衡點(diǎn)和一個穩(wěn)定的極限環(huán)，不同的初始條件會使得助推器的姿態(tài)收斂到不同的吸引子，驗(yàn)證了基于演化搜索算法求解的平衡狀態(tài)的正確性。

將不同初始狀態(tài)對應(yīng)的相軌跡在t-ζ˙平面上進(jìn)行投影，所得結(jié)果如圖9所示。

圖9 解軌線在t-ζ˙平面投影Fig.9 Solve for the Projection onto t-ζ˙ Plane

由圖9可知，在Ma=0.8的工況下，在進(jìn)動角速度存在初始偏差的情況下，均能收斂到平衡點(diǎn)對應(yīng)的穩(wěn)定狀態(tài)（平衡點(diǎn)2和4）或在平衡點(diǎn)對應(yīng)穩(wěn)定狀態(tài)附近震蕩（平衡點(diǎn)3），其中平衡點(diǎn)2對應(yīng)的穩(wěn)定進(jìn)動角速度ζ˙為7.006 21，平衡點(diǎn)3和4對應(yīng)ζ˙為0。再次驗(yàn)證了平衡點(diǎn)解算的正確性。

5 結(jié) 論

本文采用基于演化搜索策略的郭濤算法，解決了助推器再入飛行模態(tài)的預(yù)示問題。首先推導(dǎo)了以經(jīng)典歐拉角描述的姿態(tài)運(yùn)動模型，使研究結(jié)論更直觀，模型更簡潔。利用該模型，進(jìn)一步推導(dǎo)了五維姿態(tài)動力學(xué)平衡點(diǎn)方程組，并采用基于演化搜索算法的改進(jìn)郭濤算法解決了非線性方程組的多解問題。進(jìn)一步，利用Lyapunov第一法對平衡點(diǎn)穩(wěn)定性進(jìn)行了分析。最后通過了時域分析的驗(yàn)證。結(jié)果表明在再入過程中的特定馬赫數(shù)域內(nèi)（Ma=（0.76，0.84））存在3組飛行模態(tài)，一組模態(tài)對應(yīng)高速旋轉(zhuǎn)運(yùn)動，其余兩組模態(tài)對應(yīng)姿態(tài)穩(wěn)定運(yùn)動。進(jìn)入該飛行速域時的不同初始狀態(tài)將使助推器姿態(tài)運(yùn)動收斂到不同的飛行模態(tài)。

本文提出了助推器無控飛行模態(tài)的預(yù)示方法，獲得了再入飛行速域內(nèi)助推器可能出現(xiàn)的多個穩(wěn)定飛行模態(tài)，為進(jìn)一步設(shè)計(jì)助推器回收系統(tǒng)或精細(xì)化設(shè)計(jì)安控區(qū)提供了有力支撐。但本文并未對各模態(tài)的穩(wěn)定域進(jìn)行分析，如何解決高維非線性模型穩(wěn)定域的確定將是下一步研究的重點(diǎn)。