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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 實例驗證

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

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

計算出t0，t1時刻REF≥41 dBz 時多邊形的幾何中心點，中心點的計算公式如下

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

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

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

3 實驗結(jié)果

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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