基于可達集方法的結冰飛機著陸階段安全風險評估

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

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

飛機著陸是整個飛行過程中飛行員參與控制的重要環節之一，對于考慮人為因素的安全性分析有重要意義。波音公司統計數據顯示，2003—2012年全球全重27 216 kg以上的商用噴氣飛機相關事故中，23%的嚴重事故和17%的事故死亡人數發生在著陸階段[1]。著陸階段飛機的飛行高度較低，如果飛機狀態出現不穩定，飛機調整狀態的時間非常短，另外由于飛機著陸過程中會遭遇風切變和微暴流，更易使飛機發生危險，因此飛機著陸階段的安全性問題尤為突出。

飛機的安全飛行是航空工業發展中不可忽視的重要問題，而造成飛機失事的原因有很多，其中，結冰就是主要原因之一。自從固定航線出現以來，結冰造成的災難性后果就不斷困擾著人們。飛機著陸時高度下降，會穿過云層等易發生結冰的環境。一旦機翼或尾翼結冰，將有可能造成飛機失速墜毀，此類事故曾經多次發生并且造成慘重的后果。2001年1月4日，兩架運-8飛機在著陸過程中因尾翼結冰相繼墜毀，造成12名機組成員和6名地面人員死亡，2人受傷。2009年美國科爾根航空公司3407號航班的DHC-8飛機在夜間儀表進近著陸過程中，由于飛機飛行速度低于結冰后的失速速度，駕駛員操縱不當引發重大飛行事故，造成50人遇難[2]。

目前各種大飛機都使用防/除冰裝置來應對結冰現象，然而由于防/除冰裝置體積和質量比較大，防/除冰效果有限，難以從根本上解決結冰的影響。另外，根據飛機的動力學特性和結冰后飛行特性的改變，從控制穩定上應對結冰后飛機失穩的方法正在逐漸發展，并有望成為解決結冰對飛行影響的有效手段。ATR-72等飛機為應對結冰對飛行性能的影響，通過改進飛行邊界控制保護系統，考慮最嚴重結冰情況并制定對應的邊界保護限制，當飛機遭遇結冰，結冰保護系統就開始工作[3-4]。美國伊利諾伊州立大學的Bragg等提出了“飛機智能防冰系統”的設想[5-6]。NASA Glenn研究中心啟動智能防冰計劃，開展控制與傳感、安全性、空氣動力學等學科組合的系統研究，并通過飛行試驗進行了驗證[7]。國內一些大學和科研機構也針對飛機結冰后穩定包線的改變進行研究，對結冰模型和結冰影響進行分析，提出了有關結冰參數辨識的研究方案。

可達集方法最初用于研究變化的流體外形，近年來利用可達集方法開展飛行安全分析得到一定的應用。文獻[8]利用水平集方法對飛機的自動著陸過程進行安全性分析，文獻[9]以飛機的筋斗動作為研究對象分析飛機機動過程中的可達集包線，文獻[10]將可達集方法用于四旋翼飛行器的后翻機動安全操縱范圍的確定中，分析不同模式下狀態參數的約束條件及其對應的可達集。

本文建立了飛機結冰前后的縱向動力學質點模型，運用可達集理論分析不同結冰程度下的四維可達集，得到飛機的飛行安全包線。從理論上說明結冰影響飛行安全的機理，對結冰飛機的著陸操縱提出指導，達到飛機在結冰狀態下安全著陸的目的。通過統計學方法確定極值理論所需的關鍵參數，建立飛機風險概率計算模型，根據可達集結果提取的航跡傾角和速度組成的安全包線，分析不同結冰條件下的風險概率。根據結冰風險概率對駕駛員操縱提出指導，創新地提出利用可達集評估結冰風險，所得的風險評估結果對于研究結冰引起的飛行安全和適航性問題具有重要意義。

1 飛機縱向質點動力學模型

本文選擇雙水獺飛機為研究對象，其模型數據參考文獻[11]。設飛機為剛體，為簡化研究，在不考慮縱向和橫航向耦合基礎下，提取飛機縱向動力學質點方程，以縱向的速度、航跡傾角、俯仰角速度和迎角4個參數展示飛機縱向運動的變化，其動力學方程[12]可寫為

(1a)

(1b)

(1c)

(1d)

由于結冰對飛機動力學的影響非常復雜并且結冰條件多種多樣，運用飛行試驗獲得飛機結冰條件下的氣動參數難以實現，而采用結冰影響模型的方法評估飛機結冰嚴重程度具有一定科學性并且簡單實用。因此本文采用結冰影響模型分析結冰前后飛機氣動導數的變化。

飛機結冰前后的氣動導數變化可表示為[13]

CA(iced)=(1+ηkiced)C(A)

(2)

式中：C(A)為飛機干凈翼型下的某個氣動導數值；CA(iced)為飛機發生結冰之后的氣動導數值；kiced為某個氣動參數受結冰影響的程度；η為結冰程度的大小，不結冰時其值為零，η隨結冰程度增加而增加。

2 可達集計算

飛行器的動態特性可用常微分方程描述為

(3)

式中：x∈Rn表示n維狀態變量；u∈U表示系統的輸入變量。

系統的初始狀態受控制量u∈U的作用，如果能夠在一定時間t∈[0,τ]內到達目標集，這些所有初始狀態的集合就是反向可達集[14-15]。

在分析系統方程時，由上述目標集的說明可知，可達集內的狀態都能夠在某一控制下經過給定時間后進入目標集，而可達集外的狀態則無論經過什么控制，都不能在給定時間內進入目標集。如圖1所示，B點和C點都是可達集內的狀態點，在某一控制量作用下，經過一定時間最終會進入目標集并停留，而可達集外的狀態點A無論在怎樣的控制作用下都無法進入目標集。

圖1 目標集與可達集的關系Fig.1 Relationship between object set and reachability set

水平集方程可以表示為

(4)

式中：φ(x,t)為水平集函數。

目標集G0可由水平集函數表示為

G0={x∈Rn|φ(x,0)≤0}

(5)

通過計算如下Hamilton-Jacobi方程的黏性解能夠得到目標集G0對應的可達集

(6)

式中：Hamilton函數H(x,p)為

(7)

u*(x,p)=arg maxpTf(x,t,u)

(8)

關于可達集理論的介紹，以及可達集的求解過程和步驟分析可參考文獻[9]。

3 飛機著陸階段分析

3.1 飛機著陸參數

飛機在著陸過程中，發動機推力T為飛機著陸的一個輸入控制量，假定發動機推力方向與飛機航跡方向一致；另一個輸入控制量是飛機的俯仰舵偏角。發動機推力和飛機俯仰舵偏角都有取值限制，由于結冰程度的影響，迎角最大值即失速迎角會改變。

在下滑階段，飛機沿著下滑斜率飛行，并且下滑角度必須保持在理想下滑傾角的變化范圍內，因此航跡傾角的取值范圍是[γmin,γmax]，并且γmin=γ0-dγ，γmax=γ0+dγ，γ0為理想的航跡角。在著陸過程中的速度也有限制，著陸速度超過最大限制會對飛機結構造成損壞，可能引起飛機沖出跑道的危險；而著陸速度超出最小限制會造成失速，失速速度可表示為

(9)

式中：ρ為空氣密度；CLmax為最大升力系數。速度范圍取為不同結冰程度下的安全速度限制，高度則根據實際著陸情況，確定在一定正值范圍內。

綜上所述，將著陸即將結束時刻的狀態范圍定義為目標集。如表1所示，不同的結冰程度對應飛機著陸不同的失速速度和失速迎角，根據飛機著陸可承受的沖擊得到著陸允許的最大速度。

表1 參數取值范圍Table 1 Range of parameters

3.2 飛機著陸可達集

由式(7)得到動力學模型下的Hamilton函數為

(10)

式中：p1、p2、p3、p4分別為φ(x,t)對狀態V、γ、q和α的偏導。最優控制量的取值是通過將輸入的不同取值組合代入到系統中計算可達集，并對可達集大小進行比較得到的。最優控制量的選取可以通過計算機的仿真試驗進行，同時進行多組輸入得到最大的可達集，最大的可達集對應的輸入被視為最安全的操縱指令。這一指令指導駕駛員操縱飛機，使飛機達到最安全的飛行狀態。

目標集是根據速度、航跡傾角、俯仰角速度和迎角的合理范圍確定的，而可達集由目標集擴展而成，飛機著陸狀態在這一范圍內是穩定和安全的。可達集的狀態經過一定的時間會進入目標集，一旦狀態偏離可達集范圍，飛機狀態在任何輸入條件下都無法進入目標集，著陸將不能達到安全狀態。

為將可達集結果可視化顯示，在著陸安全包線的計算分析過程中，選取狀態的3項參數作為性能指標，構建一個三維可視圖形，根據圖形的區域范圍確定安全包線的范圍。本文將狀態參數迎角切片進行數據分析，選取速度、航跡傾角和俯仰角速度作為求解目標函數。

首先計算結冰程度為0時給定目標集下的反向可達集，其對應的失速迎角為17.5°，結果如圖2 所示，網格內部的藍色部分為給定的目標集，黃色網格所包圍的區域是對應的可達集。以著陸下滑結束段的安全狀態為起始狀態，對飛機著陸階段的飛行狀態模型進行反時間方向的求解，所得的可達集為飛機的飛行安全包線。

在計算所得的可達集范圍內，任何一個狀態點都能夠通過合理的操縱，進入理想的目標集內。當飛機狀態偏離可達集范圍時，飛機在任意操縱下均無法到達理想的著陸飛行狀態。因此，駕駛員操縱或機載計算機的控制律應該確保飛機狀態始終位于可達集內。從圖2中看出，未結冰時，飛機的安全飛行狀態空間較大，飛機的安全裕度較大，駕駛員的安全操縱范圍也較大。

可達集刻畫的結果是存在誤差的，誤差主要由計算網格數目決定。網格數目越多，計算精度越高，誤差相對較小，計算機的計算時間增長。為使計算時間縮短，同時計算精度誤差相對于飛行安全范圍能夠忽略，應選取合適的網格數目進行計算。

圖2 結冰程度為0時目標集與可達集Fig.2 Object set and reachability set when icing degree at 0

結冰程度不同時，對應的失速迎角等變量限制也不同。結冰程度是0.1時的失速迎角為16.9°，結冰程度是0.2時的失速迎角為15.25°，因此需要在不同控制和不同的目標集下進行計算。另外，結冰會導致飛機的氣動導數發生變化，因此不同結冰程度的可達集有明顯的差異，其仿真結果如圖3和圖4所示。

圖3中，黃色網格部分為結冰程度為0.1時的可達集，紅色部分為不結冰飛機的可達集。可以看出，結冰后飛機的安全包線縮小，可達集范圍收縮，駕駛員可以安全操縱的區域變小，安全裕度降低。圖4中，黃色網格部分為結冰程度為0.2時的可達集，紅色部分為結冰程度為0.1時飛機的可達集。可以看出，隨著結冰程度的增大，可達集縮小，安全駕駛范圍縮小，飛機的操縱需要更加謹慎，發生危險的概率增加。

圖3 不結冰與結冰程度為0.1時的可達集對比Fig.3 Comparison of reachability sets between no ice and icing degree at 0.1

圖4 結冰程度為0.1和0.2時的可達集對比Fig.4 Comparison of reachability sets between icing degree at 0.1 and 0.2

飛行在安全包線內可以保證安全飛行，可達集計算的范圍可以表征飛行包線的范圍。駕駛員的操縱和外界瞬時擾動使飛機狀態的變動在這一范圍內時，飛機可以通過操縱使其安全著陸。

4 基于可達集的飛機結冰風險評估

第3節中的可達集是基于四維動力學方程計算的，根據動力學方程，可以將速度與航跡傾角隨時間的變化率提取出來進行計算[16]，不考慮俯仰角速度和迎角對飛行風險影響，分析飛機結冰飛行風險，將方程取為二維的動力學方程計算可達集。

本節將二維可達集作為安全評估的標準，分析不同結冰程度時的風險概率[17-18]。由于結冰飛行試驗風險很大，因此一般的風險評估使用地面模擬飛行的數據，在此基礎上分析風險概率。

由第3節著陸安全可達集分析可知，著陸時，飛機的安全穩定狀態應該處于速度、航跡傾角的合理范圍內，因此將飛行狀態超出可達集范圍作為危險發生的判據。速度大于或小于可達集限制認為飛機可能發生危險，航跡傾角超出給定范圍認為飛機將發生危險。用概率表達式表示為[19]

Pd=1

VVmax，γ<γminorγ>γmax

(11)

式中：Pd為發生危險的概率。

4.1 參數提取介紹與極值理論

為計算飛機結冰條件下的飛行風險概率，首先設定初始飛行狀態：初始飛行的高度為100 m，速度為80 m/s，航跡傾角為-3°。將駕駛員操縱模型作為系統的輸入，進行蒙特卡羅仿真，得到多組仿真結果。駕駛員操縱模型根據文獻[20]所提供的模型確定，其模型可以表示為

(12)

式中：x0為操縱量初始值；σ為模型參數；u(x-x0)為單位階躍函數。

根據統計學和概率論的知識可知，仿真次數越多，即得到的樣本容量越大，統計結果越接近真實情況，因此仿真次數需要足夠多。同時，為減少計算機的計算任務量，本文取1 000次仿真結果進行分析。

飛行風險可以通過關鍵飛行參數超出其邊界值的概率來進行評估[21-22]，在著陸階段，飛行狀態是關鍵飛行參數，當飛行狀態超出安全包線時，可認為飛行風險事件發生。飛行風險是小概率事件(事件發生概率<10-9)，很難通過計算所有發生事件次數除以總的仿真次數來得到風險事件發生的概率。而此類事件一旦發生就會帶來嚴重危害，例如金融風險、巨額保險的賠付、重大的自然災害、重大人為事故等。

極值理論在計算這種低頻高危風險事件的概率方面具有其特有的優勢[23-24]。極值理論認為，在不需要知道獨立同分布隨機變量的累積概率分布條件下，就可以得到極值的分布函數，這是因為隨著樣本容量的增加，極值的分布漸近地趨于一個確定的分布函數。

如果存在常數序列{an>0}和bn，使得當n→∞時：

(13)

式中：ξ為非退化分布函數，那么G必屬于廣義極值(Generalized Extreme Value，GEV)分布族，即

(14)

式中：-∞<μ<∞，σ>0，-∞<ξ<∞。

當ξ>0時，G(z)表示Fréchet分布。

當ξ<0時，G(z)表示Weibull分布。

4.2 風險概率計算

將飛機質點動力學方程降維成二維方程，僅對速度和航跡傾角進行可達集的求解，可以得到飛機的速度和航跡傾角的安全范圍。將這一范圍作為飛機風險評估的判據，狀態超出安全范圍認為飛機將發生危險。在計算風險概率時，風險的定義是飛機失控，即接地之前或者接地瞬間飛機失去控制，不包括飛機可控狀態下飛機著陸后沖出跑道或者飛機結構受損。因此不將飛機速度較大并且不適合著陸的情況定義為風險，僅考慮航跡傾角或速度超過最小限制的情形。

不結冰時得到的可達集結果如圖5所示。根據結果，航跡傾角或速度的較小極值在紅色虛線以內，紅色虛線所在直線為

γ+0.018 9V=0.727 4

(15)

安全包線所在區域應該滿足以下條件:

(16)

為計算風險概率，定義變量d，其值為

d=-(γ+0.018 9V)

(17)

基于提取的參數，本文對根據速度V和航跡傾角γ的參數極值計算所得的變量d進行參數辨識，分別可以得到變量d的極值分布GEV模型。其中對于變量d來說，μd=-1.608 3，σd=0.042 5。概率計算方程可表示為

(18)

由式(18)計算可得，干凈外形條件下，d=-0.727 4，飛行風險概率是

P=1-G(d)=9.962 1×10-10

(19)

結冰程度為0.1時，所得到的可達集結果如圖6所示。根據結果，仍然將航跡傾角或速度的較小極值限制在紅色虛線以內，紅色虛線所在

圖5 結冰程度為0時的安全包線Fig.5 Safety envelope when icing degree is 0

直線為

γ+0.018 9V=0.854 6

(20)

安全包線所在區域應該滿足以下條件：

(21)

結冰程度為0.1時，氣動力和氣動特性發生改變，安全包線收縮，飛行狀態超出可控范圍的可能性大大提高，風險概率增大，此時d=-0.854 6，根據Gumbel模型計算此時的風險概率為

P=1-G(d)=1.986 9×10-8

(22)

結冰程度為0.2時,所得可達集結果如圖7所示。根據上述計算過程，航跡傾角或速度的較小極值限制在紅色虛線以內，紅色虛線所在直線為

γ+0.018 9V=0.985 2

(23)

安全包線所在區域應該滿足以下條件：

(24)

當結冰程度為0.2時，由于氣動特性改變進一步加劇，安全包線進一步收縮，此時d=-0.985 2，風險概率由概率模型計算：

P=1-G(d)=4.292 7×10-7

(25)

當風險概率小于10-8時，認為結冰嚴重程度屬于“極小的”范疇，可認為飛機不會發生風險。當風險概率范圍為10-8≤P<10-6時，駕駛人員應當立即開啟除冰裝置。當風險概率為P≥10-6時，駕駛員應當立即改變飛行路線，使飛機駛離結冰區。當速度或航跡傾角大于限制值時，駕駛員不應繼續著陸，應當改出著陸階段，調整狀態后重新進行著陸。

圖6 結冰程度為0.1時的安全包線Fig.6 Safety envelope when icing degree is 0.1

平尾涉及俯仰舵面的操縱，一旦結冰，飛機縱向的俯仰操縱將受到影響。根據文獻[25]的研究，平尾結冰相對于機翼結冰來說，對飛機的縱向升力和阻力產生的影響較小。因此在相同程度的結冰條件下，平尾結冰對應的反向可達集結果范圍比機翼結冰對應的可達集范圍更大。仿真結果如圖8所示，結冰程度同為0.1時，相比于機翼結冰，平尾結冰對應的可達集結果在低速區域和負航跡傾角區域范圍更大，而在高速區域范圍有略微的收縮，這與實際情況相對應[26]。

圖7 結冰程度為0.2時的安全包線Fig.7 Safety envelope when icing degree is 0.2

圖8 平尾結冰與機翼結冰對應的安全包線Fig.8 Safety envelopes when tailplane icing and wing icing

5 結 論

基于可達集理論，以結冰條件下的飛機質點動力學模型為分析對象，利用可達集方法對其著陸階段的安全包線進行求解，根據安全包線和極值理論計算不同結冰程度下的著陸風險概率，并對駕駛員操縱提出指導，主要結論如下：

1) 可達集適合用于計算飛機飛行的安全包線，計算結果直觀可視。結冰飛機著陸時使用可達集計算安全包線能夠指導駕駛員安全著陸，飛機狀態一旦偏離可達集，駕駛員應該立即改變飛機的操縱，改出著陸階段，避免發生惡劣后果。

2) 隨著結冰程度的加劇，飛行的可達集范圍縮小，穩定裕度降低，更容易發生危險。其主要原因是飛機結冰后的各項氣動導數發生了較大變化，動力學特性變差，操縱難度增加，導致飛機的穩定性降低，風險系數增加。

3) 以可達集為安全包線的結冰飛行風險評估方法，能夠得到飛機超出結冰后安全包線的風險概率，該方法能夠用來定量地計算不同結冰程度下的飛行風險，能夠指導駕駛員進行合理的操縱，從而規避危險。

水平集方法可求解出飛機著陸的安全包線，風險評估方法結合可達集計算結果能評估出飛機結冰后的飛行風險。然而著陸過程中飛行風險涉及的因素有很多，本文利用可達集計算時采用的模型為簡化的模型，所計算的結果仍然需要進一步的研究才能更加準確。后續將增加飛機系統維數，建立更加準確的仿真框架和模型，并對仿真模型進行有效驗證。