基于吸引域與二元極值理論的結(jié)冰飛機(jī)飛行風(fēng)險(xiǎn)量化評(píng)估

伍強(qiáng)，徐浩軍，裴彬彬,*，魏揚(yáng)，張楊

1. 空軍工程大學(xué) 航空工程學(xué)院，西安 710038

2. 中國(guó)人民解放軍 95633部隊(duì)，成都 611530

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

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

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

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

1 模型建立

1.1 飛機(jī)縱向動(dòng)力學(xué)模型

飛機(jī)縱向非線性動(dòng)力學(xué)模型可表示為

(1)

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

(2)

1.2 結(jié)冰影響模型

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

=(1+)

(3)

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

1.3 飛機(jī)本體模型

飛機(jī)本體非線性動(dòng)力學(xué)模型可表示為

(4)

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

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

(5)

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

=[,,,]

(6)

運(yùn)用四元數(shù)法構(gòu)建飛機(jī)動(dòng)力學(xué)模型具體過程見附錄A。

1.4 駕駛員操縱模型

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

(7)

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

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

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

考慮多項(xiàng)式非線性自治動(dòng)力學(xué)系統(tǒng)：

(8)

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

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

(9)

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

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

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

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

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

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

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

:={∈|()≤}

(10)

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

(11)

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

:={∈|()=}

(12)

在進(jìn)行估算時(shí)，可假定吸引域?yàn)橐粰E球體，對(duì)于二維問題，可假設(shè)：

()=N=()+()

(13)

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

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

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

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

將建立起的各多項(xiàng)式模型代入常規(guī)動(dòng)力學(xué)模型中，進(jìn)而推導(dǎo)出飛機(jī)動(dòng)力學(xué)模型的多項(xiàng)式表達(dá)式。

為便于計(jì)算，通過坐標(biāo)變換將飛機(jī)多項(xiàng)式動(dòng)力學(xué)模型的平衡點(diǎn)移動(dòng)到原點(diǎn)，并去除不必要的項(xiàng)(系數(shù)很小或階數(shù)很高)。

基于SOS理論計(jì)算出飛機(jī)在原點(diǎn)附近的吸引域。

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

2.2 飛機(jī)縱向多項(xiàng)式動(dòng)力學(xué)模型

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

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

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

=++

(14)

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

=7153 1×10-

2447 4×10+0027 4

(15)

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

(16)

式中：()是只與迎角有關(guān)的基本氣動(dòng)參數(shù)，選用一元四階多項(xiàng)式進(jìn)行擬合，得到其多項(xiàng)式表達(dá)形式為

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

5137 2+0068 9

(17)

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

0034 6+0454 7+0002 5

(18)

(19)

升力系數(shù)各項(xiàng)的多項(xiàng)式擬合結(jié)果與真實(shí)值之間的對(duì)比如圖1所示。

圖1 升力系數(shù)各項(xiàng)擬合值與真實(shí)值之間的對(duì)比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ù)的多項(xiàng)式表達(dá)式。利用同樣的方法，可得到阻力系數(shù)、俯仰力矩系數(shù)的多項(xiàng)式表達(dá)式。在進(jìn)行擬合計(jì)算時(shí)，非線性動(dòng)力學(xué)模型中各角度均以(°)為單位。同時(shí)，需要指出的是，該多項(xiàng)式表達(dá)式成立的參數(shù)范圍需滿足飛機(jī)氣動(dòng)參數(shù)插值表中給出的插值范圍限制。對(duì)于該型飛機(jī)而言，各參數(shù)的有效范圍為：∈[-10°, 60°],∈[-30, 30](°)/s,∈[-30°, 20°]，∈[80,150] m/s。

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

(21)

通過上述對(duì)常規(guī)動(dòng)力學(xué)模型中非線性項(xiàng)的多項(xiàng)式轉(zhuǎn)化，可將常規(guī)的飛機(jī)縱向非線性動(dòng)力學(xué)模型轉(zhuǎn)化為多項(xiàng)式形式的微分方程組。

2.3 基于SOS理論的吸引域計(jì)算

2.3.1 多項(xiàng)式動(dòng)力學(xué)模型簡(jiǎn)化

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

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

最終將常規(guī)的飛機(jī)非線性動(dòng)力學(xué)模型(式(1))按照上述方法轉(zhuǎn)化得到的多項(xiàng)式短周期動(dòng)力學(xué)模型為

0004 3+0056 5-0537 8+0891 2

(22)

1423 2-0783 5

(23)

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

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

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

()=8207+3647 6

(24)

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

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

0538 9≤0580 5}

(25)

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

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

0538 9=0580 5}

(26)

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

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

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

(27)

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

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

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

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

3.1 一維極值分布模型

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

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

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

(28)

(29)

(30)

(31)

(32)

(33)

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

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

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

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

圖3 QQ檢驗(yàn)圖Fig.3 QQ inspection chart

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

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

3.3 二元飛行安全關(guān)鍵參數(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)

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

(40)

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

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

表5 Copula函數(shù)擬合優(yōu)度檢驗(yàn)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模型的飛行風(fēng)險(xiǎn)概率計(jì)算

綜上所述，背景飛機(jī)結(jié)冰條件下飛行風(fēng)險(xiǎn)量化評(píng)估流程如圖5所示。

圖5 結(jié)冰條件下飛行風(fēng)險(xiǎn)量化評(píng)估流程Fig.5 Process of flight risk quantitative assessment under icing conditions

4.1 典型飛行想定

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

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

4.2 不同結(jié)冰嚴(yán)重程度下的飛行風(fēng)險(xiǎn)概率

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

(41)

對(duì)于背景飛機(jī)，假定飛機(jī)以=2 000 m，=120 m/s的初始狀態(tài)保持水平直線飛行，當(dāng)結(jié)冰嚴(yán)重程度因子=0.1時(shí)，吸引域1為

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

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

(42)

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

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

圖6 結(jié)冰嚴(yán)重程度η=0.1時(shí)α-q相平面圖與吸引域Fig.6 Phase plane plot of α-q and region of attraction when icing severity η=0.1

將式(42)代入式(41)中，計(jì)算出當(dāng)前結(jié)冰嚴(yán)重程度條件下的飛行風(fēng)險(xiǎn)概率為

=887×10

(43)

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

當(dāng)結(jié)冰嚴(yán)重程度因子=0.2時(shí)，吸引域2為

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

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

(44)

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

圖7 結(jié)冰嚴(yán)重程度η=0.2時(shí)α-q相平面圖與吸引域Fig.7 Phase plane plot of α-q and region of attraction when icing severity η=0.2

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

當(dāng)前結(jié)冰嚴(yán)重程度下的飛行風(fēng)險(xiǎn)概率為

=126×10

(45)

該飛機(jī)處于輕度結(jié)冰(Light)狀態(tài)，按照航空信息手冊(cè)中的說明，飛機(jī)長(zhǎng)時(shí)間在目標(biāo)航跡區(qū)遭遇該嚴(yán)重程度的結(jié)冰時(shí)，駕駛員需要間歇開啟防/除冰設(shè)備。

當(dāng)結(jié)冰嚴(yán)重程度因子=0.3時(shí)，吸引域3為

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

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

(46)

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

圖8 結(jié)冰嚴(yán)重程度η=0.3時(shí)α-q相平面圖與吸引域Fig.8 Phase plane plot of α-q and region of attraction when icing severity η=0.3

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

當(dāng)前結(jié)冰嚴(yán)重程度條件下的飛行風(fēng)險(xiǎn)概率為

=742×10

(47)

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

由于在進(jìn)行多項(xiàng)式擬合、風(fēng)險(xiǎn)概率計(jì)算的每一步都會(huì)引入誤差，最終的計(jì)算結(jié)果很難用一個(gè)準(zhǔn)確的誤差范圍去衡量,本文通過對(duì)每一步計(jì)算結(jié)果的檢驗(yàn)來進(jìn)行誤差控制。誤差的引入主要有以下兩個(gè)方面：

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

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

5 結(jié) 論

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

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

2) 基于吸引域和二元極值理論的結(jié)冰飛機(jī)飛行風(fēng)險(xiǎn)量化評(píng)估方法，能夠用來定量地計(jì)算不同結(jié)冰嚴(yán)重程度下的飛行風(fēng)險(xiǎn)，進(jìn)而指導(dǎo)駕駛員進(jìn)行合理的操縱，從而規(guī)避飛行危險(xiǎn)。

附錄A

運(yùn)用四元數(shù)法構(gòu)建飛機(jī)動(dòng)力學(xué)模型：

(A1)

(A2)

(A3)

(A4)

同時(shí)：

(A5)

(A6)

(A7)

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

(A8)

=arcsin[2(-)]

(A9)

(A10)

附錄B

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

(B1)

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

1)Ⅰ型分布

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

(B2)

2) Ⅱ型分布

(B3)

3) Ⅲ型分布

(B4)

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

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

(B5)

附錄C

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

序號(hào)樣本值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 (°)

序號(hào)樣本值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

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

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

(D1)

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