CN115422779B - 一种基于常微分方程组的新安江模型的构建方法及其水文预报方法 - Google Patents

Abstract

本发明公开了一种基于常微分方程组的新安江模型的构建方法及其水文预报方法，构建步骤包含：推导常微分方程形式的产流、分水源、坡面汇流、河网汇流方程；以向量形式组装推导的常微分方程，形成基于常微分方程组的新安江模型；确定基于常微分方程组的新安江模型的输入和参数，采用高阶精度的数值求解方法进行计算，得到流域实际蒸散发强度和出口断面流量过程。本发明推导了新安江模型连续形式的数学方程，引入了稳定且计算精度高的数值求解方法，有效地降低和控制了新安江模型的数值误差，解决了原有新安江模型存在较大的数值误差的问题，同时为模型的实现、在水文预报中的应用和进一步改进提供便利，具有较强的应用前景。

Description

一种基于常微分方程组的新安江模型的构建方法及其水文预报方法

技术领域

本发明属于水文模型技术领域，具体涉及一种基于常微分方程组的新安江模型的构建方法及其水文预报方法。

背景技术

作为重要的防洪非工程措施，水文预报得到广泛地重视和研究，是“预报、预警、预演、预案”四预体系的起点和重要组成内容。水文预报的核心是水文模型，按照发展阶段可将水文模型分为黑箱模型、集总式模型和分布式模型。由于模型结构、参数物理意义明确且对数据要求低，集总式模型是目前水文预报主要使用的模型。新安江模型是由河海大学赵人俊教授团队于1973年提出的集总式模型，是中国水文模型研究的代表成果并得到广泛的应用。

尽管新安江模型提出已接近50年，但其在计算过程中使用时段初值代替时段均值，存在较大的数值误差，由于模型使用的方程为直接建立在时段上的离散形式的计算方程，而非连续形式的数学方程，因此无法使用高阶精度、稳定的数值求解方法以减少数值误差，同时也给模型的实现和进一步改进带来不便。新安江模型的数值误差会随着降雨强度的增大而增大，在极端降水事件频发的背景下，会进一步限制模型的模拟精度和实际应用，是新安江模型研究的重要问题。

发明内容

基于以上技术问题，本发明的一个目的是提供一种基于常微分方程组的新安江模型的构建方法，包括如下步骤：

步骤1、产流模块构建：基于三层蒸散发模型和张力水蓄水容量曲线，推导常微分方程形式的产流方程，根据实测流域平均降雨强度、实测水面蒸发强度计算流域实际蒸散发强度、产流强度和产流面积比例；

步骤2、分水源模块构建：基于自由水蓄水容量曲线，推导常微分方程形式的分水源方程，根据步骤1得到的产流面积比例将步骤1得到的产流强度划分为地表产流强度、壤中流产流强度和地下产流强度；

步骤3、坡面汇流模块构建：基于线性水库，推导常微分方程形式的坡面汇流方程，根据步骤2得到的壤中流产流强度、地下产流强度计算由坡面进入河道的壤中流径流强度、地下径流强度，地表产流直接进入河道形成地表径流，地表径流强度、壤中流径流强度及地下径流强度之和为河网入流强度；

步骤4、河网汇流模块构建：基于Nash串联线性水库，推导常微分方程形式的河网汇流模块，根据步骤3得到的河网的入流强度计算流域出口流量；

步骤5、基于常微分方程组的新安江模型构建：以向量形式组装步骤1、2、3、4得到的常微分方程，形成以常微分方程组描述的新安江模型。

进一步的，所述步骤1中，推导的常微分方程形式的产流方程为：

式(1)、(2)和(3)中，W_u、W_l和W_d分别为上层、下层和深层土壤的流域平均张力水蓄量，t为时间，P_n为净雨强度，E_u、E_l和E_d分别为上层、下层和深层土壤的实际蒸散发强度，I_u为上层土壤蓄满后对下层土壤的补给强度，I_l为下层土壤蓄满后对深层土壤的补给强度，R为产流强度。P_n对应方程为：

P_n＝max(P-E_p,0)                        (4)

E_n＝max(E_p-P,0)                        (5)

式(4)和(5)中，max为取最大值函数，P为实测流域平均降雨强度，E_n为净蒸发强度，E_p为流域蒸散发能力。E_p对应方程为：

E_p＝K_e·E_obs                              (6)

式(6)中，K_e为蒸散发系数，E_obs为实测水面蒸发强度。E_u、E_l和E_d对应方程为：

式(7)、(8)和(9)中，min为取最小值函数，W_lm是下层土壤张力水容量(W_l的最大值)，c是深层土壤蒸散发系数。流域实际蒸散发强度E_t对应方程为：

E_t＝E_u+E_l+E_d                          (10)

I_u和I_l对应方程为：

式(11)中，W_um是上层土壤张力水容量(W_u的最大值)。R对应方程为：

R＝P_nf_w                                 (13)

式(13)中，f_w为流域产流面积比例，其对应方程为：

式(14)中，A_imp为流域不透水面积比例，b为张力水蓄水容量曲线指数，W₀为流域平均张力水蓄量，W_m为流域平均张力水容量(W₀的最大值)。W₀和W_m对应方程为：

W₀＝W_u+W_l+W_d                          (15)

W_m＝W_um+W_lm+W_dm                       (16)

式(16)中，W_dm是深层土壤张力水容量(W_d的最大值)。

进一步的，所述步骤2中，推导的常微分方程形式的分水源方程为：

式(17)中，S₀是流域平均自由水蓄量，R_pa是透水面积上的产流强度，R_s是透水面积上的地表产流强度，R_i是壤中流产流强度，R_g是地下产流强度。R_pa对应方程为：

R_pa＝P_n(f_w-A_imp)                          (18)

R_s、R_i和R_g对应方程为：

R_i＝K_i(f_w-A_imp)S₀                        (20)

R_g＝K_g(f_w-A_imp)S₀                        (21)

式(19)、(20)和(21)中，S_m是流域平均自由水容量(S₀的最大值)，ex是自由水蓄水容量曲线指数，K_i和K_g分别为壤中流和地下径流出流系数。地表产流强度R_t对应方程为：

式(22)中，R_ia为不透水面积上的地表产流强度。

进一步的，所述步骤3中，推导的常微分方程形式的坡面汇流方程为：

式(23)和(24)中，O_i和O_g分别为壤中流和地下径流线性水库的蓄量，Q_i和Q_g分别为壤中流径流强度和地下径流强度。Q_i和Q_g对应方程为：

Q_i＝-O_ilnC_i/Δt                          (25)

Q_g＝-O_glnC_g/Δt                          (26)

式(25)和(26)中，C_i和C_g分别为壤中流和地下径流的消退系数，Δt为时段长。河网入流强度Q_t的对应方程为：

Q_t＝R_t+Q_i+Q_g                            (27)

进一步的，所述步骤4中，推导的常微分方程形式的河网汇流方程为：

式(28)和(29)中，F_i是第i个线性水库的蓄量(i＝1,2,…,n)，n是Nash串联线性水库个数，Q_i是第i个线性水库的出流强度(i＝1,2,…,n)，其对应方程为：

Q_i＝F_i/K_f                              (30)

式(30)中，K_f是Nash串联线性水库出流系数。第n个线性水库的出流强度Q_n即为流域出口的出流强度，量纲为L/T，需要根据流域面积A换算到常用单位m³/s，进而得到流域出口的瞬时流量Q_instant，对于一个特定的时间t＝T，流域出口流量Q对应方程为：

进一步的，所述步骤5中，描述新安江模型的常微分方程组对应方程为：

本发明的另一个目的是提供一种采用所述的基于常微分方程组的新安江模型进行水文预报的方法，基于常微分方程组的新安江模型求解：给定模型输入即步骤1中所述的实测流域平均降雨强度、实测水面蒸发强度，确定模型参数，采用高阶精度的常微分方程数值求解方法对步骤5得到的常微分方程组进行求解，模型的输出为步骤1所述的流域实际蒸散发强度、步骤4所述的流域出口流量。

进一步的，基于常微分方程组的新安江模型包含2个固定参数和15个可调参数，固定参数为流域面积A和时段长Δt，流域面积根据数字高程模型(DEM)和流域出口站点经纬度可进行计算，时段长根据输入资料时间分辨率和模拟目的进行确定，15个可调参数的物理含义及取值范围为：

参数	物理含义	取值范围
			<![CDATA[K<sub>e</sub>]]>	蒸散发系数	[0.6,1.5]
c	深层土壤蒸散发系数	[0.05,2]
			<![CDATA[W<sub>um</sub>]]>	上层土壤张力水容量	[5,30]
<![CDATA[W<sub>lm</sub>]]>	下层土壤张力水容量	[60,90]
			<![CDATA[W<sub>dm</sub>]]>	深层土壤张力水容量	[15,60]
<![CDATA[A<sub>imp</sub>]]>	流域不透水面积比例	[0.01,0.2]
			b	张力水蓄水容量曲线指数	[0.1,0.4]
<![CDATA[S<sub>m</sub>]]>	流域平均自由水容量	[10,50]
			ex	自由水蓄水容量曲线指数	[1,1.5]
<![CDATA[K<sub>i</sub>]]>	壤中流出流系数	[0.1,0.55]
			<![CDATA[K<sub>g</sub>]]>	地下径流出流系数	<![CDATA[0.7-K<sub>i</sub>]]>
<![CDATA[C<sub>i</sub>]]>	壤中流消退系数	[0.5,0.9]
			<![CDATA[C<sub>g</sub>]]>	地下径流消退系数	[0.98,0.998]
n	Nash串联线性水库个数	[1,20]
			<![CDATA[K<sub>f</sub>]]>	Nash串联线性水库出流系数	[0.01,20]

本发明所达到的有益效果：本发明提供的基于常微分方程组的新安江模型，基于三层蒸散发模型、张力水蓄水容量曲线、自由水蓄水容量曲线、线性水库和Nash串联线性水库，推导常微分方程形式的产流、分水源、坡面汇流、河网汇流方程，构建了以常微分方程组描述的新安江模型，得到了新安江模型连续形式的数学方程，引入了高阶精度、稳定的数值求解方法进行求解，有效地降低和控制了新安江模型的数值误差，同时为模型的实现、在水文预报中的应用和进一步改进提供便利，具有较强的工程意义。

附图说明

图1是本发明基于常微分方程组的新安江模型的流程图；

图2是基于常微分方程组的新安江模型的模型结构示意图；

图3是某流域基于常微分方程组的新安江模型的模拟结果图。

具体实施方式

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

实施例1

如图1所示，一种基于常微分方程组的新安江模型的构建方法，包括以下步骤：

步骤1、基于三层蒸散发模型和张力水蓄水容量曲线的思想，推导常微分方程形式的产流方程，根据实测流域平均降雨强度、实测水面蒸发强度计算流域实际蒸散发强度、产流强度和产流面积比例，对应常微分方程形式的产流方程为：

式(1)、(2)和(3)中，W_u、W_l和W_d分别为上层、下层和深层土壤的流域平均张力水蓄量(mm)，t为时间(day)，P_n为净雨强度(mm/day)，E_u、E_l和E_d分别为上层、下层和深层土壤的实际蒸散发强度(mm/day)，I_u为上层土壤蓄满后对下层土壤的补给强度(mm/day)，I_l为下层土壤蓄满后对深层土壤的补给强度(mm/day)，R为产流强度(mm/day)。P_n对应方程为：

P_n＝max(P-E_p,0)                        (4)

E_n＝max(E_p-P,0)                       (5)

式(4)和(5)中，max为取最大值函数，P为实测流域平均降雨强度(mm/day)，E_n为净蒸发强度(mm/day)，E_p为流域蒸散发能力(mm/day)。E_p对应方程为：

E_p＝K_e·E_obs                             (6)

式(6)中，K_e为蒸散发系数，E_obs为实测水面蒸发强度(mm/day)。E_u、E_l和E_d对应方程为：

式(7)、(8)和(9)中，min为取最小值函数，W_lm是下层土壤张力水容量(W_l的最大值，mm)，c是深层土壤蒸散发系数。流域实际蒸散发强度E_t(mm/day)对应方程为：

E_t＝E_u+E_l+E_d                           (10)

I_u和I_l对应方程为：

式(11)中，W_um是上层土壤张力水容量(W_u的最大值，mm)。R对应方程为：

R＝P_nf_w                                 (13)

式(13)中，f_w为流域产流面积比例，其对应方程为：

式(14)中，A_imp为流域不透水面积比例，b为张力水蓄水容量曲线指数，W₀为流域平均张力水蓄量(mm)，W_m为流域平均张力水容量(W₀的最大值，mm)。W₀和W_m对应方程为：

W₀＝W_u+W_l+W_d                          (15)

W_m＝W_um+W_lm+W_dm                      (16)

式(16)中，W_dm是深层土壤张力水容量(W_d的最大值，mm)；

步骤2、基于自由水蓄水容量曲线思想，推导常微分方程形式的分水源方程，根据步骤1得到的产流面积比例将步骤1得到的产流强度划分为地表产流强度、壤中流产流强度和地下产流强度，推导的常微分方程形式的分水源方程为：

式(17)中，S₀是流域平均自由水蓄量(mm)，R_pa是透水面积上的产流强度(mm/day)，R_s是透水面积上的地表产流强度(mm/day)，R_i是壤中流产流强度(mm/day)，R_g是地下产流强度(mm/day)。R_pa对应方程为：

R_pa＝P_n(f_w-A_imp)                        (18)

R_s、R_i和R_g对应方程为：

R_i＝K_i(f_w-A_imp)S₀                       (20)

R_g＝K_g(f_w-A_imp)S₀                      (21)

式(19)、(20)和(21)中，S_m是流域平均自由水容量(S₀的最大值，mm)，ex是自由水蓄水容量曲线指数，K_i和K_g分别为壤中流和地下径流出流系数(day^-1)。地表产流强度R_t(mm/day)对应方程为：

式(22)中，R_ia为不透水面积上的地表产流强度(mm/day)；

步骤3、基于线性水库的思想，推导常微分方程形式的坡面汇流方程，根据步骤2得到的壤中流产流强度、地下产流强度计算由坡面进入河道的壤中流径流强度、地下径流强度，地表产流直接进入河道形成地表径流，地表径流强度、壤中流径流强度及地下径流强度之和为河网入流强度，推导的常微分方程形式的坡面汇流方程为：

式(23)和(24)中，O_i和O_g分别为壤中流和地下径流线性水库的蓄量(mm)，Q_i和Q_g分别为壤中流径流强度和地下径流强度(mm/day)。Q_i和Q_g对应方程为：

Q_i＝-O_ilnC_i/Δt                          (25)

Q_g＝-O_glnC_g/Δt                        (26)

式(25)和(26)中，C_i和C_g分别为壤中流和地下径流的消退系数，Δt为时段长(day)。河网入流强度Q_t(mm/day)的对应方程为：

Q_t＝R_t+Q_i+Q_g                          (27)

步骤4、基于Nash串联线性水库思想，推导常微分方程形式的河网汇流模块，根据步骤3得到的河网的入流强度计算流域出口流量，推导的常微分方程形式的河网汇流方程为：

式(28)和(29)中，F_i是第i个线性水库的蓄量(i＝1,2,…,n)，其单位为mm，n是Nash串联线性水库个数，Q_i是第i个线性水库的出流强度(i＝1,2,…,n)，其单位为mm/day，其对应方程为：

Q_i＝F_i/K_f                               (30)

式(30)中，K_f是Nash串联线性水库出流系数(day)。第n个线性水库的出流强度Q_n即为流域出口的出流强度，量纲为L/T，需要根据流域面积A换算到常用单位m³/s，进而得到流域出口的瞬时流量Q_instant(m³/s)，对于一个特定的时间t＝T，流域出口流量Q(m³/s)对应方程为：

步骤5、以向量形式组装步骤1、2、3、4得到的常微分方程，形成以常微分方程组描述的新安江模型，模型结构示意如图2所示：

图中W为单点张力水容量(mm)，W_mm为流域内最大单点张力水容量(mm)，a_w为W₀(W_u+W_l+W_d)对应的单点张力水蓄量(mm)，这三个符号参与公式(13)和(14)的推导。

S为单点自由水容量(mm)、S_mm为流域内最大单点自由水容量(mm)，a_s为S₀对应的单点自由水蓄量(mm)，这三个符号参与公式(18)和(19)的推导。

描述新安江模型的常微分方程组对应方程为：

实施例2

采用实施例1构建的模型进行水文预报：

对于一个具体流域，收集模拟期内流域内雨量站实测降雨资料、蒸发站实测水面蒸发资料，通过流域面平均降雨计算和时间尺度转换得到日尺度的实测流域平均降雨强度、实测水面蒸发强度资料序列，作为模型输入；常微分方程形式的新安江模型包含2个固定参数和15个可调参数，固定参数为流域面积A和时段长Δt，通过收集数字高程模型(DEM)和流域出口水文站点经纬度并通过流域提取计算得到流域面积为2690km²，考虑输入资料时间分辨率和建模目的，构建日模型，时段长为1day，15个可调参数通过人工调整得到，其具体取值为：

参数	取值
		<![CDATA[K<sub>e</sub>]]>	1.35
c	0.12
		<![CDATA[W<sub>um</sub>]]>	21.2
<![CDATA[W<sub>lm</sub>]]>	74.1
		<![CDATA[W<sub>dm</sub>]]>	22.1
<![CDATA[A<sub>imp</sub>]]>	0.15
		b	0.16
<![CDATA[S<sub>m</sub>]]>	37.3
		ex	1.36
<![CDATA[K<sub>i</sub>]]>	0.50
		<![CDATA[K<sub>g</sub>]]>	0.20
<![CDATA[C<sub>i</sub>]]>	0.53
		<![CDATA[C<sub>g</sub>]]>	0.992
n	3
		<![CDATA[K<sub>f</sub>]]>	0.27

模型采用四阶显式可变步长的Runge-Kutta方法进行数值求解，输出为日尺度的流域实际蒸散发强度和流域出口流量，模型模拟结果如图3所示。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (4)

1.一种基于常微分方程组的新安江模型的构建方法，其特征在于，包括如下步骤：

步骤1、产流模块构建：推导常微分方程形式的产流方程，根据实测流域平均降雨强度、实测水面蒸发强度计算流域实际蒸散发强度、产流强度和产流面积比例，产流面积比例对应方程为：

式(14)中，f_w为流域产流面积比例，A_imp为流域不透水面积比例，b为张力水蓄水容量曲线指数，W₀为流域平均张力水蓄量，W_m为流域平均张力水容量，为W₀的最大值；

步骤2、分水源模块构建：推导常微分方程形式的分水源方程，根据步骤1得到的产流面积比例将步骤1得到的产流强度划分为地表产流强度、壤中流产流强度和地下产流强度，地表产流强度对应方程为：

式(22)中，R_t为地表产流强度，R_ia为不透水面积上的地表产流强度，R_s是透水面积上的地表产流强度，P_n为净雨强度，S₀是流域平均自由水蓄量，S_m是流域平均自由水容量，为S₀的最大值，ex是自由水蓄水容量曲线指数；

步骤3、坡面汇流模块构建：推导常微分方程形式的坡面汇流方程，根据步骤2得到的壤中流产流强度、地下产流强度计算由坡面进入河道的壤中流径流强度、地下径流强度，地表产流直接进入河道形成地表径流，地表径流强度、壤中流径流强度及地下径流强度之和为河网入流强度；

步骤4、河网汇流模块构建：推导常微分方程形式的河网汇流模块，根据步骤3得到的河网的入流强度计算流域出口流量；

步骤5、基于常微分方程组的新安江模型构建：以向量形式组装步骤1、2、3、4得到的常微分方程，形成基于常微分方程组的新安江模型；

所述步骤1中，推导的常微分方程形式的产流方程为：

式(1)、(2)和(3)中，W_u、W_l和W_d分别为上层、下层和深层土壤的流域平均张力水蓄量，t为时间，E_u、E_l和E_d分别为上层、下层和深层土壤的实际蒸散发强度，I_u为上层土壤蓄满后对下层土壤的补给强度，I_l为下层土壤蓄满后对深层土壤的补给强度，R为产流强度；P_n对应方程为：

P_n＝max(P-E_p,0)                         (4)

E_n＝max(E_p-P,0)                        (5)

式(4)和(5)中，max为取最大值函数，P为实测流域平均降雨强度，E_n为净蒸发强度，E_p为流域蒸散发能力；E_p对应方程为：

E_p＝K_e·E_obs                              (6)

式(6)中，K_e为蒸散发系数，E_obs为实测水面蒸发强度；E_u、E_l和E_d对应方程为：

式(7)、(8)和(9)中，min为取最小值函数，W_lm是下层土壤张力水容量，为W_l的最大值，c是深层土壤蒸散发系数；流域实际蒸散发强度E_t对应方程为：

E_t＝E_u+E_l+E_d                           (10)

I_u和I_l对应方程为：

式(11)中，W_um是上层土壤张力水容量，为W_u的最大值；R对应方程为：

R＝P_nf_w                 (13)

W₀和W_m对应方程为：

W₀＝W_u+W_l+W_d                    (15)

W_m＝W_um+W_lm+W_dm                (16)

式(16)中，W_dm是深层土壤张力水容量，为W_d的最大值；

所述步骤2中，推导的常微分方程形式的分水源方程为：

式(17)中，R_pa是透水面积上的产流强度，R_s是透水面积上的地表产流强度，R_i是壤中流产流强度，R_g是地下产流强度；R_pa对应方程为：

R_pa＝P_n(f_w-A_imp)                   (18)

R_s、R_i和R_g对应方程为：

R_i＝K_i(f_w-A_imp)S₀                   (20)

R_g＝K_g(f_w-A_imp)S₀                    (21)

式(19)、(20)和(21)中，K_i和K_g分别为壤中流和地下径流出流系数；

所述步骤3中，推导的常微分方程形式的坡面汇流方程为：

式(23)和(24)中，O_i和O_g分别为壤中流和地下径流线性水库的蓄量，Q_i和Q_g分别为壤中流径流强度和地下径流强度；Q_i和Q_g对应方程为：

Q_i＝-O_i lnC_i/Δt                          (25)

Q_g＝-O_g lnC_g/Δt                         (26)

式(25)和(26)中，C_i和C_g分别为壤中流和地下径流的消退系数,Δt为时段长；

河网入流强度Q_t的对应方程为：

Q_t＝R_t+Q_i+Q_g                           (27)

所述步骤4中，推导的常微分方程形式的河网汇流方程为：

式(28)和(29)中，F_i是第i个线性水库的蓄量，i＝1,2,…,n，n是Nash串联线性水库个数，Q_i是第i个线性水库的出流强度，i＝1,2,…,n，其对应方程为：

Q_i＝F_i/K_f                               (30)

式(30)中，K_f是Nash串联线性水库出流系数；第n个线性水库的出流强度Q_n即为流域出口的出流强度，量纲为L/T，需要根据流域面积A换算到常用单位m³/s，进而得到流域出口的瞬时流量Q_instant，对于一个特定的时间t＝T，流域出口流量Q对应方程为：

2.根据权利要求1所述的基于常微分方程组的新安江模型的构建方法，其特征在于：所述步骤5中，描述新安江模型的常微分方程组对应方程为：

3.一种采用权利要求1或2所述的基于常微分方程组的新安江模型进行水文预报的方法，其特征在于：基于常微分方程组的新安江模型求解：给定模型输入即步骤1中所述的实测流域平均降雨强度、实测水面蒸发强度，确定模型参数，采用高阶精度的常微分方程数值求解方法对步骤5得到的常微分方程组进行求解，模型的输出为步骤1所述的流域实际蒸散发强度、步骤4所述的流域出口流量。

4.根据权利要求3所述的水文预报方法，其特征在于：基于常微分方程组的新安江模型包含2个固定参数和15个可调参数，固定参数为流域面积A和时段长Δt，流域面积根据数字高程模型和流域出口站点经纬度通过计算获得，时段长根据输入资料时间分辨率和建模目的进行确定，15个固定参数均有物理含义及取值范围。

Images