薛 源 徐浩軍 李 強(qiáng)

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

侯世芳

(中航工業(yè)陜西飛機(jī)工業(yè)集團(tuán)有限公司 飛行試驗(yàn)中心，漢中723000)

劉仕杰

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

當(dāng)前，國(guó)內(nèi)外對(duì)尾流飛行風(fēng)險(xiǎn)的研究多為尾流機(jī)理特征分析及尾流場(chǎng)內(nèi)的飛行動(dòng)力學(xué)計(jì)算.例如:利用仿真及實(shí)驗(yàn)手段分析了尾流的產(chǎn)生機(jī)理、動(dòng)力學(xué)特性及其在外部環(huán)境下的耗散及衰減過程，并在此基礎(chǔ)上提出了多種適用性較強(qiáng)的尾流模型[1－6];對(duì)尾流遭遇情形下的飛行狀況進(jìn)行了仿真分析[7－10].與此同時(shí) NASA，F(xiàn)AA，ICAO 等機(jī)構(gòu)為改進(jìn)機(jī)場(chǎng)的尾流安全間隔標(biāo)準(zhǔn)亦做了一些工作[11－14]，取得的研究成果已經(jīng)開始逐步應(yīng)用于航空管制和機(jī)場(chǎng)起降調(diào)度.在國(guó)內(nèi)，尾流問題也日益引起重視，一些研究人員對(duì)尾流問題均進(jìn)行過相關(guān)的研究[15－18].但當(dāng)前的研究在尾流飛行風(fēng)險(xiǎn)發(fā)生概率的評(píng)估與預(yù)測(cè)方面缺乏相關(guān)的理論與方法，鮮有構(gòu)建尾流風(fēng)險(xiǎn)概率拓?fù)淇臻g結(jié)構(gòu)的相關(guān)工作，而尾流風(fēng)險(xiǎn)概率拓?fù)涞臉?gòu)建對(duì)于環(huán)境風(fēng)險(xiǎn)的可視化具有重要的意義(NASA等多家機(jī)構(gòu)目前將環(huán)境風(fēng)險(xiǎn)可視化列為未來(lái)航空技術(shù)發(fā)展的重要方向之一).因此有必要在對(duì)飛行風(fēng)險(xiǎn)概率進(jìn)行量化評(píng)估的基礎(chǔ)上構(gòu)建尾流風(fēng)險(xiǎn)概率拓?fù)?

要構(gòu)建尾流風(fēng)險(xiǎn)概率拓?fù)涫紫软毲蟮梦擦鲌?chǎng)中每個(gè)點(diǎn)上的飛行風(fēng)險(xiǎn)量化概率，此概率是一個(gè)重要的參考指標(biāo)，對(duì)于飛機(jī)的適航性與飛行安全具有重要的意義，但怎樣求得尾流遭遇情形下的飛行風(fēng)險(xiǎn)概率尚缺乏相關(guān)的理論與方法.對(duì)于評(píng)估由于硬件故障導(dǎo)致飛行事故概率，在 ARP－4761[19]，ARP－4754[20]，MIL－HDBK－516B[21]，MILSTD－882D[22]等安全性規(guī)范中有明確的根據(jù)系統(tǒng)硬件故障實(shí)驗(yàn)樣本進(jìn)行量化計(jì)算的思路及方法，但對(duì)于由外部環(huán)境(尤其是在低空近距尾流遭遇的情況下)導(dǎo)致的飛行事故及飛行風(fēng)險(xiǎn)較難進(jìn)行量化概率指標(biāo)的評(píng)估預(yù)測(cè).因此如何利用有限的樣本參數(shù)評(píng)估計(jì)算尾流導(dǎo)致的飛行風(fēng)險(xiǎn)概率是構(gòu)建尾流風(fēng)險(xiǎn)概率拓?fù)鋱D首先要解決的問題.

利用樣本參數(shù)評(píng)估飛行風(fēng)險(xiǎn)概率首先需建立描述此樣本分布的理論模型.由于尾流飛行風(fēng)險(xiǎn)屬于低頻高危事件的范疇，其判定風(fēng)險(xiǎn)發(fā)生所提取的極值樣本一般具有厚尾特性，這種分布形式在其他的低頻高危事件(如:地震、海嘯、金融風(fēng)險(xiǎn)、飛行事故)中較為常見[23－26]，針對(duì)此種分布形式的描述目前較有效的方法是采用極值理論[27].但文中尾流飛行風(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ì)多維空間的分類，其實(shí)質(zhì)就是構(gòu)造相關(guān)性核函數(shù).文中使用Copula理論[28－32]描述多維尾流極值間的相關(guān)結(jié)構(gòu)，由于Copula理論是專為評(píng)估極值分布而提出的，故在對(duì)極值相關(guān)性的描述上，Copula極值模型能較好地反映極值參數(shù)之間的聯(lián)系和發(fā)展趨勢(shì)，比其他方法具有較高的精度，較適合應(yīng)用于本文當(dāng)中.

綜上所述，本文擬采用人－機(jī)－環(huán)復(fù)雜系統(tǒng)仿真與多元極值的風(fēng)險(xiǎn)評(píng)估方法相結(jié)合的手段，探索基于多元Copula的飛行風(fēng)險(xiǎn)概率量化評(píng)估與預(yù)測(cè)方法，在有效解決一維極值模型局限性的基礎(chǔ)上利用多維極值分布求得尾流空間中的飛行風(fēng)險(xiǎn)概率;而后從基于風(fēng)險(xiǎn)發(fā)生概率的數(shù)學(xué)層面構(gòu)建尾流發(fā)展不同階段的飛行風(fēng)險(xiǎn)概率的拓?fù)浣Y(jié)構(gòu)，并對(duì)其進(jìn)行可視化處理.研究成果擬在尾流風(fēng)險(xiǎn)定量評(píng)估和風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)可視化方面有所創(chuàng)新，可為尾流風(fēng)險(xiǎn)的預(yù)測(cè)、規(guī)避、控制提供參考;同時(shí)為尾流事故預(yù)防提供分析和檢驗(yàn)依據(jù).

1 基于蒙特卡洛法的極值提取流程

尾流遭遇情形下內(nèi)外部影響因素具有復(fù)雜隨機(jī)的特性，主要體現(xiàn)在復(fù)雜的非線性物理特征以及復(fù)雜的隨機(jī)特性和不確定性.由于文中研究?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)條件苛刻;因此文中基于蒙特卡洛法考慮尾流遭遇風(fēng)險(xiǎn)的隨機(jī)性，對(duì)三維空間內(nèi)的每個(gè)尾流遭遇點(diǎn)進(jìn)行大數(shù)據(jù)量的計(jì)算從而提取多元極值參數(shù)，有效地對(duì)尾流遭遇情形的隨機(jī)性和復(fù)雜性進(jìn)行描述.

首先在尾流場(chǎng)實(shí)驗(yàn)數(shù)據(jù)和數(shù)值計(jì)算數(shù)據(jù)的基礎(chǔ)上建立尾流模型庫(kù).數(shù)據(jù)庫(kù)以經(jīng)過驗(yàn)證的尾流理論模型[6－7]為基礎(chǔ)，以特定機(jī)翼流場(chǎng)后的數(shù)值計(jì)算結(jié)果作為常用飛機(jī)尾流場(chǎng)數(shù)據(jù)的補(bǔ)充.檢索輸入為前方產(chǎn)生尾流飛機(jī)的型號(hào)與尾流場(chǎng)內(nèi)的三維坐標(biāo)，輸出為該三維空間點(diǎn)上的尾流場(chǎng)參數(shù).尾流數(shù)據(jù)庫(kù)將大氣紊流程度、溫度分層、離散突風(fēng)、Crow 鏈接因子等[1－3，6，8－10]對(duì)尾流場(chǎng)狀態(tài)有影響的內(nèi)外部因素設(shè)置為影響變量，在分析其對(duì)尾流場(chǎng)的各典型參數(shù)(衰減率、渦強(qiáng)、渦核內(nèi)徑、速度場(chǎng))和拓?fù)湫螤钣绊懙幕A(chǔ)上，將內(nèi)外部因素作為隨機(jī)影響項(xiàng)按照其出現(xiàn)概率補(bǔ)充進(jìn)數(shù)據(jù)庫(kù)，使得尾流數(shù)據(jù)庫(kù)可以在蒙特卡洛仿真過程中的每一次檢索都能反映內(nèi)外部環(huán)境影響下的尾流場(chǎng)隨機(jī)特性和動(dòng)態(tài)特性.

而后構(gòu)建飛行員補(bǔ)償操縱行為特征數(shù)據(jù)庫(kù).由于飛行員的心理素質(zhì)和操縱習(xí)慣不一樣，對(duì)新環(huán)境的適應(yīng)能力和感覺能力也不盡相同，無(wú)論是在常規(guī)操作時(shí)還是處理某些特殊情況，飛行員的操縱行為都具有隨機(jī)性的特點(diǎn).通過研究人的心理活動(dòng)、生理活動(dòng)，對(duì)飛行員操縱模型進(jìn)行研究是非常有效的.第四軍醫(yī)大學(xué)航空航天醫(yī)學(xué)系的裝備教研室對(duì)飛行員進(jìn)行了大量的實(shí)驗(yàn)采樣，本文根據(jù)飛行員的實(shí)驗(yàn)采樣數(shù)據(jù)辨識(shí)了駕駛員模型中的操縱行為參數(shù)，建立了反映飛行員補(bǔ)償操縱隨機(jī)性和不確定性的行為特征數(shù)據(jù)庫(kù).

在計(jì)算飛參數(shù)據(jù)時(shí)，利用蒙特卡洛方法將內(nèi)外部影響條件變量按照其出現(xiàn)頻率進(jìn)行隨機(jī)抽樣，需抽樣的條件變量如表1所示.將抽樣的變量數(shù)值動(dòng)態(tài)代入到實(shí)驗(yàn)系統(tǒng)的計(jì)算迭代中，以對(duì)相關(guān)的氣動(dòng)參數(shù)及操縱信號(hào)產(chǎn)生量化影響，從而反映真實(shí)內(nèi)外部環(huán)境影響下的隨機(jī)性與不確定性.極值參數(shù)提取實(shí)驗(yàn)系統(tǒng)基于某型飛機(jī)鐵鳥地面實(shí)驗(yàn)系統(tǒng)改造而成，實(shí)物如圖1所示.

圖1 提取極值參數(shù)實(shí)驗(yàn)系統(tǒng)Fig.1 Simulation system for extracting extreme parameters

表1 需用蒙特卡洛法抽樣的影響條件變量Table 1 Condition variables need to be sampled using Monte Carlo method

在單個(gè)仿真實(shí)驗(yàn)特定點(diǎn)上提取極值參數(shù)的步驟如下:

1)設(shè)置前機(jī)和后機(jī)的機(jī)體參數(shù)及前機(jī)的初始飛行狀態(tài);

2)設(shè)置仿真實(shí)驗(yàn)特定點(diǎn);

3)設(shè)置i=1;

4)根據(jù)各個(gè)影響條件變量的統(tǒng)計(jì)頻率，利用蒙特卡洛法抽樣表1中的條件變量值;

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

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

7)i=i+1，回到步驟4)，循環(huán)迭代，直至i=n.

其中n為設(shè)置的仿真特定點(diǎn)上的蒙特卡洛計(jì)算次數(shù).根據(jù)不同計(jì)算次數(shù)下的仿真結(jié)果，n越大，越能反映極值的統(tǒng)計(jì)特性，但提高n亦會(huì)增加系統(tǒng)的計(jì)算負(fù)擔(dān)，經(jīng)驗(yàn)證在n＞2 000時(shí)，計(jì)算結(jié)果趨于相對(duì)穩(wěn)定，因此根據(jù)系統(tǒng)計(jì)算資源的實(shí)際使用狀況，設(shè)置n=3000.

2 所提取極值參數(shù)的可信度驗(yàn)證

設(shè)置前機(jī)速度為90 m/s，飛行高度為100 m，前機(jī)的初始渦量為430m2/s(相當(dāng)于A340的尾渦初始強(qiáng)度)，后機(jī)進(jìn)入尾流場(chǎng)時(shí)的襟翼角度為0°.以t=17 s時(shí)的尾流左渦核中心偏左5 m處為仿真實(shí)驗(yàn)特定點(diǎn)，由于前機(jī)速度設(shè)置為90 m/s，故坐標(biāo)x=1 530 m.對(duì)此狀態(tài)通過上文中的步驟進(jìn)行n=3000次的計(jì)算迭代，圖2顯示了i=61的仿真過程中，出現(xiàn)的飛參極值超限從而導(dǎo)致發(fā)生文中所定義的飛行風(fēng)險(xiǎn).可以看出最初進(jìn)入尾流場(chǎng)時(shí)由于尾流中上升氣流的作用，高度和迎角均略有增加，隨后由于滾轉(zhuǎn)力矩的作用，飛機(jī)急劇滾轉(zhuǎn)并掉高度.確定對(duì)尾流飛行風(fēng)險(xiǎn)發(fā)生影響最大的3個(gè)飛行參數(shù)(滾轉(zhuǎn)角φ、下降高度ΔH、迎角α).本次迭代的極值參數(shù).因本文選取的近地近距尾流遭遇情形屬于高風(fēng)險(xiǎn)科目，故不可能采用試飛驗(yàn)證，因此采用飛行員在回路的飛機(jī)地面鐵鳥臺(tái)實(shí)驗(yàn)數(shù)據(jù)作為驗(yàn)?zāi)?shù)據(jù).某型飛機(jī)鐵鳥臺(tái)是飛機(jī)的地面實(shí)驗(yàn)臺(tái)，軟硬件結(jié)構(gòu)均與真實(shí)飛機(jī)相同，使用經(jīng)過風(fēng)洞實(shí)驗(yàn)與試飛驗(yàn)證后的氣動(dòng)參數(shù)，在各個(gè)飛行科目下與真實(shí)試飛數(shù)據(jù)的誤差不超過12%.相比于其他類型的飛機(jī)地面模擬器具有較高的精度，可以認(rèn)為鐵鳥臺(tái)系統(tǒng)所模擬的科目符合真實(shí)試飛狀況.

圖2 第61次迭代中的飛行參數(shù)Fig.2 Flight parameters of the 61th simulation process

列出在i≤75時(shí)提取的前75個(gè)極值參數(shù)如圖3所示，圖3中箭頭所指示即為i=61時(shí)的極值參數(shù).圖4為使用上文方法提取的前75次極值參數(shù)以及在相同條件下進(jìn)行的75次飛行員在回路鐵鳥臺(tái)實(shí)驗(yàn)數(shù)據(jù)的QQ(Quantile－Quantile)圖，可發(fā)現(xiàn)3種極值參數(shù) φmax，ΔHmax，αmax的 QQ 圖均近似為直線，說明本文提取的極值參數(shù)和實(shí)驗(yàn)數(shù)據(jù)屬于同一種分布類型.Kolmogorov－Smirnov(KS)檢驗(yàn)的結(jié)果亦表明3種極值參數(shù)的K－S值均小于0.1，而P值均大于0.25(即在比95%的置信水平低得多的情況下亦能通過檢驗(yàn))，故可認(rèn)為本文方法得到的數(shù)據(jù)與鐵鳥臺(tái)實(shí)驗(yàn)數(shù)據(jù)具有相同的分布類型.因此文中方法所提取極值參數(shù)的可信度較高，可以接受其作為評(píng)估飛行風(fēng)險(xiǎn)概率的樣本.從圖3中可以看出極值參數(shù)的分布存在較明顯的厚尾特性，又因?yàn)槲擦黠L(fēng)險(xiǎn)屬于低頻高危事件的范疇，故采用極值理論對(duì)此種分布形式進(jìn)行描述.利用極值理論能夠有效地對(duì)隨機(jī)序列最值概率分布的尾部進(jìn)行建模，用于描述極值樣本數(shù)據(jù)序列分布的尾部特征.

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

圖3 通過蒙特卡洛仿真和飛機(jī)地面實(shí)驗(yàn)分別得出的前75次三維極值參數(shù)Fig.3 Front 75 groups of extreme parameters through Monte Carlo simulation and aircraft ground test

圖4 極值參數(shù)和實(shí)驗(yàn)數(shù)據(jù)的QQ圖Fig.4 QQ plots of simulated extreme parameters and test extreme parameters

對(duì)3個(gè)極值參數(shù)進(jìn)行歸一化處理.查某型飛機(jī)氣動(dòng)數(shù)據(jù)，極值參數(shù)迎角的臨界值與Ma有關(guān)，例如在襟翼偏角δf=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/85.重心下降高度的極值參數(shù)為ΔHmax，以機(jī)翼翼尖剛好觸地時(shí)的狀態(tài)作為風(fēng)險(xiǎn)發(fā)生臨界點(diǎn)，極值參數(shù)ΔHmax的歸一化公式如式(1)所示.

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

3 基于Copula模型的三維極值分布

多維隨機(jī)變量的聯(lián)合分布除了與各分量的分布有關(guān)之外，更重要的是與變量之間的相關(guān)結(jié)構(gòu)有關(guān).當(dāng)隨機(jī)變量的維數(shù)比較大時(shí)，單個(gè)分量的極值行為未必含有整個(gè)向量的聯(lián)合極值行為.

3.1 提出FPAVW Copula模型

假設(shè)函數(shù)F是多元極值分布函數(shù).根據(jù)Copula的相關(guān)定理[28]:如果F是多元極值分布函數(shù)，則F的一維邊緣分布必然屬于廣義極值(GEV)分布族(如式(3)).

其中，ξ∈R;μ∈R;σ ＞0;1+ξ(x－ μ)/σ ＞0.

設(shè)文中極值向量(φmax，ΔHmax，αmax)的分布函數(shù)為 F(φmax，ΔHmax，αmax)，邊緣分布函數(shù)為F1(φmax)，F(xiàn)2(ΔHmax)和 F3(αmax)，經(jīng)驗(yàn)證，F(xiàn)1(φmax)，F(xiàn)2(ΔHmax)和 F3(αmax)在比 95% 的置信水平低得多的情況下亦能通過分布形式為GEV的假設(shè)檢驗(yàn)，故其符合極值理論中的GEV分布，故對(duì)于任意的(φmax，ΔHmax，αmax)∈Rd，一定存在一個(gè)Copula C，使得

文中的 F1(φmax)，F(xiàn)2(ΔHmax)和 F3(αmax)都是連續(xù)分布函數(shù)，故C是唯一的.由式(4)定義的函數(shù)C是一個(gè)邊緣分布為GEV形式的三維聯(lián)合分布函數(shù).對(duì)于本文中三維極值參數(shù)的Copula模型選擇，首先根據(jù)常用的三維Copula來(lái)構(gòu)建，其通用形式如式(5)所示.主要有Gumbel Copula模型(如式(6))、Frank Copula模型(如式(7))、Clayton Copula模型(如式(8))、GS Copula模型(如式(9))、Joe Copula模型(如式(10)).

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

3.2 辨識(shí)Copula模型中的未知參數(shù)

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

1)根據(jù)定義，Copula的邊緣分布函數(shù)即為一維極值的GEV分布函數(shù)，故將每組三維極值樣本點(diǎn)分別代入式(3)中，辨識(shí)出的GEV分布模型中未知參數(shù):

得到極值參數(shù)最終的邊緣累積概率如式(12)～式(14)所示.

2)根據(jù)式(15)求出Copula的密度函數(shù).

3)根據(jù)變量ui，vi，wi的數(shù)值構(gòu)建目標(biāo)函數(shù)如式(16)所示.

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

4)對(duì)目標(biāo)函數(shù)進(jìn)行辨識(shí)，在對(duì)最小二乘法、極大似然法、模式搜索算法、遺傳算法(GA，Genetic Algorithm)等辨識(shí)方式進(jìn)行對(duì)比測(cè)試的基礎(chǔ)上，結(jié)合目標(biāo)函數(shù)(式(16))構(gòu)造復(fù)雜及計(jì)算量較大的特點(diǎn).選用較適合本文的優(yōu)化環(huán)境的局部搜索能力強(qiáng)且收斂速度較快的自適應(yīng)區(qū)間粒子群(ARPSO，Adaptive Range Particle Swarm Optimization)算法[33].利用 ARPSO 辨識(shí)出目標(biāo)函數(shù)式(16)的未知參數(shù)如表2所示.

表2 不同Copula模型的參數(shù)辨識(shí)結(jié)果Table 2 Identified weight parameters of different Copula models

根據(jù)表2中6種Copula的辨識(shí)參數(shù)構(gòu)建其概率密度如圖5所示，以此更直觀地觀察Copula模型的密度函數(shù)特征，圖5中U和V分別表示ui和vi的連續(xù)形式.

從圖 5中亦可以看出 Gumbel Copula，Joe Copula，F(xiàn)PAVW 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é)果如表3所示.

圖5 w=0.85時(shí)不同Copula模型的概率分布圖Fig.5 Probability density maps of different Copula models when w=0.85

表3 不同Copula模型的描述優(yōu)度檢驗(yàn)Table 3 Fitting test for different Copula models

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

最終，根據(jù)構(gòu)建的FPAVW Copula模型求出風(fēng)險(xiǎn)概率，如式(17)所示.

根據(jù)表2中辨識(shí)出的FPAVW Copula未知參數(shù)求得在第1節(jié)提到的仿真特定點(diǎn)上的飛行風(fēng)險(xiǎn)概率Pr=0.051 2.需注意的是，由于飛行事故的發(fā)生是一個(gè)多因素影響的不確定過程，不可能將所有內(nèi)外部隨機(jī)因素考慮完全，因此4761[19]與882D[22]中的事故率很大程度上是一個(gè)參考值.本文得到的飛行風(fēng)險(xiǎn)量化概率值在多數(shù)狀況下亦是一個(gè)參考值，與真實(shí)值必然有一定的誤差.但其在不同狀況下飛行風(fēng)險(xiǎn)的橫向?qū)Ρ确治?、風(fēng)險(xiǎn)程度的歸類劃分中具有積極的意義.如:不同尾流遭遇情形下的飛行風(fēng)險(xiǎn)定量比較，不同惡劣環(huán)境條件下或不同硬件故障條件下的風(fēng)險(xiǎn)大小對(duì)比，預(yù)定科目或任務(wù)的風(fēng)險(xiǎn)程度預(yù)測(cè)比較等.

4 構(gòu)建三維尾流飛行風(fēng)險(xiǎn)概率拓?fù)鋱D

上文以尾流三維空間內(nèi)的單個(gè)特定點(diǎn)為例闡述了求取尾流風(fēng)險(xiǎn)概率的方法.下面將利用上文方法對(duì)整個(gè)尾流場(chǎng)三維空間進(jìn)行飛行風(fēng)險(xiǎn)的量化概率計(jì)算，而后構(gòu)建尾流不同發(fā)展階段的風(fēng)險(xiǎn)概率拓?fù)鋱D.首先對(duì)不同發(fā)展階段尾流場(chǎng)的空間結(jié)構(gòu)進(jìn)行網(wǎng)格化，每個(gè)網(wǎng)格節(jié)點(diǎn)對(duì)應(yīng)一個(gè)數(shù)據(jù)采樣點(diǎn)，并依據(jù)尾流場(chǎng)強(qiáng)度的大小對(duì)網(wǎng)格進(jìn)行逐級(jí)加密.將這些采樣點(diǎn)存入SQLSERVER數(shù)據(jù)庫(kù)，供服務(wù)器實(shí)時(shí)調(diào)用.風(fēng)險(xiǎn)拓?fù)浣Y(jié)構(gòu)圖的構(gòu)建流程如圖6所示.

圖6 尾流風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)圖的構(gòu)建Fig.6 Construction process of flight risk probability distribution maps in wakes

按照尾流場(chǎng)的發(fā)展變化規(guī)律，為風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)分為3個(gè)階段.

4.1 第1階段:尾流風(fēng)險(xiǎn)概率拓?fù)涑跏茧A段

在尾流的初始衰減階段，繪制尾流風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)圖如圖7所示.x軸為前機(jī)飛離的距離.圖7c 為 Schwarz，Hahn，Holz?pfel，Gerz 等人提出的尾流場(chǎng)內(nèi)簡(jiǎn)化風(fēng)險(xiǎn)區(qū)間(SHA，Simplified Hazard Area)[34－37]示意圖.

圖7 1.35 km≤x≤1.65 km時(shí)的風(fēng)險(xiǎn)概率拓?fù)鋱DFig.7 Flight risk probability topology maps when 1.35 km≤x≤1.65 km

從圖7b中可以看出，本文構(gòu)建的二維風(fēng)險(xiǎn)概率拓?fù)鋱D在形狀上與圖7c中的SHA有相似之處.但是由于本文考慮了隨機(jī)性和不確定性的影響，二維風(fēng)險(xiǎn)概率拓?fù)鋱D呈現(xiàn)不規(guī)則且邊緣粗糙的蝴蝶狀，更加符合真實(shí)的情況.與SHA規(guī)則、光滑且左右對(duì)稱的形狀是有顯著區(qū)別的.另一方面，SHA是通過簡(jiǎn)單計(jì)算尾流誘導(dǎo)滾轉(zhuǎn)力矩與最大可操縱滾轉(zhuǎn)能力的比值來(lái)確定的.而文中的風(fēng)險(xiǎn)概率拓?fù)鋱D則是基于定量的飛行風(fēng)險(xiǎn)概率構(gòu)建的，尾流場(chǎng)中的每一個(gè)網(wǎng)格點(diǎn)都有其對(duì)應(yīng)的概率值.這些定量的概率值對(duì)于尾流遭遇情形下的飛行風(fēng)險(xiǎn)評(píng)估與預(yù)測(cè)具有嚴(yán)格且科學(xué)的統(tǒng)計(jì)學(xué)意義，可以對(duì)現(xiàn)有飛行安全準(zhǔn)則中的飛行事故率評(píng)估方法提供有效的補(bǔ)充，這也是SHA所不具有的.

分析圖7可以看出，在左右主風(fēng)險(xiǎn)區(qū)域的兩側(cè)各有一個(gè)較小的次級(jí)風(fēng)險(xiǎn)區(qū)域.在1.35 km≤x≤1.65 km時(shí)尾渦中心的飛行風(fēng)險(xiǎn)概率大于13%，此為高度危險(xiǎn)狀態(tài);圖7b淺色區(qū)域內(nèi)的飛行風(fēng)險(xiǎn)發(fā)生概率均大于2.5%，此概率值亦遠(yuǎn)遠(yuǎn)高于 ARP－4761[19]所規(guī)定的系統(tǒng)安全性水平中的A級(jí)災(zāi)難性和B級(jí)危險(xiǎn)性失效狀態(tài)，因此無(wú)論何種情況，都應(yīng)避免進(jìn)入圖7b中的淺色區(qū)域.分析其拓?fù)浣Y(jié)構(gòu)亦可發(fā)現(xiàn)隨著x的增加，飛行風(fēng)險(xiǎn)有逐漸減小的趨勢(shì).

4.2 第2階段:風(fēng)險(xiǎn)概率拓?fù)涞拈L(zhǎng)波振蕩階段

隨著尾流出現(xiàn)長(zhǎng)波不穩(wěn)定現(xiàn)象，尾流風(fēng)險(xiǎn)概率拓?fù)渲械闹黠L(fēng)險(xiǎn)區(qū)域和左右次級(jí)風(fēng)險(xiǎn)區(qū)域呈現(xiàn)周期性的振蕩，如圖8、圖9所示.從圖8和圖10可以看出在震蕩過程中左右風(fēng)險(xiǎn)區(qū)域相互靠攏后又分離，其離得越近，風(fēng)險(xiǎn)程度越高.風(fēng)險(xiǎn)概率在核心處達(dá)到8% ～8.5%，相比第1階段有所減小.

圖8 5.4 km≤x≤5.7 km時(shí)的三維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.8 Three－dimengional flight risk probability topology map when 5.40 km≤x≤5.70 km

圖9 11.70≤x≤12.0 km時(shí)的三維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.9 Three－dimensional flight risk probability topology map when 11.7 km≤x≤12.0 km

從圖9和圖11可以看出隨著x的進(jìn)一步增大，尾流風(fēng)險(xiǎn)概率拓?fù)涞恼袷幏狄嚯S著變大.在x=11.8 km左右時(shí)，幅值達(dá)到了40 m以上.在11.7 km≤x≤12.0 km時(shí)，主風(fēng)險(xiǎn)核心分化出較多的次級(jí)風(fēng)險(xiǎn)區(qū)域，在波谷處又重新融合為一個(gè)較大的主風(fēng)險(xiǎn)核心(如圖11).隨著尾渦向著融合的方向發(fā)展，尾流風(fēng)險(xiǎn)拓?fù)浣Y(jié)構(gòu)中的左右風(fēng)險(xiǎn)區(qū)域亦開始相互融合與鏈接，圖11顯示了在左右風(fēng)險(xiǎn)區(qū)域靠近的過程中，風(fēng)險(xiǎn)核心被顯著拉長(zhǎng).

圖11中的尾流風(fēng)險(xiǎn)區(qū)域的最大風(fēng)險(xiǎn)概率下降到4%以下且風(fēng)險(xiǎn)核心的覆蓋范圍變小，相比于第1階段大于13%的風(fēng)險(xiǎn)概率有了顯著的下降.這和尾渦發(fā)展過程中的長(zhǎng)波不穩(wěn)定性(Crow不穩(wěn)定性)剛好是吻合的，在這個(gè)階段，尾渦強(qiáng)度受到外部環(huán)境和自身長(zhǎng)波不穩(wěn)定性的影響迅速衰減，從而亦導(dǎo)致了風(fēng)險(xiǎn)概率值的驟降.

圖10 5.4 km＜x＜5.7 km時(shí)不同位置處的二維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.10 Two－dimensional flight risk probability topology maps at different positoins when 5.4 km≤x≤5.7 km

圖11 11.7 km≤x≤12.0 km時(shí)不同位置處的二維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.11 Two－dimensional flight risk probability topology maps at different positions when 11.7 km≤x≤12.0 km

4.3 第3階段:尾流風(fēng)險(xiǎn)概率拓?fù)涞沫h(huán)狀階段

在尾流左右風(fēng)險(xiǎn)區(qū)域鏈接之后，風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)開始呈現(xiàn)不規(guī)則的環(huán)形(如圖12、圖13)，并分化出3個(gè)以上的次級(jí)風(fēng)險(xiǎn)區(qū)域，環(huán)繞在主風(fēng)險(xiǎn)區(qū)域周圍，呈現(xiàn)脊椎狀，三維風(fēng)險(xiǎn)拓?fù)洵h(huán)隨著x的增加越來(lái)越扁平且沿y軸被拉長(zhǎng).

圖12 15.3 km≤x≤16.1 km時(shí)的三維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.12 Three－dimensional flight risk probability topology map when 15.3 km≤x≤16.1 km

風(fēng)險(xiǎn)區(qū)域的拓?fù)湫螤詈偷?階段與第2階段相比進(jìn)一步拉伸變長(zhǎng)，核心區(qū)域進(jìn)一步減小.

在圖13的風(fēng)險(xiǎn)環(huán)拓?fù)鋬?nèi)部，風(fēng)險(xiǎn)概率二維拓?fù)浣Y(jié)構(gòu)呈現(xiàn)長(zhǎng)條狀，主風(fēng)險(xiǎn)區(qū)域與次級(jí)風(fēng)險(xiǎn)區(qū)域的劃分開始模糊，擁有多個(gè)風(fēng)險(xiǎn)核心區(qū)域.在風(fēng)險(xiǎn)環(huán)的兩端有上下分層的趨勢(shì).而在兩個(gè)風(fēng)險(xiǎn)環(huán)交界的中間過渡區(qū)域，由于偏離尾流風(fēng)險(xiǎn)概率拓?fù)涞挠绊懛秶L(fēng)險(xiǎn)概率值極小，平均在0.01%以下，可認(rèn)為風(fēng)險(xiǎn)環(huán)間的過渡區(qū)域是相對(duì)無(wú)風(fēng)險(xiǎn)的.

在風(fēng)險(xiǎn)概率拓?fù)洵h(huán)狀結(jié)構(gòu)形成后，風(fēng)險(xiǎn)概率減小的速率與第2階段相比有所下降.在尾流風(fēng)險(xiǎn)拓?fù)浒l(fā)展的第3階段后期，環(huán)狀拓?fù)浣Y(jié)構(gòu)內(nèi)的風(fēng)險(xiǎn)概率最大值在0.4%以下，且核心風(fēng)險(xiǎn)區(qū)較小，后機(jī)進(jìn)入此風(fēng)險(xiǎn)區(qū)域可以被認(rèn)為是安全的，但仍應(yīng)當(dāng)在操縱余度充裕的情況下盡量避免進(jìn)入該區(qū)域.

在尾流風(fēng)險(xiǎn)概率拓?fù)浒l(fā)展的3個(gè)階段中，前機(jī)右機(jī)翼范圍內(nèi)飛行風(fēng)險(xiǎn)概率值的強(qiáng)度與覆蓋區(qū)域比左翼的略微大一些.分析其原因應(yīng)為左翼范圍內(nèi)一般引起飛機(jī)的右滾轉(zhuǎn)效應(yīng)，而向左壓桿比向右拉桿更符合飛行員的操縱習(xí)慣，故飛行員在前機(jī)左翼范圍內(nèi)的補(bǔ)償操縱更為靈敏有效，因此左翼范圍內(nèi)的飛行風(fēng)險(xiǎn)較右翼小.而風(fēng)險(xiǎn)核心區(qū)下方的風(fēng)險(xiǎn)范圍比上方的風(fēng)險(xiǎn)范圍普遍大一些，這是由于隨著z軸的減小離地高度亦減小，從而留給飛行員的反應(yīng)時(shí)間減小，因此更容易引發(fā)飛行風(fēng)險(xiǎn).

圖13 15.3 km≤x≤16.1 km時(shí)不同位置處的二維風(fēng)險(xiǎn)概率拓?fù)鋱DFig.13 Two－dimensional flight risk probability topology maps at different positions when 15.3 km≤x≤16.1 km

5 結(jié)論

1)提出了一種量化評(píng)估飛行風(fēng)險(xiǎn)概率的新思路及新方法，并將其用到尾流場(chǎng)風(fēng)險(xiǎn)概率的評(píng)估中;在對(duì)尾流場(chǎng)空間內(nèi)的每個(gè)網(wǎng)格節(jié)點(diǎn)進(jìn)行飛行風(fēng)險(xiǎn)概率計(jì)算的基礎(chǔ)上，構(gòu)建了3個(gè)不同階段的二維及三維尾流風(fēng)險(xiǎn)概率拓?fù)浣Y(jié)構(gòu)圖，完成了三維空間內(nèi)的尾流飛行風(fēng)險(xiǎn)概率指標(biāo)的可視化，量化地描述了尾流飛行風(fēng)險(xiǎn)概率在尾流發(fā)展不同階段的動(dòng)態(tài)變化情況.對(duì)于尾流風(fēng)險(xiǎn)預(yù)警和規(guī)避控制技術(shù)的研究具有一定的參考價(jià)值，有助于提高航空器的運(yùn)行安全性.

2)文中風(fēng)險(xiǎn)概率拓?fù)鋱D的構(gòu)建方法是對(duì)現(xiàn)有各類飛行安全規(guī)范[19－22]中風(fēng)險(xiǎn)評(píng)估理論的有效補(bǔ)充，對(duì)于飛行安全與適航性管理具有積極的作用.

3)本文思路及方法不僅局限于尾流飛行風(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)等等.

References)

[1]Holz?pfel F.Probabilistic two－phase aircraft wake－vortex model:further development and assessment[J].Journal of Aircraft，2006，43(3):700 －708

[2]Holz?pfel F，Kladetzke J.Assessment of wake－vortex encounter probabilities for crosswind departure scenarios[J].Journal of Aircraft，2011，48(3):812 －822

[3]Sarpkaya T，Robins R E，Delisi D P.Wake－vortex eddy－dissipation model predictions compared with observations[J].Journal of Aircraft，2001，38(4):687 － 692

[4]Jurkovich M S.CFD prediction of the flow field behind the KC－135R tanker[R].AIAA 2011－3510，2011

[5]Harris M，Vaughan J M，Huenecke K E A.Aircraft wake vortices:a comparison of wind－tunnel data with field trial measurements by laser radar[J].Aerospace Science and Technology，2000，4(5):363 －370

[6]Gerz T，Holz?pfel F，Darracq D.Commercial aircraft wake vortices[J].Progress in Aerospace Sciences，2002，38(3):181 － 208

[7]Rossow V J.Lift－generated vortex wakes of subsonic transport aicraft[J].Progress in Aerospace Sciences，1999，35(6):507 －660

[8]Dogan A，Kim E，Blake W.Control and simulation of relative motion for aerial refueling in racetrack maneuvers[J].Journal of Guidance，Control and Dynamics，2007，30(5):1551 －1557

[9]Hahn K U，F(xiàn)ischenberg D，Niedermeier D，et al.Wake encounter flight control assistance based on forward－looking measurement processing[R].AIAA－2010－7680，2010

[10]Rossow V J，James K D.Overview of wake－vortex hazards during cruise[J].Journal of Aircraft，2000，37(6):960 － 975

[11]Sharon M J，Reveley M.An overview of the NASA aviation safety program assessment progress[R].AIAA－2003－6706，2003

[12]Luckner R，H?hne G，F(xiàn)uhrmann M.Hazard criteria for wake vortex encounters during approach[J].Aerospace Science and Technology，2004，8(8):673 －687

[13]De Bruin A.Wake vortex evolution and ENCounter(WAVENC)[J].Air＆ Space Europe，2000，2(5):84 －87

[14]Visscher D，Winckelmansy G，Lonfils T，et al.The WAKE4D simulation platform for predicting aircraft wake vortex transport and decay:description and examples of application[R].AIAA－2010－7994，2010

[15]周彬.飛機(jī)尾流的微結(jié)構(gòu)特征及散射特性研究[D].長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué)，2009

Zhou Bin.Study on the microstructure and scattering characteristics of aircraft wake vortices[D].Changsha:National University of Defense Technology，2009(in Chinese)

[16]李健兵，王雪松，王濤.基于馬蹄渦模型的飛機(jī)尾流介電常數(shù)分布的快速建模方法[J].航空學(xué)報(bào)，2009，30(7):1214－1218

Li Jianbing，Wang Xuesong，Wang Tao.Rapid modeling of distribution of aircraft wake dielectric constant based on horse－shoe Vortex model[J].Acta Aeronoutica et Astronautica Sinica，2009，30(7):1214 －1218(in Chinese)

[17]魏志強(qiáng).尾渦流場(chǎng)及安全間隔的建模與仿真計(jì)算研究[D].天津:中國(guó)民航大學(xué)，2008

Wei Zhiqiang.The research on modeling and simulation of flow field and safety spacing for wake vor－tex[D].Tianjin:Civil Aviation University of China，2008(in Chinese)

[18]李大偉，王宏倫.自動(dòng)空中加油階段加油機(jī)尾渦流場(chǎng)建模與仿真[J].北京航空航天大學(xué)學(xué)報(bào)，2010，36(7):776 －780

Li Dawei，Wang Honglun.Wake vortex effect modeling and simulation in automated aerial refueling[J].Journal of Beijing University of Aeronautics and Astronautics，2010，36(7):776 －780(in Chinese)

[19]SAE ARP 4761—1996 Guidelines and methods for conducting the safety assessment process on civil airborne systems and equipment[S]

[20]SAE ARP 4754—1996 Certification considerations for highintegrated or complex aircraft systems[S]

[21]DOD.MIL－HDBK－516B—2005 Airworthiness certification criteria[S]

[22]DOD.MIL－STD－882D—2000 Standard practice for system safety[S]

[23]石先武，方偉華，林偉，等.基于極值理論的中國(guó)臺(tái)風(fēng)降水分布不確定性分析[J].北京師范大學(xué)學(xué)報(bào)，2011，47(5):493－501

Shi Xianwu，F(xiàn)ang Weihua，Lin Wei，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(in Chinese)

[24]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

[25]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):249 － 275

[26]武占科，趙林，葛耀君.上海地區(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

Wu Zhanke，Zhao Lin，Ge Yaojun.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(in Chinese)

[27]Stuart C.An introduction to statistical modeling of extreme value[M].London:Springer，2007

[28]Nelsen R B.An introduction to copulas[M].2nd ed.New York:Springer，2006

[29]Joe H.Asymptotic efficiency of the two－stage estimation method for Copula－based models[J].Journal of Multivariate Analysis，2005，94(2):401 －419

[30]王麗芳.基于Copula理論的分布估計(jì)算法研究[D].蘭州:蘭州理工大學(xué)，2011

Wang Lifang.The research on estimation of distribution algorithm based on Copula theory[D].Lanzhou:Lanzhou University of Technology，2011(in Chinese)

[31]趙麗琴.基于Copula函數(shù)的金融風(fēng)險(xiǎn)度量研究[D].廈門:廈門大學(xué)，2009

Zhao Liqin.The study of financial risk measurement based on copula function[D].Xiamen:Xiamen University，2009(in Chinese)

[32]張英芝，鄭銳，申桂香，等.基于Copula理論的數(shù)控裝備故障相關(guān)性[J].吉林大學(xué)學(xué)報(bào):工學(xué)版，2011，41(6):1636 －1640

Zhang Yingzhi，Zheng Rui，Shen Guixiang，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(in Chinese)

[33]Parsopoulos K E，Vrahatis M N.Recent approaches to global optimization problems through particle swarm optimization[J].Natural Computing，2002，1(1):235 －306

[34]Schwarz C W，Hahn K U.Full－flight simulator study for wake vortex hazard area investigation[J].Aerospace Science and Technology，2006，10(2):136 －143

[35]Hahn K U，Schwarz C，F(xiàn)riehmelt H.A simplified hazard area prediction(SHAPe)model for wake vortex encounter avoidance[C]//24th International Congress of Aeronautical Sciences(ICAS)，2004

[36]Holz?pfel F，Gerz T，F(xiàn)rech M，et al.The wake vortex prediction and monitoring system WSVBS part I:design[J].Air Traffic Control Quarterly，2009，17(4):301 －322

[37]Gerz T，Holz?pfel F，Gerling W，et al.The wake vortex prediction and monitoring system WSVBS part II:performance and ATC integration at frankfurt airport[J].Air Traffic Control Quarterly，2009，17(4):323 －346