CN115357847A - 一种基于误差分解的日尺度星地降水融合方法 - Google Patents

Abstract

本发明公开了一种基于误差分解的日尺度星地降水融合方法，借助误差分解模型将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差，针对各项误差选取地理加权回归模型和相似矩阵法计算空间误差估计值和降水背景值，最后通过“加法模型”框架、“乘法模型”框架和回归残差值估计进行降水融合，得到各网格点的降水融合值。本发明更细致地针对站点处的各误差组分进行空间降水误差估计，减少降水融合模型的不确定性；地理加权回归模型反映了降水误差与其影响因素之间的非平稳关系，相似矩阵法修正了降水背景场中的错误信息，充分利用各误差组分的有效信息，能够获得更好的降水空间估计数据，误差小于常规的降水插值方法。

Description

一种基于误差分解的日尺度星地降水融合方法

技术领域

本发明属于水文与气象技术领域，具体涉及一种基于误差分解的日尺度星地降水融合方法。

背景技术

目前用于星地降水数据融合的方法如平均偏差修正方法、概率密度匹配法、双核平滑方法、最优插值法和地理加权回归法等。如李运龙等在文献基于地理加权回归克里金的降水数据融合及其在水文预报中的应用[J].长江流域资源与环境,2017,26(09):1359-1368中结合地理加权回归模型与克里金插值方法估算降水背景误差，并应用加法模型框架生成最终的降水融合值。将卫星降水数据的误差进行分解可以更好地描述其误差的各个组成部分和来源，因此误差分解方法常用于卫星降水精度验证与评估的研究中。如李新焕等在文献卫星降水产品在川渝地区的精度验证和误差分解[J].人民长江,2022,53(01):97-103.中采用多种评价指标和误差分解方法对GSMaP和IMERG进行对比分析，以验证卫星降水产品在川渝地区的精度。

现有的融合方法从降水的统计意义出发构建降水融合模型，通常采用回归模型构建地面降水与环境变量和卫星降水间的关系，实现卫星降水与地面降水的融合，或是计算站点处地面降水和卫星降水间误差并基于此估算整个空间的降水误差以实现降水融合，而少有方法针对降水各项误差进行剖析而进行降水融合。

发明内容

本发明提供一种基于误差分解的日尺度星地降水融合方法，通过借助误差分解模型，将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差，并针对各项误差选取地理加权回归模型和相似矩阵法计算空间误差估计值和降水背景值，通过“加法模型”框架、“乘法模型”框架和回归残差值估计进行降水融合，能有效提高卫星遥感反演降水的精度，为水文模型的输入提供数据支撑。

为了实现上述目的，本发明所采用的技术方案为：

一种基于误差分解的日尺度星地降水融合方法，包括如下步骤：

步骤1、获取雨量站观测降水、卫星遥感反演降水以及与降水空间估计相关的地理因子、气候因子和关联因子数据，并对降水数据和各因子数据进行预处理；

步骤2、借助误差分解模型，将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差；

步骤3、针对命中误差和漏报误差，利用地理加权回归模型得到这两种误差在各网格点的空间误差估计值，针对误报误差，利用相似矩阵法修正降水背景值；

步骤4、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，利用IDW法和Kriging法估计回归残差值，并将回归残差值还原到融合初始结果，基于雨量站观测降水选择精度更高的模型框架，得到各网格点的降水融合值；

进一步的，所述步骤1中对降水数据和各因子数据进行预处理包括：

步骤11、采用箱线图法对卫星遥感反演降水数据中的异常值进行处理，将异常值修正为当日卫星降水平均值；

步骤12、卫星降水数据的Drizzles处理，将0.1mm/d作为区分降水事件与无降水事件的阈值，将0.1mm/d以下的卫星遥感反演降水量设为0；

步骤13、将与降水空间估计相关的地理因子、气候因子和关联因子数据重采样到与卫星遥感反演降水数据相同的空间分辨率(0.1°)，并对各因子数据进行归一化处理。

进一步的，所述步骤2中将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差，具体包括如下步骤：

步骤21、计算雨量站观测降水与卫星遥感反演降水之间的误差；

步骤22、将计算的误差分为三种独立的误差组分：命中误差、漏报误差和误报误差，三种误差的定义如下：

式中，e代表雨量站观测降水与卫星遥感反演降水之间的误差，作为降水背景误差；e_Hit代表命中误差；e_Miss代表漏报误差；e_False代表误报误差；P_G代表雨量站观测降水；P_S代表卫星遥感反演降水。

进一步的，所述步骤3中针对命中误差和漏报误差，利用地理加权回归模型得到这两种误差在各网格点的空间估计值，针对误报误差，利用相似矩阵法修正降水背景场，具体包括如下步骤：

步骤31、针对命中误差和漏报误差，以选取的地理因子、气候因子和关联因子数据为自变量，以雨量站观测降水与卫星遥感反演降水之间的误差作为因变量，采用逐步回归法选择对降水贡献量大的因子，选取Bi-square函数作为空间权函数，构建地理加权回归估计空间误差值和误差比例因子：

式中，w_ij代表空间权函数矩阵，由站点所在网格i与无测站网格j之间的距离d_ij来确定；d_ij代表样本点与估计点间的欧式空间距离；b代表带宽，表征权重与距离函数关系的非负衰减参数；

步骤32、针对误报误差，根据站点所在网格和无测站网格之间的相似性划分相似网格来修正降水背景场，将卫星遥感反演降水数据设为降水背景场，与出现误报误差的网格相似的无测站网格处的卫星降水量设置为0，站点所在网格和无测站网格之间的相似性如下式所示：

式中，SM代表站点所在网格和无测站网格之间的相似性，由站点所在网格和无测站网格的高程、经度和纬度三个变量来确定；

代表站点所在网格的第i个变量；

代表无测站网格的第i个变量；ΔP_i代表站点所在网格中变量的最小值和最大值之间的差值。

进一步的，所述步骤4中基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，利用IDW法和Kriging法估计回归残差值，并将回归残差值还原到融合初始结果，基于雨量站观测降水选择精度更高的模型框架，得到各网格点的降水融合值，具体包括如下步骤：

步骤41、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，“加法模型”框架和“乘法模型”框架的公式分别为：

P_a＝P_bc+e_b

P_m＝P_bc*r_b

式中，P_a代表基于“加法模型”框架得到的降水融合初始结果；P_m代表基于“乘法模型”框架得到的降水融合初始结果；P_bc代表经相似矩阵法修正过的降水背景场；e_b代表根据命中误差和漏报误差利用地理加权回归估计的空间误差值；r_b代表根据命中误差和漏报误差利用地理加权回归估计的误差比例因子；

步骤42、分别利用IDW法和Kriging法估计回归残差值，将回归残差值还原到融合初始结果；

步骤43、基于雨量站观测降水，计算精度指标，在“加法模型”框架和“乘法模型”框架中选择精度更高的模型框架，在IDW法和Kriging法中选择精度更高的，重复步骤4得到各网格点的降水融合值。

进一步的，步骤4之后还可以包括：

步骤5、分别从统计精度评估与水文模拟精度评估两个角度对比评估融合降水数据的精度，具体包括如下步骤：

步骤51、采用留一交叉验证法，分别计算统计指标进行精度评估，并绘制混淆矩阵图对融合前后的误差组分进行可视化分析；

步骤52、以雨量站观测降水数据、原始卫星降水数据和基于误差分解的融合降水数据为降水输入驱动水文模型，分别采用静态参数法和动态参数法进行参数率定，对比评估融合降水数据的水文模拟精度。

与现有技术相比，本发明具有如下显著优点：

本发明更细致地针对站点处的各误差组分进行空间降水误差估计，减少降水融合模型的不确定性；地理加权回归模型反映了降水误差与其影响因素之间的非平稳关系，相似矩阵法修正了降水背景场中的错误信息，充分利用各误差组分的有效信息；易于实施，能够获得更好的降水空间估计数据，误差小于常规的降水插值方法，可为水文模型等提供精度更高的降水输入源。

附图说明

图1为本发明方法流程示意图；

图2为GSMaP-MVK降水产品在2016年6月1日降水命中误差和漏报误差基于地理加权回归的误差估计值空间分布图；

图3为GSMaP-MVK降水产品在2016年6月1日降水误报误差基于相似网格法的降水背景值空间分布图；

图4为GSMaP-MVK降水产品融合前后的混淆矩阵图；

图5为GSMaP-MVK融合前后分别驱动VIC模型的日径流过程图，率定期为2014-2017年，验证期为2017-2018年。

具体实施方式

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

本发明所述的一种基于误差分解的日尺度星地降水融合方法，基于雨量站观测降水与卫星遥感反演降水数据进行融合；以黄河源区作为研究区，分别融合GSMaP-MVK、GSMaP-Gauge和IMERG-F卫星降水产品，三种降水产品的空间尺度均为0.1°，时间尺度为日尺度。流程图如图1所示，具体实现包括如下步骤：

步骤1、获取雨量站观测降水、卫星遥感反演降水以及与降水空间估计相关的地理因子、气候因子和关联因子数据，并对降水数据和各因子数据进行预处理，包括以下步骤：

步骤11、从中国气象数据网获取所选流域内雨量站观测降水数据，从各卫星降水产品官网上下载同一时段的GSMaP-MVK、GSMaP-Gauge和IMERG-F三种卫星产品的降水数据，空间分辨率分别为0.1°，并通过累加逐小时降水整理成日尺度降水，采用箱线图法对各卫星遥感反演降水数据中的异常值进行处理，将异常值修正为当日该卫星降水产品的降水平均值；

步骤12、卫星降水数据的Drizzles处理，将0.1mm/d作为区分降水事件与无降水事件的阈值，将0.1mm/d以下的卫星遥感反演降水量设为0；

步骤13、下载研究区DEM数据，提取经纬度、高程和坡度坡向等地理因子信息，从欧洲中期天气预报中心的再分析数据集ERA5-Land获取温度、大气压和风速风向等气候因子信息，从美国国家航空航天局(NASA)的SMAP数据集获取土壤湿度数据，从NASA的MOD13A3数据集获取归一化植被指数(NDVI)数据，将搜集的各数据采用双线性内插法重采样到与卫星遥感反演降水数据相同的空间分辨率，即0.1°×0.1°，并对各因子数据进行归一化处理。

步骤2、借助误差分解模型，将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差，包括以下步骤：

步骤21、计算雨量站观测降水与卫星遥感反演降水之间的误差；

步骤22、将计算的误差分为三种独立的误差组分：命中误差、漏报误差和误报误差，三种误差的定义如下：

式中，e代表雨量站观测降水与卫星遥感反演降水之间的误差，作为降水背景误差；e_Hit代表命中误差；e_Miss代表漏报误差；e_False代表误报误差；P_G代表雨量站观测降水；P_S代表卫星遥感反演降水。

步骤3、针对命中误差和漏报误差，利用地理加权回归模型得到这两种误差在各网格点的空间估计值，针对误报误差，利用相似矩阵法修正降水背景值，包括以下步骤：

步骤31、针对命中误差和漏报误差，以GSMaP-MVK降水产品在2016年6月1日降水为例，图2为命中误差和漏报误差基于地理加权回归的误差估计值空间分布图，以选取的地理因子、气候因子和关联因子数据为自变量，以雨量站观测降水与卫星遥感反演降水之间的误差作为因变量，采用逐步回归法选择对降水贡献量大的因子，选取Bi-square函数作为空间权函数构建地理加权回归估计空间误差值和误差比例因子：

式中，w_ij代表空间权函数矩阵，由站点所在网格i与无测站网格j之间的距离d_ij来确定；d_ij代表样本点与估计点间的欧式空间距离；b代表带宽，表征权重与距离函数关系的非负衰减参数；

步骤32、针对误报误差，以GSMaP-MVK降水产品在2016年6月1日降水为例，图3为误报误差基于相似网格法的降水背景值空间分布图，根据站点所在网格和无测站网格之间的相似性划分相似网格来修正降水背景场，将卫星遥感反演降水数据设为降水背景场，与出现误报误差的网格相似的无测站网格处的卫星降水量设置为0，站点所在网格和无测站网格之间的相似性如下式所示：

式中，SM代表站点所在网格和无测站网格之间的相似性，由站点所在网格和无测站网格的高程、经度和纬度三个变量来确定；

代表站点所在网格的第i个变量；

代表无测站网格的第i个变量；ΔP_i代表站点所在网格中变量的最小值和最大值之间的差值。

步骤4、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，利用IDW法和Kriging法估计回归残差值，并将回归残差值还原到融合初始结果，基于雨量站观测降水，选择“加法模型”框架和“乘法模型”框架中精度更高的模型框架，得到各网格点的降水融合值，包含如下步骤：

步骤41、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，“加法模型”框架和“乘法模型”框架的公式分别为：

P_a＝P_bc+e_b

P_m＝P_bc*r_b

式中，P_a代表基于“加法模型”框架得到的降水融合初始结果；P_m代表基于“乘法模型”框架得到的降水融合初始结果；P_bc代表经相似矩阵法修正过的降水背景场；e_b代表根据命中误差和漏报误差利用地理加权回归估计的空间误差值；r_b代表根据命中误差和漏报误差利用地理加权回归估计的误差比例因子；

步骤42、分别利用IDW法和Kriging法估计回归残差值，将回归残差值还原到融合初始结果；

步骤43、基于雨量站观测降水，计算皮尔逊相关系数(CC)、平均绝对误差(MAE)和均方根误差(RMSE)这三种精度指标，比较两种框架得到的降水融合结果的精度，选择“加法模型”框架作为空间误差估计值和降水背景值耦合的框架，比较两种回归残差值估计方法，选择IDW法作为回归残差值的插值方法，并重复步骤4得到各网格点的降水融合值，三种卫星降水产品分别记作EC-GSMaP-MVK、EC-GSMaP-Gauge和EC-IMERG-F。

步骤5、分别从统计精度评估与水文模拟精度评估两个角度对比评估融合降水数据的精度。

步骤51、根据步骤4获得的融合降水数据绘制混淆矩阵图对融合前后的误差组分进行可视化分析，以GSMaP-MVK为例，分别绘制融合前后的混淆矩阵图(图4)。混淆矩阵本质上是描述测试数据与真值匹配性能的数值表，混淆矩阵可以用来评估卫星降水产品融合前后对不同量级降水事件的捕获情况。与融合前相比，融合后的EC-GSMaP-MVK在大于20mm的降水事件探测能力上略微提升，说明融合后的降水产品对强降水事件有相对较好的探测能力。

此外，采用留一交叉验证法以五个统计指标来评价融合产品的性能，分别是皮尔逊相关系数(CC)、平均误差(ME)、平均绝对误差(MAE)、均方根误差(RMSE)、相对偏差(BIAS)，各指标的计算公式如下：

式中：n表示精度评估中所用到的总样本数；G_i表示雨量站观测降水值；

表示雨量站观测降水的平均值；S_i表示卫星降水估计值；

表示卫星降水估计量的平均值。

基于误差分解的融合降水产品的评价结果见表1，从表1可以看出，相比于原始卫星降水产品，融合降水产品精度指标均有不同程度的提升。比较融合降水产品和两种常用的降水空间插值估计方法反距离加权法(IDW)和普通克里金法(OK)可以发现，EC-GSMaP-Gauge和EC-IMERG-F均较降水空间插值估计方法精度更高，这表明经过与地面观测降水数据的融合，EC-GSMaP-Gauge和EC-IMERG-F在日尺度径流模拟方面具有替代地面站点插值降水数据的潜力。EC-GSMaP-Gauge的相关系数CC值由原始卫星降水数据的0.69升到了0.72，平均绝对误差MAE与均方根误差值RMSE均下降，表明GSMaP-Gauge原始卫星降水数据经过融合校正处理后，融合降水数据误差更低和数据精度更高。EC-GSMaP-MVK的相关系数提高了0.12，ME、MAE、RMSE和BIAS的值均显著地下降，这表明经过与地面观测降水数据的融合，EC-GSMaP-MVK降水数据对日降水的过度高估得到了明显改善。EC-IMERG-F较IMERG-F具有更高的CC和较低的MAE和RMSE。

表1黄河源区2014～2018年不同降水数据日尺度精度评价指标

步骤52、以雨量站观测降水数据、原始卫星降水数据和基于误差分解的融合降水数据为降水输入驱动VIC水文模型，分别采用静态参数法和动态参数法进行参数率定，对比评估融合降水数据的水文模拟精度。以GSMaP-MVK采用静态参数法进行参数率定为例，图5为GSMaP-MVK融合前后分别驱动VIC模型的日径流过程图，率定期为2014-2017年，验证期为2017-2018年。从图中可以看出，融合前的GSMaP-MVK的模拟结果较差，经过误差分解融合模型得到的EC-GSMaP-MVK能显著提高径流模拟的精度，虽然模拟结果然较实测值仍然偏多，但径流的变化趋势能与实测径流保持一致，这表明误差分解融合模型能够有效修正原始GSMaP-MVK降水，提高其在水文模拟应用中的效果。综上所述，本实施例中误差分解的融合降水产品相比于原始卫星降水产品，精度指标均有不同程度的提升，可认为本发明提出的一种基于误差分解的日尺度星地降水融合方法是合理有效的。

以上所述仅是本发明的优选实施例而已，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (6)

1.一种基于误差分解的日尺度星地降水融合方法，其特征在于，包括如下步骤：

步骤1、获取雨量站观测降水、卫星遥感反演降水以及与降水空间估计相关的地理因子、气候因子和关联因子数据，并对降水数据和各因子数据进行预处理；

步骤2、借助误差分解模型，将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差；

步骤3、针对命中误差和漏报误差，利用地理加权回归模型得到这两种误差在各网格点的空间估计值，针对误报误差，利用相似矩阵法修正降水背景场；

步骤4、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，利用IDW法和Kriging法估计回归残差值，并将回归残差值还原到融合初始结果，基于雨量站观测降水选择精度更高的模型框架和回归残差值估计方法，得到各网格点的降水融合值。

2.根据权利要求1所述的一种基于误差分解的日尺度星地降水融合方法，其特征在于：所述步骤1中对降水数据和各因子数据进行预处理，具体包括如下步骤：

步骤11、采用箱线图法对卫星遥感反演降水数据中的异常值进行处理，将异常值修正为当日该卫星降水平均值；

步骤12、卫星降水数据的Drizzles处理，将0.1mm/d作为区分降水事件与无降水事件的阈值，将0.1mm/d以下的卫星遥感反演降水量设为0；

步骤13、将与降水空间估计相关的地理因子、气候因子和关联因子数据重采样到与卫星遥感反演降水数据相同的空间分辨率(0.1°)，并对各因子数据进行归一化处理。

3.根据权利要求1所述的一种基于误差分解的日尺度星地降水融合方法，其特征在于：所述步骤2中将雨量站观测降水与卫星遥感反演降水之间的误差分解为命中误差、漏报误差和误报误差，具体包括如下步骤：

步骤21、计算雨量站观测降水与卫星遥感反演降水之间的误差；

步骤22、将计算的误差分为三种独立的误差组分：命中误差、漏报误差和误报误差，三种误差的定义如下：

式中，e代表雨量站观测降水与卫星遥感反演降水之间的误差，作为降水背景误差；e_Hit代表命中误差；e_Miss代表漏报误差；e_False代表误报误差；P_G代表雨量站观测降水；P_S代表卫星遥感反演降水。

4.根据权利要求3所述的一种基于误差分解的日尺度星地降水融合方法，其特征在于：所述步骤3具体包括如下步骤：

步骤31、针对命中误差和漏报误差，以选取的地理因子、气候因子和关联因子数据为自变量，以雨量站观测降水与卫星遥感反演降水之间的误差作为因变量，采用逐步回归法选择对降水贡献量大的因子，选取Bi-square函数作为空间权函数构建地理加权回归估计空间误差值和误差比例因子：

式中，w_ij代表空间权函数矩阵，由站点所在网格i与无测站网格j之间的距离d_ij来确定；d_ij代表样本点与估计点间的欧式空间距离；b代表带宽，表征权重与距离函数关系的非负衰减参数；

步骤32、针对误报误差，根据站点所在网格和无测站网格之间的相似性划分相似网格来修正降水背景场，将卫星遥感反演降水数据设为降水背景场，与出现误报误差的网格相似的无测站网格处的卫星降水量设置为0，站点所在网格和无测站网格之间的相似性如下式所示：

式中，SM代表站点所在网格和无测站网格之间的相似性，由站点所在网格和无测站网格的高程、经度和纬度三个变量来确定；

代表站点所在网格的第i个变量；

代表无测站网格的第i个变量；ΔP_i代表站点所在网格中变量的最小值和最大值之间的差值。

5.根据权利要求1所述的一种基于误差分解的日尺度星地降水融合方法，其特征在于：所述步骤4具体包括如下步骤：

步骤41、基于“加法模型”框架和“乘法模型”框架分别将空间误差估计值和降水背景值耦合得到降水融合初始结果，“加法模型”框架和“乘法模型”框架的公式分别为：

P_a＝P_bc+e_b

P_m＝P_bc*r_b

式中，P_a代表基于“加法模型”框架得到的降水融合初始结果；P_m代表基于“乘法模型”框架得到的降水融合初始结果；P_bc代表经相似矩阵法修正过的降水背景场；e_b代表根据命中误差和漏报误差利用地理加权回归估计的空间误差值；r_b代表根据命中误差和漏报误差利用地理加权回归估计的误差比例因子；

步骤42、分别利用IDW法和Kriging法估计回归残差值，将回归残差值还原到融合初始结果；

步骤43、基于雨量站观测降水，计算精度指标，在“加法模型”框架和“乘法模型”框架中选择精度更高的模型框架，在IDW法和Kriging法中选择精度更高的回归残差值估计方法，重复步骤4得到各网格点的降水融合值。

6.根据权利要求1所述的一种基于误差分解的日尺度星地降水融合方法，其特征在于：步骤4之后还包括：

步骤5、分别从统计精度评估与水文模拟精度评估两个角度对比评估融合降水数据的精度，具体包括如下步骤：

步骤51、采用留一交叉验证法，分别计算统计指标进行精度评估，并绘制混淆矩阵图对融合前后的误差组分进行可视化分析；

步骤52、以雨量站观测降水数据、原始卫星降水数据和基于误差分解的融合降水数据为降水输入驱动水文模型，分别采用静态参数法和动态参数法进行参数率定，对比评估融合降水数据的水文模拟精度。

Images