CN114429089A - 一种岩溶区分布式非线性水文模拟方法 - Google Patents

Abstract

本发明提供一种岩溶区分布式非线性水文模拟方法，首先获取流域数字高程模型数据并在此基础上计算栅格的地形指数，获取流域岩溶发育地质图和石漠化地质图，判断流域内各栅格的岩溶发育和石漠化情况；然后构建表层岩溶带自由水蓄水库；之后计算表层岩溶带自由水蓄水容量，表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数，表层岩溶带自由水蓄水库对地下水的日出流系数；最后通过栅格内坡面汇流和栅格间马斯京根法汇流将各栅格的产流汇到流域出口。本发明构建表层岩溶带自由水蓄水库和地下溶隙线性水库计算岩溶流域的径流，模拟岩溶流域表层岩溶带和地下溶隙对径流的调蓄作用，为岩溶流域水文模拟提供方法，促进分布式水文模型的发展。

Description

一种岩溶区分布式非线性水文模拟方法

技术领域

本发明属于水文模拟与预报领域，具体涉及一种岩溶区分布式非线性水文模拟方法。

背景技术

全球岩溶面积约为220万km²，居住人口约10亿人。中国西南部是世界上岩溶发育最典型的地区之一，具有集中连片的岩溶地貌，分布面积约5.4×10⁵km²。由于石漠化、植被稀少等原因，岩溶流域的土层厚度、土地覆盖类型等下垫面条件空间异质性大。流域内岩溶裂隙发育程度不同，导致了流域产汇流机制复杂，洪水模拟难度大。岩溶流域中的表层岩溶带以及地下深层溶隙、溶洞等都是地表水入渗、地下水赋存和运移的良好介质，在产汇流过程中起重要作用，探究地表水-岩溶水-土壤水的转化和产汇流规律十分必要。岩溶流域下垫面空间异质性较大，其洪水模拟难度大。集总式的岩溶水文模型难以刻划流域下垫面条件的复杂性；水动力岩溶模型在实验小流域模拟效果好，但是需要精细的数据资料，且对于面积更大的中小河流流域，模型的洪水模拟计算量庞大。分布式和半分布式概念性岩溶模型，可以在考虑下垫面异质性的同时简化计算，但是目前分布式岩溶模型大部分被用于日径流和月径流的模拟，用于场次洪水模拟的分布式岩溶模型较少。需要设计一种用于岩溶流域面积较大的模拟场次洪水的水文模拟方法。

发明内容

本发明的目的是提供一种岩溶区分布式非线性水文模拟方法，构建表层岩溶带自由水蓄水库模拟表层岩溶带的调蓄作用，在坡面汇流的过程中构建消退系数不同的地下径流线性水库模拟岩溶发育区域的地下径流，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明具体采用以下技术方案：

一种岩溶区分布式非线性水文模拟方法，针对包含岩溶区的目标流域，执行以下步骤，实现目标流域出口断面流量的预测：

步骤1，基于目标流域数字高程模型，计算目标流域内各栅格的地形指数；

基于目标流域岩溶发育地质图和目标流域石漠化地质图，得到目标流域内各栅格的岩溶发育情况和石漠化情况；

步骤2，分别针对目标流域内各栅格，基于栅格新安江模型的模拟方法，得到各栅格土壤自由水蓄水库各径流的径流深；分别针对目标流域岩溶区的各岩溶栅格，结合各岩溶栅格的石漠化情况，构建各岩溶栅格的表层岩溶带自由水蓄水库；

步骤3，针对各岩溶栅格构建的表层岩溶带自由水蓄水库，结合目标流域内各栅格的地形指数，计算得到各岩溶栅格表层岩溶带自由水蓄水库容量，以及得到各岩溶栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数、表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数，进而得到各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深；

步骤4，针对目标流域，基于各栅格土壤自由水蓄水库各径流的径流深、各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深、以及各栅格的岩溶发育情况，通过各栅格内各径流的坡面汇流量，结合栅格间使用马斯京根法得到目标流域出口断面各径流的流量，进而实现目标流域出口断面流量的预测。

作为本发明的一种优选技术方案，所述步骤1包括：

步骤1-1，基于目标流域数字高程模型计算目标流域内各栅格的地形指数：

式中，Slope为栅格坡度，TI为栅格地形指数，

为栅格在经度方向的高程变化率，

为栅格在纬度方向的高程变化率，a为栅格单宽累积面积；

步骤1-2，基于目标流域石漠化地质图，将目标流域内的各栅格划分为3种类型：(1)非岩溶栅格，(2)未石漠化岩溶栅格，(3)石漠化岩溶栅格。

作为本发明的一种优选技术方案，所述步骤2包括：

步骤2-1，目标流域内的各栅格基于栅格新安江模型的模拟方法，在降水经过植被截流、蒸发、满足张力水蓄水容量后，获得的剩余雨量R称为栅格总径流深：

R＝P-I_CA-E-W_M (3)

获得的剩余雨量流入该栅格的土壤自由水蓄水库；式中：R为栅格总径流深，P为栅格降雨量，I_CA为栅格植被截流量，E为栅格蒸发量，W_M为栅格张力水蓄水容量；

步骤2-2，基于栅格新安江模型的模拟方法，分别对目标流域内的各非岩溶栅格，采用土壤自由水蓄水库进行分水源调蓄，定义土壤自由水蓄水库蓄水容量为S_M，判断总径流深R与当前时段初始土壤自由水蓄水库蓄水量S₀之和与S_M的关系，其中，R+S₀大于S_M产生地表径流，R+S₀小于S_M不产生地表径流；土壤自由水蓄水库有两个出口，一个旁侧出口形成壤中流，一个向下出口形成地下径流：

S_0t+1＝S-K_I×S-K_G×S (7)

式中：R为栅格总径流深，R_S为栅格地表径流深，R_I为栅格壤中流径流深，R_G为非岩溶栅格地下径流深，S_M为栅格土壤自由水蓄水库蓄水容量，K_I为栅格土壤自由水蓄水库对壤中流的日出流系数；K_G为栅格土壤自由水蓄水库对地下径流的日出流系数，S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，S为栅格当前时段土壤自由水蓄水库蓄水量，S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量；

步骤2-3，对于各未石漠化岩溶栅格，在土壤自由水蓄水库的下方构建表层岩溶带自由水蓄水库；上层土壤自由水蓄水库与下层表层岩溶带自由水蓄水库串联，定义表层岩溶带自由水蓄水库容量为S_MEK；

通过格林-安普特公式计算上层土壤自由水蓄水库的自由水向下层表层岩溶带自由水蓄水库渗透的过程：

式中：f(t)为栅格土壤自由水库下渗能力，K为栅格饱和水力传导度，Ψ为栅格湿润锋处土壤吸力，Δθ为栅格土壤饱和含水率与初始含水率之差，F(t)为栅格累计渗漏量，t为当前时段；

实际下渗量计算公式为：

f_a＝f(t)×Δt (9)

式中：f_a为实际下渗量，f(t)为栅格土壤自由水库下渗能力，Δt为时间步长；

基于栅格新安江模型的模拟方法，分别对于各未石漠化岩溶栅格土壤自由水蓄水库分水源调蓄的计算过程如下：

R_I＝K_I×S (12)

S_0t+1＝S-K_I×S-f_a (13)

式中：R为栅格总径流深，R_S为栅格地表径流深，R_I为栅格壤中流径流深，S_M为栅格土壤自由水蓄水库蓄水容量，S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，S为栅格当前时段土壤自由水蓄水库蓄水量，f_a为栅格实际下渗量，S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量；

下层表层岩溶带自由水蓄水库有两个出口，一个旁侧出口形成表层岩溶带侧向径流，一个向下出口形成岩溶地下径流：

S_EK＝S_EK0+f_a (14)

S_EK0t+1＝S_EK-K_EKI×S_EK-K_EKG×S_EK (16)

式中：S_EK0为栅格当前时段初始表层岩溶带自由水蓄水库蓄水量，S_EK为栅格当前时段表层岩溶带自由水蓄水库蓄水量，S_EK0t+1为栅格下一时段初始表层岩溶带自由水蓄水库蓄水量，R_EKI为栅格表层岩溶带侧向径流深，R_GK为栅格岩溶地下径流深，K_EKI为栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数；K_EKG为栅格表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数；

步骤2-4，分别对于各石漠化岩溶栅格，该栅格的土壤自由水蓄水库只有一个向下的出口，构建的表层岩溶带自由水蓄水库与未石漠化岩溶栅格相同。

作为本发明的一种优选技术方案，所述步骤3包括：

步骤3-1，针对各岩溶栅格，建立各岩溶栅格表层岩溶带自由水蓄水库容量S_MEK与各岩溶栅格地形指数的线性关系，假设各岩溶栅格中地形指数最大的栅格对应的表层岩溶带自由水蓄水库容量最大，地形指数最小的栅格对应的表层岩溶带自由水蓄水库容量最小：

式中：TI_max为岩溶区域内各栅格地形指数中的最大值；TI_min为岩溶区域内各栅格地形指数中的最小值；S_MEKmax为表层岩溶带自由水蓄水库最大值，S_MEKmin为表层岩溶带自由水蓄水库最小值，ξ_a与ξ_b为线性关系的两个系数；

通过参数率定确定S_MEKmax和S_MEKmin，从而得到ξ_a与ξ_b两个系数的值；

进而基于步骤1得到的每个岩溶栅格的地形指数经下式，得到每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK：

ξ_a×TI+ξ_b＝S_MEK (18)

式中：TI为各栅格的地形指数；S_MEK为岩溶区域内各栅格表层岩溶带蓄水容量，ξ_a与ξ_b为线性关系的两个系数；

步骤3-2，假设K_EKI与K_EKG之和与表层岩溶带自由水蓄水库容量S_MEK呈线性关系，表层岩溶带蓄水容量最大的栅格对应的K_EKI与K_EKG之和最小，表层岩溶带蓄水容量最小的栅格对应的K_EKI与K_EKG之和最大；

通过参数率定确定K_EKI与K_EKG之和的最大值与最小值，从而得到a与b两个系数的值；

进而基于公式(18)得到的每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK，经以下公式确定每个岩溶栅格K_EKI与K_EKG之和的值：

a×S_MEK+b＝(K_EKI+K_EKG) (20)

式中：S_MEK为岩溶区域内各栅格表层岩溶带自由水蓄水库容量，a与b为线性关系的两个系数；S_MEKmax、S_MEKmin为流域内各栅格表层岩溶带蓄水容量的最大值和最小值；

步骤3-3，假设K_EKI与K_EKG的比例为一个预设常数，通过参数率定确定该常数，即可确定各岩溶栅格的K_EKI和K_EKG值；

步骤3-4，基于各岩溶栅格对应的K_EKI、K_EKG值、以及表层岩溶带自由水蓄水库容量S_MEK，结合步骤2各岩溶栅格构建的表层岩溶带自由水蓄水库，得到各岩溶栅格的表层岩溶带侧向径流深R_EKI、岩溶地下径流深R_GK。

作为本发明的一种优选技术方案，所述步骤4包括：

步骤4-1，针对目标流域内各栅格采用线性水库法计算各径流坡面汇流量，

对于各非岩溶栅格：

式中：Q_S为栅格地表径流流量，Q_I为栅格壤中流流量，Q_G为非岩溶栅格地下径流流量，C_S为栅格地面径流消退系数；C_I为栅格壤中流消退系数；C_G为非岩溶栅格地下径流消退系数；t为当前时段，t-1为上一时段，U为单位换算系数，U＝栅格面积(km²)/3.6Δt，Δt为计算时间步长；

对于各岩溶栅格：

未石漠化岩溶栅格

石漠化岩溶栅格

式中：Q_EKI为表层岩溶带侧向出流流量，Q_GK为岩溶栅格地下径流流量，C_EKI为表层岩溶带侧向出流消退系数；C_GK为岩溶区地下径流消退系数；

步骤4-2，栅格间使用马斯京根法汇流将各栅格的流量向目标流域出口演算，第i个栅格流向第i+1个栅格的马斯京根法计算为：

Q(t)_i+1＝C₀Q(t)_i+C₁Q(t-1)_i+C₂Q(t-1)_i+1 (24)

式中：

X_e和K_e为马斯京根法的两个参数，Δt为计算时间步长，t为当前时段，t-1为上一时段；

目标流域内各栅格对应的Q_S、Q_I、Q_G、Q_EKI分别通过栅格间马斯京根法向目标流域出口演算，其中如果是岩溶栅格流向非岩溶栅格，则Q_EKI和Q_GK作为Q_G进行演算，最终得到目标流域出口断面的流量。

作为本发明的一种优选技术方案，所述步骤4-1中针对岩溶区地下径流消退系数C_GK，结合目标流域各栅格的岩溶发育情况，建立多个地下径流消退系数C_GK{C_GK1、C_GK2、…、C_GKn}，各地下径流消退系数分别与各岩溶发育程度一一对应，不同岩溶发育程度的栅格取不同的地下径流消退系数C_GKn。

本发明的有益效果是：本发明提供一种岩溶区分布式非线性水文模拟方法，首先获取流域数字高程模型数据并在此基础上计算栅格的地形指数，获取流域岩溶发育地质图和石漠化地质图，判断流域内各栅格的岩溶发育和石漠化情况；然后构建表层岩溶带自由水蓄水库；之后计算表层岩溶带自由水蓄水容量，表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数，表层岩溶带自由水蓄水库对地下水的日出流系数；最后通过栅格内坡面汇流和栅格间马斯京根法汇流将各栅格的产流汇到流域出口。本发明通过构建表层岩溶带自由水蓄水库和地下溶隙线性水库计算岩溶流域的径流，模拟岩溶流域表层岩溶带和地下溶隙对径流的调蓄作用，可以从空间上反映岩溶流域表层岩溶带和地下溶隙对产汇流过程的影响，且计算效率较高，可以适用于较大流域的水文模拟，为岩溶流域水文模拟提供方法，促进分布式水文模型的发展。

附图说明

图1是本发明提供的一种岩溶区分布式非线性水文模拟方法流程示意图；

图2是具体实施例中普厅河流域流域下垫面的TI值分布图、岩溶发育图和石漠化情况图；

图3是具体实施例中表层岩溶带自由水蓄水库结构图；

图4是具体实施例中普厅河流域表层岩溶带自由水蓄水容量、表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数、表层岩溶带自由水蓄水库对地下水的日出流系数空间分布图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步描述。

应当理解的是，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

本发明的目的是提供一种岩溶区分布式非线性水文模拟方法，构建表层岩溶带自由水蓄水库模拟表层岩溶带的调蓄作用，在坡面汇流的过程中构建消退系数不同的地下径流线性水库模拟岩溶发育区域的地下径流，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明具体采用以下技术方案：

如图1所示，一种岩溶区分布式非线性水文模拟方法，针对包含岩溶区的目标流域，执行以下步骤，实现目标流域出口断面流量的预测：

步骤1，基于目标流域数字高程模型(DEM)，计算目标流域内各栅格的地形指数；

基于目标流域岩溶发育地质图和目标流域石漠化地质图，得到目标流域内各栅格的岩溶发育情况和石漠化情况；

具体步骤包括：

步骤1-1，基于目标流域数字高程模型计算目标流域内各栅格的地形指数：

式中，Slope为栅格坡度，TI为栅格地形指数，

为栅格在经度方向的高程变化率，

为栅格在纬度方向的高程变化率，a为栅格单宽累积面积；

步骤1-2，基于目标流域石漠化地质图，将目标流域内的各栅格划分为3种类型：(1)非岩溶栅格，(2)未石漠化岩溶栅格，(3)石漠化岩溶栅格。

如图2所示为云南普厅河流域地形指数分布图，普厅河流域岩溶发育情况图和石漠化情况图。

步骤2，分别针对目标流域内各栅格，基于栅格新安江模型的模拟方法，得到各栅格土壤自由水蓄水库各径流的径流深；分别针对目标流域岩溶区的各岩溶栅格，结合各岩溶栅格的石漠化情况，构建各岩溶栅格的表层岩溶带自由水蓄水库；

具体步骤包括：

步骤2-1，目标流域内的各栅格基于栅格新安江模型的模拟方法，在降水经过植被截流、蒸发、满足张力水蓄水容量后，获得的剩余雨量R称为栅格总径流深：

R＝P-I_CA-E-W_M (3)

获得的剩余雨量流入该栅格的土壤自由水蓄水库；式中：R为栅格总径流深，mm；P为栅格降雨量，mm；I_CA为栅格植被截流量，mm；E为栅格蒸发量，mm；W_M为栅格张力水蓄水容量，mm；

步骤2-2，基于栅格新安江模型的模拟方法，分别对目标流域内的各非岩溶栅格，采用土壤自由水蓄水库进行分水源调蓄，定义土壤自由水蓄水库蓄水容量为S_M，判断总径流深R与当前时段初始土壤自由水蓄水库蓄水量S₀之和与S_M的关系，其中，R+S₀大于S_M产生地表径流，R+S₀小于S_M不产生地表径流；土壤自由水蓄水库有两个出口，一个旁侧出口形成壤中流，一个向下出口形成地下径流：

S_0t+1＝S-K_I×S-K_G×S (7)

式中：R为栅格总径流深，mm；R_S为栅格地表径流深，mm；R_I为栅格壤中流径流深，mm；R_G为栅格地下径流深，mm；S_M为栅格土壤自由水蓄水库蓄水容量，mm；K_I为栅格土壤自由水蓄水库对壤中流的日出流系数；K_G为栅格土壤自由水蓄水库对地下径流的日出流系数；S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，mm；S为栅格当前时段土壤自由水蓄水库蓄水量，mm；S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量，mm；

步骤2-3，对于各未石漠化岩溶栅格，在土壤自由水蓄水库的下方构建表层岩溶带自由水蓄水库；上层土壤自由水蓄水库与下层表层岩溶带自由水蓄水库串联，定义表层岩溶带自由水蓄水库容量为S_MEK；

通过格林-安普特公式计算上层土壤自由水蓄水库的自由水向下层表层岩溶带自由水蓄水库渗透的过程，使用格林-安普特下渗公式替代原始土壤自由水蓄水库的向下出流计算公式：

式中：f(t)为栅格土壤自由水库下渗能力，mm/h；K为栅格饱和水力传导度，mm/h；Ψ为栅格湿润锋处土壤吸力，mm；Δθ为栅格土壤饱和含水率与初始含水率之差，％；F(t)为栅格累计渗漏量，mm；t为当前时段，h；

实际下渗量计算公式为：

f_a＝f(t)×Δt (9)

式中：f_a为实际下渗量，mm；f(t)为栅格土壤自由水库下渗能力，mm/h；Δt为时间步长，h；

基于栅格新安江模型的模拟方法，分别对于各未石漠化岩溶栅格土壤自由水蓄水库分水源调蓄的计算过程如下：

R_I＝K_I×S (12)

S_0t+1＝S-K_I×S-f_a (13)

式中：R为栅格总径流深，mm；R_S为栅格地表径流深，mm；R_I为栅格壤中流径流深，mm；S_M为栅格土壤自由水蓄水库蓄水容量，mm；S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，mm；S为栅格当前时段土壤自由水蓄水库蓄水量，mm；f_a为栅格实际下渗量，mm；S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量，mm；

下层表层岩溶带自由水蓄水库有两个出口，一个旁侧出口形成表层岩溶带侧向径流，一个向下出口形成岩溶地下径流：

S_EK＝S_EK0+f_a (14)

S_EK0t+1＝S_EK-K_EKI×S_EK-K_EKG×S_EK (16)

式中：S_EK0为栅格当前时段初始表层岩溶带自由水蓄水库蓄水量，mm；S_EK为栅格当前时段表层岩溶带自由水蓄水库蓄水量，mm；S_EK0t+1为栅格下一时段初始表层岩溶带自由水蓄水库蓄水量，mm；R_EKI为栅格表层岩溶带侧向径流深，mm；R_GK为栅格岩溶地下径流深，mm；K_EKI为栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数；K_EKG为栅格表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数；

步骤2-4，分别对于各石漠化岩溶栅格，对于各石漠化岩溶栅格，由于土壤贫瘠，认为不存在壤中流，认为该栅格的土壤自由水蓄水库只有一个向下的出口，构建的表层岩溶带自由水蓄水库与未石漠化岩溶栅格相同，其余计算与未石漠化岩溶栅格相同。

则针对各石漠化岩溶栅格，基于栅格新安江模型的模拟方法，分别对于各石漠化岩溶栅格土壤自由水蓄水库分水源调蓄的计算过程如下：

S_0t+1＝S-f_a

各石漠化岩溶栅格下层表层岩溶带自由水蓄水库有两个出口，一个旁侧出口形成表层岩溶带侧向径流，一个向下出口形成岩溶地下径流：

S_EK＝S_EK0+f_a

S_EK0t+1＝S_EK-K_EKI×S_EK-K_EKG×S_EK

如图3所示为表层岩溶带自由水蓄水库结构图。

步骤3，针对各岩溶栅格构建的表层岩溶带自由水蓄水库，结合目标流域内各栅格的地形指数，计算得到各岩溶栅格表层岩溶带自由水蓄水库容量，以及得到各岩溶栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数、表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数，进而得到各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深；

具体步骤包括：

步骤3-1，针对各岩溶栅格，建立各岩溶栅格表层岩溶带自由水蓄水库容量S_MEK与各岩溶栅格地形指数的线性关系，假设各岩溶栅格中地形指数最大的栅格对应的表层岩溶带自由水蓄水库容量最大，地形指数最小的栅格对应的表层岩溶带自由水蓄水库容量最小：

式中：TI_max为岩溶区域内各栅格地形指数中的最大值；TI_min为岩溶区域内各栅格地形指数中的最小值；S_MEKmax为表层岩溶带自由水蓄水库最大值，mm；S_MEKmin为表层岩溶带自由水蓄水库最小值，mm；ξ_a与ξ_b为线性关系的两个系数；

通过参数率定确定S_MEKmax和S_MEKmin，从而得到ξ_a与ξ_b两个系数的值；

进而基于步骤1得到的每个岩溶栅格的地形指数经下式，得到每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK：

ξ_a×TI+ξ_b＝S_MEK (18)

式中：TI为各栅格的地形指数；S_MEK为岩溶区域内各栅格表层岩溶带蓄水容量，mm；ξ_a与ξ_b为线性关系的两个系数；最终计算的普厅河流域表层岩溶带自由水蓄水容量S_MEK如图4。

步骤3-2，假设K_EKI与K_EKG之和与表层岩溶带自由水蓄水库容量S_MEK呈线性关系，由于表层岩溶带越深渗透性越弱，所以表层岩溶带浅的区域整体上渗透性大，K_EKI与K_EKG之和越大渗透性越高，故假设表层岩溶带蓄水容量最大的栅格对应的K_EKI与K_EKG之和最小，表层岩溶带蓄水容量最小的栅格对应的K_EKI与K_EKG之和最大。通过率定K_EKI与K_EKG之和的最大值与最小值，即可确定每个岩溶栅格K_EKI与K_EKG之和的值。表层岩溶带蓄水容量最大的栅格对应的K_EKI与K_EKG之和最小，表层岩溶带蓄水容量最小的栅格对应的K_EKI与K_EKG之和最大；

通过参数率定确定K_EKI与K_EKG之和的最大值与最小值，从而得到a与b两个系数的值；

进而基于公式(18)得到的每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK，经以下公式确定每个岩溶栅格K_EKI与K_EKG之和的值：

a×S_MEK+b＝(K_EKI+K_EKG) (20)

式中：S_MEK为岩溶区域内各栅格表层岩溶带自由水蓄水库容量，mm；a与b为线性关系的两个系数；S_MEKmax、S_MEKmin为流域内各栅格表层岩溶带蓄水容量的最大值和最小值，mm；普厅河流域3个参数S_MEK、K_EKI、K_EKG计算结果的空间分布如图4所示。

步骤3-3，假设K_EKI与K_EKG的比例为一个预设常数，通过参数率定确定该常数，即可确定各岩溶栅格的K_EKI和K_EKG值；

步骤3-4，基于各岩溶栅格对应的K_EKI、K_EKG值、以及表层岩溶带自由水蓄水库容量S_MEK，结合步骤2各岩溶栅格构建的表层岩溶带自由水蓄水库，得到各岩溶栅格的表层岩溶带侧向径流深R_EKI、岩溶地下径流深R_GK。

步骤4，针对目标流域，基于各栅格土壤自由水蓄水库各径流的径流深、各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深、以及各栅格的岩溶发育情况，通过各栅格内各径流的坡面汇流量，结合栅格间使用马斯京根法得到目标流域出口断面各径流的流量，进而实现目标流域出口断面流量的预测；

具体步骤包括：

步骤4-1，针对目标流域内各栅格采用线性水库法计算各径流坡面汇流量，

对于各非岩溶栅格：

式中：Q_S为栅格地表径流流量，m³/s；Q_I为栅格壤中流流量，m³/s；Q_G为栅格地下径流流量，m³/s；C_S为栅格地面径流消退系数；C_I为栅格壤中流消退系数；C_G为栅格地下径流消退系数；t为当前时段，h；t-1为上一时段，h；U为单位换算系数，U＝栅格面积(km²)/3.6Δt，Δt为计算时间步长，h；

对于各岩溶栅格：

未石漠化岩溶栅格

石漠化岩溶栅格

式中：Q_EKI为表层岩溶带侧向出流流量，m³/s；Q_GK为岩溶栅格地下径流流量，m³/s；C_EKI为表层岩溶带侧向出流消退系数；C_GK为岩溶区地下径流消退系数；

针对岩溶区地下径流消退系数C_GK，结合目标流域各栅格的岩溶发育情况，建立多个地下径流消退系数C_GK{C_GK1、C_GK2、…、C_GKn}，各地下径流消退系数分别与各岩溶发育程度一一对应，不同岩溶发育程度的栅格取不同的地下径流消退系数C_GKn。

如图2所示，普厅河流域的有两种岩溶发育情况，“岩溶发育区”和“岩溶洞穴区”，所以在普厅河流域选取C_GK1、C_GK2两个参数。C_GK1用于一般的“岩溶发育区”，C_GK2用于“岩溶洞穴区”，C_GKn值越小汇流时间越短，“岩溶洞穴区”汇流时间短，所以认为C_GK1>C_GK2。

步骤4-2，栅格间使用马斯京根法汇流将各栅格的流量向目标流域出口演算，第i个栅格流向第i+1个栅格的马斯京根法计算为：

Q(t)_i+1＝C₀Q(t)_i+C₁Q(t-1)_i+C₂Q(t-1)_i+1 (24)

式中：

X_e和K_e为马斯京根法的两个参数，Δt为计算时间步长，h；t为当前时段，h；t-1为上一时段，h；

目标流域内各栅格对应的Q_S、Q_I、Q_G、Q_EKI分别通过栅格间马斯京根法向目标流域出口演算，其中如果是岩溶栅格流向非岩溶栅格，则Q_EKI和Q_GK作为Q_G进行演算，目标流域内各栅格向目标流域出口通过栅格间马斯京根法运算最终得到的Q_S、Q_I、Q_G、Q_EKI为各径流的流量，各径流流量的总和为最终目标流域出口断面的流量。

上述技术方案所设计提供一种岩溶区分布式非线性水文模拟方法，首先获取流域数字高程模型数据并在此基础上计算栅格的地形指数，获取流域岩溶发育地质图和石漠化地质图，判断流域内各栅格的岩溶发育和石漠化情况；然后构建表层岩溶带自由水蓄水库；之后计算表层岩溶带自由水蓄水容量，表层岩溶带自由水蓄水库对表层岩溶带侧向出流的日出流系数，表层岩溶带自由水蓄水库对地下水的日出流系数；最后通过栅格内坡面汇流和栅格间马斯京根法汇流将各栅格的产流汇到流域出口。本发明通过构建表层岩溶带自由水蓄水库和地下溶隙线性水库计算岩溶流域的径流，可以从空间上反映岩溶流域表层岩溶带和地下溶隙对产汇流过程的影响，且计算效率较高，可以适用于较大流域的水文模拟。

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (6)

1.一种岩溶区分布式非线性水文模拟方法，其特征在于：针对包含岩溶区的目标流域，执行以下步骤，实现目标流域出口断面流量的预测：

步骤1，基于目标流域数字高程模型，计算目标流域内各栅格的地形指数；

基于目标流域岩溶发育地质图和目标流域石漠化地质图，得到目标流域内各栅格的岩溶发育情况和石漠化情况；

步骤2，分别针对目标流域内各栅格，基于栅格新安江模型的模拟方法，得到各栅格土壤自由水蓄水库各径流的径流深；分别针对目标流域岩溶区的各岩溶栅格，结合各岩溶栅格的石漠化情况，构建各岩溶栅格的表层岩溶带自由水蓄水库；

步骤3，针对各岩溶栅格构建的表层岩溶带自由水蓄水库，结合目标流域内各栅格的地形指数，计算得到各岩溶栅格表层岩溶带自由水蓄水库容量，以及得到各岩溶栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数、表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数，进而得到各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深；

步骤4，针对目标流域，基于各栅格土壤自由水蓄水库各径流的径流深、各岩溶栅格的表层岩溶带自由水蓄水库各径流的径流深、以及各栅格的岩溶发育情况，通过各栅格内各径流的坡面汇流量，结合栅格间使用马斯京根法得到目标流域出口断面各径流的流量，进而实现目标流域出口断面流量的预测。

2.根据权利要求1所述的一种岩溶区分布式非线性水文模拟方法，其特征在于，所述步骤1包括：

步骤1-1，基于目标流域数字高程模型计算目标流域内各栅格的地形指数：

式中，Slope为栅格坡度，TI为栅格地形指数，

为栅格在经度方向的高程变化率，

为栅格在纬度方向的高程变化率，a为栅格单宽累积面积；

步骤1-2，基于目标流域石漠化地质图，将目标流域内的各栅格划分为3种类型：(1)非岩溶栅格，(2)未石漠化岩溶栅格，(3)石漠化岩溶栅格。

3.根据权利要求2所述的一种岩溶区分布式非线性水文模拟方法，其特征在于，所述步骤2包括：

步骤2-1，目标流域内的各栅格基于栅格新安江模型的模拟方法，在降水经过植被截流、蒸发、满足土壤张力水蓄水容量后，获得的剩余雨量R称为栅格总径流深：

R＝P-I_CA-E-W_M (3)

获得的剩余雨量流入该栅格的土壤自由水蓄水库；式中：R为栅格总径流深，P为栅格降雨量，I_CA为栅格植被截流量，E为栅格蒸发量，W_M为栅格张力水蓄水容量；

步骤2-2，基于栅格新安江模型的模拟方法，分别对目标流域内的各非岩溶栅格，采用土壤自由水蓄水库进行分水源调蓄，定义土壤自由水蓄水库蓄水容量为S_M，判断总径流深R与当前时段初始土壤自由水蓄水库蓄水量S₀之和与S_M的关系，其中，R+S₀大于S_M产生地表径流，R+S₀小于S_M不产生地表径流；

土壤自由水蓄水库有两个出口，一个旁侧出口形成壤中流，一个向下出口形成地下径流：

S_0t+1＝S-K_I×S-K_G×S (7)

式中：R为栅格总径流深，R_S为栅格地表径流深，R_I为栅格壤中流径流深，R_G为非岩溶栅格地下径流深，S_M为栅格土壤自由水蓄水库蓄水容量，K_I为栅格土壤自由水蓄水库对壤中流的日出流系数；K_G为栅格土壤自由水蓄水库对地下径流的日出流系数；S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，S为栅格当前时段土壤自由水蓄水库蓄水量，S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量；

步骤2-3，对于各未石漠化岩溶栅格，在土壤自由水蓄水库的下方构建表层岩溶带自由水蓄水库；上层土壤自由水蓄水库与下层表层岩溶带自由水蓄水库串联，定义表层岩溶带自由水蓄水库容量为S_MEK；

通过格林-安普特公式计算上层土壤自由水蓄水库的自由水向下层表层岩溶带自由水蓄水库渗透的过程：

式中：f(t)为栅格土壤自由水库下渗能力，K为栅格饱和水力传导度，Ψ为栅格湿润锋处土壤吸力，Δθ为栅格土壤饱和含水率与初始含水率之差，F(t)为栅格累计渗漏量，t为当前时段；

实际下渗量计算公式为：

f_a＝f(t)×Δt (9)

式中：f_a为实际下渗量，f(t)为栅格土壤自由水库下渗能力，Δt为时间步长；

基于栅格新安江模型的模拟方法，分别对于各未石漠化岩溶栅格土壤自由水蓄水库分水源调蓄的计算过程如下：

R_I＝K_I×S (12)

S_0t+1＝S-K_I×S-f_a (13)

式中：R为栅格总径流深，R_S为栅格地表径流深，R_I为栅格壤中流径流深，S_M为栅格土壤自由水蓄水库蓄水容量，S₀为栅格当前时段初始土壤自由水蓄水库蓄水量，S为栅格当前时段土壤自由水蓄水库蓄水量，f_a为栅格实际下渗量，S_0t+1为栅格下一时段初始土壤自由水蓄水库蓄水量；

下层表层岩溶带自由水蓄水库有两个出口，一个旁侧出口形成表层岩溶带侧向径流，一个向下出口形成岩溶地下径流：

S_EK＝S_EK0+f_a (14)

S_EK0t+1＝S_EK-K_EKI×S_EK-K_EKG×S_EK (16)

式中：S_EK0为栅格当前时段初始表层岩溶带自由水蓄水库蓄水量，S_EK为栅格当前时段表层岩溶带自由水蓄水库蓄水量，S_EK0t+1为栅格下一时段初始表层岩溶带自由水蓄水库蓄水量，R_EKI为栅格表层岩溶带侧向径流深，R_GK为岩溶栅格地下径流深，K_EKI为栅格表层岩溶带自由水蓄水库对表层岩溶带侧向径流的日出流系数；K_EKG为栅格表层岩溶带自由水蓄水库对岩溶地下径流的日出流系数；

步骤2-4，分别对于各石漠化岩溶栅格，该栅格的土壤自由水蓄水库只有一个向下的出口，构建的表层岩溶带自由水蓄水库与未石漠化岩溶栅格相同。

4.根据权利要求3所述的一种岩溶区分布式非线性水文模拟方法，其特征在于，所述步骤3包括：

步骤3-1，针对各岩溶栅格，建立各岩溶栅格表层岩溶带自由水蓄水库容量S_MEK与各岩溶栅格地形指数的线性关系，假设各岩溶栅格中地形指数最大的栅格对应的表层岩溶带自由水蓄水库容量最大，地形指数最小的栅格对应的表层岩溶带自由水蓄水库容量最小：

式中：TI_max为岩溶区域内各栅格地形指数中的最大值；TI_min为岩溶区域内各栅格地形指数中的最小值；S_MEKmax为表层岩溶带自由水蓄水库最大值，S_MEKmin为表层岩溶带自由水蓄水库最小值，ξ_a与ξ_b为线性关系的两个系数；

通过参数率定确定S_MEKmax和S_MEKmin，从而得到ξ_a与ξ_b两个系数的值；

进而基于步骤1得到的每个岩溶栅格的地形指数经下式，得到每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK：

ξ_a×TI+ξ_b＝S_MEK (18)

式中：TI为各栅格的地形指数；S_MEK为岩溶区域内各栅格表层岩溶带蓄水容量，ξ_a与ξ_b为线性关系的两个系数；

步骤3-2，假设K_EKI与K_EKG之和与表层岩溶带自由水蓄水库容量S_MEK呈线性关系，表层岩溶带蓄水容量最大的栅格对应的K_EKI与K_EKG之和最小，表层岩溶带蓄水容量最小的栅格对应的K_EKI与K_EKG之和最大；

通过参数率定确定K_EKI与K_EKG之和的最大值与最小值，从而得到a与b两个系数的值；

进而基于公式(18)得到的每个岩溶栅格的表层岩溶带自由水蓄水库容量S_MEK，经以下公式确定每个岩溶栅格K_EKI与K_EKG之和的值：

a×S_MEK+b＝(K_EKI+K_EKG) (20)

式中：S_MEK为岩溶区域内各栅格表层岩溶带自由水蓄水库容量，a与b为线性关系的两个系数；S_MEKmax、S_MEKmin为流域内各栅格表层岩溶带蓄水容量的最大值和最小值；

步骤3-3，假设K_EKI与K_EKG的比例为一个预设常数，通过参数率定确定该常数，即可确定各岩溶栅格的K_EKI和K_EKG值；

步骤3-4，基于各岩溶栅格对应的K_EKI、K_EKG值、以及表层岩溶带自由水蓄水库容量S_MEK，结合步骤2各岩溶栅格构建的表层岩溶带自由水蓄水库，得到各岩溶栅格的表层岩溶带侧向径流深R_EKI、岩溶地下径流深R_GK。

5.根据权利要求4所述的一种岩溶区分布式非线性水文模拟方法，其特征在于，所述步骤4包括：

步骤4-1，针对目标流域内各栅格采用线性水库法计算各径流坡面汇流量，

对于各非岩溶栅格：

式中：Q_S为栅格地表径流流量，Q_I为栅格壤中流流量，Q_G为非岩溶栅格地下径流流量，C_S为栅格地面径流消退系数；C_I为栅格壤中流消退系数；C_G为非岩溶栅格地下径流消退系数；t为当前时段，t-1为上一时段，U为单位换算系数，U＝栅格面积(km²)/3.6Δt，Δt为计算时间步长；

对于各岩溶栅格：

未石漠化岩溶栅格

石漠化岩溶栅格

式中：Q_EKI为表层岩溶带侧向出流流量，Q_GK为岩溶栅格地下径流流量，C_EKI为表层岩溶带侧向出流消退系数；C_GK为岩溶区地下径流消退系数；

步骤4-2，栅格间使用马斯京根法汇流将各栅格的流量向目标流域出口演算，第i个栅格流向第i+1个栅格的马斯京根法计算为：

Q(t)_i+1＝C₀Q(t)_i+C₁Q(t-1)_i+C₂Q(t-1)_i+1 (24)

式中：

X_e和K_e为马斯京根法的两个参数，Δt为计算时间步长，t为当前时段，t-1为上一时段；

目标流域内各栅格对应的Q_S、Q_I、Q_G、Q_EKI分别通过栅格间马斯京根法向目标流域出口演算，其中如果是岩溶栅格流向非岩溶栅格，则Q_EKI和Q_GK作为Q_G进行演算，最终得到目标流域出口断面的流量。

6.根据权利要求5所述的一种岩溶区分布式非线性水文模拟方法，其特征在于，所述步骤4-1中针对岩溶区地下径流消退系数C_GK，结合目标流域各栅格的岩溶发育情况，建立多个地下径流消退系数C_GK{C_GK1、C_GK2、…、C_GKn}，各地下径流消退系数分别与各岩溶发育程度一一对应，不同岩溶发育程度的栅格取不同的地下径流消退系数C_GKn。

Images