伍強，徐浩軍，裴彬彬,*，魏揚，張楊

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

2. 中國人民解放軍 95633部隊，成都 611530

在一定氣象條件下，大氣中的過冷水滴與飛機迎風撞擊就有可能導致飛機表面結冰，波音公司統(tǒng)計結果表明，飛機結冰是誘發(fā)飛行失控(Loss of Control, LOC)嚴重事故的3大因素之一，飛機結冰問題長期以來一直受到人們的高度關注。2004年11月21日，東方航空公司一架CRJ-200飛機由于停放后除冰不徹底導致起飛過程中墜毀，機上55人全部遇難。2006年6月3日，某型飛機在多次穿越云層后結冰導致飛機失控，機上40人全部遇難。因此，有必要深入研究結冰條件下的飛機飛行風險，探索結冰對飛機飛行安全的威脅程度，從而制定合理的結冰風險規(guī)避策略。

研究結冰飛行風險問題必然涉及到飛行風險評估方法，傳統(tǒng)的飛行風險評估方法以定性分析為主，定量分析也僅依據(jù)靜態(tài)安全性指標。Mohaghegh等將系統(tǒng)動力學理論、貝葉斯網(wǎng)絡、事件圖和故障樹等理論應用于飛行風險管控中，提出了航空維修系統(tǒng)風險分析方法。Brooker將馬爾科夫鏈、故障樹、事件樹等方法成功應用于飛行風險評估中。傳統(tǒng)的飛行風險分析方法主要針對元部件故障或人為失誤等靜態(tài)因素，無法研究飛機飛行過程中的動態(tài)飛行風險，因此國內外已經(jīng)將研究重點放在運用仿真手段開展復雜外部條件下飛行風險量化評估，de Mendon?a等提出了基于飛行模擬仿真的飛機安全性分析方法，并證明了運用飛行模擬開展飛行安全保障工作的必要性。Blum等討論了基于蒙特卡羅模擬開展飛機安全性分析的可行性。蒙特卡羅法是在已知隨機變量分布規(guī)律的條件下，依據(jù)變量概率分布隨機提取變量樣本數(shù)據(jù)。該方法能夠有效模擬隨機變量的動態(tài)變化過程，準確反映隨機變量的參數(shù)分布規(guī)律，解決復雜系統(tǒng)中的隨機性問題，適用于提取飛行仿真中隨機變量參數(shù)樣本。飛行安全關鍵參數(shù)極值樣本均具有明顯的厚尾分布特性，這種分布形式在低頻高危風險事件(如金融風險、自然災害等)中較為常見。武朋瑋等基于可達集方法和極值理論對結冰條件下的飛機著陸風險進行了量化評估，其研究對于指導飛行員操縱具有一定的應用參考價值，但是其采用的是一元極值分布模型，沒有考慮各個參數(shù)之間相關性。多元極值理論與Copula函數(shù)模型快速發(fā)展，已經(jīng)成為研究高危低頻風險事件的有力研究工具，李哲、魏揚構建了適用于飛行風險量化評估的多元極值Copula模型，其選擇的飛行參數(shù)主要包括速度、迎角和滾轉角，對于飛機的姿態(tài)角及相應的角速度之間的相關性沒有進行深入的研究。

運用蒙特卡羅方法和極值理論進行復雜條件下飛行風險量化評估的前提和關鍵是要設定準確的事故發(fā)生判定條件。目前的飛行風險判定標準主要針對某一飛行參數(shù)是否超過其極限值，由于各飛行參數(shù)極限值是依據(jù)經(jīng)驗和試飛得到一維數(shù)據(jù)，沒有考慮到各個飛行參數(shù)之間具有的強烈相關性，僅僅根據(jù)一維極限值作為判斷標準開展飛行風險評估，難以全面有效的描述飛行安全特性。而吸引域方法可用于計算非線性系統(tǒng)的吸引域及穩(wěn)定邊界，進而作為飛行風險判定條件進行飛行風險的定量計算。

本文首先建立了典型運輸類飛機縱向動力學模型和結冰影響模型，基于吸引域方法得到了不同結冰程度條件下的飛機基于迎角和俯仰角速率參數(shù)的二維吸引域及穩(wěn)定邊界。建立了典型的人-機-環(huán)系統(tǒng)模型，運用蒙特卡羅仿真方法，以典型的駕駛員操縱——升降舵脈沖作為輸入得到了迎角和俯仰角速率等關鍵參數(shù)極值。根據(jù)極值理論，建立了適用于飛行風險量化評估的二元極值Copula模型，通過遺傳算法辨識模型參數(shù)，根據(jù)多種擬合優(yōu)度檢驗方法確定了最優(yōu)分布模型，以飛行狀態(tài)超出穩(wěn)定邊界作為發(fā)生飛行事故的判定條件，計算不同結冰嚴重程度下的風險概率，進而根據(jù)結冰風險概率對駕駛員操縱提出指導。本文創(chuàng)新地提出利用吸引域和二元極值理論評估結冰風險，所得風險評估結果對于研究結冰引起的飛行安全和適航性問題具有重要意義。

1 模型建立

1.1 飛機縱向動力學模型

飛機縱向非線性動力學模型可表示為

(1)

式中：、、、分別表示空速、迎角、機體俯仰角、俯仰角速率；、分別為推力在機體坐標系、軸上的投影；是由推力產(chǎn)生的俯仰力矩；表示飛機的轉動慣量；、、分別表示飛機受到的升力、阻力與俯仰力矩，其計算方式為

(2)

1.2 結冰影響模型

由于開展結冰飛行試驗難度大、成本高，本文通過建立結冰影響模型來分析結冰后氣動參數(shù)的變化。Hui等利用最大似然法獲取飛機結冰前后的氣動參數(shù)，初步建立了結冰外界環(huán)境無量綱參數(shù)與退化系數(shù)相關的結冰影響模型， Lampton與Valasek提出了一種基于試飛數(shù)據(jù)定義每個氣動參數(shù)對應的退化因子的簡化方法，進而對不同結冰嚴重程度及結冰分布的動力學響應進行分析；此外，他們還研究了非對稱結冰影響模型，用于分析除冰系統(tǒng)單側故障時的動力學特性，當前常用的是Bragg等提出的結冰參量影響模型，模型可表示為

=(1+)

(3)

式中：與分別表示任意結冰前后飛機的性能、穩(wěn)定性與控制參數(shù)或其導數(shù)；表示結冰對飛機氣動參數(shù)的影響參數(shù)，它與飛機自身尺寸、飛行狀態(tài)或與飛機受結冰影響的敏感性相關；為飛機結冰程度參數(shù)，只與氣象條件有關。

1.3 飛機本體模型

飛機本體非線性動力學模型可表示為

(4)

式中：為狀態(tài)向量，包含飛行速度、迎角、側滑角、四元數(shù)、俯仰角速率、滾轉角速率、偏航角速率和空間位置參數(shù)，即

=[,,,,,,,,,,,,]

(5)

為控制向量，包括油門偏度指令、升降舵偏度指令、副翼偏度指令和方向舵偏度指令，即

=[,,,]

(6)

運用四元數(shù)法構建飛機動力學模型具體過程見附錄A。

1.4 駕駛員操縱模型

進行蒙特卡羅仿真時，需要將駕駛員操縱作為輸入，考慮到駕駛員操縱的隨機性，認為駕駛員的隨機性操縱參數(shù)服從參數(shù)為的瑞利分布：

(7)

式中：(-)為單位階躍函數(shù);為操縱量初始值。

2 非線性動力學系統(tǒng)吸引域分析

2.1 基于平方和理論的多項式非線性系統(tǒng)吸引域及穩(wěn)定邊界估計

考慮多項式非線性自治動力學系統(tǒng)：

(8)

式中：()為狀態(tài)向量。

假定=是系統(tǒng)的局部漸近穩(wěn)定平衡點，那么平衡點附近的吸引域(Region of Attraction)定義為

(9)

式中：(,)為初始狀態(tài)為時動力學系統(tǒng)的解。所有初始狀態(tài)位于吸引域內的點最終都將收斂于平衡點。

若存在一個連續(xù)可微的函數(shù)：→使得

1)()=0且?≠,()>0。

2):={∈|()≤}有界。

3)?{∈|▽()()<0}∪{}。

則對于所有的∈，方程(1)的解存在且滿足lim()∈。從而，是方程(1)吸引域的子集。

滿足定理1的函數(shù)稱作李亞普諾夫函數(shù)(Lyapunov Function)，而可作為吸引域的一個估計。為最大化，可定義一個可變區(qū)域()，在約束∈的條件下最大化。

:={∈|()≤}

(10)

式中：是大于零的值；()為正定多項式。利用代數(shù)幾何中的Positivstellensatz定理，可將該最優(yōu)化問題轉化成下面的平方和規(guī)劃(Sum-of-Squares Programming, SOSP)問題：

(11)

式中：為所有個變量的多項式集合；表示個變量的平方和多項式集合。結合-迭代算法，可對上述最優(yōu)化問題進行優(yōu)化求解，即可估計出多項式非線性動力學系統(tǒng)的吸引域，吸引域的邊界也就是該動力學系統(tǒng)的穩(wěn)定邊界：

:={∈|()=}

(12)

在進行估算時，可假定吸引域為一橢球體，對于二維問題，可假設：

()=N=()+()

(13)

式中：=diag(,)為形狀函數(shù)。

利用SOS理論對非線性動力學系統(tǒng)吸引域求解的前提是建立動力學系統(tǒng)的多項式模型。為此，首先要在合理的區(qū)間范圍內將飛行動力學模型轉化為多項式模型，然后在此基礎上計算飛機在平衡點附近的穩(wěn)定域。

非線性動力學系統(tǒng)吸引域的計算過程如下：

利用多項式擬合對每一項氣動系數(shù)、推力特性、狀態(tài)變量的逆等建立起相應的多項式模型。

將建立起的各多項式模型代入常規(guī)動力學模型中，進而推導出飛機動力學模型的多項式表達式。

為便于計算，通過坐標變換將飛機多項式動力學模型的平衡點移動到原點，并去除不必要的項(系數(shù)很小或階數(shù)很高)。

基于SOS理論計算出飛機在原點附近的吸引域。

通過坐標反變換將計算出來的吸引域映射到以原平衡點為中心點的區(qū)域，該區(qū)域即為飛機狀態(tài)參數(shù)的局部漸近穩(wěn)定區(qū)域,即吸引域。

2.2 飛機縱向多項式動力學模型

進行非線性動力學系統(tǒng)吸引域的計算時，首先需要將常規(guī)的飛機非線性動力學模型轉化為多項式模型。常規(guī)的飛機非線性動力學模型表達式中存在著一些非多項式項，主要包括三角函數(shù)項、氣動力(矩)系數(shù)、速度的逆等。在轉換為多項式的過程中，多項式的階數(shù)越高，計算精度也就越高，但同時會帶來計算量的增加。因此，需要確立合適的多項式階數(shù)。

對于三角函數(shù)的轉化可采用泰勒級數(shù)展開方法，一般展開至三階即可。在[(45°,45°]的范圍內，sin函數(shù)與cos函數(shù)的誤差都比較小，例如30°時，sin函數(shù)與cos函數(shù)估算值與實際值之間的誤差分別為0.065%與0.35%

常規(guī)的飛機非線性動力學模型(如式(1)所示)中含速度的逆，即1/。由于理想情況下，動力學系統(tǒng)的多項式模型不應該含有指數(shù)為負的項,為此，在合適的速度區(qū)間上用二次多項式來擬合速度的逆，從而得到近似的多項式表達式：

=++

(14)

式中：、、為待定常系數(shù)，通過曲線擬合來計算。飛機結冰一般在速度較低的情形下發(fā)生，考慮到飛機的最小平飛速度限制，在對速度進行擬合時，擬合區(qū)間選定為80～150 m/s，在此區(qū)間內，擬合誤差均在10數(shù)量級，例如當速度為100 m/s 時，誤差為0.791%，擬合結果為

=7153 1×10-

2447 4×10+0027 4

(15)

在飛機動力學模型中，空氣動力(矩)以及發(fā)動機的推力與飛機狀態(tài)參數(shù)、操縱參數(shù)有關，一般用插值算法對氣動參數(shù)表進行計算。本文以通用運輸類模型(Generic Transport Model, GTM)風洞試驗得出的數(shù)據(jù)為基準，經(jīng)過相似準則變換后形成TCM(Transport Class Model)模型進行分析計算。為便于開展吸引域分析，本文基于最小二乘擬合算法來計算氣動參數(shù)(包括升力系數(shù)、阻力系數(shù)、俯仰力矩系數(shù))的多項式表達式。這里以背景飛機在干凈構型下的升力系數(shù)為例進行說明，升力系數(shù)可分解為

(16)

式中：()是只與迎角有關的基本氣動參數(shù)，選用一元四階多項式進行擬合，得到其多項式表達形式為

()=-0791 6+4388 6-7790 9+

5137 2+0068 9

(17)

(,)是與迎角和升降舵操縱量有關的項，這里選定用二元二階多項式進行擬合，得到其多項式表達形式為

0034 6+0454 7+0002 5

(18)

(19)

升力系數(shù)各項的多項式擬合結果與真實值之間的對比如圖1所示。

圖1 升力系數(shù)各項擬合值與真實值之間的對比Fig.1 Comparison between fitted and true values of lift coefficient

10-2316 4×10-

0000 6+5587 4×10+

0000 9+0080 5-0022 5+

0022 9-0097 7-4415 1-

0016 7+2464 7+0003 7

(20)

將式(17)、式(18)、式(20)代入式(16)便可得到升力系數(shù)的多項式表達式。利用同樣的方法，可得到阻力系數(shù)、俯仰力矩系數(shù)的多項式表達式。在進行擬合計算時，非線性動力學模型中各角度均以(°)為單位。同時，需要指出的是，該多項式表達式成立的參數(shù)范圍需滿足飛機氣動參數(shù)插值表中給出的插值范圍限制。對于該型飛機而言，各參數(shù)的有效范圍為：∈[-10°, 60°],∈[-30, 30](°)/s,∈[-30°, 20°]，∈[80,150] m/s。

在進行飛機穩(wěn)定性分析時，通常不考慮發(fā)動機推力的變化，而是把發(fā)動機的推力維持在配平狀態(tài)，這里發(fā)動機的推力采用簡化的推力模型進行計算，認為發(fā)動機的推力只與油門偏度有關，油門偏度的取值范圍為0～100%。根據(jù)發(fā)動機推力值與油門偏度之間的差值數(shù)據(jù)，同樣可用多項式擬合的方法，得到發(fā)動機推力的多項式模型：

(21)

通過上述對常規(guī)動力學模型中非線性項的多項式轉化，可將常規(guī)的飛機縱向非線性動力學模型轉化為多項式形式的微分方程組。

2.3 基于SOS理論的吸引域計算

2.3.1 多項式動力學模型簡化

飛機縱向擾動有兩種典型的模態(tài)：長周期模態(tài)與短周期模態(tài)。飛機受到擾動后，首先是短周期模態(tài)參數(shù)、的變化較為迅速，一旦發(fā)散對飛行安全造成很大影響；長周期模態(tài)參數(shù)、的變化較為緩慢，在短時間內可視為不變，并且駕駛員具有較長的時間對其進行修正。因此，短周期多項式模型與全量多項式模型在短時間內的動態(tài)響應差別不大，其主要的區(qū)別在于擾動后飛行軌跡趨于平衡點階段，短周期模型能很快地趨于平衡點，而全量模型由于其中長周期模態(tài)的存在，狀態(tài)參數(shù)趨于平衡點過程中存在著低頻小幅震蕩，這些差別對于吸引域的計算影響不大。同時考慮到狀態(tài)變量增多引起的計算量增大的問題，本文在計算過程中將、視為定值(即配平值)，將全量多項式模型簡化為只含有短周期模態(tài)參數(shù)、的多項式短周期模型。

基于上述思想，飛機的初始飛行狀態(tài)設定為在=2 000 m高度上以=120 m/s的速度做水平直線飛行，多項式動力學模型的平衡點為：=120 m/s，=626°，=0 (°)/s，此時的輸入變量為：=2225，=049°。將全量多項式模型中的、、、分別替換為它們各自的配平值，便可得到代表飛機短周期模態(tài)的多項式模型。按照吸引域的定義，在進行多項式系統(tǒng)的吸引域估計時，平衡點為零點。為滿足這一條件，將模型的平衡點移動至原點，可以通過將、分別替換為+、+來進行計算。同時，為簡化計算，可略去系數(shù)小于10的項。

最終將常規(guī)的飛機非線性動力學模型(式(1))按照上述方法轉化得到的多項式短周期動力學模型為

0004 3+0056 5-0537 8+0891 2

(22)

1423 2-0783 5

(23)

需要再次指出的是，上述轉化的多項式動力學模型的有效范圍同樣必須滿足2.2節(jié)中多項式模型推導過程中的約束條件：對于迎角而言，其有效范圍為∈[-10°,45°]，俯仰角速率有效范圍為∈[-30,30] (°)/s，升降舵偏角有效范圍為∈[-30°,20°]，飛行速度有效范圍為∈[80,150] m/s。

2.3.2 吸引域及穩(wěn)定邊界計算

將式(13)中的形狀函數(shù)定義為：=diag(20°,30(°)/s)，便可得到正定多項式：

()=8207+3647 6

(24)

基于SOS理論，利用-迭代算法，便可得到飛機在原點附近的吸引域：

′={,∈|0603 6-0279 7+

0538 9≤0580 5}

(25)

相應的穩(wěn)定邊界為

′={,∈|0603 6-0279 7+

0538 9=0580 5}

(26)

上述計算吸引域的過程中，需要將平衡點從配平狀態(tài)點移至坐標原點，所以在獲得原點附近的吸引域′后，將吸引域中心平移至原配平狀態(tài)點，即可得到飛機在配平狀態(tài)點附近的吸引域：

={,∈|0603 6(-0109 3)-

0279 7(-0109 3)+0538 9≤0580 5}

(27)

吸引域的物理意義為：當飛機受到陣風或其他不利擾動后，只要其狀態(tài)參數(shù)位于吸引域內，飛機就能回到原來的平衡狀態(tài)。

多項式模型轉化過程中擬合函數(shù)的誤差會導致多項式模型的動態(tài)響應與常規(guī)的非線性動力學模型的響應存在誤差，進而會導致基于SOS理論計算得到的吸引域與飛機實際的穩(wěn)定域存在一定的出入。為驗證多項式模型轉化及SOS理論計算吸引域方法的有效性，將背景飛機常規(guī)動力學模型中的短周期模態(tài)參數(shù)在各個初始狀態(tài)點上的動態(tài)響應投影到-相平面上,與所計算出的吸引域進行對比分析，如圖2所示。每一條相軌跡始于平面內設定的初始狀態(tài)點矩陣，紅色曲線表示初始狀態(tài)點最終演化至發(fā)散的軌跡，藍色曲線表示初始狀態(tài)點最終演化至平衡點的軌,綠色曲線圍成的區(qū)域為計算出的吸引域。顯然，吸引域位于平衡點附近的穩(wěn)定范圍內且占據(jù)了穩(wěn)定范圍大部分的面積，始于吸引域內每一個狀態(tài)點均能回到平衡點，說明采用SOS理論計算得出的穩(wěn)定邊界雖然偏于保守，不能完全包含整個穩(wěn)定范圍，但基本上能夠反映出飛機在平衡點附近的穩(wěn)定范圍。

圖2 干凈外形時α-q的相平面圖Fig.2 Phase plane plot of α-q in clean shape

3 基于極值理論的飛行安全關鍵參數(shù)分布模型

3.1 一維極值分布模型

一維極值分布模型能夠有效地描述結冰情形下飛行安全關鍵參數(shù)極值樣本的概率分布情況和樣本數(shù)據(jù)邊際概率分布尾部特征。當前，其他高危低頻事件中應用最為廣泛的一維極值分布模型是廣義極值分布模型(GEV)。GEV模型由Fisher-Tippett定理推導所得，推導過程見附錄B。

3.2 一維飛行安全關鍵參數(shù)分布模型

針對飛行安全關鍵參數(shù)具有的厚尾特性，選取能夠有效描述尾部分布規(guī)律的5種分布模型，分別是極值理論中的極值分布(EV)、廣義極值分布(GEV)以及能夠描述尾部分布特性的對數(shù)正態(tài)分布(Logn)、指數(shù)分布(Exp)和威布爾分布(Wei)，同時選擇正態(tài)分布(Norm)用做對比分析。上述模型的分布函數(shù)分別為

(28)

(29)

(30)

(31)

(32)

(33)

在運用蒙特卡羅仿真方法獲取迎角及俯仰角速率極值樣本(見附錄C)的基礎上，通過遺傳算法辨識上述模型中的未知參數(shù)變量，辨識結果如表1所示。

表1 一維極值模型參數(shù)辨識結果Table 1 Parameter identification results of one-dimensional extreme value model

首先通過QQ圖檢驗法初步判定6種分布的擬合精度。在已知分布函數(shù)的未知參數(shù)取值后，可依據(jù)原樣本極值參數(shù)繪制QQ圖。從理論角度分析，極值參數(shù)的分布函數(shù)為(;,,…,)時，QQ圖中的圖像應近似為一條直線，若圖形偏離直線較大，則認為該種形式的分布不滿足樣本觀測值分布特性。上述目標樣本分布模型的QQ圖如圖3所示。

由圖3(a)可知，對于迎角極值，GEV分布模型的QQ圖最接近于直線，而其他分布模型的QQ圖均不同程度的有偏離趨勢，其中，Exp和Logn分布模型的偏離程度較低，而Norm、Wei和EV分布模型偏離較為嚴重。同理，由圖3(b)可知，對于俯仰角速率極值，GEV分布模型的線性程度最高。

圖3 QQ檢驗圖Fig.3 QQ inspection chart

表2 迎角極值擬合優(yōu)度檢驗Table 2 Goodness-of-fit test for extreme angle of attack

表3 俯仰角速率角極值擬合優(yōu)度檢驗Table 3 Goodness-of-fit test for extreme pitch rate

3.3 二元飛行安全關鍵參數(shù)Copula模型

通過Sklar定理可以看出，Copula函數(shù)能夠有效地分析多維聯(lián)合分布(詳見附錄D)。常見的二元Copula模型主要有Gumbel Copula模型(Gum)、Frank Copula模型(Frank)、Clayton Copula模型(Clay)及Joe Copula模型(Joe)。上述模型的分布函數(shù)分別為

(34)

(35)

(36)

Joe:(,)=1-[(1-)+(1-)-

(37)

式中：、分別為迎角極值和俯仰角速率極值的邊緣分布。

=()

(38)

=()

(39)

相應的概率密度函數(shù)(,)可表示為

(40)

通過遺傳算法辨識上述模型中的未知參數(shù)，辨識結果如表4所示。

表4 Copula模型參數(shù)辨識結果Table 4 Parameter identification results of Copula model

表5 Copula函數(shù)擬合優(yōu)度檢驗Table 5 Goodness-of-fit test for Copula function

圖4 Gumbel Copula模型概率密度Fig.4 Probability density of Gumbel Copula model

將式(38)和式(39)代入到式(40)中，便可得到迎角極值和俯仰角速率極值聯(lián)合分布概率密度函數(shù)(,)。

4 基于吸引域理論和二元極值Copula模型的飛行風險概率計算

綜上所述，背景飛機結冰條件下飛行風險量化評估流程如圖5所示。

圖5 結冰條件下飛行風險量化評估流程Fig.5 Process of flight risk quantitative assessment under icing conditions

4.1 典型飛行想定

結冰相關安全關鍵飛行參數(shù)有迎角與飛行速度等，當駕駛員沒有意識到結冰的嚴重程度時，一旦操縱不當導致這些運動參數(shù)超過其安全邊界，便有可能發(fā)生飛行事故。將所有可能的結冰致災因素考慮到風險評估中會提高風險預測的準確性，然而這會使工作變得異常繁瑣，而且沒有必要?？紤]到迎角通常是邊界保護問題研究中首要考慮的因素，本文在縱向通道上，以典型的駕駛員操縱——升降舵脈沖作為輸入，將飛行風險事件的發(fā)生定義為飛機狀態(tài)超出基于迎角和俯仰角速率的二維穩(wěn)定邊界，以此為基準，對結冰后的飛行動力學仿真結果進行分析評估。

之所以選擇升降舵脈沖，一方面該操縱被廣泛用于結冰后的參數(shù)辨識中，能夠充分激發(fā)飛機的動力學響應,另一方面在駕駛員意識到飛機因結冰高度降低時，由于操縱者對結冰帶來的飛行包線萎縮并沒有充分的認識，很可能會采用近似于升降舵脈沖的操縱方式來修正高度。

4.2 不同結冰嚴重程度下的飛行風險概率

以飛行狀態(tài)超出穩(wěn)定邊界為飛行風險判定條件，采用Gumbel Copula模型計算結冰條件下的飛行風險概率：

(41)

對于背景飛機，假定飛機以=2 000 m，=120 m/s的初始狀態(tài)保持水平直線飛行，當結冰嚴重程度因子=0.1時，吸引域1為

1={,∈|0621 8(-0119 7)-

0274 4(-0119 7)+0608 9≤0046}

(42)

將1及背景飛機常規(guī)動力學模型的穩(wěn)定與不穩(wěn)定相軌跡投影到-相平面圖上，如圖6所示。辨識Gumbel Copula模型中的未知參數(shù)，結果如表6所示。

表6 結冰嚴重程度η=0.1時Gumbel模型參數(shù)辨識結果Table 6 Parameter identification results of Gumbel model when icing severity η=0.1

圖6 結冰嚴重程度η=0.1時α-q相平面圖與吸引域Fig.6 Phase plane plot of α-q and region of attraction when icing severity η=0.1

將式(42)代入式(41)中，計算出當前結冰嚴重程度條件下的飛行風險概率為

=887×10

(43)

該飛機處于弱結冰(Trace)狀態(tài)，按照航空信息手冊中的說明，除非飛機暴露于結冰區(qū)中較長時間(大于1 h)，駕駛員在穿越結冰區(qū)時無需開啟防/除冰設備。

當結冰嚴重程度因子=0.2時，吸引域2為

2={,∈|0601 2(-0127 9)-

0301 7(-0127 9)+0638 3≤0028 2}

(44)

同樣，將2及背景飛機常規(guī)動力學模型的穩(wěn)定與不穩(wěn)定相軌跡投影到-相平面圖上，如圖7所示。辨識Gumbel Copula模型中的未知參數(shù)，結果如表7所示。

圖7 結冰嚴重程度η=0.2時α-q相平面圖與吸引域Fig.7 Phase plane plot of α-q and region of attraction when icing severity η=0.2

表7 結冰嚴重程度η=0.2時Gumbel模型參數(shù)辨識結果Table 7 Parameter identification results of Gumbel model when icing severity η=0.2

當前結冰嚴重程度下的飛行風險概率為

=126×10

(45)

該飛機處于輕度結冰(Light)狀態(tài)，按照航空信息手冊中的說明，飛機長時間在目標航跡區(qū)遭遇該嚴重程度的結冰時，駕駛員需要間歇開啟防/除冰設備。

當結冰嚴重程度因子=0.3時，吸引域3為

3={,∈|0542 1(-0147 9)-

0252 7(-0147 9)+0857 8≤0014 4}

(46)

同樣，將3及背景飛機常規(guī)動力學模型的穩(wěn)定與不穩(wěn)定相軌跡投影到-相平面圖上，如圖8所示。辨識Gumbel Copula模型中的未知參數(shù)，結果如表8所示。

圖8 結冰嚴重程度η=0.3時α-q相平面圖與吸引域Fig.8 Phase plane plot of α-q and region of attraction when icing severity η=0.3

表8 結冰嚴重程度η=0.3時Gumbel模型參數(shù)辨識結果Table 8 Parameter identification results of Gumbel model when icing severity η=0.3

當前結冰嚴重程度條件下的飛行風險概率為

=742×10

(47)

該飛機處于中度結冰(Moderate)狀態(tài)，按照航空信息手冊中的說明，飛機一旦在目標航跡區(qū)遭遇該嚴重程度的結冰，發(fā)生飛行風險的可能性較高，在運行過程中，駕駛員必須一直開啟防/除冰設備或者改變路線駛離結冰區(qū)域。

由于在進行多項式擬合、風險概率計算的每一步都會引入誤差，最終的計算結果很難用一個準確的誤差范圍去衡量,本文通過對每一步計算結果的檢驗來進行誤差控制。誤差的引入主要有以下兩個方面：

1) 多項式擬合帶來的誤差。雖然多項式模型轉化過程中引入了誤差，但圖2所示的相平面對比圖說明誤差對結果影響不大，究其原因，基于平方和理論的吸引域計算方法是一種偏保守的方法，存在一定的安全裕度，該方法在解決F/A-18飛機控制器對落葉飄模態(tài)敏感性分析上得到了應用，驗證了該方法的有效性。

2) 應用極值理論帶來的誤差。雖然極值參數(shù)分布的隨機性帶來的誤差不可避免，但只要通過文中所述的檢驗方法，即可認為誤差是可接受的。

5 結 論

本文首先建立了飛機縱向動力學模型和結冰影響模型，基于吸引域理論，得到了不同結冰程度條件下的飛機二維吸引域及穩(wěn)定邊界；建立典型的人-機-環(huán)系統(tǒng)模型，運用蒙特卡羅仿真方法，得到了迎角極值和俯仰角速率極值；根據(jù)極值理論，建立并確定了適用于結冰條件下飛行風險量化評估的二元極值Copula模型，計算了不同結冰嚴重程度下的風險概率并對駕駛員操縱提出指導，主要結論如下：

1) 吸引域方法適用于計算飛機飛行的穩(wěn)定邊界，計算結果直觀明了。隨著結冰嚴重程度的增加，飛行的穩(wěn)定邊界縮小，導致發(fā)生飛行事故的風險增加。

2) 基于吸引域和二元極值理論的結冰飛機飛行風險量化評估方法，能夠用來定量地計算不同結冰嚴重程度下的飛行風險，進而指導駕駛員進行合理的操縱，從而規(guī)避飛行危險。

附錄A

運用四元數(shù)法構建飛機動力學模型：

(A1)

(A2)

(A3)

(A4)

同時：

(A5)

(A6)

(A7)

采用四元數(shù)法求解運動方程，能夠有效避免出現(xiàn)奇點，同時能夠減少三角函數(shù)的計算，提高運行效率。四元數(shù)滿足平方和為1，但其物理意義不夠明確，可依據(jù)四元數(shù)計算出歐拉角，飛機姿態(tài)角、、可通過式(A8)～式(A10)計算。在反求歐拉角的過程中，可能出現(xiàn)分母為零的情況，可以通過巧妙編程解決，不影響運動方程求解連續(xù)性，同時為減小迭代計算產(chǎn)生的累積誤差，每步仿真均將四元數(shù)進行歸一化處理。

(A8)

=arcsin[2(-)]

(A9)

(A10)

附錄B

-設序列,,…,為獨立同分布的隨機變量，當存在常數(shù)列{>0}、{}時，滿足：

(B1)

則()是非退化分布函數(shù)，且()必屬于下列3種類型分布之一：

1)Ⅰ型分布

()=exp(-e-) -∞<<+∞

(B2)

2) Ⅱ型分布

(B3)

3) Ⅲ型分布

(B4)

上述3種分布分別稱為Gumbel分布、Frechet分布、Weibull分布，可統(tǒng)稱為極值分布,和稱為規(guī)范化常數(shù)。當=1時，(,1)和(,1)分別稱為標準Frechet分布和標準Weibull分布。

進一步引入位置參數(shù)和尺度參數(shù)，則可以將上述3種類型的極值分布形式統(tǒng)一為

(B5)

附錄C

表C1 迎角極值樣本
Table C1 Sample of extreme angle of attack (°)

序號樣本值1-58.7788.8948.9888.9979.0186-109.0219.0399.0589.0829.09311-159.1199.1539.1829.2019.22816-209.2519.2729.3059.3179.33221-259.3339.3349.3419.3469.35626-309.3629.3889.4039.4329.45031-359.4749.4989.5279.6019.64336-409.6829.7019.7239.8329.89541-459.96210.00810.05910.16510.19746-5010.23110.38910.41810.48310.532

表C2 俯仰角極值樣本
Table C2 Sample of extreme pitch angle (°)

序號樣本值1-55.7965.8065.9045.9485.9776-105.9815.9856.0116.0176.02511-156.0466.0756.1056.1516.16316-206.2216.2256.2346.2416.25721-256.2596.2606.2666.2856.29126-306.3176.3306.3726.3746.39831-356.4336.4886.5136.5966.60136-406.6236.6896.7116.7936.84241-456.9086.9717.1327.1747.19846-507.2557.3497.3567.4677.896

附錄D

假設聯(lián)合分布函數(shù)具有邊緣分布函數(shù)()，()，…，()，則?-Copula函數(shù)，滿足：

(,,…,)=((),(),…,())

(D1)

若()，()，…，()連續(xù)，則唯一確定。