基于雷達(dá)數(shù)據(jù)的飛行受限區(qū)形變及方位預(yù)測(cè)

樊立艷，趙鶴宇，常錦才,3，李印鳳

（1.華北理工大學(xué)理學(xué)院，河北 唐山 063210；2.華北理工大學(xué)建筑工程學(xué)院，河北 唐山 063210；3.華北理工大學(xué)河北省數(shù)據(jù)科學(xué)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，河北 唐山 063210）

隨著科技的發(fā)展和人民生活水平的提升，人們使用航空器的次數(shù)及貨物運(yùn)輸?shù)念l率也在日益增多，同時(shí)空域資源也越來越緊張，一旦出現(xiàn)不良?xì)庀蠡蚴侵卮笫录沼蛸Y源的利用率將會(huì)急劇下降，進(jìn)而造成飛機(jī)晚點(diǎn)、乘客滯留甚至航班取消，這不但會(huì)造成巨大的經(jīng)濟(jì)損失，也會(huì)形成很大的安全隱患[1]。強(qiáng)對(duì)流天氣對(duì)航班飛行的影響是制約中國民航運(yùn)輸持續(xù)發(fā)展的主要原因之一，據(jù)民航局公布的數(shù)據(jù)顯示，構(gòu)成航班延誤的各類因素中天氣原因占比最高為56.8%[2]。

飛行限制區(qū)的劃設(shè)方法主要分為靜態(tài)和動(dòng)態(tài)，由于氣象是不斷變化的，所以目前在航空領(lǐng)域的氣象預(yù)測(cè)主要以動(dòng)態(tài)預(yù)測(cè)為主[3-4]。Krozel[5]在對(duì)不良?xì)夂蛴绊懴碌娘w行受限區(qū)進(jìn)行劃設(shè)時(shí)，給出了采用多邊形對(duì)飛行受限區(qū)進(jìn)行劃設(shè)的方式。Bokadia 等[6]提出用平移劃設(shè)法來預(yù)測(cè)強(qiáng)對(duì)流天氣中的飛行受限區(qū)。Nilim 等[7]通過馬爾可夫預(yù)測(cè)模型對(duì)惡劣天氣進(jìn)行預(yù)測(cè)，將受惡劣天氣影響的扇區(qū)分為各個(gè)小區(qū)域，并通過預(yù)測(cè)雷暴的范圍來對(duì)飛行限制區(qū)進(jìn)行動(dòng)態(tài)劃設(shè)。謝春生等[8]提出了基于Graham 算法的初始飛行受限區(qū)劃設(shè)方法，得到了一種不規(guī)則的多邊形限制區(qū)，在劃設(shè)受限區(qū)時(shí)加入了動(dòng)態(tài)因素例如雷雨的移動(dòng)、雷雨大小的變化以及軍事活動(dòng)等，將靜態(tài)劃設(shè)方法擴(kuò)展到動(dòng)態(tài)。楊惠東等[9]結(jié)合雷達(dá)回波特征，將氣象數(shù)據(jù)進(jìn)行“比較賦值”運(yùn)算，對(duì)處理后的數(shù)據(jù)圖進(jìn)行灰度化、二值化以便進(jìn)行飛行受限區(qū)的邊緣探測(cè)，運(yùn)用Graham 算法將由上述邊緣探測(cè)得到的邊界點(diǎn)圍成凸多邊形，在氣象預(yù)報(bào)的基礎(chǔ)上提出了動(dòng)態(tài)飛行受限區(qū)的劃設(shè)方法。陳可嘉等[10]通過對(duì)影響飛機(jī)飛行的氣象屬性進(jìn)行研究并且找到最關(guān)鍵的因素，之后在Graham 算法的基礎(chǔ)上結(jié)合灰色關(guān)聯(lián)度分析的預(yù)測(cè)模型對(duì)飛行受限區(qū)進(jìn)行動(dòng)態(tài)劃設(shè)。蔣昕等[11]通過使用線性回歸的預(yù)測(cè)方法對(duì)飛行受限區(qū)域的邊界進(jìn)行實(shí)時(shí)預(yù)測(cè)。陳金良等[12]通過分析飛機(jī)對(duì)空間的利用時(shí)長和范圍，利用平移法對(duì)原有的飛行限制區(qū)預(yù)測(cè)，再根據(jù)外推法得出動(dòng)態(tài)飛行限制區(qū)，使用數(shù)學(xué)模型來描述飛行限制區(qū)的邊界變動(dòng)情形。由于氣象變化是隨機(jī)的且沒有固定的變化規(guī)律，所以采用常規(guī)的預(yù)測(cè)方法對(duì)動(dòng)態(tài)飛行限制區(qū)進(jìn)行預(yù)測(cè)會(huì)產(chǎn)生較大的誤差，且該難點(diǎn)較難克服，因此本文是在現(xiàn)有方法的基礎(chǔ)上綜合了雷暴的實(shí)時(shí)變化對(duì)飛行限制區(qū)進(jìn)行動(dòng)態(tài)預(yù)測(cè)。

各個(gè)地區(qū)采用不同設(shè)備對(duì)氣象數(shù)據(jù)進(jìn)行采集會(huì)導(dǎo)致氣象數(shù)據(jù)的時(shí)間分辨率不同，時(shí)間分辨率的高低對(duì)預(yù)測(cè)結(jié)果也有一定的影響[13]，因此針對(duì)不同時(shí)間分辨率的氣象數(shù)據(jù)應(yīng)該采用不同的方法進(jìn)行動(dòng)態(tài)飛行限制區(qū)的預(yù)測(cè)。本文所使用的數(shù)據(jù)為每小時(shí)更新一次，主要針對(duì)時(shí)間分辨率較低的氣象數(shù)據(jù)進(jìn)行動(dòng)態(tài)飛行限制區(qū)的預(yù)測(cè)方法研究。

1 動(dòng)態(tài)飛行限制區(qū)的預(yù)測(cè)

1.1 靜態(tài)飛行限制區(qū)的劃設(shè)

據(jù)反射率因子（REF）的大小可以確定雷達(dá)的回波等級(jí)，并依此將空中區(qū)域劃分為不同的范圍。依據(jù)反射率信息可以將雷達(dá)回波分為7 個(gè)等級(jí), 并給出在對(duì)應(yīng)回波等級(jí)下航空器的飛行狀態(tài)信息[14],如表1 所示。

表1 雷達(dá)回波等級(jí)表Tab.1 Diagram of radar echoes

首先對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，表2 為收集到的部分原始雷達(dá)數(shù)據(jù)。

表2 部分原始雷達(dá)數(shù)據(jù)Tab.2 Partial primary radar data

提取REF≥41 dBz 的回波強(qiáng)度，并對(duì)提取出的數(shù)據(jù)用Graham 算法劃設(shè)出靜態(tài)飛行限制區(qū)的初始多邊形。Graham 掃描法是用來尋找多邊形比較普遍且方便高效的一個(gè)方法，它的原理是先找到多邊形上的一個(gè)點(diǎn)，然后將其他點(diǎn)與該點(diǎn)進(jìn)行極角排序，按逆時(shí)針方向去找多邊形上的每一個(gè)點(diǎn)，符合條件的就存在多邊形點(diǎn)集里，否則就舍棄，具體的算法步驟如下。

輸入：平面直角坐標(biāo)系中散亂的點(diǎn)集；

輸出：點(diǎn)集構(gòu)成的凸包。

Step 1 將點(diǎn)集內(nèi)的所有點(diǎn)放在平面直角坐標(biāo)系中，找到最左下方的點(diǎn)，即橫縱坐標(biāo)最小的點(diǎn)，由幾何知識(shí)可知該點(diǎn)一定是多邊形上的點(diǎn)，如圖3 中的紅色點(diǎn)即為最左下方的點(diǎn)。

Step 2 然后以這點(diǎn)為極點(diǎn), 計(jì)算出其他所有點(diǎn)與這點(diǎn)的極角并從小到大進(jìn)行排序，如果出現(xiàn)多個(gè)點(diǎn)與該點(diǎn)的極角相同則根據(jù)其他點(diǎn)到該點(diǎn)的距離進(jìn)行排序，距離近的排在前面。

Step 3 由幾何知識(shí)可知點(diǎn)1 在多邊形上，將其加入多邊形的集合中，此時(shí)多邊形集合中有0 和1兩點(diǎn)，接下來找第三點(diǎn)，此時(shí)點(diǎn)2 稱為當(dāng)前點(diǎn)，連接點(diǎn)1 和當(dāng)前點(diǎn)得到直線L，若點(diǎn)3 在直線L 的左側(cè)則保留當(dāng)前點(diǎn)到多邊形集合中，若在右側(cè)則當(dāng)前點(diǎn)不在多邊形集合中，該點(diǎn)集中的點(diǎn)2 在多邊形集合內(nèi)。

Step 4 接下來找第四個(gè)點(diǎn)，此時(shí)點(diǎn)3 為當(dāng)前點(diǎn)，連接點(diǎn)2 和當(dāng)前點(diǎn)得到直線M，因?yàn)辄c(diǎn)4 在直線M 的右側(cè)所以當(dāng)前點(diǎn)不在多邊形集合中，所以連接點(diǎn)2 和點(diǎn)4 得到直線N，此時(shí)點(diǎn)4 為當(dāng)前點(diǎn)，因?yàn)辄c(diǎn)5 在直線N 的左側(cè)，所以當(dāng)前點(diǎn)在多邊形集合內(nèi)，又因?yàn)辄c(diǎn)6 在直線N 的右側(cè)，所以當(dāng)前點(diǎn)不在多邊形集合內(nèi)，點(diǎn)6 在多邊形集合內(nèi)，依次對(duì)各個(gè)點(diǎn)進(jìn)行判斷。

Step 5 依次對(duì)每個(gè)點(diǎn)進(jìn)行判斷，直到點(diǎn)集中的最后一個(gè)點(diǎn)時(shí)，連接點(diǎn)0 與最后一個(gè)點(diǎn)。

Step 6 將以上多邊形集里的點(diǎn)依次連接即可得到該點(diǎn)集所對(duì)應(yīng)的凸包，如圖1 為平面直角坐標(biāo)系中一些散亂的點(diǎn)構(gòu)成的凸包。

圖1 點(diǎn)集構(gòu)成的凸包Fig.1 Convex package composed of the above point set

1.2 動(dòng)態(tài)飛行限制區(qū)的預(yù)測(cè)

對(duì)飛行限制區(qū)的動(dòng)態(tài)預(yù)測(cè)主要分為4 個(gè)過程：首先是采用Graham 算法劃設(shè)出靜態(tài)飛行限制區(qū)的初始多邊形；然后由歷史數(shù)據(jù)得到的距離均值的方法對(duì)飛行限制區(qū)的形狀進(jìn)行預(yù)測(cè)；其次引入Markov思想通過類狀態(tài)轉(zhuǎn)移矩陣對(duì)飛行限制區(qū)的中心點(diǎn)位置進(jìn)行預(yù)測(cè)；最后根據(jù)歷史數(shù)據(jù)得到的角度增量來預(yù)測(cè)飛行限制區(qū)中心點(diǎn)的角度變化。

在對(duì)飛行限制區(qū)進(jìn)行劃設(shè)時(shí)，除了考慮雷暴中心的移動(dòng)趨勢(shì)，雷暴本身的變化也是影響飛機(jī)飛行的一個(gè)重要因素，研究雷暴變化的過程中，REF=35 dBz 時(shí)表明該天氣狀況需要引起注意，REF=41 dBz 時(shí)已經(jīng)嚴(yán)重影響到航空交通。本文通過對(duì)歷史數(shù)據(jù)進(jìn)行分析，找到了飛行受限區(qū)劃設(shè)區(qū)域35 dBz≤REF≤41 dBz 的多邊形形變關(guān)系，并基于此給出了一種對(duì)飛行限制區(qū)進(jìn)行動(dòng)態(tài)預(yù)測(cè)的方法即距離均值的方法。

首先通過Graham 算法分別畫出當(dāng)前時(shí)刻REF≥41 dBz 和REF≥35 dBz 時(shí)的飛行限制區(qū)初始多邊形。假設(shè)圖2 為當(dāng)前REF≥41 dBz 和REF≥35 dBz 時(shí)的多邊形ABCDE 和A'B'C'D'E'F'，計(jì)算出多邊形的中心點(diǎn)O。由歷史氣象數(shù)據(jù)可得到下一時(shí)刻REF≥41 dBz 時(shí)畫出的多邊形形狀與當(dāng)前時(shí)刻REF≥35 dBz 時(shí)畫出的多邊形形狀的關(guān)系,分別計(jì)算出REF≥41 dBz 和REF≥35 dBz 時(shí)多邊形上各頂點(diǎn)與中心點(diǎn)的距離之和

圖2 飛行受限區(qū)同一時(shí)刻不同范圍REF 示意圖Fig.2 Schematic diagram of different ranges of REF in flight restricted areas at the same time

氣象變化是隨機(jī)的沒有固定的變化規(guī)律，所以認(rèn)為其下一時(shí)刻的氣象狀況只與上一時(shí)刻的氣象狀況有關(guān)，而與之前更早的狀況無關(guān)，即符合Markov 的無后效性，要對(duì)整個(gè)飛行限制區(qū)進(jìn)行動(dòng)態(tài)預(yù)測(cè)則需要找到符合變化規(guī)律的類Markov 狀態(tài)轉(zhuǎn)移矩陣。令兩個(gè)中心點(diǎn)分別為A（x1，y1）與B（x2，y2），且Markov 類狀態(tài)轉(zhuǎn)移矩陣為P，則可得

對(duì)中心點(diǎn)進(jìn)行預(yù)測(cè)時(shí)將每個(gè)點(diǎn)分為橫向預(yù)測(cè)和縱向預(yù)測(cè)，如圖3 的點(diǎn)A，在橫向可以向前即M點(diǎn)移動(dòng)，也可以向后即N 點(diǎn)移動(dòng)，縱向可以向上即D 點(diǎn)移動(dòng)，也可以向下即E 點(diǎn)移動(dòng)。

圖3 飛行限制區(qū)中心點(diǎn)位移示意圖Fig.3 Schematic diagram of the displacement of the center point of the flight restriction zone

在基于距離均值雷暴形變預(yù)測(cè)的基礎(chǔ)上考慮了雷暴運(yùn)動(dòng)的角度變化，由上一時(shí)刻的中心點(diǎn)與當(dāng)前時(shí)刻的中心點(diǎn)的角度變化來確定整個(gè)飛行受限區(qū)移動(dòng)的角度，在計(jì)算時(shí)找一個(gè)飛行限制區(qū)多邊形以外的點(diǎn)O（x0，y0）作為參考點(diǎn)。

圖4 飛行受限區(qū)角度變化示意圖Fig.4 Schematic diagram of the change in the angle of the flight restricted area

2 實(shí)例驗(yàn)證

記當(dāng)前時(shí)刻為t1，上一時(shí)刻為t0，下一時(shí)刻為t2，則第n 時(shí)刻為tn-1。首先提取出t1時(shí)刻REF≥41 dBz和REF≥35 dBz 的數(shù)據(jù)，根據(jù)航空規(guī)定危險(xiǎn)天氣分布之間大于20 km 時(shí)，可以將飛行受限區(qū)劃分為多個(gè)更小的區(qū)域，所以在劃設(shè)靜態(tài)飛行限制區(qū)時(shí)對(duì)原始數(shù)據(jù)畫出的靜態(tài)飛行限制區(qū)進(jìn)行修改，畫出t1時(shí)刻的初始多邊形。

由此可以得到靜態(tài)飛行限制區(qū)的初始形狀，根據(jù)規(guī)定我國航路寬度為20 km，即航路中心線兩側(cè)各10 km，所以最終的飛行受限區(qū)應(yīng)將多邊形的邊界向外擴(kuò)10 km。接下來需要求出多邊形的幾何中心點(diǎn)，而是否增加向外擴(kuò)的部分對(duì)幾何中心點(diǎn)都沒有影響，在不考慮外擴(kuò)的情況下對(duì)靜態(tài)飛行限制區(qū)進(jìn)行動(dòng)態(tài)預(yù)測(cè)。

計(jì)算出t0，t1時(shí)刻REF≥41 dBz 時(shí)多邊形的幾何中心點(diǎn)，中心點(diǎn)的計(jì)算公式如下

可得到t0時(shí)刻和t1時(shí)刻多邊形的中心坐標(biāo)。

由以上步驟計(jì)算出的t1時(shí)刻REF≥41 dBz 和REF≥35 dBz 時(shí)多邊形頂點(diǎn)坐標(biāo)以及t1時(shí)刻REF≥41 dBz 時(shí)多邊形的中心點(diǎn)坐標(biāo)，計(jì)算出中心點(diǎn)到兩個(gè)多邊形各個(gè)頂點(diǎn)的距離，令REF≥41 dBz 時(shí)中心點(diǎn)到多邊形各頂點(diǎn)的距離和為

通過觀察t0、t1時(shí)刻以及之后更多時(shí)刻REF≥41 dBz 和REF≥35 dBz 時(shí)多邊形頂點(diǎn)坐標(biāo)，選取了O （x0，y0）點(diǎn)作為本次預(yù)測(cè)的參考點(diǎn)，令t0、t1時(shí)刻REF≥41 dBz 時(shí)兩個(gè)多邊形中心點(diǎn)分別為A（x1，y1）與B（x2，y2），則點(diǎn)A 與點(diǎn)B 之間的夾角分別為α 和β，其中

3 實(shí)驗(yàn)結(jié)果

由上文飛行限制區(qū)的動(dòng)態(tài)預(yù)測(cè)步驟可以得到t1時(shí)刻REF≥41 dBz 時(shí)飛行限制區(qū)的預(yù)測(cè)多邊形，將其與t2時(shí)刻REF≥41 dBz 時(shí)實(shí)時(shí)的飛行限制區(qū)進(jìn)行對(duì)比，如表3 所示。

表3 飛行受限區(qū)多邊形頂點(diǎn)坐標(biāo)對(duì)比Tab.3 Vertice coordinate comparison of flight restricted area polygons

畫出t2時(shí)刻REF≥41 dBz 時(shí)實(shí)時(shí)飛行限制區(qū)的多邊形以及預(yù)測(cè)多邊形并進(jìn)行對(duì)比，如圖5所示。

圖5 t2 時(shí)刻飛行限制區(qū)動(dòng)態(tài)預(yù)測(cè)對(duì)比圖Fig.5 Comparison chart of dynamic predictions in flight restriction zones at t2

采用以上提出的方法對(duì)另一組氣象數(shù)據(jù)進(jìn)行預(yù)測(cè)并得到預(yù)測(cè)結(jié)果如圖6 所示。

由圖5 和圖6 可以看出劃設(shè)的動(dòng)態(tài)飛行限制區(qū)幾乎覆蓋了t2時(shí)刻REF≥41 dBz 時(shí)飛行限制區(qū)。通過精確度和偏離率對(duì)危險(xiǎn)天氣飛行限制區(qū)的動(dòng)態(tài)預(yù)測(cè)作出評(píng)價(jià)，精確度主要用來評(píng)估所劃設(shè)動(dòng)態(tài)飛越限制區(qū)的準(zhǔn)確率，即實(shí)際飛行限制區(qū)與所劃設(shè)的動(dòng)態(tài)飛行限制區(qū)重疊部分的比率，精確度越高則說明預(yù)測(cè)的飛行限制區(qū)越準(zhǔn)確，用Da表示；偏差率表示在動(dòng)態(tài)飛行限制區(qū)域漏劃面積和實(shí)際飛越限制區(qū)域面積之間的比率，偏差率越小則說明劃設(shè)的動(dòng)態(tài)飛行限制區(qū)覆蓋的實(shí)際飛行限制區(qū)的區(qū)域越大，則對(duì)飛機(jī)飛行越有安全保障，用Db表示，其中Da與Db的計(jì)算公式如下

圖6 t2 時(shí)刻飛行限制區(qū)動(dòng)態(tài)預(yù)測(cè)對(duì)比圖Fig.6 Comparison chart of dynamic predictions in flight restriction zones at t2

式中，SAFRA為實(shí)際飛行受限區(qū)面積，SDFRA為劃設(shè)的動(dòng)態(tài)飛行限制區(qū)的面積，本文實(shí)驗(yàn)得到的結(jié)果如表4。

表4 預(yù)測(cè)結(jié)果對(duì)比Tab.4 Comparison of prediction results

通過參考文獻(xiàn)所提供的預(yù)測(cè)方法并與之對(duì)比，例如蔣昕等[11]中的一元線性回歸預(yù)測(cè)方法所給出的結(jié)果，其精確度均值為69.41%，偏差度均值為24.00%，而本文精確度為74.72%，偏差度為10.11%，可知本文所提出的方法精確度較高，偏差度較低，有較大的應(yīng)用價(jià)值與發(fā)展前景。

4 結(jié)論

本文在Graham 算法和Markov 轉(zhuǎn)移矩陣的基礎(chǔ)上提出了距離均值和角度增量的方法對(duì)飛行限制區(qū)進(jìn)行動(dòng)態(tài)預(yù)測(cè)，不僅考慮了雷雨天氣雷暴中心點(diǎn)的位置變化，對(duì)其自身的形變也進(jìn)行了預(yù)測(cè)，并得到了以下的結(jié)論：

1）本文所提出的預(yù)測(cè)方法在預(yù)測(cè)時(shí)間精度較低的氣象數(shù)據(jù)時(shí)精確度較高，能有效覆蓋實(shí)際的飛行限制區(qū)，偏差率較低，預(yù)測(cè)的誤差較小，對(duì)于氣象預(yù)測(cè)以及航空安全方面有較大的應(yīng)用價(jià)值。

2）本文所提出的預(yù)測(cè)方法在對(duì)氣象進(jìn)行預(yù)測(cè)時(shí)不僅考慮了雷暴中心點(diǎn)的位置移動(dòng)情況，而且也對(duì)雷暴本身的形變進(jìn)行了預(yù)測(cè)，使得預(yù)測(cè)結(jié)果更貼近實(shí)際情況。

3）本文所提出的預(yù)測(cè)方法在對(duì)氣象進(jìn)行預(yù)測(cè)時(shí)可以做到實(shí)時(shí)更新預(yù)測(cè)區(qū)域，具有實(shí)際且廣闊的應(yīng)用前景。