CN112651189B - 一种基于自然子流域的通用流域水循环模拟计算方法 - Google Patents

一种基于自然子流域的通用流域水循环模拟计算方法 Download PDF

Info

Publication number
CN112651189B
CN112651189B CN202011463721.8A CN202011463721A CN112651189B CN 112651189 B CN112651189 B CN 112651189B CN 202011463721 A CN202011463721 A CN 202011463721A CN 112651189 B CN112651189 B CN 112651189B
Authority
CN
China
Prior art keywords
basin
sub
river
formula
calculation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011463721.8A
Other languages
English (en)
Other versions
CN112651189A (zh
Inventor
杨开斌
夏建荣
韩兵
刘阳容
卢鹏
周鹏程
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
PowerChina Kunming Engineering Corp Ltd
Original Assignee
PowerChina Kunming Engineering Corp Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by PowerChina Kunming Engineering Corp Ltd filed Critical PowerChina Kunming Engineering Corp Ltd
Priority to CN202011463721.8A priority Critical patent/CN112651189B/zh
Publication of CN112651189A publication Critical patent/CN112651189A/zh
Application granted granted Critical
Publication of CN112651189B publication Critical patent/CN112651189B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Sewage (AREA)

Abstract

本发明涉及一种基于自然子流域的通用流域水循环模拟计算方法，本发明以自然子流域作为产流、坡面汇流的计算单元，各个子流域的面雨量输入采用本次发明的高容错面雨量计算方法计算得到。而各个子流域可采用不同的产汇流模型以充分反映气候、下垫面特性，子流域之间通过拓扑关系模型及水动力学模型实现河网汇流，最终计算得到各个子流域出口断面的流量过程。对于跨流域引调水工程影响，则概化为上游流量过程输入边界下，单一河道流量演进对于最邻近下游子流域出口出流的增减影响。本发明可为流域水资源分析计算、流域水资源优化配置提供有力的理论及技术支撑。

Description

一种基于自然子流域的通用流域水循环模拟计算方法
技术领域
本发明属于流域水循环模拟分析领域，更具体地，涉及一种模仿流域天然产汇流结构特性，可适应流域不同气候、下垫面分布条件，可定量分析跨流域引调水工程影响的流域水循环模拟计算方法。
背景技术
目前在流域水循环模拟分析领域，基于雨量站计算单元进行流域水循环模拟计算方法得到广泛应用。该方法根据流域地形、地貌条件及布设的雨量站网，用泰森多边形法将研究流域划分为若干个雨量站计算单元，每个计算单元对应于相应雨量站的泰森多边形，计算面积即为泰森多边形所占面积，并认为单元内降雨输入可由相应雨量站雨量过程代表。对每个计算单元采用水文模型分别进行蒸散发计算、产流计算及坡面汇流计算，得到计算单元的出流过程；将单元面积的出流过程用马斯京根分段连续演算法进行出口以下的河道汇流演算，求得单元面积在流域出口的流量过程线；将每个单元面积在流域出口的流量过程线性叠加，即得到研究流域出口断面的流量过程。
该方法计算逻辑简单，计算量较小，成果也能有一定精度，因此在水文预报、水循环模拟计算等领域得到广泛应用。
但该方法存在以下无法忽视的问题：
(1)所有雨量站计算单元均采用同一套产流、坡面汇流计算方法计算单元出口出流过程，但实际上单元出口位置无法得知，无法得到流域内有意义的分区流量成果，更重要的是每个雨量站计算单元仅具有空间统计意义，而缺乏集水单元的物理意义。对于研究流域内气候、下垫面条件存在明显差异的情况，没有理论基础也无法针对特定雨量站计算单元选择不同的产流、坡面汇流计算方法来适应流域这种空间差异性，影响了水循环模拟计算的精度，也影响了方法的通用性、实用性。
(2)此外，对于各个雨量站计算单元的出流过程，采用马斯京根法进行河道汇流，叠加得到研究流域出口出流的方法，存在着因为使用的是虚拟河道而不是实际河道从而导致的河网汇流物理意义不足的问题，也进一步影响了水循环模拟计算的精度。
(3)再有，该方法由于无法进行汇水节点及实际河道的有效表达，无法处理当流域内存在跨流域引调水工程时引调水影响的定量分析计算问题，也进一步影响了方法的通用性、实用性。
发明内容
针对现有基于雨量站计算单元进行流域水循环模拟计算方法在实际应用中的实用性、通用性不高，计算精度有限的问题，提出一种基于自然子流域的通用流域水循环模拟计算方法，该方法以自然子流域作为产流、坡面汇流的计算单元，各个子流域的面雨量输入采用本次发明的高容错面雨量计算方法计算得到。而各个子流域可采用不同的产汇流模型以充分反映气候、下垫面特性，子流域之间通过拓扑关系模型及水力学模型实现河网汇流，最终计算得到各个子流域出口断面的流量过程。对于跨流域引调水工程影响，则概化为上游流量过程输入边界下，单一河道流量演进对于最邻近下游子流域出口出流的增减影响。本发明解决了原有方法无法针对性考虑流域内不同水文气候条件的子单元产汇流计算、河网汇流缺乏物理意义导致计算精度不高、无法定量化分析跨流域来水/出流对流域水循环的影响的三大问题，可为流域水资源分析计算、流域水资源优化配置提供有力的理论及技术支撑。
具体技术解决方案如下：
(1)高容错面雨量计算方法解决自然子流域面雨量计算问题
对各自然子流域划分正方形网格，各网格以有效面积占自然子流域面积的比例为权重，自然子流域面雨量采用下式计算得到。
Figure GDA0003798180280000021
式中：n为子流域内网格个数；A网格i为第i个网格有效面积；A子流域为子流域面积；P网格i为第i个网格雨量。
对于单一网格的雨量推求，通过排序算法搜索最接近网格中心点的N个雨量站点，根据站点雨量值判断(例如***设定-999为缺测或传输不及时，那检测到-999即判断该站点该时段的雨量缺失，此外还可加入其它雨量值合理性判断方法)、挑选出资料完备、准确的m(m小于N)个站点，基于m个站点的雨量资料采用距离平方倒数法插值算得，距离平方倒数法公式如下所示：
Figure GDA0003798180280000022
Figure GDA0003798180280000023
式中：P站点j为m个站点中第j个站点的雨量；dij为第j个站点距第i个网格中心点的距离；lat为纬度；lon为经度。
通过以上设定，方法具备以下优点：
①不以存在“变数”的雨量站为加权基本，而是选择以网格雨量按面积加权求取子流域面雨量，计算稳定，成果精度高，网格生成及面积权重计算均为一次性工作，计算复杂度不高。
②筛选、挑选站点进行距离平方倒数插值网格雨量的过程，对资料质量要求不高，容错能力强，计算精度高且速度快，雨量站变动后只需更新站网位置信息，维护简单。
本次发明的高容错面雨量计算方法可稳定、高效地解决雨量站资料偶有缺失、传输不及时情况下的各自然子流域的面雨量计算问题。
(2)构建产汇流计算方法库解决不同水文气候条件的子流域产汇流计算问题
经过归纳，子流域不同水文气候条件可基本分为以下情况：
1)气候干旱，土层厚度大，超渗产流显著；
2)气候湿润，地表易蓄满产流；
3)子流域缺乏代表性的蒸发资料输入；
4)子流域高寒高海拔，存在融雪径流；
在建立产汇流方法库的前提下，算法可通过方法编号来根据需求配置各子流域适合的产汇流方法。
对于第一种情况，算法可采用霍顿下渗曲线或菲利普下渗曲线进行超渗产流计算；
1、霍顿下渗曲线
霍顿下渗曲线的相关公式为：
f＝fc+(f0-fc)e-kt (4)
Figure GDA0003798180280000031
联立式(4)和式(5)，有：
Figure GDA0003798180280000032
式中：f为子流域时段内平均下渗率，mm/h；f0、fc分别为子流域平均最大、最小下渗能力，mm/h；t为历时，h；k为下渗能力衰减系数，h-1；W为土壤含水量，mm。
需用迭代求出f～W的关系。迭代过程为：以T＝W/f0作为t的第一次近似值，即可由式(5)计算得到W的第一个近似值ST，如果|ST-W|>允许误差e，则由式(4)计算出f的第一个近似值U，然后T＝T+(W-ST)/U，迭代多次后，直到|ST-W|≤e，即可求得所需的f值。
2、菲利普下渗曲线
菲利普下渗曲线的相关公式为：
Figure GDA0003798180280000041
Figure GDA0003798180280000042
Figure GDA0003798180280000043
式中：B和A为两个待定参数；其他参数含义如上。
菲利普下渗曲线，在给定一组系数A、B参数值，便可直接求得f～W的关系。
在得到f～W关系后，为充分考虑下垫面土壤下渗的变异性，假定流域内各点的下渗能力在流域上的分布为抛物线，根据下式计算超渗产流量：
Figure GDA0003798180280000044
F＝f△t (11)
PE＝P-E (12)
fmm＝f(1+BX) (13)
式中：RIE为超渗产流量，mm；F为时段下渗量，mm；fmm为流域平均下渗能力为f时流域最大的点下渗能力；BX为指数系数。
对于第二种情况，算法可采用基于地形指数或蓄水容量～面积曲线的蓄满产流计算方法。
1、基于地形指数蓄满产流计算方法
①蒸发量计算
Figure GDA0003798180280000045
式中：Ea,i为点i处实际蒸发量，m；EP为蒸散发能力，m；Srz,i为植被根系区缺水量，m；Srmax,i为植被根系区最大蓄水容量，m。
②产流量计算
Figure GDA0003798180280000046
Figure GDA0003798180280000047
式中：ai为点i处单宽集水面积，m2；tanβi为点i处的地表坡度；zi为点i处地下水距地表深度，m；
Figure GDA0003798180280000051
为饱和地下水水面平均深度，m；Szm为非饱和区最大蓄水深度，m。
若zi为负值，饱和地下水将漫出地面，形成地表径流。
点i处下渗率计算公式为：
Figure GDA0003798180280000052
式中：Suz,i为点i处非饱和区土壤含水量，m；SDi为非饱和区土壤蓄水能力，m；td为时间参数，h。
整个流域的下渗率为：
Figure GDA0003798180280000053
式中：Ai为地形指数数值相同的各处面积之和，m2。
Figure GDA0003798180280000054
式中：T0为饱和导水率，m2/h。
饱和地下水水面平均深度
Figure GDA0003798180280000055
的计算公式为：
Figure GDA0003798180280000056
2、基于蓄水容量～面积曲线的蓄满产流计算方法
①蒸发量计算
当WU+P≥EP时
EU＝EP EL＝0 ED＝0 (21)
当WU+P<EP，WL≥C·WLM时
EU＝WU+P EL＝(EP-EU)WL/WLM ED＝0 (22)
当WU+P<EP，C(EP-EU)≤WL<C·WLM时
EU＝WU+P EL＝C(EP-EU) ED＝0 (23)
当WU+P<EP，WL<C(EP-EU)时
EU＝WU+P EL＝WL ED＝C(EP-EU)-EL (24)
E＝EU+EL+ED (25)
式中：EP为蒸散发能力；P为降雨量；WL为下层土壤含水量；WU为上层土壤含水量；WLM为下层土壤含水容量；C为蒸发扩散系数；EU为上层土壤蒸发量；EL为下层土壤蒸发量；ED为深层土壤蒸发量；E为总蒸发量。
②蓄满产流量计算
当a+PE≤WMM(流域部分面积产流)：
Figure GDA0003798180280000061
当a+PE>WMM(全流域产流)：
R＝PE-(WM-W) (27)
其中：
Figure GDA0003798180280000062
式中：R为蓄满产流量，mm；a为与流域初始平均蓄水量W相应的流域最大蓄水量，mm；b为抛物线指数；WM为流域平均蓄水容量；WMM为流域最大蓄水容量；PE为扣除蒸发后的净雨。
③产流分配计算
当PE+AU<Smm：
Figure GDA0003798180280000063
当PE+AU≥Smm：
Figure GDA0003798180280000064
其中：
Figure GDA0003798180280000065
FR＝R/PE (32)
时段自由水蓄水量为：
Figure GDA0003798180280000066
壤中产流量和地下产流量分别为：
RI＝KI·S·FR (34)
RG＝KG·S·FR (35)
式中：Smm为流域最大自由水蓄水容量，mm；Sm为流域平均自由水蓄水容量，mm；EX为抛物线指数；S1为时段初始流域平均自由水蓄水量，mm；AU为与S1相应的流域最大自由水蓄水容量；FR1和、FR分别为上一时段及本时段的产流面积比例；PE为扣除蒸发后的净雨；R为蓄满产流量；RS为地表产流量；RI为壤中产流量；RG为地下产流量。
对于第三种情况，算法推荐采用基于增益因子的产流计算方法。
在基于增益因子的产流计算方法中，产流R表现为降雨P和增益G的乘积：
R(t)＝G(t)P(t) (36)
增益G与前期土壤含水量W有关，可表示为：
Figure GDA0003798180280000071
经泰勒展开后，有：
G(t)＝g1+g2W(t) (38)
R(t)＝g1P(t)+g2W(t)P(t) (39)
对于第四种情况，算法可采用度日因子计算方法进行积融雪计算。
JRsnow＝D(Tt-Tc) (40)
式中：JRsnow为正时表示为融雪量，为负值时表示为积雪量，mm；D为度日因子，mm/d；Tt为日平均气温；Tc为临界气温，一般可设为0℃。
若时段降雨为P，前期积雪深为S，则计算时段末积雪深为S-JRsnow(限制不小于0)，时段净雨为P+JRsnow(JRsnow<S)或P+S(JRsnow≥S)。
此外，对于地表径流汇流计算，可采用线性水库或地貌单位线方法。地下径流汇流计算，采用线性水库计算方法。
1、线性水库法
线性水库法公式为：
Qt+1＝Rt+1(1-C)U+Qt*C (41)
U＝AREA/(△t*3.6) (42)
式中：Qt+1、Qt分别为t+1、t时刻的流量，m3/s；Rt+1为t+1时刻产流量，mm；C为消退系数；U为单位转换系数；AREA为流域面积，km2；△t为时段长，h；
2、地貌单位线
地貌单位线基本形式如下：
Figure GDA0003798180280000081
式中：N为反映流域调蓄能力的参数，K为线性水库的蓄泄系数；Γ(N)为Γ函数，即
Figure GDA0003798180280000082
在参数N及参数K推求上，可根据霍顿地貌几何率(面积比、河长比、分叉比)推求：
Figure GDA0003798180280000083
式中：RB,RL,RA—流域水系的分叉比、河长比和面积比，可基于斯特拉勒级别通过DEM资料求得。
推求K参数的问题实质是如何根据地形资料确定流域平均汇流时间。根据不同级别河流的流速主要依赖于其地形坡度的事实，有如下关系式：
τ＝1-(1-λ)(1-ρ) (45)
其中：
Figure GDA0003798180280000084
进一步分析还可得到如下关系式：
τ＝λ1-mλ (47)
由式(45)和式(47)还可推导出：
Figure GDA0003798180280000085
用霍顿河长定律，可推导出：
Figure GDA0003798180280000086
以上式中：τ为净雨质点自河源至下游某一断面的平均汇流时间与河源至流域出口断面的平均汇流时间的比值；ρ为与河流长度及河底比降有关的参数；n为自河源至下游某断面的子河段数；N为自河源至流域出口断面的子河段数；△lj为自河源开始划分的第j子河段长；pj为第j子河段的平均坡度；m为反映河道纵剖面特性的综合参数；Ω为河系最高级别河流的级数；VΩ为流域出口断面的流速，一般由出口断面洪水过程线的涨洪段的平均流速给定；α为流域形心至流域出口断面的距离与流域长度的比值。
m参数可认为是反映河道纵剖面特性的综合参数，目前需结合干支流的实际资料，先计算τ和λ，并点绘τ～λ图进行分析得到，较为繁琐。
根据霍顿河长定律及比降定律，针对ρ参数计算，可构造河长比降比RLS这一概念，其为各级河流lp-0.5值的平均比值，则有：
Figure GDA0003798180280000091
在此基础上，可联立式(45)及式(47)，通过迭代求解，较方便的计算参数m。
相比现有基于雨量站计算单元进行流域水循环模拟计算方法，通过构建针对不同水文气候条件的产汇流计算方法库，可灵活配置各自然子流域的方法编号，可在充分适应各自然子流域的水文气候特性的基础上，提高计算子流域自身产汇流反映在子流域出口的流量过程的精度。
(3)拓扑关系模型解决子流域河道汇流计算顺序问题
自然子流域拓扑关系为二叉树结构。
经分析可知，子流域拓扑关系图，也是有向无环图，则可构建自然子流域拓扑关系模型的算法为：
1、统计所有子流域的入度(入度即为汇流至某一子流域的相邻子流域个数)；
2、分离出对于入度为0的子流域，将入度为0的子流域汇流至的相邻子流域的入度减1；
3、重复步骤2直至所有子流域被分离出来，完成子流域河道汇流计算顺序的排序计算。
(4)水力学模型解决子流域河道汇流计算、跨流域引调水工程影响定量分析问题
水力学模型采用一维非恒定流模型，根据二叉树结构的拓扑关系分析可知，对于有n个子流域的研究流域，需建立(n-1)/2个河段的单河道一维非恒定流模型。模型建立的主要依据为一维圣维南方程组，包括连续方程和动量守恒方程。
1)连续方程
Figure GDA0003798180280000092
2)动量守恒方程
Figure GDA0003798180280000093
在有限差分中采用四点隐式差分格式，其中n和j分别表示时间和空间(沿河道)的离散，n时段的水位和流量已知，待求的是n+1时段的水位和流量。差分形式为：
Figure GDA0003798180280000101
Figure GDA0003798180280000102
Figure GDA0003798180280000103
按此形式写出各个项的差分并代入连续方程和动量守恒方程。对于方程中的非线性项进行适当线性化，例如有：
Figure GDA0003798180280000104
经过整理可得：
A1j△Qj+B1j△Zj+C1j△Qj+1+D1j△Zj+1＝E1j (57)
A2j△Qj+B2j△Zj+C2j△Qj+1+D2j△Zj+1＝E2j (58)
式中A1j、B1j、C1j、D1j、E1j、A2j、B2j、C2j、D2j、E2j均为系数，如
Figure GDA0003798180280000105
Figure GDA0003798180280000106
上边界采用上游子流域入流，下边界采用相应单河道最下游断面的水位流量关系，可由曼宁公式计算。
Figure GDA0003798180280000107
式中：Q为流量；n为糙率；A为断面过水面积；R为断面水力半径；S为断面位置水面比降，以底坡代替。
发明的效果
与现有技术相比，本发明具有以下优点和有益效果：
1、本发明通过耦合拓扑关系模型、水文模型、水力学模型，对原有基于雨量站单元“先演后合”的水循环计算方法进行了改进，自然子流域分别进行产汇流计算后，再根据实际的河网结构进行有序的河网汇流，改进后模型结构更符合流域天然产汇流特性，模型参数也更具物理意义，提高了模拟精度；
2、相比原方法，本发明方法可基于拓扑关系模型及单河道一维非恒定流模型，灵活分析跨流域引调水工程对流域水循环的影响。
3、本发明通过构建产汇流计算方法库，可适应流域不同气候、下垫面分布条件，针对性的进行流域内不同水文气候条件子单元的产汇流计算。基于汇水节点及实际河道的有效表达，除流域出口断面流量过程外，方法还可获得各个子流域的流量过程，在实际应用中具有较高价值。
总体来说，本发明方法相比原方法更为符合流域天然产汇流特性，提高了水循环模拟精度，方法通用性更强，应用面更广，可为流域水循环模拟、水文预报等研究工作提供有力的理论与技术支撑。
附图说明
图1为本发明基于自然子流域的通用流域水循环模拟计算方法实施流程图；
图2为本发明研究流域-虎山水文站控制流域子流域分布图；
图3为本发明研究流域-虎山水文站控制流域土壤类别分布图；
图4为本发明研究流域-虎山水文站控制流域土地覆被类别分布图；
图5为本发明研究流域-虎山水文站控制流域坡度类别分布图；
图6为本发明方法计算得到的研究流域-虎山水文站控制流域各子流域出口流量过程图；
图7为本发明方法与原有方法在研究流域-虎山水文站控制流域出口流量过程的对比图。
图3中：(图例中的3位字母字符是在FAO90标准下的土壤索引代码。其中，LVh代指弱发育淋溶土，CMd代指不饱和始成土，FLe代指饱和冲积土，ATc代指土垫旱耕人为土，ACh代指典型低活性强淋溶土，ACn代指典型腐殖质强淋溶土，ALn代指典型高活性强淋溶土，WR代指水体。数字代表流域内存在土种相同但质地不同的亚类)。
图4中：(RICE代指稻田；AGRL代指耕地；FRST代指林地；PAST代指草地；URHD代指高密度居民区；BARR代指裸土；WATR代指水体)
具体实施方式
以江西饶河虎山水文站控制流域为研究对象进行基于自然子流域的流域水循环模拟分析。
本发明主要实施步骤如下(主要流程见图1)：
步骤一：基于DEM地形资料采用ArcGIS软件对研究流域-虎山水文站控制流域进行水文分析，将研究流域划分为13个子流域，如图2所示，提取各子流域内河道。明确各子流域相应下游子流域编号，采用拓扑关系模型分析出子流域河道汇流计算顺序。
表1基于拓扑关系模型分析出的虎山水文站控制流域子流域河道汇流计算顺序
序号 子流域编号 下游子流域编号 计算顺序号
1 1 3 2
2 2 3 1
3 3 5 4
4 4 5 3
5 5 7 6
6 6 7 5
7 7 10 8
8 8 8 13
9 9 10 7
10 10 11 10
11 11 8 12
12 12 8 11
13 13 11 9
步骤二：基于虎山水文站控制流域内及周边44个雨量站2009年～2013年逐日降雨资料，采用高容错面雨量计算方法计算得到各子流域2009年～2013年面雨量过程。选取三都站2009年～2013年逐日蒸发过程作为流域代表性蒸发过程。
步骤三：根据地形资料识别研究流域属于山区、丘陵地貌。根据柯本气候分类识别研究流域内均属于CFA(温带-无干季-夏季炎热)气候。基于土壤、土地覆被及地形资料，绘制研究流域土壤类别分布、土地覆被类别分布及坡度类别分布图，见图3～图5。在此基础上，分析子流域5、7、8、10、11、12存在较多人为活动，产汇流计算方法选择适用于山区地形可反映一定人为活动的新安江模型方法，而其他子流域考虑采用结构简单、计算效率较高的水箱模型方法。
步骤四：根据各子流域河道流经区域的土地利用情况，经验确定河道糙率值。耦合拓扑关系模型、水文(产汇流计算)模型及水力学模型，逐时段进行联合计算，以模拟得到的研究流域出口断面流量过程与实测流量过程的Nash效率系数作为参数优化对象，优化产汇流计算的参数，得到各子流域出流过程，见图6。本发明方法与原有基于雨量站计算单元进行流域水循环模拟计算方法在研究流域出口断面模拟流量过程与实测的对比情况见图7。
表2本发明方法与原有方法在虎山水文站控制流域出口流量过程的指标对比
指标 本发明方法 原有方法 是否本发明方法指标表现更优
Nash效率系数 0.968 0.954 是
水量平衡系数 0.998 1.038 是
步骤五：在子流域13增加调水过程，通过概化为子流域13出流过程的减少，根据步骤四优化得到的参数及其他已确定参数，计算各子流域出流过程，分析调水过程对流域水循环的影响。
以上所述的仅是本发明的部分具体实施例(本发明的技术方案的实施例不能穷举，本发明所记载的保护范围以本发明的数值范围和其他技术要点范围为准)，方案中公知的具体内容或常识在此未作过多描述。应当指出，上述实施例不以任何方式限制本发明，对于本领域的技术人员来说，凡是采用等同替换或等效变换的方式获得的技术方案均落在本发明的保护范围内。本申请要求的保护范围应当以其权利要求的内容为准，说明书中的具体实施方式等记载可以用于解释权利要求的内容。

Claims (1)

1.一种基于自然子流域的通用流域水循环模拟计算方法，其特征在于，包括以下步骤：
步骤一：流域水文分析及子流域河道汇流计算顺序生成
基于DEM地形资料对研究流域进行水文分析，将研究流域划分为若干个子流域，明确各子流域相应下游子流域编号，采用拓扑关系模型分析出子流域河道汇流计算顺序；
步骤二：计算子流域面雨量过程、选择流域代表性蒸发过程
基于研究流域内及周边雨量站降雨资料，采用高容错面雨量计算方法计算得到各子流域面雨量过程；选取流域内代表站蒸发过程作为流域代表性蒸发过程：
对各自然子流域划分正方形网格，各网格以有效面积占自然子流域面积的比例为权重，自然子流域面雨量采用下式计算得到；
Figure FDA0003798180270000011
式中：n为子流域内网格个数；A网格i为第i个网格有效面积；A子流域为子流域面积；P网格i为第i个网格雨量；
对于单一网格的雨量推求，通过排序算法搜索最接近网格中心点的N个雨量站点，根据站点雨量值判断、挑选出资料完备、准确的m个站点，m小于N，基于m个站点的雨量资料采用距离平方倒数法插值得，距离平方倒数法公式如下所示：
Figure FDA0003798180270000012
Figure FDA0003798180270000013
式中：
Figure FDA0003798180270000021
为m个站点中第j个站点的雨量；dij为第j个站点距第i个网格中心点的距离；lat为纬度；lon为经度；
步骤三：根据气候及下垫面特性，选择各子流域产汇流计算方法
根据地形资料识别研究流域所属地形地貌；根据柯本气候分类识别研究流域内气候类别分布；绘制研究流域下垫面特性：土壤类别分布、土地覆被类别分布及坡度类别分布图；结合气候类别分布及下垫面特性分布，设定各子流域选取的产汇流计算方法；
步骤三：产汇流计算方法针对的子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；或2)气候湿润，地表易蓄满产流；或3)子流域缺乏代表性的蒸发资料输入；或4)子流域高寒高海拔，存在融雪径流；
针对子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；算法采用霍顿下渗曲线进行超渗产流计算；
霍顿下渗曲线的相关公式为：
f＝fc+(f0-fc)e-kt (4)
Figure FDA0003798180270000022
联立式(4)和式(5)，有：
Figure FDA0003798180270000023
式中：f为子流域时段内平均下渗率，mm/h；f0、fc分别为子流域平均最大、最小下渗能力，mm/h；t为历时，h；k为下渗能力衰减系数，h-1；W为土壤含水量，mm；
需用迭代求出f～W的关系；迭代过程为：以T＝W/f0作为t的第一次近似值，由式(5)计算得到W的第一个近似值ST，如果|ST-W|>允许误差e，则由式(4)计算出f的第一个近似值U，然后T＝T+(W-ST)/U，迭代多次后，直到|ST-W|≤e，求得所需的f值；
针对子流域不同水文气候条件为：1)气候干旱，土层厚度大，超渗产流显著；算法采用菲利普下渗曲线进行超渗产流计算；
菲利普下渗曲线的相关公式为：
Figure FDA0003798180270000031
Figure FDA0003798180270000032
Figure FDA0003798180270000033
式中：B和A为两个待定参数；其他参数含义如上；
菲利普下渗曲线，在给定一组系数A、B参数值，求得f～W的关系；
在得到f～W关系后，根据下式计算超渗产流量：
Figure FDA0003798180270000034
F＝f△t (11)
PE＝P-E (12)
fmm＝f(1+BX) (13)
式中：RIE为超渗产流量，mm；F为时段下渗量，mm；fmm为流域平均下渗能力为f时流域最大的点下渗能力；BX为指数系数；
针对子流域不同水文气候条件为：2)气候湿润，地表易蓄满产流；算法采用基于地形指数的蓄满产流计算方法；
①蒸发量计算
Figure FDA0003798180270000035
式中：Ea,i为点i处实际蒸发量，m；EP为蒸散发能力，m；Srz,i为植被根系区缺水量，m；Srmax,i为植被根系区最大蓄水容量，m；
②产流量计算
Figure FDA0003798180270000041
Figure FDA0003798180270000042
式中：ai为点i处单宽集水面积，m2；tanβi为点i处的地表坡度；zi为点i处地下水距地表深度，m；z为饱和地下水水面平均深度，m；Szm为非饱和区最大蓄水深度，m；
若zi为负值，饱和地下水将漫出地面，形成地表径流；
点i处下渗率计算公式为：
Figure FDA0003798180270000043
式中：Suz,i为点i处非饱和区土壤含水量，m；SDi为非饱和区土壤蓄水能力，m；td为时间参数，h；
整个流域的下渗率为：
Figure FDA0003798180270000044
式中：Ai为地形指数数值相同的各处面积之和，m2；
Qb＝AT0exp(-λ)exp(-z/Szm) (19)
式中：T0为饱和导水率，m2/h；
饱和地下水水面平均深度z的计算公式为：
Figure FDA0003798180270000045
针对子流域不同水文气候条件为：2)气候湿润，地表易蓄满产流；算法采用蓄水容量～面积曲线的蓄满产流计算方法；
①蒸发量计算
当WU+P≥EP时
EU＝EP EL＝0 ED＝0 (21)
当WU+P<EP，WL≥C·WLM时
EU＝WU+P EL＝(EP-EU)WL/WLM ED＝0 (22)
当WU+P<EP，C(EP-EU)≤WL<C·WLM时
EU＝WU+P EL＝C(EP-EU) ED＝0 (23)
当WU+P<EP，WL<C(EP-EU)时
EU＝WU+P EL＝WL ED＝C(EP-EU)-EL (24)
E＝EU+EL+ED (25)
式中：EP为蒸散发能力；P为降雨量；WL为下层土壤含水量；WU为上层土壤含水量；WLM为下层土壤含水容量；C为蒸发扩散系数；EU为上层土壤蒸发量；EL为下层土壤蒸发量；ED为深层土壤蒸发量；E为总蒸发量；
②蓄满产流量计算
当a+PE≤WMM流域部分面积产流：
Figure FDA0003798180270000051
当a+PE>WMM全流域产流：
R＝PE-(WM-W) (27)
其中：
Figure FDA0003798180270000052
式中：R为蓄满产流量，mm；a为与流域初始平均蓄水量W相应的流域最大蓄水量，mm；b为抛物线指数；WM为流域平均蓄水容量；WMM为流域最大蓄水容量；PE为扣除蒸发后的净雨；
③产流分配计算
当PE+AU<Smm：
Figure FDA0003798180270000061
当PE+AU≥Smm：
Figure FDA0003798180270000062
其中：
Figure FDA0003798180270000063
FR＝R/PE (32)
时段自由水蓄水量为：
Figure FDA0003798180270000064
壤中产流量和地下产流量分别为：
RI＝KI·S·FR (34)
RG＝KG·S·FR (35)
式中：Smm为流域最大自由水蓄水容量，mm；Sm为流域平均自由水蓄水容量，mm；EX为抛物线指数；S1为时段初始流域平均自由水蓄水量，mm；AU为与S1相应的流域最大自由水蓄水容量；FR1和、FR分别为上一时段及本时段的产流面积比例；PE为扣除蒸发后的净雨；R为蓄满产流量；RS为地表产流量；RI为壤中产流量；RG为地下产流量；
针对子流域不同水文气候条件为：3)子流域缺乏代表性的蒸发资料输入；算法推荐采用基于增益因子的产流计算方法：
在基于增益因子的产流计算方法中，产流R表现为降雨P和增益G的乘积：
R(t)＝G(t)P(t) (36)
增益G与前期土壤含水量W有关，为：
Figure FDA0003798180270000065
经泰勒展开后：
G(t)＝g1+g2W(t) (38)
R(t)＝g1P(t)+g2W(t)P(t) (39)；
针对子流域不同水文气候条件为：4)子流域高寒高海拔，存在融雪径流；算法推荐采用基于增益因子的产流计算方法：
算法采用度日因子计算方法进行积融雪计算；
JRsnow＝D(Tt-Tc) (40)
式中：JRsnow为正时表示为融雪量，为负值时表示为积雪量，mm；D为度日因子，mm/d；Tt为日平均气温；Tc为临界气温，设为0℃；
若时段降雨为P，前期积雪深为S，则计算时段末积雪深为S-JRsnow限制不小于0，时段净雨为P+JRsnow(JRsnow<S)或P+S(JRsnow≥S)；
此外，对于地表径流汇流计算，采用线性水库或地貌单位线方法；地下径流汇流计算，采用线性水库计算方法；
1)线性水库法
线性水库法公式为：
Qt+1＝Rt+1(1-C)U+Qt*C (41)
U＝AREA/(△t*3.6) (42)
式中：Qt+1、Qt分别为t+1、t时刻的流量，m3/s；Rt+1为t+1时刻产流量，mm；C为消退系数；U为单位转换系数；AREA为流域面积，km2；△t为时段长，h；
2)地貌单位线
地貌单位线形式如下：
Figure FDA0003798180270000071
式中：N为反映流域调蓄能力的参数，K为线性水库的蓄泄系数；Γ(N)为Γ函数，即
Figure FDA0003798180270000081
在参数N及参数K推求上，根据霍顿地貌几何率：面积比、河长比、分叉比推求：
Figure FDA0003798180270000082
式中：RB,RL,RA—流域水系的分叉比、河长比和面积比，基于斯特拉勒级别通过DEM资料求得；
使用如下关系式：
τ＝1-(1-λ)(1-ρ) (45)
其中：
Figure FDA0003798180270000083
得到如下关系式：
τ＝λ1-mλ (47)
由式(45)和式(47)推导出：
Figure FDA0003798180270000084
用霍顿河长定律，推导出：
Figure FDA0003798180270000085
以上式中：τ为净雨质点自河源至下游某一断面的平均汇流时间与河源至流域出口断面的平均汇流时间的比值；ρ为与河流长度及河底比降有关的参数；n为自河源至下游某断面的子河段数；N为自河源至流域出口断面的子河段数；△lj为自河源开始划分的第j子河段长；pj为第j子河段的平均坡度；m为反映河道纵剖面特性的综合参数；Ω为河系最高级别河流的级数；VΩ为流域出口断面的流速，一般由出口断面洪水过程线的涨洪段的平均流速给定；α为流域形心至流域出口断面的距离与流域长度的比值；
m参数认为是反映河道纵剖面特性的综合参数；
根据霍顿河长定律及比降定律，针对ρ参数计算，构造河长比降比RLS这一概念，其为各级河流lp-0.5值的平均比值，则有：
Figure FDA0003798180270000091
在此基础上，联立式(45)及式(47)，通过迭代求解，计算参数m；
步骤四：耦合拓扑关系模型、水文模型及水力模型，优化产汇流参数
根据各子流域河道流经区域的土地利用情况，确定河道糙率值；耦合拓扑关系模型、水文产汇流计算模型及水力学模型，逐时段进行联合计算，以模拟得到的研究流域出口断面流量过程与实测流量过程的Nash效率系数作为参数优化对象，优化产汇流计算的参数，得到各子流域出流过程；
拓扑关系模型为自然子流域拓扑关系模型；算法为：步骤1)统计所有子流域的入度，入度即为汇流至某一子流域的相邻子流域个数；步骤2)分离出对于入度为0的子流域，将入度为0的子流域汇流至的相邻子流域的入度减1；步骤3)重复步骤2直至所有子流域被分离出来，完成子流域河道汇流计算顺序的排序计算；
水力学模型采用一维非恒定流模型，对于有n个子流域的研究流域，需建立(n-1)/2个河段的单河道一维非恒定流模型；模型建立包括连续方程和动量守恒方程；
1)连续方程
Figure FDA0003798180270000101
2)动量守恒方程
Figure FDA0003798180270000102
在有限差分中采用四点隐式差分格式，其中n和j分别表示时间和空间的离散，n时段的水位和流量已知，待求的是n+1时段的水位和流量；差分形式为：
Figure FDA0003798180270000103
Figure FDA0003798180270000104
Figure FDA0003798180270000105
按此形式写出各个项的差分并代入连续方程和动量守恒方程；对于方程中的非线性项进行适当线性化，有：
Figure FDA0003798180270000106
经过整理得：
A1j△Qj+B1j△Zj+C1j△Qj+1+D1j△Zj+1＝E1j (57)
A2j△Qj+B2j△Zj+C2j△Qj+1+D2j△Zj+1＝E2j (58)
式中A1j、B1j、C1j、D1j、E1j、A2j、B2j、C2j、D2j、E2j均为系数，如
Figure FDA0003798180270000107
Figure FDA0003798180270000108
上边界采用上游子流域入流，下边界采用相应单河道最下游断面的水位流量关系，由曼宁公式计算；
Figure FDA0003798180270000111
式中：Q为流量；n为糙率；A为断面过水面积；R为断面水力半径；S为断面位置水面比降，以底坡代替；
步骤五：若存在引水/调水过程，则分析对流域水循环的影响若特定子流域内存在引水/调水过程，通过概化为该子流域出流过程的增加/减少，根据步骤四优化得到的参数及其他已确定参数，计算各子流域出流过程，分析引水/调水过程对流域水循环的影响。
CN202011463721.8A 2020-12-11 2020-12-11 一种基于自然子流域的通用流域水循环模拟计算方法 Active CN112651189B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011463721.8A CN112651189B (zh) 2020-12-11 2020-12-11 一种基于自然子流域的通用流域水循环模拟计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011463721.8A CN112651189B (zh) 2020-12-11 2020-12-11 一种基于自然子流域的通用流域水循环模拟计算方法

Publications (2)

Publication Number Publication Date
CN112651189A CN112651189A (zh) 2021-04-13
CN112651189B true CN112651189B (zh) 2023-03-10

Family

ID=75353856

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011463721.8A Active CN112651189B (zh) 2020-12-11 2020-12-11 一种基于自然子流域的通用流域水循环模拟计算方法

Country Status (1)

Country Link
CN (1) CN112651189B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113128067A (zh) * 2021-05-06 2021-07-16 大连理工大学 一种基于分布式时变地貌单位线的山丘区小流域洪水预报方法
CN115034160B (zh) * 2022-06-28 2023-04-21 河海大学 一种基于马斯京根法参数转换为等效河道的计算方法
CN115659686B (zh) * 2022-11-06 2023-11-21 长江水利委员会水文局 一种复杂河型断面法河道沿程槽蓄量批量计算方法
CN115618687B (zh) * 2022-11-09 2023-04-28 中国水利水电科学研究院 一种基于子流域和一维有限体积单元的降雨径流模拟方法
CN116090744B (zh) * 2022-12-02 2023-11-07 广东省水利水电科学研究院 山地丘陵区小型灌区灌溉用水配置方法、计算机装置及存储介质

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5557082B2 (ja) * 2009-02-25 2014-07-23 独立行政法人防災科学技術研究所 降水分布の推定システムおよび降水分布の推定方法
CN102034001A (zh) * 2010-12-16 2011-04-27 南京大学 一种以栅格为模拟单元的分布式水文模型设计方法
CN102567635B (zh) * 2011-12-23 2015-06-03 中国水利水电科学研究院 一种定量区分水循环演变过程中不同因素贡献的方法
CN104281780B (zh) * 2014-10-11 2016-03-23 水利部交通运输部国家能源局南京水利科学研究院 线性水库滞留汇流及嵌套流域（多子流域）汇流方法
CN108416049B (zh) * 2018-03-19 2020-07-17 河海大学 一种高寒山区流域雨雪混合产流计算方法
CN110570517B (zh) * 2019-08-07 2020-10-16 河海大学 一种基于下垫面特征的重配置产流模拟方法
CN111914432B (zh) * 2020-08-14 2022-11-29 贵州东方世纪科技股份有限公司 一种基于大数据的水文预报方法

Also Published As

Publication number Publication date
CN112651189A (zh) 2021-04-13

Similar Documents

Publication Publication Date Title
CN112651189B (zh) 一种基于自然子流域的通用流域水循环模拟计算方法
Abdulla et al. Application of a macroscale hydrologic model to estimate the water balance of the Arkansas‐Red River Basin
Dash et al. A SWAT-Copula based approach for monitoring and assessment of drought propagation in an irrigation command
Karvonen et al. A hydrological model for predicting runoff from different land use areas
Batelaan et al. GIS-based recharge estimation by coupling surface–subsurface water balances
Zhang et al. Predicting hydrologic response to climate change in the Luohe River basin using the SWAT model
CN102722909B (zh) 一种基于自适应尺度dem的流水线拓扑网络动态模拟方法
CN110928965B (zh) 一种基于流域精细分类的多模型灵活架构的模拟方法
CN113742910A (zh) 基于中小流域洪水预报的水库来水量预警预报方法及***
Wu et al. Reuse of return flows and its scale effect in irrigation systems based on modified SWAT model
Carpenter et al. Continuous streamflow simulation with the HRCDHM distributed hydrologic model
Abdullahi et al. Assessment of water availability in the Sokoto Rima River Basin
Ahmadzadeh et al. Assessment of water demand reliability using SWAT and RIBASIM models with respect to climate change and operational water projects
Klaas et al. Evaluating the impact of grid cell properties in spatial discretization of groundwater model for a tropical karst catchment in Rote Island, Indonesia
CN114091277B (zh) 一种考虑初始状态变量影响的新安江模型参数率定方法
Sang et al. Prediction of water resources change trend in the Three Gorges Reservoir Area under future climate change
CN114969655A (zh) 流域输沙量的模拟估算方法
Kaur et al. Impact of climate change on groundwater levels in Sirhind Canal Tract of Punjab, India
Bemporad et al. A distributed approach for sediment yield evaluation in Alpine regions
Fu et al. Water resource availability assessment through hydrological simulation under climate change in the Huangshui watershed of the qinghai–tibet plateau
Abraham et al. Simulation of chain of tanks to augment water supply: a case study from Tamil Nadu
Kuok et al. Stage-Storage and Flood Risk Assessments of Upgraded Batu Kitang Submersible Weir
Mirza Hydrologic modeling approaches for climate impact assessment in South Asia
CN114611290B (zh) 一种基于量变参数水文不确定性处理器的场次洪水水文模型实时预报方法
Wijaya et al. The Hydraulic Modelling of Capacity of Water Pool in Universitas Jambi

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant