考虑降水时间相关性的地面观测-雷达-卫星遥感逐时降水融合方法研究

阮惠华，张钧民 ，许剑辉，戴晓爱

(1.广东省气象探测数据中心，广东 广州 510641；2.广东省科学院广州地理研究所/广东省遥感与地理信息应用重点实验室/广东省地理空间信息技术与应用公共实验室，广东 广州 510070；3.成都理工大学地球科学学院，四川 成都 610059)

1 引 言

随着气象观测系统的迅猛发展，通过地面自动气象站、雷达、卫星遥感等手段获取的不同时空分辨率的降水观测数据越来越多[1]；各行业对高质量的时空连续降水数据产品要求越来越高，特别是针对极端天气情况下高时空分辨率的降水数据，对洪涝灾害及滑坡监测等具有极其重要的意义[2-4]。

综合气象站-雷达-卫星遥感-地表辅助参量等多种资料，多源数据融合与数据同化技术已成为获取高精度、高质量、时空连续的多时空尺度降水融合网格数据产品的有效手段[1,5]。根据地面观测-雷达、地面观测-卫星遥感、雷达-卫星遥感、地面观测-雷达-卫星遥感等不同源的降水数据组合，相关学者开发了包括地理加权回归方法(Geographically Weighted Regression，GWR)[6-7]、回归克里金插值方法[8-9]、贝叶斯模式平均方法(Bayesian Model Averaging，BMA)[10-12]、最 优 插值[13-14]、多尺度模式融合方法[15]、随机森林与XGBoost 机器学习方法[16-19]等不同的多源降水数据融合模型，并在子午河盆地、淮河流域、扬子江流域、华南地区等区域进行应用分析。为获取更高时间分辨率的降水数据产品，相关学者利用“概率密度函数+贝叶斯模型平均(BMA)+最优插值”组合方法有效地融合地面观测降水-雷达定量降水估测产品-卫星遥感降水产品，研制了高精度的1 km 分辨率的逐时降水融合产品[1,20-21]，并对这些逐时降水融合产品在四川、黄土高原、长江流域、华南地区等区域的适用性进行对比评估[22-24]。结果表明，在华南区域，融合降水产品与实际降水量的偏差较小，对实际降水量的再现能力强。

尽管现有的逐时降水数据融合研究取得一定进展，针对龙舟水、台风水等不同类型的暴雨过程，逐时多源降水融合方法的适用性仍需进一步研究。此外，现有的逐时降水融合方法较少考虑暴雨过程降水的时间相关性。受到热带与中纬度天气系统的共同影响，华南地区成为我国暴雨最频繁的地区之一。基于此，本文利用广东省北部山区2018 年汛期强降水时段4 月23—28 日(北京时间，下同)龙舟水、5 月7—11 日龙舟水和9 月16—17 日台风“山竹”3 次暴雨过程250 个气象站逐时降水数据，以及雷达降水、卫星遥感降水、地形与海岸线距离等辅助变量，在对不同类型的逐时降水数据进行质量控制处理基础上，基于机器学习算法与地统计学理论，开展空间分辨率为1 km逐时气象站点-雷达-卫星遥感降水的融合试验，充分考虑相邻时刻降水的时间相关性，以期获取复杂地形区域暴雨过程高精度的降水融合数据，并对逐时降水融合试验结果进行分析评估。

2 数据及预处理

2.1 研究区概况

研究区主要位于粤北山区(图1)。研究区海拔最低6 m，最高1 427 m，高海拔区域主要集中在东西两边，中间区域海拔整体相对较低，主要集中在清远东南部的英德、清新，清城境内的北江河谷。研究区属亚热带季风气候，雨水资源丰富，平均年降水量1 631.4～2 149.3 mm，年平均降水日(日降水量≥0.1 mm/d)为160～173 天，是广东省一个典型的“雨窝”[25]。

图1 研究区域及气象站点分布 a.广东省边界范围；b.研究区域及气象站点分布。

2.2 研究数据

本文以广东省气象局提供的逐时降水数据作为地面基准数据，250个气象站点的空间分布如图1a 所示。随机选择空间均匀分布的200 个气象站点(占80%)作为训练集站点，50 个气象站点(占20%)为测试集站点并用来评价不同降水融合模型的精度(图1b)。选取2018 年4 月23—27 日(过程I)、5 月7—10 日(过程II)、9 月16—17 日(过程III)3次暴雨过程的降水数据进行逐时降水数据融合试验，相关数据介绍如表1所示。

表1 相关数据来源及介绍

雷达降水：由1 km 分辨率的逐6 min 雷达定量降水估测产品计算得到，具体的计算过程参考文献[25]。

与海岸线距离、DEM：这两个地表辅助参量的处理过程主要参考文献[25]。

IMERG 和GSMaP 降 水：IMERG(Integrated multi-satellite retrievals for GPM) 和 GSMaP(Global satellite mapping of precipitation)是新一代全球降雨观测计划(Global precipitation measurement, GPM)广泛应用的卫星降水产品，具有覆盖范围更广(扩展到全球)、时空分辨率更高的优势(IMERG：空间分辨率为0.1 °、时间分辨率为0.5 h；GSMaP：空间分辨率为0.1 °、时间分辨率为1 h)[26-28]。本文利用面到点克里金插值方法[27]分别对IMERG 和GSMaP 降水产品进行空间降尺度，得到空间分辨率为1 km 的逐时IMERG 和GSMaP降水产品。

3 方 法

3.1 考虑降水时间相关性的XGBoost 逐时降水数据融合方法

本文提出的考虑时间相关性的XGBoost逐时降水融合方法技术路线如图2所示。

图2 考虑时间相关性的多源逐时降水融合方法技术流程图

(1) 利用普通克里金插值算法对前两个时刻的地面气象站逐时降水进行插值，得到空间分辨率为1 km的逐时网格降水数据(PreSiteOK)；前两个时刻的降水插值结果(PreSiteOK)作为降水时间相关性引入到XGBoost模型中。

(2) 采用地统计学的面到点克里金插值(Area-to-point Kriging，ATPK)[27-29]分别对空间分辨率为0.1 °的逐时IMERG 和GSMaP 降水产品进行空间降尺度。面到点克里金插值方法是利用已知面对未知点进行插值估计的空间降尺度方法，其原理与普通克里金类似，即未知点值为其所在面以及附近面数据的线性加权求和[27]：

权重λx的求解如下：

式中，x为高空间分辨率(1 km)的待插值网格点，K为低空间分辨率(0.1 °)的网格数目，vi为低空间分辨率的网格，表示面与点协方差函数，表示面与面协方差函数为空间分辨率为0.1 °网格的IMERG 和GSMaP 降水产品,λx为权重，μx为拉格朗日乘子,为降尺度的空间分辨率为1 km的IMERG和GSMaP降水产品。

(3) XGBoost 是一种基于从GBDT(梯度提升树)算法改进和扩展而来的经过优化的集成学习算法[30]，本文拟采用XGBoost 算法建立当前时刻气象站-前两个时刻网格插值-雷达-卫星遥感多源降水预测模型：

通过步骤(3)获取每时刻的XGBoost降水融合模型以及对应的逐时降水残差。将空间分辨率为1 km 的逐时雷达降水数据、前两个时刻网格降水数据、空间降尺度的卫星遥感降水产品和地表辅助变量作为输入，利用构建的XGBoost 逐时降水融合模型预测1 km分辨率的逐时降水。

(4) 对XGBoost 降水融合模型的逐时降水残差进行普通克里金插值，得到1 km 分辨率的逐时降水残差，并与模型预测结果进行相加得到1 km分辨率的高质量逐时降水融合结果。

3.2 融合试验设计

为了进一步验证不同逐时降水融合模型的性能，设计了“考虑降水时间相关性的XGBoost逐时降水数据融合(本文提出方法，记为方案I)”、“不考虑降水时间相关性的XGBoost逐时降水数据融合(XGBoost 方法，记为方案II)”、“不考虑降水时间相关性的随机森林逐时降水数据融合(RF方法，记为方案III)”三组逐时降水融合对比试验方案。具体的试验设计为：方案I，利用普通克里金插值算法对前两个时刻的地面气象站逐时降水进行插值处理，并作为辅助变量，采用XGBoost算法进行逐时气象站点降水-前两个时刻逐时网格降水-雷达降水-卫星遥感降水融合分析；方案II，不考虑前两个时刻逐时网格降水，采用XGBoost 算法进行逐时气象站点降水-雷达降水-卫星遥感降水融合分析；方案III，不考虑前两个时刻逐时网格降水，采用RF算法进行逐时气象站点降水-雷达降水-卫星遥感降水融合分析。

3.3 评价方法

结合检验站点的降水观测数据，采用主要的评价指标对逐时降水融合结果进行定量评价，包括均方根误差(RMSE)、平均绝对误差(MAE)和决定系数(R2)，其计算公式如下：

式中，n表示用于检验的站点数，Ok表示第k个气象站点降水观测数据，表示第k个气象站点对应的逐时降水融合结果，Oˉ表示气象站点降水观测数据的平均值。

4 结果与讨论

结合研究区的50个地面站点逐小时降水观测数据，分别利用决定系数(R2)、均方根误差(RMSE)和平均绝对误差(MAE)3 个统计指标评价2018 年4 月23—28 日(龙舟水)、5 月7—11 日(龙舟水)和9月16—17 日(台风“山竹”)3 次暴雨过程的小时尺度多源降水数据融合结果的精度。

图3 对应2018 年3 次暴雨过程的方案I、方案II和方案III逐时降水融合结果。结果显示，3次暴雨过程中，方案I融合结果比方案II和方案III融合结果的精度高，R2均大于0.609。9 月暴雨过程方案I 融合结果的RMSE 和MAE 最小，分别为2.198和1.065 mm；方案I 降水融合的R2达到了0.78，高于方案II和方案III降水融合的R2。这可能是因为方案I不仅考虑了当前时刻地面降水与雷达降水、卫星遥感降水、地表辅助变量间的非线性关系，还将前两个时刻地面降水作为自变量引入到模型中，考虑了降水数据的时间相关性。此外，图3 也显示了3 次暴雨过程的方案II 逐时降水融合都优于方案III 逐时降水融合，3 次暴雨过程的最大逐时降水量逐渐增加，暴雨过程II 的最大逐时降水量超过50 mm；暴雨过程I的逐时降水融合结果比过程II 和过程III 的逐时降水融合结果的精度低。结果表明，在降水量相对较大的情况下，三种方案都能获取较高精度的逐时降水融合数据；然而，在降水量相对较小的情况下，三种方案的降水融合结果精度相对较低。这可能是因为在降水量相对较小的情况下，卫星遥感降水产品的估测存在不足[31]，逐时降水融合模型引入卫星遥感降水产品导致降水融合结果存在一定的误差。

图3 2018年过程I(a)、过程II(b)和过程III(c)三种方案的逐时降水融合与站点降水观测比较

图4 显示了3 次暴雨过程不同融合方法的逐时降水融合的R2、RMSE 和MAE 箱图。对比3 次暴雨过程的逐时降水融合，方案I 的3 次暴雨过程降水融合的中位数R2分别为0.415、0.466和0.734，都高于其他两种方法降水融合的中位数R2；除了5月暴雨过程方案I 降水融合的中位数RMSE 略大于方案II 降水融合的中位数RMSE，方案I 的暴雨过程I 和III 降水融合中位数RMSE 都小于其他两种方法降水融合的中位数RMSE，暴雨过程I和III降水融合的中位数RMSE 分别为0.418 mm 和1.735 mm；方案I 的3 次暴雨过程降水融合的中位数MAE 都高于其他两种方法降水融合的中位MAE。结果表明，方案I降水融合结果优于方案II和方案III降水融合结果。这进一步表明在逐时降水融合过程中引入多个相邻时刻的降水信息能有效提高降水融合模型的精度。

图4 不同融合方法的逐时降水融合的R2(a)、RMSE(b)和MAE(c)箱图

比较方案II 和方案III 降水融合发现，暴雨过程I 和过程II 的方案II 降水融合的中位数R2小于方案III 降水融合的中位数R2，但3 次暴雨过程方案II 降水融合的R2均值都大于方案III 降水融合的R2均值；3次暴雨过程的方案II降水融合的中位数RMSE 和MAE 整体上小于方案III 降水融合的中位数RMSE 和MAE。这表明，方案II 降水融合的精度整体上优于方案III降水融合。从图4 也可看出，方案I、方案II 和方案III 三种融合方案在过程III 都能获取较高精度的逐时降水融合结果，其中位数R2分别达到了0.734、0.722和0.643，远高于暴雨过程I和过程II逐时降水融合的中位数R2。

表2 和图5～7 分别显示了3 次暴雨过程方案I、方案II和方案III的R2表现最好的时刻占总暴雨过程比例情况、逐时降水融合的R2、RMSE和MAE时间序列以及对应时间50个地面站点逐时降水观测的平均值。从表2 可看出，3 次暴雨过程方案I获得更高精度的逐时降水融合，优于方案II 和方案III 降水融合，R2占比全部大于53.82%，远高于其他两种降水融合方案。

表2 不同降水融合方案R2表现最好的时刻占总暴雨过程的比例情况

图5 暴雨过程I逐时降水融合R2(a)、RMSE(b)、MAE(c)精度评价

从图5 看出，暴雨过程I 的方案I 降水融合整体上优于方案II 和方案III 降水融合(表2)。但也出现方案I 降水融合精度比方案II 和方案III 降水融合精度低的情况，比如暴雨过程I 的23 日23 时方案I 降水融合的R2远小于方案II 和方案III 降水融合的R2。这可能是因为暴雨过程I的23 日21—22 时降水量比较大，而到了23 时降水量显著减少，出现有的区域可能没有降水的情况，方案I 引入21—22时两个时刻相对较大降水信息反而增加降水融合误差。从图5也可看出，尽管方案II降水融合精度整体上稍微优于方案III 降水融合，但也有不少时刻出现方案II降水融合精度差于方案III降水融合。

从图6看出，暴雨过程II的方案I降水融合R2、RMSE 和MAE 整体上优于方案II 和方案III 降水融合，特别是在持续降水量比较大的暴雨过程II的9 日02 时—10 日12 时这段时间内，方案I 降水融合精度明显优于方案II和方案III降水融合。在这段时间内，部分方案II 降水融合的R2低于方案III 降水融合的R2，但这两种方法的RMSE 和MAE相差不大。

图6 暴雨过程II逐时降水融合R2(a)、RMSE(b)、MAE(c)精度评价

图7 显示了暴雨过程III 的方案I、方案II 和方案III 降水融合的R2、RMSE 和MAE 时间序列，尽管出现方案II和方案III降水融合的R2略优于方案I降水融合的R2，如暴雨过程III的16日10时；但整体上，方案I 降水融合的R2、RMSE 和MAE 明显优于方案II 和方案III 降水融合的R2(表2)、RMSE 和MAE。这表明在这次暴雨过程中，考虑前两个时刻降水信息能更好地改善逐时降水融合模型精度。此外，图7 也显示了，方案II 降水融合精度整体上优于方案III 降水融合精度，因为方案II 降水融合具有更高的R2和更低的RMSE和MAE。

图7 暴雨过程III的逐时降水融合R2(a)、RMSE(b)、MAE(c)精度评价

图8 显示了2018 年3 次暴雨过程不同融合方法预测的累计降水空间分布，在3次暴雨过程方案I、方案II与方案III方法的降水融合结果具有相似的空间特征。

图8 不同融合方法的累计降水融合结果空间分布图 a.暴雨过程I；b.暴雨过程II；c.暴雨过程III。

暴雨过程I的累积最大降水量不超过110 mm，暴雨过程II 的累积最大降水量增加到180 mm，到了暴雨过程III，累积最大降水量超过了210 mm。在累积降水量比较高的区域，方案II 降水融合结果的高值范围整体上比方案I 和方案III 降水融合结果的高值范围大；在累积降水量比较低的区域，方案II 降水融合结果的低值范围整体上比方案I和方案III 降水融合结果的低值范围大。此外，方案III 在极大与极小降水融合的性能表现不太好，在降水量相对较高和较低的区域，方案III 降水融合不能较好地刻画降水分布的空间细节。总体上来看，方案I 降水融合结果获得较满意的结果，具有更多的空间分布特征。

总的来说，相比方案II 和方案III 逐时降水融合模型，方案I 获得较高精度的逐时降水融合结果，逐时降水融合具有更丰富的空间分布特征。然而，方案I 逐时降水融合结果仍存在一定的误差，这可能与引入的空间分辨率为0.1 °的IMERG和GsMaP卫星遥感降水产品有关。粗分辨率的卫星遥感降水产品在高值和低值的降水估测仍存在一定的误差[31-33]。未来的研究应该引入更多时间序列暴雨过程的逐时降水信息，通过长短期记忆人工神经网络(LSTM)挖掘暴雨过程降水的时间特征信息[34]，并以此优化逐时降水融合模型，进一步提高逐时降水融合精度。此外，结合遥感数据获取的云量、水汽等辅助参数[35]，在充分挖掘降水量与云量、水汽之间时空特征信息基础上，引入到逐时降水数据融合模型中提高多源降水数据融合模型的精度。

5 结 论

本文以广东省北部山区2018年汛期强降水时段4 月23—28 日龙舟水、5 月7—11 日龙舟水和9月16—17 日台风“山竹”3 次典型暴雨过程的逐时降水为研究对象，建立基于XGBoost 与克里金插值算法的地面观测-雷达-卫星遥感多源逐时降水融合模型，充分考虑相邻时刻降水的时间相关性，开展了空间分辨率为1 km 的逐时降水融合试验，并与不考虑降水时间相关性的XGBoost和随机森林(RF)算法的逐时降水融合模型进行对比，得到如下结论。

(1) 考虑降水时间相关性的方案I逐时降水融合模型融合结果精度整体上高于方案II和方案III逐时降水融合结果，与气象站逐时降水数据最接近，因为方案I通过引入前两个时刻逐时网格降水数据充分考虑了相邻时刻降水的时间相关性，使得降水融合结果能够更好地刻画空间细节。在出现降水量极大和极小值的区域，不考虑降水时间相关性的方案II和方案III逐时降水融合精度有所降低。

(2) 与方案II 和方案III 逐时降水融合结果相比，方案I逐时降水融合结果在不同暴雨过程的准确性均有明显改进，3次暴雨过程的RMSE分别降低了6.0%和9.7%、6.3%和9.5%、10.5%和30.0%。

(3) 方案II 的精度整体上优于方案III，3 次暴雨过程的RMSE 分别降低了3.9%、3.4%和21.8%，表明XGBoost算法在刻画逐时气象站降水与雷达降水、卫星遥感降水、地表辅助参量间的非线性关系能力方面比RF算法更有优势。