基于雷達數據的飛行受限區形變及方位預測

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

（1.華北理工大學理學院，河北 唐山 063210；2.華北理工大學建筑工程學院，河北 唐山 063210；3.華北理工大學河北省數據科學與應用重點實驗室，河北 唐山 063210）

隨著科技的發展和人民生活水平的提升，人們使用航空器的次數及貨物運輸的頻率也在日益增多，同時空域資源也越來越緊張，一旦出現不良氣象或是重大事件，空域資源的利用率將會急劇下降，進而造成飛機晚點、乘客滯留甚至航班取消，這不但會造成巨大的經濟損失，也會形成很大的安全隱患[1]。強對流天氣對航班飛行的影響是制約中國民航運輸持續發展的主要原因之一，據民航局公布的數據顯示，構成航班延誤的各類因素中天氣原因占比最高為56.8%[2]。

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

各個地區采用不同設備對氣象數據進行采集會導致氣象數據的時間分辨率不同，時間分辨率的高低對預測結果也有一定的影響[13]，因此針對不同時間分辨率的氣象數據應該采用不同的方法進行動態飛行限制區的預測。本文所使用的數據為每小時更新一次，主要針對時間分辨率較低的氣象數據進行動態飛行限制區的預測方法研究。

1 動態飛行限制區的預測

1.1 靜態飛行限制區的劃設

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

表1 雷達回波等級表Tab.1 Diagram of radar echoes

首先對數據進行預處理，表2 為收集到的部分原始雷達數據。

表2 部分原始雷達數據Tab.2 Partial primary radar data

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

輸入：平面直角坐標系中散亂的點集；

輸出：點集構成的凸包。

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

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

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

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

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

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

圖1 點集構成的凸包Fig.1 Convex package composed of the above point set

1.2 動態飛行限制區的預測

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

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

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

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

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

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

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

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

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

2 實例驗證

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

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

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

可得到t0時刻和t1時刻多邊形的中心坐標。

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

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

3 實驗結果

由上文飛行限制區的動態預測步驟可以得到t1時刻REF≥41 dBz 時飛行限制區的預測多邊形，將其與t2時刻REF≥41 dBz 時實時的飛行限制區進行對比，如表3 所示。

表3 飛行受限區多邊形頂點坐標對比Tab.3 Vertice coordinate comparison of flight restricted area polygons

畫出t2時刻REF≥41 dBz 時實時飛行限制區的多邊形以及預測多邊形并進行對比，如圖5所示。

圖5 t2 時刻飛行限制區動態預測對比圖Fig.5 Comparison chart of dynamic predictions in flight restriction zones at t2

采用以上提出的方法對另一組氣象數據進行預測并得到預測結果如圖6 所示。

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

圖6 t2 時刻飛行限制區動態預測對比圖Fig.6 Comparison chart of dynamic predictions in flight restriction zones at t2

式中，SAFRA為實際飛行受限區面積，SDFRA為劃設的動態飛行限制區的面積，本文實驗得到的結果如表4。

表4 預測結果對比Tab.4 Comparison of prediction results

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

4 結論

本文在Graham 算法和Markov 轉移矩陣的基礎上提出了距離均值和角度增量的方法對飛行限制區進行動態預測，不僅考慮了雷雨天氣雷暴中心點的位置變化，對其自身的形變也進行了預測，并得到了以下的結論：

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

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

3）本文所提出的預測方法在對氣象進行預測時可以做到實時更新預測區域，具有實際且廣闊的應用前景。