CN112861072B - 一种星地多源降水自适应动态融合方法 - Google Patents

一种星地多源降水自适应动态融合方法 Download PDF

Info

Publication number
CN112861072B
CN112861072B CN202110175249.6A CN202110175249A CN112861072B CN 112861072 B CN112861072 B CN 112861072B CN 202110175249 A CN202110175249 A CN 202110175249A CN 112861072 B CN112861072 B CN 112861072B
Authority
CN
China
Prior art keywords
precipitation
time
satellite
fine resolution
resolution grid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110175249.6A
Other languages
English (en)
Other versions
CN112861072A (zh
Inventor
张珂
晁丽君
刘林鑫
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202110175249.6A priority Critical patent/CN112861072B/zh
Publication of CN112861072A publication Critical patent/CN112861072A/zh
Application granted granted Critical
Publication of CN112861072B publication Critical patent/CN112861072B/zh
Priority to NL2030246A priority patent/NL2030246A/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/10—Devices for predicting weather conditions
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G01W1/14—Rainfall or precipitation gauges
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Environmental & Geological Engineering (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Atmospheric Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Ecology (AREA)
  • Environmental Sciences (AREA)
  • Computing Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Hydrology & Water Resources (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • Probability & Statistics with Applications (AREA)
  • Databases & Information Systems (AREA)
  • Radio Relay Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种星地多源降水自适应动态融合方法，包括：获取流域粗分辨率和细分辨率的时不变因子和时变因子；时不变因子包括：DEM、坡度、Hand指数及粗糙度；时变因子包括：NDVI、风速、地表温度、散射指数及微波极化指数；对每颗卫星进行降尺度，得到每颗卫星的细分辨率降水；构建自适应动态贝叶斯框架，通过降水误差的分布转换、降水概率密度函数参数优化来计算每颗卫星在实测站点所在网格的最优权重；计算每颗卫星细分辨率网格的最优权重，每颗卫星细分辨率网格降水加权求和得到融合降水。本发明通过自适应动态融合方法融合站点降水和多颗卫星降水，得到细分辨率的融合降水产品，为精细化水文模型的输入提供数据支撑。

Description

一种星地多源降水自适应动态融合方法
技术领域
本发明属于水文和气象技术领域，具体涉及基于获取时间和空间上的动态权重进而实现星地多源降水融合的方法。
背景技术
中国是一个洪涝灾害频发的国家，水文模型是解决洪涝灾害的重要手段。降水作为水文模型最关键的输入源，降水的精度和时效性直接影响模拟结果的精度和可靠性。
目前降水数据的获取方式主要有地面观测、卫星和雷达定量降水估计、模式定量降水预报。长期以来，降水的常规观测主要依赖于布设于地表的观测站点，采用有限的观测结果代表周边几十甚至几百平方公里范围内的真实降水。实际降水的大小、类型等具有显著的时空变异性，地面站点存在以点代面的问题，特别是站点稀少的区域观测降水不能有效反映空间降水的空间变异性，降水观测的空间局限性成为水文研究中的难点。雷达定量降水估计具有空间分辨率高、实时性强的优点，但是因容易受覆盖物的影响，其覆盖范围有限。伴随着国内外卫星遥感技术的发展，基于天气雷达与卫星的遥感降水观测得以不断完善，弥补了地面站点空间分布的不足，也为降水的监测提供了新的手段。目前，卫星遥感在获取时空变化的全球降水方面具有独特的优势，提供了前所未有的卫星降水产品如GPM、TRMM、COMRPH、PERSIANN、PERSIANN_CDR、FY-3B、FY-3C等。卫星定量降水估计具有覆盖面广、观测时间较连续的优势，但是由于遥感探测仪器、反演算法等限制，卫星降水产品的精度相对较低，并且对固态降水的反演能力有限。
为了更好的描述降水的空间变化，充分发挥不同降水来源的优点，提高降水产品的质量，构建高精度、细时空分辨率的水文模型驱动场，采用混合加权逐步回归方法对原始卫星降水进行降尺度、自适应动态贝叶斯融合方法融合多颗卫星降尺度降水与站点降水。
发明内容
为了获取细时空分辨率的降水产品，本发明提供了一种星地多源降水自适应动态融合方法，利用细时空分辨率的时不变因子(地形数据)和时变因子(植被数据、气象数据、辐射指数)，融合站点实测降水和多颗卫星降水产品，得到细分辨率的融合降水产品，为精细化水文模型的输入提供数据支撑。
为了获取细分辨率的降水融合产品，本发明具体采用以下技术方案：
一种星地多源降水自适应动态融合方法，其特征在于，包括以下步骤：
步骤1，获取流域粗分辨率(Coarse Resolution，CR)和细分辨率(FineResolution，FR)的时不变因子，主要包括：DEM(DEMCR、DEMFR)、坡度(SlopeCR、SlopeFR)、Hand指数(HandCR、HandFR)、粗糙度(SRCR、SRFR)；粗分辨率和细分辨率的时变因子分别为：NDVI(NDVICR、NDVIFR)、风速(WSCR、WSFR)、地表温度(LSTCR、LSTFR)、散射指数(SICR、SIFR)、微波极化指数(PDCR、PDFR)；卫星编号用m表示，m＝1,2，……，5分别代表GPM、TRMM、CMORPH、PERSIANN、PERSIANN_CDR；粗分辨率网格编号用i表示，细分辨率网格编号用j表示，站点个数用n表示；
步骤2，基于混合地理加权逐步回归模型对每颗卫星进行降尺度，得到细分辨率降水；
步骤3，构建自适应动态贝叶斯框架，通过降水误差的分布转换、降水概率密度函数参数优化来计算每颗卫星在降水实测站点所在细分辨率网格的最优权重；
步骤4，基于克里金插值通过实测站点所在细分辨率网格求得的权重，计算每颗卫星每个细分辨网格的最优权重，每颗卫星细分辨率降水加权求和得到融合降水。
前述的一种星地多源降水自适应动态融合方法，其特征在于，所述步骤1中获取粗分辨率、细分辨率的时不变和时变因子。
前述的一种星地多源降水自适应动态融合方法，其特征在于，所述步骤2中基于混合地理加权逐步回归模型对多颗卫星进行降尺度，得到每颗卫星细分辨率降水，具体包括以下步骤：
步骤2.1，基于粗分辨率的时不变因子和时变因子，构建混合地理加权逐步回归模型(以某个时刻某颗卫星单个粗分辨率网格为例)：
Figure BDA0002940491190000021
式中，i为粗分辨率网格编号，m为卫星编号，t表示时间，CR表示粗分辨率，
Figure BDA0002940491190000022
为t时刻m卫星第i个粗分辨率网格的降水；
Figure BDA0002940491190000023
为第i个粗分辨率网格的高程，
Figure BDA0002940491190000024
为
Figure BDA0002940491190000025
对应的时不变系数；
Figure BDA0002940491190000026
为第i个粗分辨率网格的坡度，
Figure BDA0002940491190000027
为
Figure BDA0002940491190000028
对应的时不变系数；
Figure BDA0002940491190000029
为第i个粗分辨率网格的Hand指数，
Figure BDA00029404911900000210
为
Figure BDA00029404911900000211
对应的时不变系数；
Figure BDA00029404911900000212
为第i个粗分辨率网格的地表粗糙度，
Figure BDA00029404911900000213
为
Figure BDA00029404911900000214
对应的时不变系数；
Figure BDA00029404911900000215
为t时刻第i个粗分辨率网格的植被归一化指数，
Figure BDA00029404911900000216
为
Figure BDA00029404911900000217
对应的时变系数；
Figure BDA00029404911900000218
为t时刻第i个粗分辨率网格的风速，
Figure BDA00029404911900000219
为
Figure BDA00029404911900000220
对应的时变系数；
Figure BDA00029404911900000221
为t时刻第i个粗分辨率网格的地表温度，
Figure BDA00029404911900000222
为
Figure BDA00029404911900000223
对应的时变系数；
Figure BDA00029404911900000224
为t时刻第i个粗分辨率网格的散射指数，
Figure BDA00029404911900000225
为
Figure BDA00029404911900000226
对应的时变系数；
Figure BDA00029404911900000227
为t时刻第i个粗分辨率网格的微波极化指数，
Figure BDA00029404911900000228
为
Figure BDA00029404911900000229
对应的时变因子；
Figure BDA00029404911900000230
为t时刻第i个粗分辨率网格的残差。
步骤2.2，方程(1)进行逐步回归，逐步引入对降水贡献量大的因子：
逐步引入9个粗分辨率因子，包括：4个时不变因子
Figure BDA00029404911900000231
和5个时变因子
Figure BDA00029404911900000232
进行逐步回归分析和对应的回归系数进行F检验，若检验通过给定显著性水平α，则将该因子保留，逐步回归完成以后留下的因子记为
Figure BDA0002940491190000031
Figure BDA0002940491190000032
为4个时不变因子和5个时变因子的子集，由两个部分组成分别为：
Figure BDA0002940491190000033
(时不变因子子集)、
Figure BDA0002940491190000034
(时变因子子集)。
步骤2.3，假设上步所得
Figure BDA0002940491190000035
由n1个时不变因子和n2个时变因子组成，某个时刻单个细分辨率网格的混合地理加权逐步回归模型：
Figure BDA0002940491190000036
式中，j为细分辨率网格编号，m为卫星编号，t表示时间，FR表示细分辨率，
Figure BDA0002940491190000037
为t时刻第m颗卫星第j个细分辨率网格的降水；
Figure BDA0002940491190000038
为第j个细分辨率网格的第1个时不变因子，
Figure BDA0002940491190000039
为其对应的系数；
Figure BDA00029404911900000310
为第j个细分辨率网格的第n1个时不变因子，
Figure BDA00029404911900000311
为其对应的系数；
Figure BDA00029404911900000312
为第j个细分辨率网格的第1个时变因子，
Figure BDA00029404911900000313
为其对应的系数；
Figure BDA00029404911900000314
为第j个细分辨率网格的第n2个时变因子，
Figure BDA00029404911900000315
为其对应的系数；
Figure BDA00029404911900000316
为t时刻第j个细分辨率网格的残差。
对于t时刻细分辨率单个网格系数求解如下：
Figure BDA00029404911900000317
式中，
Figure BDA00029404911900000318
为时不变因子系数矩阵；
Figure BDA00029404911900000319
为时变因子系数矩阵；
Figure BDA00029404911900000320
为时不变因子矩阵，
Figure BDA00029404911900000321
为时变因子矩阵，I为单位矩阵，W为权重矩阵，ym,CR为第m颗原始卫星降水矩阵，
Figure BDA00029404911900000322
计算t时刻第m颗卫星第j个细分辨率网格的降水
Figure BDA00029404911900000323
(即降尺度后的降水)：
Figure BDA00029404911900000324
式中，sum定义为矩阵的逐项相加函数，
Figure BDA00029404911900000325
作为系统残差可取常数，
Figure BDA00029404911900000326
分别为：
Figure BDA00029404911900000327
通过j次循环得到t时刻第m颗卫星降尺度之后的降水，通过m次循环得到每颗卫星降尺度降水即细分辨率时空一致的降水。
前述的一种星地多源降水自适应动态融合方法，其特征在于，所述步骤3通过降水分布转换、降水概率密度函数参数优化来计算单颗卫星的最优权重，具体包括以下步骤：
步骤3.1，基于不同卫星细分辨率降水、站点降水条件下融合降水的后验分布，根据全概率公式，建立自适应动态贝叶斯框架，融合降水的概率密度可表示为：
Figure BDA0002940491190000041
式中，t1为训练时间段(t1的取值范围为一个动态范围，可为同一个季节内的0-90天，获取t1+1时刻的权重以t1时间段训练数据为基础)；obs为t1时段内实测站点降水矩阵，
Figure BDA0002940491190000042
其中
Figure BDA0002940491190000043
代表t1时段第一个时刻的n个站点降水列向量，
Figure BDA0002940491190000044
代表t1时段第二个时刻的n个站点降水列向量，
Figure BDA0002940491190000045
代表t1时段第t1时刻的n个站点降水列向量，具体表示如下：
Figure BDA0002940491190000046
式中，
Figure BDA0002940491190000047
为t1时段第t1时刻第n个站点的降水值。
Pm,FR为t1时段内站点对应的第m颗卫星细分辨率降水，表示如下：
Figure BDA0002940491190000048
式中，
Figure BDA0002940491190000049
为t1时段第t1时刻第n个站点所在细分辨率网格第m颗卫星的降水值。
pm(PFusion|Pm,FR,obs)为在不同卫星细分辨率降水、站点降水发生的条件下融合降水的后验概率，p(PFusion|obs)为站点降水发生的条件下细分辨率降水发生的概率；ωm为第m颗卫星细分辨率降水的权重；PFusion为t1时段内站点所在细分辨率网格的融合降水，具体表示为：
Figure BDA0002940491190000051
式中，
Figure BDA0002940491190000052
为t1时段第t1时刻第n个站点所在细分辨率网格的融合降水值。
步骤3.2，在自适应动态贝叶斯框架下，采用最大期望方法求权重的前提是数据集符合正态分布，将单颗卫星的细分辨率降水数据集转化为正态分布。
对上式求期望和方差：
Figure BDA0002940491190000053
Figure BDA0002940491190000054
式中，
Figure BDA0002940491190000055
为第m颗卫星细分辨率降水的方差。
自适应动态贝叶斯的权重ωm与Pm,FR直接相关，假设每颗卫星的初始权重都相等即
Figure BDA0002940491190000056
通过最大期望算法来求解权重ωm。最大期望算法的基本假设是pm(PFusion|Pm,FR,obs)符合正态分布，因此需要通过Box-Cox变换进行转换：
Figure BDA0002940491190000057
式中，Pm,γ为转换之后的降尺度降水数据，γ为转换参数(一般用极大似然方法确定γ的取值)。
步骤3.3，利用最大期望算法计算权重ωm，得到单颗卫星细分辨率降水所在站点的最优权重，自适应动态贝叶斯框架的参数集可表示为：
Figure BDA0002940491190000058
式中，δ为参数集，ωm为第m颗卫星细分辨率降水的最优权重，
Figure BDA0002940491190000059
为第m颗卫星细分辨降水的方差。
采用极大似然方法求解，对上述公式(5)取对数，并求导同时导数等于0即期望最大，反推求出权重ωm：
Figure BDA00029404911900000510
Figure BDA00029404911900000511
根据公式(11)反推得到ωm。
前述的一种星地多源降水自适应动态融合方法，其特征在于，所述步骤4中基于克里金插值计算每颗卫星每个细分辨率网格的最优权重，每颗卫星细分辨率降水加权求和得到融合降水，具体步骤包括：
步骤4.1，用克里金插值方法将每颗卫星在不同站点的最优权重(ωm)插值到每个网格，得到t1+1时刻每颗卫星每个细分辨率网格的权重，继续上述步骤，得到t1+2时刻每颗卫星每个细分辨率网格的权重；
步骤4.2，通过加权求和得到融合降水产品。
本发明的有益效果：本发明提供的一种星地多源降水自适应动态融合方法，首先构建混合地理加权逐步回归模型对卫星降水进行降尺度，得到时空一致的细分辨率卫星降水；构建自适应动态贝叶斯框架，获取在每个时刻每颗卫星每个细分辨率网格的权重，乘以细分辨率卫星降水相加得到融合降水产品。本方法以多颗卫星降水、细分辨率的时不变因子、时变因子为基础数据，数据来源稳定可靠，方法中变量之间的函数关系明确，建立了卫星降水与多种因子之间的关系和自适应动态权重求解方法，保证了结果的客观合理性；弥补了卫星降水和站点降水各自的缺点，获取了细分辨率的融合降水产品，作为精细化水文模型的输入。
附图说明
图1是本发明的计算流程示意图；
图2是粗分辨和细分辨率因子分布图(以风速为例)；
图3是每颗卫星降水空间分布图；
图4是每颗卫星细分辨降水分布图；
图5是融合降水分布图。
具体实施方式
下面结合附图和具体实施例对本发明作进一步描述。
如图1所示，本发明提供的一种星地多源降水自适应动态融合方法，包括以下步骤：
步骤1，获取流域粗分辨率(Coarse Resolution，CR)和细分辨率(FineResolution，FR)的时不变因子，主要包括：DEM(DEMCR、DEMFR)、坡度(SlopeCR、SlopeFR)、Hand指数(HandCR、HandFR)、粗糙度(SRCR、SRFR)；粗分辨率和细分辨率的时变因子分别为：NDVI(NDVICR、NDVIFR)、风速(WSCR、WSFR)、地表温度(LSTCR、LSTFR)、散射指数(SICR、SIFR)、微波极化指数(PDCR、PDFR)；卫星编号用m表示，m＝1,2，……，5分别代表GPM、TRMM、CMORPH、PERSIANN、PERSIANN_CDR；粗分辨率网格编号用i表示，细分辨率网格编号用j表示，站点个数用n表示；如图2所示。
步骤2中基于混合地理加权逐步回归模型对多颗卫星进行降尺度，得到每颗卫星细分辨率降水，具体包括以下步骤：
步骤2.1，基于粗分辨率的时不变因子和时变因子，构建混合地理加权逐步回归模型(以某个时刻某颗卫星单个粗分辨率网格为例)：
Figure BDA0002940491190000071
式中，i为粗分辨率网格编号，m为卫星编号，t表示时间，CR表示粗分辨率，
Figure BDA0002940491190000072
为t时刻m卫星的第i个粗分辨率网格的降水(如图3所示)；
Figure BDA0002940491190000073
为第i个粗分辨率网格的高程，
Figure BDA0002940491190000074
为
Figure BDA0002940491190000075
对应的时不变系数；
Figure BDA0002940491190000076
为第i个粗分辨率网格的坡度，
Figure BDA0002940491190000077
为
Figure BDA0002940491190000078
对应的时不变系数；
Figure BDA0002940491190000079
为第i个粗分辨率网格的Hand指数，
Figure BDA00029404911900000710
为
Figure BDA00029404911900000711
对应的时不变系数；
Figure BDA00029404911900000712
为第i个粗分辨率网格的地表粗糙度，
Figure BDA00029404911900000713
为
Figure BDA00029404911900000714
对应的时不变系数；
Figure BDA00029404911900000715
为t时刻第i个粗分辨率网格的植被归一化指数，
Figure BDA00029404911900000716
为
Figure BDA00029404911900000717
对应的时变系数；
Figure BDA00029404911900000718
为t时刻第i个粗分辨率网格的风速，
Figure BDA00029404911900000719
为
Figure BDA00029404911900000720
对应的时变系数；
Figure BDA00029404911900000721
为t时刻第i个粗分辨率网格的地表温度，
Figure BDA00029404911900000722
为
Figure BDA00029404911900000723
对应的时变系数；
Figure BDA00029404911900000724
为t时刻第i个粗分辨率网格的散射指数，
Figure BDA00029404911900000725
为
Figure BDA00029404911900000726
对应的时变系数；
Figure BDA00029404911900000727
为t时刻第i个粗分辨率网格的微波极化指数，
Figure BDA00029404911900000728
为
Figure BDA00029404911900000729
对应的时变因子；
Figure BDA00029404911900000730
为t时刻第i个粗分辨率网格的残差。
步骤2.2，方程(1)进行逐步回归，逐步引入对降水贡献量大的因子：
逐步引入9个粗分辨率因子，包括：4个时不变因子
Figure BDA00029404911900000731
和5个时变因子
Figure BDA00029404911900000732
进行逐步回归分析和对应的回归系数进行F检验，若检验通过给定显著性水平α，则将该因子保留，逐步回归完成以后留下的因子记为
Figure BDA00029404911900000733
Figure BDA00029404911900000734
为4个时不变因子和5个时变因子的子集，由两个部分组成分别为：
Figure BDA00029404911900000735
(时不变因子子集)、
Figure BDA00029404911900000736
(时变因子子集)。
步骤2.3，假设上步所得
Figure BDA00029404911900000737
由n1个时不变因子和n2个时变因子组成，某个时刻单个细分辨率网格的混合地理加权逐步回归模型为：
Figure BDA00029404911900000738
式中，j为细分辨率网格编号，m为卫星编号，t表示时间，FR表示细分辨率，
Figure BDA00029404911900000739
为t时刻第m颗卫星的第j个细分辨率网格的降水；
Figure BDA00029404911900000740
为第j个细分辨率网格的第1个时不变因子，
Figure BDA00029404911900000741
为其对应的系数；
Figure BDA00029404911900000742
为第j个细分辨率网格的第n1个时不变因子，
Figure BDA00029404911900000743
为其对应的系数；
Figure BDA00029404911900000744
为第j个细分辨率网格的第1个时变因子，
Figure BDA0002940491190000081
为其对应的系数；
Figure BDA0002940491190000082
为第j个细分辨率网格的第n2个时变因子，
Figure BDA0002940491190000083
为其对应的系数；
Figure BDA0002940491190000084
为t时刻第j个细分辨率网格的残差。
对于t时刻细分辨率单个网格系数求解如下：
Figure BDA0002940491190000085
式中，
Figure BDA0002940491190000086
为时不变因子系数矩阵；
Figure BDA0002940491190000087
为时变因子系数矩阵；
Figure BDA0002940491190000088
为时不变因子矩阵，
Figure BDA0002940491190000089
为时变因子矩阵，I为单位矩阵，W为权重矩阵，ym,CR为第m颗原始卫星降水矩阵，
Figure BDA00029404911900000810
计算t时刻第m颗卫星第j个细分辨率网格的降水
Figure BDA00029404911900000811
(即降尺度后的降水)：
Figure BDA00029404911900000812
式中，sum定义为矩阵的逐项相加函数，
Figure BDA00029404911900000813
作为系统残差可取常数，
Figure BDA00029404911900000814
分别为：
Figure BDA00029404911900000815
通过j次循环得到t时刻第m颗卫星降尺度之后的降水，通过m次循环得到每颗卫星降尺度降水即细分辨率时空一致的降水，如图4所示。
步骤3，通过降水分布转换、降水概率密度函数参数优化来计算单颗卫星每个站点所在细分辨率网格的最优权重，具体包括以下步骤：
步骤3.1，基于不同卫星细分辨率降水、站点降水条件下融合降水的后验分布，根据全概率公式，建立自适应动态贝叶斯框架，融合降水的概率密度可表示为：
Figure BDA00029404911900000816
式中，t1为训练时间段(t1的取值范围为一个动态范围，可为同一个季节内的0-90天，获取t1+1时刻的权重以t1时间段训练数据为基础)；obs为t1时段内实测站点降水矩阵，
Figure BDA00029404911900000817
其中
Figure BDA00029404911900000818
代表t1时段第一个时刻的n个站点降水列向量，
Figure BDA00029404911900000819
代表t1时段第二个时刻的n个站点降水列向量，
Figure BDA00029404911900000820
代表t1时段第t1时刻的n个站点降水列向量，具体表示如下：
Figure BDA0002940491190000091
式中，
Figure BDA0002940491190000092
为t1时段第t1时刻第n个站点的降水值。
Pm,FR为t1时段内站点对应的第m颗卫星细分辨率降水，表示如下：
Figure BDA0002940491190000093
式中，
Figure BDA0002940491190000094
为t1时段第t1时刻第n个站点所在细分辨率网格第m颗卫星的降水值。
pm(PFusion|Pm,FR,obs)为在不同卫星细分辨率降水、站点降水发生的条件下融合降水的后验概率，p(PFusion|obs)为站点降水发生的条件下细分辨率降水发生的概率；ωm为第m颗卫星细分辨率降水的权重；PFusion为t1时段内站点所在细分辨率网格的融合降水，具体表示为：
Figure BDA0002940491190000095
式中，
Figure BDA0002940491190000096
为t1时段第t1时刻第n个站点所在细分辨率网格的融合降水值。
步骤3.2，在自适应动态贝叶斯框架下，采用最大期望方法求权重的前提是数据集符合正态分布，将单颗卫星的细分辨率降水数据集转化为正态分布。
对上式求期望和方差：
Figure BDA0002940491190000097
Figure BDA0002940491190000098
式中，
Figure BDA0002940491190000099
为第m颗卫星细分辨率降水的方差。
自适应动态贝叶斯的权重ωm与Pm,FR直接相关，假设每颗卫星的初始权重都相等即
Figure BDA0002940491190000101
通过最大期望算法来求解权重ωm。最大期望算法的基本假设是pm(PFusion|Pm,FR,obs)符合正态分布，因此需要通过Box-Cox变换进行转换：
Figure BDA0002940491190000102
式中，Pm,γ为转换之后的降尺度降水数据，γ为转换参数(一般用极大似然方法确定γ的取值)。
步骤3.3，利用最大期望算法计算权重ωm，得到单颗卫星细分辨率降水所在站点的最优权重，自适应动态贝叶斯框架的参数集可表示为：
Figure BDA0002940491190000103
式中，δ为参数集，ωm为第m颗卫星细分辨率降水的最优权重，
Figure BDA0002940491190000104
为第m颗卫星细分辨降水的方差。
采用极大似然方法求解，对上述公式(5)取对数，并求导同时导数等于0即期望最大，反推求出权重ωm：
Figure BDA0002940491190000105
Figure BDA0002940491190000106
根据公式(11)反推得到ωm。
步骤4，基于克里金插值计算每颗卫星每个细分辨率网格的最优权重，每颗卫星细分辨率降水加权求和得到融合降水，具体步骤包括：
步骤4.1，用克里金插值方法将每颗卫星在不同站点的最优权重(ωm)插值到每个网格，得到t1+1时刻每颗卫星每个细分辨率网格的权重，继续上述步骤，得到t1+2时刻每颗卫星每个细分辨率网格的权重；
步骤4.2，通过加权求和得到融合降水产品，如图5所示。
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (2)

1.一种星地多源降水自适应动态融合方法，其特征在于，包括以下步骤：
步骤1，获取流域粗分辨率和细分辨率的时不变因子和时变因子；时不变因子包括：DEM、坡度、Hand指数、地表粗糙度；时变因子包括：NDVI、风速、地表温度、散射指数及微波极化指数；
步骤2，基于混合地理加权逐步回归模型对卫星降水进行降尺度，得到每颗卫星细分辨率降水；
步骤2.1，基于粗分辨率的时不变因子和时变因子，构建某个时刻某颗卫星单个粗分辨率网格的混合地理加权逐步回归模型：
Figure FDA0003195744820000011
式中，i为粗分辨率网格编号，m为卫星编号，t表示时间，CR表示粗分辨率，
Figure FDA0003195744820000012
为t时刻第m颗卫星第i个粗分辨率网格的降水；
Figure FDA0003195744820000013
为第i个粗分辨率网格的高程，
Figure FDA0003195744820000014
为
Figure FDA0003195744820000015
对应的时不变系数；
Figure FDA0003195744820000016
为第i个粗分辨率网格的坡度，
Figure FDA0003195744820000017
为
Figure FDA0003195744820000018
对应的时不变系数；
Figure FDA0003195744820000019
为第i个粗分辨率网格的Hand指数，
Figure FDA00031957448200000110
为
Figure FDA00031957448200000111
对应的时不变系数；
Figure FDA00031957448200000112
为第i个粗分辨率网格的地表粗糙度，
Figure FDA00031957448200000113
为
Figure FDA00031957448200000114
对应的时不变系数；
Figure FDA00031957448200000115
为t时刻第i个粗分辨率网格的植被归一化指数，
Figure FDA00031957448200000116
为
Figure FDA00031957448200000117
对应的时变系数；
Figure FDA00031957448200000118
为t时刻第i个粗分辨率网格的风速，
Figure FDA00031957448200000119
为
Figure FDA00031957448200000120
对应的时变系数；
Figure FDA00031957448200000121
为t时刻第i个粗分辨率网格的地表温度，
Figure FDA00031957448200000122
为
Figure FDA00031957448200000123
对应的时变系数；
Figure FDA00031957448200000124
为t时刻第i个粗分辨率网格的散射指数，
Figure FDA00031957448200000125
为
Figure FDA00031957448200000126
对应的时变系数；
Figure FDA00031957448200000127
为t时刻第i个粗分辨率网格的微波极化指数，
Figure FDA00031957448200000128
为
Figure FDA00031957448200000129
对应的时变因子；
Figure FDA00031957448200000130
为t时刻第i个粗分辨率网格的残差；
步骤2.2，对上步构建的混合地理加权逐步回归模型进行逐步回归，逐步引入对降水贡献量大的因子：
逐步引入9个粗分辨率因子，包括：4个时不变因子和5个时变因子，4个时不变因子为
Figure FDA00031957448200000131
5个时变因子为
Figure FDA00031957448200000132
Figure FDA0003195744820000021
进行逐步回归分析和对应的回归系数进行F检验，若检验通过给定显著性水平α，则将该因子保留，逐步回归完成以后留下的因子记为
Figure FDA0003195744820000022
Figure FDA0003195744820000023
为4个时不变因子和5个时变因子的子集，由两个部分组成分别为：时不变因子子集
Figure FDA0003195744820000024
时变因子子集
Figure FDA0003195744820000025
步骤2.3，假设上步所得
Figure FDA0003195744820000026
由n1个时不变因子和n2个时变因子组成，某个时刻单个细分辨率网格的混合地理加权逐步回归模型为：
Figure FDA0003195744820000027
式中，j为细分辨率网格编号，m为卫星编号，t表示时间，FR表示细分辨率，
Figure FDA0003195744820000028
为t时刻第m颗卫星第j个细分辨率网格的降水；
Figure FDA0003195744820000029
为第j个细分辨率网格的第1个时不变因子，
Figure FDA00031957448200000210
为其对应的系数；
Figure FDA00031957448200000211
为第j个细分辨率网格的第n1个时不变因子，
Figure FDA00031957448200000212
为其对应的系数；
Figure FDA00031957448200000213
为第j个细分辨率网格的第1个时变因子，
Figure FDA00031957448200000214
为其对应的系数；
Figure FDA00031957448200000215
为第j个细分辨率网格的第n2个时变因子，
Figure FDA00031957448200000216
为其对应的系数；
Figure FDA00031957448200000217
为t时刻第j个细分辨率网格的残差；
对于t时刻细分辨率单个网格系数求解如下：
Figure FDA00031957448200000218
式中，
Figure FDA00031957448200000219
为时不变因子系数矩阵；
Figure FDA00031957448200000220
为时变因子系数矩阵；
Figure FDA00031957448200000221
为时不变因子矩阵，
Figure FDA00031957448200000222
为时变因子矩阵，I为单位矩阵，W为权重矩阵，ym,CR为第m颗卫星降水矩阵，
Figure FDA00031957448200000223
计算t时刻第m颗卫星第j个细分辨率网格的降水
Figure FDA00031957448200000224
Figure FDA00031957448200000225
式中，sum定义为矩阵的逐项相加函数，
Figure FDA00031957448200000226
作为系统残差可取常数，
Figure FDA00031957448200000227
Figure FDA0003195744820000031
分别为：
Figure FDA0003195744820000032
通过对每个细分辨率网格j逐一循环得到t时刻第m颗卫星降尺度之后的降水，通过对每个卫星m逐一循环得到每颗卫星降尺度降水即细分辨率时空一致的降水；
步骤3，构建自适应动态贝叶斯框架，通过降水误差的分布转换、降水概率密度函数参数优化来计算每颗卫星在降水观测站点所在细分辨率网格的最优权重；
步骤3.1，基于不同卫星细分辨率降水、降水观测站点降水条件下融合降水的后验分布，根据全概率公式，建立自适应动态贝叶斯框架，融合降水的概率密度表示为：
Figure FDA0003195744820000033
式中，t1为训练时间段；obs为t1时段内实测降水观测站点降水矩阵，
Figure FDA0003195744820000034
其中
Figure FDA0003195744820000035
代表t1时段第一个时刻的n个降水观测站点降水列向量，
Figure FDA0003195744820000036
代表t1时段第二个时刻的n个降水观测站点降水列向量，
Figure FDA0003195744820000037
代表t1时段第t1时刻的n个降水观测站点降水列向量，具体表示如下：
Figure FDA0003195744820000038
式中，
Figure FDA0003195744820000039
为t1时段第t1时刻第n个降水观测站点的降水值；
Pm,FR为t1时段内降水观测站点对应的第m颗卫星细分辨率降水，表示如下：
Figure FDA0003195744820000041
式中，
Figure FDA0003195744820000042
为t1时段第t1时刻第n个降水观测站点所在细分辨率网格第m颗卫星的降水值；pm(PFusion|Pm,FR,obs)为在不同卫星细分辨率降水、降水观测站点降水发生的条件下融合降水的后验概率，p(PFusion|obs)为降水观测站点降水发生的条件下细分辨率降水发生的概率；ωm为第m颗卫星在降水观测站点所在细分辨率网格的最优权重；PFusion为t1时段内降水观测站点所在细分辨率网格的融合降水，具体表示为：
Figure FDA0003195744820000043
式中，
Figure FDA0003195744820000044
为t1时段第t1时刻第n个降水观测站点所在细分辨率网格的融合降水值；
步骤3.2，对降水观测站点降水发生的条件下不同卫星细分辨率降水发生的概率求期望和方差：
Figure FDA0003195744820000045
Figure FDA0003195744820000046
式中，
Figure FDA0003195744820000047
为第m颗卫星细分辨率降水的方差；
通过Box-Cox变换进行转换：
Figure FDA0003195744820000051
式中，Pm,γ为转换之后的降尺度降水数据，γ为转换参数；
步骤3.3，利用最大期望算法计算得到单颗卫星降水所在降水观测站点的最优权重ωm：
Figure FDA0003195744820000052
式中，δ为自适应动态贝叶斯框架的参数集；
Figure FDA0003195744820000053
式中，ωm为第m颗卫星在降水观测站点所在细分辨率网格的最优权重，
Figure FDA0003195744820000054
为第m颗卫星细分辨率降水的方差；
步骤4，基于克里金插值通过降水观测站点所在细分辨率网格求得的最优权重，计算每颗卫星每个细分辨率网格的最优权重，每颗卫星细分辨率降水加权求和得到融合降水。
2.根据权利要求1所述的一种星地多源降水自适应动态融合方法，其特征在于，所述步骤4包括：
步骤4.1，用克里金插值方法将第m颗卫星在不同降水观测站点所在细分辨率网格的最优权重ωm插值到每个网格，得到t1+1时刻每颗卫星每个细分辨率网格的最优权重，继续用克里金插值方法将t1+1时刻第m颗卫星在不同降水观测站点所在细分辨率网格的最优权重ωm插值到每个网格，得到t1+2时刻每颗卫星每个细分辨率网格的最优权重；
步骤4.2，通过加权求和得到融合降水：该方法有效的将卫星降尺度方法与多源降水融合方法进行结合，最终产生了细分辨率的降水融合产品。
CN202110175249.6A 2021-02-09 2021-02-09 一种星地多源降水自适应动态融合方法 Active CN112861072B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110175249.6A CN112861072B (zh) 2021-02-09 2021-02-09 一种星地多源降水自适应动态融合方法
NL2030246A NL2030246A (en) 2021-02-09 2021-12-22 Satellite-ground multi-source precipitation adaptive dynamic fusion method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110175249.6A CN112861072B (zh) 2021-02-09 2021-02-09 一种星地多源降水自适应动态融合方法

Publications (2)

Publication Number Publication Date
CN112861072A CN112861072A (zh) 2021-05-28
CN112861072B true CN112861072B (zh) 2021-10-19

Family

ID=75989390

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110175249.6A Active CN112861072B (zh) 2021-02-09 2021-02-09 一种星地多源降水自适应动态融合方法

Country Status (2)

Country Link
CN (1) CN112861072B (zh)
NL (1) NL2030246A (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115203639B (zh) * 2022-06-21 2023-03-10 中国长江三峡集团有限公司 一种基于矩阵运算的不规则网格面雨量计算方法及系统
CN115357847B (zh) * 2022-08-28 2023-09-19 水利部交通运输部国家能源局南京水利科学研究院 一种基于误差分解的日尺度星地降水融合方法
CN115755221B (zh) * 2022-10-22 2024-08-23 天津大学 一种基于数理不确定性分析的多源降水数据融合方法
CN115795399B (zh) * 2023-01-31 2023-04-28 中国科学院地理科学与资源研究所 一种多源遥感降水数据自适应融合方法和系统
CN117349795B (zh) * 2023-12-04 2024-02-02 水利部交通运输部国家能源局南京水利科学研究院 基于ann和gwr耦合的降水融合方法和系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107918166A (zh) * 2016-10-09 2018-04-17 清华大学 多卫星融合降水方法及系统
CN107918165A (zh) * 2016-10-09 2018-04-17 清华大学 基于空间插值的多卫星融合降水预测方法及系统
CN108647740A (zh) * 2018-05-16 2018-10-12 河海大学 利用高分辨率地形和气象因子进行多源降水融合的方法
CN108681639A (zh) * 2018-05-16 2018-10-19 河海大学 一种考虑局部变量和全局变量的动态降水降尺度方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107918166A (zh) * 2016-10-09 2018-04-17 清华大学 多卫星融合降水方法及系统
CN107918165A (zh) * 2016-10-09 2018-04-17 清华大学 基于空间插值的多卫星融合降水预测方法及系统
CN108647740A (zh) * 2018-05-16 2018-10-12 河海大学 利用高分辨率地形和气象因子进行多源降水融合的方法
CN108681639A (zh) * 2018-05-16 2018-10-19 河海大学 一种考虑局部变量和全局变量的动态降水降尺度方法

Also Published As

Publication number Publication date
CN112861072A (zh) 2021-05-28
NL2030246A (en) 2022-04-29

Similar Documents

Publication Publication Date Title
CN112861072B (zh) 一种星地多源降水自适应动态融合方法
Xu et al. Do ERA5 and ERA5-land precipitation estimates outperform satellite-based precipitation products? A comprehensive comparison between state-of-the-art model-based and satellite-based precipitation products over mainland China
US11333796B2 (en) Spatial autocorrelation machine learning-based downscaling method and system of satellite precipitation data
CN111401602B (zh) 基于神经网络的卫星以及地面降水测量值同化方法
CN110222911B (zh) 一种卫星遥感与地面数据协同的雨量站网优化布局方法
Conti et al. Evaluation and comparison of satellite precipitation estimates with reference to a local area in the Mediterranean Sea
Yang et al. Merging high‐resolution satellite‐based precipitation fields and point‐scale rain gauge measurements—A case study in Chile
Soltani et al. Review of assimilating GRACE terrestrial water storage data into hydrological models: Advances, challenges and opportunities
Aryee et al. Development of high spatial resolution rainfall data for Ghana
Islam Statistical comparison of satellite-retrieved precipitation products with rain gauge observations over Bangladesh
Wang et al. Performance of three reanalysis precipitation datasets over the Qinling‐Daba Mountains, eastern fringe of Tibetan Plateau, China
CN115795399B (zh) 一种多源遥感降水数据自适应融合方法和系统
CN114417646B (zh) 一种高维异构降水数据融合方法及系统
CN112329265A (zh) 卫星遥感降水精细化空间估算方法及系统
Chen et al. Sub-daily soil moisture estimate using dynamic Bayesian model averaging
CN114461983B (zh) 一种基于水量平衡原理的卫星降水产品空间降尺度方法
Fu et al. Implications of diurnal variations in land surface temperature to data assimilation using MODIS LST data
Zhang et al. Merging multisatellite precipitation products using stacking method and the censored-shifted gamma ensemble model output statistics in china's Beimiaoji basin
Jing et al. Cascaded downscaling–calibration networks for satellite precipitation estimation
Hisam et al. Comprehensive evaluation of Satellite-Based and reanalysis precipitation products over the Mediterranean region in Turkey
Chen et al. Establishment and assessment of hourly high‐resolution gridded air temperature data sets in Zhejiang, China
CN117493475A (zh) 一种基于机器学习重构缺资料地区月径流的方法及系统
Zhou et al. A station-data-based model residual machine learning method for fine-grained meteorological grid prediction
Du et al. Real-time tropospheric delay map retrieval using sparse GNSS stations
CN116523189A (zh) 一种考虑水文特征的墒情站点规划方法、装置及存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant