Python 水文预报三水源新安江模型构建

这篇具有很好参考价值的文章主要介绍了Python 水文预报三水源新安江模型构建。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

1 前言

新安江模型是极为重要的水文模型,此拙作只作为简单的引导,希望大家在阅读的时候能够有自己的思考,例如引入智能进化算法率定参数、程序自动绘制精美的图表、程序自动输出评价指标、分析模型对参数的敏感程序、思考模型产生误差的原因、实现改进的新安江模型、提升自己的编程能力。

2 新安江模型介绍

新安江模型是华东水利学院(现河海大学)赵人俊教授团队提出的一个水文模型,是中国少有的一个具有世界影响力的水文模型。新安江模型模型适用于湿润地区和半湿润地区。它将整个流域划分为多个块状单元流域,并对每个单元流域进行产汇流计算,得出单元流域的出口流量过程。接着进行河道洪水演算,求得流域出口的流量过程。将每个单元流域的出流过程相加,即可得到流域的总出流过程。

2.1 模型结构

为了考虑降水和流域下垫面分布不均匀的影响,新安江模型的结构设计为分散性的,分为蒸散发计算产流计算分水源计算汇流计算四个层次结构。每块单元流域的计算流程如下图所示。
新安江模型python代码,Hydrology,python

2.1.1 蒸散发计算

在新安江模型中,流域蒸散发计算没有考虑流域内土壤含水量在面上分布的不均匀性,而是按土壤垂向分布的不均匀性将土层分为三层,用三层蒸散发模型计算蒸散发量,计算公式如下:
新安江模型python代码,Hydrology,python

2.1.2 产流计算

产流计算采用蓄满产流机制。蓄满是指包气带的土壤含水量达到田间持水量。蓄满产流是指:降水在满足田间持水量以前不产流,所有的降水都被土壤所吸收而成为张力水;降水在满足田间持水量以后,所有的降水扣除蒸发量都产流。
按照蓄满产流的概念,采用蓄水容量-面积分配曲线来考虑土壤缺水量分布不均匀的问题,以此来计算产流量 R。

  • 计算流域平均的初始土壤含水量对应的纵坐标 a:

新安江模型python代码,Hydrology,python

  • 如果 P - E > 0 ,计算径流量 R:

新安江模型python代码,Hydrology,python

2.1.3 三水源划分

根据径流实验观测和径流形成原理,将径流划分为地面径流 RS、壤中流 RI、地下径流 RG。与蓄满产流模型相类似,由于下垫面的不均匀性,自由水蓄量也存在空间分布不均匀性,因此用自由水蓄水容量曲线来考虑这点,以此来进行三水源划分。

  • 计算流域平均的初始自由水蓄量对应的纵坐标 AU:

新安江模型python代码,Hydrology,python

  • 计算地面径流 RS:

新安江模型python代码,Hydrology,python

  • 本时段的自由水蓄量 S:

新安江模型python代码,Hydrology,python

  • 相应的壤中流 RI 和地下径流 RG:

新安江模型python代码,Hydrology,python
新安江模型python代码,Hydrology,python

  • 下一时段初的自由水蓄量 S1:

新安江模型python代码,Hydrology,python

2.1.4 三水源汇流计算

  • 地表径流汇流,采用线性水库法:
    新安江模型python代码,Hydrology,python
  • 壤中流汇流,表层自由水侧向流动,出流后成为表层壤中流进入河网。若土层较厚,表层自由水还可以渗入到深层土,经过深层土的调蓄作用才进入河网。壤中流汇流也可采用线性水库法:
    新安江模型python代码,Hydrology,python
  • 地下径流汇流,采用线性水库法:
    新安江模型python代码,Hydrology,python
  • 单元面积河网总入流,单元面积河网总入流 QT 为地表径流量、壤中流与地下径流量之和:
    新安江模型python代码,Hydrology,python
  • 单元面积河网汇流,单位面积的河流采用滞后演算法:
    新安江模型python代码,Hydrology,python
  • 河道汇流,河道汇流采用马斯京根分段流量演算法,将演算河段分为 n 个单元河段,用马斯京根法连续进行 n 次演算,以求得出流过程:
    新安江模型python代码,Hydrology,python

2.2 模型参数

新安江模型中各参数意义如下(还包括初始土壤含水量):

参数 参数物理意义
K 蒸散发折算系数
B 流域蓄水容量分布曲线指数
C 深层散发系数
WM 张力水容量
WUM 上层张力水容量
WLM 下层张力水容量
IM 不透水面积比例
SM 自由水容量
EX 流域自由水容量分布曲线指数
KG 地下水出流系数
KI 壤中流出流系数
CG 地下水消退系数
CI 壤中流消退系数
CS 河网水流消退系数
L 河网汇流滞时
X 马斯京根法参数

3 模型代码

新安江模型有两种时间尺度:日模型(24 h)、次洪模型(1 h)。
日模型的结果可为次洪模型提供初始值,这里只给出次洪模型相应阶段的代码。如果汇流时间小于24 h,日模型可以简化不考虑河道汇流。

3.1 代码说明

参考:Python 建立流域三层蒸发和蓄满产流模型(二水源划分)
需要用到pandas库和math库,不会的自行CSDN。

3.2 更新土壤含水量

初次计算需要设置初始土壤含水量,后每次计算需要先更新土壤含水量。

def calculate_WU_WL_WD(df, i):  # 更新土壤含水量
    if i == 0:
        # 设置初始值
        df.loc[i, ["WU"]] = WUM_init
        df.loc[i, ["WL"]] = WLM_init
        df.loc[i, ["WD"]] = WDM_init
    else:
        infi = df["PE"][i - 1] - df["R"][i - 1]  # 上层蓄水量不为0,则只存在EU;且上层蓄水量与下渗的净降雨有关
        df.loc[i, ["WU"]] = df["WU"][i - 1] + infi
        df.loc[i, ["WL"]] = df["WL"][i - 1]  # 假设扣除损失量后上层蓄水量仍大于或等于0,则中、深层蓄水量不会变化
        df.loc[i, ["WD"]] = df["WD"][i - 1]
        if df["WU"][i] < 0:  # 若上层蓄水量小于损失量
            df.loc[i, ["WL"]] = df["WL"][i - 1] + df["WU"][i]  # 则蒸发损失认为发生中层
            df.loc[i, ["WU"]] = 0
            df.loc[i, ["WD"]] = df["WD"][i - 1]  # 假设扣除损失量后中层蓄水量仍大于或等于0,则深层蓄水量不会变化
            if df["WL"][i] < 0:  # 若中层蓄水量小于损失量
                df.loc[i, ["WD"]] = df["WD"][i - 1] + df["WL"][i]  # 则蒸发损失认为发生深层
                df.loc[i, ["WL"]] = 0
                df.loc[i, ["WU"]] = 0
                if df["WD"][i] < 0:  # 若深层蓄水量小于损失量
                    df.loc[i, ["WD"]] = 0
                    df.loc[i, ["WL"]] = 0
                    df.loc[i, ["WU"]] = 0

        # 设置蓄水量上限,不应超过蓄水容量
        if df["WU"][i] > WUM:  # 若扣除损失后入渗量足以补充上层蓄水容量,则多余的会逐层向下补充
            df.loc[i, ["WL"]] = df["WU"][i] - WUM + df["WL"][i - 1]  # 首先补充中层
            df.loc[i, ["WU"]] = WUM
            df.loc[i, ["WD"]] = df["WD"][i - 1]
            if df["WL"][i] > WLM:  # 若扣除损失后入渗量足以补充中层蓄水容量,则多余的会逐层向下补充
                df.loc[i, ["WD"]] = df["WL"][i] - WLM + df["WD"][i - 1]  # 补充深层
                df.loc[i, ["WU"]] = WUM
                df.loc[i, ["WL"]] = WLM
                if df["WD"][i] > WDM:  # 完全补充
                    df.loc[i, ["WD"]] = WDM
                    df.loc[i, ["WU"]] = WUM
                    df.loc[i, ["WL"]] = WLM

3.3 三层蒸发模型计算

四个条件判断,需要先计算流域蒸发能力。

def calculate_EU_EL_ED(df, i):  # 三层蒸发模型计算
    df.loc[i, ["EP"]] = df["E0"][i] * K  # 计算流域蒸发能力

    """四个条件判断"""
    if df["WU"][i] + df["P"][i] >= df["EP"][i]:
        df.loc[i, ["EU"]] = df["EP"][i]
        df.loc[i, ["EL"]] = 0
        df.loc[i, ["ED"]] = 0
    if df["WU"][i] + df["P"][i] < df["EP"][i] and df["WL"][i] >= C * WLM:
        df.loc[i, ["EU"]] = df["WU"][i] + df["P"][i]
        df.loc[i, ["EL"]] = (df["EP"][i] - df["EU"][i]) * (df["WL"][i] / WLM)
        df.loc[i, ["ED"]] = 0
    if df["WU"][i] + df["P"][i] < df["EP"][i] and C * (df["EP"][i] - df["EU"][i]) <= df["WL"][i] and df["WL"][
        i] < C * WLM:
        df.loc[i, ["EU"]] = df["WU"][i] + df["P"][i]
        df.loc[i, ["EL"]] = C * (df["EP"][i] - df["EU"][i])
        df.loc[i, ["ED"]] = 0
    if df["WU"][i] + df["P"][i] < df["EP"][i] and df["WL"][i] < C * (df["EP"][i] - df["EU"][i]):
        df.loc[i, ["EU"]] = df["WU"][i] + df["P"][i]
        df.loc[i, ["EL"]] = df["WL"][i]
        df.loc[i, ["ED"]] = C * (df["EP"][i] - df["EU"][i]) - df["EL"][i]

3.4 蒸发量、净降水量、土壤总蓄水量计算

需要计算蒸发量E、净降水量PE、土壤总蓄水量计算W,为后续计算产流量R做准备。

def calculate_E_PE_W(df, i):
    df.loc[i, ["E"]] = df["EU"][i] + df["EL"][i] + df["ED"][i]  # 计算总的蒸发量
    df.loc[i, ["PE"]] = df["P"][i] - df["E"][i]  # 计算降水扣除蒸发量
    df.loc[i, ["W"]] = df["WU"][i] + df["WL"][i] + df["WD"][i]  # 计算土壤总蓄水量

3.5 蓄满产流模型计算产流

采用蓄满产流模型,认为当PE > 0时才有可能产生径流R,否则不产生。

def calculate_R(df, i):  # 蓄满产流模型计算产流
    """计算公式"""
    if df["PE"][i] > 0:  # PE>0才有可能产生径流
        a = WMM * (1 - math.pow((1 - df["W"][i] / WM), (1 / (1 + b))))
        if a + df["PE"][i] <= WMM:
            df.loc[i, ["R"]] = df["PE"][i] + df["W"][i] - WM + WM * math.pow((1 - (df["PE"][i] + a) / WMM), (b + 1))
        else:
            df.loc[i, ["R"]] = df["PE"][i] - (WM - df["W"][i])
    else:
        df.loc[i, ["R"]] = 0

3.6 流域产流面积计算

根据蓄水容量-面积分配曲线计算流域产流面积。

def calculate_FR(df, i):  # 计算流域产流面积
    """三水源划分产流面积迭代"""
    if df["R"][i] > 0:
        df.loc[i, "FR"] = df["R"][i] / df["PE"][i]
        if df["FR"][i] > 1:
            df.loc[i, "FR"] = 1
    else:
        if i == 0:
            df.loc[i, "FR"] = FR1
        else:
            df.loc[i, "FR"] = df["FR"][i - 1]

3.7 三水源划分

将产流R划分为RS、RG、RI,可参照三水源划分公式,首个时段设置初始的自由水蓄量。

def calculate_RS_RG_RI(df, i):  # 三水源划分判断函数
    """计算公式"""
    if i == 0:  # 第一时段计算
        df.loc[i, "S1"] = S1
        if df["PE"][i] > 0:
            AU = SMM * (1 - math.pow((1 - (((df["S1"][i] * FR1) / df["FR"][i]) / Sm)), 1 / (1 + EX)))
            if df["PE"][i] + AU < SMM:
                df.loc[i, "RS"] = df["FR"][i] * (df["PE"][i] + (df["S1"][i] * FR1) / df["FR"][i] - Sm +
                                                 Sm * math.pow((1 - (df["PE"][i] + AU) / SMM), 1 + EX))
            if df["PE"][i] + AU >= SMM:
                df.loc[i, "RS"] = df["FR"][i] * (df["PE"][i] + (df["S1"][i] * FR1) / df["FR"][i] - Sm)
            S = (df["S1"][i] * FR1) / df["FR"][i] + (df["R"][i] - df["RS"][i]) / df["FR"][i]
            df.loc[i, "RI"] = KI * S * df["FR"][i]
            df.loc[i, "RG"] = KG * S * df["FR"][i]
            df.loc[i + 1, "S1"] = S * (1 - KI - KG)
        else:
            S = (df["S1"][i] * FR1) / df["FR"][i]
            df.loc[i + 1, "S1"] = S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KI * S * df["FR"][i]
            df.loc[i, "RI"] = KG * S * df["FR"][i]
    else:  # 其余时段计算
        if df["PE"][i] > 0:
            AU = SMM * (1 - math.pow((1 - (((df["S1"][i] * df["FR"][i - 1]) / df["FR"][i]) / Sm)),
                                     1 / (1 + EX)))
            if df["PE"][i] + AU < SMM:
                df.loc[i, "RS"] = df["FR"][i] * (
                        df["PE"][i] + (df["S1"][i] * df["FR"][i - 1]) / df["FR"][i] - Sm +
                        Sm * math.pow((1 - (df["PE"][i] + AU) / SMM), 1 + EX))
            if df["PE"][i] + AU >= SMM:
                df.loc[i, "RS"] = df["FR"][i] * (
                        df.loc[i, "PE"] + (df["S1"][i] * df["FR"][i - 1]) / df["FR"][i] - Sm)
            S = (df["S1"][i] * df["FR"][i - 1]) / df["FR"][i] + (df["R"][i] - df["RS"][i]) / df["FR"][i]
            df.loc[i, "RI"] = KI * S * df["FR"][i]
            df.loc[i, "RG"] = KG * S * df["FR"][i]
            df.loc[i + 1, "S1"] = S * (1 - KI - KG)
        else:
            S = (df["S1"][i] * df["FR"][i - 1]) / df["FR"][i]
            df.loc[i + 1, "S1"] = S * (1 - KG - KI)
            df.loc[i, "RS"] = 0
            df.loc[i, "RG"] = KG * S * df["FR"][i]
            df.loc[i, "RI"] = KI * S * df["FR"][i]

3.8 三水源汇流计算

根据线性水库、滞后演算公式,对RS、RG、RI三种径流成分进行汇流计算。

def calculate_Q(df, i, U):  # 三水源汇流计算,单元面积河网汇流
    """计算公式"""
    df.loc[i, "QS"] = df["RS"][i] * U  # 地面径流的坡地汇流
    if i == 0:  # 设置初始壤中流、地下径流
        df.loc[i, "QI"] = 1 / 3 * Q
        df.loc[i, "QG"] = 1 / 3 * Q
    else:
        df.loc[i, "QI"] = CI * df["QI"][i - 1] + (1 - CI) * df["RI"][i] * U  # 壤中流汇流,消退系数CI
        df.loc[i, "QG"] = CG * df["QG"][i - 1] + (1 - CG) * df["RG"][i] * U  # 地下径流汇流,消退系数CG

    df.loc[i, "QT"] = df["QS"][i] + df["QI"][i] + df["QG"][i]  # 单元面积河网汇流,滞后演算法

    if 0 <= i <= L:
        df.loc[i, "Qt"] = Q
    elif i > L:
        df.loc[i, "Qt"] = df["Qt"][i - 1] * CS + (1 - CS) * df["QT"][i - L]

3.9 河道汇流演算

采用马斯京根法进行河道汇流演算,需要提前计算好子河段数量n。

def calculate_Qi(df, i, n):  # 河道汇流演算,马斯京根法
    df.loc[i, "Q1"] = df["Qt"][i]
    if n > 0:
        """单元河段参数"""
        K_l = T
        x_l = 0.5 - n * (1 - 2 * X) / 2
        """计算C0、C1、C2"""
        C0 = (0.5 * T - K_l * x_l) / (0.5 * T + K_l - K_l * x_l)
        C1 = (0.5 * T + K_l * x_l) / (0.5 * T + K_l - K_l * x_l)
        C2 = 1 - C0 - C1
        # df.loc[i, "Q1"] = df["Qt"][i]
        for p in range(n):
            """进行分段计算"""
            I2 = df[f"Q{1 + p}"][i]
            if i == 0:
                I1 = Q
                Q1 = Q
            else:
                I1 = df[f"Q{1 + p}"][i - 1]
                Q1 = df[f"Q{2 + p}"][i - 1]
            df.loc[i, f"Q{2 + p}"] = C0 * I2 + C1 * I1 + C2 * Q1

3.10 计算顺序

需要注意的是,这里是对每个子流域分别(循环)计算,每个sheet表对应一个子流域的降水,最后需要加在一起才是流域出口断面流量。

if __name__ == '__main__':
   	...
   	忽略数据读取代码和保存至excel代码
   	...
    for k in range(sheet_count):
        print(f"第{k + 1}次运行ing")
        df = read_data(k) # 这个是自编的函数
        for i in range(len(df)):  # 循环
            calculate_WU_WL_WD(df, i)
            calculate_EU_EL_ED(df, i)
            calculate_E_PE_W(df, i)
            calculate_R(df, i)
            calculate_FR(df, i)
            calculate_RS_RG_RI(df, i)
            U = Area[k] / (3.6 * T)  # F单位为km^2,T单位为h
            calculate_Q(df, i, U)
            calculate_Qi(df, i, N[k])
       ...
   	忽略保存至excel代码
   	...

这里没有给出完整的读取数据、处理数据、绘制图表代码,需要根据自己的数据来写,不断思考和寻找问题答案的过程才是编程!。

以下是用于计算的Excel格式
新安江模型python代码,Hydrology,python

4 总结

新安江模型是极为重要的水文模型,此拙作只作为简单的引导,希望大家在阅读的时候能够有自己的思考。

如有问题,欢迎探讨,邮箱:mymail_zo@qq.com文章来源地址https://www.toymoban.com/news/detail-861148.html

到了这里,关于Python 水文预报三水源新安江模型构建的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 【SWAT水文模型】SwatWeather软件使用教程

    当气象数据部分缺失时,SWAT模型构建的天气发生器可对缺测数据进行插补和完善。天气发生器参变量体系的计算主要依靠流域内部及临近站点的观测数据,经长序列分析,得到各关键水文气象要素的月平均值及标准差/偏态系数,包括最低气温、最高气温、降水量、干日日数

    2024年02月03日
    浏览(51)
  • 实操:用Flutter构建一个简单的微信天气预报小程序

    ​ 微信小程序是一种快速、高效的开发方式,Flutter则是一款强大的跨平台开发框架。结合二者,可以轻松地开发出功能丰富、用户体验良好的微信小程序。 这里将介绍如何使用Flutter开发一个简单的天气预报小程序,并提供相应的代码示例。 1. 准备工作 在开始之前,确保你

    2024年02月16日
    浏览(41)
  • 水文水动力模型在城市内涝、城市排水、海绵城市规划设计中深度应用

    随着计算机的广泛应用和各类模型软件的发展,将排水系统模型作为城市洪灾评价与防治的技术手段已经成为防洪防灾的重要技术途径。将聚焦于综合利用GIS及CAD等工具高效地进行大规模城市排水系统水力模型的建立,利用SWMM实现排水系统水力模拟。讲解SWMM深度应用方法,

    2024年02月15日
    浏览(37)
  • R语言的水文、水环境模型优化技术及快速率定方法与多模型案例实践

    在水利、环境、生态、机械以及航天等领域中,数学模型已经成为一种常用的技术手段。同时,为了提高模型的性能,减小模型误用带来的风险;模型的优化技术也被广泛用于模型的使用过程。模型参数的快速优化技术不但涉及到优化本身而且涉及采样设计、不确定性估计等

    2024年02月16日
    浏览(34)
  • Flutter开发微信小程序实战:构建一个简单的天气预报小程序

    微信小程序是一种快速、高效的开发方式,Flutter则是一款强大的跨平台开发框架。结合二者,可以轻松地开发出功能丰富、用户体验良好的微信小程序。 这里将介绍如何使用Flutter开发一个简单的天气预报小程序,并提供相应的代码示例。 在开始之前,确保你已经安装了Fl

    2024年02月12日
    浏览(52)
  • 水文模型有哪些?SWAT模型、VIC模型、HEC模型、HSPF模型、HYPE模型、SWMM模型、FVCOM模型、Delft3D模型等应用

    目录 ㈠ 从小白到精通SWAT模型学习建模方法、实例应用、高级进阶 ㈡ R+VIC模型融合实践技术应用及未来气候变化模型预测 ㈢ HEC-RAS一维、二维建模方法及实践技术应用 ㈣ HEC-HMS水文模型实践技术应用 ㈤ HSPF 模型应用 ㈥ HYPE分布式水文模型建模方法与案例分析 ㈦ 基于SWMM及

    2023年04月08日
    浏览(59)
  • Python机器学习、深度学习技术提升气象、海洋、水文领域实践应用

    Python是功能强大、免费、开源,实现面向对象的编程语言,在数据处理、科学计算、数学建模、数据挖掘和数据可视化方面具备优异的性能,这些优势使得Python在气象、海洋、地理、气候、水文和生态等地学领域的科研和工程项目中得到广泛应用。可以预见未来Python将成为的

    2023年04月23日
    浏览(49)
  • 基于Python机器学习、深度学习技术提升气象、海洋、水文领域实践应用能力

    查看原文 基于Python机器学习、深度学习技术提升气象、海洋、水文领域实践应用能力 目录 专题一、Python软件的安装及入门 专题二、气象常用科学计算库 专题三、气象海洋常用可视化库 专题四、爬虫和气象海洋数据 专题五、气象海洋常用插值方法 专题六、机器学习基础理

    2023年04月21日
    浏览(62)
  • Python| 水文 |批量读取NC文件降水数据并导出为Excel相应格式

    ①首先需要在网上下载相应whl的安装包 (注意:一定要保证版本与Python解释器版本一致,比如我的解释器是3.11版本,电脑是windows64位的系统,那么就选择netCDF4-1.6.5-cp311-cp311-win_amd64.whl进行下载) 在网上看到过一些下载地址,但是好多没法用,这里提供一个,不需要梯子:

    2024年04月28日
    浏览(37)
  • Python构建简单线性回归模型教程

    本文介绍如何构建简单线性回归模型及计算其准确率,最后介绍如何持久化模型。 线性回归表示发现函数使用线性组合表示输入变量。简单线性回归很容易理解,使用了基本的回归技术,一旦理解了这些基本概念,可以更好地学习其他类型的回归模型。 回归用于发现输入变

    2023年04月21日
    浏览(39)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包