[发明专利]一种基于自然子流域的通用流域水循环模拟计算方法

权利要求书

1.一种基于自然子流域的通用流域水循环模拟计算方法，其特征在于，包括以下步骤：

步骤一：流域水文分析及子流域河道汇流计算顺序生成

基于DEM地形资料对研究流域进行水文分析，将研究流域划分为若干个子流域，明确各子流域相应下游子流域编号，采用拓扑关系模型分析出子流域河道汇流计算顺序；

步骤二：计算子流域面雨量过程、选择流域代表性蒸发过程

基于研究流域内及周边雨量站降雨资料，采用高容错面雨量计算方法计算得到各子流域面雨量过程；选取流域内代表站蒸发过程作为流域代表性蒸发过程：

对各自然子流域划分正方形网格，各网格以有效面积占自然子流域面积的比例为权重，自然子流域面雨量采用下式计算得到；

式中：n为子流域内网格个数；A_网格i为第i个网格有效面积；A_子流域为子流域面积；P_网格i为第i个网格雨量；

对于单一网格的雨量推求，通过排序算法搜索最接近网格中心点的N个雨量站点，根据站点雨量值判断、挑选出资料完备、准确的m个站点，m小于N，基于m个站点的雨量资料采用距离平方倒数法插值得，距离平方倒数法公式如下所示：

式中：为m个站点中第j个站点的雨量；d_ij为第j个站点距第i个网格中心点的距离；lat为纬度；lon为经度；

步骤三：根据气候及下垫面特性，选择各子流域产汇流计算方法

根据地形资料识别研究流域所属地形地貌；根据柯本气候分类识别研究流域内气候类别分布；绘制研究流域下垫面特性：土壤类别分布、土地覆被类别分布及坡度类别分布图；结合气候类别分布及下垫面特性分布，设定各子流域选取的产汇流计算方法；

步骤三：产汇流计算方法针对的子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；或2)气候湿润，地表易蓄满产流；或3)子流域缺乏代表性的蒸发资料输入；或4)子流域高寒高海拔，存在融雪径流；

针对子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；算法采用霍顿下渗曲线进行超渗产流计算；

霍顿下渗曲线的相关公式为：

f＝f_c+(f₀-f_c)e^-kt (4)

联立式(4)和式(5)，有：

式中：f为子流域时段内平均下渗率，mm/h；f₀、f_c分别为子流域平均最大、最小下渗能力，mm/h；t为历时，h；k为下渗能力衰减系数，h^-1；W为土壤含水量，mm；

需用迭代求出f～W的关系；迭代过程为：以T＝W/f₀作为t的第一次近似值，由式(5)计算得到W的第一个近似值ST，如果|ST-W|允许误差e，则由式(4)计算出f的第一个近似值U，然后T＝T+(W-ST)/U，迭代多次后，直到|ST-W|≤e，求得所需的f值；

针对子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；算法采用菲利普下渗曲线进行超渗产流计算；

菲利普下渗曲线的相关公式为：

式中：B和A为两个待定参数；其他参数含义如上；

菲利普下渗曲线，在给定一组系数A、B参数值，求得f～W的关系；

在得到f～W关系后，根据下式计算超渗产流量：

F＝f△t (11)

PE＝P-E (12)

f_mm＝f(1+BX) (13)

式中：R_IE为超渗产流量，mm；F为时段下渗量，mm；f_mm为流域平均下渗能力为f时流域最大的点下渗能力；BX为指数系数；

针对子流域不同水文气候条件为：2)气候湿润，地表易蓄满产流；算法采用基于地形指数的蓄满产流计算方法；

①蒸发量计算

式中：E_a,i为点i处实际蒸发量，m；E_P为蒸散发能力，m；S_rz,i为植被根系区缺水量，m；S_rmax,i为植被根系区最大蓄水容量，m；

②产流量计算

式中：a_i为点i处单宽集水面积，m²；tanβ_i为点i处的地表坡度；z_i为点i处地下水距地表深度，m；z为饱和地下水水面平均深度，m；S_zm为非饱和区最大蓄水深度，m；

若zi为负值，饱和地下水将漫出地面，形成地表径流；

点i处下渗率计算公式为：

式中：S_uz,i为点i处非饱和区土壤含水量，m；SD_i为非饱和区土壤蓄水能力，m；t_d为时间参数，h；

整个流域的下渗率为：

式中：A_i为地形指数数值相同的各处面积之和，m²；

Q_b＝AT₀exp(-λ)exp(-z/S_zm) (19)

式中：T₀为饱和导水率，m²/h；

饱和地下水水面平均深度z的计算公式为：

针对子流域不同水文气候条件为：2)气候湿润，地表易蓄满产流；算法采用蓄水容量～面积曲线的蓄满产流计算方法；

①蒸发量计算

当WU+P≥E_P时

E_U＝E_P E_L＝0 E_D＝0 (21)

当WU+PE_P，WL≥C·WLM时

E_U＝WU+P E_L＝(E_P-E_U)WL/WLM E_D＝0 (22)

当WU+PE_P，C(E_P-E_U)≤WLC·WLM时

E_U＝WU+P E_L＝C(E_P-E_U) E_D＝0 (23)

当WU+PE_P，WLC(E_P-E_U)时

E_U＝WU+P E_L＝WL E_D＝C(E_P-E_U)-E_L (24)

E＝E_U+E_L+E_D (25)

式中：E_P为蒸散发能力；P为降雨量；WL为下层土壤含水量；WU为上层土壤含水量；WLM为下层土壤含水容量；C为蒸发扩散系数；E_U为上层土壤蒸发量；E_L为下层土壤蒸发量；E_D为深层土壤蒸发量；E为总蒸发量；

②蓄满产流量计算

当a+PE≤WMM流域部分面积产流：

当a+PEWMM全流域产流：

R＝PE-(WM-W) (27)

其中：

式中：R为蓄满产流量，mm；a为与流域初始平均蓄水量W相应的流域最大蓄水量，mm；b为抛物线指数；WM为流域平均蓄水容量；WMM为流域最大蓄水容量；PE为扣除蒸发后的净雨；

③产流分配计算

当PE+AUS_mm：

当PE+AU≥S_mm：

其中：

FR＝R/PE (32)

时段自由水蓄水量为：

壤中产流量和地下产流量分别为：

RI＝KI·S·FR (34)

RG＝KG·S·FR (35)

式中：S_mm为流域最大自由水蓄水容量，mm；S_m为流域平均自由水蓄水容量，mm；EX为抛物线指数；S₁为时段初始流域平均自由水蓄水量，mm；AU为与S₁相应的流域最大自由水蓄水容量；FR₁和、FR分别为上一时段及本时段的产流面积比例；PE为扣除蒸发后的净雨；R为蓄满产流量；RS为地表产流量；RI为壤中产流量；RG为地下产流量；

针对子流域不同水文气候条件为：3)子流域缺乏代表性的蒸发资料输入；算法推荐采用基于增益因子的产流计算方法：

在基于增益因子的产流计算方法中，产流R表现为降雨P和增益G的乘积：

R(t)＝G(t)P(t) (36)

增益G与前期土壤含水量W有关，为：

经泰勒展开后：

G(t)＝g₁+g₂W(t) (38)

R(t)＝g₁P(t)+g₂W(t)P(t) (39)；

针对子流域不同水文气候条件为：4)子流域高寒高海拔，存在融雪径流；算法推荐采用基于增益因子的产流计算方法：

算法采用度日因子计算方法进行积融雪计算；

JR_snow＝D(T_t-T_c) (40)

式中：JR_snow为正时表示为融雪量，为负值时表示为积雪量，mm；D为度日因子，mm/d；T_t为日平均气温；T_c为临界气温，设为0℃；

若时段降雨为P，前期积雪深为S，则计算时段末积雪深为S-JR_snow限制不小于0，时段净雨为P+JR_snow(JR_snowS)或P+S(JR_snow≥S)；

此外，对于地表径流汇流计算，采用线性水库或地貌单位线方法；地下径流汇流计算，采用线性水库计算方法；

1)线性水库法

线性水库法公式为：

Q_t+1＝R_t+1(1-C)U+Q_t*C (41)

U＝AREA/(△t*3.6) (42)

式中：Q_t+1、Q_t分别为t+1、t时刻的流量，m³/s；R_t+1为t+1时刻产流量，mm；C为消退系数；U为单位转换系数；AREA为流域面积，km²；△t为时段长，h；

2)地貌单位线

地貌单位线形式如下：

式中：N为反映流域调蓄能力的参数，K为线性水库的蓄泄系数；Γ(N)为Γ函数，即

在参数N及参数K推求上，根据霍顿地貌几何率：面积比、河长比、分叉比推求：

式中：R_B,R_L,R_A—流域水系的分叉比、河长比和面积比，基于斯特拉勒级别通过DEM资料求得；

使用如下关系式：

τ＝1-(1-λ)(1-ρ) (45)

其中：

得到如下关系式：

τ＝λ^1-mλ (47)

由式(45)和式(47)推导出：

用霍顿河长定律，推导出：

以上式中：τ为净雨质点自河源至下游某一断面的平均汇流时间与河源至流域出口断面的平均汇流时间的比值；ρ为与河流长度及河底比降有关的参数；n为自河源至下游某断面的子河段数；N为自河源至流域出口断面的子河段数；△l_j为自河源开始划分的第j子河段长；p_j为第j子河段的平均坡度；m为反映河道纵剖面特性的综合参数；Ω为河系最高级别河流的级数；V_Ω为流域出口断面的流速，一般由出口断面洪水过程线的涨洪段的平均流速给定；α为流域形心至流域出口断面的距离与流域长度的比值；

m参数认为是反映河道纵剖面特性的综合参数；

根据霍顿河长定律及比降定律，针对ρ参数计算，构造河长比降比R_LS这一概念，其为各级河流lp^-0.5值的平均比值，则有：

在此基础上，联立式(45)及式(47)，通过迭代求解，计算参数m；

步骤四：耦合拓扑关系模型、水文模型及水力模型，优化产汇流参数

根据各子流域河道流经区域的土地利用情况，确定河道糙率值；耦合拓扑关系模型、水文产汇流计算模型及水力学模型，逐时段进行联合计算，以模拟得到的研究流域出口断面流量过程与实测流量过程的Nash效率系数作为参数优化对象，优化产汇流计算的参数，得到各子流域出流过程；

拓扑关系模型为自然子流域拓扑关系模型；算法为：步骤1)统计所有子流域的入度，入度即为汇流至某一子流域的相邻子流域个数；步骤2)分离出对于入度为0的子流域，将入度为0的子流域汇流至的相邻子流域的入度减1；步骤3)重复步骤2直至所有子流域被分离出来，完成子流域河道汇流计算顺序的排序计算；

水力学模型采用一维非恒定流模型，对于有n个子流域的研究流域，需建立(n-1)/2个河段的单河道一维非恒定流模型；模型建立包括连续方程和动量守恒方程；

1)连续方程

2)动量守恒方程

在有限差分中采用四点隐式差分格式，其中n和j分别表示时间和空间的离散，n时段的水位和流量已知，待求的是n+1时段的水位和流量；差分形式为：

按此形式写出各个项的差分并代入连续方程和动量守恒方程；对于方程中的非线性项进行适当线性化，有：

经过整理得：

A_1j△Q_j+B_1j△Z_j+C_1j△Q_j+1+D_1j△Z_j+1＝E_1j (57)

A_2j△Q_j+B_2j△Z_j+C_2j△Q_j+1+D_2j△Z_j+1＝E_2j (58)

式中A_1j、B_1j、C_1j、D_1j、E_1j、A_2j、B_2j、C_2j、D_2j、E_2j均为系数，如

上边界采用上游子流域入流，下边界采用相应单河道最下游断面的水位流量关系，由曼宁公式计算；

式中：Q为流量；n为糙率；A为断面过水面积；R为断面水力半径；S为断面位置水面比降，以底坡代替；

步骤五：若存在引水/调水过程，则分析对流域水循环的影响若特定子流域内存在引水/调水过程，通过概化为该子流域出流过程的增加/减少，根据步骤四优化得到的参数及其他已确定参数，计算各子流域出流过程，分析引水/调水过程对流域水循环的影响。