周馳， 李穎暉， 鄭無計， 武朋瑋， 董澤洪

(空軍工程大學(xué)航空工程學(xué)院, 西安 710038)

結(jié)冰不僅會惡化飛機的動力學(xué)特性，也會嚴重影響駕駛員正常操縱，因而結(jié)冰飛行事故持續(xù)發(fā)生[1]。駕駛員在結(jié)冰情況下仍繼續(xù)保持常規(guī)操縱，飛行狀態(tài)極易超出安全邊界[2]。NASA對1976—1994年間的飛行事故進行了數(shù)據(jù)統(tǒng)計，其中由結(jié)冰引發(fā)的飛行事故就有16起，共造成139人死亡[3]。結(jié)冰不僅使飛機升力減小，阻力增大，失速迎角也會對應(yīng)減小，這導(dǎo)致飛機更易失速。1994年，美國鷹航公司的一架ATR72-212飛機因機翼嚴重結(jié)冰在印第安納州墜毀，當時飛行迎角僅為5°，遠遠低于飛控系統(tǒng)給出的限制值18.1°，但事故卻仍因飛機失速而發(fā)生[4]。

對于剛體飛機而言，一旦控制結(jié)構(gòu)確定下來，整個飛控系統(tǒng)也被固定。同時飛行指令也存在一個限制，只要機組下達的飛行指令超出這個限制，飛機將很難回到原有的平衡狀態(tài)。當飛機的動力學(xué)特性發(fā)生變化，就可能出現(xiàn)抖動甚至失速。文獻[5]指出通過分岔分析方法能夠清晰地對環(huán)境影響下的飛機動力學(xué)狀態(tài)變化進行觀測，此時駕駛員可通過操縱將飛行狀態(tài)調(diào)整到新的穩(wěn)定平衡狀態(tài)，這也是一種常用的邊界恢復(fù)方法。分岔分析方法目前正被廣泛應(yīng)用于研究飛機的動力學(xué)特性變化，并用于指導(dǎo)飛行操縱[6-7]。Engelbrecht等[8]指出飛行系統(tǒng)的穩(wěn)定性與分岔點有關(guān)。操縱過度，分岔就會出現(xiàn)，系統(tǒng)就可能失穩(wěn)。Harry和Jean[9]利用分岔分析方法對飛機動力學(xué)特性進行了研究，但其只是研究了對飛行指令的限制，雖然提出了當飛機失穩(wěn)后將飛行狀態(tài)調(diào)整到某一穩(wěn)定平衡點(Stable Equilibrium Point，SEP)的策略，但是該狀態(tài)點附近穩(wěn)定邊界的范圍卻沒有定義，如果找到的穩(wěn)定平衡點附近的穩(wěn)定邊界較小，即使將飛行狀態(tài)調(diào)整到該穩(wěn)定狀態(tài)，一個輕微的擾動都能使飛行狀態(tài)重新失穩(wěn)。為提升飛行安全，應(yīng)該在此基礎(chǔ)上找到可靠的邊界確定方法。

近年來，國內(nèi)外對安全邊界的確定方法已經(jīng)進行了大量研究。傳統(tǒng)的邊界確定方法是通過數(shù)值計算或飛行試驗，給出迎角、側(cè)滑角、俯仰角、飛行速度、過載等飛行參量的限制值。但由于傳統(tǒng)方法無法將環(huán)境因素考慮在內(nèi)，因此就不再適用于結(jié)冰條件安全邊界確定。近年來，吸引域(Region of Attraction，ROA)方法[10]被用于安全邊界的確定，但是ROA方法主要是基于Lyapunov能量函數(shù)[11]確定邊界，精確性依賴于所選擇的Lyapunov函數(shù)，因此保守性較強。Weekly和Tinka[12]提出用可達集方法確定飛行安全邊界，并將最大可控不變集作為安全邊界，但可達集方法基于最優(yōu)控制，在不同控制下可能確定的邊界也不同，保守性較強，同時該方法不是基于穩(wěn)定平衡狀態(tài)確定邊界，因此在指導(dǎo)飛行操縱上存在一定難度。

本文提出了利用微分流形理論對結(jié)冰條件下的飛行系統(tǒng)進行非線性穩(wěn)定域確定，該方法在穩(wěn)定域的確定上主要是基于穩(wěn)定平衡點。穩(wěn)定域的確定是通過計算穩(wěn)定平衡點穩(wěn)定邊界上的不穩(wěn)定平衡點(Unstable Equilibrium Point，UEP)的穩(wěn)定流形的并集得到，并將其作為飛行系統(tǒng)的安全邊界[13-14]。只要將系統(tǒng)的飛行狀態(tài)控制在穩(wěn)定域內(nèi)，飛行系統(tǒng)就不會有失穩(wěn)風險，而在穩(wěn)定域外的點則會振蕩發(fā)散[15]。本文首先利用分岔分析方法對不同結(jié)冰程度條件下的飛機的動力學(xué)特性變化進行觀測，給出了不同條件下駕駛員的可操縱安全范圍，并得到了飛機的穩(wěn)定平衡點；然后在選定的穩(wěn)定平衡點處基于微分流形理論進行安全邊界確定。通過將2種方法相結(jié)合，實現(xiàn)了結(jié)冰條件下，即使飛行狀態(tài)超出原安全邊界，仍然可以通過操縱指導(dǎo)使其到達新的安全邊界。

1 結(jié)冰飛機非線性模型

1.1 飛機縱向通道動力學(xué)模型

結(jié)冰飛機非線性動力學(xué)方程可表示成如下微分方程形式：

(1)

式中：狀態(tài)變量x中包含飛行迎角α、俯仰角θ以及俯仰角速度q；η為結(jié)冰因子，反映結(jié)冰程度。

建立的飛機三維非線性縱向通道動力學(xué)模型如下：

(2)

式中：m為飛機質(zhì)量；L為飛機升力；V為飛行速度；M為飛機縱向通道力矩；Iy為飛機縱向通道轉(zhuǎn)動慣量；g為重力加速度。

升力和力矩可以表示為

(3)

式中：Tx和Tz分別為推力在機體軸上的分量，T為飛機推力；Fx和Fz分別為氣動力在機體軸上的分量；ρ為大氣密度；S為飛機參考面積；c為平均空氣動力學(xué)弦長；Cx、Cz、Cm為力及力矩系數(shù)，具體見參考文獻[9]。

1.2 結(jié)冰飛機多項式系數(shù)模型

對飛行數(shù)據(jù)進行多項式參數(shù)擬合能夠精確地得到氣動數(shù)據(jù)隨飛行迎角α、俯仰角速度q以及升降舵偏角δe變化的多項式形式的解析表達式[16]，如下：

(4)

式中：多項式系數(shù)xi(i=1, 2, 3, 4) 和zi(i=1, 2, 3, 4)為結(jié)冰因子η的函數(shù)。

可通過如下表達式體現(xiàn)結(jié)冰程度對飛機氣動特性的影響[2]：

C(A)iced=(1+ηKC(A))C(A)

(5)

式中：KC(A)表示飛行參數(shù)的變化；C(A)為結(jié)冰前力或力矩系數(shù)。結(jié)冰因子η是一個時變量，隨結(jié)冰程度變化，變化范圍從0到ηmax。η=0表示飛機不受結(jié)冰影響，η=ηmax表示飛機完全結(jié)冰。另外，式(5)所用的結(jié)冰程度影響模型是由Bragg提出的初步模型，雖然模型較為簡單，但在一定程度上能夠反映結(jié)冰的趨勢，根據(jù)公開發(fā)表的文獻[4,17-18]，可知該方法得到的結(jié)果對于本文使用的GTM飛機具有較強的通用性，且比較精確可靠。

1.3 增穩(wěn)控制器設(shè)計

為了提升飛機的穩(wěn)定性，對飛行系統(tǒng)進行反饋增穩(wěn)控制。選取飛行迎角α、俯仰角θ以及俯仰角速度q作為控制變量，對整個飛行系統(tǒng)進行增穩(wěn)設(shè)計，具體結(jié)構(gòu)如圖1所示。圖中：θc為駕駛員操縱指令輸入；Kα、Kq為相應(yīng)的反饋控制律參數(shù)；θe為指令輸入誤差；KP為指令放大系數(shù)。

圖1 飛機縱向通道控制器結(jié)構(gòu)Fig.1 Structure for longitudinal channel controller of aircraft

2 安全性分析方法

本文的主要研究內(nèi)容是將微分流形理論與分岔分析方法相結(jié)合，并共同用于指導(dǎo)結(jié)冰條件下的駕駛員操縱。下面對這2種方法的基礎(chǔ)理論進行介紹。

2.1 分岔分析方法

分岔分析方法可用于如下的非線性微分系統(tǒng)：

(6)

式中：f(·)為非線性微分系統(tǒng)的函數(shù)；x為狀態(tài)變量；u為控制輸入。分岔分析方法可以用于展示全維飛行狀態(tài)，但對于本文研究的飛機縱向通道而言，主要展示(V，α，θ，q)這4維。在分岔分析中，分岔點的形成主要是由于系統(tǒng)平衡點的穩(wěn)定性質(zhì)發(fā)生變化所導(dǎo)致。當出現(xiàn)結(jié)冰情況，飛機的動力學(xué)系統(tǒng)將被破壞，這也會對平衡點的穩(wěn)定性質(zhì)造成影響。伴隨著結(jié)冰程度的不斷增加，原本穩(wěn)定的平衡點也可能變成不穩(wěn)定平衡點，這樣就形成了分岔點。同樣的，當駕駛員操縱過猛也會造成動力學(xué)系統(tǒng)失穩(wěn)，特別是對于結(jié)冰飛機，駕駛員的操縱會變得更加敏感。因此，本文分別將操縱指令輸入δe及結(jié)冰因子η作為分岔參數(shù)進行研究。

2.2 微分流形理論

穩(wěn)定性是衡量飛行安全的一項主要指標，只要能將飛行狀態(tài)保持在穩(wěn)定域內(nèi)，飛機就不會有失穩(wěn)的風險[19-20]。本文提出的基于微分流形理論方法求穩(wěn)定域是通過以穩(wěn)定平衡點穩(wěn)定邊界上的不穩(wěn)定平衡點的穩(wěn)定流形作為安全邊界，該方法能精確地刻畫出飛行系統(tǒng)的穩(wěn)定域。

2.2.1 基本理論

非線性動力學(xué)系統(tǒng)可表示為以下形式：

(7)

對于系統(tǒng)式(7)，若存在點x0，使f(x0)=0，則x0稱為系統(tǒng)的平衡點。若x0處的Jacobian矩陣的特征值實部全不為零，則x0為雙曲平衡點。若特征值實部全為負，則x0為穩(wěn)定平衡點；反之，若其中存在實部大于零的特征值，則x0是不穩(wěn)定平衡點。對于不穩(wěn)定平衡點，實部小于零的部分的特征向量{v1,v2,…,vl}張成穩(wěn)定的特征空間Es，大于零的部分對應(yīng)的特征向量{vl+1,vl+2,…,vn}張成不穩(wěn)定的特征空間Eu，全空間表示為Rn=Es+Eu。

雙曲平衡點x0的穩(wěn)定流形Ws(x0)和不穩(wěn)定流形Wu(x0)可表示如下：

(8)

式中：φ(t,x)為系統(tǒng)式(7)始于x的解。

2.2.2 非線性系穩(wěn)定域具體構(gòu)造方法

本文的穩(wěn)定域構(gòu)造主要是基于軌道弧長法[21]，流程如圖2所示。具體步驟如下：

步驟1確定初始圓。利用不穩(wěn)定平衡點xi所對應(yīng)的穩(wěn)定特征向量組成特征子空間(即切平面)，然后在切平面上以xi為中心，r為半徑作初始圓，并在圓上均勻取N個點{p1,1,p1,2,…,p1,N}。

步驟2求解軌線。以點集{p1，1,p1，2,…,p1，N}為初始點，分別對系統(tǒng)進行反時間積分，當軌線長度達到設(shè)定值L0時停止，得到第一代軌線{T1,1,T1,2,…,T1,N}，各條軌線的終點記為{p2,1,p2,2,…,p2,N}。

步驟3軌線疏密程度判斷。首先檢查各條軌線間的距離，若2條軌線間的距離大于Dmax，則在這2條軌線的初始點間重新插入一個新初始點。若2條軌線間的距離小于Dmin，則刪除其中一條軌線的初始點。最后重新調(diào)整初始點集，返回步驟2，直到滿足要求，進入步驟4。

圖2 軌道弧長法構(gòu)造穩(wěn)定邊界流程Fig.2 Flowchart of stability boundary with track arc length method

步驟4以步驟3中確定的軌線終點為二代初始點，重復(fù)步驟2和步驟3，迭代達到Zmax次時停止。

步驟5連接相鄰代的軌線。將始于相鄰初始點的軌線連接，形成邊界面。

步驟6將穩(wěn)定邊界上所有的不穩(wěn)定平衡點所確定的穩(wěn)定流形組合在一起形成系統(tǒng)的穩(wěn)定域。

3 結(jié)冰飛機案例分析

本文以NASA的GTM為研究對象，考慮增穩(wěn)控制結(jié)構(gòu)的影響，利用1.1節(jié)給出的結(jié)冰飛機縱向通道模型，對飛行系統(tǒng)的動力學(xué)特性及非線性穩(wěn)定域進行分析。飛機的相關(guān)動力學(xué)參數(shù)如下：飛機質(zhì)量m=90 703 kg，大氣密度ρ=0.538 4 kg/m3，重力加速度g=9.806 3 m/s2，飛行高度H=6 000 m，俯仰軸的轉(zhuǎn)動慣量Iy=14 679 000 kg/m2，機翼參考面積S=185.25 m2，平均空氣動力學(xué)弦長c=4.87 m，飛機具體的氣動力及力矩參數(shù)見參考文獻[9]。

3.1 分岔分析

本文中分岔分析方法主要是用于研究結(jié)冰后飛機的動力學(xué)特性變化，通過展示飛行狀態(tài)平衡點的形式揭示。為了便于分析，研究過程中各個控制參數(shù)保持不變，而駕駛員操縱指令作為變化量，其中平衡點的計算是通過計算飛機的動力學(xué)方程得到。仿真分別針對不同結(jié)冰條件給出(V，α，θ，q)的分岔分析圖，如圖3～圖5所示。

圖3 無結(jié)冰(η=0)分岔圖Fig.3 Bifurcation diagram without icing (η=0)

圖4 結(jié)冰因子η=0.1分岔圖Fig.4 Bifurcation diagram for icing factor η=0.1

圖5 結(jié)冰因子η=0.3分岔圖Fig.5 Bifurcation diagram for icing factor η=0.3

圖3為飛機不結(jié)冰(η=0)時的分岔圖，其中紅色和藍色部分分別表示穩(wěn)定平衡點和不穩(wěn)定平衡點，指令輸入δe的范圍為-10°～10°之間變化?？梢钥闯觯旓w機不結(jié)冰時指令輸入δe在-2.2°～10°之間都存在穩(wěn)定的平衡點，當操縱指令值低于2.2°，飛機就會存在失速的風險。文獻[9]研究了一旦飛機失穩(wěn)，通常是將駕駛桿控制到對應(yīng)的指令為δe=0°時的穩(wěn)定平衡狀態(tài)附近進行邊界恢復(fù)，這也是一種最直接有效的方法。

圖4為飛機輕度結(jié)冰(η=0.1)時的分岔圖。通過與圖3對比可知，當飛機遭受結(jié)冰影響，飛機動力學(xué)模型結(jié)構(gòu)也發(fā)生變化，從而平衡點分布也被改變。相對于不結(jié)冰情形，穩(wěn)定平衡點的范圍變化不大，駕駛員指令δe在-1.8°～10°之間都存在穩(wěn)定的平衡點，因此輕度結(jié)冰對飛機的動力學(xué)特性的影響較小。

圖5為飛機重度結(jié)冰(η=0.3)時的分岔圖。可以看出，隨著結(jié)冰程度增加，在同一操縱指令下不僅平衡點的分布發(fā)生巨大變化，穩(wěn)定平衡點的穩(wěn)定性質(zhì)也發(fā)生較大變化，原本穩(wěn)定的平衡點也因結(jié)冰的影響變得不穩(wěn)定。同時結(jié)冰程度的加劇也使駕駛員的可操縱范圍相應(yīng)減小，駕駛員指令δe只有在0°～6.8°范圍內(nèi)存在穩(wěn)定平衡點，從而造成駕駛員的操縱變得更加敏感，輕易推拉桿都容易導(dǎo)致飛機失穩(wěn)。

3.2 微分流形理論確定安全邊界

3.2.1 精確性驗證

本節(jié)中蒙特卡羅方法被用于對微分流形方法進行精確性驗證。蒙特卡羅方法作為一種傳統(tǒng)方法，已經(jīng)廣泛運用于飛行仿真驗證。其原理是：首先在工作點(任意取的某一穩(wěn)定平衡點)所確定的狀態(tài)空間內(nèi)大量取初始點，然后在初始點處進行正向積分，最后將能夠收斂于工作點的初始點保留下來構(gòu)成穩(wěn)定域。該方法的特點是初始點取的越多，精度越高。但是初始點取得過多，計算所消耗的時間也越長。下面分別利用微分流形方法及蒙特卡羅方法計算穩(wěn)定域，仿真環(huán)境在Intel (R) Core (TM) i7-4790 處理器下進行，其中微分流形方法選擇的圓弧向外增長長度(步長)為0.01，向外擴展圈數(shù)為10圈。蒙特卡羅方法在滿足精度要求的前提下，將其初始點數(shù)控制在50×50×50。分別用蒙特卡羅方法和微分流形方法所計算的穩(wěn)定域如圖6所示。

圖6 流形方法精確性驗證Fig.6 Accuracy verification for manifold method

圖6中，藍色點為蒙特卡羅方法構(gòu)成的穩(wěn)定域，彩色曲面為基于微分流形理論確定的安全邊界。可以看出，用微分流形方法確定的穩(wěn)定域能夠與蒙特卡羅方法高度吻合，因此證明了基于微分流形理論構(gòu)造穩(wěn)定域具有較高的精度。另外在仿真時間上微分流形方法耗時較短僅需2 min，而蒙特卡羅方法雖然同樣能夠得到精確的穩(wěn)定域，但耗時卻超過了1 h。

3.2.2 結(jié)冰飛機安全邊界確定及操縱保護

本節(jié)主要研究基于微分流形理論確定結(jié)冰飛機非線性穩(wěn)定域，并將其作為飛機的安全邊界。對于式(2)的模型，從3.1節(jié)的分岔分析圖中選駕駛桿操縱值為δe=5°時的穩(wěn)定平衡點(0.01，-0.35,0)作為工作點在(α，θ，q)3維進行穩(wěn)定域確定。針對不同的結(jié)冰程度，所確定的穩(wěn)定域如圖7所示。

圖7中紅色、綠色、黑色區(qū)域分別代表無結(jié)冰(η=0)、輕度結(jié)冰(η=0.1)以及重度結(jié)冰(η=0.3)時飛機的穩(wěn)定域。可以看出，當飛機只是輕度結(jié)冰時，穩(wěn)定域變化不大，飛行狀態(tài)仍然容易維持在穩(wěn)定域內(nèi)。但當飛機發(fā)生重度結(jié)冰，穩(wěn)定域已經(jīng)嚴重收縮，原先在穩(wěn)定域內(nèi)的平衡狀態(tài)可能會超出穩(wěn)定域，這給飛機帶來失穩(wěn)的風險，下面給出重度結(jié)冰時的飛機操縱保護策略。

當飛機遭遇重度結(jié)冰影響，可能會導(dǎo)致飛行狀態(tài)超出原有安全邊界，由分岔分析可知，可通過改變駕駛桿，使飛行狀態(tài)進入新的安全邊界，如圖8所示。圖中黑色區(qū)域為重度結(jié)冰(η=0.3)，駕駛桿操縱量δe=5°時飛機的穩(wěn)定域。A(0.17，0.08，0.06) rad為原先穩(wěn)定域內(nèi)的一點，但當飛機重度結(jié)冰，飛行狀態(tài)已經(jīng)超出原穩(wěn)定域，如圖中黑色區(qū)域所示。在點A處進行時域仿真，解軌線為圖中藍色曲線1所示，從解軌線可看出，最終飛行狀態(tài)振蕩發(fā)散。通過拉桿將操縱桿拉回至δe=0°處，此時狀態(tài)點A進入了新的穩(wěn)定域內(nèi)，如圖中紅色區(qū)域所示。在A點處進行時域仿真，解軌線為圖中綠色曲線2，從解軌線可看出，軌線最后將會收斂于新的穩(wěn)定平衡點，最終飛行狀態(tài)達到穩(wěn)定。

圖7 不同結(jié)冰程度穩(wěn)定域Fig.7 Stability region for different degrees of icing

圖8 不同操縱下穩(wěn)定域Fig.8 Stability region for different manipulations

4 結(jié) 論

本文以GTM為案例飛機，將飛機氣動參數(shù)進行多項式擬合并結(jié)合結(jié)冰因子模型，建立了結(jié)冰條件下的飛機縱向通道動力學(xué)模型。提出了基于微分流形理論確定飛機的安全邊界，然后結(jié)合分岔分析方法對結(jié)冰后飛機的穩(wěn)定性質(zhì)變化進行了研究，最后對于重度結(jié)冰情形給出了駕駛員操縱策略。主要結(jié)論如下：

1) 分岔分析可以得到飛機在不同結(jié)冰條件下的平衡點，并給出安全操縱范圍，而微分流形理論正是基于穩(wěn)定平衡點確定安全邊界，相比在狀態(tài)空間內(nèi)進行大量搜索尋點，簡化了工作量的同時能夠更加清晰直觀地對飛行狀態(tài)變化進行觀測。

2) 不同結(jié)冰程度對于飛機的操穩(wěn)性能影響不同。對于輕度結(jié)冰，穩(wěn)定域向內(nèi)收縮程度低于10%，此時在同一穩(wěn)定平衡狀態(tài)下飛機仍具有較強抗干擾能力，同時對于駕駛員的可操縱范圍影響較小，操縱范圍相較不結(jié)冰情形只減少5%左右。但對于嚴重結(jié)冰，穩(wěn)定域向內(nèi)收縮程度已經(jīng)超過50%，飛機的抗擾動能力也嚴重降低，甚至原先的穩(wěn)定平衡狀態(tài)將變得不穩(wěn)定，同時駕駛員的可操縱范圍也減少了45%左右，這給飛行安全帶來巨大風險。

3) 微分流形理論在不同程度結(jié)冰條件下都能夠精確地確定飛機高維穩(wěn)定域，將其與分岔分析方法相結(jié)合能夠?qū)崿F(xiàn)即使駕駛員的飛行狀態(tài)在穩(wěn)定域以外，仍能夠通過指導(dǎo)操縱使其到達一個新的穩(wěn)定域，最終使飛行狀態(tài)達到穩(wěn)定。

另外針對結(jié)冰因子模型的局限性，將進一步完善飛機結(jié)冰后的動力學(xué)非線性模型建立，并針對不同冰型、不同部位結(jié)冰以及不同飛機結(jié)冰進行研究，同時對飛機結(jié)冰保護控制方法進行更加深入研究。