基于遥感土壤湿度数据的水文模型参数估计方法及设备

1.本发明实施例涉及水文模型技术领域，尤其涉及一种基于遥感土壤湿度数据的水文模型参数估计方法及设备。


背景技术：

2.流域水文模型是一种用数学物理方法定量表达流域内降雨-径流关系的工具，常被用于水文预报、水库调度和水资源规划与管理等领域的科学研究和生产实践中。模型模拟精度与模型参数具有密切关系，传统的模型参数优选方法依赖于实测水文资料。但是，目前在世界范围内仍存在着大量无水文观测资料或缺资料的流域，导致流域水文模型很难在无资料地区取得满意的模拟精度。因此，开发一种基于遥感土壤湿度数据的水文模型参数估计方法及设备，可以有效克服上述相关技术中的缺陷，就成为业界亟待解决的技术问题。


技术实现要素：

3.针对现有技术存在的上述问题，本发明实施例提供了一种基于遥感土壤湿度数据的水文模型参数估计方法及设备。
4.第一方面，本发明的实施例提供了一种基于遥感土壤湿度数据的水文模型参数估计方法，包括：s1：获取参考流域范围内的遥感土壤湿度数据，遥感反演土壤湿度数据仅包含土壤表层范围内的土壤含水量信息，需要计算土壤湿度指数以反映整个土壤剖面的含水量情况；s2：根据遥感土壤湿度数据，计算参考流域和目标流域范围内各栅格处的土壤蓄水容量，绘制其频率分布直方图，选择概率分布线型作为蓄水容量曲线的线型，并进行k-s检验；s3：采用矩估计法，计算出蓄水容量曲线的概率分布参数的前l阶矩μ1，
…
，μ
l
，并联立μ1，
…
，μ
l
求出蓄水容量曲线概率分布参数的估计值；s4：对蓄水容量曲线参数采用估计值，集总式水文模型的其余参数采用实测水文气象数据进行参数自动率定，驱动模型对参考流域的降雨-径流过程进行模拟，并与通过常规率定方法得到的模拟结果进行对比；s5：在目标流域与参考流域之间，建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系，基于所述关系计算目标流域的蓄水容量曲线概率分布的参数；s6：对目标流域的蓄水容量曲线参数采用计算值，集总式水文模型的其余参数直接移植参考流域的参数值，并与移植参考流域全部参数值的移植方案得到的模拟结果进行对比。
5.在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s1中的土壤湿度指数的计算方式为：
[0006][0007]
式(1)中，为tn时刻i栅格位置的土壤湿度指数；为tn时刻i栅格位置的表层土壤湿度的观测值；kn为取值范围为0～1的参数，计算公式为：
[0008][0009]
式(2)中，t0代表土壤湿度变化的时间尺度，以天为单位，当预测0～100cm深度的土壤湿度值时，取t0＝20。
[0010]
在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s2中：
[0011]
a1、土壤蓄水容量的计算方式为：
[0012]
定义i栅格处的土壤蓄水容量为sri，通过以下方式计算：
[0013][0014][0015][0016]
(3)、(4)、(5)三式中，θi(t)均为流域内i栅格处根系层的遥感土壤湿度数据时间序列；max为取最大值；min为取最小值；t为时间；t为序列总长度；l为根系层在垂直方向的深度；
[0017]
a2、蓄水容量曲线：
[0018]
流域蓄水容量曲线是对流域饱和缺水量空间不均匀性的一种统计学描述，是将流域内各地点包气带的蓄水容量，按从小到大的顺序排列得到的一条蓄水容量与相应面积关系的统计曲线，
[0019]
a3、选择合适的概率分布线型：
[0020]
概率分布线型的选择基于三点根据：理论分布曲线的线型要尽量接近计算得到的蓄水容量分布特点；根据蓄水容量的性质，理论分布函数的自变量取值范围大于或大于等于0；由于理论分布的累积分布函数需要代入模型中参与计算，故理论分布函数的参数不宜过多；
[0021]
a4、k-s检验：
[0022]
kolmogorov-smirnov检验是一种非参数检验方法，计算经验累积分布与参考的理论累积分布函数之间的绝对差，绝对差的最大值就是k-s检验的统计量d
ks
，数学表达式如下：
[0023][0024]
式(6)中，f0为假定的理论累积分布函数；f
data
为样本的经验累积分布函数；n为样本总量，
[0025]
当k-s检验的统计量d
ks
符合(7)式时，将拒绝原假设，即样本不服从假定的理论分布：
[0026][0027]
式(7)中，α为显著性水平，取0.05；c(α)通过α、n查表得出；n为样本总量；对k-s检验的结果将使用p值法，若计算出的p值小于指定的显著性水平0.05，则拒绝原假设。
[0028]
在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s3中的矩估计法：
[0029]
对于服从某概率分布的连续型随机变量x，定义其单变量概率密度函数为f(x；u1，
…
，u
l
)，其累积分布函数为f(x；u1，
…
，u
l
)，u1，
…
，u
l
为待估参数，随机变量x的k阶矩μk的计算方式为：
[0030]
μk＝e(xk)＝∫xkdf(x；u1，
…
，u
l
)＝∫xkf(x；u1，
…
，u
l
)dx
ꢀꢀꢀ
(8)
[0031]
由式(8)可以得到前l阶矩的表达式：
[0032][0033]
联立求解出l个待估参数的表达式：
[0034][0035]
设(x1，x2，...，xn)为x的样本，样本的k阶矩的计算方式为：
[0036][0037]
式(11)中，n为样本总量，
[0038]
将样本的k阶矩ak作为μk的估计量，代入式(10)的fi中，得到ui的估计量。
[0039]
在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s4中：
[0040]
b1、常规率定方法：
[0041]
常规率定方法指，采用参数优化算法对集总式水文模型的所有参数进行率定，
[0042]
b2、结果对比：
[0043]
计算两种方法的nse效率系数和kge效率系数，
[0044]
nse的计算方式为：
[0045][0046]
式(12)中，q
sim，j
和q
obs，j
分别为第j时段的径流模拟值和实测值；为整个计算时段内径流实测值的均值；n为时段总长度；
[0047]
kge的计算方式为：
[0048][0049]
式(13)中，r为径流模拟值序列相对于观测值序列的pearson相关系数；a为径流模拟值序列的标准差与观测值序列的标准差的比值；b为径流模拟值序列的平均值与观测值序列的平均值的比值；
[0050]
根据nse值和kge值，分别定义两个模型稳健性指标mrc
nse
和mrc
kge
，计算方式分别为：
[0051][0052][0053]
式(14)、(15)中，nsev(n)和kgev(n)分别为检验期第n组参数值计算的nse结果和kge结果；nsec(n)和kgec(n)分别为率定期第n组参数值计算的nse结果和kge结果；n为参数组合的总数，mrc的结果可能为正也可能为负，数值越大模型的稳健性越好。
[0054]
在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s5中的建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系：
[0055]
对于流域蓄水容量的计算，满足(16)式：
[0056]
wm＝∫
0+∞
(1-f(c))dc＝e(c)
ꢀꢀꢀꢀꢀ
(16)
[0057]
式(16)中，wm为流域总蓄水量；c为流域内各点的蓄水容量；f(c)为蓄水容量曲线的累积分布函数；e(c)为流域内每个点的蓄水容量的统计期望；
[0058]
从式(16)可以看出，在数值上，流域的总蓄水量wm等于流域内每个点的蓄水容量c的统计期望，基于这种关系，对于目标流域和参考流域之间，都可以建立如下关系：
[0059]
[0060]
式(17)中，wma为目标流域的总蓄水量；wmb为参考流域的总蓄水量；和分别为目标流域和参考流域根据遥感土壤湿度数据计算得到的土壤蓄水容量均值；通过式(17)中的关系，利用遥感土壤湿度数据将参考流域的参数移植到目标流域。
[0061]
在上述方法实施例内容的基础上，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s6中结果对比：计算两种参数移植方案的nse值和kge值，计算方式分别参考式(12)和式(13)。
[0062]
第二方面，本发明的实施例提供了一种基于遥感土壤湿度数据的水文模型参数估计装置，包括：第一主模块，用于获取参考流域范围内的遥感土壤湿度数据，遥感反演土壤湿度数据仅包含土壤表层范围内的土壤含水量信息，需要计算土壤湿度指数以反映整个土壤剖面的含水量情况；第二主模块，用于根据遥感土壤湿度数据，计算参考流域和目标流域范围内各栅格处的土壤蓄水容量，绘制其频率分布直方图，选择概率分布线型作为蓄水容量曲线的线型，并进行k-s检验；第三主模块，用于采用矩估计法，计算出蓄水容量曲线的概率分布参数的前l阶矩μ1，
…
，μ
l
，并联立μ1，
…
，μ
l
求出蓄水容量曲线概率分布参数的估计值；第四主模块，用于对蓄水容量曲线参数采用估计值，集总式水文模型的其余参数采用实测水文气象数据进行参数自动率定，驱动模型对参考流域的降雨-径流过程进行模拟，并与通过常规率定方法得到的模拟结果进行对比；第五主模块，用于在目标流域与参考流域之间，建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系，基于所述关系计算目标流域的蓄水容量曲线概率分布的参数；第六主模块，用于对目标流域的蓄水容量曲线参数采用计算值，集总式水文模型的其余参数直接移植参考流域的参数值，并与移植参考流域全部参数值的移植方案得到的模拟结果进行对比。
[0063]
第三方面，本发明的实施例提供了一种电子设备，包括：
[0064]
至少一个处理器；以及
[0065]
与处理器通信连接的至少一个存储器，其中：
[0066]
存储器存储有可被处理器执行的程序指令，处理器调用程序指令能够执行第一方面的各种实现方式中任一种实现方式所提供的基于遥感土壤湿度数据的水文模型参数估计方法。
[0067]
第四方面，本发明的实施例提供了一种非暂态计算机可读存储介质，非暂态计算机可读存储介质存储计算机指令，计算机指令使计算机执行第一方面的各种实现方式中任一种实现方式所提供的基于遥感土壤湿度数据的水文模型参数估计方法。
[0068]
本发明实施例提供的基于遥感土壤湿度数据的水文模型参数估计方法及设备，对于参考流域，需要具有较长时间序列的实测水文资料；对于目标流域，需要在地形地貌、植被、土地利用和水文气候特征等方面与参考流域具有一定相似性，为水利行业从事无资料流域的水文预报、流域水资源利用与管理决策提供支撑，具有很强的实用性和广泛的适用性。
附图说明
[0069]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做一简单的介绍，显而易见地，下面描述中的附图是本发
明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0070]
图1为本发明实施例提供的基于遥感土壤湿度数据的水文模型参数估计方法流程图；
[0071]
图2为本发明实施例提供的基于遥感土壤湿度数据的水文模型参数估计装置结构示意图；
[0072]
图3为本发明实施例提供的电子设备的实体结构示意图；
[0073]
图4为本发明实施例提供的赣江流域(左侧)、湘江流域(右侧)的土壤蓄水容量空间分布结果示意图；
[0074]
图5为本发明实施例提供的赣江流域(左侧)、湘江流域(右侧)的蓄水容量的频率分布直方图与两种理论分布曲线的拟合结果示意图；
[0075]
图6为本发明实施例提供的hd_p1、hd_p2、hd_ln1和hd_ln2这四个模型在率定期的赣江流域径流模拟结果的nse和kge的箱线图；
[0076]
图7为本发明实施例提供的hd_p1、hd_p2、hd_ln1和hd_ln2这四个模型在检验期的赣江流域径流模拟结果的nse和kge的箱线图；
[0077]
图8为本发明实施例提供的湘江流域在情景1和情景2下对径流过程的模拟结果对比示意图；
[0078]
图9为本发明实施例提供的湘江流域在情景3和情景4下对径流过程的模拟结果对比示意图；
[0079]
图10为本发明实施例提供的湘江流域在四种情景下的nse和kge的箱线图。
具体实施方式
[0080]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。另外，本发明提供的各个实施例或单个实施例中的技术特征可以相互任意结合，以形成可行的技术方案，这种结合不受步骤先后次序和/或结构组成模式的约束，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时，应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0081]
本发明实施例提供了一种基于遥感土壤湿度数据的水文模型参数估计方法，参见图1，该方法包括：s1：获取参考流域范围内的遥感土壤湿度数据，遥感反演土壤湿度数据仅包含土壤表层(地表以下0～5cm)范围内的土壤含水量信息，需要计算土壤湿度指数以反映整个土壤剖面的含水量情况；s2：根据遥感土壤湿度数据，计算参考流域和目标流域范围内各栅格处的土壤蓄水容量，绘制其频率分布直方图，选择概率分布线型作为蓄水容量曲线的线型，并进行k-s检验(即kolmogorov-smirnov检验)；s3：采用矩估计法，计算出蓄水容量曲线的概率分布参数的前l阶矩μ1，
…
，μ
l
，并联立μ1，
…
，μ
l
求出蓄水容量曲线概率分布参数的估计值；s4：对蓄水容量曲线参数采用估计值，集总式水文模型的其余参数采用实测水文气象数据进行参数自动率定，驱动模型对参考流域的降雨-径流过程进行模拟，并与通过常
规率定方法得到的模拟结果进行对比；s5：在目标流域与参考流域之间，建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系，基于所述关系计算目标流域的蓄水容量曲线概率分布的参数；s6：对目标流域的蓄水容量曲线参数采用计算值，集总式水文模型的其余参数直接移植参考流域的参数值，并与移植参考流域全部参数值的移植方案得到的模拟结果进行对比。
[0082]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s1中的土壤湿度指数的计算方式为：
[0083][0084]
式(1)中，为tn时刻i栅格位置的土壤湿度指数；为tn时刻i栅格位置的表层土壤湿度的观测值；kn为取值范围为0～1的参数，计算公式为：
[0085][0086]
式(2)中，t0代表土壤湿度变化的时间尺度，以天为单位，当预测0～100cm深度的土壤湿度值时，取t0＝20。
[0087]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s2中：
[0088]
a1、土壤蓄水容量的计算方式为：
[0089]
定义i栅格处的土壤蓄水容量为sri，通过以下方式计算：
[0090][0091][0092][0093]
(3)、(4)、(5)三式中，θi(t)均为流域内i栅格处根系层的遥感土壤湿度数据时间序列；max为取最大值；min为取最小值；t为时间；t为序列总长度；l为根系层在垂直方向的深度；
[0094]
a2、蓄水容量曲线：
[0095]
流域蓄水容量曲线是对流域饱和缺水量空间不均匀性的一种统计学描述，是将流域内各地点包气带的蓄水容量，按从小到大的顺序排列得到的一条蓄水容量与相应面积关系的统计曲线，
[0096]
a3、选择合适的概率分布线型：
[0097]
概率分布线型的选择基于三点根据：理论分布曲线的线型要尽量接近计算得到的蓄水容量分布特点；根据蓄水容量的性质，理论分布函数的自变量取值范围大于或大于等于0；由于理论分布的累积分布函数需要代入模型中参与计算，故理论分布函数的参数不宜
过多；
[0098]
a4、k-s检验：
[0099]
kolmogorov-smirnov检验是一种非参数检验方法，计算经验累积分布与参考的理论累积分布函数之间的绝对差，绝对差的最大值就是k-s检验的统计量d
ks
，数学表达式如下：
[0100][0101]
式(6)中，f0为假定的理论累积分布函数；f
data
为样本的经验累积分布函数；n为样本总量，
[0102]
当k-s检验的统计量d
ks
符合(7)式时，将拒绝原假设，即样本不服从假定的理论分布：
[0103][0104]
式(7)中，α为显著性水平，取0.05；c(α)通过α、n查表得出；n为样本总量；对k-s检验的结果将使用p值法，若计算出的p值小于指定的显著性水平0.05，则拒绝原假设。
[0105]
具体如图4所示(横竖坐标为经纬度)，根据smap土壤湿度数据，计算出赣江流域、湘江流域的土壤蓄水容量空间分布结果。
[0106]
如图5所示，绘制了赣江流域(图5左侧)、湘江流域(图5右侧)的蓄水容量频率分布直方图对两种理论分布曲线的拟合结果。观察直方图，可以看出无论是赣江流域还是湘江流域，其蓄水容量都是向中间集中，而两端较少。这种分布情况不再符合hymod模型中原始的蓄水容量曲线的pareto分布。由选择线型的三点根据，认为lognormal分布更适合作为两个流域的蓄水容量曲线的线型。
[0107]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s3中的矩估计法：
[0108]
对于服从某概率分布的连续型随机变量x，定义其单变量概率密度函数为f(x；u1，
…
，u
l
)，其累积分布函数为f(x；u1，
…
，u
l
)，u1，
…
，ui为待估参数，随机变量x的k阶矩μk的计算方式为：
[0109]
μk＝e(xk)＝∫xkdf(x；u1，
…
，u
l
)＝∫xkf(x；u1，
…
，u
l
)dx
ꢀꢀꢀ
(8)
[0110]
由式(8)可以得到前l阶矩的表达式：
[0111][0112]
联立求解出l个待估参数的表达式：
[0113][0114]
设(x1，x2，...，xn)为x的样本，样本的k阶矩的计算方式为：
[0115][0116]
式(11)中，n为样本总量，
[0117]
将样本的k阶矩ak作为μk的估计量，代入式(10)的fi中，得到ui的估计量。
[0118]
具体地，pareto分布有两个参数，分别为sm和β。根据式(8)计算出pareto分布的一阶矩和二阶矩并联立μ1、μ2求出线型参数的估计值根据式(5)，得到β的估计值结果
[0119]
lognormal分布有两个参数，分别为μ和σ。根据式(8)计算出lognormal分布的一阶矩和二阶矩μ2＝exp(2μ+2σ2)，并联立μ1、μ2求出线型参数的估计值根据式(5)，得到σ的估计值结果
[0120]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s4中：
[0121]
b1、常规率定方法：
[0122]
常规率定方法指，采用参数优化算法对集总式水文模型的所有参数进行率定，
[0123]
b2、结果对比：
[0124]
计算两种方法的nse效率系数和kge效率系数，
[0125]
nse的计算方式为：
[0126][0127]
式(12)中，q
sim，j
和q
obs，j
分别为第j时段的径流模拟值和实测值；为整个计算时段内径流实测值的均值；n为时段总长度；
[0128]
kge的计算方式为：
[0129][0130]
式(13)中，r为径流模拟值序列相对于观测值序列的pearson相关系数；a为径流模拟值序列的标准差与观测值序列的标准差的比值；b为径流模拟值序列的平均值与观测值序列的平均值的比值；
[0131]
根据nse值和kge值，分别定义两个模型稳健性指标mrc
nse
和mrc
kge
，计算方式分别为：
[0132][0133][0134]
式(14)、(15)中，nsev(n)和kgev(n)分别为检验期第n组参数值计算的nse结果和kge结果；nsec(n)和kgec(n)分别为率定期第n组参数值计算的nse结果和kge结果；n为参数组合的总数，mrc的结果可能为正也可能为负，数值越大模型的稳健性越好。
[0135]
具体地，参数优化算法采用dream算法，以流量的均方根误差rmse作为目标函数。模型参数的先验分布选择基于其取值上下界的均匀分布，在dream算法结束后，从参数的后验分布中选取最后的2000组采样结果作为参数后验分布的样本，以考虑模型参数的不确定性。从后验分布中选择一组在率定期目标函数值最小的一组参数值作为最优参数。
[0136]
对于采用原始的蓄水容量曲线(即pareto分布)的hymod模型，所有参数全部由dream算法率定得到的hymod模型记为hd_p1；参数β采用估计值，其余参数都由dream算法率定得到的hymod模型记为hd_p2。
[0137]
对于采用lognormal分布作为蓄水容量曲线的hymod模型，全部参数都由dream算法率定得到的hymod模型记为hd_ln1；参数σ采用估计值，其余参数都由dream算法率定得到的hymod模型记为hd_ln2。
[0138]
图6展示了hd_p1、hd_p2、hd_ln1和hd_ln2这四个模型在率定期的赣江流域径流模拟结果的nse和kge的箱线图。由图可以看出，在率定期，对于hd_p1，使用了估计得到的参数
后，模拟效果下降并不明显；对于hd_ln1，使用了估计得到的参数后，模拟效果有所提升。
[0139]
图7展示了hd_p1、hd_p2、hd_ln1和hd_ln2这四个模型在检验期的赣江流域径流模拟结果的nse和kge的箱线图。可以看出在检验期，对于hd_p1，使用了估计得到的参数后，模拟效果下降并不明显；对于hd_ln1，使用了估计得到的参数后，模拟效果有所提升。
[0140]
可以看出，在率定期和检验期，hd_p1、hd_p2、hd_ln1和hd_ln2这四个模型均能较好地模拟出赣江流域的日径流过程。hd_p1和hd_p2的模拟效果和不确定性范围相差不大；hd_ln2在率定期和检验期的模拟效果均比hd_ln1更优，不确定性范围也更小。
[0141]
根据(14)、(15)两式，计算出hd_ln1和hd_ln2的两种稳健性指标：hd_ln1的mrc
nse
为1.93％，mrc
kge
为-0.83％；hd_ln2的mrc
nse
为2.39％，mrc
kge
为-0.08％。可以看出，本发明的参数估计方法能够提高水文模型的稳健性。
[0142]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计方法，步骤s5中的建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系：
[0143]
对于流域蓄水容量的计算，满足(16)式：
[0144]
wm＝∫
0+∞
(1-f(c))dc＝e(c)
ꢀꢀꢀꢀꢀꢀ
(16)
[0145]
式(16)中，wm为流域总蓄水量；c为流域内各点的蓄水容量；f(c)为蓄水容量曲线的累积分布函数；e(c)为流域内每个点的蓄水容量的统计期望；
[0146]
从式(16)可以看出，在数值上，流域的总蓄水量wm等于流域内每个点的蓄水容量c的统计期望，基于这种关系，对于目标流域和参考流域之间，都可以建立如下关系：
[0147][0148]
式(17)中，wma为目标流域的总蓄水量；wmb为参考流域的总蓄水量；和分别为目标流域和参考流域根据遥感土壤湿度数据计算得到的土壤蓄水容量均值；通过式(17)中的关系，利用遥感土壤湿度数据将参考流域的参数移植到目标流域。
[0149]
具体地，对于使用原始的蓄水容量曲线(即pareto分布)的hymod模型，根据式(17)可以构建目标流域与参考流域之间的总蓄水量与蓄水容量均值之间的关系为经过整理，可以得到sma的计算式为
[0150]
对于蓄水容量曲线使用lognormal分布的hymod模型，根据式(17)可以构建目标流域与参考流域之间的总蓄水量与蓄水容量均值之间的关系为经过整理，可以得到μa的计算式为
[0151]
基于上述方法实施例的内容，作为一种可选的实施例，本发明实施例中提供的基
于遥感土壤湿度数据的水文模型参数估计方法，步骤s6中结果对比：计算两种参数移植方案的nse值和kge值，计算方式分别参考式(12)和式(13)。
[0152]
具体地，对于无资料流域的水文模拟，不存在参数率定过程，因此模拟时间设置为整个有资料时期(即2011年至2017年)，其中2011年作为“预热期”，不参与结果对比。
[0153]
将hd_p2模型在赣江流域率定得到的参数，全部移植到湘江流域，记为情景1；利用hd_p2模型在赣江流域率定得到的参数计算出湘江流域的参数sm，再将hd_p2模型的剩余参数全部移植至湘江流域，记为情景2；将hd_ln2模型在赣江流域的率定得到参数，全部移植到湘江流域，记为情景3；利用hd_ln2模型在赣江流域率定得到的参数计算出湘江流域的参数μ，再将hd_ln2模型的剩余参数全部移植至湘江流域，记为情景4。
[0154]
如图8所示，展示了湘江流域在情景1和情景2下对径流过程的模拟结果对比。可以看出，情景2的模拟序列在低水部分和高水部分都优于情景1，且两个情景的不确定性区间的带宽接近。
[0155]
如图9所示，展示了湘江流域在情景3和情景4下对径流过程的模拟结果对比。可以看出，情景4的不确定性区间带宽明显大于情景3。
[0156]
如图10所示，展示了湘江流域在四种情景下的nse和kge的箱线图。可以看出，无论是对于nse还是kge，情景2的模拟结果都比情景1更优，不确定性范围变化不大；情景4比情景3的总体模拟结果更优，但不确定性范围明显增大。
[0157]
本发明实施例提供的基于遥感土壤湿度数据的水文模型参数估计方法，对于参考流域，需要具有较长时间序列的实测水文资料；对于目标流域，需要在地形地貌、植被、土地利用和水文气候特征等方面与参考流域具有一定相似性，为水利行业从事无资料流域的水文预报、流域水资源利用与管理决策提供支撑，具有很强的实用性和广泛的适用性。
[0158]
本发明各个实施例的实现基础是通过具有处理器功能的设备进行程序化的处理实现的。因此在工程实际中，可以将本发明各个实施例的技术方案及其功能封装成各种模块。基于这种现实情况，在上述各实施例的基础上，本发明的实施例提供了一种基于遥感土壤湿度数据的水文模型参数估计装置，该装置用于执行上述方法实施例中的基于遥感土壤湿度数据的水文模型参数估计方法。参见图2，该装置包括：第一主模块，用于获取参考流域范围内的遥感土壤湿度数据，遥感反演土壤湿度数据仅包含土壤表层范围内的土壤含水量信息，需要计算土壤湿度指数以反映整个土壤剖面的含水量情况；第二主模块，用于根据遥感土壤湿度数据，计算参考流域和目标流域范围内各栅格处的土壤蓄水容量，绘制其频率分布直方图，选择概率分布线型作为蓄水容量曲线的线型，并进行k-s检验；第三主模块，用于采用矩估计法，计算出蓄水容量曲线的概率分布参数的前l阶矩μ1，
…
，μ
l
，并联立μ1，
…
，μ
l
求出蓄水容量曲线概率分布参数的估计值；第四主模块，用于对蓄水容量曲线参数采用估计值，集总式水文模型的其余参数采用实测水文气象数据进行参数自动率定，驱动模型对参考流域的降雨-径流过程进行模拟，并与通过常规率定方法得到的模拟结果进行对比；第五主模块，用于在目标流域与参考流域之间，建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系，基于所述关系计算目标流域的蓄水容量曲线概率分布的参数；第六主模块，用于对目标流域的蓄水容量曲线参数采用计算值，集总式水文模型的其余参数直接移植参考流域的参数值，并与移植参考流域全部参数值的移植方案得到的模拟结果进行对比。
[0159]
本发明实施例提供的基于遥感土壤湿度数据的水文模型参数估计装置，采用图2中的若干模块，对于参考流域，需要具有较长时间序列的实测水文资料；对于目标流域，需要在地形地貌、植被、土地利用和水文气候特征等方面与参考流域具有一定相似性，为水利行业从事无资料流域的水文预报、流域水资源利用与管理决策提供支撑，具有很强的实用性和广泛的适用性。
[0160]
需要说明的是，本发明提供的装置实施例中的装置，除了可以用于实现上述方法实施例中的方法外，还可以用于实现本发明提供的其他方法实施例中的方法，区别仅仅在于设置相应的功能模块，其原理与本发明提供的上述装置实施例的原理基本相同，只要本领域技术人员在上述装置实施例的基础上，参考其他方法实施例中的具体技术方案，通过组合技术特征获得相应的技术手段，以及由这些技术手段构成的技术方案，在保证技术方案具备实用性的前提下，就可以对上述装置实施例中的装置进行改进，从而得到相应的装置类实施例，用于实现其他方法类实施例中的方法。例如：
[0161]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第一子模块，用于实现步骤s1中的土壤湿度指数的计算方式为：
[0162][0163]
式(1)中，为tn时刻i栅格位置的土壤湿度指数；为tn时刻i栅格位置的表层土壤湿度的观测值；kn为取值范围为0～1的参数，计算公式为：
[0164][0165]
式(2)中，t0代表土壤湿度变化的时间尺度，以天为单位，当预测0～100cm深度的土壤湿度值时，取t0＝20。
[0166]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第二子模块，用于实现步骤s2中：
[0167]
a1、土壤蓄水容量的计算方式为：
[0168]
定义i栅格处的土壤蓄水容量为sri，通过以下方式计算：
[0169][0170][0171][0172]
(3)、(4)、(5)三式中，θi(t)均为流域内i栅格处根系层的遥感土壤湿度数据时间序列；max为取最大值；min为取最小值；t为时间；t为序列总长度；l为根系层在垂直方向的深度；
[0173]
a2、蓄水容量曲线：
[0174]
流域蓄水容量曲线是对流域饱和缺水量空间不均匀性的一种统计学描述，是将流域内各地点包气带的蓄水容量，按从小到大的顺序排列得到的一条蓄水容量与相应面积关系的统计曲线，
[0175]
a3、选择合适的概率分布线型：
[0176]
概率分布线型的选择基于三点根据：理论分布曲线的线型要尽量接近计算得到的蓄水容量分布特点；根据蓄水容量的性质，理论分布函数的自变量取值范围大于或大于等于0；由于理论分布的累积分布函数需要代入模型中参与计算，故理论分布函数的参数不宜过多；
[0177]
a4、k-s检验：
[0178]
kolmogorov-smirnov检验是一种非参数检验方法，计算经验累积分布与参考的理论累积分布函数之间的绝对差，绝对差的最大值就是k-s检验的统计量d
ks
，数学表达式如下：
[0179][0180]
式(6)中，f0为假定的理论累积分布函数；f
data
为样本的经验累积分布函数；n为样本总量，
[0181]
当k-s检验的统计量d
ks
符合(7)式时，将拒绝原假设，即样本不服从假定的理论分布：
[0182][0183]
式(7)中，α为显著性水平，取0.05；c(α)通过α、n查表得出；n为样本总量；对k-s检验的结果将使用p值法，若计算出的p值小于指定的显著性水平0.05，则拒绝原假设。
[0184]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第三子模块，用于实现步骤s3中的矩估计法：
[0185]
对于服从某概率分布的连续型随机变量x，定义其单变量概率密度函数为f(x；u1，
…
，u
l
)，其累积分布函数为f(x；u1，
…
，u
l
)，u1，
…
，u
l
为待估参数，随机变量x的k阶矩μk的计算方式为：
[0186]
μk＝e(xk)＝∫xkdf(x；u1，
…
，u
l
)＝∫xkf(x；u1，
…
，u
l
)dx
ꢀꢀꢀꢀ
(8)
[0187]
由式(8)可以得到前l阶矩的表达式：
[0188][0189]
联立求解出l个待估参数的表达式：
[0190][0191]
设(x1，x2，...，xn)为x的样本，样本的k阶矩的计算方式为：
[0192][0193]
式(11)中，n为样本总量，
[0194]
将样本的k阶矩ak作为μk的估计量，代入式(10)的fi中，得到ui的估计量。
[0195]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第四子模块，用于实现步骤s4中：
[0196]
b1、常规率定方法：
[0197]
常规率定方法指，采用参数优化算法对集总式水文模型的所有参数进行率定，
[0198]
b2、结果对比：
[0199]
计算两种方法的nse效率系数和kge效率系数，
[0200]
nse的计算方式为：
[0201][0202]
式(12)中，q
sim，j
和q
obs，j
分别为第j时段的径流模拟值和实测值；为整个计算时段内径流实测值的均值；n为时段总长度；
[0203]
kge的计算方式为：
[0204][0205]
式(13)中，r为径流模拟值序列相对于观测值序列的pearson相关系数；a为径流模拟值序列的标准差与观测值序列的标准差的比值；b为径流模拟值序列的平均值与观测值序列的平均值的比值；
[0206]
根据nse值和kge值，分别定义两个模型稳健性指标mrc
nse
和mrc
kge
，计算方式分别为：
[0207]
[0208][0209]
式(14)、(15)中，nsev(n)和kgev(n)分别为检验期第n组参数值计算的nse结果和kge结果；nsec(n)和kgec(n)分别为率定期第n组参数值计算的nse结果和kge结果；n为参数组合的总数，mrc的结果可能为正也可能为负，数值越大模型的稳健性越好。
[0210]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第五子模块，用于实现步骤s5中的建立总蓄水量和通过遥感土壤湿度数据计算得到的土壤蓄水容量均值的关系：
[0211]
对于流域蓄水容量的计算，满足(16)式：
[0212]
wm＝∫
0+∞
(1-f(c))dc＝e(c)
ꢀꢀꢀꢀ
(16)
[0213]
式(16)中，wm为流域总蓄水量；c为流域内各点的蓄水容量；f(c)为蓄水容量曲线的累积分布函数；e(c)为流域内每个点的蓄水容量的统计期望；
[0214]
从式(16)可以看出，在数值上，流域的总蓄水量wm等于流域内每个点的蓄水容量c的统计期望，基于这种关系，对于目标流域和参考流域之间，都可以建立如下关系：
[0215][0216]
式(17)中，wma为目标流域的总蓄水量；wmb为参考流域的总蓄水量；和分别为目标流域和参考流域根据遥感土壤湿度数据计算得到的土壤蓄水容量均值；通过式(17)中的关系，利用遥感土壤湿度数据将参考流域的参数移植到目标流域。
[0217]
基于上述装置实施例的内容，作为一种可选的实施例，本发明实施例中提供的基于遥感土壤湿度数据的水文模型参数估计装置，还包括：第六子模块，用于实现步骤s6中结果对比：计算两种参数移植方案的nse值和kge值，计算方式分别参考式(12)和式(13)。
[0218]
本发明实施例的方法是依托电子设备实现的，因此对相关的电子设备有必要做一下介绍。基于此目的，本发明的实施例提供了一种电子设备，如图3所示，该电子设备包括：至少一个处理器(processor)、通信接口(communications interface)、至少一个存储器(memory)和通信总线，其中，至少一个处理器，通信接口，至少一个存储器通过通信总线完成相互间的通信。至少一个处理器可以调用至少一个存储器中的逻辑指令，以执行前述各个方法实施例提供的方法的全部或部分步骤。
[0219]
此外，上述的至少一个存储器中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个方法实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read-only memory)、随机存取存储器(ram，random access memory)、磁碟或者光盘等各种可以存储程序代码的介质。
[0220]
以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。
[0221]
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件实现。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如rom/ram、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。
[0222]
附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。基于这种认识，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
[0223]
在本专利中，术语
″
包括
″
、
″
包含
″
或者其任何其它变体意在涵盖非排它性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其它要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句
″
包括
……″
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
[0224]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。