何光勤,魯 力,胡敬玉,馬 昕
(中國民用航空飛行學(xué)院空中交通管理學(xué)院,四川 廣漢 618300)
飛機在雷暴天氣中飛行是十分危險的,不僅會產(chǎn)生嚴重顛簸影響旅客的舒適度,還有可能受到雷擊以及積冰和冰雹的影響,對飛機電子設(shè)備和結(jié)構(gòu)安全造成巨大的威脅[1]。1994年2月4日,美國一家航空公司的DC-9-31飛機在飛行中進入雷暴區(qū)時,兩臺發(fā)動機吸入大量水和冰雹,發(fā)動機壓縮器損傷導(dǎo)致失速墜毀。因此,預(yù)測雷暴的位置,有效避免飛機飛行中受到雷暴天氣的影響有著重要意義。早在1842年,奧地利物理學(xué)家J.Doppler從運動聲源受到啟發(fā)得出多普勒效應(yīng),根據(jù)此效應(yīng)研發(fā)的多普勒氣象雷達對民航業(yè)的運行安全發(fā)揮了重要作用,它通過發(fā)射高頻脈沖電磁波信號探測搜索范圍內(nèi)的氣象情況,比如湍流、雷暴等天氣,再根據(jù)反射信號的頻差、相位差、衰減等特性分析前方是何種天氣狀況,并用不同的顏色反映到雷達回波圖上[2]。近些年,國內(nèi)很多學(xué)者根據(jù)氣象雷達回波圖,采用最小二乘法進行了雷暴走勢預(yù)測,如陳嵐峰等[3]采用Matlab軟件對曲線進行了最小二乘法擬合仿真研究;梁平等[4]運用最小二乘法和地理加權(quán)模型對漢江流域土地利用類型與水質(zhì)關(guān)系進行了評估;肖云等[5]利用空域最小二乘法分析了重力衛(wèi)星誤差;陳曉倩等[6]采用最小二乘法預(yù)測了無人機路徑。對于航空器碰撞風(fēng)險,國內(nèi)一些學(xué)者對此也做了一些研究,如呂宗平等[7]采用CREAM和IDAC模型對航空器碰撞風(fēng)險進行了評估;薛宇敬陽等[8]采用飛行路徑算法預(yù)測判斷了飛機碰撞風(fēng)險。另外,還有一些學(xué)者采用各種算法模型對民航事故征候進行了分析和預(yù)測[9-14]。
基于上述研究,本文基于氣象雷達回波圖分析,采用最小二乘法進行了雷暴走勢預(yù)測,并判斷飛機正常飛行是否會受到雷暴天氣的影響,若受影響飛機則需要改變原飛行航跡,以避免穿越雷暴區(qū)而引發(fā)的嚴重事故。
本文利用兩運動物體相撞模型,判斷飛機沿預(yù)計飛行航跡運動是否會在未來某時刻受雷暴天氣的影響。對于飛機的運動,知道運動速度和航向,可計算飛機運動的軌跡方程。但對于雷暴的運動,運動方向和速度在不斷的變化,需借助軟件找到雷暴運動的軌跡方程。
物體運動軌跡在二維坐標系中可分解到x軸和y軸上,用x(t)和y(t)表示物體相對兩坐標軸關(guān)于時間t的運動方程。飛機沿預(yù)計飛行航跡的運動速度和航向是已知的,可找到飛機運動的軌跡方程,記為x0(t)和y0(t)。
在氣象雷達回波圖中,雷暴顏色越濃(見圖1),雷暴等級越大。CAD軟件具有分析圖像以及準確、方便地找出坐標數(shù)據(jù)的功能。因此,利用氣象雷達回波圖得到的某時刻雷暴的坐標數(shù)據(jù)[15-16],并采用最小二乘法可擬合得到雷暴運動的軌跡方程x′(t)和y′(t)。
圖1 氣象雷達回波圖Fig.1 Meteorological radar echo map
(1)
其中,x′(ti)可以是ti的一次、二次或三次方程,由Q對參數(shù)求導(dǎo)等于0即可得到Q的最小值,則x′(ti)則可由參數(shù)a、b、c表達。
本文選取多個方程,但找出最符合雷暴運動的軌跡方程則需要進行方差分析。假設(shè)在未來t時刻已知雷暴的坐標數(shù)據(jù)為(x0,y0),再利用Matlab軟件擬合得到的雷暴運動軌跡方程計算在t時刻雷暴的坐標x′(t)和y′(t),其方差計算公式如下:
D=(x′(t)-x0)2+(y′(t)-y0)2
(2)
找出方差最小的方程即為最匹配雷暴運動的軌跡方程。
將飛機和雷暴當成兩運動物體,由于兩物體的運動軌跡均可視為沿x、y軸運動相對于時間t的函數(shù)方程,則兩運動物體之間的距離可寫成關(guān)于時間t的函數(shù),即:
(3)
S對t求導(dǎo)可得極值,由此可得出在未來何時兩物體間的距離最小,并將其最小距離與飛機運行標準相比,若該最小距離小于飛機運行標準的要求,則飛機需要提前改變飛行航跡以避免飛機與雷暴相遇。
本文以2018年6月9日首都航空JD5690次航班經(jīng)由宜昌三峽機場飛往杭州蕭山機場的飛行過程為例,根據(jù)氣象雷達回波圖可知當日在此航班飛行航路上有雷暴天氣,利用兩物體運動相撞模型,評估飛機在終端區(qū)若沿預(yù)計飛行航跡運動是否會受雷暴天氣的影響。本文基于氣象雷達回波圖分析,并采用最小二乘法進行雷暴走勢預(yù)測,該方法的具體分析過程如下:
第一步:找出飛機運動的軌跡方程。
首先,根據(jù)飛機的運動速度和航向,以及氣象雷達回波圖并結(jié)合Auto CAD軟件,可獲取JD5690次航班不同時刻在氣象雷達回波圖上的坐標位置,見圖2。
圖2 JD5690次航班不同時刻在氣象雷達回波圖上的 坐標位置Fig.2 Coordinate position of flight JD5690 on the meteorological radar echo wave at different times
本節(jié)給出的飛機運動軌跡坐標值的比例尺均為1∶25 000,單位為m。
然后,根據(jù)飛機在17∶00 pm和17∶04 pm時刻的坐標位置分別為(0.821 7,5.733 8)、(2.036 3,6.401 8),可得到飛機運動的軌跡方程為
(4)
最后,使用Matlab軟件作出方程(4)的曲線,即飛機沿預(yù)計飛行航線運動的軌跡圖,見圖3。
圖3 飛機沿預(yù)計飛行航線運動的軌跡圖Fig.3 Trajectory map of the aircraft along the expected track movement
第二步:找出雷暴運動的軌跡方程。
雷暴的運動速度和方向都是時刻變化的,首先根據(jù)氣象雷達回波圖并使用Auto CAD軟件,找出已知時刻雷暴中心的坐標位置,見表1。
表1 已知時刻雷暴中心的坐標位置數(shù)據(jù)列表
然后,利用最小二乘法并使用Matlab軟件,找出符合已知條件的雷暴運動軌跡擬合方程如下:
(5)
(6)
(7)
第三步:根據(jù)方差找出最匹配雷暴運動軌跡的擬合方程。
采取16∶31pm時刻實際雷暴中心的坐標位置(4.345,7.751 8)和第二步中三個擬合方程[公式(5)、(6)、(7)]預(yù)測得到的在該時刻雷暴中心坐標位置的方差進行分析,并判斷哪個方程的擬合效果最佳[17]。設(shè)三個擬合方程的方差計算結(jié)果為D1、D2、D3,則三個擬合方程的方差值見表2。
表2 三個擬合方程的方差值列表
由表2可知,擬合方程(5)的方差值最小,故選取該方程作為雷暴運動的軌跡方程,并使用Matlab軟件作出擬合方程(5)的曲線,即雷暴運動的軌跡圖,見圖4。
圖4 雷暴運動的軌跡圖Fig.4 Motion trajectory map of thunderstorms
第四步:在規(guī)定的時間內(nèi)計算飛機與雷暴之間的最小距離。
由兩運動物體之間距離S(t)的計算公式(3),可計算未來20 min飛機與雷暴之間的最小距離。將飛機和雷暴運動的軌跡使用Matlab軟件繪制到一張圖中,由此可以直觀地看出何時兩者之間的距離最小,見圖5。
由圖5可見,大約在飛機運動軌跡的第8個軌跡點即約16 min時刻,飛機與雷暴之間的距離最小。
利用公式(3)具體求解過程中,應(yīng)提前將飛機與雷暴運動的軌跡方程式在時間上相對應(yīng),再采用S(t)對t進行求導(dǎo),即可求得在t0時的極小值S(t0),此步驟可在Matlab軟件中編程求解,具體編程代碼如下:
symst
x=@(t)(1.285 1*(10^(-4))*(t+113)^2+0.014 5*(t+113)+2.115 0-0.303 7*t-0.821 7)^2+((-3.828 4)*(10^(-5))*(t+113)^2+0.011 6*(t+113)+6.803 9-0.167*t-5.733 8)^2;<此代碼表示飛機與雷暴之間的間隔>
[tmin]=fminbnd(x,0,20)
tmin=
15.900 4
[xtmin]=double(subs(x,t,[tmin]))
Matlab軟件運行得出結(jié)果如下:
tmin
15.9004
xtmin=
0.746 9
圖5 飛機和雷暴運動的軌跡圖Fig.5 Motion trajectory map of the aircraft and thunderstorms注:點與點的時間間隔為2 min。
本文運用數(shù)理統(tǒng)計學(xué)中的區(qū)間估計方法來驗證該方法的可行性。首先利用中心極限定理可近視認為S(t0)服從N(μ,σ2)的正態(tài)分布,則μ=21,σ2=D(x)=0.070 4,而根據(jù)切比雪夫不等式:
(8)
然后根據(jù)飛機運行標準可知,飛機與雷暴之間的最小距離應(yīng)當不低于25 km,而本文計算得到的飛機與雷暴之間的距離最小值為21 km,那么可設(shè)ε=4,使得在t0=15.9 min時刻飛機與雷暴之間的距離不大于25 km,計算得到其概率P≥0.996 5,即該事件為大概率事件,表明此方法是可行的。
相應(yīng)地,t0為服從于N(μ,σ2)的正態(tài)分布,則μ=15.9,σ2=D(x)=0.070 4,在知道該時刻飛機與雷暴之間有最小距離21 km的前提下,給定置信水平為0.9,估計t0的置信區(qū)間,根據(jù)該置信區(qū)間可找出飛機需要改變飛行航跡的最早時間點,其置信區(qū)間為
(9)
由于預(yù)計雷暴運動軌跡方程使用的采樣點有4個,故n=4,zα/2為標準正態(tài)分布的上分位點(見圖6),經(jīng)查表計算可得該置信區(qū)間為(15.842,15.958),則可視為飛機應(yīng)在15.8 min時刻之前改變飛行航跡,以避免雷暴天氣的影響[18-20]。
圖6 標準正態(tài)分布上α分位點示意圖Fig.6 Diagram of quantil on standard normal
根據(jù)第2.2節(jié)的區(qū)間估計方法,JD5690次航班已在10 min后改變了飛行航跡(見圖7),這與本文利用最小二乘法推導(dǎo)出的結(jié)果基本一致。
圖7 10 min后JD5690次航班改航圖Fig.7 JD5690 flight diversion chart after 10min
2018年5月20號CA1860次航班由柳州飛往北京,本沿預(yù)計航路SJG飛往P59航路點,中途經(jīng)過雷暴區(qū)(見圖8,雷暴區(qū)顏色為紅色),雷暴移動方向由西北方向向P59航路點移動。本文利用最小二乘法對雷暴運動的軌跡進行了預(yù)測,其預(yù)測結(jié)果如下。
雷暴中心在0 min
(10)
飛機在0 min
(11)
根據(jù)第2.2節(jié)的區(qū)間估計方法,可得出飛機預(yù)計在(1.982,2.342)區(qū)間與雷暴發(fā)生碰撞,因此飛機必須通過改變飛行航跡實行避讓雷暴。如圖8中的主飛行儀表(Primary Flight Display,PFD)所示,CA1860次航班的航向已改為280°,向西飛行,以避讓雷暴,從而驗證了本文提出的最小二乘法預(yù)測雷暴走勢的準確性。
圖8 CA1860次航班途經(jīng)雷暴區(qū)Fig.8 CA1860 flight passing through the thunderstorm area
本文先基于氣象雷達回波圖并結(jié)合Auto CAD進行圖像分析提取坐標位置數(shù)據(jù),再采用最小二乘法進行線性擬合獲得飛機和雷暴運動的軌跡方程,并通過兩者運動軌跡方程的聯(lián)立計算并判斷飛機是否會進入雷暴區(qū)且受其影響,最后通過案例驗證表明,實際結(jié)果處于最小二乘法預(yù)測的范圍之內(nèi),說明該方法具有一定的準確性。本文在對雷暴運動的軌跡方程進行選取時運用了方差公式和切比雪夫不等式等數(shù)據(jù)統(tǒng)計學(xué)方法進行可行性分析,并采用置信區(qū)間評估飛機與雷暴之間最小距離時刻的時間容差,使該方法具有實際預(yù)測意義。但本研究在選取雷暴中心的坐標位置數(shù)據(jù)時有一定誤差,且飛機在實際飛行中運動速度和航向也會發(fā)生變化,在以后的研究中,應(yīng)對這些變化進行處理,使預(yù)測結(jié)果更加精確。