基于雷达数据的飞行受限区形变及方位预测

樊立艳，赵鹤宇，常锦才,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）本文所提出的预测方法在对气象进行预测时可以做到实时更新预测区域，具有实际且广阔的应用前景。