基于虚拟站点的格状遥测雨量和点状实测雨量的融合框架

1.本发明属于少资料地区降雨估计技术领域，涉及一种基于虚拟站点的格状遥测雨量和点状实测雨量的融合框架。


背景技术：

2.降水在水文、农业、气象等领域至关重要，但降雨场较强的空间异质性导致对其精确估计成为一项挑战。目前降水监测主要包括两种方式，一种是基于地面雨量站的点状直接测量方式，另一种是基于雷达、卫星和数值气象模型等的网格状间接测量方式。直接测量方式在站点位置处测量精准，但受投资资金、地形气候以及设备维护等条件限制，在站点空间覆盖稀疏的地区，雨量站点难以捕捉降水的全部空间分布信息，例如漏测降雨中心。得力于遥感技术快速发展，间接监测方法的空间连续估计可反映降雨场的空间分布，特别是对于空间降雨最值的捕捉；然而，间接测量方法对降雨量的数值估计存在误差和较多不确定性，特别是对于近实时或实时降水产品。因此，亟需发明一项融合格状遥测雨量和点状实测雨量的降雨估计方法，从而更准确地估计降雨场。
3.根据两种测量方式的特性，融合思路主要分为两类，一类是利用点状实测雨量校正格状遥测雨量；另一类是利用格状遥测雨量补充点状实测雨量缺测的降雨空间分布信息。许多研究在融合格状遥测雨量和点状实测雨量方面做出尝试，融合方法主要包含三类。第一类，利用统计分析模型以点状实测雨量为基准校正格状遥测雨量；第二类，给定多源降雨估计以不同权重从而融合多源降雨估计；第三类，基于机器学习方法融合多源降雨估计。以上融合方法的步骤可以概括为：首先，在有雨量站的格点，以地面站点的实测雨量为基准，利用不同回归分析方法，训练或校准不同的转化模型；其次，在无雨量站的格点，利用训练后的转化模型估计降雨。以上方法采用了逐格融合方式，而未考虑格与格间的空间关系，可以看作非空间类融合方法。因此，这三类方法在融合过程中只遵循了第一条融合思路，而忽略了第二条融合思路。
4.相较于仅用点状实测雨量进行空间插值估计降雨场，上述融合方法虽然在雨量站点分布极稀疏的地区可以提高降雨估计精度，但由于第二类融合思路的缺失，随着估计降雨场的目标区域的雨量站点密度增加，这种精度的提高越来越弱，当站点密度超过某一阈值时，融合两类降雨的性能会劣于仅用点状实测雨量空间插值。在wu等人的研究中，该阈值约为60站/106km2，与世界气象组织推荐的不少于1000站/106km2相比，站点仍极为稀疏，难以准确捕捉降水空间分布。因此，亟需发明一项在任何雨量站点密度条件下都适用的格状遥测雨量和点状实测雨量融合方法。


技术实现要素：

5.针对现有技术存在的问题，本发明提供一种基于虚拟站点的格状遥测雨量和点状实测雨量的融合框架。
6.本发明的技术方案如下：
7.一种基于虚拟站点的格状遥测雨量和点状实测雨量的融合框架，包括以下步骤：
8.第一步，根据雨量站点集合g＝[g1,g2,
…
]，将待估计降雨场的对象地区划分为多个子区域。
[0009]
首先，将对象区域按选定分辨率划分为网格；其次，将由同一组雨量站点控制的网格划分为一个子区域，其中雨量站点的可控距离在10～50km选取。如图1所示，以空间中任意两个雨量站点g1和g2为例，dc是雨量站点的可控距离，图中以不同颜色区分不同子区域，g
x
代表控制同一子区域的雨量站，它是所有雨量站点的一个子集。
[0010]
对于一个子区域，任一格点x的估计降雨的理论取值范围均是：
[0011][0012]
其中，
[0013][0014][0015]
其中，是实际雨量站点gi的实测雨量值，gi是格点x控制距离内的实际雨量站点。
[0016]
第二步，针对每一个子区域，推荐所需虚拟站点的数量和位置。其中，虚拟站点的数量和位置由多源格状遥测雨量产品与点状实测雨量做空间比较后推荐。
[0017]
实际雨量超出第一步中所述的理论取值范围的格点均可选取为待补充的虚拟站点。为有效扩宽的理论取值范围，选取超出理论取值范围的最大雨量和最小雨量分别所在的两个格点作为该子区域的虚拟站点。值得注意的是，某子区域可能存在的情况包括：存在部分格点实际雨量小于p
min
，同时存在部分格点实际雨量大于p
max
，则该子区域需在最大雨量和最小雨量分别所在的两个格点补充虚拟站；不存在格点的实际雨量小于p
min
，存在部分格点实际雨量大于p
max
，则该子区域需在最大雨量所在格点补充一个虚拟站；存在格点的实际雨量小于p
min
，不存在格点的实际雨量大于p
max
，则该子区域需在最小雨量所在格点补充一个虚拟站；所有格点的实际雨量均在[p
min
,p
max
]范围内，则该子区域不需补充虚拟站。因此，在每个子区域内需补充虚拟站数量不超过两个。
[0018]
其中，因实际雨量未知，实际雨量与理论取值范围边界的比较需借助多源格状遥测雨量产品，具体步骤如下：
[0019]
将每个格状遥测雨量产品在单个时间步长估计的降雨场做0-1归一化处理，得到归一化的降雨场则在点x：
[0020][0021]
式中，是在单个时间步长内产品s估计的降雨场，是在单个时间步长内产品s在格点x估计的降雨量。
[0022]
对于一个子区域，如果，
[0023][0024]
那么，的位置就是产品s推荐的最大雨量虚拟站点的位置。反之，产品s未
推荐最大雨量虚拟站点。
[0025]
式中，是一个子区域位置处的x＝[x1,x2,
…
,xm]是子区域所包含的格点，是控制该子区域的站点位置处的a
max
是一个增强不等式显著性的参数，建议取值为0.1-0.2之间。
[0026]
同理，如果，
[0027][0028]
那么，的位置就是产品s推荐的最小雨量虚拟站点的位置。反之，产品s未推荐最小雨量虚拟站点。
[0029]
式中，a
min
是一个增强不等式显著性的参数，建议取值为0-0.1之间。
[0030]
所有降水产品均完成推荐虚拟站点的步骤后，若超过半数产品推荐的虚拟站点位置相同或者相邻，则该位置为最终确定的虚拟站点位置，反之，该子区域不需补充虚拟站点。
[0031]
第三步，利用现有的非空间类融合方法估计虚拟站点处的降雨量，所述的非空间类融合方法的思路可概括为：
[0032]
2.1)训练模型
[0033]
基于各类回归分析建立一个因变量与自变量的转化模型，其中，因变量为站点实测雨量，自变量为多源降雨产品在实际站点位置的雨量估计值以及其他相关的协变量。
[0034][0035]
式中，f(
·
)是所训练的转化模型，pg是实际站点g的实测雨量，是降雨产品si在g处的雨量估计值，ag是与g相关的协变量集合，如地形等信息，ε是模型误差。
[0036]
2.2)估计降雨
[0037]
基于训练后的转化模型估计虚拟站点处的降雨量因变量为虚拟站点处的降雨量，自变量为多源降雨产品在虚拟站点位置的雨量估计值以及其他相关的协变量。
[0038][0039]
式中，降雨产品si在虚拟站点v处的雨量估计值，av是与v相关的协变量。
[0040]
第四步，判断目标区域的所有格点距最近雨量站点的距离是否小于第一步中选定的雨量站点可控距离，其中雨量站点是实际站点和虚拟站点的总集合g＝[g1,g2,
…
,v1,v2,
…
]。若是，则转到第五步；若否，则以站点总集合作为更新后的实际站点集合并转到第一步，如图3所示。
[0041]
第五步，基于现有的空间插值方法，利用实际站点和虚拟站点进行空间插值，得到降雨场。
[0042]
本发明的效果和益处为：
[0043]
本发明利用多源格状遥测雨量产品推测虚拟站点的数量、位置和雨量，进而整合实际站点和虚拟站点进行雨量空间插值，从而补充了雨量站点缺测的降雨空间分布信息，并削弱了格状遥测产品的不确定性带来的干扰。本发明兼并多源产品融合以及站点空间插值两种方法的优势，为获得更高精度的降水空间分布提供了新思路。
附图说明
[0044]
图1是本发明涉及的子区域划分流程图，图中字母含义与公式中涉及字母的含义相同。
[0045]
图2是本发明涉及的推断虚拟站点位置的流程。
[0046]
图3是本发明的运行框架。
具体实施方式
[0047]
本发明在多源降水产品融合与站点空间插值的基础上，提出一种基于虚拟站点的格状遥测雨量和点状实测雨量的融合框架。
[0048]
下面通过实施例，对本发明做进一步说明。
[0049]
选用中国北部某区域作为实例应用区域，面积为31352km2，海拔为3-2400米。该区域具有典型的温带季风气候，冬季干燥寒冷，夏季多雨炎热，春秋季短。年降水量约为1100mm，降雨量年际变化较大，年内分配不均，降水集中在5～9月。该地区雨量站点密集，共有624个雨量站，几乎每0.1
°
网格至少存在一个雨量站，因此，能提供每天的降雨场真值作为精度评价的依据。
[0050]
具体步骤如下：
[0051]
第一步，根据雨量站点将待估计降雨场的对象地区划分为多个子区域。
[0052]
首先，将对象区域划分为0.1
°
网格。对象区域按照0.1
°
分辨率可划分为322个网格。若多个雨量站点位于同一网格，取多雨量站实测值的平均值作为该格点雨量真值。若网格内无雨量站，取该网格周围最临近的四个网格的雨量平均值作为该格点的雨量真值。若网格内无雨量站且最邻近网格内也无雨量站，则该网格雨量记为缺测。在对象区域内选择10个缺测率小于2％的格点假定为少资料地区的实际雨量站点，其他格点值用于评价降雨场的估计精度，则实际雨量站点集合记为g＝[g1,g2,
…g10
]。
[0053]
其次，将由同一组雨量站点控制的网格划分为一个子区域，其中雨量站点的可控距离50km。如图1所示，以空间中任意两个雨量站点g1和g2为例，dc是雨量站点的可控距离，图中以不同颜色区分不同子区域，g
x
代表控制同一子区域的雨量站，它是所有雨量站点的一个子集。
[0054]
对于一个子区域，任一格点x的估计降雨的理论取值范围均是：
[0055][0056]
其中，
[0057][0058][0059]
其中，是实际雨量站点gi的实测雨量值，gi是格点x控制距离内的实际雨量站点。
[0060]
第二步，针对每一个子区域，推荐所需虚拟站点的数量和位置。其中，虚拟站点的数量和位置由多源格状遥测雨量产品与点状实测雨量做空间比较后推荐。
[0061]
选用5种基于卫星的降水产品和3种基于数值气象模型的降水产品为实例驱模型，即gsmap-v6-nrt(gsmap-n)、gsmap-v6-gauge-nrt(gsmap-gn)、merg-v6-early(imerg-e)、
imerg-v6-late(imerg-l)、persiann-ccs、ecmwf、jma以及ncep，表1为选用的降水产品的基本信息。本发明实施实例的降水产品统一采用0.1
°
的空间分辨率，对于初始空间分辨率不同的降水产品，采用重采样的方式将其精度转化为0.1
°
。
[0062]
表1选用的格状遥测降水产品的基本信息
[0063][0064][0065]
将每个格状遥测雨量产品在单个时间步长估计的降雨场做0-1归一化处理，得到归一化的降雨场则在点x：
[0066][0067]
式中，是在单个时间步长内产品s估计的降雨场，是在单个时间步长内产品s在格点x估计的降雨量。
[0068]
对于一个子区域，如果，
[0069][0070]
那么，的位置就是产品s推荐的最大雨量虚拟站点的位置。反之，产品s未推荐最大雨量虚拟站点。
[0071]
式中，是一个子区域位置处的x＝[x1,x2,
…
,xm]是子区域所包含的格点，是控制该子区域的站点位置处的a
max
是一个增强不等式显著性的参数，本实例取0.15。
[0072]
同理，如果，
[0073][0074]
那么，的位置就是产品s推荐的最小雨量虚拟站点的位置。反之，产品s未推荐最小雨量虚拟站点。
[0075]
式中，a
min
是一个增强不等式显著性的参数，本实例取0.05。
[0076]
所有降水产品均完成推荐虚拟站点的步骤后，若超过半数产品推荐的虚拟站点位置相同或者相邻，则该位置为最终确定的虚拟站点位置，反之，该子区域不需补充虚拟站点。如图2所示，三种产品在图示子区域的雨量分布均显示白点处为的位置，且该
位置雨量大于实际雨量站所在的黑点的雨量，即因此图中白点被多源降水产品推荐为最大雨量虚拟站点。
[0077]
第三步，利用随机森林(rf)作为确定虚拟站降雨量的非空间类融合方法估计虚拟站点处的降雨量。
[0078]
首先，在实际站点位置训练随机森林模型f(
·
)。
[0079][0080]
式中，f(
·
)是所训练的随机森林模型，pg是实际站点g的实测雨量，是降雨产品si在g处的雨量估计值，本实例共8种产品，ag是与g相关的协变量集合，本实例选取距离站点g最近的两个实际站点的雨量以及这两个最近的实际雨站点到目标格点的距离，ε是随机森林模型的误差。
[0081]
其次，基于训练后的f(
·
)估计虚拟站点v处的降雨量
[0082][0083]
式中，降雨产品si在v处的雨量估计值，本实例共8种产品，av是与v相关的协变量集合，本实例选取距离站点v最近的两个实际站点的雨量以及这两个最近的实际雨站点到目标格点的距离。
[0084]
第四步，判断对象区域的所有格点距最近雨量站点的距离是否小于第一步中选定的雨量站点可控距离，其中雨量站点是实际站点和虚拟站点的总集合g＝[g1,g2,
…
,g
10
,v1,v2,
…
]。若是，则转到第五步；若否，则以站点总集合作为更新后的实际站点集合并转到第一步，如图3所示。
[0085]
第五步，基于反距离权重(idw)法，利用实际站点和虚拟站点进行空间插值，得到降雨场。
[0086]
最后，选用皮尔逊相关系数和均方根误差评价估计的降雨场的精度，其表达式如下：
[0087][0088][0089]
其中，
[0090][0091][0092]
式中，t为参与评价的时间步长数，n为参与评价的格点数量，指的是格点x处的降雨量。
[0093]
以上为10个实际站点的少资料地区的实施步骤。在研究区域内选择2、3、4、6、14、22、30个缺测率小于2％的格点方案假定分别为少资料地区的实际雨量站点，其他格点值用于评价降雨场的估计精度。基于上述选取的7种雨量站点方案分别驱动本发明。
[0094]
站点密度不同时，反距离权重法(idw)、随机森林法(rf)以及本发明估计雨量场的精度(and)如表2所示。可以看出，在站点极为稀疏的情况下，随机森林法(rf)对雨量场的估计比反距离权重法(idw)好；当站点密度超过约300站/106km2时，反距离权重法(idw)对雨量场的估计比随机森林法(rf)好。但无论是哪种站点密度方案，本发明对雨量场的估计都是精度最高的。
[0095]
表2不同方法在不同站点密度雨量场估计综合表现
[0096][0097]
以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。