梯度提升决策树在雷达定量降水估测中的应用研究

2019-11-15陈训来陈元昭赵春阳张珂

Advances in Meteorological Science and Technology 2019年3期

陈训来陈元昭赵春阳张珂

（1 深圳市气象局，深圳 518040；2 深圳南方强天气研究重点实验室，深圳 518040；3 阿里云计算有限公司，杭州 310099）

0 引言

定量降水估测（QPE）是现代天气业务体系的重要组成部分，一直以来是业务预报的难点和重点，也是定量降水预报、强降水短临预警的基础，对城市内涝、地质灾害精细化预报起重要支撑作用。随着新一代多普勒天气雷达观测网的逐步建成和业务运行，应用多普勒天气雷达联合自动气象站雨量开展定量降水估测已成为一种较为成熟可靠的实现手段[1]。

利用多普勒天气雷达定量降水估测，相对于用自动气象站雨量观测，具有测量范围广、时空分辨率高、可及时取得大面积定量降水资料的优点。由于雷达回波（反射率因子）反映的是空中水滴、冰晶等对雷达电磁波的反射或散射，因此雷达回波强度（Z）与降水强度（I）之间存在密切的指数型关系，即Z=A×Ib[2]。通过Z与I之间的对应关系就可以估测出雷达覆盖范围内的降水强度。降水定量估测精准度在很大程度上取决于Z-I关系式中系数A和b，降雨常具有时空特性强、变化快的特点，且不同季节、不同地点的降水性质不同，其Z-I关系式也不同，常用来建立Z-I关系的方法有最优化法[3]、概率密度法[4-5]、动态关系法[6]、遗传算法[7]和动态分级法[8]等。为了得到更高精度的降水估计值，通常利用自动气象站雨量资料对雷达降水估测值做进一步订正，包括应用平均校准法[9]、变分同化法[10]、客观插值分析法[11]和卡尔曼滤波校准法[12]等。

随着人工神经网络（ANN）等深度学习技术的发展，很多专家学者利用深度学习的算法开展雷达定量降水估测的研究，由于深度学习是一种非参数运算，且具有自适应学习方法，能有效减少误差，提高了降水估测的准确性[13-16]。Liu等[13]研究发现基于人工神经网络的定量降水估测比用Z-I关系更精确。Chiang等[14]利用台湾六次台风降雨过程研究表明，采用动态人工神经网络比Z-I关系有更高的精度和稳定性。邵月红等[15]通过BP神经网络方法估测定量降雨，精度要优于Z-I关系式，估测的降雨量与站点实测雨量吻合性较好。傅德胜等[16]基于径向基函数神经网络（RBF）建立的雷达定量估测降水模型，其对降水估测精度和稳定性高于以变分法得到的Z-I关系。

目前针对华南地区开展基于深度学习的定量降水估测算法的研究还比较少，本文主要应用三维网格雷达回波资料和自动气象站数据开展基于梯度提升决策树（GBDT）算法在雷达定量估测降水中的研究和应用，提高雷达定量降水估测的精度。

1 数据和预处理

1.1 数据

本研究使用的雷达资料为广东省境内广州、深圳、韶关、清远、阳江、河源、汕尾、梅州、湛江和肇庆等10部S波段多普勒雷达在VCP21探测模式下，每隔6 min对9个仰角进行三维体扫描，为确保所有雷达能同步扫描，采用了北京敏视达公司的雷达同步控制和数据采集系统时间同步软件进行时间同步。为计算方便，在雷达三维体扫描原始数据基础上，利用Cressman插值法将极坐标格式的雷达探测资料插值到三维直角坐标系中，并采用了兰红平等[17]提出的“膨胀-侵蚀”算法对回波缺陷进行质量控制。由于在雷达组网拼图时，若低于1.5 km，任意两部雷达之间几乎没有共同的覆盖区域，观测的范围太小；若高于5 km，雷达覆盖的范围扩大了，然而却难以形成降水[18]。因此，本研究中用到经过雷达回波质量控制和拼图处理后得到1.5、2.5、3.5和4.5 km共四层高度上CAPPI（constant altitude plan position indicator）数据，空间分辨率为0.01°×0.01°，水平网格点为约为500 km×300 km，数据覆盖为2015—2017年的每年4—9月的降水过程。

另外，本研究中还用到该范围内区域自动气象站逐小时雨量观测资料，自动气象站分布如图1。

图1 自动气象站分布图（黑点） Fig. 1 Distribution of automatic weather stations(AWSs) (black dots)

1.2 预处理

由于雷达回波为空中水滴、冰晶等对雷达电磁波的反射或散射，自动气象站记录的是落到地面的水滴，在选取雷达的回波强度数据时，需要综合考虑多种因素的影响。梁建茵等[19]认为，对低纬度地区，自动气象站点上降水强度与其上空约5 km×5 km区域内的平均雷达回波强度关系较好，使用平均雷达回波强度的方法可以在一定程度上减少降水估测的误差。因此，为了更有效的利用三维雷达回波强度数据进行定量降水估测，本文建立以自动气象站为中心、1 h内逐6 min的5 km×5 km格点范围的雷达回波强度数据（图2）。研究中用到的自动气象站降水资料为整点的1 h雨量，与雷达资料的匹配方式为：T时刻记录的过去1 h雨量（T-1到T的1 h累计雨量）与T-1到T的逐6 min、10个时次5 km×5 km格点的雷达回波反射率相匹配，形成一条数据样本，作为降雨定量估测模型的输入数据。本文共建立100万组训练数据集和10万组测试数据集。

2 基于GBDT的雷达定量估测降水算法

梯度提升决策树是一个应用广泛、效果较好的监督式机器学习模型，它于2001年由Friedman[20]提出的组合决策树算法，由决策树（decision tree）[21]和梯度提升（gradient boosting）[22]组合而成，提高回归预测的精度。

GBDT是一种迭代的决策树算法，其迭代的目的是通过计算上一次模型的负梯度来改进模型，每一次新模型的建立都减少了上一个模型的残差，然后在残差减少的梯度方向上建立新的决策树；每个新模型的输入都是上一个模型输出，输出是下一个模型的输入，预测模型的最终结果为之前所有模型结果的相加。该算法由多个弱预测模型（通常是决策树）通过集合的形式组合产生预测模型，所有树的结论加起来作为最终答案。该方法的优点在于可以尽量避免过拟合现象，每迭代一棵树，残差代表的都是全局最优的绝对方向。由于GBDT算法考虑方差和偏差之间的权衡，具有较好的泛化能力，几乎可以解决所有的数据回归问题并具有高检测精度的特性，且在非线性数据处理方面具有良好的效果，被广泛应用消息检测[23]、径流预测[24]等研究领域。

本文将每一个数据样本是具有一组特征属性的雷达回波三维网格数据和自动气象站小时雨量数据组成，定义训练数据集：

图2 以自动气象站为中心的雷达回波反射率数据 Fig. 2 Three dimensional data of radar reflectivity for GBDT as input (The gauge location is at the center of the grid)

假设通过递归迭代分割将特征空间（即输入变量）划分为m个子集，则GBDT模型表示为m个模型函数相加来预测输出：

基于GBDT的定量降水估测QPE模型算法的详细步骤如下：

1）使用常量值初始化模型，估计使损失函数极小化的常数值，它是只有一个根节点的树，即β是一个常数值：

迭代训练模型，在上一次模型损失函数的梯度下降方向上建立模型，从m=1到M（M为迭代次数）

计算损失函数的负梯度在当前模型的值，将它作为残差ri的估计值：

2）将式（5）中得到的估计残差ri作为输入，拟合到一颗回归树，求得回归树的叶节点区

3 个例分析

利用建立的训练数据集进行GBDT算法的建模和训练，得到基于GBDT算法的雷达定量降水估测模型，并用测试数据集进行检验和误差分析。作为比较，本文还利用测试数据集建立了固定Z-I关系法和动态Z-I关系法进行定量降水估测。其中固定Z-I关系法应用汪瑛等[8]通过统计历史上强降水日的雷达反射率因子和1 h降水资料，得到适合广东地区的系数A和b的值，即Z=220.8I2.23，这个统计结果可直接用于雷达定量降水估测，且计算速度快。动态Z-I关系法主要利用自动气象站雨量计降水资料，基于Z-I关系算法，实时计算当前时刻最优系数A和b的值，使雷达定量估测降水与区域自动站观测降水的误差最小，作为下一时次的Z-I关系计算多普勒雷达估测降水，该方法的优点是前期不需要收集大量的资料样本做统计准备，只需要统计上一个小时的最优Z-I关系。

图3给出了2017年4月21日07时的影响广东区域一条飑线雷达定量估测降水和06—07时的1 h雨量实况。该小时降水的主要特点是有一条东北—西南向覆盖广东中部地区、小时雨量大于10 mm的降雨带，在广州和茂名存在两个强降水中心（小时雨量在40 mm以上），其中最大小时雨量出现在广州，为65.3 mm（图3d）。图3a为固定Z-I关系法的定量降水估测，与自动气象站实况雨量相比，降水分布和落区较为一致，但是降水量级估测偏小很多，对两个强降水中心的降水估测较实况比偏小20～30 mm。动态Z-I关系法对降水的估测比固定Z-I关系法有所改善（图3b），对位于广州的强降水中心估测在40 mm以上，虽然仍比实况偏小，但在降水量级上比固定Z-I关系法已有10 mm的提高。图3c是基于GBDT算法反演的雷达定量估测降水，不同等级降水量的面积和强度相对固定Z-I关系法和动态Z-I关系法都有所改进，尤其对位于广州的强降水中心的降水量在50～55 mm，与实况更接近，改善了降水量偏小的问题。

4 雷达定量降水结果和误差分析

为了描述定量估测降水的误差和精度特征，在检验和分析时，使用平均误差（ME）、相对误差（RE）和均方根误差（RMSE）等进行统计评价，具体计算公式如下：

图3 2017年4月21日07时雷达1 h定量降水结果和自动气象站实测降水（a）固定Z-I关系法；（b）动态Z-I关系法；（c）GBDT算法；（d）雨量实况 Fig. 3 Radar-derived 1h QPE and gauge spatial distributions on 06:00 BT 21 April 2017(a) 1 h QPE of Z-I relations; (b) 1 h QPE of dynamic Z-I relations; (c)1 h QPE of GBDT scheme; (d) observation

式中，n为总样本数，i为QPE雨量计的序列，Pi和Gi分别为雷达估测雨量和自动气象站实测雨量。

用平均误差评估数据的整体差异，其值越接近0则数据整体差异越小等评价方法；使用相对误差评估数据的可信度，其值越小数据可信度越高；用均方根差评估雷达估测雨量与自动气象站观测雨量的离散程度，其值越小分布越集中；从整体上评估不同方法对降水估测精度的高低，所有误差都在站点上计算。为了使评估的结果具有客观性，剔掉没有通过质量控制的自动气象站观测雨量数据，将1 h雨量划分为1～5、5～10、10～20、20～30、30～40、40～50 mm和50 mm以上的7个降水量级，分析这三种方法对降水的反演效果。

图4 给出了这三种方法的雷达定量降水和自动气象站观测雨量散点图。越接近直线点表明雷达定量降水的效果越好，高于直线说明高估，低于直线点为低估。可以看出，固定Z-I关系法对雷达定量降水的反演存在低估，随着降水量级的增大，低估程度趋于明显，与实测降水的相关系数为0.71；动态Z-I关系法对雷达定量降水的估测准确性有所改进，雷达估测降水和实测降水的相关系数略有提高，为0.74；应用GBDT算法定量估测的降水与实况观测相关性较好，相关系数提高到0.80，散点图显示集中性较好，改善了降水估测的精度，但是仍然存在低量级降水高估、高量级降水低估的现象。

表1是对不同量级降水的雷达定量降水误差分析，可以看出，固定Z-I关系法对10 mm/h以下的降水存在高估，而对10 mm/h以上的降水明显低估，对50 mm/h以上的强降水低估平均误差达到34.14 mm，相对误差为56.95%。动态Z-I关系法对雷达定量降水估测效果优于固定Z-I关系法，各量级的平均误差和相对误差较固定Z-I关系法有所减小，但仍然存在对10 mm/h以下量级降水高估，而10 mm/h以上量级降水低估的现象。这可能是由于动态Z-I关系法虽然根据自动站观测雨量和雷达回波关系不断调整A和b值，提高雷达定量降水估测能力，但仍然存在对同一时刻的降水类型没有进行分类，所有降水类型都用同一个A、b值的Z-I关系的不足。

图4 实测雨量和雷达定量估测降水散点图：（a）固定Z-I关系法；（b）动态Z-I关系法；（c）GBDT算法 Fig. 4 The scatter diagram of hourly rainfall estimation: (a)Z-I relations; (b) dynamic Z-I relations; (c) GBDT scheme

基于GBDT算法反演的降水平均误差和相对误差都小于固定Z-I关系法和动态Z-I关系法的结果，尤其对30 mm/h以上量级降水的估测，相对误差比动态Z-I关系法都提高了5%以上，对降水反演的系统误差更小，改善了降水估测的精度。RMSE的含义是雷达估测雨量与自动气象观测的降雨量之间的离散程度，对5～20 mm/h量级的降水估测，动态Z-I关系法RMSE最小，但对于20 mm/h以上的强降水，GBDT算法估测降雨量的RMSE都小于固定Z-I关系法和动态Z-I关系法，表明GBDT算法对强降水的雨量估测拟合度较高。