基于多維極值參數(shù)的飛行風(fēng)險(xiǎn)量化評(píng)估方法

薛 源，徐浩軍，裴彬彬，陳怡然

（1．空軍工程大學(xué)航空航天工程學(xué)院，陜西西安710038；2．空軍北京航空裝備訓(xùn)練基地，北京100076）

基于多維極值參數(shù)的飛行風(fēng)險(xiǎn)量化評(píng)估方法

薛 源1，徐浩軍1，裴彬彬1，陳怡然2

（1．空軍工程大學(xué)航空航天工程學(xué)院，陜西西安710038；2．空軍北京航空裝備訓(xùn)練基地，北京100076）

結(jié)合極值理論與多因素耦合系統(tǒng)建模仿真思路，提出了一種基于Copula的多維極值風(fēng)險(xiǎn)評(píng)估方法。針對(duì)飛行過(guò)程中的復(fù)雜隨機(jī)性，基于蒙特卡羅法提取所需要的三維極值參數(shù)，驗(yàn)證了所提取極值參數(shù)具有和試驗(yàn)數(shù)據(jù)相同的分布形式，并構(gòu)建了飛行風(fēng)險(xiǎn)發(fā)生的判定條件。在對(duì)一維極值參數(shù)符合廣義極值分布的假設(shè)進(jìn)行證明的基礎(chǔ)上，提出了三維極值參數(shù)的四參數(shù)變權(quán)重（four adaptive weight parameters，F(xiàn)AWP），Copula模型利用自適應(yīng)粒子群算法對(duì)一維和三維目標(biāo)函數(shù)中的未知參數(shù)進(jìn)行了辨識(shí)，對(duì)多種Copula辨識(shí)出的三維極值分布進(jìn)行了擬合優(yōu)度檢驗(yàn)，結(jié)果表明FAWP Copula對(duì)三維極值參數(shù)分布形式的描述最為精確。利用FAWP Copula模型對(duì)尾流遭遇情形下的飛行風(fēng)險(xiǎn)概率進(jìn)行了量化計(jì)算，所得指標(biāo)可用來(lái)研究尾流場(chǎng)內(nèi)的風(fēng)險(xiǎn)規(guī)避策略及算法。

多維極值參數(shù)；廣義極值分布；Copula模型；自適應(yīng)粒子群算法；飛行風(fēng)險(xiǎn)概率

0 引 言

飛行風(fēng)險(xiǎn)的量化概率是一個(gè)重要的參考指標(biāo)，對(duì)于飛機(jī)的適航性與飛行安全具有重要的意義。對(duì)由于硬件故障導(dǎo)致的飛行風(fēng)險(xiǎn)概率在ARP-4761［1］，ARP-4754［2］，MILHDBK-516B__［3］，STD-882D［4］及國(guó)軍標(biāo)GJB900-90［5］等國(guó)內(nèi)外安全性規(guī)范中有明確的量化指標(biāo)和計(jì)算方法，但大多是根據(jù)系統(tǒng)硬件故障實(shí)驗(yàn)樣本進(jìn)行可靠性指標(biāo)的計(jì)算與分配，其考慮飛行過(guò)程中的不確定性和隨機(jī)性因素較少，對(duì)于如何求出外部環(huán)境導(dǎo)致的飛行事故及飛行風(fēng)險(xiǎn)概率指標(biāo)尚缺乏相關(guān)的理論與方法。而飛行風(fēng)險(xiǎn)一般是由多個(gè)因素耦合引起的。據(jù)統(tǒng)計(jì)，92%的飛行事故是由多個(gè)因子導(dǎo)致，平均每個(gè)事故有4．39個(gè)誘發(fā)因素，多的可達(dá)20個(gè)［6］。許多歐美科學(xué)家［6-9］經(jīng)大量研究認(rèn)為：事故的發(fā)生通常是由偶然的、耦合作用的不安全因素累積導(dǎo)致的。在上述這種多因素耦合、強(qiáng)非線性的情況下，現(xiàn)有國(guó)內(nèi)外安全性規(guī)范［15］中的飛行風(fēng)險(xiǎn)評(píng)估理論與方法較難進(jìn)行飛行風(fēng)險(xiǎn)的量化預(yù)測(cè)，因此如何考慮多因素耦合情形提取風(fēng)險(xiǎn)評(píng)估樣本，并利用有限的樣本數(shù)據(jù)對(duì)飛行風(fēng)險(xiǎn)進(jìn)行量化評(píng)估是飛行安全研究的難點(diǎn)問(wèn)題。

利用樣本參數(shù)評(píng)估飛行風(fēng)險(xiǎn)概率首先需建立描述此樣本分布的理論模型。由于飛行風(fēng)險(xiǎn)屬于低頻高危事件（如：地震、海嘯、金融風(fēng)險(xiǎn)、飛行事故）［10-13］的范疇，所提取的極值樣本一般具有厚尾特性，針對(duì)此種分布形式的描述目前較有效的方法為采用極值理論［14］。但多因素耦合復(fù)雜飛行情形下的風(fēng)險(xiǎn)概率的評(píng)估牽扯到多維極值參數(shù)及其相關(guān)性特征，一維極值理論的計(jì)算結(jié)論并不能平行推廣到多維情形，因此需探索對(duì)多維極值參數(shù)空間進(jìn)行描述的方法。目前對(duì)二維及以上的參數(shù)空間進(jìn)行評(píng)估較常用的有效方法是構(gòu)造參數(shù)間的相關(guān)性結(jié)構(gòu)，如近些年比較流行的支持向量機(jī)對(duì)多維空間的分類(lèi)，其實(shí)質(zhì)就是構(gòu)造相關(guān)性核函數(shù)。文中使用Copula理論［15-19］描述多維極值間的相關(guān)結(jié)構(gòu)，由于Copula理論是專為評(píng)估極值分布而提出的，故在對(duì)極值相關(guān)性的描述上，Copula極值模型能較好地反映極值參數(shù)之間的聯(lián)系和發(fā)展趨勢(shì)，具有較高的精度，較適合應(yīng)用于本文當(dāng)中。

鑒于此，本文擬在考慮多因素耦合情形下的不確定性和隨機(jī)性的基礎(chǔ)上，提取多維極值參數(shù)。探索基于多元Copula的飛行風(fēng)險(xiǎn)概率量化評(píng)估與預(yù)測(cè)方法，擬解決一維極值模型的局限性，并利用Copula對(duì)尾流遭遇情形下三維極值參數(shù)的相關(guān)性進(jìn)行分析。研究成果擬在飛行風(fēng)險(xiǎn)定量評(píng)估方法上有所創(chuàng)新，對(duì)現(xiàn)有的安全性規(guī)范進(jìn)行補(bǔ)充，同時(shí)為飛行事故的預(yù)測(cè)與預(yù)防提供分析和檢驗(yàn)的新方法。

1 基于蒙特卡羅法的三維極值參數(shù)提取過(guò)程

研究?jī)?nèi)外部因素的復(fù)雜隨機(jī)性所需要的數(shù)據(jù)量較大，試飛數(shù)據(jù)與人在回路地面實(shí)驗(yàn)數(shù)據(jù)無(wú)法滿足數(shù)據(jù)量的要求；加之模擬飛行風(fēng)險(xiǎn)具有較大的安全隱患且需要考慮條件較多，導(dǎo)致實(shí)驗(yàn)條件較為苛刻；因此文中基于蒙特卡羅法對(duì)飛行風(fēng)險(xiǎn)科目進(jìn)行多次仿真實(shí)驗(yàn)以提取評(píng)估飛行風(fēng)險(xiǎn)所需要的飛行參數(shù)極值。極值參數(shù)提取實(shí)驗(yàn)系統(tǒng)基于某型飛機(jī)鐵鳥(niǎo)地面試驗(yàn)系統(tǒng)改造而成，實(shí)物如圖1所示。某型飛機(jī)鐵鳥(niǎo)臺(tái)是飛機(jī)的地面試驗(yàn)臺(tái)，軟硬件結(jié)構(gòu)均與真實(shí)飛機(jī)相同，使用經(jīng)過(guò)風(fēng)洞實(shí)驗(yàn)與試飛驗(yàn)證后的氣動(dòng)參數(shù)，在各個(gè)科目下與真實(shí)試飛數(shù)據(jù)的誤差不超過(guò)12%。

文中評(píng)估的飛行風(fēng)險(xiǎn)案例引自GJB 625A-2001［20］中固定翼飛機(jī)復(fù)雜科目里的第22項(xiàng)：著陸進(jìn)場(chǎng)滾轉(zhuǎn)。設(shè)置風(fēng)險(xiǎn)條件為：飛機(jī)在近地面的情況下進(jìn)入前方飛機(jī)所產(chǎn)生的尾渦中，導(dǎo)致了飛行狀態(tài)的突變，從而引發(fā)飛行風(fēng)險(xiǎn)的概率。在每次計(jì)算飛參數(shù)據(jù)之前，首先利用蒙特卡羅方法將內(nèi)外部影響條件變量按照其出現(xiàn)頻率進(jìn)行隨機(jī)抽樣，需抽樣的條件變量如表1所示。將抽樣的變量數(shù)值作為檢索條件從數(shù)據(jù)庫(kù)中提取尾流數(shù)據(jù)、飛行員操縱行為特性參數(shù)以及其他影響飛行狀況的條件數(shù)據(jù)，從而對(duì)每次計(jì)算迭代過(guò)程中所使用的參數(shù)產(chǎn)生影響，以此反映真實(shí)內(nèi)外部環(huán)境影響下的隨機(jī)性與不確定性。在對(duì)全部隨機(jī)參數(shù)進(jìn)行蒙特卡羅抽樣后，將其動(dòng)態(tài)代入到實(shí)驗(yàn)系統(tǒng)的飛行仿真計(jì)算迭代中，以對(duì)相關(guān)的氣動(dòng)參數(shù)及操縱信號(hào)產(chǎn)生具有多因素耦合效應(yīng)的量化影響。飛機(jī)本體方程為基于四元數(shù)法的六自由度方程，微分算法為四階龍格庫(kù)塔算法。使用（real time workshop，RTW）將Simulink搭建的控制系統(tǒng)轉(zhuǎn)化為實(shí)時(shí)系統(tǒng)Vx Works支持的C代碼，將其下載到實(shí)時(shí)仿真機(jī)，時(shí)間步長(zhǎng)為20 ms。

圖1 提取極值參數(shù)實(shí)驗(yàn)系統(tǒng)

表1 需用蒙特卡羅法抽樣的影響條件變量___________

2 三維極值參數(shù)的可信度驗(yàn)證

圖2顯示了在第61次的蒙特卡羅抽樣仿真過(guò)程中，出現(xiàn)的飛參極值超限從而導(dǎo)致發(fā)生下文中定義的飛行風(fēng)險(xiǎn)?？梢钥闯鲎畛踹M(jìn)入尾流場(chǎng)時(shí)由于尾流中上升氣流的作用，高度和迎角均略有增加，隨后由于滾轉(zhuǎn)力矩的作用，飛機(jī)急劇滾轉(zhuǎn)并高度下降。確定對(duì)尾流飛行風(fēng)險(xiǎn)發(fā)生影響最大的3個(gè)飛行參數(shù)（滾轉(zhuǎn)角φ、下降高度ΔH、迎角α）。本次迭代的極值參數(shù)φimax＝75.61°，ΔHimax＝87.84 m，αimax＝7.91°。

圖2 第61次迭代中的飛參圖

因本文選取的近地近距尾流遭遇情形屬于高風(fēng)險(xiǎn)科目，故不可能采用試飛驗(yàn)證，因此采用飛行員在回路的飛機(jī)地面鐵鳥(niǎo)臺(tái)實(shí)驗(yàn)數(shù)據(jù)作為驗(yàn)?zāi)?shù)據(jù)。列出在i≤75時(shí)提取的前75個(gè)極值參數(shù)如圖3所示，圖3中紅色標(biāo)示即為i＝61時(shí)的極值參數(shù)。圖4為使用上文方法提取的前75次極值參數(shù)以及在相同條件下進(jìn)行的75次飛行員在回路鐵鳥(niǎo)臺(tái)實(shí)驗(yàn)數(shù)據(jù)的分位數(shù)圖（quantiles-quantiles，QQ）圖，可發(fā)現(xiàn)3種極值參數(shù)φmax、ΔHmax、αmax的QQ圖均近似為直線，說(shuō)明本文提取的極值參數(shù)和實(shí)驗(yàn)數(shù)據(jù)屬于同一種分布類(lèi)型。Kolmogorov-Smirnov（K-S）檢驗(yàn)的結(jié)果亦表明3種極值參數(shù)的K-S值均小于0.1，而P值均大于0.25（即在比95%的置信水平低得多的情況下亦能通過(guò)檢驗(yàn)），故可認(rèn)為本文方法得到的數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)具有相同的分布類(lèi)型，文中方法所提取極值參數(shù)的可信度較高，可以接受其作為評(píng)估飛行風(fēng)險(xiǎn)概率的樣本。

圖3 前75次的三維極值參數(shù)圖

圖4 極值參數(shù)的QQ圖

3 定義尾流飛行風(fēng)險(xiǎn)發(fā)生的條件

在對(duì)分布模型進(jìn)行辨識(shí)之前，首先對(duì)文中所涉及到的飛行風(fēng)險(xiǎn)進(jìn)行定義如下：以超過(guò)95%的概率極易引起STD-882D［4］中所定義的風(fēng)險(xiǎn)范疇中評(píng)估值為1～5的災(zāi)難性飛行事故。即不能安全飛行和著陸的失效情況，引起飛機(jī)結(jié)構(gòu)損傷并導(dǎo)致至少一人的傷亡。

對(duì)3個(gè)極值參數(shù)進(jìn)行歸一化處理。查某型飛機(jī)氣動(dòng)數(shù)據(jù)，極值參數(shù)迎角的臨界值與Ma有關(guān)，例如在襟翼0°度時(shí)，當(dāng)Ma＝0．2時(shí)，臨界迎角αc為20．50°；而當(dāng)Ma＝0．7，αc僅為10.90°。根據(jù)氣動(dòng)數(shù)據(jù)和提取極值參數(shù)時(shí)的Ma進(jìn)行差值處理，得到歸一化的極值迎角參數(shù)為αmax／αc（δf，Ma）。根據(jù)氣動(dòng)手冊(cè)，滾轉(zhuǎn)角的臨界風(fēng)險(xiǎn)極值為φc＝85°，歸一化的極值滾轉(zhuǎn)角參數(shù)為φmax／8。重心下降高度的極值參數(shù)為ΔHmax，以機(jī)翼翼尖剛好觸地時(shí)的狀態(tài)作為風(fēng)險(xiǎn)發(fā)生臨界點(diǎn)，極值參數(shù)ΔHmax的歸一化公式為

式中，b為機(jī)翼展長(zhǎng)，取值為38 m；φ和θ為極值參數(shù)ΔHmax對(duì)應(yīng)時(shí)間點(diǎn)上的滾轉(zhuǎn)角和俯仰角。給出判定文中定義的尾流飛行風(fēng)險(xiǎn)是否發(fā)生的公式為

4 驗(yàn)證一維極值參數(shù)符合廣義極值分布

從圖3和圖4中可以看出極值參數(shù)的分布存在較明顯的厚尾特性，又因?yàn)槲擦黠w行風(fēng)險(xiǎn)屬于低頻高危事件的范疇，故適合采用極值理論對(duì)此種分布形式進(jìn)行描述。極值理論是關(guān)于隨機(jī)變量序列最值漸近分布的理論，利用極值理論能夠有效地對(duì)隨機(jī)序列最值概率分布的尾部進(jìn)行建模，用于描述極值樣本數(shù)據(jù)序列分布的尾部特征。一維極值分布有固定的解析模型可以參照，但多維極值的聯(lián)合分布除了與各分量的分布有關(guān)之外，更重要的是與變量之間的相關(guān)性有關(guān)。當(dāng)極值參數(shù)的個(gè)數(shù)比較大時(shí)，單個(gè)分量的極值行為未必含有整個(gè)向量的聯(lián)合極值行為。因此，本文采用適用于多維極值分布問(wèn)題的Copula理論研究三維極值參數(shù)。

根據(jù)Copula的相關(guān)定理［15］：如果函數(shù)F是多元極值分布函數(shù)，則函數(shù)F的一維邊緣分布必然屬于廣義極值（generalized extreme value，GEV）分布族。因此在對(duì)三維極值參數(shù)的分布形式與相關(guān)結(jié)構(gòu)進(jìn)行研究之前，須驗(yàn)證式（2）中涉及到的一維極值參數(shù)φmax、ΔHmax、αmax符合極值理論中的GEV分布形式（如式（3））。

4.1 辨識(shí)一維極值參數(shù)的分布形式

為驗(yàn)證式（2）涉及到的一維極值參數(shù)φmax、ΔHmax、αmax是否符合極值理論中的GEV分布形式（如式（3）），首先采用不同的分布模型來(lái)描述一維極值參數(shù)，較具有代表性的有廣義Pareto（generalized Pareto，GP）分布、GEV分布、正態(tài)（normal，norm）分布、威布（Weibull distribution，Weibull）分布、指數(shù)（Exponential，Exp）分布、泊松（Poisson）分布。

式中，ξ∈R；μ∈R；σ＞0；1＋ξ（x－μ）／σ＞0。

設(shè)分布族的統(tǒng)一形式為F（x；θ1，θ2，…，θm），f（x；θ1，θ2，…，θm）為其密度函數(shù)，其中θ為未知參數(shù)。將一維極值參數(shù)的子樣值升序排列得到（x1，x2，…，xn），構(gòu)建目標(biāo)函數(shù)

下一步需對(duì)目標(biāo)函數(shù)進(jìn)行辨識(shí)，在對(duì)最小二乘法、極大

（2）計(jì)算每個(gè)粒子對(duì)應(yīng)的目標(biāo)函數(shù)值。

式中，r1和r2是0～1之間的隨機(jī)值；c1和c2是正的常數(shù)，c1＋c2≤4，一般情況下取c1＝c2＝2。

（5）由下式更新權(quán)重w：似然法、模式搜索算法、遺傳算法（genetic algorithm，GA）等辨識(shí)方式進(jìn)行對(duì)比測(cè)試的基礎(chǔ)上，結(jié)合目標(biāo)函數(shù)（如式（4））構(gòu)造復(fù)雜及計(jì)算量較大的特點(diǎn)。選用較適合本文的優(yōu)化環(huán)境的局部搜索能力強(qiáng)且收斂速度較快的粒子群算法［21］對(duì)目標(biāo)函數(shù)進(jìn)行辨識(shí)。但對(duì)于粒子群算法來(lái)說(shuō)，它最初階段給定的搜索范圍通常在以后的整個(gè)搜索迭代過(guò)程中是固定的。隨著迭代過(guò)程的進(jìn)行，最初的搜索區(qū)間變得相對(duì)過(guò)大從而影響了找到最優(yōu)解的速度和精度。而粒子群算法的改進(jìn)算法自適應(yīng)區(qū)間粒子群（adaptive range particle swarm optimization，ARPSO）可以動(dòng)態(tài)地改變搜索區(qū)間，它的具體思路是根據(jù)迭代過(guò)程中變量的變化來(lái)減小動(dòng)態(tài)搜索區(qū)間，從而使粒子的聚集區(qū)間越來(lái)越小進(jìn)而獲得高精度的全局最優(yōu)值。

利用ARPSO對(duì)目標(biāo)函數(shù)的辨識(shí)流程如下：

（1）將目標(biāo)函數(shù)（如式（4））的未知參數(shù)θ＝（θ1，θ2，…，θm）視為一個(gè)m維的粒子，設(shè)置最初的搜索區(qū)間、粒子的數(shù)量和最大搜索迭代次數(shù)kmax。初始化迭代次數(shù)k為1。在此后每次迭代的搜索區(qū)間中隨機(jī)地給定每個(gè)粒子最初的位置和速度。

式中，f表示粒子當(dāng)前迭代的目標(biāo)函數(shù)值；favg和fmin分別表示當(dāng)前迭代中所有粒子的目標(biāo)函數(shù)平均值和最小目標(biāo)函數(shù)值。

（6）更新迭代次數(shù)k為k＝k＋1。計(jì)算變量θik＋1的平均值μki＋1和標(biāo)準(zhǔn)差σik＋1。設(shè)定標(biāo)準(zhǔn)差為＝σ＝（i＝1，2，…，m）。

（7）修正標(biāo)準(zhǔn)差

（8）隨著迭代搜索的進(jìn)行許多隨機(jī)地分布在搜索區(qū)間內(nèi)的粒子會(huì)往全局最優(yōu)點(diǎn)方向移動(dòng)，因此動(dòng)態(tài)搜索區(qū)間會(huì)隨著迭代搜索過(guò)程的進(jìn)行而減小。利用下式來(lái)更新系統(tǒng)參數(shù)a：

式中，amin可以設(shè)置為趨近于0的非負(fù)值（如amin＝1.0× 10－5）；amax可由下式求得：

采用式（10）后，可以在迭代開(kāi)始階段設(shè)置比較大的動(dòng)態(tài)搜索區(qū)間，而隨著迭代的進(jìn)行可以減小動(dòng)態(tài)搜索區(qū)間的范圍，有效地提高了收斂的速度和最優(yōu)解的精度。

（9）根據(jù)下式設(shè)置動(dòng)態(tài)搜索區(qū)間：

1，L和＋1，R分別代表變量θ＋1左右兩邊的標(biāo)準(zhǔn)差。

（10）如果全局最優(yōu)點(diǎn)Pg不在動(dòng)態(tài)搜索區(qū)間內(nèi)，根據(jù)下式調(diào)整動(dòng)態(tài)搜索區(qū)間：

（11）如果超出邊界約束，根據(jù)下式計(jì)算得出標(biāo)準(zhǔn)差，從而設(shè)置更新的動(dòng)態(tài)搜索區(qū)間：

表2 極值參數(shù)φmax的分布模型辨識(shí)結(jié)果

表3 極值參數(shù)ΔHmax的分布模型辨識(shí)結(jié)果

表4 極值參數(shù)αmax的分布模型辨識(shí)結(jié)果

4.2 GEV描述一維極值參數(shù)的準(zhǔn)確性驗(yàn)證

在已知分布函數(shù)的未知參數(shù)后，根據(jù)原采樣極值參數(shù)對(duì)現(xiàn)有分布函數(shù)進(jìn)行擬合優(yōu)度檢驗(yàn)并作出概率密度圖（見(jiàn)圖5）。使用K-S方法的擬合優(yōu)度檢驗(yàn)結(jié)果如表5所示。

圖5 極值參數(shù)的QQ圖與概率密度圖

______表5 極值參數(shù)的擬合優(yōu)度檢驗(yàn)________

觀察圖5可以看出，3個(gè)一維極值參數(shù)GEV分布的QQ圖接近線性，其他分布的QQ圖都不同程度地偏離線性趨勢(shì)，其中Norm分布和GP分布偏離線性較大，從極值參數(shù)的概率密度圖中亦可看出GEV對(duì)極值樣本的描述最為準(zhǔn)確。分析表5可以看出，GEV的K-S值小于其他模型，而P值亦遠(yuǎn)遠(yuǎn)大于其他分布模型，分析P值可知GEV在比95%的置信水平低得多的情況下亦能通過(guò)檢驗(yàn)，而其他的分布模型甚至在99%的置信水平下都未能通過(guò)檢驗(yàn)。

綜上所述，GEV對(duì)極值參數(shù)的描述是極為準(zhǔn)確的，尾流遭遇情形下一維極值參數(shù)的分布符合GEV分布。

5 利用Copula求取飛行風(fēng)險(xiǎn)概率

5.1 提出FAWP Copula

在驗(yàn)證了一維邊緣極值分布符合GEV分布的基礎(chǔ)上，繼續(xù)研究三維極值參數(shù)的分布結(jié)構(gòu)與相關(guān)性。設(shè)極值隨機(jī)向量（φmax，ΔHmax，αmax）的分布函數(shù)為F（φmax，ΔHmax，αmax），邊緣分布函數(shù)分別為符合GEV分布的F1（φmax）、F2（ΔHmax）和F3（αmax），根據(jù)Copula的相關(guān)定理［15］，則對(duì)于任意的（φmax，ΔHmax，αmax）∈Rd，一定存在一個(gè)Copula C，使得

文中的F1（φmax）、F2（ΔHmax）和F3（αmax）都是連續(xù)分布函數(shù)，故C是唯一的。由式（17）定義的函數(shù)C是一個(gè)邊緣分布為GEV形式的三維聯(lián)合分布函數(shù)。對(duì)于本文中三維極值參數(shù)的Copula模型選擇，首先根據(jù)常用的二維Copula來(lái)構(gòu)建，其通用形式如式（18）所示。利用式（18）構(gòu)建的Copula模型主要有：Gumbel Copula（見(jiàn)式（19））、Frank Copula（見(jiàn)式（20））、Clayton Copula（見(jiàn)式（21））、GS Copula（見(jiàn)式（22））、Joe Copula（見(jiàn)式（23））。

根據(jù)三維極值的分布可以初步判定對(duì)上尾變化敏感的Copula模型能較好地反映本文中極值的分布情況。上文中的Gumbel Copula和Joe Copula均對(duì)上尾的變化較敏感，但其未知參數(shù)僅有2個(gè)，這使得在描述三維變量對(duì)相關(guān)性的各自影響程度時(shí)具有一定的局限性，故本文在Gumbel模型的基礎(chǔ)上提出一種四參數(shù)變權(quán)重（four adaptive weight parameters-FAWP）Copula模型，如式（24）所示。

5.2 Copula模型辨識(shí)與描述精度分析

根據(jù)提取的極值參數(shù)對(duì)上文中Copula模型的未知參數(shù)進(jìn)行辨識(shí)，具體步驟如下。

步驟1 根據(jù)定義，Copula的邊緣分布函數(shù)即為一維極值的GEV分布函數(shù)，故將每組三維極值樣本點(diǎn)x、Δ分別代入上文中已辨識(shí)出未知參數(shù)（未知參數(shù)見(jiàn)表2～表4）的式（3）中，極值參數(shù)最終的邊緣累積概率為

記Ui；ξ，μ，σ）為變量ui，Vi（ΔHimax；ξ，μ，σ）為變量vi，Wi（；ξ，μ，σ）為變量wi。

步驟2 根據(jù)下式求出Copula的密度函數(shù)：

步驟3 根據(jù)變量ui、vi、wi的數(shù)值構(gòu)建目標(biāo)函數(shù)

H（u1，u2，…，un；v1，v2，…，vn；w1，w2，…，wn；θ1，θ2，θ3，θ4）＝

式中，Ai表示極值樣本點(diǎn)在區(qū)間（u≤ui，v≤vi，w≤wi）內(nèi)的個(gè)數(shù)。

步驟4 利用ARPSO算法辨識(shí)出目標(biāo)函數(shù)（見(jiàn)式（29））的未知參數(shù)。

根據(jù)表6中6種Copula的辨識(shí)參數(shù)構(gòu)建其三維的概率密度如圖6所示，以此更直觀地觀察Copula模型的密度函數(shù)特征。從圖中亦可以看出Gumbel Copula，Joe Copula，F(xiàn)AWP Copula對(duì)三維極值參數(shù)的厚尾特性描述較好，對(duì)比w＝0．85時(shí)極值參數(shù)（ui，vi）的邊緣分布等高線圖，亦可以看出FAWP Copula對(duì)三維極值參數(shù)分布的描述是合適的。對(duì)于本文中涉及到的Copula模型，分別應(yīng)用AIC（Akaike information criteria）準(zhǔn)則、BIC（Bayesian information criteria）準(zhǔn)則、χ2檢驗(yàn)法、K-S檢驗(yàn)法評(píng)價(jià)其描述極值分布的準(zhǔn)確程度，結(jié)果如表7所示。

表6 Copula模型的參數(shù)辨識(shí)表

圖6 w＝0．85時(shí)Copula模型對(duì)參數(shù)（ui，vi）的描述

表7 Copula模型對(duì)極值參數(shù)的描述優(yōu)度檢驗(yàn)

從表7中可以看出，F(xiàn)rank Copula，Gumbel Copula，Joe Copula，F(xiàn)AWP Copula的P值大于顯著性水平0．01，0．02，0．05，即這4種Copula在99%，98%，95%的置信水平下均能通過(guò)檢驗(yàn)；而Clayton Copula甚至在99%的置信水平下亦未能通過(guò)檢驗(yàn)；Gumbel Copula，Joe Copula，F(xiàn)AWP Copula在遠(yuǎn)小于95%的置信水平下亦能通過(guò)檢驗(yàn)，3種Copula的AIC值、BIC值亦比較小，故完全可以接受這3種Copula作為三維極值分布的描述模型；但同時(shí)FAWP Copula的P值最大，χ2和K-S值最小，說(shuō)明其對(duì)極值參數(shù)相關(guān)結(jié)構(gòu)的描述更為準(zhǔn)確。

5.3 極值參數(shù)相關(guān)性分析

利用FAWP Copula的二維形式構(gòu)建極值參數(shù)見(jiàn)的相關(guān)函數(shù)如下：

式中，A（ω）是凸的；A（0）＝A（1）＝1；max｛ω，1－ω｝≤A（ω）≤1（ω∈［0，1］）；如果A（ω）＝1，則2個(gè)隨機(jī)變量相互獨(dú)立；如果A（ω）＝max｛ω，1－ω｝，則完全相關(guān)。

在變量ω從0～1變化的過(guò)程中，極值參數(shù)（φmax，ΔHmax）的FAWP Copula相關(guān)函數(shù)A（ω）的值在0.85～1之間變動(dòng)，但始終不等于1，計(jì)算得知其相關(guān)程度在0.3左右；極值參數(shù)（φmax，αmax）的FAWP Copula相關(guān)函數(shù)A（ω）的值在0.93～1之間變動(dòng)，其相關(guān)程度較弱，計(jì)算得知其相關(guān)程度在0.15左右；極值參數(shù)（ΔHmax，αmax）的FAWP Copula相關(guān)函數(shù)A（ω）的值在0.96～1之間變動(dòng)，其相關(guān)程度亦較弱，計(jì)算得知其相關(guān)程度在0.1左右。

5.4 評(píng)估飛行風(fēng)險(xiǎn)概率

最終，根據(jù)構(gòu)建的三維極值參數(shù)的FAWP Copula模型求出風(fēng)險(xiǎn)概率，如

根據(jù)表6辨識(shí)出的未知參數(shù)構(gòu)建FAWP Copula的解析模型，求得在第一節(jié)提到的仿真特定點(diǎn)上的飛行風(fēng)險(xiǎn)概率Pr為0.0512。需注意的是，由于飛行事故的發(fā)生是一個(gè)多因素影響的不確定過(guò)程，不可能將所有內(nèi)外部隨機(jī)因素考慮完全，因此文獻(xiàn)［1］與文獻(xiàn)［4］中的事故率很大程度上是一個(gè)參考值。文中得到的飛行風(fēng)險(xiǎn)量化概率值在多數(shù)狀況下亦是一個(gè)參考值，與真實(shí)值必然有一定的誤差。但其在不同狀況下飛行風(fēng)險(xiǎn)的橫向?qū)Ρ确治?、風(fēng)險(xiǎn)程度的歸類(lèi)劃分中具有積極的意義。如：不同惡劣環(huán)境條件下或不同硬件故障條件下的風(fēng)險(xiǎn)大小對(duì)比，預(yù)定科目或任務(wù)的風(fēng)險(xiǎn)程度預(yù)測(cè)比較等。

6 結(jié) 論

（1）結(jié)合人-機(jī)-環(huán)復(fù)雜系統(tǒng)建模方法與多元極值理論提出了一種量化評(píng)估飛行風(fēng)險(xiǎn)概率的新思路及新方法，并將其用到尾流場(chǎng)的風(fēng)險(xiǎn)概率量化評(píng)估中。首先考慮尾流遭遇情形下的不確定性與隨機(jī)性因素，利用蒙特卡羅法對(duì)飛行情況進(jìn)行大數(shù)據(jù)量的仿真實(shí)驗(yàn)從而提取多尾極值參數(shù)；而后在驗(yàn)證了一維的極值參數(shù)符合GEV分布的基礎(chǔ)上，基于極值Copula模型描述了多元極值參數(shù)的分布結(jié)構(gòu)和相關(guān)性，提出了三維極值的FAWP Copula模型，驗(yàn)證了FAWP Copula在辨識(shí)具有后尾特性的多維極值參數(shù)分布時(shí)相比于其他Copula模型有更高的精度。

（2）利用FAWP Copula量化地描述了尾流飛行風(fēng)險(xiǎn)概率，對(duì)于尾流場(chǎng)內(nèi)的導(dǎo)航控制與風(fēng)險(xiǎn)規(guī)避，機(jī)場(chǎng)起降尾流安全間隔改進(jìn)，環(huán)境風(fēng)險(xiǎn)可視化等研究方向有一定的參考價(jià)值，有助于提高航空器的運(yùn)行安全性。

（3）文中思路亦可對(duì)其他內(nèi)外部環(huán)境因素影響下的飛行風(fēng)險(xiǎn)概率評(píng)估提供參考。其中的飛行風(fēng)險(xiǎn)概率評(píng)估方法是對(duì)現(xiàn)有各類(lèi)飛行安全規(guī)范［15］中風(fēng)險(xiǎn)評(píng)估理論的有效補(bǔ)充，對(duì)于飛行安全與適航性管理具有積極的作用。本文思路及方法不僅局限于尾流飛行風(fēng)險(xiǎn)的定量評(píng)估，也可以用來(lái)評(píng)估其他有飛參極值數(shù)據(jù)的情況，比如：危險(xiǎn)科目下的試飛風(fēng)險(xiǎn)、復(fù)雜外部環(huán)境下的飛行風(fēng)險(xiǎn)、飛機(jī)軟件或硬件故障下的飛行風(fēng)險(xiǎn)等等。

［1］SAE ARP4761．Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment［R］．1996．

［2］SAE ARP4754．Certification considerations for high-integrated or complex aircraft systems［R］．1996．

［3］DOD．Airworthiness certification criteria［R］．2005．

［4］DOD．Standard practice for system safety［R］．2000．

［5］GJB900-90．General program for system safety［S］．1990．（GJB900-90．系統(tǒng)安全性通用大綱，1990．）

［6］Julien S，Dimitri N M，Ivan Y B．Use of flight simulation in early design：formulation and application of the virtual testing and evaluation methodology［R］．AIAA-2000- 5590，2000．

［7］Ivan Y B，Daniel A D，Dimitri N M．Modeling and simulation of airworthiness requirements for an HSCT prototype in early design［R］．AIAA-1998- 4936，1998．

［8］Dimitri N M，Patrick T B，Neil R W．Advanced design of complex systems using the collaborative visualization environment（Co VE）［R］．AIAA-2005- 126，2005．

［9］Belcastro C M，Celeste M．On the validation of safety critical aircraft systems-partⅠ：analytical＆simulation methods［C］∥Proc.of the Conference on AIAA Guidance，Navigaion and Control，2003．

［10］Shi X W，F(xiàn)ang W H，Lin W，et al．Uncertainty of china typhoon rainfall probability estimation with different extreme value models［J］．Journal of Beijing Normal University，2011，47（5）：493- 501．（石先武，方偉華，林偉，等．基于極值理論的中國(guó)臺(tái)風(fēng)降水分布不確定性分析［J］．北京師范大學(xué)學(xué)報(bào)，2011，47（5）：493- 501．）

［11］Liu D F，Wen S Q，Wang L P．Poison-gumble mixed compound distribution and its application［J］．Chinese Science Bulletin，2002，47（22）：1901- 1906．

［12］Ho L C，Burridge P，Caddle J，et al．Value-at-risk：applying the extreme value approach to Asian markets in the recent financial turmoil［J］．Pacific-Basin Finance Journal，2000，8（2）：249- 275．

［13］Wu Z K，Zhao L，Ge Y J．Statistical analysis of wind velocity and rainfall intensity joint probability distribution of Shanghai area in typhoon condition［J］．Acta Aerodynamica Sinica，2010，28（4）：393- 399．（武占科，趙林，葛耀君．上海地區(qū)臺(tái)風(fēng)條件風(fēng)速和雨強(qiáng)聯(lián)合概率分布統(tǒng)計(jì)［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28（4）：393- 399．）

［14］Stuart C．An introduction to statistical modeling of extreme value［M］．London：Springer，2007．

［15］Nelsen R B．An introduction to copulas［M］．2nd ed．New York：Springer，2006．

［16］Joe H．Asymptotic efficiency of the two-stage estimation method for Copula-based models［J］．Journal of Multivariate Analysis，2005，94（2）：401- 419．

［17］Wang L F．The research on estimation of distribution algorithm based on Copula theory［D］．Lanzhou：Lanzhou University of Technology，2011．（王麗芳．基于Copula理論的分布估計(jì)算法研究［D］．蘭州：蘭州理工大學(xué)，2011．）

［18］Zhao L Q．The study of financial risk measurement based on copula function［D］．Xiamen：Xiamen University，2009．（趙麗琴．基于Copula函數(shù)的金融風(fēng)險(xiǎn)度量研究［D］．廈門(mén)：廈門(mén)大學(xué)，2009．）

［19］Zhang Y Z，Zhen R，Sheng G X，et al．Failure dependency of CNC equipment based on copula theory［J］．Journal of Jilin University（Engineering and Technology Edition），2011，41（6）：1636- 1640．（張英芝，鄭銳，申桂香，等．基于Copula理論的數(shù)控裝備故障相關(guān)性［J］．吉林大學(xué)學(xué)報(bào)（工學(xué)版），2011，41（6）：1636- 1640．）

［20］GJB 625A-2006．Complex subjects of flight test research for military airplane［S］．2006（GJB 625A-2006．軍用飛機(jī)科研試飛復(fù)雜科目，2006．）

［21］Parsopoulos K E，Vrahatis M N．Recent approaches to global optimization problems through particle swarm optimization［J］．Natural Computing，2002，1（1）：235- 306．

Quantitative flight risk evaluation method based on multi-dimensional extreme parameters

XUE Yuan1，XU Hao-jun1，PEI Bin-bin1，CHEN Yi-ran2
（1.Aeronautics and Astronautics Engineering College，Air Force Engineering University，Xi’an 710038，China；2.Air Force Aviation Equipment Training Base，Beijing 100076，China）

A new flight risk assessment approach based on multidimensional extreme Copula is proposed using multivariate extreme value theory and coupled system modeling ideas．First，we extract three-dimensional wake extreme parameters required for assessing the risk using Monte Carlo method，verify the extracted extreme parameters and the test data has the same distribution form，then build a flight risk determination condition；Second，we propose the four adaptive weight parameters（FAWP）for three-dimensional extreme parameters based on the result that the one-dimensional extreme parameters meet generalized extreme value distribution；Third，adaptive range particle swarm optimization algorithm is used to identify unknown parameters of the one-dimensional and three-dimensional objective function．The results of fitting test show FAWP Copula model has higher accuracy than the other Copula models，so it is the most suitable model to describe the thick tail formed by multi-dimensional extreme values．At last，the risk probability in the situation of near-ground wake encounter is evaluated using FAWPCopula，and it has some certain reference values for research directions such as wake navigation control and risk aversion．

multi-dimensional extreme parameters；generalized extreme value distribution；Copula model；adaptive range particle swarm optimization algorithm；flight risk probability

O 213．2；V 328．5

A

10．3969／j．issn．1001-506X．2015．01．18

薛 源（1986-），男，博士研究生，主要研究方向?yàn)轱w行仿真、飛行安全。

E-mail：wowszxy＠163．com

徐浩軍（1965-），男，教授，碩士，主要研究方向?yàn)轱w行安全、飛行風(fēng)險(xiǎn)。

E-mail：xuhaojun＠xjtu．edu．cn

裴彬彬（1990-），男，碩士研究生，主要研究方向?yàn)轱w行力學(xué)建模。E-mail：szxy1986＠126．com

陳怡然（1986-），女，講師，碩士，主要研究方向?yàn)轱w參數(shù)據(jù)采集。

E-mail：chenyiran＠126．com

1001-506X（2015）01-0109-08

網(wǎng)址：www．sys-ele．com

2013- 09- 26；

2014- 04- 22；網(wǎng)絡(luò)優(yōu)先出版日期：2014- 06- 17。

網(wǎng)絡(luò)優(yōu)先出版地址：http：∥w ww．cnki．net／kcms／detail／11．2422．TN．20140617．1639．009．html

國(guó)家自然科學(xué)基金（U1333131，61374145）資助課題