CN115755221A - 一种基于数理不确定性分析的多源降水数据融合方法 - Google Patents

Abstract

本发明公开了一种基于数理不确定性分析的多源降水数据融合方法，包括：S1.降水数据的时空分辨率统一；S2，降水数据的系统偏差校正；S3，构建多源降水数据的误差矩阵；S4，多源降水自适应融合；S5，融合降水数据集的晴雨分类误差校正；S6，降水融合数据后处理。本发明最大限度地克服了多源降水数据在不确定性分析与融合过程中对地面观测的依赖，完善了地面观测站点稀疏区的降水融合技术，提升了大空间范围的多源降水数据融合稳定性，显著降低了降水事件的误报率，避免了卫星数据噪音的过拟合问题，最终得到更为准确、可靠的大范围降水数据产品，为进一步提升区域水资源管理与自然灾害预警和监测能力提供有力支撑。

Description

一种基于数理不确定性分析的多源降水数据融合方法

技术领域

本发明涉及遥感水文气象数据处理与遥感信息不确定分析领域，尤其涉及一种基于数理不确定性分析的多源降水数据融合方法。

背景技术

降水是水循环的重要组成部分，降水资料的可靠程度在一定程度上决定了径流与蒸散发等陆面水文过程的模拟精度，同时也是影响洪水预报、干旱监测与预警及农作物产量估算精度的关键因素。卫星遥感和再分析系统是目前获取大范围、高时空分辨率降水信息的主要途径。然而，单一的卫星或再分析降水数据通常伴随着较大的不确定性。因此，综合利用地面观测、多种卫星与再分析降水资料的多源数据融合是降低降水数据不确定性、提升降水估算精度的重要途径，也是当今遥感水文的研究热点。

随机误差与晴雨分类误差是卫星与再分析降水数据的主要误差来源。现有的多源降水数据融合研究中，大多通过地面观测资料对多源降水数据进行校正。高密度地面监测网可有效地降低上述两种误差，提升多源数据融合结果的精度，但此方法对于广泛分布的地面站点稀疏区存在较大的局限性。数理不确定性分析方法可以通过采用多组彼此独立的数据集，经数学建模求解卫星降水数据误差，从而解决地面观测稀疏区的数据不确定性评估问题。然而目前已有的大部分数据集彼此间通常存在一定的误差相关性，不完全满足彼此独立的假设前提，因此给传统数理不确定性分析方法的应用带来了一系列困扰。此外，目前的融合方法会累积子数据集中所有的误报降水事件，造成较高的误报率。

发明内容

本发明的目的是提供一种基于数理不确定性分析的多源降水数据融合方法，以解决现有遥感多源降水数据误差的全面量化与融合技术中对地面观测的依赖，导致无法客观估算随机误差、晴雨分类误差及各数据集间的互相关性误差，使多源降水数据融合结果偏离最优值的问题。

为此，本发明采用以下技术方案：

一种基于数理不确定性分析的多源降水数据融合方法，包括以下步骤：

S1.降水数据的时空分辨率统一：

对常用的至少三种降水数据集进行时空分辨率统一的预处理，通过算数平均与最近距离插值法将多源降水数据的时空分辨率分别统一至≥0.5小时和≥0.05°，得到多个时空分辨率统一的降水数据集；

S2，降水数据的系统偏差校正：

基于全球月尺度降水量的参考值，对步骤S1得到的多个时空分辨率统一的降水数据集进行逐月校正，得到校正后的多源降水数据集；

S3，构建多源降水数据的误差矩阵：

基于三重排列分析方法TC和降水数据集之间的误差独立性，建立基于四重排列分析方法QC的多源降水不确定性分析方案；使用QC方法计算步骤S2校正后的多源降水数据集的降水误差及各数据集间的误差互相关性，得到多源降水数据误差矩阵；

S4，多源降水自适应融合：

使用最小二乘法与步骤S3构建的多源降水数据误差矩阵，求解各个卫星降水数据集的最优融合权重，当权重出现负值时，自动剔除不确定性最大的降水数据，重新计算融合权重，直至所有选取的卫星数据融合权重为正值，得到自适应的优化融合降水数据集；

S5，融合降水数据集的晴雨分类误差校正：

通过分类变量融合算法，逐格点计算最优晴雨时间序列，使用该最优晴雨时间序列来纠正步骤S4得到的优化融合降水数据集的晴雨分类误差，得到晴雨分类误差校正后的降水融合数据集；

S6，降水融合数据后处理：

再次使用全球月尺度降水量的参考值，对步骤S5所生成的降水融合数据集进行系统偏差校正，确保消除融合过程中可能引入的系统偏差，完成多源降水数据融合。

优选的是，步骤S1中，将所述时空分辨率分别统一为日尺度和0.1°。对于空间分辨率高于0.1°的降水数据，将对0.1°内所有网格内的降水强度进行平均；对于空间分辨率较低的降水数据集，采用降水融合中最近距离法将数据插值到0.1°空间分辨率；对于高时间分辨率的降水数据，通过时段内累加的方法对时间分辨率进行统一。

步骤S2中对多源卫星降水数据进行逐月校正的公式为：

式中，

和

分别为步骤S2中的多源降水数据集在M月的校正和原始日降水量；S^M为M月的月尺度参考降水量，对

低于5mm的情况不进行校正。

步骤S3构建多源降水数据的误差矩阵的步骤如下：

1)基于TC的降水数据误差独立性判别：

通过TC方法，对多源数据中的任意三组降水数据集x₁、x₂、x₃的精度进行计算：

其中，COV代表协方差；

与

分别代表x₁、x₂、x₃三组降水数据所对应的精度；当至少两个数据集x_s1与x_s2的精度

及

不随第三组数据集x_i的变化而变化时，判定数据集x_s1与x_s2独立于其它降水数据；

2)基于QC的降水数据误差信息估算：

根据上述TC分析所选定的独立数据集，构建QC不确定性分析方案；使用QC方法，对S2中多源降水数据的误差结构进行数学建模与联合求解，计算得到多源降水数据的误差与各数据集间的互相关性误差；

3)构建多源降水数据误差矩阵：

使用QC方法，通过变更S2中不同的降水数据集，循环求解所有候选降水数据集的误差及互相关性误差，构建误差矩阵：

其中C为误差矩阵，当E的下标中两个数字相同时，其表示某数据集的误差方差(如E₁₁)，当其下标中两个数字不同时，表示两个数据集间的误差协方差(如E₁₂)。

步骤S4进行多源降水自适应融合的具体方法为：

根据最小二乘法，多源降水数据融合表达式为：

式中，x_m为融合后的降水数据集；λ_i为融合权重；其中，各个数据集的最优融合权重表示为：

其中

为式(7)计算的降水数据集误差矩阵的逆矩阵中的不同元素；

当降水数据的随机误差或互相关性误差较高时，使用迭代程序剔除基于QC的误差方差最高的单个降水产品，重新计算各个降水数据集的融合权重，直到所有保留的融合权重λi值全部为正值。

步骤S5的具体方法为：

根据CTC-M算法的原理，首先将降水数据转换为晴雨二进制时间序列，设三组独立的晴雨分类时间序列分别为d₁，d₂与d₃，其协方差表示为：

Q₁₂＝COV(d₁，d₂)＝f(P)(2π₁-1)(2π₂-1) (10)

Q₁₃＝COV(d₁，d₃)＝f(P)(2π₁-1)(2π₃-1) (11)

Q₂₃＝COV(d₂，d₃)＝f(P)(2π₂-1)(2π₃-1) (12)

其中，Q为不同晴雨二进制时间序列的协方差；P为晴雨时间序列真值；f(P)是P统计属性的函数；π为不同降水数据集的晴雨分类精度；

根据式(10)到(12)，计算得到与各个晴雨分类精度正相关的统计量v：

据此优化融合这三组卫星数据，得到精度更高的晴雨分类时间序列：

d_m＝sign(w₁d₁+w₂d₂+w₃d₃) (16)

其中，d₁、d₂、d₃为三种产品的晴雨时间序列；w₁、w₂、w₃为对应的融合权重；w_i为x_i产品对应的融合权重；

使用优化融合的晴雨分类时间序列(d_m)，校正步骤S4所生成的优化融合的降水数据集(即x_m)。

本发明基于最小二乘的优化加权平均理论，通过四重排列分析估算降水数据的随机误差方差及数据集间互相关性误差协方差，建立可以自动遴选并融合优势降水数据集的自适应融合方法，降低了降水估算值的随机误差；通过分类变量不确定性分析与融合，结合晴雨分类误差校正方法，进一步提升融合结果的有、无降水辨识精度。

与现有技术相比，本发明具有以下有益效果：

1.本发明最大限度地克服了多源降水数据在不确定性分析与融合过程中对地面观测的依赖，完善了地面观测站点稀疏区的降水融合技术；

2.本发明开发了符合多源降水数据误差客观特征的优化融合方法，通过自适应融合方法提升了大空间范围的多源降水数据融合稳定性；

3.本发明单独考虑了晴雨分类误差，显著降低了降水事件的误报率；

4.本发明的方法考虑了不同降水数据间的互相关误差，避免了卫星数据噪音的过拟合问题，最终得到更为准确、可靠的大范围降水数据产品，为进一步提升区域水资源管理与自然灾害预警和监测能力提供有力支撑。

附图说明

图1是本发明的多源降水数据融合方法的流程图；

图2是本发明中多源降水数据误差信息分析方法的流程图。

具体实施方式

以下结合附图和实施例对本发明的基于数理不确定性分析的多源降水数据融合方法进行详细说明。

实施例1

鉴于目前主流降水数据的分辨率及陆面水文过程模拟的需求，本发明选取7组卫星与再分析降水数据(ERA5、SM2Rain、IMERG、TRMM 3B42RT、CMORPH、PERSIANN-CDR、CHIRPS)作为评价实例，以数理不确定性分析方法为基础，对降水数据中的多重不确定性信息进行评估，构建“自适应”多源融合方法来降低降水估算值中的随机误差，并根据CTC-M算法的原理，构建全球晴雨时间序列数据集，对降水数据进行晴雨误差校正，最终形成完整的多源降水数融合体系。

参见图1，本发明的一种基于数理不确定性分析的多源降水数据融合方法包括以下步骤：

S1，降水数据的时空分辨率统一：

统一多源降水数据的时空辨率，本实施例中，分别统一为日尺度和0.1°。对于空间分辨率高于0.1°的降水数据，将对0.1°内所有网格内的降水强度进行平均；对于空间分辨率较低的降水数据集，采用降水融合中广泛应用的“最近距离法”，将数据插值到0.1°空间分辨率。对于高时间分辨率(如0.5h、1h等)的降水数据，通过时段内累加的方法对时间分辨率进行统一，得到多个时空分辨率统一的降水数据集。

S2，降水数据的系统偏差校正：

使用全球月尺度降水量的参考值(如Climate Hazards group Precipitationclimatology，CHPclim)，对多源卫星降水数据进行逐月校正，以纠正降水产品的系统偏差：

式中，

和

分别为步骤S2中的多源降水数据集在M月的校正和原始日降水量，S^M为M月的月尺度参考降水量(如CHPclim)，对

低于5mm的情况不进行校正。

S3，构建多源降水数据的误差矩阵：

参见图2，综合使用三重排列分析(Triple collocation analysis，TC)与四重排列分析(Quadruple collocation analysis,QC)方法来构建步骤S2校正后的多源降水数据集的误差矩阵，进而用于多源降水数据融合，具体步骤包括：

1)基于TC的降水数据误差独立性判别：

通过TC方法，对多源数据中的任意三组降水数据集(例如x₁,x₂,x₃)的精度进行计算：

其中COV代表协方差，

与

分别代表x₁,x₂,x₃三组降水数据集所对应的精度。使用TC计算所有排列组合情况下的所有降水数据精度。当某些数据集(如x_s1与x_s2)的精度不随第三组数据集的变化而变化时，即TC(x_s1,x_s2,x_i)所估算的

及

不随x_i的变化而变化时，可判定x_s1与x_s2独立于其它降水数据。

2)基于QC的降水数据误差信息估算：

在QC方法中，考虑了四个降水产品(例如，x_s1,x_s2,x_i及x_j)，其中x_s1和x_s2独立于x_i与x_j，由前述TC分析中所选定。QC允许x_i和x_j为非独立的数据集，即存在互相关误差的可能性。QC方法可以同时估算上述四种数据集的误差A及互相关性误差b：

根据A与b，可以求解：

a＝(A^TA)^-1A^Tb. (6)

其中，a为使用QC计算所得的列向量，其包括了各降水数据的误差及误差协方差信息；A^T为A的转置矩阵。

3)构建多源降水数据误差矩阵：

将式(5)与(6)定义为QC算子，即QC(x_s1,x_s2,x_i,x_j)，通过变更x_i与x_j所对应的数据集，循环求解所有候选降水数据集的误差及互相关性误差，构建误差矩阵：

其中C为误差矩阵，当E的下标中两个数字相同时，其表示某数据集的误差方差(如E₁₁)，当其下标中两个数字不同时，表示两个数据集间的误差协方差(如E₁₂)。

S4，多源降水数据自适应融合：

根据最小二乘法，多源降水数据融合表达式为：

式中，x_m为融合后的降水数据集；λ_i为融合权重。其中，各个数据集的最优融合权重可表示为：

其中

为式(7)计算的降水数据集误差矩阵的逆矩阵中的不同元素。当降水数据的随机误差或互相关性误差较高时，λ_i可能超出[0，1]合理区间。

为此，使用迭代程序剔除基于QC的误差方差最高的单个降水产品，重新计算各个降水数据集的融合权重，直到所有保留的降水融合权重(λi值)全部为正值，进而形成“自适应”的多源降水融合方法。

S5，融合降水数据集的晴雨分类误差校正：

通过对每个像素应用分类变量融合算法(Categorical Triple Collocationbased Merging，CTC-M)来计算最优晴、雨时间序列。根据CTC-M算法的原理，首先利用三组相对独立的卫星数据构建晴雨二进制时间序列(d₁，d₂，d₃)，再通过联合求解的方式估算各个降水数据集的晴雨分类精度。

CTC-M首先将降水数据转换为晴、雨二进制时间序列，以三组独立的晴雨分类时间序列，如d₁，d₂与d₃为例，其协方差可表示为：

Q₁₂＝COV(d₁，d₂)＝f(P)(2π₁-1)(2π₂-1) (10)

Q₁₃＝COV(d₁，d₃)＝f(P)(2π₁-1)(2π₃-1) (11)

Q₂₃＝COV(d₂，d₃)＝f(P)(2π₂-1)(2π₃-1) (12)

其中，Q为不同晴雨二进制时间序列的协方差，P为晴雨时间序列真值，f(P)是P统计属性的函数，π为不同降水数据集的晴雨分类精度。根据式(10)到(12)，可以计算得到与各个晴雨分类精度正相关的统计量(v)：

据此可以优化融合这三组卫星数据，从而得到精度更高的晴雨时间序列：

d_m＝sign(w₁d₁+w₂d₂+w₃d₃) (16)

其中，d₁、d₂、d₃为三种产品的晴雨时间序列；w₁、w₂、w₃为对应的融合权重；w_i为x_i产品对应的融合权重。使用优化融合的晴雨分类时间序列(d_m)，校正S4所生成的优化融合的降水数据集(即x_m)。

S6，降水融合数据后处理：

再次使用公式(1)，利用CHPclim数据对步骤S5中产生的降水时间序列逐月进行系统偏差的校正，以确保消除融合过程中所可能引入的系统偏差，最终完成多源降水数据融合。

对于本发明的方法最终融合的数据集，可使用均方根误差、皮尔逊相关系数、误报率、探测率、HSS技巧指数等多个指标进行评估。

Claims (7)

1.一种基于数理不确定性分析的多源降水数据融合方法，包括以下步骤：

S1.降水数据的时空分辨率统一：

对常用的至少三种降水数据集进行时空分辨率统一的预处理，通过算数平均与最近距离插值法将多源降水数据集的时空分辨率分别统一至≥0.5小时和≥0.05°，得到多个时空分辨率统一的降水数据集；

S2，降水数据集的系统偏差校正：

基于全球月尺度降水量的参考值，对步骤S1得到的多个时空分辨率统一的降水数据集进行逐月校正，得到校正后的多源降水数据集；

S3，构建多源降水数据的误差矩阵：

基于三重排列分析方法TC和降水数据集之间的误差独立性，建立基于四重排列分析方法QC的多源降水不确定性分析方案；使用QC方法计算步骤S2校正后的多源降水数据集的降水误差及各数据集间的误差互相关性，得到多源降水数据的误差矩阵；

S4，多源降水自适应融合：

使用最小二乘法与步骤S3构建的多源降水数据的误差矩阵，求解各个卫星降水数据集的最优融合权重，当权重出现负值时，自动剔除不确定性最大的降水数据，重新计算融合权重，直至所有选取的卫星数据融合权重为正值，得到自适应的优化融合降水数据集；

S5，融合降水数据集的晴雨分类误差校正：

通过分类变量融合算法，逐格点计算最优晴雨分类时间序列，使用该最优晴雨分类时间序列来纠正步骤S4得到的优化融合降水数据集的晴雨分类误差，得到晴雨分类误差校正后的降水融合数据集；

S6，降水融合数据后处理：

再次使用全球月尺度降水量的参考值，对步骤S5所生成的降水融合数据集进行系统偏差校正，确保消除融合过程中可能引入的系统偏差，完成多源降水数据融合。

2.根据权利要求1所述的多源降水数据融合方法，其特征在于：步骤S1中，将所述时空分辨率分别统一为日尺度和0.1°。

3.根据权利要求2所述的多源降水数据融合方法，其特征在于：对于空间分辨率高于0.1°的降水数据，将对0.1°内所有网格内的降水强度进行平均；对于空间分辨率较低的降水数据集，采用降水融合中最近距离法将数据插值到0.1°空间分辨率；对于高时间分辨率的降水数据，通过时段内累加的方法对时间分辨率进行统一。

4.根据权利要求1所述的多源降水数据融合方法，其特征在于，步骤S2中对多源卫星降水数据进行逐月校正的公式为：

式中，

和

分别为步骤S2中的多源降水数据集在M月的校正和原始日降水量；S^M为M月的月尺度参考降水量，对

低于5mm的情况不进行校正。

5.根据权利要求1所述的多源降水数据融合方法，其特征在于，步骤S3构建多源降水数据的误差矩阵的步骤如下：

1)基于TC的降水数据误差独立性判别：

通过TC方法，对多源数据中的任意三组降水数据集x₁、x₂、x₃的精度进行计算：

其中，COV代表协方差；

与

分别代表x₁、x₂、x₃三组降水数据所对应的精度；当至少两个数据集x_s1与x_s1的精度

及

不随第三组数据集x_i的变化而变化时，判定数据集x_s1与x_s2独立于其它降水数据；

2)基于QC的降水数据误差信息估算：

根据上述TC分析所选定的独立数据集(x_s1与x_s2)，构建QC不确定性分析方案；使用QC方法，对S2中多源降水数据的误差结构进行数学建模与联合求解，计算得到多源降水数据的误差与各数据集间的互相关性误差；

3)构建多源降水数据误差矩阵：

使用QC方法，通过变更S2中不同的降水数据集，循环求解所有候选降水数据集的误差及互相关性误差，构建误差矩阵：

其中C为误差矩阵，当E的下标中两个数字相同时，其表示某数据集的误差方差(如E₁₁)，当其下标中两个数字不同时，表示两个数据集间的误差协方差(如E₁₂)。

6.根据权利要求5所述的多源降水数据融合方法，其特征在于，步骤S4进行多源降水自适应融合的具体方法为：

根据最小二乘法，多源降水数据融合表达式为：

式中，x_m为融合后的降水数据集；λ_i为融合权重；其中，各个数据集的最优融合权重表示为：

其中

为式(7)计算的降水数据集误差矩阵的逆矩阵中的不同元素；

当降水数据的随机误差或互相关性误差较高时，使用迭代程序剔除基于QC的误差方差最高的降水数据集，重新计算剩余降水数据集的融合权重，直到所有保留的融合权重λ_i值全部为正值。

7.根据权利要求6所述的多源降水数据融合方法，其特征在于，步骤S5的具体方法为：

根据CTC-M算法的原理，首先将降水数据转换为晴雨二进制时间序列，设三组独立的晴雨分类时间序列分别为d₁，d₂与d₃，其协方差表示为：

Q₁₂＝COV(d₁,d₂)＝f(P)(2π₁-1)(2π₂-1) (10)

Q₁₃＝COV(d₁,d₃)＝f(P)(2π₁-1)(2π₃-1) (11)

Q₂₃＝COV(d₂,d₃)＝f(P)(2π₂-1)(2π₃-1) (12)

其中，Q为不同晴雨二进制时间序列的协方差；P为晴雨时间序列真值；f(P)是P统计属性的函数；π为不同降水数据集的晴雨分类精度；

根据式(10)到(12)，计算得到与各个晴雨分类精度正相关的统计量v：

据此优化融合这三组卫星数据，得到精度更高的晴雨分类时间序列：

d_m＝sign(w₁d₁+w₂d₂+w₃d₃) (16)

其中，d₁、d₂、d₃为三种产品的晴雨时间序列；w₁、w₂、w₃为对应的融合权重；w_i为x_i产品对应的融合权重；

使用优化融合的晴雨分类时间序列d_m，校正步骤S4所生成的最终的优化融合降水数据集x_m。

Images