水文气象要素空间插值最优分辨率的确定方法

摘要

本发明公开了一种水文气象要素空间插值最优分辨率的确定方法，包括以下步骤：在同一地区观测站点数目相同情况下，采用Kriging插值法对水文气象要素观测数据在不同分辨率下进行空间插值展布，得到空间插值结果图像；计算各分辨率下空间插值结果图像的空间插值信息量、空间插值噪声量及空间插值信噪比；对各分辨率下的空间插值信噪比进行排序，空间插值信噪比最大的分辨率为空间差值最优分辨率。本发明通过定义空间插值信噪比的概念，作为一个衡量空间插值结果信息质量的指标，使空间插值质量具体量化，并通过计算各个分辨率下的空间插值信噪比，找到最优分辨率；在现有样本数据资料下，通过选择最优分辨率，使获得的空间插值结果图像包含的有效数据最大化，更具参考使用价值。

说明书

技术领域

本发明涉及水文、水资源研究领域，具体涉及水文气象要素空间插值最优分辨率的确定方法。

背景技术

空间插值是在观测站点基础上获取水文气象要素的空间展布信息的重要方式，这种以专题图像为载体的信息方式，能为我们提供观测站点无法提供的连续空间信息。这种通过空间插值得到的图像所包含的信息精度由原始观测站点的数目和图像空间分辨率大小两方面决定。

通常，人们在关注插值结果的时候，着眼点主要在原始观测站点数据的空间采样特性，即站点样本容量、站点位置、数据空间密度、数据空间分布和空间变化等对插值精度的影响。实际上，除了原始观测站点特性，空间插值中所选的分辨率不同带来的不确定性对插值精度的影响同样值得关注。Gottschalk等人(Gottschalk L.,Leblois,E andJ.Correlation and covariance of runoff revisited.Journal of Hydrology,398(2011)76-90,doi:10.1016/j.jhydrol.2010.12.011)在研究地理对象相关关系时，就开始注意到栅格尺度变化(空间分辨率变化)给插值结果可靠性带来了影响，在其它条件不变的情况下，增加分辨率，减少栅格尺寸，虽然图像中可用信息量变大，但是像元间的协方差增加，相应的根据Kriging插值算法，单个单元的插值误差就越大。

可见，对于空间插值得到的图像所包含的信息精度而言，并非分辨率越大，精度越高，在给定观测站点数目的情况下，如何确定适当分辨率使水文气象要素空间插值最优，为我们提供精度更高的信息量，变得至关重要。

发明内容

本发明所要解决的技术问题是在给定观测站点数目的情况下，如何确定适当分辨率使水文气象要素空间插值最优的问题。

为了解决上述技术问题，本发明所采用的技术方案是提供一种水文气象要素空间插值最优分辨率的确定方法，水文气象要素空间插值最优分辨率的确定方法，包括以下步骤：

在同一地区观测站点数目相同情况下，采用Kriging插值法对水文气象要素观测数据在不同分辨率下进行空间插值展布，得到空间插值结果图像；

计算各分辨率下空间插值结果图像的空间插值信息量、空间插值噪声量及空间插值信噪比；

对各分辨率下的空间插值信噪比进行排序，空间插值信噪比最大的分辨率为空间差值最优分辨率。

在上述方法中，所述Kriging插值法公式为：

$>Z_V^{*}=\underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i{Z}_i;$>

其中，为目标单元V的估计值；Z_i为第i个观测样本；n是观测样本Z_i的数目；λ_i是每一个观测样本Z_i的权重。

在上述方法中，每一个观测样本的权重λ_i通过采用拉格朗日算法求解λ_i的方程组求得，方程组为：

$>\left(\begin{array}{cc}\underset{j=1}{\overset{n}{\Sigma }}{\lambda }_{j}C(x_i,x_j)-\mu =\bar{C}(x_i,V)& i=1,2,\mathrm{...},n\\ \underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i=1& \end{array}\right)$>

其中，V代表目标单元，x观测站点，C(x_i,x_j)代表观测站点之间的协方差函数；代表观测站点与目标单元之间的协方差函数。

在上述方法中，

所述空间插值信息量值为：采用Kriging插值法得到的空间插值结果图像上所有栅格作为样本计算得到的方差：

$>{\mathrm{Signal}}_{map}=D\left(Z^{*}\right)=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}{(Z_s^{*}-\overline{Z^{*}})}^2$>

其中，m为目标单元的数目，为第s个目标单元的估计值，为所有单元估计值的均值；

所述空间插值噪声量值为：采用Kriging插值法得到的空间插值结果图像上所有栅格估计方差的平均值：

$>{\mathrm{Noise}}_{map}=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}{\sigma }_s^2=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}\left(\bar{C}\right(V,V)+\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{\lambda }_i{\lambda }_{j}C(x_i,x_j)-2\underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i\bar{C}(V_s,x_i\left)\right)$>

其中，为Kriging插值法计算得到的空间插值结果图像上每个目标单元的插值方差，为目标单元与目标单元之间的协方差函数；

所述空间插值信噪比SNR_map为：当前分辨率下，空间插值信息量Signal_map值与空间插值噪声量Noise_map值之比，公示表达为：

$>{\mathrm{SNR}}_{map}=\frac{{\mathrm{Signal}}_{map}}{{\mathrm{Noise}}_{map}}.$>

在上述方法中，在对各个分辨率下的空间插值信噪比进行排序后，采用绘制空间插值信噪比与分辨率的曲线图方式得到最优分辨率，曲线图中最高点对应的分辨率为最优分辨率。

本发明通过定义空间插值信噪比的概念，作为一个衡量空间插值结果信息质量的指标，使空间插值质量具体量化，并通过计算各个分辨率下的空间插值信噪比，找到最优分辨率；在现有样本数据资料下，通过选择最优分辨率，使获得的空间插值结果图像包含的有效数据最大化，更具参考使用价值。

附图说明

图1为本发明提供的水文气象要素空间插值最优分辨率的确定方法流程图；

图2为本发明实施例中空间插值信噪比与分辨率的曲线图；

图3为本发明实施例一中降雨量的不同分辨率空间插值结果图；

图4为本发明实施例一中降雨量空间插值信息量与分辨率的曲线图；

图5为本发明实施例一中降雨量空间插值噪声量与分辨率的曲线图；

图6为本发明实施例一中降雨量空间插值信噪比与分辨率的曲线图；

图7为本发明实施例二中实际蒸发量的不同分辨率空间插值结果图；

图8为本发明实施例二中实际蒸发量空间插值信息量与分辨率的曲线图；

图9为本发明实施例二中实际蒸发量空间插值噪声量与分辨率的曲线图；

图10为本发明实施例二中实际蒸发量空间插值信噪比与分辨率的曲线图；

图11为本发明实施例三中径流深的不同分辨率空间插值结果图；

图12为本发明实施例三中径流深空间插值信息量与分辨率的曲线图；

图13为本发明实施例三中径流深空间插值噪声量与分辨率的曲线图；

图14为本发明实施例三中径流深空间插值信噪比与分辨率的曲线图。

具体实施方式

本发明借鉴了信号学和图像处理技术中的“信噪比”概念，并在此基础上提出空间插值信噪比的概念，揭示空间分辨率变化给插值结果带来的规律性特征，从而得到在给定观测站点数目的情况下，确定使水文气象要素空间插值最优的分辨率。其中，

信噪比SNR(Signal to Noise Ratio)，又称为“讯噪比”；在信号学中，信噪比SNR为有用信号功率(Power of Signal)P_signal与噪声功率(Power>noise的比，即为用信号幅度A_signal(Amplitude>noise(Amplitude>

$>SNR=\frac{P_{signal}}{P_{noise}}={\left(\frac{A_{signal}}{A_{noise}}\right)}^2;$>

信噪比SNR的单位一般使用分贝，其值为十倍对数信号与噪声功率比：

$>SNR\left(dB\right)=10{\log }_{10}\left(\frac{P_{signal}}{P_{noise}}\right)=20{\log }_{10}\left(\frac{A_{signal}}{A_{noise}}\right);$>

由信号学信噪比定义可知，信噪比越大，说明混在信号里的噪声越小，信号质量越高，否则相反。

下面结合说明书附图和具体实施例对本发明做出详细的说明。

如图1所示，本发明提供的水文气象要素空间插值最优分辨率的确定方法，包括以下步骤：

步骤S10、在同一地区观测站点数目相同情况下，采用Kriging插值法对水文气象要素观测数据在不同分辨率下进行空间插值展布，得到空间插值结果图像。

其中，Kriging插值法包括以Kriging插值法为内核的所有插值法，具体可以根据文气象要素数据结构特点选择最为适当的插值法，Kriging插值法公式为：

$>Z_V^{*}=\underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i{Z}_i;$>

其中，为目标单元V的估计值；Z_i为第i个观测样本；n是观测样本Z_i的数目；λ_i是每一个观测样本Z_i的权重。

每一个观测样本的权重λ_i通过采用拉格朗日算法求解λ_i的方程组求得，方程组为：

$>\left(\begin{array}{cc}\underset{j=1}{\overset{n}{\Sigma }}{\lambda }_{j}C(x_i,x_j)-\mu =\bar{C}(x_i,V)& i=1,2,\mathrm{...},n\\ \underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i=1& \end{array}\right)$>

其中，V代表目标单元，x观测站点，C(x_i,x_j)代表观测站点之间的协方差函数；代表观测站点与目标单元之间的协方差函数。

步骤S20、计算各分辨率下空间插值结果图像的空间插值信息量和空间插值噪声量及空间插值信噪比。

根据信号学信噪比的概念，定义空间插值信息量为插值结果图像所有栅格值的总体方差，反映的是图像上所有栅格值的差异程度，空间插值噪声量为空间插值结果图像所有栅格的插值估计方差的平均值，反映的是图像上所有栅格插值误差的平均值，所以：

空间插值信息量Signal_map值为：采用Kriging插值法得到的空间插值结果图像上所有栅格作为样本计算得到的方差：

$>{\mathrm{Signal}}_{map}=D\left(Z^{*}\right)=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}{(Z_s^{*}-\overline{Z^{*}})}^2$>

其中，m为目标单元的数目，为第s个目标单元的估计值；为所有单元估计值的均值。

空间插值结果图像噪声量Noise_map值为：采用Kriging插值法得到的空间插值结果图像上所有栅格估计方差的平均值：

$>{\mathrm{Noise}}_{map}=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}{\sigma }_s^2=\frac{1}{m}\underset{s=1}{\overset{m}{\Sigma }}\left(\bar{C}\right(V,V)+\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{\lambda }_i{\lambda }_{j}C(x_i,x_j)-2\underset{i=1}{\overset{n}{\Sigma }}{\lambda }_i\bar{C}(V_s,x_i\left)\right)$>

其中，m为目标单元的数目，为Kriging插值法计算得到的空间插值结果图像上每个目标单元的插值方差，为目标单元与目标单元之间的协方差函数；λ_i为采用拉格朗日算法求解的到的第i个观测值的插值权重；代表单元与点之间的协方差函数。

空间插值信噪比SNR_map为：当前分辨率下，空间插值信息量Signal_map值与空间插值噪声量Noise_map值之比，即

$>{\mathrm{SNR}}_{map}=\frac{{\mathrm{Signal}}_{map}}{{\mathrm{Noise}}_{map}}.$>

由空间插值信噪比定义可知，其可作为一个衡量空间插值结果信息质量的指标。

步骤S30、对各分辨率下的空间插值信噪比进行排序，选取信噪比最大的分辨率作为空间差值最优分辨率；为了更为直观，也可在对各个分辨率下的空间插值信噪比进行排序后，绘制空间插值信噪比与分辨率的曲线图(如图2所示)，图中最高点对应的分辨率即为最优分辨率。

下面以水文循环三要素降雨量、实际蒸发量和径流深为具体实施例对本发明进行详细说明。

实施例一：降雨量。

(1)在同一研究区，选择相同观测站点数目的降雨数据，采用Kriging插值法对降雨量观测数据在不同分辨率下进行空间插值展布，其中，观测站点数目为145个，输出分辨率分别选取5km、10km、15km、20km、25km、30km、35km、40km、45km、50km，输出降雨量的不同分辨率空间插值结果图如图3所示，从图中可以明显看出随着栅格尺寸的减小，降雨空间插值结果的空间特征均有明显的变化。

(2)计算降雨量空间插值信噪比。

以上述输出的降雨量的不同分辨率插值结果图像为基础，利用Kriging插值法及信噪比公式，计算各个分辨率下降雨量的空间插值信息量、噪声量和信噪比，如表2所示。

表2：不同分辨率下输出的降雨量空间插值信息量、噪声量及信噪比。

分辨率信息量噪声量信噪比547306.2587705.3940991042527.673295.802648153724960356.1721622026846.8250085.3607872518251.3536634.982624301428933434.2743033510397.4225824.026886406489.19823542.756668153013.22622321.3500125084120000.4205

(3)对各个分辨率下的空间插值信噪比进行排序后，绘制空间插值信息量、噪声量及信噪比与分辨率的曲线图，分别如图4、5、6所示；从图4、5、6可知，当输出栅格大小依次选取5km、10km、15km、20km、25km、30km、35km、40km、45km、50km时，随着分辨率的降低，对于降雨量空差值结果来说，无论是信息量还是噪声量都在减少，而各个分辨率下的空间插值信噪比则出现了先增后减的趋势，且在15km分辨率下，降雨空间插值信噪比达到最大，降雨空间插值信噪比SNR_P_MAX＝6.1，即在此分辨率下，相邻栅格的差异程度是栅格插值平均误差的6.1倍。

实施例二：实际蒸发量。

(1)在同一研究区，选择相同观测站点数目的实际蒸发量数据，采用Kriging插值法对实际蒸发量观测数据在不同分辨率下进行空间插值展布，其中，观测站点数目为13个，输出分辨率分别选取20km、40km、60km、80km、100km、120km、140km、160km、180km、200km，输出实际蒸发量的不同分辨率空间插值结果图，如图7所示，从图中可以明显看出随着栅格尺寸的减小，实际蒸发空间插值结果的空间特征均有明显的变化。

(2)计算实际蒸发量空间插值信噪比。

以上述输出的实际蒸发量的不同分辨率空间插值结果为基础，利用Kriging插值法及信噪比公式，计算各个分辨率下实际蒸发量的空间插值信息量、噪声量和信噪比，如表3所示。

表3：不同输出分辨率下实际蒸发量空间插值信息量、噪声量及信噪比。

分辨率信息量噪声量信噪比201768836504.846027401494330254.9398356011756.7622115.317393808197.21313975.867726100492011174.4046551202714.9588833.0746981401844.9077002.635581160955.10435001.9102091804842001.42001211001.21

对各个分辨率下的空间插值信噪比进行排序后，绘制空间插值信息量、噪声量及信噪比与分辨率的曲线图，分别如图8、9、10所示，从图8、9、10中可见，当输出栅格大小依次选取20km、40km、60km、80km、100km、120km、140km、160km、180km、200km时，随着分辨率的降低，对于实际蒸发量来说，无论是信息量还是噪声量都在减少，而各个分辨率下的空间插值信噪比则出现了先增后减的趋势，且在80km分辨率下，实际蒸发量空间插值信噪比达到最大，实际蒸发量空间插值信噪比SNR_E_MAX＝5.9，即在此分辨率下，相邻栅格的差异程度是栅格插值平均误差的5.9倍。

实施例三：径流深。

(1)在同一研究区，选择相同观测站点数目的径流深数据，采用Kriging插值法对径流深观测数据在不同分辨率下进行空间插值展布，其中，观测站点数目为20个，输出分辨率分别选取5km、10km、15km、20km、25km、30km、35km、40km、45km、50km，输出径流深的不同分辨率空间插值结果图，如图11所示，从图中可以明显看出随着栅格尺寸的减小，径流深空间插值结果的空间特征均有明显的变化。

(2)计算实径流深空间插值信噪比。

以上述输出的径流深的不同分辨率空间插值结果为基础，借助Kriging插值法及信噪比公式，计算各个分辨率下径流深的空间插值信息量、噪声量和信噪比，如表4所示。

表4：不同输出分辨率下径流空间插值信息量、噪声量及信噪比。

分辨率信息量噪声量信噪比5326367161.54.55714610308896218.44.9673551525575.034585.95.57688420214963495.56.1496212515087.3623726.3606073096041368.37.018929356025.96742.38.117958403441.7786045.698308151482.253004.940833504842002.42

对各个分辨率下的空间插值信噪比进行排序后，绘制空间插值信息量、噪声量及信噪比与分辨率的曲线图，分别如图12、13、14所示，从图12、13、14中可见，当输出栅格大小依次选取5km、10km、15km、20km、25km、30km、35km、40km、45km、50km时，随着分辨率的降低，对于年平均径流深来说，无论是信息量还是噪声量都在减少，而各个分辨率下的空间插值信噪比则出现了先增后减的趋势，且在35km分辨率下，径流空间插值信噪比达到最大，径流空间插值信噪比SNR_Q_MAX＝8，即在此分辨率下，相邻栅格的差异程度是栅格插值平均误差的8倍。

结合实施例一、二、三可知空间插值信噪比的一般变化规律，即对于水文循环三要素降雨、蒸发和径流来说，都存在着相同的情况，即随着分辨率的降低空间插值信息量越少，同时空间插值噪声量也在减少，但由于两者的减少幅度、速度不同，使得最终的空间插值信噪比出现了起伏的趋势，且存在一个极大的峰值。而空间插值信噪比的峰值反映的是在此分辨率下，栅格的变化幅度(信息量)与栅格的平均插值误差(噪声量)之间的比值最大，此时的空间插值信噪比则为“最大信噪比”，此时的输出分辨率为“最优分辨率”。

另一方面，实施例一、二、三进行空间插值所用到的站点数目分别为145、13、20个，而它们相对应的空间分辨率阈值则分别出现在15km、80km、35km，由此我们可以猜测最优空间分辨率与站点数目有关，即同一地区站点越少，其所能得到的最优分辨率就越低；但并非仅由站点数目的反映，所以本发明名仅适用在同一地区站点数目相同下，确定最优分辨率。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。