具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下文中将结合附图对本发明的实施例进行详细说明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互任意组合。

图1为本发明实施例复杂地形区域的降水数据空间插值方法的流程图，如图1所示，包括：将包括多个站点的待处理区域划分为两个以上网格，以其中一个网格作为目标网格进行如下处理；

步骤101，将距离目标网格预设范围内的站点确定为辅助站点；

在一种示例性实例中，以长江上游流域作为待处理区域为例，共使用396个站点，其中数字高程模型(Digital Elevation Model)数据，使用了原始精度90m的SRTM(ShuttleRadar Topography Mission)数据。对于待处理区域的每一个网格，可选择距离目标网格最近的8个站点为辅助站点。

步骤102，根据辅助站点与目标网格的距离分别获得每个辅助站点的距离权重；

在一种示例性实例中，步骤102可以通过公式(1)获得：

公式(1)中，x为目标网格(m,n)与站点i的距离；根据长江上游流域站点的分布情况，本申请实施例中，假设x₀取值为400km；s为校正系数，按照实际应用场景，比如可以取值为4。

步骤103，根据辅助站点中的站点与站点间以目标网格为中心所形成的角度，计算距离方向权重；

在一种示例性实例中，步骤103可以通过公式(2)获得：

w_(m,n),i＝w_0(m,n),i(1+a_(m,n),i) (2)

公式(2)中的调整系数可以通过公式(3)获得：

公式(3)中，cosθ_m,n(i,j)表示第i、j个站点以目标网格(m,n)为中心所形成的角度，如附图3所示；nog表示所用站点个数；w_0(m,n),j表示第j个站点的距离权重。

步骤104，根据辅助站点的多年平均月降水数据及海拔高程数据，拟合得到降水高程梯度；

在一种示例性实例中，根据多年平均月降水数据与海拔高程数据的相关程度，选择距离目标网格最近的5个站点，拟合得到降水高程梯度。

步骤105，根据拟合得到的降水高程梯度，计算辅助站点到目标网格路径上的第一平均降水高程梯度；

在一种示例性实例中，步骤105可以通过公式(4)获得：

公式(4)中，ele(u_m,n)和ele(u_g)分别表示目标网格u_m,n和站点g的海拔高程数据，Nop表示路径上的网格数目。附图3标示了站点s2到目标网格(m,n)的路径。slope_j为路径上第j个网格的降水高程梯度。

步骤106，利用距离方向权重和第一平均降水高程梯度对目标网格的多年平均月降水数据进行插值，得到目标网格的多年平均月降水数据的插值。

在一种示例性实例中，步骤106可以通过公式(5)获得：

公式(5)中，w_g(u_m,n)表示插值所用的第g个站点到目标网格u_m,n的距离方向权重；Nog表示插值所用的站点数目；S(u_g,u_m,n)表示第g个站点到网格u_m,n平均降水高程梯度，见计算公式(4)。

在一种示例性实例中，目标网格的多年平均7月份降水数据的插值参见附图6a。

本发明实施例根据降水-高程梯度关系对距离方向权重进行修正，可以估计降水数据随地形的变化的情况，避免分布密集站点权重影响太大的问题，提高了降水空间插值的精度。

在一种示例性实例中，本发明实施例方法还包括：

根据所述辅助站点到所述目标网格的距离计算多年平均月降水数据插值的不确定性；

根据所述多年平均月降水数据插值的不确定性，计算多年平均月降水数据插值的平均不确定性。

在一种示例性实例中，所述计算多年平均月降水数据插值的平均不确定性，包括：

根据所述降水高程梯度，计算辅助站点中的站点到站点的第二平均降水高程梯度。

所述第二平均降水高程梯度，可以通过公式(6)获得：

公式(6)中，Nop表示路径上的网格数目；ele(u_i+h)和ele(u_i)表示距离为h的两个站点的海拔高程；slope_j为路径上第j个网格的降水高程梯度。

多年平均月降水数据插值的平均不确定性，可以通过公式(7)获得：

公式(7)中，n(h)为相互距离为h的站点对的个数；P(u_i+h)和P(u_i)为多年平均月降水数据。

对于地形复杂的区域，可以进一步考虑站点的高程变化对降水空间变异性的影响。在这种情况下，利用所述第二平均降水高程梯度计算多年平均月降水数据插值的平均不确定性，可以通过公式(8)获得：

公式(8)中，ele(u_i+h)及ele(u_i)均表示海拔高程；S(u_i,u_i+h)表示第二平均降水高程梯度，见计算公式(6)。

在一种示例性实例中，所述根据辅助站点到目标网格的距离计算多年平均月降水数据插值的不确定性，包括：

采用预先设置的数学模型拟合半变异函数与相应的距离的关系，绘制多年平均各月份的拟合曲线(如附图5中曲线所示)。

在一种示例性实例中，利用长江上游396个站点的多年平均月降水数据，分别在1至12月份计算半变异函数值。根据所用站点的分布情况，以100km为单位划分区间，共分为22个区间，每个区间最终的值为区间内所有的算术平均值，便可求得共22个值及相应的距离h，绘制每对离散点(如附图5中离散点所示)。通常使用带有块金效应的球形模型或指数模型等理论变异函数模型分析降水量的空间变异。

在一种示例性实例中，所述数学模型为球形模型或指数模型。

在一种示例性实例中，球形模型和指数模型拟合半变异函数与距离的关系，分别可以通过公式(9)和(10)获得：

公式(9)和(10)中，C₀表示块金值，反映因为采样密度不足从而无法刻画的变异特征，或者来源于观测的误差；a表示变程值，代表空间相关尺度，一般如果采样点间距小于该距离，则说明观测数据的相关性较强；C₀+C表示基台值，C与C₀的相对大小可以体现变量的空间变异性是由空间相关过程引起或是随机效应引起。

分别用带有块金效应的球形模型或指数模型拟合半变异函数值及相应的距离h的关系，在长江上游地区，指数模型拟合得更好，因此优选指数模型拟合。通过目测样本变异函数值的分布来估计指数模型中的初始参数值，并通过试算的方式进一步确定更为合理的指数模型参数初始值。按照实际应用场景，块金值C0可以取值0.5，变程值a可以取值1900，基台值C可以取值400。

在一种示例性实例中，本发明实施例方法还包括：

对所述多年平均月降水数据插值的平均不确定性进行量化。

在一种示例性实例中，所述对多年平均月降水数据插值的平均不确定性进行量化，包括：

计算目标网格的多年平均月降水数据插值的方差。

所述方差可以通过公式(11)获得：

多年平均月降水数据插值结果的不确定性量化结果参见附图6b。

在一种示例性实例中，本发明实施例方法还包括：

计算目标网格的日降水观测值与日降水气候值的比值；

在一种示例性实例中，在日尺度上，降水数据的时空变异性更大，不同年份同一天的降水数据可能有明显不同的降水-高程梯度。所以，用多年平均各月份(1-12)降水数据插值的结果作为研究区域各月份降水的气候值，同时假设各月份内日尺度的气候值等于各月份降水气候值除以相应天数，进一步计算目标网格的日降水观测值与日降水气候值的比值。

所述比值可以通过公式(12)获得：

公式(12)中，x_cli(u_i)表示站点u_i处的日降水气候值；x(u_i)表示日降水观测值。

采用距离方向权重插值方法对所述比值进行空间插值，得到比值的空间插值。

在一种示例性实例中，鉴于日尺度降水气候值的空间分布已包含了地形变化对降水的作用，所以将比值进行空间插值时不再考虑高程。

所述比值的空间插值可以通过公式(13)获得：

根据目标网格的日降水气候值及所述比值的空间插值，计算目标网格的日降水数据的插值。

所述目标网格的日降水数据的插值，可以通过公式(14)获得：

x_ADW(u_m,n)＝x_cli(u_m,n)×ratio_ADW(u_m,n) (14)

在一种示例性实例中，目标网格2010年10月1日降水数据的插值参见附图8a。

在一种示例性实例中，本发明实施例方法还包括：

利用所述目标网格的日降水观测值与日降水气候值的比值，计算目标网格日降水数据插值的平均不确定性。

在一种示例性实例中，半变异函数可以通过公式(15)获得：

公式(15)中，ratio(u_i+h)和ratio(u_i)分别表示相互距离为h的两个站点的ratio值，日尺度的经验表达式的拟合方式与月尺度一致。

在日尺度上，根据日尺度上降水插值的方法以及误差量化方法，在获得每一天的半变异函数后，由于研究区域有些站点在某些时段有缺测情况，因此对某一月份的所有天数的半变异函数值取平均值的方法得到多年平均该月份每一天的半变异函数值。同样分别用带有块金效应的球形模型或指数模型拟合半变异函数值及相应的距离h的关系，指数模型拟合更佳，因此优选指数模型拟合。通过目测样本变异函数值的分布来估计指数模型中的初始参数值，并通过试算的方式进一步确定更为合理的指数模型初始参数值。按照实际应用场景，块金值C0可以取值0，变程值a可以取值20，基台值C可以取值400。

对目标网格日降水数据插值的平均不确定性进行量化，可以通过公式(16)获得：，

公式(16)中：为当站点g到目标网格u_m,n的距离为h时，Ratio的插值结果的半变异函数经验拟合值。

在一种示例性实例中，目标网格日降水数据插值的平均不确定性量化结果参见附图8b。

本发明实施例还提供一种计算机存储介质，所述计算机存储介质中存储有计算机程序，所述计算机程序被处理器执行时实现上述降水数据空间插值方法。

本发明实施例对比了考虑高程及不考虑高程两种情况下插值结果的误差。具体做法为使用留一法交叉检验方法，针对长江上游396个站点所在的网格，每个网格处的降水分别在不考虑高程要素和考虑高程要素的情况下由距离其最近的其他8个站点通过ADW插值方法得到。并基于站点实测数据，对两种结果进行了检验，检验指标为月均绝对误差(MAAE)。图9a、9b分别展示了两种情况下的MAAE交叉验证结果。总体上，在研究区域的大部分区域，两种插值结果MAAE均较小，具体表现为在所用的396个站点所在网格中，当未考虑高程时，MAAE小于5mm的站点占145个，MAAE小于15mm的站点数共为351个；当考虑高程时，MAAE小于5mm的站点占159个，MAAE小于15mm的站点数为354个。在考虑高程的情况下，有60.6％的站点所在网格的插值效果优于未考虑高程要素的情况。改善情况如图10所示，可直观发现MAAE得到改善的站点主要分布在山区。说明在地形复杂的山区，考虑高程要素后的插值方法更为合理。表1统计了各高程区间考虑高程要素后的MAAE得到改善的站点数所占比重。

表1不同高程区间内考虑高程要素插值后误差的改善情况

在低海拔地区改善的比重较小，但随着高程增加，改善的比重有增大趋势。在海拔2500～3500m的高程区间内，改善的比重最大。综合以上分析，考虑高程要素后的ADW插值方法更为合理，插值结果的精度也更高。

本领域普通技术人员可以理解，上文中所公开方法中的全部或某些步骤、系统、装置中的功能模块/单元可以被实施为软件、固件、硬件及其适当的组合。在硬件实施方式中，在以上描述中提及的功能模块/单元之间的划分不一定对应于物理组件的划分；例如，一个物理组件可以具有多个功能，或者一个功能或步骤可以由若干物理组件合作执行。某些组件或所有组件可以被实施为由处理器，如数字信号处理器或微处理器执行的软件，或者被实施为硬件，或者被实施为集成电路，如专用集成电路。这样的软件可以分布在计算机可读介质上，计算机可读介质可以包括计算机存储介质(或非暂时性介质)和通信介质(或暂时性介质)。如本领域普通技术人员公知的，术语计算机存储介质包括在用于存储信息(诸如计算机可读指令、数据结构、程序模块或其他数据)的任何方法或技术中实施的易失性和非易失性、可移除和不可移除介质。计算机存储介质包括但不限于RAM、ROM、EEPROM、闪存或其他存储器技术、CD-ROM、数字多功能盘(DVD)或其他光盘存储、磁盒、磁带、磁盘存储或其他磁存储装置、或者可以用于存储期望的信息并且可以被计算机访问的任何其他的介质。此外，本领域普通技术人员公知的是，通信介质通常包含计算机可读指令、数据结构、程序模块或者诸如载波或其他传输机制之类的调制数据信号中的其他数据，并且可包括任何信息递送介质。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。