王健名， 徐浩軍， 薛源， 王小龍， 李哲

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

基于極值理論的平尾結(jié)冰飛行風(fēng)險(xiǎn)評估

王健名， 徐浩軍*， 薛源， 王小龍， 李哲

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

提出了結(jié)合極值理論與Copula模型來量化評估平尾結(jié)冰條件下飛行風(fēng)險(xiǎn)概率的方法。通過建立人-機(jī)-環(huán)復(fù)雜系統(tǒng)模型，對平尾在進(jìn)近與著陸過程中的結(jié)冰情形進(jìn)行仿真，采用蒙特卡羅法提取平尾結(jié)冰極值參數(shù)，驗(yàn)證了所提取極值參數(shù)符合一維廣義極值(GEV)分布，根據(jù)飛行風(fēng)險(xiǎn)的定義和相關(guān)安全性準(zhǔn)則，建立了平尾結(jié)冰飛行風(fēng)險(xiǎn)發(fā)生的判定條件，計(jì)算得出一維極值飛行風(fēng)險(xiǎn)概率；在此基礎(chǔ)上選取Copula模型來描述二維極值參數(shù)的相關(guān)性，對多種Copula模型的未知參數(shù)進(jìn)行辨識，通過擬合優(yōu)度檢驗(yàn)對精度進(jìn)行驗(yàn)證，得出Joe Copula模型對二維極值分布的描述最為準(zhǔn)確，運(yùn)用Joe Copula模型計(jì)算出二維極值飛行風(fēng)險(xiǎn)概率，有效解決了一維極值具有的局限性。所提方法對飛行安全評估等理論有一定參考價(jià)值，能為平尾結(jié)冰飛行事故的預(yù)防提供分析和檢驗(yàn)依據(jù)。

極值理論； Copula模型； 平尾結(jié)冰； 飛行風(fēng)險(xiǎn)概率； 蒙特卡羅法； 參數(shù)辨識

飛機(jī)飛行時(shí)，常常會產(chǎn)生平尾結(jié)冰現(xiàn)象，導(dǎo)致平尾氣動(dòng)特性被破壞，平尾負(fù)升力減小，配平能力下降，升降舵操縱效率降低。當(dāng)平尾結(jié)冰發(fā)生在進(jìn)近與著陸過程中時(shí)，飛機(jī)極易發(fā)生平尾失速現(xiàn)象，甚至導(dǎo)致飛行事故發(fā)生。目前，平尾結(jié)冰已成為世界較為關(guān)注的安全隱患問題，評估在平尾結(jié)冰條件下的飛行風(fēng)險(xiǎn)概率顯得極其重要。

國外最早對飛機(jī)平尾結(jié)冰現(xiàn)象進(jìn)行研究。包括風(fēng)洞試驗(yàn)和飛行測試，如1994年NASA和FAA聯(lián)合研究的平尾結(jié)冰項(xiàng)目Tailplane Icing Program(TIP)，其以DHC-6為試驗(yàn)機(jī)型，進(jìn)行了一系列平尾結(jié)冰飛行試驗(yàn)和風(fēng)洞試驗(yàn)，得到了大量的試驗(yàn)數(shù)據(jù)[1-2]；通過計(jì)算流體力學(xué)(CFD)的應(yīng)用，進(jìn)行水滴運(yùn)動(dòng)和撞擊特性計(jì)算、冰形預(yù)測以及分析平尾結(jié)冰后的氣動(dòng)參數(shù)變化[3-5]；對平尾結(jié)冰進(jìn)行飛行仿真，分析平尾結(jié)冰對飛機(jī)動(dòng)力學(xué)特性和飛行品質(zhì)的影響[6-7]；進(jìn)行相應(yīng)的結(jié)冰控制和防、除冰方法與裝置的研究[8-10]，如1998年首次提出了飛機(jī)智能防冰系統(tǒng)(Smart Icing System， SIS)等。

與國外相比，國內(nèi)關(guān)于平尾結(jié)冰方面的研究起步較晚，主要在以下方面進(jìn)行了一定的相關(guān)性研究：對平尾結(jié)冰飛機(jī)縱向氣動(dòng)參數(shù)進(jìn)行辨識，定量地分析了平尾結(jié)冰對飛機(jī)縱向氣動(dòng)參數(shù)的影響[11]；通過研究平尾結(jié)冰導(dǎo)致的飛行事故，分析了平尾結(jié)冰對平尾氣動(dòng)特性的影響，給出預(yù)防平尾失速的方法[12]；預(yù)測結(jié)冰外形、結(jié)冰邊界保護(hù)以及防除冰[13-15]等。

但綜合國內(nèi)外來看，關(guān)于平尾結(jié)冰條件下的飛行風(fēng)險(xiǎn)概率評估方面的研究極少。而國軍標(biāo)GJB900-90[16]、SAE-ARP-4761[17]、MIL-HDBK-516B[18]及MIL-STD-882D[19]等安全性規(guī)范考慮飛行過程中的不確定性和隨機(jī)性因素較少，較難對多因素耦合、非線性較強(qiáng)的飛行條件進(jìn)行飛行風(fēng)險(xiǎn)量化預(yù)測，尤其是在進(jìn)近與著陸過程中平尾結(jié)冰的條件下。因此，迫切需要建立有效的理論與方法來定量評估平尾結(jié)冰條件下的飛行風(fēng)險(xiǎn)概率。鑒于此，本文提出了一種基于人-機(jī)-環(huán)復(fù)雜系統(tǒng)建模仿真和極值理論[20]提取相應(yīng)極值參數(shù)，利用Copula模型來量化評估飛行風(fēng)險(xiǎn)概率的方法。此方法對研究平尾結(jié)冰引起的飛行安全和適航性問題具有重要意義。圖1為本文飛行風(fēng)險(xiǎn)評估流程圖。

圖1 飛行風(fēng)險(xiǎn)評估流程Fig.1 Flowchart of flight risk evaluation

1 平尾結(jié)冰后人-機(jī)-環(huán)復(fù)雜系統(tǒng)建模

平尾結(jié)冰的內(nèi)外部影響因素具有復(fù)雜隨機(jī)性的特點(diǎn)：一是復(fù)雜的非線性物理特性，主要指駕駛員和飛機(jī)在解析模型中的非線性；二是復(fù)雜的隨機(jī)性，主要指外部天氣變化、駕駛員操縱行為差異以及航空設(shè)備故障發(fā)生等。考慮到駕駛員、飛機(jī)和平尾結(jié)冰環(huán)境之間的相互耦合及相互作用的復(fù)雜關(guān)系，本文分別建立了駕駛員模型、飛機(jī)運(yùn)動(dòng)模型和平尾結(jié)冰參量模型，實(shí)現(xiàn)平尾結(jié)冰動(dòng)力學(xué)仿真。

1.1 駕駛員模型

平尾結(jié)冰后，駕駛員的操縱對飛機(jī)影響極大，因此駕駛員模型是評估平尾結(jié)冰飛行風(fēng)險(xiǎn)需考慮的重要因素。本文將駕駛員的操作分為3個(gè)階段。當(dāng)t

圖2 駕駛員補(bǔ)償模型Fig.2 Compensation model of pilot

則駕駛員模型的輸入輸出關(guān)系可以表示為

(1)

平尾結(jié)冰會導(dǎo)致平尾負(fù)升力減小，升降舵效率降低，使得飛機(jī)配平能力減弱，進(jìn)而會使飛機(jī)轉(zhuǎn)入俯沖狀態(tài)，若駕駛員不及時(shí)調(diào)整升降舵，飛機(jī)可能會發(fā)生飛行事故。以升降舵為例，取俯仰角為信號，建立升降舵運(yùn)動(dòng)的非線性模型(副翼運(yùn)動(dòng)的非線性模型與之相似，在此不再贅述)：

(2)

式中：t0為駕駛員開始操縱升降舵的時(shí)刻；Δt為偏轉(zhuǎn)所需時(shí)間；δe0為初始升降舵偏角；Vδe為升降舵偏轉(zhuǎn)速率；tp為延遲時(shí)間，主要包括駕駛員操縱過程中的固有延時(shí)特性、肌肉神經(jīng)滯后時(shí)間以及飛行參數(shù)傳遞和處理過程中的滯后時(shí)間等。

假設(shè)駕駛員操縱升降舵的偏轉(zhuǎn)速率Vδe服從均勻分布，其概率密度函數(shù)為

(3)

式中：a和b分別為升降舵偏轉(zhuǎn)速率的最小值和最大值。

1.2 基于四元數(shù)法的飛機(jī)運(yùn)動(dòng)模型

基于四元數(shù)法構(gòu)建飛機(jī)本體的六自由度運(yùn)動(dòng)模型為

(4)

(5)

(6)

(7)

(8)

(9)

c1=(Iyy-Izz)Izz-IxzIxz/λI

c2=(Ixx-Iyy+Izz)Ixz/λI

c3=Izz/λI

c4=Ixz/λI

c5=(Izz-Ixx)/Iyy

c6=Ixz/Iyy

c7=1/Iyy

c8=[(Ixx-Iyy)Ixx+IxzIxz]/λI

c9=Ixx/λI

λI=IxxIzz-IxzIxz

1.3 平尾結(jié)冰參量模型

采用平尾結(jié)冰后的氣動(dòng)參數(shù)模型[6]：

(10)

(11)

式中：KCA為氣動(dòng)導(dǎo)數(shù)對結(jié)冰的敏感性，只與飛機(jī)有關(guān)而與氣象條件無關(guān)，對于給定的飛機(jī)是常值。η和ηice相似，只是采用實(shí)際飛機(jī)平尾氣動(dòng)弦長和相應(yīng)飛行速度進(jìn)行計(jì)算。它是關(guān)于大氣溫度、液態(tài)水含量、水滴直徑、結(jié)冰遭遇時(shí)間等的函數(shù)，只與氣象條件有關(guān)，與具體飛機(jī)無關(guān)。

由式(10)和式(11)可得

C(A)iced=1+ηKCAC(A)

(12)

此結(jié)冰氣動(dòng)導(dǎo)數(shù)的計(jì)算公式是通過NASA對雙水獺飛機(jī)進(jìn)行結(jié)冰飛行試驗(yàn)研究，由Bragg等提出的，通過結(jié)冰后氣動(dòng)導(dǎo)數(shù)下降比例和結(jié)冰嚴(yán)重影響程度η表征飛機(jī)受結(jié)冰影響。文獻(xiàn)[21]通過試飛驗(yàn)證了此模型在Cessna 208B飛機(jī)上同樣適用，因此模型具有較高的適用性。本文亦采用此模型。

2 基于蒙特卡羅法仿真的二維極值參數(shù)選取

本文進(jìn)行的平尾結(jié)冰飛行風(fēng)險(xiǎn)概率評估是基于飛機(jī)在進(jìn)近與著陸過程中，平尾發(fā)生結(jié)冰導(dǎo)致平尾臨界負(fù)失速迎角減小(平尾臨界失速迎角增大)，平尾負(fù)升力降低，飛機(jī)進(jìn)入俯沖狀態(tài)，進(jìn)而引起飛行風(fēng)險(xiǎn)來展開研究的。

研究此科目內(nèi)外部影響因素需要大量數(shù)據(jù)，試飛與人在回路地面試驗(yàn)都無法達(dá)到此數(shù)據(jù)量；并且飛機(jī)在進(jìn)近與著陸過程中的平尾結(jié)冰飛行風(fēng)險(xiǎn)評估屬于高風(fēng)險(xiǎn)科目，飛行風(fēng)險(xiǎn)大，需模擬的外部環(huán)境等條件又十分困難。因此，如何利用有限數(shù)據(jù)量對平尾結(jié)冰飛行風(fēng)險(xiǎn)進(jìn)行量化評估是一個(gè)難點(diǎn)。本文通過建立平尾結(jié)冰后的人-機(jī)-環(huán)復(fù)雜系統(tǒng)模型，利用極值理論，基于蒙特卡羅法提取相應(yīng)飛行參數(shù)極值，能較好地解決上述問題。

2.1 飛行極值參數(shù)提取

仿真共有3個(gè)數(shù)據(jù)庫，分別是駕駛員模型數(shù)據(jù)庫、平尾結(jié)冰數(shù)據(jù)庫和外部環(huán)境數(shù)據(jù)庫。駕駛員模型數(shù)據(jù)庫：駕駛員個(gè)體具有不同的操縱習(xí)慣、心理狀態(tài)和身體素質(zhì)，使得駕駛員在應(yīng)對新環(huán)境時(shí)表現(xiàn)出操縱行為的差異性，通過采集飛行員的實(shí)驗(yàn)操縱數(shù)據(jù)，辨識了駕駛員模型中的未知參數(shù)，主要包括駕駛員操縱過程中的固有延時(shí)特性，肌肉神經(jīng)滯后時(shí)間以及飛行參數(shù)傳遞和處理過程中的滯后時(shí)間等，從而建立駕駛員操縱行為數(shù)據(jù)庫；平尾結(jié)冰數(shù)據(jù)庫：考慮到本文飛機(jī)氣動(dòng)結(jié)構(gòu)與DHC-6相近，故仿真利用DHC-6平尾結(jié)冰數(shù)據(jù)；外部環(huán)境數(shù)據(jù)庫：主要包括空氣密度、當(dāng)?shù)貧鈮?、離散突風(fēng)、風(fēng)切變狀況、雨雪狀況和空氣濕度等。

基于蒙特卡羅法提取出相應(yīng)飛行極值參數(shù)過程如圖3所示，具體步驟如下：

1) 設(shè)定飛機(jī)初始飛行狀態(tài)。初始飛行高度、速度、姿態(tài)角、結(jié)冰時(shí)間等。

2) 將由天氣狀況得到的云層液態(tài)水含量(LWC)、平均有效水滴直徑(MVD)、溫度T，在平尾結(jié)冰數(shù)據(jù)庫中進(jìn)行蒙特卡羅抽樣，提取平尾結(jié)冰數(shù)據(jù)的隨機(jī)變量值，將其作為檢索條件，結(jié)合設(shè)定的結(jié)冰時(shí)間，依據(jù)統(tǒng)計(jì)數(shù)據(jù)從平尾結(jié)冰數(shù)據(jù)庫中提取出平尾結(jié)冰嚴(yán)重程度η，其值是一個(gè)常值。

3) 從外部環(huán)境數(shù)據(jù)庫中進(jìn)行蒙特卡羅抽樣，提取外部環(huán)境影響因素的隨機(jī)變量值。

4) 從駕駛員模型數(shù)據(jù)庫中進(jìn)行蒙特卡羅抽樣，提取駕駛員操縱行為參數(shù)的隨機(jī)變量值。

5) 進(jìn)行第i次飛行仿真計(jì)算。

飛機(jī)本體方程為基于四元數(shù)的六自由度方程，采用四階龍格庫塔微分算法。實(shí)時(shí)仿真計(jì)算機(jī)通過相應(yīng)飛控系統(tǒng)轉(zhuǎn)化的C代碼進(jìn)行實(shí)時(shí)仿真，仿真時(shí)間步長為0.02 s。

6) 記錄第i次計(jì)算結(jié)果中的飛行參數(shù)，提取飛行參數(shù)極值，存入數(shù)據(jù)庫。

7)i=i+1，返回到步驟2)，循環(huán)迭代到i=n結(jié)束。

n為蒙特卡羅計(jì)算次數(shù)，n越大越能反應(yīng)極值的統(tǒng)計(jì)特性，但n太大會又會導(dǎo)致計(jì)算任務(wù)增加。經(jīng)驗(yàn)證當(dāng)n>2 000時(shí)，計(jì)算結(jié)果趨于穩(wěn)定，因此本文選取n=2 000。

圖3 飛行極值參數(shù)的提取過程Fig.3 Process of extracting flight extreme parameters

圖4以第56次飛行仿真計(jì)算為例，表征說明飛行仿真及參數(shù)極值提取的過程，仿真是在初始高度H=1 000 m(機(jī)場高度500 m)，速度V=75 m/s，飛機(jī)遭遇平尾結(jié)冰時(shí)間為5 min，襟翼角為35° 時(shí)進(jìn)行的，得到的平尾結(jié)冰嚴(yán)重程度η=0.4。從仿真曲線可以看出，飛機(jī)平尾發(fā)生結(jié)冰后，法向過載急劇增大，易超出其飛行包線；平尾負(fù)迎角(Negative Tailplane Angle of Attack，Negative Tail AOA)也增大，易達(dá)到平尾臨界負(fù)迎角，甚至導(dǎo)致平尾失速的發(fā)生，因此確定對平尾結(jié)冰飛行風(fēng)險(xiǎn)影響最大的兩個(gè)參數(shù)為平尾負(fù)迎角和法向過載。提取本次仿真中的平尾負(fù)迎角極值αtailmax=8.307 0°，法向過載極值nzmax=1.410 1。

圖4 第56次迭代的飛行參數(shù)Fig.4 Flight parameters in the 56th iteration

2.2 極值參數(shù)可信度驗(yàn)證

圖5為i≤100所提取的極值樣本與地面試驗(yàn)數(shù)據(jù)駕駛員在回路地面飛行實(shí)時(shí)仿真試驗(yàn)數(shù)據(jù)的散點(diǎn)圖和QQ圖(若呈一條直線，則兩極值樣本分布近似相同)。從QQ圖可以看出平尾負(fù)迎角和法向過載的兩種極值樣本分位數(shù)曲線接近為直線，從散點(diǎn)圖可以看出兩種極值樣本的數(shù)據(jù)點(diǎn)分布趨于一致，表明所提取的仿真數(shù)據(jù)和駕駛員試驗(yàn)數(shù)據(jù)的分布類型近似相同。表1為相關(guān)性系數(shù)和K-S檢驗(yàn)計(jì)算結(jié)果。K-S值均小于0.1，P值大于0.05，即在95%的置信水平下能通過檢驗(yàn)；相關(guān)性系數(shù)均大于0.9，表明平尾負(fù)迎角和法向過載的兩種極值樣本有強(qiáng)線性相關(guān)性。因此，認(rèn)為所提取的仿真數(shù)據(jù)和駕駛員試驗(yàn)數(shù)據(jù)屬于同一種分布類型，具有較高的準(zhǔn)確性和可信度，可以當(dāng)作平尾結(jié)冰風(fēng)險(xiǎn)評估的樣本數(shù)據(jù)。

圖5 提取極值樣本與地面試驗(yàn)數(shù)據(jù)的QQ圖和散點(diǎn)圖Fig.5 QQ and scatter plots of simulated and ground test extreme value samples

表1 極值樣本擬合優(yōu)度檢驗(yàn)Table 1 Goodness-of-fit test of extreme value samples

ExtremevaluesampleCorrelationcoefficientK?SvaluePαtailmax0．90870．07000．9610nzmax0．95720．08000．8938

2.3 平尾結(jié)冰飛行風(fēng)險(xiǎn)判據(jù)

(13)

3 一維極值參數(shù)分布

表2為參數(shù)極值的統(tǒng)計(jì)特性，可以看出峰度系數(shù)均大于3，表明所提取的極值樣本較集中，比正態(tài)分布有更長的尾部；偏度系數(shù)均大于0，表明分布類型在右側(cè)具有較長的尾部。因此，可以得出結(jié)論：選取的極值參數(shù)分布具有較明顯的厚尾特性。而極值理論適合對具有厚尾特性的低頻高危飛行風(fēng)險(xiǎn)案例進(jìn)行研究，因此針對此種分布類型，本文采用極值理論對其進(jìn)行描述。與一維極值分布不同的是，對于二維參數(shù)極值分布，需要研究參數(shù)之間的相關(guān)性，采用Copula理論[22-24]對涉及到多維極值參數(shù)的平尾結(jié)冰飛行風(fēng)險(xiǎn)概率進(jìn)行評估具有較高準(zhǔn)確性。由Copula理論多元極值分布定義[25]可知：F是多元極值分布函數(shù)，則F的一維邊緣分布屬于GEV(Generalized Extreme Value)分布族。因此，研究參數(shù)之間的相關(guān)性之前，要先驗(yàn)證一維極值參數(shù)符合GEV分布。

表2 極值樣本統(tǒng)計(jì)量Table 2 Statistics of extreme value samples

3.1 一維極值參數(shù)辨識

極值理論中較為廣泛應(yīng)用的有GEV分布(式(14))、正態(tài)(Normal)分布(式(15))、對數(shù)正態(tài)(Lognormal)分布(式(16))、威布(Weibull)分布(式(17))、指數(shù)(Exponential)分布(式(18))和EV(Extreme Value)分布(式(19))。分布函數(shù)為

(14)

(15)

(16)

(17)

F(x;μ)=1-exp-x/μ

(18)

(19)

式中：ξ、μ、σ均為分布參數(shù)。

基于所提取飛行參數(shù)極值對以上函數(shù)中的未知參數(shù)進(jìn)行辨識，結(jié)果如表3所示。

表3不同分布模型的參數(shù)辨識

Table3Identifiedparametersofdifferentdistributionmodel

ModelExtremeparameterαtailmaxnzmaxGEVξ=0．5020，μ=1．3764，σ=0．0539ξ=0．4690，μ=7．1652，σ=0．2542Normalμ=1．4472，σ=0．1368μ=7．4546，σ=0．5007Lognormalμ=0．3656，σ=0．0882μ=2．0067，σ=0．0649Weibulla=1．5147，b=9．1173a=7．7047，b=13．4147Exponentialμ=1．4472μ=7．4546EVμ=1．5247，σ=0．1732μ=7．7276，σ=0．5878

3.2 擬合優(yōu)度檢驗(yàn)

采用K-S檢驗(yàn)法、A-D檢驗(yàn)法和卡方檢驗(yàn)法對極值參數(shù)進(jìn)行擬合優(yōu)度檢驗(yàn)，來判斷飛行極值參數(shù)分布函數(shù)的準(zhǔn)確性。表4和表5為擬合優(yōu)度檢驗(yàn)結(jié)果。在平尾負(fù)迎角和法向過載的極值參數(shù)分布中，GEV分布的K-S檢驗(yàn)值小于0.1，P值大于0.05，表明在95%的置信水平下能通過檢驗(yàn)。GEV分布的A-D檢驗(yàn)值和χ2值相對其他分布非常小，A-D檢驗(yàn)的P值遠(yuǎn)大于0.05，體現(xiàn)出了較高的擬合性。因此，對比GEV和其他分布可以得出，GEV分布的辨識精度遠(yuǎn)遠(yuǎn)高于其他分布。

表4 極值參數(shù)αtailmax的擬合優(yōu)度檢驗(yàn)Table 4 Goodness-of-fit test of extreme parameter αtailmax

表5 極值參數(shù)nzmax的擬合優(yōu)度檢驗(yàn)Table 5 Goodness-of-fit test of extreme parameter nzmax

3.3 極值參數(shù)的概率密度圖和累計(jì)分布圖

圖6和圖7分別為平尾負(fù)迎角和法向過載極值參數(shù)的概率密度圖和累計(jì)分布圖。由圖6和圖7 可知，GEV分布對平尾負(fù)迎角和法向過載極值樣本的描述最為準(zhǔn)確，其他分布類型不能如實(shí)反映極值樣本特征，存在較大誤差。

圖6 平尾負(fù)迎角極值的概率密度圖和累計(jì)分布圖Fig.6 Probability density and cumulative probability maps of extreme negative tail AOA

圖7 法向過載極值的概率密度圖和累計(jì)分布圖Fig.7 Probability density and cumulative probability maps of extreme normal overload

綜上所述，平尾結(jié)冰條件下一維極值參數(shù)(平尾負(fù)迎角和法向過載)的分布符合GEV分布。

3.4 一維極值參數(shù)飛行風(fēng)險(xiǎn)概率計(jì)算

根據(jù)平尾結(jié)冰飛行風(fēng)險(xiǎn)發(fā)生的判定條件，將辨識得到的GEV未知參數(shù)ξ、μ、σ代入GEV分布式(14)，分別以平尾負(fù)迎角和法向過載作為風(fēng)險(xiǎn)評估因子，利用式(20)和式(21)求得風(fēng)險(xiǎn)概率分別為0.011 6和0.021 7，根據(jù)MIL-STD-882D標(biāo)準(zhǔn)可知，事故發(fā)生概率都達(dá)到了B水平，即“可能的”。表明平尾結(jié)冰對飛機(jī)安全性危害極大。

(20)

(21)

相同條件下獨(dú)立求解平尾負(fù)迎角和法向過載的風(fēng)險(xiǎn)概率是不同的，說明僅僅考慮某個(gè)極值參數(shù)評估飛行風(fēng)險(xiǎn)具有一定局限性，因此需要綜合考慮多個(gè)飛行極值參數(shù)。

4 二維極值Copula模型

考慮到極值參數(shù)之間的相關(guān)性問題，就需要研究邊緣分布之間的相關(guān)性，建立相關(guān)性核函數(shù)，本文利用的Copula模型就是采用相同的思路，但廣義的Copula模型無法準(zhǔn)確描述平尾結(jié)冰條件下極值參數(shù)之間的相關(guān)性，因此，選取適合平尾結(jié)冰條件下極值參數(shù)相關(guān)性的Copula函數(shù)是本節(jié)的重點(diǎn)。

設(shè)隨機(jī)向量的(αtailmax,nzmax)的分布函數(shù)為F(αtailmax,nzmax)，邊緣分布函數(shù)分別為服從GEV分布的u=F1(αtailmax)、v=F2(nzmax)。則對于任意的(αtailmax,nzmax)∈R，一定存在一個(gè)CopulaC，使得

C(u,v)=C(F1(αtailmax),F2(nzmax))=

F(αtailmax,nzmax)

(22)

4.1 二維極值參數(shù)辨識

常見的對厚尾特性較敏感的二維阿基米德Copula模型有Gumbel Copula模型(式(23))，Clayton Copula模型(式(24))，GS Copula模型(式(25))，F(xiàn)rank Copula模型(式(26))和Joe Copula模型(式(27))。Gumbel Copula模型、GS Copula模型和Joe Copula模型對飛行參數(shù)在分布上尾處的變化較為敏感，能較好地描述上尾相關(guān)性；Frank Copula模型對變量之間有近似對稱的上下尾部相關(guān)結(jié)構(gòu)有較好的描述；Clayton Copula模型對飛行參數(shù)在分布下尾處的變化十分敏感，能較好地描述下尾相關(guān)性。本文根據(jù)二維極值參數(shù)的分布規(guī)律選取了Joe Copula 模型。

C(u,v)=exp-[(-lnu)β+(-lnv)β]1/θ

(23)

C(u,v)=(u-β+v-β-1)-1/β

(24)

(25)

(26)

C(u,v)=1-[(1-u)β+(1-v)β-

(27)

式中：β為不同模型參數(shù)。

基于符合GEV的飛行極值參數(shù)對以上Copula模型中的未知參數(shù)進(jìn)行辨識，結(jié)果如表6所示。

表6 不同Copula模型的參數(shù)辨識Table 6 Identified parameters of different Copula models

4.2 擬合優(yōu)度檢驗(yàn)

從表7可以看出，5種Copula模型K-S檢驗(yàn)的P值均大于0.05，表明在95%的置信水平下，均能通過檢驗(yàn)。其中，Joe Copula模型的AIC(最小信息準(zhǔn)則)、BIC值(貝葉斯信息規(guī)則，與AIC準(zhǔn)則一樣，數(shù)值越小，擬合精度越高)、χ2值和K-S 值最小，P值最大。因此，相對于其他Copula模型，Joe Copula模型對二維極值分布的描述精度最高。

表7 不同Copula模型的擬合優(yōu)度檢驗(yàn)Table 7 Goodness-of-fit test for different Copula models

4.3 Copula模型的概率密度圖

根據(jù)表6中未知參數(shù)的辨識結(jié)果，畫出5種Copula模型的概率密度圖，如圖8所示。圖中X、Y軸分別是以平尾迎角和法向過載為極值參數(shù)的一維極值分布函數(shù)F1(αtailmax)和F2(nzmax)，Z軸是所求的二維極值飛行風(fēng)險(xiǎn)概率密度。從圖8可以看出，和其他Copula模型相比，Joe Copula模型對厚尾特性描述最好，因此選取Joe Copula模型作為二維極值分布的描述模型。

4.4 二維極值參數(shù)飛行風(fēng)險(xiǎn)概率計(jì)算

依據(jù)Joe Copula模型及平尾結(jié)冰飛行風(fēng)險(xiǎn)發(fā)生的判定條件求出相應(yīng)風(fēng)險(xiǎn)概率為

(28)

圖8 不同Copula模型的概率密度圖Fig.8 Probability density maps of different Copula models

將辨識出的未知參數(shù)代入Joe Copula模型，計(jì)算得到遭遇平尾結(jié)冰為5 min的飛行風(fēng)險(xiǎn)概率為0.026 5(隨著設(shè)定的遭遇結(jié)冰時(shí)間的增長，其值增大)?？梢钥闯觯C合考慮二維極值參數(shù)計(jì)算出的平尾結(jié)冰飛行風(fēng)險(xiǎn)概率值不是簡單地由單個(gè)極值算出的風(fēng)險(xiǎn)概率值相加，它比兩者都大，但比兩者之和小，具有更高的準(zhǔn)確性和可信度。

5 結(jié) 論

1) 確定對平尾結(jié)冰飛行風(fēng)險(xiǎn)影響最大的兩個(gè)參數(shù)為平尾負(fù)迎角和法向過載，采用蒙特卡羅法提取出相應(yīng)飛行極值參數(shù)，驗(yàn)證了所提取參數(shù)極值符合一維廣義極值分布，并建立了平尾結(jié)冰飛行風(fēng)險(xiǎn)判斷條件，計(jì)算出一維極值飛行風(fēng)險(xiǎn)概率。

2) 對多種Copula模型進(jìn)行參數(shù)辨識和擬合優(yōu)度檢驗(yàn)，驗(yàn)證了Joe Copula模型對二維極值分布的描述準(zhǔn)確度最高。依據(jù)Joe Copula模型及平尾結(jié)冰飛行風(fēng)險(xiǎn)判定條件計(jì)算出二維極值平尾結(jié)冰飛行風(fēng)險(xiǎn)概率，其值比一維極值飛行風(fēng)險(xiǎn)概率大，具有更高的可信度與準(zhǔn)確性。

3) 所提出的平尾結(jié)冰飛行風(fēng)險(xiǎn)概率評估方法量化地描述了平尾結(jié)冰對飛行風(fēng)險(xiǎn)的影響。對飛機(jī)平尾結(jié)冰情形下的風(fēng)險(xiǎn)規(guī)避、控制等問題具有一定的參考價(jià)值，是現(xiàn)有飛行適航性準(zhǔn)則及飛行安全規(guī)范在風(fēng)險(xiǎn)評估方面的有效補(bǔ)充。

4) 由于飛行事故的發(fā)生受內(nèi)外部多方面因素影響，是一個(gè)隨機(jī)的不確定過程，不可能全面準(zhǔn)確地考慮到所有因素，因此本文和SAE-ARP-4761、MIL-STD-882E中的飛行風(fēng)險(xiǎn)量化概率在某種程度上都是一個(gè)參考值，和真實(shí)值存在一定誤差。但在平尾結(jié)冰飛行風(fēng)險(xiǎn)預(yù)測、不同結(jié)冰情況下飛行風(fēng)險(xiǎn)的橫向?qū)Ρ确治龊惋w行事故發(fā)生概率的等級劃分等方面具有積極的意義。

[1] RANAUDO R J, BATTERSON J G, REEHORST A L, et al. Determination of longitudinal aerodynamic derivatives using flight data from an icing research aircraft： AIAA-1989-0754[R]. Reston: AIAA,1989.

[2] RATVASKY T P, VAN ZANTE J F. NASA/FAA tailplane icing program: flight test report: NASA/TP-2000-209908[R]. Washington, D.C.: NASA, 2000.

[3] BRAGG M B. Aircraft aerodynamic effects due to large droplet ice accretions: AIAA-1996-0932[R]. Reston: AIAA, 1996.

[4] GINGRAS D R, DICKES E G, RATVASKY T P. Modeling of in-flight icing effects for pilot training: AIAA-2002-4605[R]. Reston: AIAA, 2002.

[5] HAMMOND D, QUERO M, IVEY P, et al. Analysis and experimental aspects of the impact of supercooled water droplets into thin water films: AIAA-2005-0077[R]. Reston: AIAA, 2005.

[6] BRAGG M B, HUTCHISON T, MERRET J, et al. Effect of ice accretion on aircraft flight dynamics: AIAA-2000-0360[R]. Reston: AIAA, 2000.

[7] AMANDA L, JOHN V. Prediction of icing effects on the coupled dynamic response of light airplanes[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(3): 656-673.

[8] ANSELL P J, KERHO M F. Envelope protection for contaminant-induced adverse aerodynamics on a wing using flap hinge moment measurements: AIAA-2013-2654[R]. Reston: AIAA, 2013.

[9] THOMAS P R. Demonstration of an ice contamination effects flight training device: AIAA-2006-0677[R]. Reston: AIAA, 2006.

[10] THOMAS P R, BILLY P B, LEE S. Current methods modeling and simulating icing effects on aircraft performance, stability, control[J]. Journal of Aircraft, 2010, 47(1): 201-211.

[11] 徐忠達(dá), 蘇媛, 曹義華.平尾結(jié)冰對飛機(jī)縱向氣動(dòng)參數(shù)的影響[J].航空學(xué)報(bào), 2013,34(7): 1563-1570.

XU Z D, SU Y, CAO Y H. Effects of tailplane icing on aircraft longitudinal aerodynamic parameters[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(7): 1563-1570 (in Chinese).

[12] 史剛, 李云. Y-8飛機(jī)平尾積冰導(dǎo)致的飛行事故分析[J]. 飛行力學(xué), 2011, 29(5): 84-86.

SHI G, LI Y. Study of flight accidents in a row caused by horizontal tail icing[J]. Flight Dynamics, 2011, 29(5): 84-86 (in Chinese).

[13] 潘環(huán), 艾劍良. 飛機(jī)結(jié)冰冰形預(yù)測的建模與仿真[J]. 系統(tǒng)仿真學(xué)報(bào), 2014, 26(1): 221-224.

PAN H, AI J L. Modeling and simulation of aircraft ice shape prediction[J]. Journal of System Simulation, 2014, 26(1): 221-224 (in Chinese).

[14] 周莉, 徐浩軍, 楊哲. 飛機(jī)在結(jié)冰條件下的最優(yōu)邊界保護(hù)方法[J]. 上海交通大學(xué)學(xué)報(bào), 2013, 47(8): 1217-1221.

ZHOU L, XU H J, YANG Z. Optimal boundary protection method for aircraft under icing conditions[J]. Journal of Shanghai Jiao Tong University, 2013, 47(8): 1217-1221 (in Chinese).

[15] 陳斌, 王立文. 飛機(jī)除冰液地面除冰過程模型仿真與實(shí)驗(yàn)[J]. 系統(tǒng)仿真學(xué)報(bào), 2012, 24(3): 556-560.

CHEN B, WANG L W. Model simulation and experiment of aircraft deicing process using deicing fluids on ground[J]. Journal of System Simulation, 2012, 24(3): 556-560 (in Chinese).

[16] GJB900-90. General program for system safety[S]. 1990.

[17] ASE. Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment: SAE ARP 4761[S]. New York: SAE, 1996.

[18] DOD. Airworthiness certification criteria: MIL-HDBK-516B[S]. Washington, D.C.: DOD, 2005.

[19] DOD. Standard practice for system safety: MIL-STD-882D[S]. Washington, D.C.: DOD, 2000.

[20] STUART C. An introduction to statistical modeling of extreme value[M]. London: Springer, 2007.

[21] LAMPTON A, VALASEK J. Prediction of icing effects on the lateral/directional stability and control of light airplanes[J]. Aerospace Science and Technology, 2012: 305-311.

[22] STELIOS D B, DIMITRIS A G. Estimation of value-at-risk by extreme value and conventional methods: A comparative evaluation of their predictive performance[J]. Journal of International Financial Markets, Institutions & Money, 2006, 8: 209-228.

[23] JOE H. Asymptotic efficiency of the two-stage estimation method for Copula-based models[J]. Journal of Multivariate Analysis, 2005, 94(2): 401-419.

[24] FREY R，MENEIL A J. Copula and credit models[J]. The Risk Metrics, 2001.

[25] NELSEN R B. An introduction to copulas[M]. New York: Springer, 1999.

王健名男， 碩士研究生。主要研究方向： 飛行器飛行品質(zhì)與綜合控制。

E-mail: 1454893168@qq.com

徐浩軍男， 教授， 博士生導(dǎo)師。主要研究方向： 飛行安全與作戰(zhàn)效能。

Tel.: 029-84787637

E-mail: xuhaojun@xjtu.edu.cn

URL：www.cnki.net/kcms/detail/11.1929.V.20160114.0930.002.html

Flightriskevaluationoftailplaneicingbasedonextremevaluetheory

WANGJianming,XUHaojun*,XUEYuan,WANGXiaolong,LIZhe

AeronauticsandAstronauticsEngineeringCollege,AirForceEngineeringUniversity,Xi’an710038,China

AnewmethodcombiningextremevaluetheoryandCopulamodelsisproposedtoquantitativelyevaluatetheflightriskoftailplaneicing.Byestablishingthecomplexpilot-aircraft-environmentmodel,thesituationoftailplaneicingduringapproachingandlandingissimulated.Theflightextremeparameterswhichareprovedtofitthegeneralizedextremevalue(GEV)distributionareextractedthroughMonteCarlomethod.Accordingtothedefinitionofflightriskandrelevantsafetycriterions,theflightriskdeterminationconditionisbuilttocomputetheflightriskprobabilityofone-dimensionalextreme.Thencopulamodelsarechosetodescribethecorrelationoftwo-dimensionalextremeparameters,andunknownparametersindifferentCopulamodelsareidentified.Theresultsofgoodness-of-fittestshowthatJoeCopulamodelhasthehighestaccuracywhendescribingthedistributionoftwo-dimensionalextremeparameters.Thus,theflightriskprobabilityoftwo-dimensionalextremeparametersiscalculatedusingJoeCopula,whichsolvesthelimitationofone-dimensionalextremeparameter.Theapproachhascertainreferencevaluesforthetheoriesofflightsafetyassessment,andprovidesanalysisandteststandardforpreventingflightaccidentinthecircumstanceoftailplaneicing.

extremevaluetheory;Copulamodel;tailplaneicing;flightriskprobability;MonteCarlomethod;parameteridentification

2015-12-04；Revised2015-12-30；Accepted2016-01-07；Publishedonline2016-01-140930

s：NationalNaturalScienceFoundationofChina(61374145,61503406);NationalBasicResearchProgramofChina(2015CB755802)

.Tel.：029-84787637E-mailxuhaojun@xjtu.edu.cn

2015-12-04；退修日期2015-12-30；錄用日期2016-01-07； < class="emphasis_bold">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2016-01-140930

www.cnki.net/kcms/detail/11.1929.V.20160114.0930.002.html

國家自然科學(xué)基金 (61374145,61503406)； 國家“973”計(jì)劃 (2015CB755802)

.Tel.：029-84787637E-mailxuhaojun@xjtu.edu.cn

王健名， 徐浩軍， 薛源， 等． 基于極值理論的平尾結(jié)冰飛行風(fēng)險(xiǎn)評估J． 航空學(xué)報(bào),2016,37(10):3011-3022.WANGJM,XUHJ,XUEY,etal.FlightriskevaluationoftailplaneicingbasedonextremevaluetheoryJ.ActaAeronauticaetAstronauticaSinica,2016,37(10):3011-3022.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0011

V212

A

1000-6893(2016)10-3011-12