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

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

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

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

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

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

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

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

1 飛機縱向質(zhì)點動力學模型

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

(1a)

(1b)

(1c)

(1d)

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

飛機結(jié)冰前后的氣動導數(shù)變化可表示為[13]

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

(2)

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

2 可達集計算

飛行器的動態(tài)特性可用常微分方程描述為

(3)

式中：x∈Rn表示n維狀態(tài)變量；u∈U表示系統(tǒng)的輸入變量。

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

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

圖1 目標集與可達集的關(guān)系Fig.1 Relationship between object set and reachability set

水平集方程可以表示為

(4)

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

目標集G0可由水平集函數(shù)表示為

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

(5)

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

(6)

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

(7)

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

(8)

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

3 飛機著陸階段分析

3.1 飛機著陸參數(shù)

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

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

(9)

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

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

表1 參數(shù)取值范圍Table 1 Range of parameters

3.2 飛機著陸可達集

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

(10)

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

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

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

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

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

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

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

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

圖3中，黃色網(wǎng)格部分為結(jié)冰程度為0.1時的可達集，紅色部分為不結(jié)冰飛機的可達集?？梢钥闯觯Y(jié)冰后飛機的安全包線縮小，可達集范圍收縮，駕駛員可以安全操縱的區(qū)域變小，安全裕度降低。圖4中，黃色網(wǎng)格部分為結(jié)冰程度為0.2時的可達集，紅色部分為結(jié)冰程度為0.1時飛機的可達集?？梢钥闯觯S著結(jié)冰程度的增大，可達集縮小，安全駕駛范圍縮小，飛機的操縱需要更加謹慎，發(fā)生危險的概率增加。

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

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

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

4 基于可達集的飛機結(jié)冰風險評估

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

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

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

Pd=1

VVmax，γ<γminorγ>γmax

(11)

式中：Pd為發(fā)生危險的概率。

4.1 參數(shù)提取介紹與極值理論

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

(12)

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

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

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

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

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

(13)

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

(14)

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

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

當ξ<0時，G(z)表示W(wǎng)eibull分布。

4.2 風險概率計算

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

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

γ+0.018 9V=0.727 4

(15)

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

(16)

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

d=-(γ+0.018 9V)

(17)

基于提取的參數(shù)，本文對根據(jù)速度V和航跡傾角γ的參數(shù)極值計算所得的變量d進行參數(shù)辨識，分別可以得到變量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)

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

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

直線為

γ+0.018 9V=0.854 6

(20)

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

(21)

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

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

(22)

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

γ+0.018 9V=0.985 2

(23)

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

(24)

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

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

(25)

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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