利用空间随机森林方法提升GPM卫星遥感降水质量

胡保健，李伟，陈传法，胡占占

1.山东科技大学 测绘与空间信息学院,青岛 266590;

2.湖州中核勘测规划设计有限公司,湖州 313000

1 引 言

降水是地球能量循环与物质交换重要组成部分，是陆地生态系统重要驱动因素，在气候变化中起着决定性作用（Markonis 等，2019；Zhou 等，2020b）。高分辨率和高精度的空间降水分布对区域水资源调度管理、灾害预测防控、生态研究等领域具有重要的应用价值（Brodeur和Steinschneider，2020；Renard 等，2011；Zhang 等，2020）。然而，降水数据估算误差在上述应用中会带来很大的不确定性（Xie 和Xiong，2011）。如何提高降水空间分辨率和准确性是近年来生态水文等领域研究热点和难点（Shen等，2014b）。

降水监测数据主要来自气象观测站、地基雷达测量、卫星遥感降水反演产品等。气象站观测降水可直接测得实际降水量，但受地形限制，难以反映降水的空间连续性分布（Sharifi 等，2019；Shen 等，2014a）。地基雷达观测可以提供连续高分辨率空间降水数据，但只能反映局部区域降水，并且复杂地形区域布设困难、易受环境因素干扰多。随着遥感技术的发展，产生了大量基于星载探测技术的卫星遥感降水产品，如TRMM、PERSIANN、GPM、CHIRPS、CMORPH（Hu等，2014；Islam等，2020；Sun 等，2016）。相对于气象站点观测与雷达测量降水，卫星遥感降水产品具有覆盖范围广、时空连续变化效果好的特点（Zhou 等，2020a），在水文模型和生态环境研究中得到广泛应用（Jiang等，2021；Wei等，2020）。

然而，现有卫星遥感降水产品分辨率较低（最优仅为0.05°），而且受传感器误差与降水反演算法不确定性等影响，现有降水产品均含有一定的系统偏差，导致其不能满足局部尺度和生态流域尺度的应用研究（Immerzeel 等，2009）。因此，国内外研究者们相继提出了多种点面融合方法，以借助地面观测站降水对卫星降水产品修正，如最优插值法（OI）（卢新玉 等，2017；潘旸 等，2012）、贝叶斯模型（Fu等，2016）、地理差异分析（GDA）（Jongjin等，2016）、地理加权回归（GWR）（Lu等，2019）、K最近邻模型（KNN）（王玉丹 等，2016）等。这些方法虽然有效提高了卫星降水产品的精度，但仍有一些问题需要解决：（1）融合后的降水产品分辨率与原有的卫星降水产品相同，不利于局部地区水文模型研究；（2）一些方法只考虑了实测降水背景场，未考虑与降水有关的环境因素；（3）在进行偏差校正时，粗分辨率的卫星降水产品与基于点的观测站实测降水存在尺度差异，导致点面融合时捕捉不到局部降水细节信息。因此，在点面融合之前对遥感降水产品进行空间降尺度是十分必要的。

近年来，国内外学者们提出和发展了多种遥感降水产品降尺度方法，如指数回归（ER）（Immerzeel等，2009）、多元线性回归（MLR）（嵇涛 等，2015；Jia 等，2011；马金辉 等，2013）、GWR 方法等（胡实 等，2020）。随着研究深入，一些研究者们开始尝试使用机器学习方法对遥感降水进行降尺度研究。大量研究表明，随机森林RF（Random Forest）模型在处理多维数据中具有良好的表现，且能够很好的捕捉自变量与因变量之间非线性关系（Baez-Villanueva 等，2020；Njuki 等，2020），具有良好的应用前景（Jing 等，2016；Ma 等，2018）。然而，上述降尺度研究均将RF 模型作为一种简单的统计工具，忽略了相邻遥感（站点）降水数据的空间相关性。

基于以上讨论，本文提出了一种顾及空间相关性的双阶段卫星降水产品质量提升方法，即降尺度与点面融合相结合的思路，以提升卫星降水产品的空间分辨率与精度。该方法以RF 为基础模型，构造了顾及空间相关性的空间随机森林模型SRF（Spatial Random Forest）。与现有方法相比，新方法优势包括：（1）考虑了降水之间空间相关性；（2）降尺度与点面融合两阶段都采用了SRF模型，且均融入了环境变量影响。以四川省为研究区域，选取最新GPM IMERG V06B 遥感降水数据作为数据源，使用SRF-SRF 方法对GPM 降水数据进行降尺度与点面融合处理，并将计算结果与经典方法比较，验证新方法的可行性和高效性。

2 研究区域和数据源

2.1 研究区概况

四川省位于中国西南部（97°21′E—108°31′E，26°03′N—34°19′N），总面积约为4.86×105km2，介于青藏高原与长江中下游，是两者的过渡地带。研究区域地形地貌复杂多样，涵盖了山脉、高原、丘陵、平原和盆地等地形，地形起伏走势西高东低。受地形因素和空间位置影响，该区域降水夏季多冬季少，其中夏季降水占全年80%—90%，且受川西垂直地势起伏影响，四川盆地降水充沛，年降水量800—1600 mm，呈现由东到西降水逐渐减少的变化特点（Yang 等，2020；Lai 和Gong，2017）。图1为研究区域地形及气象站点分布。

图1 四川省气象站分布Fig.1 Distribution of weather stations in Sichuan Province

2.2 数据源

（1）气象站数据。气象站观测数据由中国气象局气象数据中心（http://data.cma.cn/［2021-01-01］）提供。选取四川省内2015 年—2019 年期间156 个逐日气象站点降水数据（图1），所有数据经过内部一致性检查、极值检查以及空间一致性检查等严格质量控制后（Song 等，2004），经累积计算得到每月降水量。

（2）卫星降水数据。GPM 降水数据从美国宇航局地球科学数据中心（https://pmm.nasa.gov/［2021-01-01］）获取，是美国NASA（美国国家航空航天局）和日本JAXA（日本宇宙航空研发机构）共同研发的继TRMM 之后全球降水测量卫星，搭载全球先进的DPR（Dual-frequency Precipitation Radar）和GPM GMI（GPM Microwave Imager）系统，对微小降雨探测能力比前一代产品TRMM 更精准，覆盖范围（60°S—60°N）更广阔。本文选用月尺度0.1°×0.1°分辨率的降水产品GPM IMERG V06B（简称GPM）数据集作为本研究所采用的数据源，时间跨度为2015年1月—2019年12月。

（3）DEM、NDVI 以及地表温度（LST）数据。研究采用的DEM 数据为SRTM（Shuttle Radar Topography Mission）DEM V4.1 数据集，从地理空间数据云（http://www.gscloud.cn/［2021-01-01］）获取，空间分辨率为90 m，经像素平均重采样生成1 km 分辨率，然后在Arcgis 中利用表面分析工具获得坡度、坡向、地形起伏度数据。NDVI（MOD13A 3）、LST（MOD11A2）数据从美国航空局（https://ladsweb.modaps.eosdis.nasa.gov/［2021-01-01］）中获取，空间分辨率为1 km，NDVI 时间分辨率为月，LST时间分辨率为8 d。其中，LST数据包括白天地表温度（LSTd）与夜间地表温度（LSTn）以及白天与夜间的温度差（LSTd-n），由每8 天LST 求平均得出月平均温度LST。

3 研究方法

本文构建了一种双阶段卫星降水产品质量提升方法。第一阶段，通过考虑降水与其他因素（如地形、NDVI、地表温度、经纬度等）之间的关系，借助SRF 将10 km GPM 卫星降水产品降尺度到1 km 空间分辨率（D_GPM）。第二阶段，将实测站点降水数据与D_GPM 及上述自变量（如地形、NDVI、地表温度、经纬度等）相结合，再次借助SRF 生成高精度高空间分辨率的降水数据。该方法的流程图如图2所示。

图2 双阶段降尺度流程图Fig.2 Two-stage downscaling flowchart

3.1 随机森林（RF）模型

随机森林（RF）算法是Breiman 于2001 年提出的一种基于bagging 的集成学习方法，通过构造多种决策树来处理自变量与因变量之间的关系，可用作数据的分类和回归预测。RF 通过构建大量的树模型，对多种特征值的重要性进行整合筛选，充分考虑不同特征值之间的重要性，选择最优样本特征值进而找到最优解，并求得所有预测值的平均值作为最终预估值。与传统回归预测方法相比，RF 可以处理复杂多维的特征值，回归预测更准确，采取树状重复放回抽样，避免过拟合，且不必考虑线性回归中的多重共线性问题，稳健性更好，在多种领域得到广泛应用。RF 模型的通用公式如下表示：

式中，P（s0）为s0处的降水预测值，Xi（s0）（i=1，2，3，…，k）为s0处的自变量，k为自变量个数，ε为预测误差。

3.2 顾及空间相关性的随机森林（SRF）模型

为充分考虑邻近降水之间的空间相关性，本文构建了顾及空间相关性的空间随机森林（SRF）模型。其中，SRF 是在RF 的基础上引入克里金估计值作为输入变量，其通用公式如下表示：

式中，Ps（s0）为s0处的克里金降水估计值。

克里金插值法是一种无偏最优估计值方法（史文娇 等，2012；Chen 和Li，2019；Kim 等，2013），对于预估点s0处的降水值PS（s0）可通过搜索该点周围邻近n个站点的降水P（x）经线性加权求得，其计算公式为

式中，PS（s0）为s0处的降水估计值，λi为克里金权重系数，n为周围邻近站点数，P（si）为站点si处的降水值。

克里金权重系数λi不仅取决于周围邻近降水站点与预估点的距离，还取决于邻近点空间分布，其值是通过求解克里金目标函数实现，即PS（s0）对真实值P（s0）的估计为无偏（式（4））且方差最小（式（5））。

由于区域性降水分布具有较大的空间相关性（Sekulić 等，2020），因此，借助克里金插值可充分考虑样本之间的空间相关性。

3.3 SRF双阶段降尺度

本研究的双阶段卫星降水产品质量提升方法大致分为SRF 降尺度与SRF 点面融合两部分，具体步骤如下：

（1）克里金插值。将10 km 分辨率的GPM 数据通过克里金插值获得10 km 分辨率的GPM10kmkri和1 km 分辨率的GPM1kmkri，对站点实测降水数据（RGS）进行克里金插值获得1 km分辨率RGS1kmkri。

（2）将1 km 空间分辨率的所有自变量X1km（包括NDVI、LSTd、LSTn、LSTd-n、DEM、坡度、坡向、地形起伏度以及经纬度）分别使用像元平均重采样至10 km空间分辨率，记为X10km。

（3）以步骤（2）中重采样自变量X10km以及步骤（1）中GPM10kmkri为自变量，原始GPM 为因变量建立SRF降尺度模型：

式中，PGPM10km(s0)为s0处的GPM降水值。

（4）将所有1 km 空间分辨率的自变量X1km输入到SRF 降尺度模型获得1 km 空间分辨率的卫星降水D_GPMsrf：

（5）对降尺度数据进行点面融合，即建立以RGS为因变量的校正模型，自变量包括以步骤（4）中的降尺度结果D_GPMsrf、1 km空间分辨率下的自变量X1km以及RGS1kmkri，即：

式中，RGS(s0)为站点s0处的实测降水。

（6）SRF 模型训练以后，将全部1 km 分辨率自变量输入到SRF 融合模型中，获取融合校正后的1 km降水数据DC_GPMsrf。

3.4 评价方法

为评价新方法（SRF-SRF）的可行性和高效性，本文选取了7 种方法与其进行比较，包括3 种传统方法，即GWR、RF、BPNN；3 种基于SRF 框架下的方法，即GPM 经双线性插值降尺度后使用SRF方法进行点面融合方法（Bi-SRF）、对GPM 使用SRF 降尺度后利用站点实测降水进行地理差异分析校准方法（SRF-GDA）、对GPM 年降水使用SRF 降尺度后按月比例分解并经SRF 校正方法（SRFdis）；最后一种方法为借助克里金直接对站点降水插值（记为Kriging）。

本文采用十折交叉验证方法验证所有方法的性能。选取的精度评价指标包括平均绝对误差（MAE）、均方根误差（RMSE）、相关系数（CC）（Lu等，2020）。3种评价指标公式如下：

式中，Poi、Psi分别代表第i个站点降水实测值和该站点对应的降水预估数据，分别代表站点实测降水的算术平均值和该站点对应降水预估数据算术平均值，n表示气象站点个数。相关系数CC表示两种数据相关性，值越接近1 表示相关程度越好。RMSE 和MAE 用于评价预估降水与实测数据之间的误差，值越小表示精度越高。

4 结果与分析

采用SRF-SRF 方法对四川省GPM 月降水数据进行降尺度与点面融合研究，并将该方法的结果与经典方法的结果在月、季、年3种时间尺度上进行比较。

4.1 月尺度分析

表1 显示了各模型在2015 年—2019 年全部月份上的计算精度。结果显示，在所有模型结果中，GWR 与BPNN 效果最差，可能是这两种方法没有有效处理降水与环境变量之间的复杂关系。基于站点的kriging比RF方法好，且所有考虑空间相关性模型（包括SRF-SRF、Bi-SRF、SRF-GDA、SRFdis）计算结果精度都优于传统的方法，说明了引入降水之间的空间相关性信息对模型精度提升尤为重要。与原始GPM相比，SRF-SRF的MAE和RMSE分别降低了19.51%、16.35%。整体而言，SRF-SRF 精度优于Bi-SRF、SRF-GDA、SRFdis方法。

表1 2015年—2019年全部月份各模型精度对比Table 1 Precision comparison of each model in all months from 2015 to 2019

各个方法在月尺度上的MAE、RMSE、CC 箱线图如图3 所示。结果显示，BPNN 方法在所有模型中效果最差，MAE、RMSE、CC 等3种指标的中位数分别为22.66 mm、30.48 mm、0.64，其次为GWR、RF、Kriging 方法，这与全部月份精度对比（表1）一致。基于SRF的4种方法结果精度都优于传统方法，其中SRF-SRF方法的MAE、RMSE、CC等3 种指标的中位数分别为15.66 mm、21.03 mm、0.81，而其他3 种基于SRF 方法的3 个指标分别为15.83—16.15 mm、21.41—22.27 mm、0.77～0.79，进一步说明了SRF-SRF的高效性。

图3 2015年—2019年月尺度精度指标箱线图Fig.3 Box plot of monthly scale accuracy indicators from 2015 to 2019

图4 为GWR、BPNN、RF、Kriging、SRFdis、SRF-SRF 在月尺度上各站点RMSE 分布图。由于受地势和气候影响（图1），来自印度洋水汽在四川中部盆地与高原相交地段受阻凝结形成降水，故雅安与乐山地区降水偏多，导致其较大的RMSE。GWR、BPNN、RF 方法在各个站点的RMSE 均高于基于SRF 的模型，尤其在四川盆地降水较多地带。这是因为在盆地区域站点分布较为均匀，考虑站点之间的相关性可以有效校正该区域的误差。由图4（e）和 图4（f）表明，与SRFdis 相比，SRF-SRF 方法在四川东北地区与四川盆中地区略好一些。

图4 2015年—2019年月尺度RMSE站点分布图Fig.4 RMSE distribution of all sites on a monthly scale from 2015 to 2019

4.2 季尺度分析

表2 显示了各模型在不同季节上的精度表现。四川地区气候普遍表现为冬暖夏热，夏季降水多，春秋冬降水少，所以在夏季降水的GPM 反演降水误差较高。具体而言，GWR 与BPNN 方法在冬季的表现要优于其他季节，可能该地区冬季气候较为温暖，降水与不同环境因子之间相关性较强。3种传统方法中，RF优于GWR、BPNN，而所有基于SRF 方法（除冬季的SRF-GDA 方法外）的精度都优于这3 种方法。SRF-GDA 在冬季的效果较差，主要是因为冬季降水较少，且在进行点面融合时，未考虑站点以外的环境变量。SRF-SRF 方法在这所有模型中方法效果最好，其中在冬季提升效果明显；相较于原始GPM，SRF-SRF 的MAE 降低了40.69%，RMSE降低了44.24%，CC提高了43.36%。

表2 2015年—2019年季尺度各模型精度对比Table 2 Accuracy comparison of various models on a seasonal scale from 2015 to 2019

4.3 年尺度分析

图5显示了各模型在年尺度上的精度表现。由于四川省从2015 年—2018 年年降水量逐年增加，且在2018 年降水量最多，所以各方法在2018 年的MAE、RMSE 最大。所有方法中，BPNN 最差，其次为GWR，其与月尺度和季尺度表现效果一致。与Kriging 与GWR、BPNN 相比，RF 略好一些；基于SRF 模型效果表现均优于传统的方法。整体而言，SRF-SRF 精度最优，而SRF-GDA 方法在所有基于SRF 模型中表现最差，这主要是因为后者在降水校正时未考虑环境变量有关。

图5 2015年—2019年年尺度精度指标Fig.5 2015—2019 annual scale accuracy indicators

4.4 降水的空间分布特征

图6 显示了各方法在2018 年7 月的降水空间分布图，其中，该月份为5 年内降水最多的月份。从降水分布可以看出，四川省地区降水分布空间差异明显，与地形走势较为接近，即降水多分布于四川中东部地势较低地区，而川西地区由于地势突起降水相对较少。原始GPM（图6（a））降水有效捕捉了川西地区的降水空间分布与降水量，而在四川中部地区却低估了降水，这可能是因为四川中部地形复杂，遥感卫星监测时受四川中部地势和潮湿气候影响，导致反演的卫星降水与实测降水具有一定的误差。经降尺度融合后，各种方法明显改善了中部地区的低估情况，空间分布整体上呈一致性。就空间分辨率而言，原始GPM（图6（a））降水数据分辨率较为粗糙，克里金插值站点结果（图6（b））虽然分辨率提升，但图像比较模糊，细节特征不明显；其他方法降水空间分布细节信息丰富，空间分辨率得到显著提升。GWR 方法（图6（c））的降水空间分布存在较多的噪声，即在部分区域呈现断崖式降低，不符合降水连续性分布的特点；RF 方法（图6（d））的降水在四川东部与四川南部呈现部分块状分布，但整体上相对于GWR 更平滑一些；SRFSRF 方法（图6（e））更能刻画出四川中部地区降水的空间分布情况，且减少了RF 预测模型结果的块状分布，降水分布在空间上更具有连续性，符合降水分布的特点。

图6 2018年7月降水空间分布图Fig.6 The spatial distribution of precipitation in July 2018

5 讨 论

5.1 RF模型变量重要性分析

降水是大气与环境相互作用的结果，因此降水与当地的地形、植被等因素存在着较高的相关性。地理空间位置与地形因素以及植被因素一直以来是研究降水降尺度的常用变量（李净和张晓，2015），然而，在一些高海拔地区，受积雪覆盖等影响，NDVI 和地形因素与降水之间的关系往往不太明显。LST在植被稀少的地区与降水之间存在着显著关系（Jing 等，2016），即降水能影响局部地区的温度变化。因此本方法考虑了空间位置、地形因素、植被因素与地表温度作为降水降尺度的环境变量。基于RF 模型的相对变量重要性（图7）分析表明，Kriging 插值降水的重要性最大，其次为SRF 降尺度值，说明研究区域站点实测降水之间具有显著的相关性。地表温度的重要性要高于NDVI，可能研究区域中高海拔地区较多，温度较低不利于植被生长，难以有效反映降水和NDVI 关系。地形因子中坡向的重要性最低，其次是坡度。而空间位置对降水的影响也很重要，这与研究区域所处的气候有很大关系。

图7 SRF模型中各变量的重要性Fig.7 The importance of variables in the SRF model

5.2 NDVI滞后性影响分析

有研究表明，NDVI 可能对降水存在3 个月以上的滞后性（胡实 等，2020；Karbalaye Ghorbanpour等，2021）。为验证这种长时间滞后性是否对降尺度融合产生影响，本研究将年均NDVI 融入到年降水降尺度中，然后按月降水占的比例获得月降水降尺度结果，最终进行SRF 点面融合（即SRFdis）。图8给出了2015 年—2019年月尺 度SRF-SRF 降水预估值与SRFdis 降水预估值的散点分布图。两种方法的离散程度与实测降水的拟合程度都接近于1∶1 线，与原始GPM（表1）相比均显著提高了GPM 降水的准确性。SRF-SRF 与SRFdis 相比，前者的MAE 和RMSE 值比后者分别降低了0.36 mm 和0.76 mm。因此，与年平均NDVI作为环境变量相比，选择月NDVI 值可以有效提高计算结果精度。

图8 SRF-SRF与SRFdis降水预估值与实测降水散点分布图Fig.8 Comparison between SRF-SRF and SRFdis

6 结 论

为提高GPM IMERG 遥感降水产品空间分辨率与准确性，本文以RF 模型为基础，构造了一种顾及空间相关性的空间随机森林模型，主要包括降尺度和点面融合两个阶段。首先使用SRF 模型并结合降水之间的空间相关性以及环境变量对GPM IMERG 降尺度，然后以实测降水产品为基准再次基于SRF 模型对降尺度后的结果点面融合，获得高空间分辨率和高精度降水产品。实验结果表明：

（1）在复杂地形区域考虑降水空间相关性可有效提升降水产品空间分布细节特征和降水的准确性，并且通过对RF 模型变量重要性分析得出kriging 插值信息重要性最大，说明了研究区域内降水之间的空间相关性不可忽视。该方法解决了研究区域内因降水产品空间分辨率粗糙和准确度低等导致难以开展精细化水文研究等难题，为降水产品精细化研究提供了技术支持。

（2）为验证本文方法的有效性，将该方法实验结果与其他经典方法结果加以对比，结果表明不同时间尺度上该方法表现效果均优于其他方法，并且本文方法的预估值与实测降水更具一致性，说明在降尺度与点面融合两阶段都考虑与降水相关的影响因素可提高结果精度。

后续研究中，将使用其他高时空分辨率降水驱动因素（如土壤湿度、地表温度、风速等）并融合不同遥感降水产品，进一步提升遥感降水产品时空分辨率（如逐日、逐小时）和精度。