CN114463616B - 一种基于Stacking与EMOS-CSG的多源卫星降水融合方法 - Google Patents

Abstract

本发明公开了一种基于Stacking与EMOS‑CSG的多源卫星降水融合方法，步骤包括：数据预处理，包括获取雨量站、原始卫星降水和辅助变量数据；借助地理加权回归模型，对每种原始卫星产品进行降尺度；采用Stacking集成框架第一层多个基学习器分别对每种降尺度卫星产品进行偏差校正；利用Stacking集成框架第二层学习器对各降尺度卫星产品的偏差校正结果进行集成；基于EMOS‑CSG方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品。本发明通过Stacking与EMOS‑CSG的融合方法将不同卫星降水产品的优势有效结合，并借助机器学习模型校正了卫星产品的偏差，提高了融合产品的精度，可为水文模型提供了高精度的降水输入。

Description

一种基于Stacking与EMOS-CSG的多源卫星降水融合方法

技术领域

本发明属于水文与气象技术领域，具体涉及一种基于Stacking与EMOS-CSG的多源卫星降水融合方法。

背景技术

高精度降水数据对流域的洪水预报、干旱监测及水资源管理具有重要意义。获取降水数据的途径通常有两种：一是通过地面雨量站直接获取，二是通过遥感手段(雷达和卫星为主)间接获取，两种途径各有优劣。雨量站数据精度高，资料系列长，但往往站点密度不高，很难精确反映降水在空间上的分布；卫星降水覆盖范围广，但受反演算法、云层性质和传感器的影响，一般存在较大的系统偏差和随机误差；雷达降水时空分辨率高，但是布设困难，维护成本高，且易受地形影响。因此，充分利用各降水产品的优势，将各种来源的降水数据有效“融合”，不仅可以获得更真实的空间分布特征，还可以提高降水估计的准确性。

传统的降水融合方法包括平均偏差校正、最优插值和克里金等，这些方法通常是针对单一卫星产品进行融合。近年来针对多卫星降水融合的方法不断涌现出来，包括地理加权岭回归、动态聚类贝叶斯模型平均和加权平均最小二乘等，这些融合方法可以结合不同来源降水的优势，提供更加可靠的降水分布信息；除上述方法外，机器学习算法因其强大的自学能力和处理非线性问题的优越性而被应用于降水融合中。常用的机器学习模型包含k最近邻方法、多元自回归样条方法、支持向量机、多项对数线性模型和人工神经网络、卷积神经网络长短时记忆网络等。

从多种机器学习算法中筛选出最优算法或者采用一种新的机器学习算法进行融合是可行的，但不同机器学习算法各有优劣，如何充分利用各机器学习算法优势在降水融合中研究较少；此外，传统融合方法通常是基于原始卫星分辨率(0.25°或0.1°)进行融合，所得的融合产品对于水文和气象应用仍然过于粗糙。

发明内容

本发明所要解决的技术问题是，提供一种既可以提高融合产品的空间分辨率，又可以结合不同机器学习算法的优势，有效提高融合产品精度的方法。

为解决上述技术问题，本发明所采用的技术方案为：

一种基于Stacking与EMOS-CSG的多源卫星降水融合方法，包括如下步骤：

步骤1、数据预处理，包括获取雨量站、原始卫星降水和辅助变量数据，辅助变量包括经纬度、高程、坡度坡向、归一化植被指数(NDVI)、土壤湿度、风向风速、温度和大气压；

步骤2、借助地理加权回归模型，对每种原始卫星产品进行降尺度；

步骤3、采用Stacking集成框架第一层6种基学习器(6种机器学习模型)分别对每种降尺度卫星产品进行偏差校正，得到每种产品的6种校正结果；其中，6种机器学习模型包括随机森林(RF)、K最近邻算法(KNN)、极端梯度提升树(XGB)、轻量梯度提升树(LGB)、CatBoost和梯度提升树(GBM)；

步骤4、利用Stacking集成框架第二层学习器(随机森林)对各降尺度卫星产品的偏差校正结果进行集成；

步骤5、基于左移截尾伽马分布的集成模型输出统计(EMOS-CSG)方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品；其中EMOS-CSG最终求解方程为：

λ＝θk[1-G_k,θ(δ)][1-G_k+1,θ(δ)]-δ[1-G_k,θ(δ)]² (6)

式(6)中，λ为左移截尾伽马分布的均值，即融合产品均值；θ为伽马分布的尺度参数；k为形状参数；G_k,θ(δ)是以k,θ为参数的原始伽马分布的累计分布函数；δ>0表示向左偏移量；G_k+1,θ(δ)是以k+1,θ为参数的原始伽马分布的累计分布函数。

进一步的，所述步骤2中借助地理加权回归模型，对每种原始卫星产品进行降尺度，包含如下步骤：

步骤21、选取与降雨关系密切的经纬度、高程、坡度坡向和NDVI作为自变量，以卫星产品作为因变量，构建地理加权回归模型：

式(1)中，P_i ^L代表第i个原始卫星产品的降雨； 分别为与第i个原始卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/>代表对应的系数；/>是常数项的系数；/>是残差项；

步骤22、采用双线性插值对式(1)中系数项和自变量因子项重采样到细分辨率；

步骤23、计算细分辨率的卫星降水：

式(2)中P_i ^H为第i个细分辨率卫星产品的降水； 分别为与第i个细分辨率卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/> 代表相应重采样到细分辨率后的系数；/>是重采样到细分辨率后常数项的系数；

步骤24、采用普通克里金法处理式(1)中的残差项，获取和细分辨率卫星产品空间分辨率一致的残差

步骤25、将步骤23得到的细分辨率的卫星降水P_i ^H与步骤24残差相加得到最终降尺度细分辨率卫星产品的降水值；

式(3)中，P_i ^HH为第i个最终降尺度细分辨率卫星产品的降水值。

进一步的，所述步骤3中采用Stacking集成框架第一层6种基学习器(6种机器学习模型)分别对每种降尺度卫星产品进行偏差校正，得到每种产品的6种校正结果，包含如下步骤：

步骤31、选择地面雨量站点所在细网格处的经纬度、土壤湿度、风速风向、温度、大气压、卫星降水因子和雨量站降水作为原始数据集，不同降尺度卫星产品的原始数据集除卫星降水不同外，其他因子均相同；

步骤32、利用双线性插值将各因子的分辨率统一到同一分辨率；

步骤33、在第一层6种基学习器中利用每种降尺度卫星产品的原始数据集进行训练，得到第一层的6种降水校正值；

步骤34、重复步骤33，依次得到每种降尺度卫星产品的6种偏差校正值。

进一步的，所述步骤4中利用Stacking集成框架第二层学习器(随机森林)对每种降尺度卫星产品的偏差校正结果进行集成，包含如下步骤：

步骤41、将每种降尺度卫星产品的6种偏差校正结果作为新的训练集，其中每种降尺度卫星产品的训练集均不相同；

步骤42、采用第二层学习器(随机森林)对步骤41得到的训练集进行训练，得到第二层的集成结果；

步骤43、重复步骤42，依次得到各降尺度卫星产品的Stacking集成结果。

进一步的，所述步骤5中基于左移截尾伽马分布的集成模型输出统计(EMOS-CSG)方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品，包含如下步骤：

步骤51、根据各卫星Stacking集成结果与伽马分布均值μ和方差σ²的关系求解gamma分布的均值和方差：

式(4)中，a₀…a_N为均值参数；b₀，b₁是方差参数；代表各卫星Stacking集成结果均值；f₁…f_N代表各卫星Stacking集成结果；其中待定参数a₀…a_N，b₀，b₁和式(6)中δ通过最小化连续等级概率评分来确定；

步骤52、根据伽马分布的均值μ、方差σ²和形状参数k和尺度参数θ的关系求解形状和尺度参数：

步骤53、通过式(6)对各卫星Stacking集成结果进行融合计算，得到高时空分辨率的降水融合产品。

本发明所达到的有益效果：本发明提供的基于Stacking与EMOS-CSG的多源卫星降水融合方法，首先借助地理加权回归方法提高了原始卫星产品的空间分辨率，然后根据Stacking集成框架将不同机器学习算法的优势结合起来，有效校正了卫星降水的偏差，最后根据集成模型输出统计方法将不同卫星产品有效融合，得到了高时空分辨率的降水数据集，为水文模型提供了高精度的降水输入。

附图说明

图1是本发明基于Stacking与EMOS-CSG的多源卫星降水融合方法的流程图。

具体实施方式

下面结合附图和具体实例对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

如图1所示，

步骤1、对目标流域进行数据预处理，包括获取流域内雨量站日降水数据，并从相应官网上下载同一时间内TMPA-3B42RT、CMORPH、GSMaP_NRT和PERSIANN四种原始卫星产品的降水数据，空间分辨率分别为0.25°、0.25°、0.1°和0.25°，并通过累计叠加整理成日尺度降雨，从地理空间数据云网站上下载DEM数据，并使用ArcGIS软件进一步提取经纬度、高程、坡度和坡向等地形因子，归一化植被指数(NDVI)来自美国国家航空航天局(NASA)的MOD13A3数据集，土壤湿度来自NASA的SMAP数据集，温度、大气压、风速和风向数据从欧洲中期天气预报中心(ECMWC)的再分析数据集ERA5-Land获取；

步骤2、借助地理加权回归模型，对每种原始卫星产品进行降尺度，得到1km的降尺度产品，包括以下步骤：

步骤21，选取与降雨关系密切的经纬度、高程、坡度坡向和NDVI作为自变量，以卫星产品作为因变量，构建地理加权回归模型：

式(1)中，P_i ^L代表第i个原始卫星产品的降雨； 分别为与第i个原始卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/>代表对应的系数；/>是常数项的系数；/>是残差项；

步骤22、采用双线性插值对式(1)中系数项和自变量因子项重采样到细分辨率1km；

步骤23、计算细分辨率的卫星降水：

式(2)中P_i ^H为第i个细分辨率卫星产品的降水； 分别为与第i个细分辨率卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/> 代表相应重采样到细分辨率后的系数；/>是重采样到细分辨率后常数项的系数；

步骤24、采用普通克里金法处理式(1)中的残差项，获取和细分辨率卫星产品空间分辨率一致的残差

步骤25、将步骤23得到的细分辨率的卫星降水P_i ^H与步骤24残差相加得到最终降尺度细分辨率卫星产品的降水值；

式(3)中，P_i ^HH为第i个最终降尺度细分辨率卫星产品的降水值。

经过降尺度后，可以获得TMPA-3B42RT₁(1km)、CMORPH₁(1km)、GSMaP_NRT₁(1km)、PERSIANN₁(1km)4种降尺度卫星产品数据集。

步骤3、采用Stacking集成框架第一层6种基学习器(6种机器学习模型)分别对每种降尺度卫星产品进行偏差校正，得到每种产品的6种校正结果，包括以下步骤：

步骤31、选择地面雨量站点所在细网格处的经纬度、土壤湿度、风速风向、温度、大气压、卫星降水因子和雨量站降水作为原始数据集，不同降尺度卫星产品的原始数据集除卫星降水不同外，其他因子均相同；

步骤32、利用双线性插值将各因子的分辨率统一到同一分辨率1km；

步骤33、在第一层6种基学习器中利用每种降尺度卫星产品的原始数据集进行训练，得到第一层的6种降水校正值；其中，6种机器学习模型包括随机森林(RF)、K最近邻算法(KNN)、极端梯度提升树(XGB)、轻量梯度提升树(LGB)、CatBoost和梯度提升树(GBM)；

步骤34、重复步骤33，依次得到每种降尺度卫星产品的6种偏差校正值。

步骤4、利用Stacking集成框架第二层学习器(随机森林)对各降尺度卫星产品的偏差校正结果进行集成，包括以下步骤：

步骤41、将每种降尺度卫星产品的6种偏差校正结果作为新的训练集，其中每种降尺度卫星产品的训练集均不相同；

步骤42、采用第二层学习器(随机森林)对步骤41得到的训练集进行训练，得到第二层的集成结果；

步骤43、重复步骤42，依次得到各降尺度卫星产品的Stacking集成结果，分别为TMPA-3B42RT₂(1km)、CMORPH₂(1km)、GSMaP_NRT₂(1km)、PERSIANN₂(1km)。

步骤5、基于左移截尾伽马分布的集成模型输出统计(EMOS-CSG)方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品，包含如下步骤：

步骤51、根据各卫星Stacking集成结果与伽马分布均值μ和方差σ²的关系求解gamma分布的均值和方差：

式(4)中，a₀…a_N为均值参数；b₀，b₁是方差参数；代表各卫星Stacking集成结果均值；f₁…f_N代表各卫星Stacking集成结果；其中待定参数a₀…a_N，b₀，b₁通过最小化连续等级概率评分来确定；

步骤52、根据伽马分布的均值μ、方差σ²和形状参数k和尺度参数θ的关系求解形状和尺度参数：

步骤53、左移截尾伽马分布的三个参数求解之后，伽马分布的累计分布函数就可以得到，可对各卫星Stacking集成结果进行融合计算，得到高时空分辨率的降水融合产品均值λ：

λ＝θk[1-G_k,θ(δ)][1-G_k+1,θ(δ)]-δ[1-G_k,θ(δ)]² (6)

式(6)中，λ为左移截尾伽马分布的均值，即降水融合产品均值；G_k,θ(δ)是以k,θ为参数的原始伽马分布的累计分布函数；δ>0表示向左偏移量，通过最小化连续等级概率评分来确定；G_k+1,θ(δ)是以k+1,θ为参数的原始伽马分布的累计分布函数。

经过步骤53之后，可以获得目标流域1km分辨率日尺度融合数据集。本实施案例采用六个指标来评价融合产品的性能，分别是均方根误差、相关系数、平均绝对误差、误报率、频率偏差指数、击中率。各指标的计算如下。

均方根误差RMSE：

相关系数CC：

平均绝对误差MAE：

误报率FAR：

频率偏差指数FBI：

击中率FH：

上述六个评价指标中：P_i表示卫星降水量；表示卫星降水量均值；O_i表示雨量站降水量；/>表示雨量站降水量均值；n表示站点总数；N₁₁表示雨量站和卫星均捕捉到降雨的次数；N₁₀表示卫星捕捉到而雨量站未捕捉到降雨的次数；N₀₁表示雨量站捕捉到而卫星未捕捉到降雨的次数。

表1融合产品的评价结果

融合产品的评价结果见表1，从表1可以看出，4种原始卫星产品均存在较大的均方根误差和平均绝对误差，相关系数较低，说明卫星降水与真实降水之间相关性较低；降尺度卫星产品的精度略有提高，虽然空间分辨率提高到1km，仍然存在较大的误差，因此有必要进一步融合地面雨量站数据对卫星产品进行偏差校正；经过Stacking和EMOS-CSG融合后，最终融合产品的精度有了显著的提高，具体表现为：均方根误差由均值12.45mm降低到4.33mm(降低了65.2％)，相关系数由均值0.20提高到0.64(提高了0.44)，平均绝对误差由均值10.83mm降低到3.13mm(降低了71.1％)，误报率由均值0.18降低到0.07(降低了0.11)，频率偏差指数由均值0.87提高到0.96(提高了0.09)，击中率由均值0.80提高到0.93(提高了0.13)。综上所述，融合产品在目标流域精度高，可认为本发明提出的基于Stacking与EMOS-CSG的多源卫星降水融合方法是合理有效的。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (9)

1.一种基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于，包括如下步骤：

步骤1、数据预处理：包括获取雨量站、原始卫星降水和辅助变量数据，所述辅助变量包括经纬度、高程、坡度坡向、NDVI、土壤湿度、风向风速、温度和大气压；

步骤2、降尺度：借助地理加权回归模型，对每种原始卫星产品进行降尺度；

步骤3、Stacking第一层基学习器校正：采用Stacking集成框架第一层6种基学习器分别对每种降尺度卫星产品进行偏差校正，得到每种产品的6种校正结果；其中，6种基学习器为随机森林、K最近邻算法、极端梯度提升树、轻量梯度提升树、CatBoost和梯度提升树；

步骤4、利用Stacking集成框架第二层学习器对各降尺度卫星产品的偏差校正结果进行集成；所述第二层学习器为随机森林；

步骤5、基于EMOS-CSG方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品；其中EMOS-CSG最终求解方程为：

λ＝θk[1-G_k,θ(δ)][1-G_k+1,θ(δ)]-δ[1-G_k,θ(δ)]² (6)

式(6)中，λ为左移截尾伽马分布的均值，即融合产品均值；θ为伽马分布的尺度参数；k为形状参数；G_k,θ(δ)是以k,θ为参数的原始伽马分布的累计分布函数；δ>0表示向左偏移量；G_k+1,θ(δ)是以k+1,θ为参数的原始伽马分布的累计分布函数；

所述步骤2中借助地理加权回归模型，对每种原始卫星产品进行降尺度，包含如下步骤：

步骤21、选取与降雨关系密切的经纬度、高程、坡度坡向和NDVI作为自变量，以卫星产品作为因变量，构建地理加权回归模型：

式(1)中，P_i ^L代表第i个原始卫星产品的降雨； 分别为与第i个原始卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/>代表对应的系数；/>是常数项的系数；/>是残差项；

步骤22、采用双线性插值对式(1)中系数项和自变量因子项重采样到细分辨率；

步骤23、计算细分辨率的卫星降水：

式(2)中P_i ^H为第i个细分辨率卫星产品的降水； 分别为与第i个细分辨率卫星产品空间分辨率一致的归一化植被指数、高程、坡度、坡向、经度、纬度值；/> 代表相应重采样到细分辨率后的系数；/>是重采样到细分辨率后常数项的系数；

步骤24、采用普通克里金法处理式(1)中的残差项，获取和细分辨率卫星产品空间分辨率一致的残差

步骤25、将步骤23得到的细分辨率的卫星降水P_i ^H与步骤24残差相加得到最终降尺度细分辨率卫星产品的降水值；

式(3)中，P_i ^HH为第i个最终降尺度细分辨率卫星产品的降水值。

2.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：步骤1原始卫星降水数据自TMPA-3B42RT、CMORPH、GSMaP_NRT和PERSIANN中任意一个或多个中获取。

3.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：步骤1中从地理空间数据云网站上下载DEM数据，并使用ArcGIS软件进一步提取经纬度、高程、坡度和坡向数据。

4.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：步骤1中NDVI、土壤湿度数据来自美国国家航空航天局。

5.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：步骤1中温度、大气压、风向风速数据从欧洲中期天气预报中心获取。

6.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：所述步骤22中采用双线性插值对式(1)中系数项和自变量因子项重采样到细分辨率1km。

7.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：所述步骤3中采用Stacking集成框架第一层6种基学习器分别对每种降尺度卫星产品进行偏差校正，得到每种产品的6种校正结果，包含如下步骤：

步骤31、选择地面雨量站点所在细网格处的经纬度、土壤湿度、风速风向、温度、大气压、卫星降水因子和雨量站降水作为原始数据集，不同降尺度卫星产品的原始数据集除卫星降水不同外，其他因子均相同；

步骤32、利用双线性插值将各因子的分辨率统一到同一分辨率；

步骤33、在第一层6种基学习器中利用每种降尺度卫星产品的原始数据集进行训练，得到第一层的6种降水校正值；

步骤34、重复步骤33，依次得到每种降尺度卫星产品的6种偏差校正值。

8.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：所述步骤4中利用Stacking集成框架第二层学习器对每种降尺度卫星产品的偏差校正结果进行集成，包含如下步骤：

步骤41、将每种降尺度卫星产品的6种偏差校正结果作为新的训练集，其中每种降尺度卫星产品的训练集均不相同；

步骤42、采用第二层学习器对步骤41得到的训练集进行训练，得到第二层的集成结果；

步骤43、重复步骤42，依次得到各降尺度卫星产品的Stacking集成结果。

9.根据权利要求1所述的基于Stacking与EMOS-CSG的多源卫星降水融合方法，其特征在于：所述步骤5中基于EMOS-CSG方法对Stacking集成结果进行融合，得到高时空分辨率的降水融合产品，包含如下步骤：

步骤51、根据各卫星Stacking集成结果与伽马分布均值μ和方差σ²的关系求解gamma分布的均值和方差：

μ＝a₀+a₁f₁…+a_Nf_N

式(4)中，a₀…a_N为均值参数；b₀，b₁是方差参数；代表各卫星Stacking集成结果均值；f₁…f_N代表各卫星Stacking集成结果；其中待定参数a₀…a_N，b₀，b₁和式(6)中δ通过最小化连续等级概率评分来确定；

步骤52、根据伽马分布的均值μ、方差σ²和形状参数k和尺度参数θ的关系求解形状和尺度参数：

步骤53、通过式(6)对各卫星Stacking集成结果进行融合计算，得到高时空分辨率的降水融合产品。