Python实现新安江三水源水文模型代码与降雨流量过程图绘制指南

Tips：

        本期向大家分享三水源新安江模型代码的编写和最终的降雨流量过程图绘制，需要注意的是：

①为便于读者理解过程，将“蒸散发、产流、分水源、汇流”4个部分的代码按顺序编写；

②各阶段的参数在每个阶段开始时进行定义赋值，方便读者理解参数含义和用途；

③因为并没有划分子流域，所以也没有使用定义函数来调用，显得更直观清晰一些。

1 新安江模型简介         

        新安江模型按照三层蒸散发模式计算流域蒸散发，按蓄满产流概念计算降雨产生的总径流量，采用流域蓄水曲线考虑下垫面不均匀对产流面积变化的影响。在径流成分划分方面，对三水源情况，按“山坡水文学”产流理论用一个具有有限容积和测孔、底孔的自由水蓄水库把总径流划分成饱和地面径流、壤中水径流和地下水径流。在汇流计算方面，单元面积的地面径流汇流一般采用单位线法，壤中水径流和地下水径流的汇流则采用线性水库法。河网汇流一般采用分段连续演算的Muskingum法或滞时演算法。

        具体的新安江模型原理的学习主要参考了“贵仁模型技术手册”、“Ai尚研修相关文章”、水文预报教材和相关论文”。（给出相关链接如下）

2.1 新安江模型 · 贵仁模型技术手册 (keepsoft.net)

【水文模型】三水源新安江模型 (qq.com)

水文预报（第5版）_百度百科 (baidu.com)

2 三水源新安江模型代码

2.1 导入相关库

import os # 导入os模块，用于文件操作
from openpyxl import load_workbook # 导入openpyxl模块，用于读取Excel文件
import numpy as np # 导入numpy模块
import matplotlib.pyplot as plt # 导入绘图模块
from matplotlib import rcParams # 设置matplotlib使用的字体
from scipy.interpolate import make_interp_spline # 折线平滑处理的方法

2.2 读取Excel数据

2.2.1 数据展示

2.2.2 读取数据

filename = 'test01_zhengsanfa_chanliu.xlsx'
if os.path.exists(filename):
    wb = load_workbook(filename)
   
    ws = wb.active
    rows = ws.rows
  
    header = []
    for s_header in  next(rows):
        header.append(s_header.value)
   
    date = []
    for row in rows:
        temprow_dict = { }
        for headercell,datecell in zip(header,row):
            temprow_dict[headercell] = datecell.value
        date.append(temprow_dict)
    

def duqu(key):
    lst=[]
    for i in range (len(date)):
        dic_1=date[i]
        # print(dic_1)
        lst.append(dic_1[key])
    return lst
   
p=duqu('P')
em=duqu('EM')
t=duqu('t')

2.3 蒸散发＋产流阶段

## ---------------蒸散发（三层蒸散发模型）+产流--------------------
# 以下为初始参数设定：
WU=3 #昨天的
WL=15 #昨天的
WD=20 #昨天的
K=0.85
WUM=20
WLM=80
WDM=20
C=0.15  # 土壤深层蒸散发系数
WM=120  # 流域平均蓄水容量
IM=0.01
B=0.4
WMM=(1+B)*WM/(1-IM)

W=WU+WL+WD #就是今天的初始土壤含水量（昨天结束的土壤含水量）
lst1=[]  #存储R
lst2=[]  #EU
lst3=[]  #EL
lst4=[]  #ED
lst5=[]  #E
lst6=[]  #WU
lst7=[]  #WL
lst8=[]  #WD
lst9=[]  #W

for i in range(len(p)):
    P=p[i]
    EM=em[i]
    a = WMM * (1 - (1 - W / WM) ** (1 / (B + 1)))
    if P - K * EM <= 0:
        R = 0
    elif P - K * EM + a < WMM:
        R = P - K * EM + W - WM + WM * (1 - (P - K * EM + a) / WMM) ** (1 + B)
    else:
        R = P - K * EM - (WM - W)
    # 计算R的W都是前一天结束时的W
    if P + WU >= K * EM:
        EU = K * EM
        EL = 0
        ED = 0
        if P<=K*EM:
            WU=WU-(K*EM-P)

        elif P+WU-K*EM-R<WUM:
            WU=P+WU-K*EM-R

        elif WL+P+WU-K*EM-R-WUM<WLM:
            WL=WL+P+WU-K*EM-R-WUM
            WU = WUM

        elif WL+WD+P+WU-K*EM-R-WUM-WLM<WDM:
                                            # （所以要保证本行和下一行计算的公式一致，这样就不会错）
            WD=WL+WD+P+WU-K*EM-R-WUM-WLM
            WL = WLM
            WU = WUM

        else:
            WU = WUM
            WL = WLM
            WD = WDM
    else:
        EU = P + WU
        WU = 0
        if WL / WLM >= C:
            EL = (K * EM - EU) * WL / WLM
            ED = 0
            WL-=EL
        else:
            if WL >= C * (K * EM - EU):
                EL = C * (K * EM - EU)
                ED = 0
                WL -= EL
            else:
                EL = WL
                ED = C * (K * EM - EU) - WL
                WL=0
                if WD>ED:
                    WD-=ED
                else:
                    WD=0
    E = EU + EL + ED  # 蒸散发量

    W=WU+WL+WD #今天结束的土壤含水量（明天的初始土壤含水量）

    lst1.append(R)
    lst2.append(EU)
    lst3.append(EL)
    lst4.append(ED)
    lst5.append(E)
    lst6.append(WU)
    lst7.append(WL)
    lst8.append(WD)
    lst9.append(W)

补充：三层蒸散发模型代码案例验证

        将代码计算结果与1984年赵人俊教授所著的《流域水文模拟——新安江模型与陕北模型》P77页表5-3进行对比，经验证结果一致，说明代码正确可行。

《流域水文模拟——新安江模型与陕北模型》

2.4 分水源阶段

        本阶段的自由水蓄量S容易搞错，上一时段末的S是本时段初的S0，在计算过程中还要利用S和产流面积比例FR的关系再做变化。

## -------------------分水源（三水源） ------------------------
# 以下为初始参数设定：
IM=0.01  #不透水面积比例
EX=1.5  #自由水蓄水容量曲线指数
SM=36  #自由水蓄水容量平均值
KI=0.13  #壤中流的日出流系数
KG=0.13  #地下径流的日出流系数
MS=SM*(1+EX)

# 以下两参数初次运行可以使用，之后将其放入相应列表
S0=0  #本时段初的自由水蓄量
FR0=0.05  #上时段产流面积比例

lst_FR=[FR0]  #存储FR
lst_Si=[]  #Si=S0*FRO/FR
lst_AU=[]  #AU
lst_RS=[]  #RS
lst_S=[]  #S
lst_RI=[]  #RI
lst_RG=[]  #RG
lst_S1=[S0]  #存放下一阶段S0，即S1
for i in range(len(p)):
    P=p[i]
    EM=em[i]
    R=lst1[i] # 产流：由上一步得到
    S0=lst_S1[i]
    FR0=lst_FR[i]
    if P-K*EM<=0:
        FR=FR0

    else:
        FR = (R - IM * (P - K * EM)) / (P - K * EM)

    Si=S0*FR0/FR
    MS=SM*(1+EX)
    AU=MS*(1-(1-(S0*FR0/FR)/SM)**(1/(EX+1)))
    if AU+P-K*EM<=MS:
        RS=FR*(P-K*EM+S0*FR0/FR-SM+SM*(1-(P-K*EM+AU)/MS)**(EX+1))
    else:
        RS=FR*(P-K*EM+S0*FR0/FR-SM)
    if R==0:
        RS=0
    S=S0*FR0/FR+(R-RS)/FR
    RI=KI*S*FR
    RG=KG*S*FR
    S1=S*(1-KI-KG)

    lst_FR.append(FR)
    lst_AU.append(AU)
    lst_RS.append(RS)
    lst_S.append(S)
    lst_RI.append(RI)
    lst_RG.append(RG)
    lst_S1.append(S1)
lst_S1.pop(0)

2.5 汇流阶段

# # ------------------------------汇流阶段-------------------------------------------
# ----------（单元面积河网总入流：线性水库法+单元面积以下河道汇流：马斯京根法）------------------
rs=lst_RS
ri=lst_RI
rg=lst_RG
# 以下为初始参数设定：
CS=0.1  #地表径流消退系数
CI=0.3  #壤中流消退系数
CG=0.9 #地下径流消退系数
U=72  #单位转换系数
# 以下三个参数初次运行可以使用，之后将其放入相应列表
QRS0=0
QRI0=0
QRG0=0
QT0=QRS0+QRI0+QRG0
lst_QRS=[QRS0]
lst_QRI=[QRI0]
lst_QRG=[QRG0]
lst_QT=[QT0]

#线性水库法
for i in range(len(rs)):
    RS = rs[i]
    RI = ri[i]
    RG = rg[i]
    QRS0 = lst_QRS[i]
    QRI0 = lst_QRI[i]
    QRG0 = lst_QRG[i]
    QRS = QRS0 * CS + RS * (1 - CS) * U # 滞后演算法
    QRI = QRI0 * CI + RI * (1 - CI) * U
    QRG = QRG0 * CG + RG * (1 - CG) * U
    QT = QRS + QRI + QRG #单元面积河网总入流
    lst_QRS.append(QRS)
    lst_QRI.append(QRI)
    lst_QRG.append(QRG)
    lst_QT.append(QT)
dt=6 # 演算时段长度
K_m=6 # 槽蓄曲线坡度
x=0.4 # 流量比重系数
n = len(lst_QT)
O2 = np.zeros((n, 1))
c0 = (0.5 * dt - K_m * x) / (0.5 * dt + K_m - K_m * x)
c1 = (0.5 * dt + K_m * x) / (0.5 * dt + K_m - K_m * x)
c2 = (-0.5 * dt + K_m - K_m * x) / (0.5 * dt + K_m - K_m * x)

# 马斯京根法
for i in range(1, n):
    I1 = lst_QT[i - 1]
    I2 = lst_QT[i]
    Q1 = O2[i - 1]
    O2[i] = c0 * I2 + c1 * I1 + c2 * Q1
O2=O2[1:] #去掉第1个元素，第1个元素是为了用于流量演算，实际中不存在。
Qout = np.array(O2) # 汇流后的流量，转化为数组形式
np.set_printoptions(suppress=True, precision=4) #将演算流量按照四位小数显示，不以科学计数法显示

3 降雨流量过程图代码

##-----------------绘图--------------------------------
time=t
flow=Qout
precipitation=p

# 设置matplotlib使用的字体
rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
rcParams['axes.unicode_minus'] = False  # 正常显示负号
# 创建图形
fig, ax1 = plt.subplots()
# 演算流量插值平滑处理
time_smooth = np.linspace(np.array(time).min(), np.array(time).max(), 300)  # np.linspace 等差数列,从x.min()到x.max()生成300个数，便于后续插值
flow_smooth = make_interp_spline(time, flow)(time_smooth)

# 创建第一个y轴并绘制流量图（左侧y轴）
ax1.set_xlabel('时间 (h)')
ax1.set_ylabel('流量 (m3/s)', color='black') # 设置y轴标签
ax1.plot(time_smooth, flow_smooth, color='#FA7F6F',linewidth=3,label='演算流量') # 绘制流量过程线
ax1.tick_params(axis='y', labelcolor='black') # 设置y轴颜色

# 创建第二个y轴并绘制倒置的降水图（右侧y轴）
width_close = dt  # 设置柱状图宽度为演算时段长度dt=6h，使柱子之间没有空隙
ax2 = ax1.twinx()
ax2.set_ylabel('降雨 (mm)', color='black')
# 将降雨量的时间轴调整为演算时段的中心，降雨图像平移时段的一半，否则降雨柱状图的中心会正对着时段开始时刻，不整齐不合理。
ax2.bar([t+0.5*dt for t in time], -np.array(precipitation), color='#8ECFC9', alpha=0.3, width=width_close,edgecolor='black',label='降雨量')  # 使用负值来倒置柱状图
ax2.tick_params(axis='y', labelcolor='black')
ax2.set_ylim(bottom=-max(precipitation), top=0)  # 调整y轴范围，使降雨柱状图在上方显示

# 将y轴标签设置为正数
y_ticks = np.linspace(0, -max(precipitation)*3, num=10)
ax2.set_yticks(y_ticks)
ax2.set_yticklabels([f'{int(-tick)}' for tick in y_ticks])  # 设置为正数

# 调整坐标轴范围,将左边的和底部的坐标轴移到原点位置
ax1.set_xlim(0, max(time))  # 设置 x 轴范围从 0 到最大时间值
ax1.set_ylim(0, max(flow)*1.5)  # 设置第一个纵轴范围从 0 到最大降水量值

# 图例,在此之前，要在plot函数中加入label参数
ax1.legend(loc='upper right',bbox_to_anchor=(1, 0.9))
ax2.legend(loc='upper right')

# 自动调整布局＋图形显示
fig.tight_layout()
plt.show()

4 结果展示

作者：摆烂老大