流域径流的模拟方法

1.本发明属于流域水文模拟的技术领域，具体涉及一种流域径流的模拟方法。


背景技术：

2.准确模拟流域降雨-径流过程和水平衡对于未来变化环境下的水资源规划和管理至关重要。集总式水文模型(lumped hydrological models,lhms)因其简单高效是研究水资源和水环境的重要工具，如澳大利亚水平衡模型(australian water balance model，awbm)，gr2m模型，simhyd模型，水箱模型(tank)。然而，lhms聚焦降雨-径流过程，忽略或过度简化了蒸发过程，影响了流域水量平衡模拟的准确性和合理性。目前，水文学者已通过纳入成熟的蒸发理论改进lhms并取得了不同程度的效果，但由于需要额外的土壤、植被等信息或极大增加模型复杂度难以广泛推广应用。


技术实现要素：

3.本发明的目的在于针对现有技术的不足之处，提供一种流域径流的模拟方法，该方法解决了集总式水文模型过度简化蒸发过程的问题，提高了集总式水文模型的径流模拟能力和模型适用性。
4.为解决上述技术问题，本发明采用如下技术方案：
5.一种流域径流的模拟方法，包括如下步骤：
6.步骤1、获取研究流域的基础水文气象数据；
7.步骤2、基于获取的基础水文气象数据和广义蒸发互补关系，构建蒸发估算；
8.步骤3、将蒸发估算与集总式水文模型进行耦合获得耦合后的水文模型；
9.步骤4、构建参数率定的目标函数，采用优化算法结合实测数据率定耦合后的水文模型的参数；
10.步骤5、利用率定后的水文模型模拟流域径流过程，并评价水文模型的模拟精度。
11.进一步地，基础水文气象数据包括月尺度的降水p、径流q、净辐射rn、近地面2m处风速u2、气温t、饱和水汽压差d。
12.进一步地，净辐射rn根据数据集中采集的向下短波辐射r
si
、向下长波辐射r
li
、地表反照率α和气温t确定；
13.饱和水汽压差根据数据集中采集的气温t和相对湿度rh确定；
14.近地面2m处风速u2根据数据集中的水平风速u和垂向风速v确定。
15.进一步地，步骤2中表观潜在蒸发e
pa
的计算方式如下：
[0016][0017]
式中：δ指给定的空气温度下饱和蒸汽压力曲线的斜率；γ指干湿表常数；rn指陆地表面的净辐射；g指地面热通量,在日尺度及更长的时间尺度上可忽略不计；d指饱和水汽压差；f(u2)指风函数，其表达式如下表示：
[0018]
f(u2)＝0.26(1+0.54u2)；
[0019]
其中，u2表示近地面2m处风速。
[0020]
进一步地，步骤2中平衡蒸发ee的计算方式如下：
[0021][0022]
式中：δ指给定的空气温度下饱和蒸汽压力曲线的斜率；γ指干湿表常数；rn指陆地表面的净辐射；g指地面热通量。
[0023]
进一步地，步骤2中构建的蒸发估算为：
[0024][0025]
式中，αc为参数。
[0026]
进一步地，集总式水文模型包括awbm、gr2m、simhyd、tank模型中的一种。
[0027]
进一步地，步骤4中构建的目标函数为kge，其公式为：
[0028][0029]
式中：q
obs
指观测径流；q
sim
指模拟径流；分别指观测径流和模拟径流序列的平均值；σo、σs分别指观测径流和模拟径流序列的标准差；为i时刻的实测径流量。
[0030]
进一步地，步骤4中率定参数的方法为：
[0031]
步骤4.1初始化阶段；输入水文气象观测资料序列，指定模型参数可行范围，在可行范围内随机搜索确定一个参数集，可得出第1次模拟的径流序列{q
sim,1
(i)}，计算实测序列与模拟序列之间的目标函数值kge1；
[0032]
步骤4.2迭代阶段；从n＝2开始，优化算法在参数范围内选定较上次更优的参数集，得到第n次模拟的径流序列{q
sim,n
(i)}，计算实测序列与模拟序列之间的目标函数值kgen，直至kgen值趋于收敛；此时，保留水文模型参数，用于流域径流模拟。
[0033]
进一步地，步骤5中采用水文评价指标kge和r2评估率定后的水文模型的径流模拟精度。
[0034]
与现有技术相比，本发明的有益效果为：本发明考虑了更具有物理意义的实际蒸发模拟方法，解决了集总式水文模型过度简化蒸发过程的问题，从而提高了集总式水文模型的径流模拟能力和模型适用性；此外，本发明具有参数少、结构简单、计算高效等特点，因此具有较好的适用性。
附图说明
[0035]
图1为本发明实施例中四个集总式水文模型的结构示意图；
[0036]
图2为本发明实施例的流程图；
[0037]
图3为本发明实施例中径流模拟表现，图中方框的底部、中部和顶部分别代表25％、50％和75％分位数，误差条的顶部和底部指离群值的分界点；其中(a、b)为原始集总式水文模型的径流模拟表现，(c、d)实施例耦合水文模型的径流模拟表现，(e、f)改善的径
流模拟表现。
具体实施方式
[0038]
下面将结合本发明实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0039]
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
[0040]
下面结合具体实施例对本发明作进一步说明，但不作为本发明的限定。
[0041]
本发明实施例公开了一种流域径流的模拟方法，包括如下步骤：
[0042]
步骤1、获取研究流域的基础水文气象数据；
[0043]
在实施例中，收集全球径流数据中心(global runoff data centre，grdc,https://www.bafg.de/grdc/)和大样本研究的集水区属性和气象学(camels)数据集的2116个流域为研究实例，其中，研究流域涵盖了赤道气候、暖温带气候、干旱气候、雪地气候和极地气候五个基本的气候类型。获取研究流域的观测径流q(mm)、降水p(mm)和常规气象数据，包括净辐射rn(mm)、近地面2m处风速u2(m/s)、气温t(℃)、饱和水汽压差d(hpa)。
[0044]
其中，净辐射利用向下短波辐射r
si
(mm)、向下长波辐射r
li
(mm)、地表反照率α和气温t(℃)确定，实现方式如下：
[0045]rn
＝r
si-r
so
+r
li-r
lo
ꢀꢀꢀꢀꢀꢀ
(1)
[0046]rso
＝αr
si
ꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0047]rlo
＝εsσt
s4
ꢀꢀꢀꢀꢀꢀ
(3)
[0048]
式中，εs为地表反射率，一般取0.97；σ为斯蒂芬-玻尔兹曼(stefan-boltzmann)常数，4.903
×
10-9
mj
·
k-4
·
m-2
·
d-1
；ts为地表温度，用气温近似表示。
[0049]
其中，饱和水汽压差利用气温t(℃)和相对湿度rh(％)确定，实现方式如下：
[0050]
d＝10(e
s-ea)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0051][0052][0053]
式中，es为饱和水汽压(kpa)；ea为实际水汽压(kpa)。
[0054]
其中，近地面2m处风速u2(m/s)利用地面上10m高处水平风速u(m/s)和垂向风速v(m/s)确定，实现方式如下：
[0055][0056]
值得注意的是，camels数据集有降水数据可直接使用；grdc数据集无配套的降水数据，来源于cru-jra再分析数据集。向下短波辐射r
si
(mm)、向下长波辐射r
li
(mm)、气温t(℃)、相对湿度rh(％)、地面上10m高处水平风速u(m/s)和垂向风速v(m/s)等常规气象数据均来源于cru-jra v2.2再分析数据集。
[0057]
步骤2、基于步骤1获得的气象数据，利用广义蒸发互补关系构建流域蒸发估算；广义蒸发互补理论解释了陆面蒸发和大气状态之间的耦合关系，通过描述实际蒸发e(mm)、表观潜在蒸发e
pa
(mm)和潜在蒸发e
po
(mm)的非线性非对称关系计算流域蒸发。即：当干热空气团经过充分供水的区域上方时，引起陆面实际蒸散发量的增加，而自身的潜在蒸散发能力会降低。其中，e
po
可表征为可调参数αc与平衡蒸发ee(mm)的乘积；该步骤包括：
[0058]
步骤2.1、基于步骤1获得的气象数据，计算蒸发输入变量表观潜在蒸发e
pa
和平衡蒸发ee，实现方式如下：
[0059][0060][0061]
式中：δ指给定的空气温度下饱和蒸汽压力曲线的斜率(kpa℃-1
)；γ指干湿表常数(kpa℃-1
)；rn指陆地表面的净辐射(mm)；g指地面热通量(mm)，在日尺度及更长的时间尺度上可略不计；d指饱和水汽压差(hpa)；f(u2)指风函数，计算方式如下：
[0062]
f(u2)＝0.26(1+0.54u2)
ꢀꢀꢀꢀꢀꢀ
(10)
[0063]
步骤2.2利用广义蒸发互补关系及步骤2.1计算的蒸发变量构建流域的蒸发估算，实现方式如下：
[0064][0065]
式中：αc指可调参数，随流域下垫面条件及其他属性而变化。e
gcr
指采用gcr方法计算的流域实际蒸发(mm)。
[0066]
步骤3、耦合构建的蒸发估算方法与集总式水文模型，获得改进后的水文模型(lhms-gcr)；集总式水文模型基于流域水平衡(p＝e+q+δs)理论建立，直接采用表观潜在蒸发参与水平衡计算，产生了不大恰当的水平衡误差和径流模拟。lhms-gcr采用了gcr方法估计的流域实际蒸发，有效解决了改问题。其中，p为观测降水、q为观测径流、e为实际蒸发，δs为流域蓄水量变化量。
[0067]
在本实施例中，如图1和图2所示，以awbm、gr2m、simhyd、tank四个集总式水文模型为例进行说明如何将蒸发估算与集总式水文模型进行耦合。在获得蒸发估算后，构建耦合模型的产汇流模块，在一个时间步长内实现方式如下。
[0068]
对awbm-gcr模型，原始模型设置蒸发为土壤的线性耗水，且默认表观潜在蒸发等于实际蒸发参与水平衡计算；耦合模型利用e
gcr
替换e
pa
，其产汇流计算实现方式如下：
[0069]
sm＝a1
×
c1+a2
×
c2+a3
×
c3
ꢀꢀꢀꢀꢀꢀꢀ
(12)
[0070]ri
＝max(si+p
i-e
gcr,i-sm,0)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0071]
qi＝bfi
×ri
×
(1-kb)+(1-bfi)
×ri
×
(1-ks)
ꢀꢀꢀꢀꢀ
(14)
[0072]si+1
＝max(si+p
i-e
gcr,i
,sm)
ꢀꢀꢀꢀꢀꢀꢀ
(15)
[0073]
式中：ri指第i时刻的流域产流量(mm)；qi指第i时刻的流域径流量(mm)；si和s
i+1
分别指第i时刻和第i+1时刻的流域土壤含水量；pi指第i时刻的流域降水量；a1、a2、a3分别指三个独立地表储水库的面积比例，满足a1+a2+a3＝1；c1、c2、c3分别指三个独立地表储水库
的容量(mm)；sm指流域最大蓄水容量(mm)；bfi指基流指数；ks指地表径流衰退系数；kb指基流衰退系数。其中，a1、a2、c1、c2、c3、bfi、ks、kb为待校准参数。
[0074]
对gr2m-gcr模型，原始模型设置土壤的降水补给和蒸发耗水为非线性关系，并利用表观潜在蒸发e
pa
参与计算；耦合模型利用e
gcr
替换e
pa
，土壤的蒸发耗水调整为线性关系，其产汇流计算实现方式如下：
[0075][0076]
p1＝pi+s
i-s1ꢀꢀꢀꢀꢀꢀꢀ
(17)
[0077]
s2＝s
1-e
gcr
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0078][0079]
p2＝s
2-s
i+1
ꢀꢀꢀꢀꢀꢀꢀ
(20)
[0080][0081]ri+1
＝x
×
(ri+p1+p2)-qiꢀꢀꢀꢀꢀꢀꢀꢀ
(22)
[0082]
式中：s1指考虑降水的土壤水含量(mm)；s2指在s1基础上在考虑蒸发的土壤水含量(mm)；p1指补充土壤水分的剩余降雨量(mm)；p2指土壤蓄水库的过降雨量(mm)；sm指流域土壤水的最大蓄水容量(mm)；x指水交换系数。其中，sm、x为待校准参数。
[0083]
对simhyd-gcr模型，原始模型设置林冠截流模块、非线性土壤蒸发模块，但仍以表观潜在蒸发e
pa
参与计算；耦合模型利用e
gcr
替换e
pa
，土壤的蒸发耗水调整为线性关系，实际蒸发作为输入约束，林冠截流模块可被简略，其产汇流计算实现方式如下：
[0084]
inp＝max(p
i-e
gcr,i
,0)
ꢀꢀꢀꢀꢀꢀꢀ
(23)
[0085][0086][0087][0088]
orun＝max(si+inf-srun-rec-sm,0)
ꢀꢀꢀꢀꢀꢀ
(27)
[0089]
qi＝kb×
(gi+rec+orun)+(inp-inf)+srun
ꢀꢀꢀꢀꢀꢀꢀꢀ
(28)
[0090]si+1
＝min(min(si+inf-srun-rec,sm)-e
gcr,i
,0)
ꢀꢀꢀꢀꢀꢀꢀ
(29)
[0091]gi+1
＝(1-kb)
×gi
ꢀꢀꢀꢀꢀꢀꢀ
(30)
[0092]
式中：inp指有效降雨量(mm)；inf指土壤下渗量(mm)；coeff指最大入渗损失(mm)；sq指下渗损失指数；srun指壤中流(mm)；sub指壤中流出流系数；rec指地下水补给量(mm)；
crak指地下水补给的比例系数；orun指土壤过饱和水量，用以补给地下水(mm)；gi和g
i+1
分别指第i时刻和第i+1时刻的地下水蓄水量(mm)。其中，coeff、sq、sm、sub、crak、kb为待校准参数。
[0093]
对tank-gcr模型，原始模型设置蒸发为土壤的线性耗水，且默认表观潜在蒸发等于实际蒸发参与水平衡计算；耦合模型利用e
gcr
替换e
pa
，以2层水箱结构为例，其产汇流计算实现方式如下：
[0094]z1,i
＝min(max(z
1,i
+p
i-e
gcr,i
,0),hs1)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(31)
[0095]qs
＝max((z
1,i-h1)
×
a1,0)
ꢀꢀꢀꢀꢀꢀꢀ
(32)
[0096]
qi＝max((z
1,i-h2)
×
a2,0)
ꢀꢀꢀꢀꢀꢀꢀ
(33)
[0097]
inf＝max((z
1,i-q
s-qi)
×
b,0)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(34)
[0098]
qb＝(z
2,i
+inf)
×
a3ꢀꢀꢀꢀꢀꢀꢀꢀ
(35)
[0099]
qi＝qs+qi+qbꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(36)
[0100]z1,i+1
＝z
1,i-q
s-q
i-inf
ꢀꢀꢀꢀꢀꢀꢀ
(37)
[0101]z2,i+1
＝min(z
2,i
+inf-qb,hs2)
ꢀꢀꢀꢀꢀꢀ
(38)
[0102]
式中：z1、z2分别指上层和下层水箱的蓄水量(mm)；hs1、hs2分别指上层和下层水箱的最大蓄水量(mm)；h1、h2指上层水箱两个侧出口的高度(mm)；
[0103]
a1、a2、a3指三个侧出口的出流系数；b指上层水箱下出口的出流系数；qs、qi、qb分别指地表径流、壤中流和基流(mm)；其中，a1、a2、a3、b、hs1、hs2、h1、h2为待校准参数。
[0104]
步骤4、构建参数率定的目标函数,采用优化算法结合实测数据率定耦合后的水文模型的参数；
[0105]
选择同时考虑模拟径流序列与实测径流序列一致性、偏差和可变性的kling
–
gupta效率系数(kge)为目标函数，其公式为：
[0106][0107]
式中：式中：q
obs
指观测径流；q
sim
指模拟径流；分别指观测径流和模拟径流序列的平均值；σo、σs分别指观测径流和模拟径流序列的标准差；为i时刻的实测径流量。
[0108]
以kge为目标函数，采用优化算法率定模型参数具体包括：
[0109]
步骤4.1初始化阶段；输入水文气象观测资料序列，指定模型参数可行范围，算法在可行范围内随机搜索确定一个参数集，可得出第1次模拟的径流序列{q
sim,1
(i)}，计算实测序列与模拟序列之间的目标函数值kge1。
[0110]
步骤4.2迭代阶段；从n＝2开始，优化算法在参数范围内选定较上次更优的参数集，得到第n次模拟的径流序列{q
sim,n
(i)}，计算实测序列与模拟序列之间的目标函数值kgen，直至kgen值趋于收敛。此时，保留水文模型参数，用于流域径流模拟。
[0111]
步骤5、利用率定后的水文模型模拟流域径流过程，并评价水文模型的模拟精度。采用水文评价指标kge和r2评估耦合模型的径流模拟表现(附图2)，kge
→
1、r2→
1表示模型越好。其中，r2由下式所求：
[0112][0113]
式中：q
obs,i
为i时刻的实测径流量；q
sim,i
为第i时刻的模拟径流量；n指水文观测资料序列的长度。
[0114]
本实施例以全球各气候类型的2116个流域为研究对象，评估耦合后的水文模型的流域模拟径流精度，结果表明(见图3)，awbm、gr2m、simhyd、tank模型的中值kge分别从0.39到0.46，从0.71到0.72，从0.03到0.74，从0.63到0.65；中值r2分别从0.38到0.35，从0.65到0.68，从0.19到0.60，从0.47到0.49，由此可见，耦合后的水文模型的流域径流模拟方法提高了径流模拟的一致性。这是因为更具有物理意义的蒸发估计方法通过流域水量平衡为流域径流模拟增加了可靠约束，本发明的水文模型实现更优的径流模拟表现和更低的不确定性。
[0115]
以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，均应当包含在本发明的保护范围内。

技术特征：
1.一种流域径流的模拟方法，其特征在于，包括如下步骤：步骤1、获取研究流域的基础水文气象数据；步骤2、基于获取的基础水文气象数据和广义蒸发互补关系，构建蒸发估算；步骤3、将蒸发估算与集总式水文模型进行耦合获得耦合后的水文模型；步骤4、构建参数率定的目标函数，采用优化算法结合实测数据率定耦合后的水文模型的参数；步骤5、利用率定后的水文模型模拟流域径流过程，并评价水文模型的模拟精度。2.根据权利要求1所述的流域径流的模拟方法，其特征在于，基础水文气象数据包括月尺度的降水p、径流q、净辐射rn、近地面2m处风速u2、气温t、饱和水汽压差d。3.根据权利要求2所述的流域径流的模拟方法，其特征在于，净辐射rn根据数据集中采集的向下短波辐射r
si
、向下长波辐射r
li
、地表反照率α和气温t确定；饱和水汽压差根据数据集中采集的气温t和相对湿度rh确定；近地面2m处风速u2根据数据集中的水平风速u和垂向风速v确定。4.根据权利要求1所述的流域径流的模拟方法，其特征在于，步骤2中表观潜在蒸发e
pa
的计算方式如下：式中：δ指给定的空气温度下饱和蒸汽压力曲线的斜率；γ指干湿表常数；r
n
指陆地表面的净辐射；g指地面热通量,在日尺度及更长的时间尺度上可忽略不计；d指饱和水汽压差；f(u2)指风函数，其表达式如下表示：f(u2)＝0.26(1+0.54u2)；其中，u2表示近地面2m处风速，由步骤1获得。5.根据权利要求1所述的流域径流的模拟方法，其特征在于，步骤2中平衡蒸发e
e
的计算方式如下：式中：δ指给定的空气温度下饱和蒸汽压力曲线的斜率；γ指干湿表常数；r
n
指陆地表面的净辐射；g指地面热通量。6.根据权利要求1所述的流域径流的模拟方法，其特征在于，步骤2中构建的蒸发估算为：式中，α
c
为参数。7.根据权利要求1所述的流域径流的模拟方法，其特征在于，集总式水文模型包括awbm、gr2m、simhyd、tank模型中的一种。8.根据权利要求1所述的流域径流的模拟方法，其特征在于，步骤4中构建的目标函数为kge，其公式为：
式中：q
obs
指观测径流；q
sim
指模拟径流；分别指观测径流和模拟径流序列的平均值；σ
o
、σ
s
分别指观测径流和模拟径流序列的标准差；为i时刻的实测径流量。9.根据权利要求8所述的流域径流的模拟方法，其特征在于，步骤4中率定参数的方法为：步骤4.1初始化阶段；输入水文气象观测资料序列，指定模型参数可行范围，在可行范围内随机搜索确定一个参数集，可得出第1次模拟的径流序列{q
sim,1
(i)}，计算实测序列与模拟序列之间的目标函数值kge1；步骤4.2迭代阶段；从n＝2开始，优化算法在参数范围内选定较上次更优的参数集，得到第n次模拟的径流序列{q
sim,n
(i)}，计算实测序列与模拟序列之间的目标函数值kge
n
，直至kge
n
值趋于收敛；此时，保留水文模型参数，用于流域径流模拟。10.根据权利要求1所述的流域径流的模拟方法，其特征在于，步骤5中采用水文评价指标kge和r2评估率定后的水文模型的径流模拟精度。

技术总结
本发明公开了一种流域径流的模拟方法，包括如下步骤：获取研究流域的基础水文气象数据；基于获取的基础水文气象数据和广义蒸发互补关系，构建蒸发估算；将蒸发估算与集总式水文模型进行耦合获得耦合后的水文模型；构建参数率定的目标函数，采用优化算法结合实测数据率定耦合后的水文模型的参数；利用率定后的水文模型模拟流域径流过程，并评价水文模型的模拟精度。本发明解决了集总式水文模型过度简化蒸发过程的问题，提高了集总式水文模型的径流模拟能力和模型适用性。模拟能力和模型适用性。模拟能力和模型适用性。


技术研发人员：雷旭昕 程磊 张晓菁 张橹 刘攀
受保护的技术使用者：武汉大学
技术研发日：2023.05.23
技术公布日：2023/9/5