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

方诗元

Python水文预报三水源新安江模型构建1前言2新安江模型介绍2.1模型结构2.1.1蒸散发计算2.1.2产流计算2.1.3三水源划分2.1.4三水源汇流计算2.2模型参数3模型代码3.1代码说明3.2更新土壤含水量3.3三层蒸发模型计算3.4蒸发量、净降水量、土壤总蓄水量计算3.5蓄满产流模型计算产流3.6流域产流面积计算3.7三水源划分3.8三水源汇流计算3.9河道汇流演算3.10计算顺序4总结1前言新安江模型是极为重要的水文模型，此拙作只作为简单的引导，希望大家在阅读的时候能够有自己的思考，例如引入智能进化算法率定参数、程序自动绘制精美的图表、程序自动输出评价指标、分析模型对参数的敏感程序、思考模型产生误差的原因、实现改进的新安江模型、提升自己的编程能力。2新安江模型介绍新安江模型是华东水利学院(现河海大学)赵人俊教授团队提出的一个水文模型，是中国少有的一个具有世界影响力的水文模型。新安江模型模型适用于湿润地区和半湿润地区。它将整个流域划分为多个块状单元流域，并对每个单元流域进行产汇流计算，得出单元流域的出口流量过程。接着进行河道洪水演算，求得流域出口的流量过程。将每个单元流域的出流过程相加，即可得到流域的总出流过程。2.1模型结构为了考虑降水和流域下垫面分布不均匀的影响，新安江模型的结构设计为分散性的，分为蒸散发计算、产流计算、分水源计算和汇流计算四个层次结构。每块单元流域的计算流程如下图所示。2.1.1蒸散发计算在新安江模型中，流域蒸散发计算没有考虑流域内土壤含水量在面上分布的不均匀性，而是按土壤垂向分布的不均匀性将土层分为三层，用三层蒸散发模型计算蒸散发量，计算公式如下：2.1.2产流计算产流计算采用蓄满产流机制。蓄满是指包气带的土壤含水量达到田间持水量。蓄满产流是指：降水在满足田间持水量以前不产流，所有的降水都被土壤所吸收而成为张力水；降水在满足田间持水量以后，所有的降水扣除蒸发量都产流。按照蓄满产流的概念，采用蓄水容量-面积分配曲线来考虑土壤缺水量分布不均匀的问题，以此来计算产流量R。计算流域平均的初始土壤含水量对应的纵坐标a：如果P-E>0，计算径流量R：2.1.3三水源划分根据径流实验观测和径流形成原理，将径流划分为地面径流RS、壤中流RI、地下径流RG。与蓄满产流模型相类似，由于下垫面的不均匀性，自由水蓄量也存在空间分布不均匀性，因此用自由水蓄水容量曲线来考虑这点，以此来进行三水源划分。计算流域平均的初始自由水蓄量对应的纵坐标AU：计算地面径流RS：本时段的自由水蓄量S：相应的壤中流RI和地下径流RG：下一时段初的自由水蓄量S1：2.1.4三水源汇流计算地表径流汇流，采用线性水库法：壤中流汇流，表层自由水侧向流动，出流后成为表层壤中流进入河网。若土层较厚，表层自由水还可以渗入到深层土，经过深层土的调蓄作用才进入河网。壤中流汇流也可采用线性水库法：地下径流汇流，采用线性水库法：单元面积河网总入流，单元面积河网总入流QT为地表径流量、壤中流与地下径流量之和：单元面积河网汇流，单位面积的河流采用滞后演算法：河道汇流，河道汇流采用马斯京根分段流量演算法，将演算河段分为n个单元河段，用马斯京根法连续进行n次演算，以求得出流过程：2.2模型参数新安江模型中各参数意义如下（还包括初始土壤含水量）：参数参数物理意义K蒸散发折算系数B流域蓄水容量分布曲线指数C深层散发系数WM张力水容量WUM上层张力水容量WLM下层张力水容量IM不透水面积比例SM自由水容量EX流域自由水容量分布曲线指数KG地下水出流系数KI壤中流出流系数CG地下水消退系数CI壤中流消退系数CS河网水流消退系数L河网汇流滞时X马斯京根法参数3模型代码新安江模型有两种时间尺度：日模型（24h）、次洪模型（1h）。日模型的结果可为次洪模型提供初始值，这里只给出次洪模型相应阶段的代码。如果汇流时间小于24h，日模型可以简化不考虑河道汇流。3.1代码说明参考：Python建立流域三层蒸发和蓄满产流模型(二水源划分)需要用到pandas库和math库，不会的自行CSDN。3.2更新土壤含水量初次计算需要设置初始土壤含水量，后每次计算需要先更新土壤含水量。defcalculate_WU_WL_WD(df,i):#更新土壤含水量ifi==0:#设置初始值df.loc[i,["WU"]]=WUM_initdf.loc[i,["WL"]]=WLM_initdf.loc[i,["WD"]]=WDM_initelse:infi=df["PE"][i-1]-df["R"][i-1]#上层蓄水量不为0，则只存在EU；且上层蓄水量与下渗的净降雨有关df.loc[i,["WU"]]=df["WU"][i-1]+infidf.loc[i,["WL"]]=df["WL"][i-1]#假设扣除损失量后上层蓄水量仍大于或等于0，则中、深层蓄水量不会变化df.loc[i,["WD"]]=df["WD"][i-1]ifdf["WU"][i]WUM:#若扣除损失后入渗量足以补充上层蓄水容量，则多余的会逐层向下补充df.loc[i,["WL"]]=df["WU"][i]-WUM+df["WL"][i-1]#首先补充中层df.loc[i,["WU"]]=WUMdf.loc[i,["WD"]]=df["WD"][i-1]ifdf["WL"][i]>WLM:#若扣除损失后入渗量足以补充中层蓄水容量，则多余的会逐层向下补充df.loc[i,["WD"]]=df["WL"][i]-WLM+df["WD"][i-1]#补充深层df.loc[i,["WU"]]=WUMdf.loc[i,["WL"]]=WLMifdf["WD"][i]>WDM:#完全补充df.loc[i,["WD"]]=WDMdf.loc[i,["WU"]]=WUMdf.loc[i,["WL"]]=WLM12345678910111213141516171819202122232425262728293031323334353637'运行运行3.3三层蒸发模型计算四个条件判断，需要先计算流域蒸发能力。defcalculate_EU_EL_ED(df,i):#三层蒸发模型计算df.loc[i,["EP"]]=df["E0"][i]*K#计算流域蒸发能力"""四个条件判断"""ifdf["WU"][i]+df["P"][i]>=df["EP"][i]:df.loc[i,["EU"]]=df["EP"][i]df.loc[i,["EL"]]=0df.loc[i,["ED"]]=0ifdf["WU"][i]+df["P"][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"]]=0ifdf["WU"][i]+df["P"][i]0时才有可能产生径流R，否则不产生。defcalculate_R(df,i):#蓄满产流模型计算产流"""计算公式"""ifdf["PE"][i]>0:#PE>0才有可能产生径流a=WMM*(1-math.pow((1-df["W"][i]/WM),(1/(1+b))))ifa+df["PE"][i]0:df.loc[i,"FR"]=df["R"][i]/df["PE"][i]ifdf["FR"][i]>1:df.loc[i,"FR"]=1else:ifi==0:df.loc[i,"FR"]=FR1else:df.loc[i,"FR"]=df["FR"][i-1]1234567891011'运行运行3.7三水源划分将产流R划分为RS、RG、RI，可参照三水源划分公式，首个时段设置初始的自由水蓄量。defcalculate_RS_RG_RI(df,i):#三水源划分判断函数"""计算公式"""ifi==0:#第一时段计算df.loc[i,"S1"]=S1ifdf["PE"][i]>0:AU=SMM*(1-math.pow((1-(((df["S1"][i]*FR1)/df["FR"][i])/Sm)),1/(1+EX)))ifdf["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"]=0df.loc[i,"RG"]=KI*S*df["FR"][i]df.loc[i,"RI"]=KG*S*df["FR"][i]else:#其余时段计算ifdf["PE"][i]>0:AU=SMM*(1-math.pow((1-(((df["S1"][i]*df["FR"][i-1])/df["FR"][i])/Sm)),1/(1+EX)))ifdf["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"]=0df.loc[i,"RG"]=KG*S*df["FR"][i]df.loc[i,"RI"]=KI*S*df["FR"][i]123456789101112131415161718192021222324252627282930313233343536373839404142'运行运行3.8三水源汇流计算根据线性水库、滞后演算公式，对RS、RG、RI三种径流成分进行汇流计算。defcalculate_Q(df,i,U):#三水源汇流计算，单元面积河网汇流"""计算公式"""df.loc[i,"QS"]=df["RS"][i]*U#地面径流的坡地汇流ifi==0:#设置初始壤中流、地下径流df.loc[i,"QI"]=1/3*Qdf.loc[i,"QG"]=1/3*Qelse:df.loc[i,"QI"]=CI*df["QI"][i-1]+(1-CI)*df["RI"][i]*U#壤中流汇流，消退系数CIdf.loc[i,"QG"]=CG*df["QG"][i-1]+(1-CG)*df["RG"][i]*U#地下径流汇流，消退系数CGdf.loc[i,"QT"]=df["QS"][i]+df["QI"][i]+df["QG"][i]#单元面积河网汇流，滞后演算法if0L:df.loc[i,"Qt"]=df["Qt"][i-1]*CS+(1-CS)*df["QT"][i-L]12345678910111213141516'运行运行3.9河道汇流演算采用马斯京根法进行河道汇流演算，需要提前计算好子河段数量n。defcalculate_Qi(df,i,n):#河道汇流演算，马斯京根法df.loc[i,"Q1"]=df["Qt"][i]ifn>0:"""单元河段参数"""K_l=Tx_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]forpinrange(n):"""进行分段计算"""I2=df[f"Q{1+p}"][i]ifi==0:I1=QQ1=Qelse: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*Q1123456789101112131415161718192021'运行运行3.10计算顺序需要注意的是，这里是对每个子流域分别（循环）计算，每个sheet表对应一个子流域的降水，最后需要加在一起才是流域出口断面流量。if__name__=='__main__': ... 忽略数据读取代码和保存至excel代码 ...forkinrange(sheet_count):print(f"第{k+1}次运行ing")df=read_data(k)#这个是自编的函数foriinrange(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单位为hcalculate_Q(df,i,U)calculate_Qi(df,i,N[k])... 忽略保存至excel代码 ...1234567891011121314151617181920这里没有给出完整的读取数据、处理数据、绘制图表代码，需要根据自己的数据来写，不断思考和寻找问题答案的过程才是编程！。以下是用于计算的Excel格式4总结新安江模型是极为重要的水文模型，此拙作只作为简单的引导，希望大家在阅读的时候能够有自己的思考。如有问题，欢迎探讨，邮箱：mymail_zo@qq.com20240819更新：实在不会的也可以从该链接自费获取源代码