基于最小二乘法的雷暴天氣下飛行改航?jīng)Q策研究

何光勤，魯 力，胡敬玉,馬 昕

(中國(guó)民用航空飛行學(xué)院空中交通管理學(xué)院，四川 廣漢 618300)

飛機(jī)在雷暴天氣中飛行是十分危險(xiǎn)的，不僅會(huì)產(chǎn)生嚴(yán)重顛簸影響旅客的舒適度，還有可能受到雷擊以及積冰和冰雹的影響，對(duì)飛機(jī)電子設(shè)備和結(jié)構(gòu)安全造成巨大的威脅[1]。1994年2月4日，美國(guó)一家航空公司的DC-9-31飛機(jī)在飛行中進(jìn)入雷暴區(qū)時(shí)，兩臺(tái)發(fā)動(dòng)機(jī)吸入大量水和冰雹，發(fā)動(dòng)機(jī)壓縮器損傷導(dǎo)致失速墜毀。因此，預(yù)測(cè)雷暴的位置，有效避免飛機(jī)飛行中受到雷暴天氣的影響有著重要意義。早在1842年，奧地利物理學(xué)家J.Doppler從運(yùn)動(dòng)聲源受到啟發(fā)得出多普勒效應(yīng)，根據(jù)此效應(yīng)研發(fā)的多普勒氣象雷達(dá)對(duì)民航業(yè)的運(yùn)行安全發(fā)揮了重要作用，它通過(guò)發(fā)射高頻脈沖電磁波信號(hào)探測(cè)搜索范圍內(nèi)的氣象情況，比如湍流、雷暴等天氣，再根據(jù)反射信號(hào)的頻差、相位差、衰減等特性分析前方是何種天氣狀況，并用不同的顏色反映到雷達(dá)回波圖上[2]。近些年，國(guó)內(nèi)很多學(xué)者根據(jù)氣象雷達(dá)回波圖，采用最小二乘法進(jìn)行了雷暴走勢(shì)預(yù)測(cè)，如陳嵐峰等[3]采用Matlab軟件對(duì)曲線進(jìn)行了最小二乘法擬合仿真研究；梁平等[4]運(yùn)用最小二乘法和地理加權(quán)模型對(duì)漢江流域土地利用類型與水質(zhì)關(guān)系進(jìn)行了評(píng)估；肖云等[5]利用空域最小二乘法分析了重力衛(wèi)星誤差；陳曉倩等[6]采用最小二乘法預(yù)測(cè)了無(wú)人機(jī)路徑。對(duì)于航空器碰撞風(fēng)險(xiǎn)，國(guó)內(nèi)一些學(xué)者對(duì)此也做了一些研究，如呂宗平等[7]采用CREAM和IDAC模型對(duì)航空器碰撞風(fēng)險(xiǎn)進(jìn)行了評(píng)估；薛宇敬陽(yáng)等[8]采用飛行路徑算法預(yù)測(cè)判斷了飛機(jī)碰撞風(fēng)險(xiǎn)。另外，還有一些學(xué)者采用各種算法模型對(duì)民航事故征候進(jìn)行了分析和預(yù)測(cè)[9-14]。

基于上述研究，本文基于氣象雷達(dá)回波圖分析，采用最小二乘法進(jìn)行了雷暴走勢(shì)預(yù)測(cè)，并判斷飛機(jī)正常飛行是否會(huì)受到雷暴天氣的影響，若受影響飛機(jī)則需要改變?cè)w行航跡，以避免穿越雷暴區(qū)而引發(fā)的嚴(yán)重事故。

1 理論與算法

本文利用兩運(yùn)動(dòng)物體相撞模型，判斷飛機(jī)沿預(yù)計(jì)飛行航跡運(yùn)動(dòng)是否會(huì)在未來(lái)某時(shí)刻受雷暴天氣的影響。對(duì)于飛機(jī)的運(yùn)動(dòng)，知道運(yùn)動(dòng)速度和航向，可計(jì)算飛機(jī)運(yùn)動(dòng)的軌跡方程。但對(duì)于雷暴的運(yùn)動(dòng)，運(yùn)動(dòng)方向和速度在不斷的變化，需借助軟件找到雷暴運(yùn)動(dòng)的軌跡方程。

1. 1 飛機(jī)運(yùn)動(dòng)的軌跡方程

物體運(yùn)動(dòng)軌跡在二維坐標(biāo)系中可分解到x軸和y軸上，用x(t)和y(t)表示物體相對(duì)兩坐標(biāo)軸關(guān)于時(shí)間t的運(yùn)動(dòng)方程。飛機(jī)沿預(yù)計(jì)飛行航跡的運(yùn)動(dòng)速度和航向是已知的，可找到飛機(jī)運(yùn)動(dòng)的軌跡方程,記為x0(t)和y0(t)。

1. 2 采用最小二乘法擬合雷暴運(yùn)動(dòng)的軌跡方程

在氣象雷達(dá)回波圖中，雷暴顏色越濃(見圖1)，雷暴等級(jí)越大。CAD軟件具有分析圖像以及準(zhǔn)確、方便地找出坐標(biāo)數(shù)據(jù)的功能。因此，利用氣象雷達(dá)回波圖得到的某時(shí)刻雷暴的坐標(biāo)數(shù)據(jù)[15-16]，并采用最小二乘法可擬合得到雷暴運(yùn)動(dòng)的軌跡方程x′(t)和y′(t)。

圖1 氣象雷達(dá)回波圖Fig.1 Meteorological radar echo map

(1)

其中，x′(ti)可以是ti的一次、二次或三次方程，由Q對(duì)參數(shù)求導(dǎo)等于0即可得到Q的最小值，則x′(ti)則可由參數(shù)a、b、c表達(dá)。

本文選取多個(gè)方程，但找出最符合雷暴運(yùn)動(dòng)的軌跡方程則需要進(jìn)行方差分析。假設(shè)在未來(lái)t時(shí)刻已知雷暴的坐標(biāo)數(shù)據(jù)為(x0,y0)，再利用Matlab軟件擬合得到的雷暴運(yùn)動(dòng)軌跡方程計(jì)算在t時(shí)刻雷暴的坐標(biāo)x′(t)和y′(t)，其方差計(jì)算公式如下：

D=(x′(t)-x0)2+(y′(t)-y0)2

(2)

找出方差最小的方程即為最匹配雷暴運(yùn)動(dòng)的軌跡方程。

1. 3 兩運(yùn)動(dòng)物體相撞判斷

將飛機(jī)和雷暴當(dāng)成兩運(yùn)動(dòng)物體，由于兩物體的運(yùn)動(dòng)軌跡均可視為沿x、y軸運(yùn)動(dòng)相對(duì)于時(shí)間t的函數(shù)方程，則兩運(yùn)動(dòng)物體之間的距離可寫成關(guān)于時(shí)間t的函數(shù)，即：

(3)

S對(duì)t求導(dǎo)可得極值，由此可得出在未來(lái)何時(shí)兩物體間的距離最小，并將其最小距離與飛機(jī)運(yùn)行標(biāo)準(zhǔn)相比，若該最小距離小于飛機(jī)運(yùn)行標(biāo)準(zhǔn)的要求，則飛機(jī)需要提前改變飛行航跡以避免飛機(jī)與雷暴相遇。

2 方法與驗(yàn)證

2. 1 方法分析過(guò)程

本文以2018年6月9日首都航空J(rèn)D5690次航班經(jīng)由宜昌三峽機(jī)場(chǎng)飛往杭州蕭山機(jī)場(chǎng)的飛行過(guò)程為例，根據(jù)氣象雷達(dá)回波圖可知當(dāng)日在此航班飛行航路上有雷暴天氣，利用兩物體運(yùn)動(dòng)相撞模型，評(píng)估飛機(jī)在終端區(qū)若沿預(yù)計(jì)飛行航跡運(yùn)動(dòng)是否會(huì)受雷暴天氣的影響。本文基于氣象雷達(dá)回波圖分析，并采用最小二乘法進(jìn)行雷暴走勢(shì)預(yù)測(cè)，該方法的具體分析過(guò)程如下：

第一步：找出飛機(jī)運(yùn)動(dòng)的軌跡方程。

首先，根據(jù)飛機(jī)的運(yùn)動(dòng)速度和航向，以及氣象雷達(dá)回波圖并結(jié)合Auto CAD軟件，可獲取JD5690次航班不同時(shí)刻在氣象雷達(dá)回波圖上的坐標(biāo)位置，見圖2。

圖2 JD5690次航班不同時(shí)刻在氣象雷達(dá)回波圖上的 坐標(biāo)位置Fig.2 Coordinate position of flight JD5690 on the meteorological radar echo wave at different times

本節(jié)給出的飛機(jī)運(yùn)動(dòng)軌跡坐標(biāo)值的比例尺均為1∶25 000，單位為m。

然后,根據(jù)飛機(jī)在17∶00 pm和17∶04 pm時(shí)刻的坐標(biāo)位置分別為(0.821 7,5.733 8)、(2.036 3,6.401 8)，可得到飛機(jī)運(yùn)動(dòng)的軌跡方程為

(4)

最后，使用Matlab軟件作出方程(4)的曲線，即飛機(jī)沿預(yù)計(jì)飛行航線運(yùn)動(dòng)的軌跡圖，見圖3。

圖3 飛機(jī)沿預(yù)計(jì)飛行航線運(yùn)動(dòng)的軌跡圖Fig.3 Trajectory map of the aircraft along the expected track movement

第二步：找出雷暴運(yùn)動(dòng)的軌跡方程。

雷暴的運(yùn)動(dòng)速度和方向都是時(shí)刻變化的，首先根據(jù)氣象雷達(dá)回波圖并使用Auto CAD軟件，找出已知時(shí)刻雷暴中心的坐標(biāo)位置，見表1。

表1 已知時(shí)刻雷暴中心的坐標(biāo)位置數(shù)據(jù)列表

然后，利用最小二乘法并使用Matlab軟件，找出符合已知條件的雷暴運(yùn)動(dòng)軌跡擬合方程如下：

(5)

(6)

(7)

第三步：根據(jù)方差找出最匹配雷暴運(yùn)動(dòng)軌跡的擬合方程。

采取16∶31pm時(shí)刻實(shí)際雷暴中心的坐標(biāo)位置(4.345,7.751 8)和第二步中三個(gè)擬合方程[公式(5)、(6)、(7)]預(yù)測(cè)得到的在該時(shí)刻雷暴中心坐標(biāo)位置的方差進(jìn)行分析，并判斷哪個(gè)方程的擬合效果最佳[17]。設(shè)三個(gè)擬合方程的方差計(jì)算結(jié)果為D1、D2、D3，則三個(gè)擬合方程的方差值見表2。

表2 三個(gè)擬合方程的方差值列表

由表2可知，擬合方程(5)的方差值最小，故選取該方程作為雷暴運(yùn)動(dòng)的軌跡方程,并使用Matlab軟件作出擬合方程(5)的曲線，即雷暴運(yùn)動(dòng)的軌跡圖，見圖4。

圖4 雷暴運(yùn)動(dòng)的軌跡圖Fig.4 Motion trajectory map of thunderstorms

第四步：在規(guī)定的時(shí)間內(nèi)計(jì)算飛機(jī)與雷暴之間的最小距離。

由兩運(yùn)動(dòng)物體之間距離S(t)的計(jì)算公式(3),可計(jì)算未來(lái)20 min飛機(jī)與雷暴之間的最小距離。將飛機(jī)和雷暴運(yùn)動(dòng)的軌跡使用Matlab軟件繪制到一張圖中，由此可以直觀地看出何時(shí)兩者之間的距離最小，見圖5。

由圖5可見，大約在飛機(jī)運(yùn)動(dòng)軌跡的第8個(gè)軌跡點(diǎn)即約16 min時(shí)刻，飛機(jī)與雷暴之間的距離最小。

利用公式(3)具體求解過(guò)程中，應(yīng)提前將飛機(jī)與雷暴運(yùn)動(dòng)的軌跡方程式在時(shí)間上相對(duì)應(yīng)，再采用S(t)對(duì)t進(jìn)行求導(dǎo)，即可求得在t0時(shí)的極小值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;<此代碼表示飛機(jī)與雷暴之間的間隔>

[tmin]=fminbnd(x,0,20)

tmin=

15.900 4

[xtmin]=double(subs(x,t,[tmin]))

Matlab軟件運(yùn)行得出結(jié)果如下：

tmin

15.9004

xtmin=

0.746 9

圖5 飛機(jī)和雷暴運(yùn)動(dòng)的軌跡圖Fig.5 Motion trajectory map of the aircraft and thunderstorms注：點(diǎn)與點(diǎn)的時(shí)間間隔為2 min。

2. 2 利用數(shù)理統(tǒng)計(jì)學(xué)方法驗(yàn)證方法的可行性

本文運(yùn)用數(shù)理統(tǒng)計(jì)學(xué)中的區(qū)間估計(jì)方法來(lái)驗(yàn)證該方法的可行性。首先利用中心極限定理可近視認(rèn)為S(t0)服從N(μ,σ2)的正態(tài)分布，則μ=21，σ2=D(x)=0.070 4，而根據(jù)切比雪夫不等式：

(8)

然后根據(jù)飛機(jī)運(yùn)行標(biāo)準(zhǔn)可知，飛機(jī)與雷暴之間的最小距離應(yīng)當(dāng)不低于25 km，而本文計(jì)算得到的飛機(jī)與雷暴之間的距離最小值為21 km，那么可設(shè)ε=4，使得在t0=15.9 min時(shí)刻飛機(jī)與雷暴之間的距離不大于25 km，計(jì)算得到其概率P≥0.996 5，即該事件為大概率事件，表明此方法是可行的。

相應(yīng)地，t0為服從于N(μ,σ2)的正態(tài)分布，則μ=15.9，σ2=D(x)=0.070 4，在知道該時(shí)刻飛機(jī)與雷暴之間有最小距離21 km的前提下，給定置信水平為0.9，估計(jì)t0的置信區(qū)間，根據(jù)該置信區(qū)間可找出飛機(jī)需要改變飛行航跡的最早時(shí)間點(diǎn)，其置信區(qū)間為

(9)

由于預(yù)計(jì)雷暴運(yùn)動(dòng)軌跡方程使用的采樣點(diǎn)有4個(gè)，故n=4，zα/2為標(biāo)準(zhǔn)正態(tài)分布的上分位點(diǎn)(見圖6)，經(jīng)查表計(jì)算可得該置信區(qū)間為(15.842,15.958)，則可視為飛機(jī)應(yīng)在15.8 min時(shí)刻之前改變飛行航跡，以避免雷暴天氣的影響[18-20]。

圖6 標(biāo)準(zhǔn)正態(tài)分布上α分位點(diǎn)示意圖Fig.6 Diagram of quantil on standard normal

根據(jù)第2.2節(jié)的區(qū)間估計(jì)方法，JD5690次航班已在10 min后改變了飛行航跡(見圖7)，這與本文利用最小二乘法推導(dǎo)出的結(jié)果基本一致。

圖7 10 min后JD5690次航班改航圖Fig.7 JD5690 flight diversion chart after 10min

2. 3 案例驗(yàn)證

2018年5月20號(hào)CA1860次航班由柳州飛往北京，本沿預(yù)計(jì)航路SJG飛往P59航路點(diǎn)，中途經(jīng)過(guò)雷暴區(qū)(見圖8，雷暴區(qū)顏色為紅色)，雷暴移動(dòng)方向由西北方向向P59航路點(diǎn)移動(dòng)。本文利用最小二乘法對(duì)雷暴運(yùn)動(dòng)的軌跡進(jìn)行了預(yù)測(cè)，其預(yù)測(cè)結(jié)果如下。

雷暴中心在0 min

(10)

飛機(jī)在0 min

(11)

根據(jù)第2.2節(jié)的區(qū)間估計(jì)方法，可得出飛機(jī)預(yù)計(jì)在(1.982,2.342)區(qū)間與雷暴發(fā)生碰撞，因此飛機(jī)必須通過(guò)改變飛行航跡實(shí)行避讓雷暴。如圖8中的主飛行儀表(Primary Flight Display，PFD)所示，CA1860次航班的航向已改為280°，向西飛行，以避讓雷暴，從而驗(yàn)證了本文提出的最小二乘法預(yù)測(cè)雷暴走勢(shì)的準(zhǔn)確性。

圖8 CA1860次航班途經(jīng)雷暴區(qū)Fig.8 CA1860 flight passing through the thunderstorm area

3 結(jié) 論

本文先基于氣象雷達(dá)回波圖并結(jié)合Auto CAD進(jìn)行圖像分析提取坐標(biāo)位置數(shù)據(jù)，再采用最小二乘法進(jìn)行線性擬合獲得飛機(jī)和雷暴運(yùn)動(dòng)的軌跡方程，并通過(guò)兩者運(yùn)動(dòng)軌跡方程的聯(lián)立計(jì)算并判斷飛機(jī)是否會(huì)進(jìn)入雷暴區(qū)且受其影響，最后通過(guò)案例驗(yàn)證表明，實(shí)際結(jié)果處于最小二乘法預(yù)測(cè)的范圍之內(nèi)，說(shuō)明該方法具有一定的準(zhǔn)確性。本文在對(duì)雷暴運(yùn)動(dòng)的軌跡方程進(jìn)行選取時(shí)運(yùn)用了方差公式和切比雪夫不等式等數(shù)據(jù)統(tǒng)計(jì)學(xué)方法進(jìn)行可行性分析，并采用置信區(qū)間評(píng)估飛機(jī)與雷暴之間最小距離時(shí)刻的時(shí)間容差，使該方法具有實(shí)際預(yù)測(cè)意義。但本研究在選取雷暴中心的坐標(biāo)位置數(shù)據(jù)時(shí)有一定誤差，且飛機(jī)在實(shí)際飛行中運(yùn)動(dòng)速度和航向也會(huì)發(fā)生變化，在以后的研究中，應(yīng)對(duì)這些變化進(jìn)行處理，使預(yù)測(cè)結(jié)果更加精確。