基于蒙特卡洛方法的著陸擦機(jī)尾風(fēng)險(xiǎn)預(yù)測(cè)*

汪 磊，楊星月

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 基于蒙特卡洛模擬方法的擦機(jī)尾風(fēng)險(xiǎn)預(yù)測(cè)模型

2.1 建立擦機(jī)尾風(fēng)險(xiǎn)評(píng)價(jià)模型

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

RPA=PPASPA

(1)

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

2.2 隨機(jī)變量的抽樣模擬

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

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

2.2.2 抽樣模擬

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

(2)

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

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

1)計(jì)算均值和方差：

(3)

(4)

(5)

(6)

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

(7)

(8)

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

(9)

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

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

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

(10)

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

(11)

2.4 擦機(jī)尾風(fēng)險(xiǎn)值的評(píng)估與預(yù)測(cè)

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

(12)

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

2.5 蒙特卡洛模擬流程

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

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

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

3 實(shí)例分析及模型驗(yàn)證

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

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

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

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

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

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

圖6 機(jī)隊(duì)2隨機(jī)生成的著陸俯仰角度數(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，說(shuō)明蒙特卡洛模擬后的數(shù)據(jù)更加集中，驗(yàn)證了模擬后的數(shù)據(jù)質(zhì)量得到進(jìn)一步提高，增加了風(fēng)險(xiǎn)預(yù)測(cè)的可靠性。

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

圖7 不同機(jī)隊(duì)擦機(jī)尾風(fēng)險(xiǎn)預(yù)測(cè)曲線Fig.7 The tail strike prediction curve of different fleet groups

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

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

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

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

4 結(jié)論

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

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

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