基于蒙特卡洛方法的著陸擦機尾風(fēng)險預(yù)測*

汪 磊，楊星月

(中國民航大學(xué) 飛行技術(shù)學(xué)院，天津 300000)

0 引言

飛行數(shù)據(jù)顯示，民機的進(jìn)近和著陸階段是整個飛行過程中的事故高發(fā)階段。雖然著陸階段僅占整個航段飛行時間的1%，但卻有23%的事故和17%的死亡人數(shù)[1]發(fā)生在這一階段。其中，擦機尾是一類頻發(fā)的民航不安全事件，一般定義為飛機起飛或著陸時機尾與跑道道面發(fā)生接觸的現(xiàn)象，著陸階段比起飛更容易發(fā)生。1991—2013年我國民航共發(fā)生擦機尾事件54起，其中著陸階段發(fā)生的擦機尾占52%[2]。飛機在著陸過程中，由于上仰角過大，機身離地面沒有足夠的裕度，導(dǎo)致機尾擦地。一般來說，著陸俯仰角越大，擦機尾發(fā)生的可能性越高。飛機擦機尾很少直接造成機毀人亡，初始階段只是增加航班運行和飛機維修的成本, 但其對飛機結(jié)構(gòu)的隱形損傷卻是不容忽視的，極有可能在一段時間之后成為飛行重大安全隱患之一。

目前，國外學(xué)者對擦機尾的研究主要從告警裝置設(shè)計及風(fēng)險監(jiān)測方面開展：Theriault等[3]研發(fā)了1種基于圖形化警告來提高飛行員警覺性的擦機尾事件目視警告系統(tǒng)；Chan[4]運用天氣預(yù)報(NWP)模型研究了風(fēng)切變對飛機擦機尾的影響；Kevin等[5]研制出1種直升機計算機程序來提醒飛行員避免擦機尾。國內(nèi)學(xué)者則對擦機尾的定性原因分析和維修研究較多：郭煥軍等[6]分析了空客飛機起飛和著陸2種情況下擦機尾的原因；舒平等[7]對飛機起飛階段擦機尾事故進(jìn)行分析并提出預(yù)防措施；毛吉星[8]就波音737飛機機尾擦地的損傷進(jìn)行級別劃分并制定對應(yīng)的維修方案；孫瑞山等[9]利用K-S檢驗?zāi)M飛機起飛擦機尾離地仰角值的可能分布，并計算具體機隊擦機尾事件的風(fēng)險概率，但選擇的樣本較少，計算的擦機尾風(fēng)險的準(zhǔn)確性有待提高。

航空公司一般根據(jù)飛機著陸時的俯仰角大小來判斷飛機是否擦機尾，當(dāng)著陸俯仰角達(dá)到某個特定值時，就會被判定為超限。為了深入挖掘接近擦機尾超限事件的俯仰角數(shù)據(jù)，利用蒙特卡洛模擬方法，突破樣本數(shù)據(jù)不足的障礙，對飛機著陸時擦機尾的風(fēng)險進(jìn)行準(zhǔn)確預(yù)測。

1 蒙特卡洛方法與飛機著陸擦機尾事件

蒙特卡洛方法[10]是1種通過生成合適的隨機數(shù)并在已知概率分布中進(jìn)行抽樣模擬試驗的統(tǒng)計方法，適用于所有隨機事件的模擬計算，在科學(xué)技術(shù)、經(jīng)濟管理、石油、金融等領(lǐng)域都有廣泛的應(yīng)用，其優(yōu)點在于能準(zhǔn)確地反映不確定性因素的影響，通過對有限的樣本數(shù)據(jù)進(jìn)行擴充實現(xiàn)對已知問題的求解。蒙特卡洛模擬分析的一般步驟如下[11-14]：

1)建立需要分析概率統(tǒng)計的原始數(shù)據(jù)模型；

2)產(chǎn)生符合給定概率分布的隨機值；

3)運用隨機抽樣的數(shù)據(jù)計算具體的估計量；

4)統(tǒng)計分析并繪制相應(yīng)的圖表。

在整個著陸階段，造成飛機著陸擦機尾的原因有：過早收油門；對油門和速度監(jiān)控不夠；接近地面時下降率太大；高度判斷不準(zhǔn)確，拉平過高；為了避免重著陸而平飄長；復(fù)飛時機晚，加油門不果斷。

以Boeing737-800(以下簡稱B737-800)飛機的著陸階段為例, 飛機從15.24 m(50 ft)高度以正常的運動軌跡、速度和正確的配平飛行，在主起落架高于跑道6.1 m(20 ft)時，逐漸增加俯仰姿態(tài)約2°～3°后開始拉平，減小飛機的下降率。開始拉平后，柔和地將油門桿收至慢車位，使飛機以正常的速度、方向和姿態(tài)在跑道入口前304.7～609.6 m(1 000～2 000 ft)處平穩(wěn)接地。

快速存取記錄器 (Quick Access Recorder，QAR)數(shù)據(jù)是航空公司有針對性地對飛行員進(jìn)行飛行訓(xùn)練并提高航班運行品質(zhì)和安全水平的重要依據(jù)。對QAR數(shù)據(jù)超限事件進(jìn)行統(tǒng)計、分析并據(jù)此開展訓(xùn)練已經(jīng)成為航空公司主要的飛行安全保障手段。

B737-800機型的QAR設(shè)備1s內(nèi)可連續(xù)記錄4次飛機著陸俯仰角的度數(shù)，如表1所示。

目前，國內(nèi)已有學(xué)者利用QAR數(shù)據(jù)對飛機重著陸風(fēng)險進(jìn)行評估預(yù)測，但其采集的航班數(shù)據(jù)有限，評估重著陸風(fēng)險有一定的局限性。因此，本文根據(jù)某航空公司飛行品質(zhì)監(jiān)控標(biāo)準(zhǔn)，將6.1 m(20 ft)離地高度至前輪接地所記錄的QAR數(shù)據(jù)作為樣本數(shù)據(jù)，并結(jié)合蒙特卡洛模擬方法，解決樣本不足的問題，建立擦機尾風(fēng)險預(yù)測模型，為航班擦機尾風(fēng)險預(yù)測提供參考。

表1 QAR原始數(shù)據(jù)示例Table 1 Example of QAR original data

2 基于蒙特卡洛模擬方法的擦機尾風(fēng)險預(yù)測模型

2.1 建立擦機尾風(fēng)險評價模型

風(fēng)險一般指某種特定的危險事件(事故或意外事件)發(fā)生的可能性與其產(chǎn)生的后果的組合。通過風(fēng)險的定義可以看出，風(fēng)險是由2個因素共同作用，一是危險發(fā)生的可能性，即危險概率；二是危險事件發(fā)生后所產(chǎn)生的后果。簡單來說，風(fēng)險可表示為某種不利事件或損失發(fā)生的概率及其后果的函數(shù)，即R=(P,S)。其中：R為風(fēng)險；P為不利事件或損失發(fā)生的概率；S為不利事件或損失發(fā)生的可能性。飛機擦機尾風(fēng)險評價模型可表示為擦機尾發(fā)生的可能性與擦機尾發(fā)生后導(dǎo)致的嚴(yán)重度的乘積，即：

RPA=PPASPA

(1)

式中：RPA為飛機擦機尾發(fā)生的風(fēng)險值；PPA為擦機尾發(fā)生的可能性；SPA為擦機尾事件的嚴(yán)重度。

2.2 隨機變量的抽樣模擬

2.2.1 正態(tài)分布隨機數(shù)的產(chǎn)生

研究[15]中對飛行QAR數(shù)據(jù)的大樣本統(tǒng)計發(fā)現(xiàn)：大部分飛行性能參數(shù)，如接地垂直載荷、著陸接地距離、接地仰角等，在較大的樣本空間(n>100)里都近似呈現(xiàn)正態(tài)分布?；诖耍妹商乜迥M方法生成服從正態(tài)分布的著陸俯仰角隨機度數(shù)，并依據(jù)相應(yīng)的算法預(yù)測擦機尾的風(fēng)險值。用蒙特卡洛模擬方法生成的隨機數(shù)序列在本質(zhì)上不是隨機的，只是很好地模仿了隨機數(shù)的性質(zhì)，但是如果隨機數(shù)滿足分布均勻性、隨機性和獨立性，并且能滿足實際問題的需要，此偽隨機數(shù)可以同真正的隨機數(shù)一樣使用。

2.2.2 抽樣模擬

抽樣是在某個隨機變量已知的情況下，通過取隨機數(shù)，實現(xiàn)該變量不斷取值的過程。滿足正態(tài)分布的隨機數(shù)生成算法很多，Box-Muller變換算法相對便捷、有效，因此，在獲得均勻分布隨機數(shù)的基礎(chǔ)上進(jìn)行Box-Muller變換, 公式如下：

(2)

式中：α為飛機著陸俯仰角的隨機度數(shù)，(°)；μ和σ分別為已知正態(tài)分布樣本數(shù)據(jù)的均值和標(biāo)準(zhǔn)差；r和t為在(0,1)值域內(nèi)服從均勻分布的獨立隨機數(shù)。

2.3 樣本數(shù)據(jù)及生成的隨機數(shù)的正態(tài)分布檢驗

1)計算均值和方差：

(3)

(4)

(5)

(6)

2)假設(shè)該樣本數(shù)據(jù)符合正態(tài)分布，則對其進(jìn)行正態(tài)檢驗，若結(jié)果表明樣本數(shù)據(jù)符合正態(tài)分布，再對其進(jìn)行方差修正得到總體方差進(jìn)行無偏估計：

(7)

(8)

3)由此可得該機型著陸時離地俯仰角的分布密度函數(shù)：

(9)

使用以期望值為對稱軸、α為矩形的區(qū)域概率PPA表示飛機擦機尾發(fā)生的可能性，如圖1所示 。

圖1 擦機尾參數(shù)的正態(tài)分布Fig.1 Normal distribution of tail strike risks

通過飛機著陸俯仰角與其分布期望值的關(guān)系來判斷擦機尾的風(fēng)險，α越靠近期望值說明擦機尾發(fā)生的可能性越小，運用蒙特卡洛模擬方法對飛機著陸俯仰角樣本進(jìn)行擴充，運用式(1)對擦機尾風(fēng)險值進(jìn)行具體化計算，可得較準(zhǔn)確的風(fēng)險預(yù)測模型。因此，PPA的計算公式可表示為：

(10)

式中：0≤PPA≤1，以樣本數(shù)據(jù)偏離均值的程度來表示飛機擦機尾事件的嚴(yán)重度，離均值越遠(yuǎn)則嚴(yán)重度越高，計算公式如下(SPA≥0)：

(11)

2.4 擦機尾風(fēng)險值的評估與預(yù)測

抽樣得到飛機著陸俯仰角在一定范圍內(nèi)的隨機度數(shù)α1,以風(fēng)險評價模型RPA=PPASPA為例,式中參數(shù)PPA和SPA的分布如式(10)、(11)所示，則對于每個預(yù)測擦機尾風(fēng)險值的估計值可由式(12)得出。飛機擦機尾風(fēng)險預(yù)測值的表達(dá)式為：

(12)

式中：r和t為服從(0，1)分布中相互獨立的隨機數(shù)；μ和σ分別為已知機隊樣本數(shù)據(jù)的均值和方差。

2.5 蒙特卡洛模擬流程

在上述預(yù)測模型基礎(chǔ)上，運用MATLAB編程進(jìn)行抽樣模擬試驗，系統(tǒng)運行5 000次，可得到5 000個擦機尾風(fēng)險預(yù)測值。蒙特卡洛模擬步驟如下：

分別取(0,1)分布上2個互相獨立的隨機數(shù)r和t，計算得到著陸俯仰角α1；將α1代入飛機擦機尾發(fā)生的可能性與后果的嚴(yán)重性計算公式中，得到PPA1和SPA1；將2個值相乘得到擦機尾風(fēng)險預(yù)測值，完成了1輪抽樣。依次進(jìn)行，如進(jìn)行5 000次抽樣，便有5 000個擦機尾風(fēng)險值R1，R2，…，R5 000。以著陸的隨機俯仰角度數(shù)為橫坐標(biāo)，預(yù)測的風(fēng)險值為縱坐標(biāo)，進(jìn)行統(tǒng)計整理得到擦機尾風(fēng)險的分布。模擬流程如圖2所示。

圖2 蒙特卡洛模擬流程示意Fig.2 The flow diagram of Monte Carlo

3 實例分析及模型驗證

本研究將某航空公司B737-800機隊的380套QAR數(shù)據(jù)作為樣本數(shù)據(jù)，其中采集機隊1共150個航班，機隊2共230個航班，運用MATLAB編程處理共380個航班的數(shù)據(jù)，得出2個不同機隊著陸時俯仰角的頻率分布直方圖，見圖3～4。

圖3 機隊1著陸俯仰角值的頻率分布Fig.3 Histogram of the landing pitch attitude distribution of the first fleet group

圖4 機隊2著陸俯仰角值的頻率分布Fig.4 Histogram of the landing pitch attitude distribution of the second fleet group

由圖3～4可知，機隊1著陸俯仰角均值為4.340 2，機隊2著陸俯仰角均值為4.317 7，整個機隊著陸俯仰角均值為4.377 7；機隊1著陸俯仰角方差為0.335 6，機隊2著陸俯仰角方差為0.332 6，整個機隊著陸俯仰角方差為0.342 8。

利用Box-Muller變換隨機生成5 000個著陸俯仰角度數(shù)，得到2個機隊著陸俯仰角隨機度數(shù)的分布直方圖，見圖5～6。

圖5 機隊1隨機生成的著陸俯仰角度數(shù)分布Fig.5 The random-generated landing pitch attitude distribution of the first fleet group

圖6 機隊2隨機生成的著陸俯仰角度數(shù)分布Fig.6 The random-generated landing pitch attitude distribution of the second fleet group

著陸俯仰角數(shù)據(jù)抽樣模擬之后，標(biāo)準(zhǔn)差數(shù)值由0.576 7變?yōu)?.558 1，說明蒙特卡洛模擬后的數(shù)據(jù)更加集中，驗證了模擬后的數(shù)據(jù)質(zhì)量得到進(jìn)一步提高，增加了風(fēng)險預(yù)測的可靠性。

取機隊2著陸俯仰角的實際數(shù)據(jù)與每組擬合數(shù)據(jù)值的差值，取絕對值計算相對誤差(相對誤差=絕對誤差÷實際值)為4.512 3，相對誤差的平均值為0.265 4，同理可得蒙特卡洛模擬后相對誤差的平均值為0.111 8，說明了蒙特卡洛模擬后數(shù)據(jù)更準(zhǔn)確、穩(wěn)定。

圖7 不同機隊擦機尾風(fēng)險預(yù)測曲線Fig.7 The tail strike prediction curve of different fleet groups

對樣本進(jìn)行Kolmogorov-Smirno(K-S)檢驗，結(jié)果符合正態(tài)分布。利用擦機尾風(fēng)險公式，計算隨機生成的每個著陸俯仰角對應(yīng)的風(fēng)險值，以著陸俯仰角為橫坐標(biāo)，擦機尾風(fēng)險值為縱坐標(biāo),得出不同機隊飛機擦機尾的風(fēng)險預(yù)測曲線，見圖7。圖7是不同機隊擦機尾風(fēng)險進(jìn)行模擬5 000次后的結(jié)果(著陸俯仰角為負(fù)值時，風(fēng)險值記為0)，從圖中可以看出，當(dāng)飛機著陸俯仰角相同時，兩機隊擦機尾的風(fēng)險值不同：當(dāng)著陸俯仰角度數(shù)小于4.5°時，兩機隊的擦機尾風(fēng)險值無太大區(qū)別，但當(dāng)著陸俯仰角大于4.5°時，機隊2擦機尾風(fēng)險明顯大于機隊1，如當(dāng)飛機著陸俯仰角為5°時，機隊1擦機尾風(fēng)險值為1.70，機隊2擦機尾風(fēng)險值為2.03，航空公司可據(jù)此加強對機隊2飛行員的訓(xùn)練,降低擦機尾的風(fēng)險。

蒙特卡洛模擬指標(biāo)匯總?cè)绫?所示。通過蒙特卡洛模擬得出預(yù)測結(jié)果的最大值、最小值和最可能值，給出了擦機尾風(fēng)險值的預(yù)測區(qū)間范圍及分布規(guī)律，證明此模擬方法能較為合理地預(yù)測風(fēng)險。

表2 蒙特卡洛模擬指標(biāo)匯總Table 2 Indicator summary of Monte Carlo

注：以上數(shù)據(jù)均是通過MATLAB模擬顯示的最終結(jié)果，其中著陸俯仰角的模擬次數(shù)為5 000次。

4 結(jié)論

1)本文利用原始QAR數(shù)據(jù)作為樣本，基于蒙特卡洛方法對飛機著陸時的擦機尾風(fēng)險值進(jìn)行預(yù)測，進(jìn)行模擬抽樣試驗5 000次，較好地解決獲取航班信息有限的問題。

2)利用蒙特卡洛算法進(jìn)行模擬，能夠比較準(zhǔn)確地估算出整個機隊在某一運行周期內(nèi)起飛或著陸階段擦機尾的風(fēng)險值，適用于單個航班擦機尾風(fēng)險值的預(yù)測，也適用于整個機隊在一段時期內(nèi)擦機尾風(fēng)險預(yù)測，還可延伸到其他超限事件的風(fēng)險評估，為航空公司的超限事件管理提供可靠性指標(biāo)參考，實現(xiàn)對機隊操作風(fēng)險的動態(tài)管理。

3)該風(fēng)險預(yù)測方法雖模擬出飛機著陸時擦機尾的風(fēng)險值，但是由于飛機在著陸時受到諸多因素影響，各因素對風(fēng)險產(chǎn)生的影響需進(jìn)一步研究。