一种逐日降水数据综合插值方法

本发明涉及水文技术领域，特别涉及一种逐日降水数据综合插值方法。

背景技术：

降水是水文循环的重要环节，降水数据是水文模型的基础驱动数据，其时空分辨率和数据质量对水文数据的后续分析与利用有关键性的影响。在我国，由于地形、人力、物力等因素的限制，能够提供最准确降水数据的地面雨量计观测面临着雨量站点分布东多西少、平原多山地少，区域代表性不足；用卫星遥感进行降水数据反演受到传感器精度和反演算法的限制，其相关降水数据产品的精度还有待改进。气象观澜部门及水文气象研究机构也制作了很多降水数据，如国家气象科学数据中心提供的中国地面降水日值0.5°×0.5°格点数据集(v2.0)，国家青藏高原科学数据中心提供的中国区域地面气象要素驱动数据集(1979-2018)，数据空间分辨率较高(0.1°×0.1°)，但这些数据在精度及数据空间分辨率方面不能满足高精度水文应用需求。因此，如何获取兼顾高质量和高时空分辨率的降水产品受到研究者的广泛关注。

近年来对降水数据插值方法的研究不断涌现，例如，shubhamtiwari等(2019)将idw方法进行改进后用于澳大利亚murray-darling流域年降水插值，结果较好。olguaydin(2018)以高程作为辅助变量，用外部漂移克里金方法对爱琴海的年降水量进行了插值，ha-gyujeong等(2020)将prism与idw结果进行平均，改善山区降水高估问题。但目前国内外对高空间分辨率的逐日降水格点数据的插值计算方法的研究还不成熟，因此需要一种能够结合多种插值方法、能设定高空间分辨率、提高逐日降水数据插值精度的基于插值方法组合与误差校正的逐日降水数据综合插值方法。

针对此现象，本申请提供一种逐日降水数据综合插值方法，克服现有单一插值方法的精度不足，提供基于多种插值方法组合与误差校正的逐日降水格点数据综合插值方法，其充分利丰富的地面站点观测数据，考虑高程对降水空间分布的影响，即保证无观测站格点的降水插值精度，又保证在有观测站格点的插值结果最大限度接近观测值，可以设定高空间分辨率待插值格点、实现逐日程序化计算，能够解决高空间分辨率逐日降水插值精度低的问题。

技术实现要素：

本发明的目的在于提供一种逐日降水数据综合插值方法，克服现有单一插值方法的精度不足，提供基于多种插值方法组合与误差校正的逐日降水格点数据综合插值方法，其充分利丰富的地面站点观测数据，考虑高程对降水空间分布的影响，即保证无观测站格点的降水插值精度，又保证在有观测站格点的插值结果最大限度接近观测值，可以设定高空间分辨率待插值格点、实现逐日程序化计算，能够解决高空间分辨率逐日降水插值精度低的问题。

本发明提供了一种逐日降水数据综合插值方法，包括以下步骤：

s1：利用高程作为辅助变量的外部漂移克里金插值方法ked及薄板样条插值方法tps，对研究区域目标内全部n个雨量站点的逐日降水观测数据p，进行逐日插值计算，得到逐日降水栅格数据p_ked和p_tps；

s2：将n个实测站点分为6组，每组包含n/6个均匀分布于全区域的站点，每次采用其中5组站点的逐日降水数据，用ked、tps方法分别进行逐日降水插值计算，并计算未参与插值计算的1组站点处的插值数据均方根误差rmse；通过6次循环，通过上述交叉检验得到全部站点采用ked、tps插值方法的rmse；

s3：采用反距离加权idw方法，对由ked与tps方法求出的站点rmse插值到全部格点，得到ked与tps方法的逐格点均方根误差r_ked和r_tps，进而根据二者数值的倒数，计算出ked和tps在整个区域每个格点上的加权平均权重w_ked、w_tps；

s4：利用权重值w_ked、w_tps，对p_ked与p_tps进行加权平均计算，得到逐日降水的加权平均插值格点数据ci；

s5：计算加权平均插值格点数据ci与实测站点数据p在各站点上的逐日误差，并对站点误差采用idw插值，得到逐格点的误差数据erro；

s6：利用逐格点的误差数据erro，对加权平均插值格点数据ci进行修正，得到最终的逐日降水综合插值格点数据cci。

进一步地，所述步骤s1包括以下步骤：

s1.1：在数据服务平台下载空间分辨率为30m×30m的数字高程数据astgtm2_dem，利用arcgis中的aggregate工具聚合为待插值格点空间分辨率；

s1.2：根据研究区域内生成的所需空间分辨率的待插值点经纬度信息，并利用arcgis中的extractvaluetopoints工具提取待插值点及n个站点对应高程；

s1.3：将n个站点的站点号、经纬度、高程、年、月、日、日降水量一一对应，得到全站点降水信息表；

s1.4：对逐日降水观测数据p进行外部漂移克里金插值方法ked和薄板样条插值方法tps逐日插值计算，得到逐日降水格点数据p_ked、p_tps。

进一步地，所述步骤s2包括以下步骤：

s2.1：将n个站点数据按纬度由小到大排序，按1-6循环标号，得到6组p1～p6含n/6个空间分布均匀的站点；

s2.2：参照步骤s1.4的方法，分6次，每次利用p1～p6中的5组站点逐日降水数据，分别采用ked和tps方法进行插值，得到6组ked插值结果ked1-ked6和6组tps值结果tps1～tps6；

s2.3：对ked1～ked6、tps1～tps6分别计算其在每次未参与插值计算的那1组站点处的均方根误差，综合后得到在全部6组共n个站点上的ked、tps插值均方根误差a_ked、a_tps。

进一步地，所述步骤s3需要对站点均方根误差进行插值，并计算加权平均所占权重，具体包括以下步骤：

用反距离加权idw方法对a_ked、a_tps进行空间插值，得到ked和tps在所有格点上的均方根误差r_ked，r_tps；

s3.2：根据r_ked，r_tps的数值倒数，逐网格计算各自权重，计算公式如下：

其中，i表示行数，j表示列数。

进一步地，所述步骤s4利用权重值w_ked与w_tps，对p_ked与p_tps进行加权平均计算，得到逐日降水的加权平均插值格点数据ci，格点处加权平均的具体计算公式如下：

其中，d表示日期，i表示行数，j表示列数。

进一步地，所述步骤s5计算逐日加权平均插值结果ci与实测降水序列p的逐日误差，具体包括以下步骤：

s5.1：根据站点经纬度提取ci在站点id处的逐日插值降水时间序列

s5.2：截取p中站点id的逐日实测降水时间序列

s5.3：计算对应站点的逐日插值误差，公式如下：

其中，id表示站点号；

s5.4：利用idw方法对逐日站点误差erro进行空间插值，得到逐格点的逐日误差空间分布erro。

进一步地，所述步骤s6对逐日加权平均插值结果ci进行插值误差校正，第d日降水数据综合插值的具体计算公式如下：

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

本发明提供了一种逐日降水数据综合插值方法，除降水数据外仅需要流域高程、雨量站站点信息，获取简单。计算方法在r语言中都有现成程序包，调用方便，利用循环实现逐日降水自动插值，采用了科学的加权权重计算方法。雨量站点较多时可实现高空间分辨率降水插值。多插值方法组合运用可减少单一插值方法误差，误差校正可实现插值结果误差最小。

附图说明

图1为本发明实施例提供的一种逐日降水数据综合插值方法的原理流程示意图；

图2为以福建省为例，通过本发明获取的外部漂移克里金和薄板样条方法加权权重组合构成图；

图3为以福建省为例，通过本发明获取的逐日降水数据产品与其他方法所得产品在未使用的40个站点处的对比结果图；

图4为以福建省为例，通过本发明获取的逐日降水数据产品与其他方法所得产品在插值过程中使用过的20个站点处的对比结果图。

具体实施方式

下面结合本发明中的附图，对本发明实施例的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

参照图1-4，本发明提供了一种逐日降水数据综合插值方法，包括以下步骤：

s1：利用高程作为辅助变量的外部漂移克里金插值方法ked及薄板样条插值方法tps，对研究区域目标内全部n个雨量站点的逐日降水观测数据p，进行逐日插值计算，得到逐日降水栅格数据p_ked和p_tps；

s2：将n个实测站点分为6组，每组包含n/6个均匀分布于全区域的站点，每次采用其中5组站点的逐日降水数据，用ked、tps方法分别进行逐日降水插值计算，并计算未参与插值计算的1组站点处的插值数据均方根误差rmse，通过6次循环，通过上述交叉检验得到全部站点采用ked、tps插值方法的逐日插值均方根误差；

s3：采用反距离加权idw方法，对由ked与tps方法求出的站点rmse插值到全部格点，得到ked与tps方法的逐格点均方根误差r_ked和r_tps，进而根据二者数值的倒数，计算出ked和tps在整个区域每个格点上的加权平均权重w_ked、w_tps；

s4：利用权重值w_ked、w_tps，对p_ked与p_tps进行加权平均计算，得到逐日降水的加权平均插值格点数据ci；

s5：计算加权平均插值格点数据ci与实测站点数据p在各站点上的逐日误差，并对站点误差采用idw插值，得到逐格点的误差数据erro；

s6：利用逐格点的误差数据erro，对加权平均插值格点数据ci进行修正，得到最终的逐日降水综合插值格点数据cci。

本发明利用福建省及其周边共480个雨量站点1979-2018年逐日降水数据，得到福建省0.05°×0.05°逐日降水格点数据产品，并进行检验，精度比现有逐日降水格点数据产品高。

实施例1

所述步骤s1包括以下步骤：

s1.1：在数据服务平台下载空间分辨率为30m×30m的数字高程数据astgtm2_dem，利用arcgis中的aggregate工具聚合为待插值格点空间分辨率(例如0.05°×0.05°)；

s1.2：根据研究区域内生成0.05°×0.05°的待插值点经纬度信息，并利用arcgis中的extractvaluetopoints工具提取待插值点及440个站点对应高程；

s1.3：将n个站点的站点号、经纬度、高程、年、月、日、日降水量一一对应，得到全站点降水信息表；

s1.4：对逐日降水观测数据p进行外部漂移克里金插值方法ked和薄板样条插值方法tps逐日插值计算，得到逐日降水格点数据p_ked、p_tps。

实施例2

所述步骤s2包括以下步骤：

s2.1：将480个站点数据按纬度由小到大排序，按1-6循环标号，分别剔除1～6号站点，得到6组含400个空间分布均匀的站点降水资料数据集p1～p6，保证每个站点都被且只被剔除一次；

s2.2：参照步骤s1.4的方法对p1～p6进行插值，分6次，每次利用p1～p6中的5组站点逐日降水数据，分别采用ked和tps方法进行插值，得到6组ked插值结果ked1-ked6和6组tps值结果tps1～tps6；

s2.3：对ked1～ked6、tps1～tps6分别计算其在每次未参与插值计算的那1组站点处的均方根误差，综合后得到在全部6组共n个站点上的ked、tps插值均方根误差a_ked、a_tps。

实施例3

参照图2，所述步骤s3需要对站点均方根误差进行插值，并计算加权平均所占权重，具体包括以下步骤：

s3.1：用反距离加权idw方法对a_ked、a_tps进行空间插值，得到ked和tps在所有格点上的均方根误差r_ked，r_tps；

s3.2：根据r_ked，r_tps的数值倒数，逐网格计算各自权重，计算公式如下：

其中，i表示行数，j表示列数。

实施例4

所述步骤s4利用权重值w_ked与w_tps，对p_ked与p_tps进行加权平均计算，得到逐日降水的加权平均插值格点数据ci，格点处加权平均的具体计算公式如下：

其中，d表示日期，i表示行数，j表示列数。

实施例5

所述步骤s5计算逐日加权平均插值结果ci与实测降水序列p的逐日误差，具体包括以下步骤：

s5.1：根据站点经纬度提取ci在站点id处的逐日插值降水时间序列

s5.2：截取p中站点id的逐日实测降水时间序列

s5.3：计算对应站点的逐日插值误差，公式如下：

其中，id表示站点号；

s5.4：利用idw方法对逐日站点误差erro进行空间插值，得到逐格点的逐日误差空间分布erro。

实施例6

所述步骤s6利用逐格点的误差数据，对逐日加权平均插值结果ci进行插值误差校正，第d日降水数据综合插值的具体计算公式如下：

在示例研究区内，利用基于插值方法组合与误差校正的逐日降水数据综合插值方法插值得到的逐日降水数据cci比常用的单一插值方法ked、tps、idw以及现有数据集中国气象数据网提供的中国地面降水日值0.5°×0.5°格点数据集(v2.0)cma具有更高的与站点实测数据相关系数、更小的均方根误差以及更小的平均偏差，参照图3为插值计算过程中未使用的40个站点处的不同方法(或产品)精度对比结果，图4为插值过程中使用过的20个站点处的不同方法(或产品)对比结果。

以上公开的仅为本发明的几个具体实施例，但是，本发明实施例并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。