1. 引言
降水预报是天气预报中的重点和难点，精准的降水预报对于防灾减灾、工农业生产和公众生活等方面具有重要意义。以华北地区为例，“63.8”海河流域暴雨、“75.8”河南暴雨、“96.8”河北暴雨均是发生在华北地区的著名暴雨过程。随着城市化进程加快，以京津冀为代表的城市群极端暴雨天气频发，如2004年7月10日北京突发性局地暴雨以及2011年6月23日、2012年7月21日北京特大暴雨天气过程，特别是2012年7月21日北京特大暴雨，降水强度大，持续时间长，造成北京市79人死亡，河北省26人死亡、20人失踪，受灾人口达到190万，直接经济损失200多亿元。因此加强降水预报研究，努力提高暴雨预报水平对减小生命财产损失有重大意义。目前降水的预报主要依赖于数值模式，其每天在世界各地主要气象业务中心运行，在灾害性天气精细化预报中发挥着不可替代的作用，但其同时具有不完整性和不精准性的特点，尤其在具有复杂下垫面条件时和面向强天气影响时预报能力明显不足。机器学习可以方便地进行空间和时间的概括，处理大量的观测数据和预报变量，并能学习和集成物理理解到模型中。在气象领域，已有大量研究使用机器学习方法来提高天气预报的准确性，如应用单个数值预报模式的多个气象要素来构建天气预报模型 [1] [2]。门晓磊等 [3] 基于中国气象局、欧洲、日本气象厅和NCEP的四个数值的近地面2 m处温度预报和实况，分别使用岭回归、随机森林和深度学习对2 m处温度进行了订正，这三种方法的预报效果都优于四个数值天气预报模式以及通常的多模式集合平均。在强对流预警方面，Czernecki等 [4] 通过随机森林方法也实现了对冰雹灾害的精准估计；Hyangsun等 [5] 在决策树、随机森林、支持向量机等机器学习技术的基础上，利用各种气象、卫星数据，实现了强对流形成过程的检测。降水订正方面，陈浩等 [6] 选择K近邻(KNN)、多元自回归样条(MARS)、支持向量机、多项对数线性模型(MLM)和人工神经网络(ANN)等方法，利用海拔、相对湿度、坡向、植被、风速、气温和坡度等因子订正了ITPCAS (Institute of Tibetan Plateau Research, Chinese Academy of Sciences)和CMORPH (Climate prediction center morphing technique)两种常用的青藏高原日降水数据集。张晨阳等 [7] 将基于地形的加权随机森林(terrain-based weighted random forests, TWRF)方法用于雷达定量降水估计，并在此基础上考虑了反射率垂直廓线的特征重要性变化，有效提升了雷达反演降水效果。黄威等 [8] 采用支持向量机回归方法，基于ECMWF的51个集合预报和NCEP的21个集合预报，建立了多模式集成的强降水预报模型。机器学习应用于专业领域的众多研究中特征提取具有关键作用 [9] [10] [11] [12]。本文通过集成机器学习中的特征工程方法，综合挖掘多源观测和数值预报等气象信息，遴选面向强降水短临预报的关键要素和时次，以加深对小时尺度强降水特征的物理理解，为强降水发生发展的机理研究和短临预报机器学习订正提供客观支撑。
2. 方法和数据
2.1. 特征工程方法概述
特征工程包括数据清洗、数据预处理、和特征选择等过程。特征工程的目的在于获得更容易表现数据性质的训练集。一般特征越好，灵活性越强，构建的模型越简单，模型的性能越出色，构建模型的计算代价越低。常用的特征工程可分为三大类，过滤式，包裹式和嵌入式(如图1所示)。其中过滤式通过设定阈值的评价准则来判识特征与目标量的关系，这种方式可理解性强，计算效率高，但由于与分类器剥离，其对弱相关特征和冗余度的判断并不理想。过滤式结合了特征和分类器算法，搜索选择预测精度较好的特征子集实现改进，但非常耗时。嵌入式直接在机器学习算法中通过权值和贡献度提取重要特征，常能获得比前两类更好的效果，如树模型和梯度提升模型(如XGBoost)。本文采用了经典的随机森林模型进行研究，并在3.3节与相关系数、互信息等特征工程进行集成分析。
Figure 1. Schematic diagram of feature engineering classification
图1. 常用特征工程分类示意图
2.2. 随机森林模型及特征重要性度量
随机森林(Random Forest，以下简称RF)是一种能够进行回归和分类的机器学习方法，由多个分类或者回归树组成(如图2所示)。RF在以决策树为基学习器构建Bagging集成的基础上，还在决策树的训练过程中引入了随机属性选择。具体来说，对基决策树的每个结点，先从该结点的属性集合中随机选择一个包含k个属性的子集，然后再从这个子集中选择一个最优属性用于划分。随机森林的基学习器的多样性不仅来自样本扰动，还来自属性扰动，因此使得最终集成的泛化性能可通过个体学习器之间的差异的增加进一步提升。随机森林简单、容易实现、计算开销小，且在很多现实任务中展现出强大的性能，被誉为“代表集成学习技术水平的方法”。
随机森林能够较好地度量每个特征的重要性。通过计算每个特征在随机森林中的每棵树上所做的贡献，取平均，就可以比较出不同特征的贡献大小。通常用基尼指数(Gini index)或者袋外数据(OOB)错误率作为评价指标来衡量。目前python中的sk-learn库已经实现了用随机森林评估特征重要性，训练好模型后，可直接调用feature_importances属性得到每个特征的重要性。
Figure 2. Schematic diagram of Random Forest Model
图2. 随机森林模型示意图
2.3. 数据集及样本平衡
本文使用的多源数据包括加密气象观测资料、高分辨率全球环流模式预报和精细化中尺度数值预报产品资料(以欧洲中期天气预报中心(以下简称EC)和华东区域模式(以下简称SMS)预报产品为例)。选取地面观测降水为目标量(因变量)，温、压、湿、风等大气要素和综合诊断的动力热力特征量为自变量，进行多源数据的时空匹配融合，构建强降水短临预报物理要素大数据集。根据天气分析预报经验和特征完备性原则，将众多物理要素归为5类特征，具体为：模式预报的目标量(降水预报，Type 0)、水汽特征(Type 1)、动力抬升特征(Type 2)、热力不稳定特征(Type 3)和维持时间特征(Type 4)。根据天气系统连续和浑沌共存的特性，同时也选取一定历史长度的气象要素作为特征，经前期试验结果选取历史6小时以内的历史观测、模式预报和目标时刻的模式预报作为特征量，预测当前时刻的累计降水。5类特征要素和时次组合共计产生376个特征。
由于强降水属于小样本事件，为获取更多的强降水样本，特选取北方4次典型的强降水事件所在日构建数据集，在个例选取是注意囊括了系统性强降水、局地短时强降水和台风降水等不同降水类型，具体采样覆盖时间为20180804~20180807，20190802~20190804，20190810~20190812。图3给出了采样时段内降水分布，降水中心的累积降水超过500 mm，京津冀地区的降水中心量值超过250 mm。基于此样本构建特征X，分别以1小时、3小时和6小时累计降水(以下简称r1、r3和r6)为目标Y，随机选取部分时次共计129万左右的样本进行特征分析，考虑到样本不均衡(即无降水和低量级降水样本数量远远多于大量级降水)，因此需要对其进行适当的下采样，得到如表1所示的样本分布。将采样后的样本放到随机森林模型中进行回归训练，获取特征重要性得分。
Table 1. Sample distribution of 3 h rainfall before and after sampling along with thresholds
表1. r3样本采样前后不同降水量阈值样本分布
Figure 3. Distribution of accumulated rainfall of selected heavy rain events
图3. 选取强降水事件的累积降水量分布
3. 试验分析
3.1. 特征提取及精简
对原始376个特征进行分类贡献度排名，分别给出以1小时/3小时/6小时累计降水为目标时5个类型特征中得分排名前十名的特征评分值(见表2和图4)。
Figure 4. Feature scores of top ten important features in five types with 6 h (orange)、3 h (green) and 1 h (blue) precipitation
图4. 原始特征情况下排名前十特征的分类分值(橙色：以r6为目标量；绿色：以r3为目标量；蓝色：以r1为目标量)
Table 2. Feature scores of top ten important features
表2. 排名前十特征评分
图4中显示，随着降水预报精细程度提升(6小时/3小时/1小时累积降水为目标量)，模式直接预报的降水(type 0)贡献度急剧下降，而相关物理特征，包括水汽特征(type 1)和热力不稳定特征(type 3)重要性显著增强。这说明，在越来越精细的小时尺度强降水短临预报中，模式直接预报的降水可信度迅速减弱，相反，物理量预报(尤其是水汽要素)对降水订正有更好的指导意义。从6 h累计降水到1 h累计降水，模式直接预报降水贡献越来越小，与大气的多尺度特征相对应，也反应了数值模式对多尺度系统的可预报性随精细程度降低。因此，应用模式预报效果较好的物理特征量订正直接的降水预报具有较大的实际意义。
同时，从得分直方图(图略)分布中可以看出，每个类型的前10特征中，“马太效应”较为明显，即排名前三的特征得分之和占据了这类型特征得分总和的50%以上，这意味着，在不牺牲太大准确率的情况下，可以较为放心的舍弃掉较多的特征，以此换得较多的计算便利。
3.2. 特征冗余剔除及分析结果
3.2.1. 剔除冗余
原始376个特征，虽然保持了很好的完备性，但存在冗余。过多的冗余特征会带来如下麻烦：1) 冗余特征使得机器学习模型更难感知到真正有益于目标识别的重要特征，因为障碍物变多了；2) 过多冗余特征使得模型不容易学习到最优的结果；3) 增大了计算成本，模型变得复杂。因此很有必要进行冗余特征的剔除和重要特征的选择。这里选择嵌入式特征选择方法，基学习器采用随机森林。嵌入式的特征选择过程流程是一个完整的模型训练过程。
Step 1：准备数据集。从所有个例时刻中随机选取十分之一的时刻的样本构建数据集。
Step 2：考虑到不同时次目标累计降水下，重要特征可能存在差异(后面的试验也验证了这一点)。分别以r1、r2和r3为目标，训练3个模型。训练之前需要做降水量级样本平衡(因为无降水样本始终占大多数)，方法为随机剔除90%的无降水样本。通过随机森林给出的feature_importance可以定量分析出不同目标下各个特征的得分情况。
Step 3：将所有特征得分按照5个类型(以下简写为type)进行排序。将r1，r2和r3中的type 1中排名前15的特征求并集，其他类型依次类推，这样共计得到5个type的106个重要特征。
通过上述步骤既保证了特征的完备性，也极大降低了特征的冗余程度。根据精简的后的5大类106个特征，可以更为仔细的分析出各个类型特征的重要性得分情况。步骤如下：
① 样本选取。从295个时刻的样本集按时间顺序每间距一选取一个时次，选出144个时次的样本集。每个时次的样本集中包含n_samples个样本。如果n_sample个样本中r3 (目标3小时累计降水) > = 1 mm的样本数量小于n_sample/50，则认为该时刻不存在有效降水。去除该样本集。通过上述步骤，将剩下的m个时刻的样本集的所有样本合并成一个数据集。
② 样本平衡。以3小时累计降水(r3)为目标，根据上述步骤获取的样本数和不同降水量的样本分布情况，可以看到非降水样本占据了80%，另外由于之前站点到格点的插值过程会产生部分0~0.1之间的样本，因此需要做样本平衡。通过随机采样1/10的无降水样本和r3在0~0.1 mm之间的样本，即可以得到相对平衡的样本分布(如表3所示)。
③ 模型训练。分别以r1，r2，r3，r6为目标通过随机森林进行回归学习训练。
Table 3. Sample distribution with different thresholds after sample balance
表3. 样本平衡处理后的不同降水量样本分布
3.2.2. 分析结果
在保持特征物理完备性的前提下，按照上述步骤进行冗余特征的精简。得到共计106个特征下，分别以r1，r2，r3，r6为目标通过随机森林获取特征量得分统计(图略)。显示结果虽然与图2略有不同，但总体信息类似：随着预报精细化程度提升(以r6/r3/r1为预报目标)，模式直接预报的降水贡献度急剧减小，而其他物理特征重要性显著增强，除了与图2类似水汽特征和热力不稳定特征贡献增大，特征精简后动力抬升特征的作用在1 h降水预报尤其突出，维持时间特征重要性也显示出增强。为具体分析高影响的特征因子，以3 h和1 h降水预报为例，给出了分类特征的前十因子排序直方图(图3和图4)。
特征重要性排序显示，物理要素比时效特征的作用更大，整体而言越临近时刻的要素预报与目标降水量的相关性越强。直接预报降水特征(Type 0，图略)显示，EC模式预报的降水量特征评分远高于其他特征，其次是地面观测的前一时刻降水和SMS中尺度模式预报的降水。以3小时累积降水为目标量(如图5所示)，在水汽特征中(Type 1)，EC预报的整层可降水量贡献明显高于其他水汽特征，其次是模式预报的低层水汽含量和水汽通量散度。动力抬升特征中(Type 2)，中低层的风场(低层辐合)作用显著，但中低层的整层抬升也是值得考虑的因子。在热力特征中(Type 3)，EC模式预报的对流不稳定条件和SMS中尺度模式预报的对流有效位能和反射率因子重要性凸显。
Figure 5. The top ten important physical factors in the five types feature for 3 h precipitation
图5. 3小时累计降水分类特征贡献前十的物理因子
Figure 6. As the same of Figure 5, but for 1 h precipitation
图6. 同图5，但为1小时累计降水分类特征
以1小时累积降水为目标量进行特征重要性分析(如图6所示)，在水汽特征中(Type 1)，EC预报的整层可降水量贡献明显高于其他水汽特征，其次是模式预报的低层水汽含量，与r3不同，当精细到逐小时降水时，与降水强度(云微物理)相关的特征因子，如零度层高度(Type1_Z0)的重要性位列第三。动力抬升特征中(Type 2)，中底层的整层抬升和近地面层的作用显著，和r3相比，低层涡度场也是值得考虑的因子。在热力特征中(Type 3)，地面观测的热力条件和中尺度模式预报的反射率因子贡献最大，其次是EC模式预报的对流有效位能和不稳定参数(Li指数和K指数)。总体而言，随着精细化程度提升，热力特征的重要性调整较大，其次是水汽特征。因此，在精细化的强降水预报机器学习订正中应着重考虑这两类特征的提取。
3.3. 特征工程集成提取
为说明上述特征提取的通用性，本节采用特征工程里最常用的5种方法进行特征提取集成，包括皮尔逊相关系数法、互信息方法(Mutual Information)、逻辑回归(Logistic Regression, LR)和梯度提升法(GBDT)，并与随机森林(RF)结果进行对比。选取各个类型的排名前15的特征，对其特征得分进行绝对值平均，如表4所示。
Table 4. Importance of the five types of physical features by five methods
表4. 五种特征工程方法给出的五类物理特征的重要度
表4中对每种特征工程方法提取出的特征中最重要的两类特征进行了加粗标注，可以看到不同的提取方法虽然在具体评分上有所差异，但整体提取的特征信息是一致的，在以r6为目标量时，模式直接预报的降水可信度很高，是最重要的参考特征，但随着精细程度的提高，其作用急剧下降，水汽特征和动力特征重要性增强，尤其到1小时降水预报，多种方法均显示出热力特征的重要性凸显。
4. 总结
近年来，随着人工智能技术的快速发展，结合物理模型和机器学习的专业领域应用得到越来越多关注。本文利用机器学习对大数据处理的优势，通过特征工程客观提取了不同尺度强降水预报的关键影响因子，获得主要结论如下：
1) 强降水受到大气运动的高维特征影响，通过特征工程客观筛选和先验专家知识判别，精简后的特征可有效去除模型冗余，在不损失机器学习精度和效果的情况下大力提升运行效率。
2) 定时、定点、定量的强降水预报是防灾减灾的重难点之一，当预报较长时间累计降水时(r6以上)，模式直接预报降水具有较好的参考价值，但当精细到3 h以内时，降水预报可信度急剧减小，而相关物理特征更具参考价值。具有共性的是越接近预报时刻，物理要素的参考价值越大。不同的是对于较长时段的累计降水(r6)，中高层的动力维持作用、整层抬升和整层含水量具有较大贡献，而随着预报精细程度提高(r1)，低层和近地面层的动力作用、热力不稳定特征更加重要，影响降水强度的云微物理相关因子(如零度层高度)作用也更加凸显。
3) 从6 h累计降水到1 h累计，模式直接预报降水贡献越来越小，而相关物理特征如水汽特征，不稳定条件等因子对降水预报具有更大的贡献。这与大气的多尺度特征相对应，也反应了数值模式对多尺度系统的可预报性随精细程度的提升而降低。因此，通过特征工程提取的重要物理特征为提升强降水预报提供了重要支撑。
基金项目
本研究由国家重点研发计划(2018YFC1507606)，国家自然科学基金项目(42030611，42075154，41505079)资助。
NOTES
*通讯作者。