基于危险天气不确定性的最小风险路径规划方法

王岩韬，赵昕颐

中国民航大学国家空管运行安全技术重点实验室，天津 300300

航路危险天气（以下简称危险天气）指在巡航过程中直接或间接影响飞行安全的天气现象，包括雷暴、飞机颠簸和飞机积冰等，是造成飞行事故和飞行延误的主要原因之一. 据统计，近年来欧美地区四分之一的航班延误是由危险天气所导致的[1-2]；在国内，天气原因更是导致了超半数以上的不正常航班[3]. 通过现有的天气观测和预报业务，提前掌握危险天气影响范围，合理规划飞行航迹，可有效提高飞行安全、减少飞行延误. 但受限于观测数据和预报模型的误差，以及大气自身的混沌特性，危险天气的预测结果并不完全准确，因此称危险天气具有一定的不确定性.

危险天气区域的规避方式从时间尺度上可以分为起飞前的飞行计划和飞行中的改航路径规划. 飞行中的改航路径规划通常发生相遇前10～60 min[4-5]，主要依靠气象雷达和卫星的实时探测，以及基于雷达回波或卫星资料的区域外推预报技术，在短时范围内预测危险天气的位置、形状和移动[6]，从而提前规避危险天气区域. 利用区域外推技术进行改航规划的典型研究包括Hauf等[7]提出的DIVMET模型，该模型将雷达和卫星检测到的雷暴区域划设为二维多边形区域，并考虑了绕飞安全裕度与危险天气的移动，可以在二维平面内给出引导飞机绕过或穿越危险天气区域的改航方案；国内王岩韬和刘锟[8-9]通过统计气象数据与历史航迹，将垂直液态水含量和基本反射率作为雷暴区域的判断指标，划设飞行受限区，并提出了一种面向航迹运行(TBO)的4D改航回归航迹规划方法. 综合以往研究来看[10]，飞行中改航路径规划的主要研究思路是使用气象雷达产品中的基本反射率等指标[11]或卫星云图将危险天气区域划设为形状规则的椭圆[12]、多边形[13]或栅格形式的飞行受限区，采用图搜索算法（Dijkstra算法、A*和D*等[14]）、随机采样算法（RRT、RRT*等[15]）等进行路径规划. 该层面的研究主要存在以下三个问题：一是多数研究仅考虑了雷暴天气，不能代表航路上全部危险天气，少数研究对积冰[16]、颠簸[17]等天气现象做了专项研究，但均未能综合考虑各种危险天气进行路径规划；二是基于区域外推的预测方法缺乏对危险天气的发生、发展和消亡的物理机制描述，无法定量描述天气发展的不确定性，从而导致模型的预测误差随时间的推移而迅速增加；三是在路径规划过程中，仅以改航路径最短、油耗最小等经济性指标为目标，缺乏安全性考量.

为了减少航班临时改航所带来的安全隐患和经济损失，飞行前对航路危险天气的提前感知和精细的飞行计划可以使航班在很大程度上避免遭遇危险天气，但这需要时间提前量更长、准确度更高的天气预报. 相较于1 h内的区域外推预报，该时间尺度的预报主要从危险天气的发生机理和所依赖的环境条件出发，根据不同诊断物理量对相应危险天气的指示意义，来进行分类天气预报[18]，可以帮助航空公司和机组人员更早地了解潜在的危险天气，提前调整飞行计划和航线. Doswell等[19]根据雷暴发生的必要条件—对流不稳定性、水汽条件和抬升动力，提出的“配料法”被广泛应用于危险天气的短时预报，并逐渐发展成为基于多种天气资料的数值天气预报（NWP）. 作为一种确定性预报方法，数值天气预报使用动力学和热力学基本方程来模拟真实大气系统的演变，由此产生对未来天气可能的最佳估计. 但由于大气的混沌特性，初始数据的观测误差和数值模拟过程中的近似误差会随时间被无限放大，因此数值天气预报同样具有一定的不确定性. 为了准确描述这种不确定性，集合预报（EPS）将多个数值模型的预测结果转化为概率密度函数，由此计算得到各类危险天气发生的概率[20]. 近年来，航空领域的许多研究是在此基础上，结合控制论、安全管理等理论方法为航班提前规划轨迹. 例如，Soler和其他学者[4,21-22]讨论了风、雷暴等对流天气的不确定性对飞行航迹的影响，使用总指数（Total totals index,TT）和对流降水（Convective precipitation，CP）的概率预报来表征风和潜在对流区的不确定性，并提出了TBO场景下具有鲁棒性的航迹规划方法；王晓亮等[23]使用集合预报中的概率天气信息构建空域概率天气模型，对通用航空的飞行路径进行了规划；Zhang[24]将雷达反射率和上升气流垂直速度两个参数作为雷暴天气的预测指标，通过计算二者的联合概率分布得到雷暴的可能性，从而进行无人机的任务规划和风险评估. 综上，概率预报在危险天气下的航迹提前规划和风险管控领域具有良好的应用前景，但仍存在两点问题：一是雷暴的不稳定性条件和水汽量评估指标单一，且仅进行了简单的联合概率分析，未能充分体现雷暴的产生机理，从单一天气要素的概率预报到雷暴天气的概率预报仍需建立更加可靠、可解释的映射关系；二是现有基于概率预报的不确定性天气研究仅考虑了雷暴等强对流天气，而对积冰、颠簸等危险天气研究不足.

鉴于现有研究存在的不足，本文针对飞行中最易遭遇的危险天气风险，基于数值预报和集成预报，在飞行计划层面提出了一种最小风险路径规划方法. 其改进之处在于：（1）为解决现有规划问题的飞行受限区只考虑单一危险天气的问题，分别对雷暴、积冰和颠簸三类危险天气区域进行了针对性的划设，并对三区进行融合，建立一种具备危险天气风险标识的栅格化地图—风险地图；（2）针对雷暴预测指标的单一的问题，选取对流有效位能和K指数、累计降水和组合雷达反射率分别描述大气的不稳定性和水汽量，考虑是否具有抬升触发条件来预测雷暴的发生，并结合C-F模型建立从单一气象要素概率预报到雷暴天气概率预报的映射关系；（3）考虑路径风险最小进行路径规划，提出一种风险最小化的A*算法和Dijkstra算法. 综上，使用传统Dijkstra、A*、RRT算法和本文提出的风险最小化算法对实际案例进行路径规划，寻找总风险最小下的最优路径.

1 多类型危险天气区域划设

1.1 基于概率预报的雷暴区域划设方法

民航规章规定签派员至少需在航班预计起飞时刻前1 h提交飞行计划，一般情况下航班会在未来2～12 h内进入稳定的高空巡航阶段（如图1所示），而雷暴等强对流天气具有突发性强、生命史短等特点，依靠区域外推预报方法无法在过长的时间提前量内预知雷暴的发生. 因此此处利用集合天气预报中的概率预报，结合雷暴的产生机理，使用“配料法”在起飞前推断航班在高空航路巡航时可能会遇到的雷暴区域，实现在较长时间尺度内预测在航路中遭遇雷暴的风险.

集合预报由多个数值预测模型组成，例如中国气象局的区域集合预报系统（CMA-REPS）由15个集合成员组成（1个控制预报+14个扰动成员），每天预报两次，分别在UTC时间0：00和12：00，单次预报时效为0、3、6、···、72 h，预报要素包括累计降水、K指数、CAPE指数等，可为突发性、危险性的天气过程预报提供依据. 该系统中的概率预报模式，通过统计预测结果超过预报阈值的集合成员占全部集合成员的比例，即可得到某预报要素大于特定阈值的概率，本文在此预报结果基础上进行雷暴天气的诊断和预测.

“配料法”的核心思想是通过建立预报要素与天气现象之间的映射关系来实现天气预报. 因此，根据雷暴的产生需要三个必要条件，即水汽条件、不稳定条件以及抬升触发机制[25]，选取6 h累计降水、组合雷达反射率、CAPE指数（对流有效位能）、K指数作为定量预报要素，是否具有抬升触发机制作为定性预报要素. 预报要素的选取与区域所在地理位置息息相关，上述要素是根据我国气候特征，综合考虑数据的可得性和可用性所确定的，多项研究表明上述要素对雷暴的预测具有良好的指示性[6,25].

在获取了上述预报要素的概率预报结果后，根据配料法思想，需依据预报要素与天气现象之间的物理关系进行预报建模. 鉴于雷暴天气过程中存在的多种不确定性，为了直观地描述不确定性在预报过程中的传播和融合，此处采用CF（Certainty factor）模型[26]，依据雷暴的产生机理建立从预报要素到雷暴天气的映射关系，以此计算雷暴发生概率.

C-F模型是一种基于可信度的不确定性推理方法，模型主要包括知识不确定性和证据不确定性两种不确定性. C-F模型的一般采用产生式规则：

其中，E代表知识的前提证据；H代表知识结论；CF(H,E)表示知识的可信度，即当证据E为真时，E对结论为H的支持程度或支撑程度；另外使用CF(E)表示证据的不确定性，即证据E的可信度. CF(H,E)、CF(E)取值范围均为[–1,1].

为了表示由于天气预报误差所导致的证据不确定性，使用概率预报中的预报结果p表示该要素在未来一段时间处于某取值范围的可信度CF(E).例如，“某格点区域未来6 h累积降水≥1 mm的概率为50%”表示“某格点区域未来6 h累积降水≥1 mm”的可能性为50%，即该天气信息证据的可信度为50%.

假设天气事件：

E1：对流有效位能CAPE≥ε1;

E2：K指数≥ε2;

E3：6 h累计降水≥ε3;

E4：组合雷达反射率≥ε4;

H1：大气不稳定;

H2：水汽充足;

H3：具有抬升触发机制;

T：雷暴发生.

则根据雷暴的发生条件[6,25]产生规则如下：

IFE1orE2THENH1(CF(H1,E12))；

IFE3orE4THENH2(CF(H2,E34))；

IFH1andH2andH3THENT(CF(T,H123)).

其中，CF(H1,E12)表示当“对流有效位能CAPE≥ε1”和“K指数≥ε2”中的一个事件为真时“大气不稳定”为真的可信度，CF(H2,E34)同理；CF(T,H123)表示当“大气不稳定”“水汽充足”“具有抬升触发机制”三者同时为真时“雷暴发生”的可信度.

根据C-F模型中证据不确定性的合取和析取规则，多证据融合时的证据不确定性计算公式如下：

其中，ε1～ε4为导致雷暴发生的预报要素阈值，其取值随区域、季节和天气型的不同有着显著差异.例如，经统计表明[27]，引发雷暴天气的CAPE指数和K指数具有明显的季节特征，冬春最小到夏季达到峰值，且具有南北差异，南方地区的阈值高于北方地区. 因此预报要素阈值应根据具体地区、时间进行设定.

在建立了雷暴产生的规则后，使用C-F模型中结论不确定性的传递算法，从不确定的初始证据出发，最终推出结论并求出结论的可信度值. 结论H的可信度的计算公式如下：

由此计算得到雷暴的发生概率P1：

1.2 积冰区域划设方法

飞机在高空中产生积冰的主要原因是云体及其周围环境中存在大量水汽和过冷水滴，一旦该区域温度达到积冰范围，则有可能出现积冰风险.因此使用数值预报产品中的相对湿度RH和温度T参数，结合国际民航组织推荐的飞机积冰指数Ic计算积冰概率，积冰指数的计算公式为：

式中，RH为相对湿度，T为温度. 公式的前半部分用相对湿度线性拟合水滴的数量、大小的增长过程，RH越接近100，取值越接近最大值100；公式的后半部分用温度的二次方来拟合水滴的增长率，当T= –7 ℃时为最大值1，T= –14 ℃和0 ℃时取最小值0. RH低于50%或者T超过0～–14 ℃范围的，水滴增长率为0，即认为无积冰发生. 因此，积冰指数Ic输出范围为0～100%，数值越大，表示积冰越强：RH=100%和T= -7 ℃时，积冰指数取得最大值100%. 具体积冰强度判据为: 当0%≤Ic＜50%，预报有轻度积冰；当50%≤Ic＜80%，预报有中度积冰;当Ic≥80%，预报有严重积冰，如表1所示.

表1 积冰指数预报等级及阈值Table 1 Forecast level and threshold value of the Ice index

该网格内积冰的概率P2：

1.3 颠簸区域划设方法

飞机在高空航路中遭遇的颠簸通常是由晴空湍流所导致的，民航规章规定使用湍流耗散率（EDR）作为颠簸的观测标准. 常见的颠簸预测指数包括Ellord指数、Dutton指数、Brown指数和理查森数等，研究表明Dutton指数对于我国中东部及沿海地区高空急流背景下的颠簸具有良好的指示能力[28]，因此使用数值预报中的风参数计算Dutton指数，以此作为颠簸的预测指标.

Dutton提出的基于水平风切变和垂直风切变的经验指数计算公式为：

式中，10.5为经验常数，SH为水平风切变（10–5s–1），SV为垂直风切变（10–3s–1），二者的计算公式如下：

根据Dutton颠簸指数预报等级，如表2所示，Dutton指数小于等于20时认为无颠簸发生，此时颠簸概率P3=0；Dutton指数大于等于50时颠簸最为严重，认定该情况下颠簸概率P3=100%；当Dutton指数处于20～50之间，使用线性插值法计算颠簸概率P3，计算公式如下：

表2 Dutton颠簸指数预报等级及阈值Table 2 Forecast level and threshold value of the Dutton index

式中，E为Dutton指数，20和50分别为无颠簸和严重颠簸的预报阈值.

1.4 风险地图的建立

使用二维风险地图（Risk map）建立自由飞行下的空域模型，如图2所示. 风险地图是一种可标识风险的栅格化地图，其中每个栅格都包含一个风险成本，用于量化航空器飞越该位置的风险，本文使用航空器遭遇危险天气的概率P表示风险成本r，计算公式如下：

图2 风险地图的建立和融合Fig.2 Establishment and integration of the risk map

假设风险地图由n×n个正方形栅格组成，每个栅格的风险成本r(x,y)，x,y∈(1,n)，取值范围为[0,1]，其中0对应无风险区，1对应飞行风险最高的区域，即禁飞区. 为了保障飞行安全，规定危险天气概率P大于δ（0≤δ≤1）的区域r值为1，δ根据飞行员对危险的容忍度和飞行安全裕度进行相应调整，本文取50%.

风险地图的建立使用了国内两种气象预报模型的预报结果，一种是用于雷暴天气预测的区域集合预报模型（CMA-REPS），一种是用于积冰和颠簸指数计算的数值预报模型（CMA-MESO）. 考虑到两类预报模型的性能差异，风险地图的分辨率、起报时间、预报时效和更新频率取两类预报模型指标的交集，具体如表3所示.

表3 预报模型性能指标Table 3 Forecast model performance index

另一方面，本文统计了CMA-MESO和CMAREPS两个预报模型所提供的参数和预报高度，并与再分析观测模型进行了对比. 结果显示（表4），国家气象科学数据中心公布的CMA-MESO数值预报产品虽然支持气压、温度等气象要素的部分高度数据，但无巡航高度（350 hpa）预报数据，且相对湿度和经/纬向风仅支持2和10 m低空范围，无法表示高空气象条件. 同样，中央气象台发布的CMAREPS集合预报产品仅支持二维的数据预报结果，无法用于三维天气计算.

表4 气象模型对比Table 4 Comparison of meteorological models

考虑到三维数据的缺失问题，本文以二维地图为例，旨在提供一种天气预测与融合的思路，在未来气象数据完备的条件下可推广至三维天气的计算. 此外，由于现有数值预报产品仅提供2和10 m高度水平的风和相对湿度数据，第3节使用中国全球大气再分析产品（CRA）代替数值预报产品进行案例分析. CRA产品具备良好的垂直分辨率和完整的气象参数，可用于过去大气状况的真实重现.

2 风险最小化飞行路径规划方法

目前常见的路径规划算法包括Dijkstra算法、A*算法、RRT算法等已被广泛应用于飞行路径规划中，但绝大部分研究是以飞行路径最短为目标，寻找避障条件下的最优路径，没有考虑非障碍物环境中的飞行风险. 因此对传统Dijkstra算法、A*算法进行改进，基于风险地图，以路径风险最小为目标函数，寻找风险最小下的最优路径.

2.1 风险最小化Dijkstra算法

传统Dijkstra算法的代价函数是赋权有向图中两相连顶点之间边的权值，代表两顶点间的距离代价. 在风险最小化Dijkstra算法中，首先将栅格化的风险地图转化为赋权有向图，每个网格作为图中的一个顶点与周围八个顶点相连，然后将相连边上的权值赋为前点到后点的距离代价与后点对应栅格的风险值之和（做归一化处理后），以此作为代价函数代入传统Dijkstra算法求解.

2.2 风险最小化A＊算法

传统A*算法基于代价函数来评估每个节点的价值，从而选择下一个节点进行扩展. 该算法同时考虑了两个因素：一是从起始节点到当前节点的实际距离，二是从当前节点到目标节点的预估距离. A*算法的距离代价公式为：

式中，G(n)表示从起点开始移动到当前位置n的实际距离代价，H(n)表示从当前位置n移动到终点的预估距离代价，F(n)表示经过过程n后从起点到终点的距离代价. 预估距离H(n)通常采用曼哈顿距离或欧几里得距离等，本文根据实际飞行情况，考虑斜线距离，采用欧几里得距离计算预估距离，其计算公式为：

式中，(x1,x2)和(y1,y2)为两点坐标，d为两点间直线距离.

在风险最小化A*算法中，代价函数不仅包括两点间距离，还包括到达下一格点的风险代价，因此设定代价函数为：

式中，d0为起点到终点的距离，用于路径代价的归一化处理；R(n)为风险代价，即风险地图上对应格点的风险值r.

3 案例分析

2023年3月31日至4月5日，全国出现大范围强降水天气，其中江南南部、华中等地多次发生强对流天气. 本文采用中国标准时间（CST）2023年4月3日中央气象台发布的强对流天气预报、数值预报CMA区域集合模式，以及国家气象科学数据中心发布的中国全球大气再分析40 a产品（逐6 h产品）进行案例分析，案例区域范围为22°N～32°N，108°E～118°E. 根据案例区域建立100×100的网格地图，每个网格大小为0.1°×0.1°.

3.1 基于危险天气预测的风险地图生成

文献[29]指出当CAPE＜400 J·kg–1时出现雷暴发生概率仅为18.3%、K指数＜30 ℃时雷暴发生概率仅为11.8%，考虑中央气象台概率预报范围，选取CAPE≥500 J·kg–1，K指数≥30 ℃作为雷暴的判定标准. 同时设定6 h累计降水量大于等于4 mm，组合雷达反射率大于等于雷达图中黄色回波范围30 dbz，表示大气中具有一定的水汽条件，具体如下所示：

E1：对流有效位能CAPE≥500 J·kg–1;

E2：K指数≥30 ℃;

E3：6 h累计降水≥4 mm;

E4：组合雷达反射率≥30 dbz.

根据2023年4月3日中央气象台发布的CMA区域集合模式预报，未来12 h内实验地区“对流有效位能≥500 J·kg–1”、“K指数≥30 ℃”、“6 h累计降水≥4 mm”、“组合雷达反射率≥30 dbz”四个事件发生的概率如图3所示，图中灰度越大表示该事件发生的概率越高.

图3 2023.04.03实验地区雷暴诊断指标概率预报图. (a)对流有效位能≥500 J·kg–1; (b) K指数≥30 ℃; (c) 6 h累计降水≥4 mm; (d)组合雷达反射率≥30 dbzFig.3 Probability forecast chart of the diagnostic index of thunderstorms in the experimental area on April 3, 2023: (a) CAPE ≥ 500 J·kg–1; (b) K index≥ 30 ℃; (c) 6 h cumulative precipitation ≥ 4 mm; (d) combined radar reflectivity ≥ 30 dbz

将事件E1～E4的发生概率转化为可信度CF(E1)～CF(E4)，并设定CF(H1,E12)=80%，CF(H2,E34)=80%，CF(T,H123)=100%. 春季冷暖气团交汇频繁，华中地区中尺度对流活跃，因此具备一定的抬升触发条件，即CF(H3)=100%. 使用C-F推理模型结合1.1节中雷暴的发生条件产生规则，计算雷暴发生的可能性P2，计算结果可视化后如图4所示.

图4 2023.04.03实验地区雷暴概率图Fig.4 Thunderstorm probability graph in the experimental area on April 3, 2023

根据公式（3）、公式（5）～（7），使用中国全球大气再分析产品中350 hpa气压高度上的湿度（%）、温度（℃）、纬向风速(m·s–1)、经向风速 (m·s–1)、垂直风切变等参数计算实验区域内的积冰指数和Dutton指数. 并根据公式（2）、（4）、（8），得到地图中各格点的雷暴概率P1、积冰概率P2、颠簸概率P3.由于所选高度为民航飞机一般巡航高度，该高度层内几乎无积冰现象产生，故P2=0. 根据公式（9），计算得到危险天气的混合概率，即地图风险值，并将危险天气的混合概率大于50%的区域设置为障碍物区，如图5所示.

图5 实验地区危险天气区域图. (a)风险地图; (b)障碍物区域（P>50%）Fig.5 Dangerous weather region in the experimental area: (a) risk map; (b) obstacle region (P>50%)

3.2 飞行路径规划

基于3.1节中危险天气区域图，设定左上角为起点，右下角为终点，采用文章第2部分提到的风险最小化A*算法和Dijkstra算法进行路径规划，并与经典路径最短条件下A*算法、Dijkstra算法和RRT算法进行对比，分别计算五者的路径长度和总风险值. 其中RRT算法重复1000次实验，选取最优解. 路径的总风险值为路线所遍历格点的风险值之和，计算结果如表5所示，路径结果如图6所示.

图6 不同算法下的路径规划图Fig.6 Path planning graphs for different algorithms

从表3中可以看出，在传统的A*和Dijkstra算法下，所得路径最短，为148.79单位长度，但由于计算方式不同，所得路径有所差异，因此总风险值不同，且数值均偏大；采用风险最小化A*算法得到的路径总风险值最小，为0，但同时也需要更大的路径代价；采用风险最小化Dijkstra算法得到的路径与最短路径相比，略有增加，但总风险值大大减少；RRT算法具有较强的随机性，因此在有限时间内不保证结果收敛到最优结果，根据以往的研究经验，该方法在高维、复杂、动态的环境中，更能发挥其优势，在处理本问题的过程中未表现出明显优势.

由此可知，①当只考虑风险最小时，应选择风险最小化A*算法下的飞行路径，此时路径安全性最高，无遭遇危险天气风险；②当综合考虑路径长度和风险时，可选择风险最小化Dijkstra算法下的飞行路径，此时路径长度和风险值均较小，且与前者相比该方案经济性更高，可节约燃油成本和时间成本.

4 总结

为解决飞行计划层面危险天气下的路径规划问题，提出了一种基于危险天气不确定性的最小风险路径规划方法，经计算讨论得到结论如下：

（1）C-F模型可以用于融合雷暴诊断中出现的多种不确定性.

（2）考虑雷暴、积冰、颠簸等多种航路危险天气，建立了基于危险天气概率的风险地图，用于航前飞行路径规划.

（3）风险最小化A*算法和Dijkstra算法可以在不同程度上以增加路径长度为代价，降低飞行风险；RRT算法更适用于高维复杂环境下的路径规划，在解决本文所述问题时无明显优势.

本文仅给出了判定危险天气的一种方法，具体参数如1.1节中导致雷暴发生的预报要素阈值ε1～ε4和1.4节中危险天气概率的判定阈值δ，还需考虑区域、季节、安全裕度等因素进行相应调整. 此外，由于相对湿度和经/纬向风高空数值预报数据的缺失，本文使用了再分析数据替代的方法来弥补数据的不足. 也希望相关气象单位能够补充数值预报、集合预报的三维数据，特别是高空预报数据，以满足航空气象业务和科研需求.