結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)的飛行風(fēng)險概率

薛源， 徐浩軍， 胡孟權(quán)

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

結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)的飛行風(fēng)險概率

薛源， 徐浩軍*， 胡孟權(quán)

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

以結(jié)冰條件下的飛行風(fēng)險量化概率為研究對象，基于蒙特卡羅飛行仿真實驗對結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)的耦合特性進(jìn)行了分析, 并獲取了飛行參數(shù)極值樣本。構(gòu)建了飛行風(fēng)險發(fā)生的判定條件；對飛行參數(shù)極值樣本進(jìn)行了統(tǒng)計特性分析，驗證了其厚尾分布特征。一維分布類型辨識結(jié)果表明廣義極值分布對相對速度和迎角極值的描述精度最高。為描述二維變量對相關(guān)性的各自影響程度，提出了一種新的雙參數(shù)變權(quán)重Copula模型；辨識結(jié)果表明該Copula模型能以較高的精度通過假設(shè)檢驗。相關(guān)性分析的結(jié)果表明相對速度和迎角同時出現(xiàn)極大值和極小值的概率較大?；诙S極值樣本的Copula分布模型求出了不同結(jié)冰程度下的飛行風(fēng)險概率值，探討了飛行風(fēng)險的非線性增長趨勢。

蒙特卡羅仿真； 結(jié)冰條件； 人-機(jī)-環(huán)系統(tǒng)； 飛行風(fēng)險； 極值分布； Copula模型

飛機(jī)結(jié)冰較復(fù)雜的隨機(jī)性和不確定性容易引起駕駛員的耦合連鎖反應(yīng)，導(dǎo)致人-機(jī)-環(huán)系統(tǒng)失去穩(wěn)定性，從而引發(fā)飛行風(fēng)險或?qū)е嘛w行事故。如何評估結(jié)冰條件下的飛行風(fēng)險一直是一個難題。自從航空器出現(xiàn)以來，由于飛機(jī)內(nèi)部部件失效、駕駛員操縱水平等內(nèi)在原因?qū)е碌娘w行風(fēng)險及飛行事故均可以被有效地解決以及避免，但結(jié)冰導(dǎo)致的飛行風(fēng)險一直伴隨著航空器的發(fā)展史，每年均有類似的飛行事故發(fā)生。根據(jù)ICAO(International Civil Aviation Organization)和美國NTSB(National Transportation Safety Board)的飛行事故數(shù)據(jù)庫統(tǒng)計，從2002—2014年，在所有氣象因素帶來的飛行事故中，12%是由于飛機(jī)結(jié)冰導(dǎo)致的，并且其中92%是在飛行中發(fā)生積冰。近年來，最嚴(yán)重的一起結(jié)冰事故發(fā)生在2009年6月 1日，法國航空447號班機(jī)的空客A330客機(jī)在大西洋上空墜毀，造成216名乘客以及12名機(jī)組人員死亡，事故的原因為結(jié)冰導(dǎo)致自動駕駛儀關(guān)閉，飛行員隨后錯誤操作導(dǎo)致失速。

結(jié)冰條件下的飛行仿真與動力學(xué)計算一直是研究外部環(huán)境風(fēng)險的重要組成部分[1]。如：Bragg等提出了物理意義明顯、結(jié)構(gòu)簡單的描述積冰影響的模型[2]，可用于估算積冰后飛機(jī)的氣動參數(shù)；Hui等通過對Convair 580飛機(jī)的建模研究指出，積冰氣象條件的微觀處理參數(shù)與積冰后氣動導(dǎo)數(shù)的變化具有特定的相關(guān)性，可以使用無量綱的參數(shù)去體現(xiàn)積冰后配平以及穩(wěn)定性和操縱性導(dǎo)數(shù)的變化[3]；Krzysztof等模擬了飛機(jī)在積冰環(huán)境的爬升過程，再現(xiàn)了該飛機(jī)失事的過程[4]；Frank等對三維機(jī)翼積冰后失速迎角與最大升力系數(shù)降低量的關(guān)系進(jìn)行了分析[5]；Thomas等利用“雙水獺”飛機(jī)積冰試飛數(shù)據(jù)，應(yīng)用多種飛行策略對多項飛行導(dǎo)數(shù)、操縱導(dǎo)數(shù)進(jìn)行了研究，并著重分析了機(jī)翼和平尾的失速問題[6]；王明豐等在突風(fēng)擾動情況下，分別對無積冰、弱積冰和嚴(yán)重積冰情況下飛機(jī)飛行參數(shù)的時間響應(yīng)進(jìn)行了計算[7]；Lampton和Valasek利用簡化的積冰影響模型，對左右機(jī)翼非對稱積冰進(jìn)行了仿真，同時還進(jìn)行了縱向、橫向和全機(jī)積冰條件的飛機(jī)仿真分析，通過比較積冰前后飛行參數(shù)的變化分析積冰對飛機(jī)穩(wěn)定性和操縱性的影響[8-10]；美國伊利諾伊州大學(xué)的Bragg等提出并領(lǐng)導(dǎo)了飛機(jī)智能積冰系統(tǒng)研究[11]；Robert和Glen開發(fā)了結(jié)冰遭遇飛行仿真模擬器[12]；David等對結(jié)冰后的飛行包線保護(hù)系統(tǒng)進(jìn)行了研究[13-15]。

上述工作具有較高的學(xué)術(shù)價值，但在結(jié)冰條件下飛行風(fēng)險發(fā)生概率的評估與預(yù)測方面尚缺乏相關(guān)的理論與方法，無法量化評估結(jié)冰條件下動態(tài)的飛行風(fēng)險概率。而結(jié)冰條件下的飛行風(fēng)險評估也正是當(dāng)前飛機(jī)安全性和適航性工作重點關(guān)注的內(nèi)容，是完善和改進(jìn)現(xiàn)有的飛行安全分析理論亟待解決的問題?，F(xiàn)有的SAE ARP 4761[16]、SAE ARP 4754A[17]、MIL-HDBK-516B[18]、MIL-STD-882E[19]和MIL-STD-1797B[20]等安全性規(guī)范與指南對于評估由于飛機(jī)內(nèi)部硬件系統(tǒng)失效導(dǎo)致飛行事故概率有明確的思路及方法。國內(nèi)外學(xué)術(shù)界對飛行安全與飛行風(fēng)險評估的理論與方法也有較多的研究工作[21-25]。但其大多是定性的分析或者是靜態(tài)的可靠性評估，此類方法在評估結(jié)冰條件下的飛行風(fēng)險概率方面具有明顯的局限性：一是由于此類方法基于確定性的模型，其不能反應(yīng)結(jié)冰條件下多因素耦合復(fù)雜情形的隨機(jī)性與不確定性；二是不能基于航空器的飛行過程評估動態(tài)的飛行風(fēng)險大小，對于結(jié)冰導(dǎo)致的飛行風(fēng)險較難進(jìn)行量化概率指標(biāo)的評估預(yù)測。

另一方面，Copula理論成為多元極值問題的研究熱點[26]，在金融、災(zāi)害預(yù)警、信號處理等方面均有應(yīng)用[27-36]。由于評估飛行風(fēng)險需研究多元飛行參數(shù)極值(飛行參數(shù)極值)的尾部相關(guān)性，故較適合將多元極值Copula模型應(yīng)用到其中。但尚未見到將這一有效理論應(yīng)用到結(jié)冰飛行風(fēng)險評估中的公開報道。鑒于此，本文首先對直觀反映結(jié)冰后運(yùn)動狀態(tài)的飛行參數(shù)進(jìn)行分析。由于結(jié)冰后飛行參數(shù)的變化可以直觀反映人-機(jī)-環(huán)系統(tǒng)的運(yùn)動狀態(tài)，飛行參數(shù)極值的超限情形即是飛行風(fēng)險發(fā)生的判定條件；因此，在對一維極值分布進(jìn)行研究的基礎(chǔ)上，對結(jié)冰條件下基于多元極值Copula模型的飛行風(fēng)險概率量化評估方法進(jìn)行探索；擬在有效解決一元極值模型局限性的基礎(chǔ)上利用多元極值分布求得結(jié)冰空間中的飛行風(fēng)險概率。

1 基于蒙特卡羅仿真實驗提取結(jié)冰條件下的飛行參數(shù)極值

一方面，嚴(yán)重結(jié)冰條件下的試飛具有較大的飛行風(fēng)險；另一方面，研究結(jié)冰條件下的隨機(jī)性與不確定性需要的數(shù)據(jù)量較大。因此文中基于蒙特卡羅法考慮不同程度結(jié)冰條件下的物理特性，對結(jié)冰后的動態(tài)過程進(jìn)行多次計算迭代，從而提取反映人-機(jī)-環(huán)系統(tǒng)在結(jié)冰后物理特性與隨機(jī)不確定性的飛行參數(shù)極值樣本點。為下文中基于多元極值Copula的飛行風(fēng)險評估方法研究提供數(shù)據(jù)輸入。飛行參數(shù)極值的提取系統(tǒng)基于某型飛機(jī)的地面試驗系統(tǒng)改造而成，蒙特卡羅仿真實驗的流程如圖1所示。

首先，設(shè)定結(jié)冰程度，提取特定結(jié)冰程度下的結(jié)冰氣動模型；然后，利用蒙特卡羅法將結(jié)冰條件下的內(nèi)外部環(huán)境變量按照其出現(xiàn)頻率進(jìn)行隨機(jī)抽樣，從而對每次計算迭代過程中所使用的內(nèi)外部環(huán)境參數(shù)產(chǎn)生影響；同時，飛行員對每一次蒙特卡羅仿真實驗的迭代計算過程進(jìn)行操縱，從而提供結(jié)冰后的舵面輸入信號，對相關(guān)的氣動參數(shù)及操縱信號產(chǎn)生量化影響。如此便可以通過多次的循環(huán)迭代反映人-機(jī)-環(huán)系統(tǒng)在結(jié)冰條件下的隨機(jī)性與不確定性。飛機(jī)本體方程為基于四元數(shù)法的六自由度方程，微分算法為四階龍格庫塔(Runge-Kutta)算法，仿真實驗的時間步長為20 ms。氣動數(shù)據(jù)與基于某型機(jī)的全機(jī)測力風(fēng)洞試驗，所得到的各氣動參數(shù)精度均滿足合格或先進(jìn)指標(biāo)的要求。

結(jié)冰對飛機(jī)氣動參數(shù)的影響模型是飛行力學(xué)建模的一個至關(guān)重要的環(huán)節(jié)。在飛行仿真計算中常用Bragg等提出的積冰影響模型：

C(A)ice=(1+η(t)KCA)C(A)

(1)

圖1 基于蒙特卡羅仿真的飛行參數(shù)極值提取過程Fig.1 Extreme flight parameters extraction process based on Monte Carlo simulation

2 結(jié)冰條件下飛行參數(shù)分析與風(fēng)險判定

蒙特卡羅仿真實驗的初始條件為：高度5 000 m，速度120 m/s，某型飛機(jī)以配平狀態(tài)進(jìn)入中度結(jié)冰區(qū)域，η(t)的變化范圍為0～0.35，開始結(jié)冰時全機(jī)質(zhì)量為49 000 kg，蒙特卡羅仿真的次數(shù)為150次。選取兩次較典型的飛行員分別在正確與錯誤操縱狀態(tài)下飛機(jī)結(jié)冰后的響應(yīng)為案例進(jìn)行研究。一次為第53次，飛行員在結(jié)冰后介入延遲較大且操縱不當(dāng)情形下飛機(jī)的響應(yīng)；一次為第91次，飛行員及時介入且操縱恰當(dāng)情形下飛機(jī)的響應(yīng)。圖2為第53次和第91次蒙特卡羅仿真實驗中飛行參數(shù)的變化情況。從圖2中可看到在第53次迭代計算中，當(dāng)飛機(jī)進(jìn)入結(jié)冰空域后，結(jié)冰情形逐漸加劇，飛行員在55.9 s時感受到過載和俯仰角的變化，開始操縱駕駛桿進(jìn)行修正以期重新回到配平狀態(tài)，但由于飛行員介入的延遲較大，加之修正時的操縱動作過猛，迎角與速度并未達(dá)到修正預(yù)期，同時引發(fā)了橫向滾轉(zhuǎn)效應(yīng)，偏航角開始改變，飛行員在拉桿的同時還需橫向壓桿以平衡姿態(tài)。在111.2 s時，飛行員注意到速度有加速減小的趨勢，開始操縱油門桿增加推力，138.7 s時，迎角超出極限值，飛機(jī)進(jìn)入失速振蕩，速度繼續(xù)減小，飛行員此時將駕駛桿拉到極限并保持油門最大位置，但仍然無法改出失速，此刻即標(biāo)志著飛行風(fēng)險的發(fā)生。在第91次迭代計算中，飛機(jī)在50 s時結(jié)冰效應(yīng)開始顯現(xiàn)，飛行員在53.7 s時及時介入，平穩(wěn)操縱駕駛桿穩(wěn)定姿態(tài)，并提前推油門以保持速度和高度，迎角和速度在短暫變化后趨于平穩(wěn)，飛機(jī)未進(jìn)入風(fēng)險狀態(tài)。第53次和第91次迭代計算所提取的迎角和速度極值[αmax,vmax]分別為[21.332 3°,63.247 8 m/s]和[6.335 9°,116.788 3 m/s]，如圖2所示。

150次蒙特卡羅仿真實驗后所提取的150組極值樣本散點分布如圖3所示，其中相對速度為飛機(jī)配平速度與結(jié)冰速度的比值，將這個比值的極大值記為vstable/vmin，文中vstable=120 m/s。如此可以將對速度極小值的研究轉(zhuǎn)換為對相對速度極大值的研究。圖3中標(biāo)注為上文中第53次和第91次蒙特卡羅仿真實驗的極值樣本。

圖2 第53次和第91次蒙特卡羅仿真實驗中飛行參數(shù)變化情況Fig.2 Flight parameter changes in the 53rd and 91st Monte Carlo simulation experiments

圖3 150組極值樣本散點分布Fig.3 Scatter distributions of 150 extreme samples

根據(jù)150次蒙特卡羅迭代計算的結(jié)果可確定迎角或者速度超限為飛行風(fēng)險發(fā)生的判定條件。對文中所涉及到的飛行風(fēng)險進(jìn)行定義如下：以超過95%的概率極易引起MIL-STD-882E[19]中所定義的風(fēng)險范疇中評估值為1-5的災(zāi)難性飛行事故。即不能安全飛行和著陸的失效情況，引起飛機(jī)結(jié)構(gòu)損傷并導(dǎo)致至少一人的傷亡。給出判斷文中定義的尾流飛行風(fēng)險是否發(fā)生的判定條件為

(2)

3 一維飛行參數(shù)極值樣本的統(tǒng)計特性

從圖3中可以看出滾轉(zhuǎn)角極值、下降高度極值和迎角極值的分布都存在明顯的厚尾特性，這種分布形式在低頻高危事件(如地震、海嘯、金融風(fēng)險、飛行事故等)中較為常見，文中涉及遭遇尾流情形下的飛行風(fēng)險亦屬于低頻高危事件的范疇。針對此種分布形式，目前較有效的描述方法為采用極值理論。下文將針對飛行參數(shù)極值樣本分布模型的不確定性，對尾流遭遇情形下一維飛行參數(shù)極值的尾部分布規(guī)律進(jìn)行研究，采用最優(yōu)化辨識算法確定分布模型中的未知參數(shù)，并利用辨識后的分布模型對一維飛行參數(shù)極值進(jìn)行擬合優(yōu)度檢驗，從而找出最合適的分布類型。在此之前，首先對飛行參數(shù)極值變量的統(tǒng)計特性進(jìn)行研究。對飛行參數(shù)極值變量進(jìn)行統(tǒng)計分析，結(jié)果如表1所示。

表1 極值樣本統(tǒng)計量Table 1 Statistics of extreme samples

觀察表1可以發(fā)現(xiàn)，相對速度極值與迎角極值的最大值均比最小值偏離均值與中位數(shù)的程度要大，說明其分布形式并不是左右對稱的。繼續(xù)分析表1可以看到三組極值樣本滾轉(zhuǎn)角、下降高度、迎角的峰度系數(shù)均大于3，說明三組極值樣本均有比正態(tài)分布更長的尾部；其偏斜度均大于0°，表明分布類型在右側(cè)具有較長尾部。

圖4為三組極值樣本的盒形圖，盒形中加粗的線為樣本的中位數(shù)，上下邊界分別為樣本正態(tài)分布范圍的25%與75%界限?？梢钥闯鰞山M極值參數(shù)在上尾均有多個樣本點超出正態(tài)分布的界限范圍。顯然，亦可得到樣本點的分布具有上厚尾特征的結(jié)論。

圖4 極值樣本盒形圖Fig.4 Box graph of extreme samples

4 一維飛行參數(shù)極值樣本的分布類型辨識

4.1 分布類型假設(shè)

針對圖3中的極值樣本，初步推斷廣義極值分布(式(3))較適合描述其厚尾分布特性。為驗證廣義極值分布的有效性，采用廣義極值分布模型和另外5種可以描述厚尾特性的分布模型來辨識一維極值參數(shù)，從而尋找對一維極值參數(shù)樣本分布的描述最準(zhǔn)確的分布模型。所采用的分布模型主要有：廣義極值分布(GEV)、正態(tài)分布模型(Normal)、對數(shù)正態(tài)分布(Lognormal)、威布分布(Weibull)、指數(shù)分布(EXP)、極值分布(EV)。其中廣義極值分布模型為

(3)

式中：ξ∈R，μ∈R，σ>0，1+ξ(x-μ)/σ>0。μ為位置參數(shù)；σ為尺度參數(shù)；ξ為形狀參數(shù)。

4.2 辨識結(jié)果

利用自適應(yīng)粒子群優(yōu)化(Adaptive Range Particle Swarm Optimization，ARPSO)算法辨識出一維極值分布模型的未知參數(shù)，如表2所示。在表2中，GEV的待辨識依次為ξ、μ和σ，Normal的待辨識參數(shù)為μ和σ，Lognormal的待辨識參數(shù)為μ和σ，Weibull的待辨識參數(shù)為a和b，EXP的待辨識參數(shù)為μ，EV的待辨識參數(shù)為μ和σ。

4.3 GEV描述一維極值分布的準(zhǔn)確性驗證

表2 極值參數(shù)的辨識結(jié)果Table 2 Identification results of extreme parameters

表3 相對速度極值的擬合優(yōu)度檢驗Table 3 Fitting test for extreme relative velocity

表4 迎角極值的擬合優(yōu)度檢驗Table 4 Fitting test for extreme angle of attack

從圖5和圖6中亦可看出GEV對極值樣本的概率密度和累積概率的擬合最為準(zhǔn)確，其他分布形式不能充分反映極值樣本的分布信息。

圖5 相對速度極值的概率密度圖與累積概率圖Fig.5 Probability density graph and cumulative probability graph of extreme relative velocity

圖6 迎角極值的概率密度圖與累積概率圖Fig.6 Probability density graph and cumulative probability graph of extreme angle of attack

(4)

(5)

根據(jù)式(2)的飛行風(fēng)險發(fā)生條件，將表2中的GEV分布參數(shù)ξ、μ和σ代入到式(6)的F2(αmax)中，再根據(jù)式(7)求得以迎角αmax為風(fēng)險評估參數(shù)時發(fā)生飛行風(fēng)險的概率為0.118 0。

(6)

Pr=1-F2(αmax/αc(δf,Ma,η)>1)

(7)

由上文可知，式(5)和式(7)求出的飛行風(fēng)險概率值是不相同的，這說明了由于飛行風(fēng)險的判定需考慮多個具有決定性作用的極值參數(shù)，利用一維極值參數(shù)與一元極值GEV分布去求飛行風(fēng)險概率不能全面地考慮其他飛行參數(shù)極值對飛行風(fēng)險的影響，具有一定的局限性。而多元極值理論可以考慮到所有對飛行風(fēng)險發(fā)生起決定作用的飛行參數(shù)極值，這也正是下面所探討的內(nèi)容。

5 結(jié)冰條件下的Copula模型構(gòu)建

上文對一維極值分布進(jìn)行了辨識，并分別給出了相對速度和迎角作為單變量臨界極值時的飛行風(fēng)險概率。但根據(jù)式(2)，文中結(jié)冰后飛行風(fēng)險的判定牽扯到二維極值參數(shù)的相關(guān)性結(jié)構(gòu)。有關(guān)單變量的結(jié)論并不能平行推廣到多元情形，其主要原因是多維隨機(jī)變量的聯(lián)合分布除了與各分量的邊緣分布有關(guān)之外，更重要的是與變量之間相關(guān)性的關(guān)系。當(dāng)隨機(jī)變量的個數(shù)比較大時，單個分量的極值行為未必含有整個向量的聯(lián)合極值行為。因此，下文基于Copula理論構(gòu)建二維極值的相關(guān)結(jié)構(gòu)，由于Copula理論是針對極值分布的評估而提出的，故在對極值相關(guān)性的描述上，Copula極值分布模型能較好地反映極值參數(shù)之間的聯(lián)系和發(fā)展趨勢。

5.1 二維分布模型假設(shè)

(8)

(9)

式中：Copula模型的未知參數(shù)ω≥1。當(dāng)ω=1時，表示極值參數(shù)u與v相互獨(dú)立；ω=+∞表示u與v完全相關(guān)。

(10)

式中：ω≠0;如果ω>0，表示u與v正相關(guān);ω→0，表示u與v相互獨(dú)立;ω<0，表示u與v負(fù)相關(guān)。

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

(11)

式中：ω∈(0,∞)，當(dāng)ω→0時，u與v趨于獨(dú)立；ω=+∞表示u與v完全相關(guān)。

(12)

式中：ω≥1，當(dāng)ω=1時，極值參數(shù)u與v相互獨(dú)立；當(dāng)ω=+∞時，極值參數(shù)u與v完全相關(guān)。

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

(1-u)ω(1-v)ω]1/ω

(13)

式中：ω≥1，當(dāng)ω=1時，極值參數(shù)u與v相互獨(dú)立；當(dāng)ω=+∞時，極值參數(shù)u與v完全相關(guān)。

5.2 提出二維極值變量的DAWP Copula模型

根據(jù)上文可以初步判定對上尾變化敏感的Gumbel Copula模型及Joe Copula模型最能反映文中極值的分布情況，但其單參數(shù)的特征使得在描述二維變量對相關(guān)性的各自影響程度時具有一定的局限性，故本文在Gumbel模型的基礎(chǔ)上提出一種新的雙參數(shù)變權(quán)重Copula模型(Double Adaptive Weight Parameter Copula,DAWP Copula)，如式(14)所示。

(14)

5.3 二維Copula模型的未知參數(shù)辨識

(15)

(16)

2) 求出Copula的密度函數(shù)為

c(ui,vi;ω1,ω2，…,ωd)=

(17)

3) 對于式(17)中未知參數(shù)ω1,ω2,…,ωd的估計，結(jié)合似然函數(shù)與χ2分布檢驗構(gòu)建目標(biāo)函數(shù)：

H(u1,u2,…,un;v1,v2,…,vn;ω1,ω2,…,ωd)=

(18)

式中：Ai為極值樣本點在區(qū)間(u≤ui,v≤vi)內(nèi)的個數(shù)。

4) 利用非線性辨識算法求出未知參數(shù)的值，如表5所示。

表5Copula模型中未知參數(shù)辨識結(jié)果

Table5IdentificationresultsofunknownparametersinCopulamodels

CopulamodelUnknownparameterGumbelω=2.3353Frankω=6.6261Claytonω=1.3951GSω=1.5058Joeω=3.0506DAWPω1=-0.9159,ω2=-0.6137

5.4 二維Copula模型的擬合優(yōu)度檢驗

對于涉及到的Copula模型，分別應(yīng)用AIC(Akaike Information Criteria)準(zhǔn)則、BIC(Bayesian Information Criteria)準(zhǔn)則、χ2檢驗、K-S檢驗評價其擬合優(yōu)度，如表6所示。

表6Copula模型對極值樣本的擬合優(yōu)度檢驗

Table6FittingtestofextremesamplefordifferentCopulamodels

CopulamodelAICBICχ2K-SPvalueGumbel-149.6779-146.66730.22330.09690.2653Frank-118.0216-115.01100.47240.21720.0764Clayton-86.7373-83.72660.36930.36470.0081GS-122.7043-119.69360.13010.20130.0235Joe-143.6251-140.61441.43800.11560.2693DAWP-149.2077-146.19710.18680.08120.3139

從表6中可以看出，Gumbel Copula、Joe Copula、DAWP Copula的P值均大于顯著性水平0.01，0.02，0.05，即這3種Copula在99%、98%、95%的置信水平下均能通過檢驗；而Clayton Copula甚至在99%的置信水平下亦未能通過檢驗。繼續(xù)比較發(fā)現(xiàn)Gumbel Copula和DAWP Copula的AIC與BIC檢驗值比其他Copula的值均要小，說明二者的辨識精度較高。而DAWP Copula的K-S值為0.081 2，比Gumbel Copula的K-S值還要小，同時其P值亦是所有 Copula 模型中最大，說明其對極值參數(shù)相關(guān)結(jié)構(gòu)的描述更為準(zhǔn)確。

5.5 相關(guān)性分析

5.6 基于DAWP Copula模型求飛行風(fēng)險概率

根據(jù)風(fēng)險判定式(2)及構(gòu)建的DAWP Copula模型求出風(fēng)險概率為

F2(αmax/αc(δf,Ma,η)>1))

(19)

CDAWP代表基于式(14)的DAWP Copula模型。根據(jù)表5中辨識出的DAWP Copula未知參數(shù)求得在第1節(jié)提到的仿真特定點上的飛行風(fēng)險概率Pr=0.130 4。對比通過一維飛行參數(shù)極值GEV模型(式(5)和式(7))求出的飛行風(fēng)險概率，可發(fā)現(xiàn)基于多元極值Copula模型求出的飛行風(fēng)險概率值比式(5)和式(7)求出的風(fēng)險概率值要大，說明其充分考慮了不同飛行參數(shù)極值超限時的飛行風(fēng)險狀況，相對于一維飛行參數(shù)極值判定飛行風(fēng)險更全面、更精確。

6 不同結(jié)冰程度下的飛行風(fēng)險

不同冰型、不同結(jié)冰位置、不同結(jié)冰厚度等會對氣動參數(shù)產(chǎn)生不同的量化影響，但其均會使得氣動參數(shù)惡化，使飛機(jī)的飛行品質(zhì)降低，操控性變差。針對文中所研究的背景飛機(jī)，在相同的高度速度以及外部環(huán)境參數(shù)下，按照不同結(jié)冰情形對氣動參數(shù)的影響程度，將無結(jié)冰狀態(tài)到極嚴(yán)重結(jié)冰狀態(tài)均分為10個等級。無結(jié)冰等級為0，此時對氣動參數(shù)無影響；極嚴(yán)重結(jié)冰等級為9，此時結(jié)冰導(dǎo)致氣動參數(shù)極具惡化，對于式(1)來說，等級9意味著參數(shù)η將大于0.3，阻力系數(shù)增加300%以上，升力系數(shù)減小30%以上。采用上文方法基于蒙特卡羅仿真實驗提取每個狀態(tài)下的迎角與速度極值樣本，利用DAWP Copula模型對樣本進(jìn)行辨識后計算出每個狀態(tài)發(fā)生飛行風(fēng)險的概率如表7和圖7所示。表7中的等級7對應(yīng)于5.6節(jié)案例中所求的飛行風(fēng)險0.130 4。

表7不同結(jié)冰等級下的飛行風(fēng)險概率

Table7Flightriskprobabilitiesunderdifferenticinglevels

IcinglevelFlightriskprobability00.00001910.00020620.00091530.00138140.00617250.01559360.05682970.13041180.25337590.360572

圖7 不同結(jié)冰等級下的飛行風(fēng)險概率曲線Fig.7 Flight risk probability curve under different icing levels

由表7和圖7可以看出，無結(jié)冰狀態(tài)發(fā)生飛行風(fēng)險的概率為0.000 019，基于MIL-STD-882E中飛行事故率的劃分可認(rèn)為此時沒有飛行風(fēng)險。從狀態(tài)0到狀態(tài)5，飛行風(fēng)險概率平穩(wěn)增大，狀態(tài)5 的飛行風(fēng)險為0.015 593，已經(jīng)高于SAE ARP-4761所規(guī)定系統(tǒng)安全性水平中的B級危險性失效狀態(tài)概率，在MIL-STD-882E從1～20的風(fēng)險評級中定為2級。從狀態(tài)5到狀態(tài)9，也就是中度結(jié)冰到極嚴(yán)重結(jié)冰這一區(qū)間，飛行風(fēng)險概率呈指數(shù)形態(tài)激增。狀態(tài)8嚴(yán)重結(jié)冰時的風(fēng)險概率已經(jīng)高達(dá)25.337 5%，此概率值遠(yuǎn)遠(yuǎn)高于MIL-STD-882E中所規(guī)定A級(經(jīng)常發(fā)生)事件10%的出現(xiàn)概率，在MIL-STD-882E從1～20的風(fēng)險評級中定為1級。這一階段需要駕駛員提前做出反應(yīng)并正確平穩(wěn)操縱飛機(jī)才能避免飛行風(fēng)險的發(fā)生。從圖7中也可以看出，從無結(jié)冰到極嚴(yán)重結(jié)冰，飛行風(fēng)險概率的增長不是一個線性的過程，其在輕度結(jié)冰時的數(shù)值較小，在中度結(jié)冰狀態(tài)之后開始激增，為一種新的飛行風(fēng)險變化趨勢。

7 結(jié) 論

1) 考慮結(jié)冰條件下的不確定性與隨機(jī)性因素，利用蒙特卡羅仿真實驗研究了人-機(jī)-環(huán)系統(tǒng)在結(jié)冰條件下的動力學(xué)特性。結(jié)果表明在中度結(jié)冰條件下某型飛機(jī)的動力學(xué)特性逐步惡化，迎角和法向過載增量較大，易于超出臨界值以致飛機(jī)失控；同時由于阻力系數(shù)的增加和升力系數(shù)的減小，飛行速度和飛行高度逐漸降低。如駕駛員在此時介入延遲較大且操縱不當(dāng)，則極易引起飛行參數(shù)極值超限從而引發(fā)飛行風(fēng)險；駕駛員需在中度結(jié)冰條件下迅速介入并平穩(wěn)操縱以脫離飛行風(fēng)險。

2) 驗證了所提取極值樣本具有厚尾分布的特性。擬合優(yōu)度檢驗的結(jié)果表明GEV模型能以較高的精度對飛行參數(shù)極值樣本的分布特性進(jìn)行描述，從而驗證了迎角極值和相對速度極值均符合GEV分布。結(jié)合多元極值Copula模型提出了一種量化評估飛行風(fēng)險概率的新方法，并將其用到結(jié)冰條件下的風(fēng)險概率量化評估中。為了解決單參數(shù)Copula模型在描述多維變量對相關(guān)性的各自影響程度時具有的局限性，提出了二維的DAWP Copula模型?；贒AWP Copula模型對某型飛機(jī)不同結(jié)冰等級下的飛行風(fēng)險進(jìn)行了計算。結(jié)果表明從無結(jié)冰到極嚴(yán)重結(jié)冰狀態(tài)，飛行風(fēng)險概率的增長不是一個線性的過程，為一種新的變化趨勢。

本文思路可對其他內(nèi)外部環(huán)境因素影響下的飛行風(fēng)險概率評估提供參考。其中的飛行風(fēng)險概率評估方法是對現(xiàn)有各類飛行安全規(guī)范中風(fēng)險評估理論的有效補(bǔ)充，對于飛行安全與適航性管理具有積極的作用。需注意的是，由于飛行事故的發(fā)生是一個多因素影響的不確定過程，不可能將所有內(nèi)外部隨機(jī)因素考慮完全。因此，本文得到的飛行風(fēng)險量化概率在多數(shù)狀況下是一種類似于飛行安全指南(如SAE ARP-4761與MIL-STD-882E)中事故率的參考值，與真實值必然有一定的誤差。但其在不同狀況下飛行風(fēng)險的橫向?qū)Ρ确治?、風(fēng)險程度的歸類劃分中具有積極的意義。同時，由于本文側(cè)重于風(fēng)險評估方法，主要考慮了結(jié)冰對氣動參數(shù)的影響，沒有全面考慮如發(fā)動機(jī)結(jié)冰、傳感器結(jié)冰以及舵面結(jié)冰卡滯等其他結(jié)冰風(fēng)險因素。因此，下一步工作將利用文中方法綜合考慮飛機(jī)結(jié)冰后多種不利因素相互耦合交聯(lián)的情形。

[1] CONNOR P, DEA A, KENNEDY Q, et al. Measuring safety climate in aviation: A review and recommendations for the future[J]. Safety Science, 2011, 49(2): 128-138.

[2] BRAGG M B, PERKINS W R, SARTER N B, et al. An interdisciplinary approach to inflight aircraft icing safety: AIAA-1998-0095[R]. Reston: AIAA, 1998.

[3] HUI K, WOLDE M, BROWN A. Flight dynamics model of turboprop transport aircraft icing effects based on preliminary flight data: AIAA-2005-1068[R]. Reston: AIAA, 2005.

[4] KRZYSZTOF S, MACIEJ L, EDYTA L, et al. Aircraft flight dynamics with simulated ice accretion: AIAA-2004-4948[R]. Reston: AIAA, 2004.

[5] FRANK T L, ADBOLLAH K. Effects of ice accretions on aircraft aerodynamics[J]. Progress in Aerospace Sciences, 2001, 37(8): 669-767.

[6] THOMAS P R, BILLY P B, SAM L. Current methods for modeling and simulating icing effects on aircraft performance stability and control: AIAA-2008-6204[R]. Reston: AIAA, 2008.

[7] 王明豐, 王立新, 黃成濤. 積冰對飛機(jī)縱向操穩(wěn)特性的量化影響[J]. 北京航空航天大學(xué)學(xué)報, 2008, 34(5): 592-595.

WANG M F, WANG L X, HUANG C T. Computational effects of ice accretion on aircraft longitudinal stability and control[J]. Journal of Beijng University of Aeronautics and Astronautics, 2008, 34(5): 592-595 (in Chinese).

[8] LAMPTON A, VALASEK J. Prediction of icing effects on the dynamic response of light airplanes[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(3): 722-732.

[9] LAMPTON A, VALASEK J. Prediction of icing effects on the coupled dynamic response of light airplanes[J]. Journal of Guidance, Control, and Dynamics, 2008, 31(3): 656-673.

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

[11] BRAGG M B, BASAR T, PERKINS W R, et al. Smart icing systems for aircraft icing safety: AIAA-2002-0813[R]. Reston: AIAA, 2002.

[12] ROBERT W D, GLEN A D. Icing encounter flight simulator[J]. Journal of Aircraft, 2006, 43(5): 1528-1537.

[13] DAVID R G, BILLY B, RICHARD R, et al. Development and implementation of a model-driven envelope protection system for in-flight ice contamination: AIAA-2010-8141[R]. Reston: AIAA, 2010.

[14] RICHARD R, BORJA M, BILLY N, et al. Piloted simulation to evaluate the utility of a real time envelope protection system for mitigating in-flight icing hazards: AIAA-2010-7987[R]. Reston: AIAA, 2010.

[15] DAVID R G. Requirements and modeling of in-flight icing effects for flight training: AIAA-2013-5075[R]. Reston: AIAA, 2013.

[16] Society of Automotive Engineers. Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment： ARP 4761[S]. Washington, D.C.: SAE,1996.

[17] Society of Automotive Engineers. Certification considerations for high-integrated or complex aircraft systems： ARP 4754[S]. Washington, D.C.: SAE, 2010.

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

[19] USA Department of Defense. Standard practice for system safety： MIL-STD-882E[S]. Washington, D.C.: DOD, 2012.

[20] USA Department of Defense. Flying qualities of piloted aircraft： MIL-STD-1797B[S]. Washington, D.C.: DOD, 2012:673-695.

[21] BROOKER P. Experts, Bayesian belief networks, rare events and aviation risk estimates[J]. Safety Science, 2011, 49(8): 1142-1155.

[22] WANG W H, JIANG X B, XIA S C. Incident tree model and incident tree analysis method for quantified risk assessment: An in-depth accident study in traffic operation[J]. Safety Science, 2010, 48(10): 1248-1262.

[23] MATTHEWS B, DAS S, BHADURI K, et al. Discovering anomalous aviation safety events using scalable data mining algorithms[J]. Journal of Aerospace Information Systems, 2013, 10(10): 467-475.

[24] OCAMPO J, MILLWATER H, SINGH G, et al. Development of a probabilistic linear damage methodology for small aircraft[J]. Journal of Aircraft, 2011， 48(6): 2090-2106.

[25] BALACHANDRAN S, ATKINS E M. A constrained Markova decision process framework for flight safety assessment and management: AIAA-2015-0115[R]. Reston: AIAA, 2015.

[26] NELSEN R B. An introduction to copulas[M]. 2rd ed. New York: Springer, 2006: 51-108.

[27] DIKS C, PANCHENKO V, SOKOLINSKIY O, et al. Comparing the accuracy of multivariate density forecasts in selected regions of the copula support[J]. Journal of Economic Dynamics & Control, 2014, 48: 79-94.

[28] SUKCHAROEN K, ZOHRABYAN T, LEATHAM D, et al. Interdependence of oil prices and stock market indices: A copula approach[J]. Energy Economics, 2014, 44: 331-339.

[30] YASMIN S, ELURU N, ABDUL R, et al. Examining driver injury severity in two vehicle crashes—A copula based approach[J]. Accident Analysis and Prevention, 2014, 66(3): 120-135.

[31] MASIN M, LAMBERTI A, ARCHETTI R. Coastal flooding: A copula based approach for estimating the joint probability of water levels and waves[J]. Coastal Engineering, 2015, 97: 37-52.

[32] BESSA R J, MIRANDA V, BOTTERUD A, et al. Time-adaptive quantile-copula for wind power probabilistic forecasting[J]. Renewable Energy, 2012, 40(1): 29-39.

[33] ERYILMAZ S. Estimation in coherent reliability systems through Copulas[J]. Reliability Engineering and System Safety, 2011, 96(5): 564-568.

[34] BERGER T. Forecasting value-at-risk using time varying copulas and EVT return distributions[J]. International Economics, 2013, 133: 93-106.

[35] MOAZAMI S, GOLIAN S, KAVIANPOUR M R, et al. Uncertainty analysis of bias from satellite rainfall estimates using copula method[J]. Atmospheric Research, 2014, 137(2): 145-166.

[36] XUE Y, XU H J, WANG X L. Build probability distribution maps of flight risk during wake encountering[J]. Journal of Aircraft, 2015, 52(3): 805-818.

薛源男， 博士， 講師。主要研究方向： 飛行仿真與飛行安全。

Tel.: 029-84787637

E-mail: xuhaojun@xjtu.edu.cn

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

Tel.: 029-84787637

E-mail: xuhaojun@xjtu.edu.cn

*Correspondingauthor.Tel.：029-84787637E-mail:xuhaojun@xjtu.edu.cn

Flightriskprobabilityofpilot-aircraft-environmentsystemundericingconditions

XUEYuan,XUHaojun*,HUMengquan

AeronauticsandAstronauticsEngineeringCollege,AirForceEngineeringUniversity,Xi’an710038,China

Thequantitativeflightriskprobabilityundericingconditionsissetastheresearchobject.BasedonMonteCarloflightsimulationexperiment,couplingcharacteristicsofpilot-aircraft-environmentsystemaftericingarestudied,andtheextremeflightparametersareextracted.Thejudgementconditionsforflightriskaregiven.Andstatisticalpropertiesoftheextremesamplesareanalyzed.Thefactthatextremeflightparametershaveheavytaildistributioncharacteristicsisverified.Identificationresultsindicatethatgeneralizedextremevaluemodelcanbestdescribethedistributioncharacteristicsofrelativevelocityandangleofattack.Inordertodescribetheinfluencelevelofeachtwo-dimensionalparameteronthecorrelations,anewCopulamodelthathastwochangeableweightsisproposed;identificationresultshowsthatthisnewCopulamodelcanpassthehypothesistestingwithhighaccuracy.Correlationanalysisresultsrevealthatrelativevelocityandangleofattackwouldbothshowupextremevalueswithgreaterprobability.FlightriskprobabilitiesunderdifferenticinglevelsarecalculatedbasedonCopuladistributionmodelwithtwo-dimensionalextremevalues.Thenonlineargrowthtrendofflightriskisalsodiscussed.

MonteCarlosimulation;icingconditions;pilot-aircraft-environmentsystem;flightrisk;extremevaluedistribution;Copulamodel

2016-01-03；Revised2016-02-22；Accepted2016-03-16；Publishedonline2016-05-041416

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

s：NationalBasicResearchProgramofChina(2015CB755802)；NationalNaturalScienceFoundationofChina(61503406)

2016-01-03；退修日期2016-02-22；錄用日期2016-03-16； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-05-041416

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

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

*

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

薛源， 徐浩軍， 胡孟權(quán)． 結(jié)冰條件下人-機(jī)-環(huán)系統(tǒng)的飛行風(fēng)險概率J． 航空學(xué)報,2016,37(11):3328-3339.XUEY,XUHJ,HUMQ.Flightriskprobabilityofpilot-aircraft-environmentsystemundericingconditionsJ.ActaAeronauticaetAstronauticaSinica,2016,37(11):3328-3339.

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

10.7527/S1000-6893.2016.0086

V212; X949

A

1000-6893(2016)11-3328-12