基于動力學邊界的結冰飛機安全預警方法

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

空軍工程大學 航空工程學院，西安 710038

飛機結冰破壞飛機的氣動特性，使升力降低、阻力增加，同時使飛機失速迎角大幅度下降，惡化飛機的飛行性能，導致傳統的邊界保護系統和操縱方法失效，因此對飛行安全造成嚴重的威脅[1]，另外由于飛機結冰導致的飛行事故具有強破壞性[2]和高發性[3]，因此引起了廣大航空領域專家學者對飛機結冰后飛行安全問題的關注。國外對飛機結冰機理[4-5]、氣動特性影響[6-8]、防/除冰方法[9]、動力學特性[10-13]等相關理論的研究較早且涉獵內容較為廣泛，而國內針對結冰飛機相關理論的研究起步相對較晚且致力于對飛機防/除冰系統[14]及致災機理方面[15-17]的研究，針對飛機結冰后動力學特性的研究相對較少，然而這部分內容卻對提高結冰飛機的飛行安全十分重要。

Bragg教授團隊[10, 18]指出從規避結冰條件和帶冰飛行兩方面入手可有效減少或防御結冰引起的飛行事故，其中通過規避結冰條件來防止飛機結冰事故是非常理想的方法，但在一定程度上是無法實現的[18-19]，因此研究飛機帶冰情況下的動力學特性是至關重要的。另外，飛機結冰使飛機飛行品質下降，導致駕駛員無法準確評估飛機的動力學特性，因此也經常發生由于駕駛員誤操縱而引發的飛行事故。例如：2002年臺灣復興航空的ATR72-200飛機由于機翼嚴重結冰而失速墜海；2006年在安徽地區同樣因為機翼嚴重結冰引發重大飛行事故。

為解決飛機帶冰飛行問題，首先必須獲得結冰飛機準確的氣動特性。鑒于此，2000年NASA[20]關于飛機結冰的報告中詳細地分析了飛機結冰對氣動特性的影響，與此同時Bragg團隊[10]提出了結冰飛機氣動參數的結冰因子影響模型，為后續對結冰飛機的動力學特性及邊界保護的研究奠定了堅實的理論基礎。在此基礎上，2008—2010年，美國田納西州大學(University of Tennessee Space Institute，UTSI)與Bihrle應用研究(Bihrle Applied Research, BAR)技術公司合作，構建了結冰污染邊界保護系統[21-22](Icing Contamination Envelope Protection system， ICEPro)，提出了結冰邊界告警與保護方法。

另外，綜合分析上述邊界保護系統及相關理論的文獻可知[13, 18, 23-25]，構建精確、有效的邊界保護系統，必須有可靠的飛行風險預警預測方法。然而現有飛機裝備的邊界保護系統很難準確評估結冰情況下的潛在危險，例如1994年的ATR72事故[26]，該飛機由于結冰導致迎角為5°時就發生了滾轉異常現象并最終導致飛行事故。而文獻[27]提出了一種基于微分流形理論確定的可考慮飛行狀態耦合情況的動力學邊界，并初步確定了飛機結冰情況的動力學邊界。因此，本文在微分流形理論確定的動力學邊界基礎上詳細分析了結冰對該動力學邊界的影響，提出基于動力學邊界進行結冰飛機飛行風險評估的方法，并以此為基礎結合飛行仿真相關理論構建結冰飛機的飛行仿真訓練系統；以著陸為測試科目對系統的正確性進行驗證，最后利用該系統初步對駕駛員的結冰風險感知能力進行訓練。結果表明該系統可對駕駛員進行結冰應對策略訓練，可為結冰飛機飛行訓練模擬器的構建提供理論支撐。

1 基于流形理論的結冰飛機動力學邊界

穩定域作為微分動力學系統穩定與不穩定的界限，充分地反映出微分動力學系統空間狀態的運動趨勢，而飛機本質為復雜微分動力學系統，其飛行包線在一定程度上也具有分割安全飛行與不安全飛行的性質，因此可利用微分動力學系統的穩定域作為飛機的某種飛行包線，因穩定域考慮了飛機的動力學特性且考慮國際上類似的概念，取名為動力學邊界(即國際上所用的Dynamic Envelope)。

1.1 流形理論

流形理論是根據微分系統的空間拓撲結構確定系統穩定平衡點吸引區的一種方法，它利用直接積分方法得到邊界上不穩定平衡點的穩定流形作為穩定邊界的一種精確估計[1]。針對穩定平衡點穩定邊界上具有多個不穩定平衡點的情況，通過依次計算邊界不穩定平衡點上的穩定流形，最終以所有穩定流形的并集所包圍的區域作為該穩定平衡點的穩定域。因此，基于流形理論計算穩定域即飛機動力學邊界的難點在于邊界上不穩定平衡點穩定流形的計算。現今針對不穩定平衡點穩定流形的計算較多[28]，且優缺點各異，其中逆軌跡法效率較高但在一定程度上限制了算法的精度，本文在此基礎上以MATLAB仿真平臺為依托對逆軌跡法終點進行校正，并利用并行算法快速計算飛機的動力學邊界，提高計算精度的同時提高了計算效率，具體計算方法如下。

考慮如下形式的非線性系統

(1)

式中：

x=[x1x2…xN]T

(2)

1) 通過求解非線性方程組f(x)=0，得到系統的平衡點，并通過f(x)在平衡點的Jacobian矩陣特征值的情況判斷平衡點的穩定性，并確定不穩定平衡點的型別。

2) 判斷不穩定平衡點是否位于穩定平衡點的邊界上。

3) 利用穩定邊界上的1型不穩定平衡點的穩定流形確定系統的穩定邊界。

下面以穩定流形的計算為例說明本文的二維流形的計算方法(對于流形而言，一維為線，二維為面，三維為體，因本文主要研究不同速度情況下，迎角、俯仰角及俯仰角速度的動力學邊界，為空間曲面，故此處應為二維流形)：

1) 計算向量場f(x)的雅克比矩陣D，得到特征值及其對應的特征向量分別為ε1、ε2、ε3和v1、v2、v3。并設ε1>0、ε2、ε3為小于零的實數或實部小于零的共軛復數。

2) 計算流形初始平面的法向量ξ：ε2、ε3為小于零實數，此時ξ=v2×v3；ε2、ε3為共軛復根，此時ξ=(v2+v3)×[(v2-v3)i]，i為復數單位。

3) 確定穩定流形初始點。以不穩定平衡點(Unstable Equilibrium Point,UEP)為原點在法向量ξ確定的平面上取圓形，并在圓周上均勻取點作為穩定流形的初始點。

4) 計算圓形軌線并進行校正。利用反時間系統計算每個初始點解軌線，并在達到特定長度l后停止本次計算，并對軌線終點進行校正，以便得到光滑的圓形軌線。校正方法如下：

① 計算終點校正時間tadd

(3)

② 計算終點的校正點xnew

xnew=xend+taddf(xend)

(4)

5) 計算相鄰解軌線終點的距離d，并作如下判斷(Δ為圓形軌線的相鄰點的距離限制，其選擇主要是為了避免軌線點過少引起的計算誤差以及軌線點過多而引起的計算資源浪費，一般取每單位軌跡圓100～200點為宜)：

① 如果0.5Δ

② 如果d>1.5Δ，則在相鄰兩點之間插入round(d/Δ)-1個點，其中round(x)函數為通過四舍五入的方法取整；

③ 如果d<0.5Δ，則舍去該點。

注：第4)、5)步采用并行算法，可同時計算圓形軌跡上的點和距離。

6) 通過上述計算得到完整的圓形軌線，并以該圓形軌線為初始軌線重復第4)、5)步過程，直到計算總長度Nl達到預定距離。

將上述N+1個圓形軌線連接成面，即為所求的二維穩定流形。

1.2 干凈飛機動力學邊界確定及算法驗證

以NASA用于研究飛機失控的GTM(Generic Transport Model)[29-30]飛機的縱平面運動模型為對象，驗證上述算法的效率及精度，最終本文所采用的動力學模型及參數為

(5)

式中：

(6)

δe=kαΔα-kqq+kθ(θref-θ)

(7)

式中：反饋參數kα=-2，kq=-1，kθ=-3；Δα=α0-α為控制誤差；θref為控制指令。

以大飛機打開15°襟翼、飛行速度Vt=85 m/s、高度H=400 m的下滑飛行任務為例，驗證上述方法的精度和效率。本文研究的飛行狀態具體如下：

(8)

式中：δth為推力系數。

計算該狀態下的平衡點分布情況，其中平衡點分布結果如表1所示，表中Type表示不穩定平衡點的型別。

表1 系統平衡點的分布Table 1 Distribution of equilibrium points for system

在不穩定平衡點添加小擾動得到動態響應曲線，并根據最終的收斂情況確定穩定平衡點邊界上的不穩定平衡點，其中不穩定平衡點的動態響應結果如圖1所示。

在UEP1添加擾動后所有的飛行狀態都是發散的，而在UEP2、UEP3和UEP4添加擾動后存在飛行狀態最終收斂于穩定平衡點的情況。結合表1最終確定邊界上的1型不穩定平衡點為UEP2和UEP4，利用1.1節方法確定該飛行情況下的動力學邊界(圖2)，并與Monte Carlo方法確定的動力學邊界對比驗證本文方法的精確性(圖3)。

圖1 動力學邊界上的不穩定平衡點確定Fig.1 Determination of UEPs on dynamic envelope

圖2 大飛機下滑階段動力學邊界Fig.2 Dynamic envelope for large aircraft during gliding phase

圖3 動力學邊界精確性驗證Fig.3 Accuracy verification of dynamic envelope

如圖2所示，為利用流形理論確定的飛機動力學邊界。其中紅色曲面為1型不穩定平衡點UEP2的穩定流形(Manifold 1)；藍色曲面為1型不穩定平衡點UEP4的穩定流形(Manifold 2)，兩者的交線為2型不穩定平衡點UEP3的穩定流形且在交線附近光滑過度。綜上說明本文方法適用于處理類似飛機這類多平衡點系統。

如圖3所示，為本文流形方法與Monte Carlo方法確定動力學邊界的對比圖。其中Monte Carlo法是在狀態空間內取大量飛行狀態點，并通過飛行狀態點的動態響應是否收斂于穩定平衡點來判斷飛行狀態的安全性，并認為收斂于穩定平衡點的飛行狀態是安全的，且所有安全飛行狀態所構成的邊界為本文所確定的動力學邊界。通過對比可知，所有安全的飛行狀態均在本文確定的動力學邊界內，且安全的飛行狀態邊界與本文方法確定的動力學邊界完美重合，因此可以證明本文方法具有較高的精度。另外根據動力學邊界內飛行狀態動態響應的收斂性可說明該動力學邊界的安全特性，即動力學邊界為動力學意義上飛行狀態安全與不安全的分界線，這在一定程度上為在線飛行風險評估與安全預警提供理論支撐。

最后，利用逆軌跡和Monte Carlo法在Intel(R) Core(TM) i7-4690 CUP、主頻3.6GHz、8 GB內存的臺式機上進行相同工況的動力學邊界仿真計算，通過對比(如表2所示)，最終可說明本文方法具有較高的計算效率。

表2 計算時間Table 2 Calculation time s

1.3 結冰飛機動力學邊界確定

飛機結冰導致氣動特性被破壞，這也是導致飛行事故的原因之一，因此研究結冰飛機的動力學特性必須有準確合理的結冰飛機氣動參數。2000年Bragg教授[10]提出的基于結冰因子的結冰飛機氣動參數確定方法，在一定程度上解決了該問題，且被廣大學者所接受，因此本文研究采用結冰因子方法確定結冰飛機的氣動模型，其表達式為

CA,iced=(1+ηfice)CA

(9)

式中：CA和CA,iced分別為結冰前、后飛機氣動導數值；fice為結冰系數，反映CA對結冰的敏感性，對于給定的飛機為常值；η為結冰因子。

利用1.2節的控制方法對帶冰飛機在高度H=400 m、打開15°襟翼，且以飛行速度Vt=85 m/s、 下滑角為2.5°穩定下滑為研究對象，確定結冰因子分別為0.1和0.3時的動力學邊界，見圖4。

如圖4所示，飛機結冰后動力學邊界收縮，尤其是結冰較為嚴重的情況下(η=0.3)，動力學邊界收縮較為明顯且邊界狀態的耦合較為嚴重，因此在這種情況下進行操縱將有較高的風險，具體的風險評估方法見第2節。

圖4 結冰程度對飛機動力學邊界的影響Fig.4 Influence of different icing degree on aircraft dynamic envelope

2 結冰飛機飛行安全預警方法

2.1 飛行風險評估及其量化方法

對飛機系統而言，本文使用的動力學邊界為該飛機微分動力系統的穩定域，因此根據穩定域的特點及上述相關論述可知，動力學邊界內部的飛行狀態將在有限時間內收斂于穩定平衡點，即飛機的工作點。且飛機飛行狀態離邊界越近，飛機越容易因外部擾動或駕駛員誤操縱而引發飛機失控現象的發生。因此可利用飛行狀態距動力學邊界的距離表征此時飛行狀態的飛行風險，即距離動力學邊界越近則相對飛行風險越高，離動力學邊界越遠則相對飛行風險越低。考慮到動力邊界的復雜性且無隱式或顯式表達，本文采用如下方法對飛行風險進行量化，并稱該方法為基于動力學邊界的飛行風險評估方法，其幾何概念圖如圖5所示。圖中點S為穩定平衡點(用xs表示)；點F為飛機實時的飛行狀態(用xf表示)；點E為SF與動力學邊界的交點(用xe表示)，其中飛行情況①為安全飛行情況，此時點E為狀態空間內射線SF與動力學邊界的交點，點E位于SF延長線上；飛行情況②為危險飛行情況，此時點E為狀態空間內SF與動力學邊界的交點，點E位于點S和點F中間。

圖5 安全系數的幾何概念Fig.5 Geometric definition of safety factor

(10)

式中:Sf為飛行狀態的安全系數(Safety factor)，用于量化飛機飛行過程中的風險等級。

圖6 飛行風險等級劃分Fig.6 Grades of flight risk

2.2 相對于傳統邊界預警方法的優越性

傳統飛行安全預警方法通常對迎角進行限制，限制值一般為略小于失速迎角的可用迎角，但該極限值不能隨飛機所處環境的變化而變化，也無法反映出飛行狀態之間的耦合關系，因此具有一定的局限性。而本文的安全預警方法可保證迎角在不超出限制的情況下，提前發出危險告警，尤其是飛機飛行狀態耦合較為明顯的情況下效果更為明顯，因此，在一定程度上可提高飛行安全，具體原因可結合圖7進行說明。

圖7所示為結冰飛機動力學邊界限制與傳統的迎角限制的對比圖。其中黃色鉛垂面為傳統意義上的迎角限制面(本文GTM飛機的失速迎角約為0.312 rad，故本文取可用迎角為0.262 rad)，即飛行狀態位于平面左側時為安全飛行狀態(α<0.262 rad)；圖中紅色、綠色以及藍色曲面分別為無冰、結冰程度0.1以及結冰程度0.3時的動力學限制，此時安全裕度設定為0.2；且當飛行狀態在動力學限制范圍內時為安全飛行狀態，超出該動力學限制后開始發出危險告警，提示駕駛員此時飛行狀態已經具有潛在的飛行風險。由圖7可知，俯仰角和俯仰角速度較大時，迎角的可用范圍收縮且低于傳統迎角的限制值，且可用迎角范圍隨結冰程度的加劇而大幅度收縮，尤其是結冰較為嚴重的情況下，動力學限制的迎角可用范圍完全位于傳統可用迎角限制值的左側，此時使用傳統的飛行風險評估及預警方法將極易導致飛行事故，而使用動力學限制的安全預警方法在一定程度上可解決類似問題的發生，下面結合圖8～圖10說明此問題。

圖7 結冰情況下動力學限制與傳統迎角限制對比Fig.7 Comparison between dynamic limit and traditional angles of attack limit under icing condition

圖8所示為結冰飛機勻速下滑階段駕駛員進行拉桿操縱后的俯仰姿態響應曲線，且在飛行狀態在55.5 s左右超出可用迎角限制；圖9所示為相同情況的三維響應，且可知飛行狀態首先超出本文使用的動力學限制，隨后才超出傳統的可用迎角限制；圖10所示為此時飛行情況下基于動力學限制的安全預警方法和傳統可用迎角限制安全預警方法的對比圖，可知飛行狀態在40 s左右超出動力學限制的安全裕度值0.2(該裕度值的設定與傳統可用迎角裕度設置相同)。

圖8 時域仿真驗證Fig.8 Verification based on real time simulation

圖9 時域仿真結果三維顯示Fig.9 Results of real-time simulation in three-dimensional viewpoint

圖10 安全風險指標對比圖Fig.10 Comparison of safety risk indicator

由圖8～圖10可知，采用基于動力學限制的安全預警方法將在仿真開始后的40 s左右發出危險告警，而使用傳統迎角限制的安全預警方法將在55.5 s左右發出危險告警，因此在該飛行情況下，本文提出的安全預警方法相比于傳統的安全預警方法將提前15.5 s左右發現飛行風險。為進一步說明本文方法的普適性及工程應用意義，后文將采用飛行模擬的方式對其在不同飛行情況下的優越性進行驗證。

3 結冰飛機著陸下滑階段飛行訓練

基于MATLAB/Simulink和Flightgear仿真軟件搭建結冰飛機的飛行仿真訓練系統(如圖11所示)，并通過大量仿真訓練說明并驗證動力學限制安全預警方法的優越性并初步實現對結冰感知能力的訓練，為說明本文設計的飛行仿真系統在分析飛機著陸過程運動特性的適用性和飛行訓練效果，駕駛員由本文作者擔任并進行如下仿真訓練(本文作者為沒有任何飛機駕駛經驗的學者，但掌握飛機著陸過程涉及的相關理論知識)。

本文采用拉飄[31]的著陸方式進行著陸，仿真從飛機位于海拔400 m高度、以2.5°航跡俯仰角和85 m/s速度穩定下滑開始，到飛機后輪接觸地面結束；著陸過程中，駕駛員在飛機位于海拔高度15 m左右時，對飛機進行拉飄著陸的相關操縱。為說明訓練效果并驗證使用動力學限制后模擬器的可用性，進行大量的著陸訓練，最終得到如下形式的著陸曲線。利用此訓練系統得到飛機未結冰情況下的某一著陸軌線如圖12所示，且此時2種安全預警方法的指示曲線如圖13所示。

圖11 結冰飛機飛行仿真訓練系統結構圖Fig.11 Structure chart of flight simulation training system for icing aircraft

由圖12(x為航向距離)可知此時飛機的著陸軌線滿足著陸標準，飛機接觸地面時的下沉速度0.561 m/s同樣滿足設計要求，且通過對飛行狀態的綜合分析可知，利用本文設計的仿真系統進行著陸過程的相關研究是可靠、有效的。

圖13所示為飛機未結冰情況下，著陸過程中本文設計的安全預警方法和傳統安全預警方法顯示結果的對比。由圖可知，2種方法均未達到安全預警限制范圍，且變化趨勢基本一致，因此說明本文提出的安全預警方法不會影響飛機的正常飛行，這也為該方法的工程應用提供前提保證。

圖12 干凈飛機著陸軌線Fig.12 Landing trajectory of clean aircraft

圖13 飛機著陸過程安全預警方法對比Fig.13 Comparison of safety warning methods during landing phase

圖14所示為飛機結冰后(結冰因子取0.3)，駕駛員根據不同的安全預警方法進行操縱應對所產成的不同著陸軌線。紅色實線為飛機結冰后駕駛員不進行任何補償操縱而形成的軌跡；綠色實線為飛機結冰后駕駛員根據傳統的迎角限制進行操縱而形成的軌跡；藍色實線為飛機結冰后駕駛員根據動力學限制安全預警方法進行操縱應對所形成的軌跡。另外根據2種安全預警方法進行的操縱應對方法是相同的，即：在保證飛機不超過安全限制的同時保證飛機速度不低于失速速度、不撞地且滿足著陸要求，且在危險情況下通過操縱改飛機下滑過程為爬升或平飛過程，以便除冰后進行二次著陸。由于本文主要進行安全預警方法的優越性說明，故不進行二次著陸的相關仿真，因此結冰飛機改出過程的仿真停止于飛機脫離飛行風險或發生事故。另外圖中藍線和綠線為駕駛員操縱后，軌線形式出現率最高的情況，即：基于動力學限制進行改出操縱訓練時也存在失敗情況，但概率相對較低，同樣基于傳統的迎角限制進行改出訓練時也存在成功的情況，同樣概率也較低。

圖14 基于不同安全預警方法的結冰飛機著陸軌跡Fig.14 Landing trajectories of icing aircraft based on different safety warning methods

通過訓練和圖14和圖15可知，基于動力學限制的安全預警方法在此種情況下將提前于迎角限制方法約4.5 s進行危險告警(2.2節理論值為15.5 s，是因為駕駛員不進行任何補償操縱所致，因為駕駛員進行補償操縱必定在一定程度上提高飛行安全，也因此減少了2種方法的差異，另外該時間針對不同的飛行情況，此處給出具體時間僅為說明文中提出方法的優越性)。且通過大量改出操縱訓練可知，在37.5 s進行結冰飛機的著陸改出的成功率(約90%)將明顯高于在40～42 s進行結冰飛機的著陸改出操縱的成功率(約30%)，且難度也明顯降低，因此可在一定程度上說明基于動力學限制的安全預警方法進行飛行訓練可提高結冰飛機的著陸飛行安全。

圖15 結冰飛機安全預警值對比Fig.15 Comparison of safety warning for icing aircraft

4 結 論

本文基于微分流形理論確定了結冰飛機的動力學邊界，并詳細分析了結冰程度對邊界的影響規律，在此基礎上利用動力學邊界的特性對飛行風險進行了量化，最后基于MATLAB/Simulink和Flightgear仿真平臺搭建了飛行仿真訓練系統，并對上述方法的優越性及可用性進行了驗證。具體結論如下：

1) 將非線性穩定域作為飛機的動力學邊界，該動力學邊界可綜合考慮飛機各狀態之間的動力學耦合情況。研究了不同結冰情況下結冰飛機動力學邊界的變化情況，結果表明，嚴重結冰導致動力學邊界明顯收縮，且耦合特性表現突出。

2) 利用動力學邊界對飛行風險進行了量化，并在此基礎上提出一種可考慮外界因素影響的動力學限制安全預警方法，該方法考慮了動力學的耦合特性，可提前于傳統迎角限制預警方法進行飛行風險告警。

3) 將基于動力學邊界的動力學限制安全預警方法應用于飛行仿真訓練系統，并對結冰飛機的著陸飛行科目進行了仿真模擬訓練，最終驗證了該方法的可行性及相比于傳統方法的優越性。