基于模糊逻辑的飞机积冰预测指数

齐 晨 金晨曦 郭文利* 甘 璐 赵德龙 卢 姁 武 帅 李黑平

1)(北京市气象服务中心， 北京 100089)2)(北京市人工影响天气办公室， 北京 100089)3)(中国人民解放军32021部队气象室， 北京 100094)4)(中国人民解放军95820部队气象台， 北京 102207)

引 言

飞机积冰是由于飞机表面迎风部位碰到云中过冷却水滴或降水中的过冷却雨滴后，在机身表面某些部位凝结成冰的现象[1]。积冰使飞机的重量和阻力增加、升力和推力减小，导致飞机的飞行效率降低，并会影响飞机的发动机性能、无线电通讯和各种仪表的正常运行，因此飞机积冰严重威胁飞行安全。

飞机积冰预测是飞机积冰研究中的重要组成部分，采用数值预报模式结合积冰诊断算法判断积冰发生情况已成为积冰预测的主流方法。目前国内常使用国际民航组织推荐的飞机积冰指数Ic[2]判断积冰的发生和强度，国外常见积冰诊断算法主要有用于积冰类型预测的RAP(Research Applications Program)算法[3]、通过探空资料判断积冰强度和类型的RAOB(Rawinsonde Observation)算法[4]以及NAWAU(National Aviation Weather Advisory Unit)积冰预报方案[5]。近年国内外学者开展了大量积冰诊断算法研究，王洪芳等[6]使用MM5(The Fifth-Generation NCAR/Penn State Mesoscale Model)中尺度数值模式对多种积冰算法进行了比较，并建立了飞机积冰预报模型；刘开宇等[7]使用世界区域预报系统(World Area Forecast System,WAFS)数值预报资料和常规数据对积冰指数进行了修正；何新党等[8]建立了基于云微物理参数的多因子积冰预测模型；McDonough等[9]基于模糊逻辑和积冰情景决策树分类方法建立了潜在积冰预报算法(the For- ecast Icing Potential Algorithm,FIP)。当前国内常用积冰判别算法以积冰相关气象条件的阈值判断的组合为核心，但缺乏关于接近阈值的气象条件对积冰影响的描述，也缺乏由非积冰环境过渡到积冰环境的过程描述,而近期逐渐发展的飞机积冰集成预报方法虽然在机理性上具有优势，但对研究样本数据和数值预报模式的输出场要求较高，难以结合区域模式进行本地化。

为建立阈值更加合理且易于本地化的飞机积冰判断指标，本研究基于模糊逻辑隶属度函数，通过对飞机积冰和非积冰个例对应气象条件的分析，定义了以气温和相对湿度为判别基础、并考虑垂直运动和过冷却液态水含量影响的积冰指数Ip(icing potential index)，综合多种气象要素判断空中有利于积冰发生的区域，为飞行员规避积冰高风险区域提供参考。

1 飞行报告积冰记录和对应气象要素信息

1.1 北京市人工影响天气办公室飞行记录积冰个例和对应机载观测数据

北京市人工影响天气办公室(简称北京人影办)飞行记录积冰个例包含2014—2017年京津冀地区共计372个积冰个例发生的日期、时间、经纬度、海拔高度和对应气象要素信息，积冰个例发生位置如图1所示。积冰个例对应气象要素信息通过机载观测仪器获取，包含积冰现象发生时段内的气温、相对湿度和垂直风速值。

图1 2014—2017年北京人影办积冰个例分布Fig.1 Distribution of icing cases from 2014 to 2017 provided by Beijing Weather Modification Office

1.2 全国空中报告积冰个例

2016年全国空中报告(PIREPs)记录了当年全国范围内飞机发生颠簸、风切变和积冰事件的日期、时间、经纬度、海拔高度和强度信息。本文采用该报告中位置确定的91个积冰记录作为积冰个例样本，同时选取了70个位置确定的颠簸或风切变记录作为非积冰个例样本，用于积冰指数Ip的建立和检验。积冰和非积冰个例的发生位置如图2所示，个例对应最临近时刻和位置的气象条件根据ERA5再分析数据确定。

1.3 ERA5再分析数据

ERA5再分析数据是欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，ECMWF)发展的第5代全球大气再分析数据(https:∥cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels)，结合模式与观测数据提供了近期气候的定量化描述。研究表明：相对于美国国家环境预报中心(National Centers for Environmental Prediction,NCEP)和日本再分析资料数据库(Japanese Re-Analysis,JRA)再分析数据，ERA数据在积冰预测方面具有准确性更高的优势[10]。相较于第4代ERA-Interim再分析数据，ERA5在时间、空间分辨率和准确性等多方面均有明显提升。本研究使用0.125°×0.125°空间分辨率的逐小时ERA5再分析数据，结合2016年全国空中报告积冰和非积冰个例的时间和位置信息，获取了91个积冰个例样本以及70个非积冰个例样本发生的最邻近时刻、位置的气温、相对湿度、垂直速度、云量和云中液态水含量信息。

图2 2016年全国空中报告积冰个例(a)和非积冰个例(b)分布Fig.2 Distribution of icing cases(a) and non-icing cases(b) from 2016 national pilot reports(PIREPs)

2 积冰指数Ip的建立

云中存在的过冷却水滴或降水中存在的过冷却雨滴(即过冷却液态水)是导致飞机发生积冰的关键，通过分析云中或降水中存在过冷却液态水的有利气象条件，判断飞机在空中出现积冰现象的可能性，建立积冰指数Ip。本研究基于模糊逻辑原理，使用北京人影办飞行记录积冰个例样本和对应机载观测气象数据建立气温和相对湿度的隶属度函数(Tmap和Rmap)，以此判断积冰事件发生的初始可能性(Ipi),在此基础上选取2016年全国空中报告中部分飞机积冰和非积冰个例样本，通过多种组合筛选，确定了表示垂直运动和过冷却液态水含量对初始积冰可能性影响的隶属度函数(Vmap和Smap)以及二者的权重系数(k1和k2)，最终计算得出积冰指数Ip，并使用2016年全国空中报告中剩余的飞机积冰和非积冰个例样本检验该指数的性能。积冰指数Ip的具体计算方法如下：

①以气温和相对湿度数据作为输入，根据气温隶属度函数(Tmap)和相对湿度隶属度函数(Rmap)计算积冰事件发生的初始可能性(Ipi)，其计算公式为

Ipi=Tmap×Rmap。

(1)

②当Ipi计算结果等于0时，积冰指数Ip等于0。

③当Ipi计算结果大于0时，在判别垂直运动方向的基础之上，分别采用不同公式计算垂直运动和过冷却液态水含量对初始积冰可能的影响：当垂直速度小于或等于0时，积冰指数Ip的计算方法为

Ip=Ipi+(1-Ipi)×(k1×Vmap+k2×Smap)；

(2)

当垂直速度大于0时，积冰指数Ip的计算方法为

Ip=Ipi+Ipi×k1×Vmap+(1-Ipi)×k2×Smap。

(3)

2.1 模糊逻辑隶属度函数

20世纪60年代Zadeh[11]使用了一种基于模糊集和模糊逻辑的方法来处理语言学中的模糊信息。此后这一方法不断发展，逐步形成了模糊逻辑系统理论，并广泛应用于气象学、医学、计算机科学和系统科学等多个领域。针对气象学领域，近年来模糊逻辑系统理论在强天气识别[12]、强对流天气监测预警[13]、雷达质量控制[14]和回波识别[15]等方面应用较多。

模糊逻辑理论可用于表达界限不清晰的定性知识和经验。基于隶属度函数的概念，使用模糊逻辑可以模仿人的不确定性判断和推理方式，描述并区分模糊集合，处理模糊关系，从而对经典逻辑难以处理的非线性、不确定问题进行定量化描述[16]。

在经典集合理论中，可以使用隶属度法定义集合，即引入等价于集合的0～1隶属度函数。为了克服经典集合理论无法描述不具有清晰边界集合的局限性，模糊集合概念应运而生。它允许隶属度函数在区间[0, 1]内取任意值描述边界模糊的集合。

本研究使用模糊逻辑隶属度函数描述气象要素对积冰可能性的非线性影响，其中气温和相对湿度的隶属度函数通过对观测数据的统计结果的拟合得出，而描述垂直运动和过冷却液态水含量对积冰可能性影响的隶属度函数则通过筛选不同隶属度函数的最优组合确定。

2.2 气温和相对湿度的隶属度函数Tmap，Rmap的确定

大量研究表明：过冷却液态水主要存在于云层或降水中，当飞机发生积冰时云量基本上是多云或者阴天[17]。同时有研究表明：云中或降水中的过冷却液态水更多地出现在环境温度接近冰点时，当气温进一步降低时由于冰晶含量增加，过冷却液态水将减少[18-20]。根据积冰个例样本对应气温和相对湿度分布规律建立相应的模糊逻辑隶属度函数，分别用于描述气温与过冷却液态水存在可能性之间以及相对湿度与云或降水存在可能性之间的相关性，综合二者判断云层或降水中存在过冷却液态水的可能性，以此判断环境的初始积冰可能性。

2.2.1 气温的隶属度函数Tmap的确定

从-35～0℃以1℃为间隔设置35个气温区间，根据北京人影办372个飞机积冰个例样本对应气温，统计每个气温区间内包含的个例样本量，如图3所示。

由图3可以看到，积冰事件对应气温主要在-15～0℃区间内，当气温低于-31℃或高于0℃时不存在积冰现象。在-15～0℃气温区间内，积冰个例样本量存在明显的先增后降趋势，并且超过半数的积冰事件发生在-8～-4℃气温区间内，而气温在-4～0℃对应的积冰个例样本量迅速减少，这与机身表面空气压缩加热有关[9]。根据图3所示的积冰个例样本量对应气温的分布规律，建立气温的模糊逻辑隶属度函数(Tmap)，用于表示不同气温条件下过冷却液态水存在的可能性，其具体形式如图4所示。

图3 各个气温区间对应积冰个例Fig.3 Sample size of icing cases in each temperature interval

图4 气温的隶属度函数Fig.4 Membership function of temperature(Tmap)

气温的隶属度函数描述了气温与过冷却液态水存在的可能性之间的相关性：当气温为-8～-4℃时Tmap达到峰值，表示当环境气温处于该气温区间时过冷却液态水存在的可能性最高；在该区间左侧，随着气温降低，云或降水中的冰晶含量逐渐升高，过冷却液态水存在的可能性逐渐降低；在该区间右侧，由于飞机表面空气的压缩加热作用，过冷却液态水存在的可能性迅速降低。

2.2.2 相对湿度的隶属度函数Rmap的确定

类似地，可确定相对湿度的隶属度函数，从0～100%间隔5%，共设置20个相对湿度区间，根据北京人影办372个积冰个例对应相对湿度，统计每个区间包含的积冰个例样本量(图5)。

图5 各相对湿度区间对应积冰个例Fig.5 Sample size of icing cases in each relative humidity interval

相对湿度的隶属度函数描述了相对湿度和云层或降水存在的可能性之间的相关性，反映了相对湿度对过冷却液态水存在的可能性的影响：当相对湿度低于30%时，很难有云层或降水存在，因此过冷却液态水存在的可能性也很低；随着相对湿度增加，云层或降水存在的可能性增加；当相对湿度大于95%时，Rmap达到峰值，表示云层或降水存在的可能性很高，环境有利于过冷却液态水存在。

综合气温和相对湿度的隶属度函数，可以判断环境的初始积冰可能性,当Tmap和Rmap同时达到峰值时，表示云层或降水中存在过冷却液态水的可能性很高。

图6 相对湿度的隶属度函数Fig.6 Membership function of relative humidity(Rmap)

2.3 垂直运动、过冷却液态水含量的隶属度函数Vmap，Smap和二者权重系数的确定

过冷却液态水的存在是飞机发生积冰现象的关键。当云层存在时，上升运动有助于液态水的产生，而下沉运动可能表示液态水的减少[21]。目前部分数值预报模式可以输出过冷却液态水含量的预测值，但将其直接应用于飞机积冰预测存在较大不确定性，有研究指出，积冰现象很可能出现在模式预测存在过冷却液态水的情况下，但模式预测过冷却液态水不存在，不能表示不发生积冰现象[9,21]。同时有研究表明，在构建积冰指标过程中考虑过冷却液态水含量的影响有助于降低虚警率[22]。基于上述前人研究成果，在气温和相对湿度的隶属度函数判别基础上，综合考虑垂直运动和过冷却液态水含量对于环境初始积冰可能性的影响(假设二者对初始积冰可能性的影响相互独立)，构建积冰指数Ip。

受使用数据限制，本研究结合观测数据、前人研究结果和经验，通过筛选多种表示垂直运动和过冷却液态水含量对初始积冰可能性影响的隶属度函数以及不同权重系数的最优组合，最终确定积冰指数Ip的计算方法。

2.3.1 筛选组合方案

本研究使用垂直速度表示垂直运动的强度和方向，并设置了3种不同的垂直速度隶属度函数形式Va，Vb和Vc(图7)。

由北京人影办积冰个例对应的垂直速度观测数据分析可知：25%的积冰个例不存在明显的垂直运动；66.4%的积冰个例对应垂直上升运动，而8.6%的积冰个例对应垂直下沉运动，且其中绝大多数积冰个例对应垂直速度为-2～2 m·s-1。这表明大多数积冰事件的发生都伴随有垂直上升运动，但较弱的下沉运动对积冰的抑制作用不明确，因此建立Va，Vb，Vc3种形式的垂直运动隶属度函数，其中Va形式为CIP(Current Icing Potential)算法中使用的表示垂直速度对积冰初始可能性影响的隶属度函数[21]，Vb和Vc分别表示较弱的下沉运动对积冰发生无影响和有促进作用。

ERA5再分析数据中云中，液态水含量(specific cloud liquid water content)和云量(cloud cover)可用于表征云中是否存在过冷却液态水,结合前人研究成果与经验综合确定[22]。Sa和Sb所采用的隶属度函数形式如图8所示。

不该禁的地方坚决不禁：各个地方要科学的划定禁养区，防止盲目扩大养殖范围，禁养区必须按照标准来划定，避免只要环境不要生产的极端做法。

图7 垂直速度的隶属度函数Fig.7 Membership functions of vertical velocity

图8 云中液态水含量Sa和云量Sb的隶属度函数Fig.8 Membership functions of specific cloud liquid water content Sa and cloud cover Sb

为确定垂直运动和过冷却液态水含量的隶属度函数在积冰指数Ip计算过程中的权重，共设置7组不同权重系数组合(表1)，用于控制垂直运动和过冷却液态水含量对于初始积冰可能性的影响程度，其中k1和k2分别表示垂直运动项和过冷却液态水含量项在积冰指数计算过程中的权重系数。

2.3.2 筛选标准与结果

随机选取2016年全国空中报告中30个积冰个例样本和25个非积冰个例样本，根据上述组合方案，对于每个积冰个例可计算输出42种组合结果。对每种组合中积冰和非积冰个例样本的积冰指数Ip计算结果，分别统计准确率A、漏报率M、虚警率F和高值率H(高值率为Ip计算结果大于0.6的样本量在积冰个例样本量中所占比例)。

表1 垂直运动项权重系数(k1)和过冷却液态水含量项权重系数(k2)的组合Table 1 Combination of weight coefficients k1 and k2

根据分类算法评估指标，假设积冰个例样本量为P，非积冰个例样本量为N，使用积冰指数Ip进行判别后：积冰个例判别为存在积冰可能性(Ip>0.1)的样本量为Tp，非积冰个例判别为不存在积冰可能性(Ip≤0.1)的样本量为Tn，积冰个例判别为不存在积冰可能性(Ip≤0.1)的样本量为Fn，非积冰个例判别为存在积冰可能性(Ip>0.1)的样本量为Fp，积冰个例判别为存在较高积冰可能性(Ip>0.6)的样本量为Q。则准确率

(4)

漏报率

(5)

虚警率

(6)

高值率

(7)

由表2可以看到，当垂直速度隶属度函数(Vmap)选择Va形式、过冷却液态水含量隶属度函数(Smap)选择Sb形式、权重系数选择k1=0.6,k2=0.4的组合时，积冰指数Ip对积冰事件的判别准确率达到最高值，漏报率和虚警率同时达到低值，同时高值率也较高，因此采用该种组合方式建立积冰指数Ip，对积冰事件的判别效果最佳。

表2 筛选结果Table 2 Screening results

综上所述，积冰指数Ip的计算步骤如下：首先以气温和相对湿度作为输入，根据其隶属度函数(Tmap和Rmap)计算得出初始积冰可能性(Ipi)；若初始积冰可能性为0，则积冰指数Ip为0。对于初始积冰可能性大于0的情况，根据垂直速度的判别，使用相应表达式和隶属度函数(Vmap和Smap)计算垂直速度和云量对于初始积冰可能性的影响，最终得出积冰指数Ip。通过对全国30个积冰个例样本和25个非积冰个例样本的分析，确定垂直运动项权重系数k1=0.6，过冷却液态水含量项权重系数k2=0.4。

根据上述方法计算得出的积冰指数Ip可以反映在一定的气象条件下，飞机在空中发生积冰现象的可能性。积冰指数Ip的计算结果与积冰可能性大小的对应关系如表3所示。

表3 积冰指数Ip的值对应的积冰可能性Table 3 Possibility of icing corresponding to the value of icing potential index

3 积冰指数Ip的检验

将上述积冰指数Ip的计算方法，用于2016年全国空中报告中未参与组合筛选的剩余61个积冰个例样本和45个非积冰个例样本，计算积冰指数Ip对积冰事件判别的准确率、漏报率和虚警率，以验证该指数的性能。

计算结果显示，积冰指数Ip判别的准确率为80.2%，其漏报率和虚警率分别为9.4%和10.4%。目前国内常使用国际民航组织推荐的飞机积冰指数Ic[2]判断积冰的发生和强度，该指数的构建为

Ic=[(R-50)×2]×[T×(T+

14)/(-49)]，

(8)

式(8)中，R为相对湿度，T为温度。使用相同的个例进行计算，结果显示其准确率、漏报率和虚警率分别为57.5%，19.8%和22.6%。这表明积冰指数Ip对积冰事件发生的判别准确率更高，且漏报率和虚警率显著降低，因此积冰指数Ip的判别效果优于积冰指数Ic。

4 结论与讨论

本文基于模糊逻辑理论，通过对积冰个例对应气象条件的分析，定义了基于气温、相对湿度、垂直速度和云量的模糊逻辑隶属度函数的积冰指数Ip，用于判断飞机在空中出现积冰现象的可能性，并使用积冰和非积冰个例样本对该积冰指数Ip进行检验，结果表明：

1) 积冰指数Ip对积冰事件的判别准确率为80.2%，漏报率和虚警率分别为9.4%和10.4%。

2) 与目前国内常用的经典积冰指数Ic相比，积冰指数Ip的判别准确率显著提升，漏报率与虚警率均显著降低，整体判别效果显著提升。

3) 计算积冰指数Ip所使用的物理量以大气温湿层结数据为主，同时也是数值预报模式常见输出量，因此该指数易于根据区域观测资料进行本地化调整，并结合高分辨率区域数值模式产品为飞机躲避积冰高风险区域提供更加精细化的参考。

需要注意的是，本文在构建积冰指数Ip的过程中使用了两组来源不同的积冰个例样本，其一为2016年全国空中报告积冰个例样本，其二为2014—2017年北京人影办积冰个例样本。由于飞行目的、计划和路线的显著差异，导致两组积冰个例样本之间存在明显的时间和三维空间分布差异。但通过对2016年全国空中报告积冰个例样本对应气象条件的分析发现，其对应气温、相对湿度和垂直速度分别表现出与北京人影办积冰个例样本对应气象要素相似的统计特征。因此两组积冰个例样本之间存在的时间、空间分布差异，对本研究中气温、相对湿度和垂直速度的隶属度函数的建立影响较小。

还需要注意的是，本文假设积冰指数Ip的计算过程中，垂直运动项和过冷却液态水含量项对初始积冰可能性的影响是独立的，并设置了互补的权重系数。垂直运动和过冷却液态水含量之间的关系本文也未讨论，有待进一步探究。

另外，本文并未讨论积冰个例样本之间由于飞机机型和飞行速度等因素不一致导致的积冰气象条件差异，不同机型在不同飞行速度下动力增温的差异会影响气温相关的隶属度函数的函数曲线[23]，因此在实际应用中应根据机型和飞行速度对气温相关曲线进行调整。

随着机载观测气象数据质量的改善[24]和数值预报系统对航空气象数据同化能力的提升[25]，以及数值预报模式对于降水云类型识别[26-27]和液态水含量预报能力的提升，积冰指数Ip在飞机积冰预测中的参考价值将进一步提升，但同时也需要大量多源观测数据对其进行验证与调整，从而更加准确判别积冰事件的发生。