Python水文三水源新安江模型代码编写+降雨流量过程图绘制

广东自由人

Tips：        本期向大家分享三水源新安江模型代码的编写和最终的降雨流量过程图绘制，需要注意的是：①为便于读者理解过程，将“蒸散发、产流、分水源、汇流”4个部分的代码按顺序编写；②各阶段的参数在每个阶段开始时进行定义赋值，方便读者理解参数含义和用途；③因为并没有划分子流域，所以也没有使用定义函数来调用，显得更直观清晰一些。1新安江模型简介             新安江模型按照三层蒸散发模式计算流域蒸散发，按蓄满产流概念计算降雨产生的总径流量，采用流域蓄水曲线考虑下垫面不均匀对产流面积变化的影响。在径流成分划分方面，对三水源情况，按“山坡水文学”产流理论用一个具有有限容积和测孔、底孔的自由水蓄水库把总径流划分成饱和地面径流、壤中水径流和地下水径流。在汇流计算方面，单元面积的地面径流汇流一般采用单位线法，壤中水径流和地下水径流的汇流则采用线性水库法。河网汇流一般采用分段连续演算的Muskingum法或滞时演算法。新安江模型原理        具体的新安江模型原理的学习主要参考了“贵仁模型技术手册”、“Ai尚研修相关文章”、水文预报教材和相关论文”。（给出相关链接如下）2.1新安江模型·贵仁模型技术手册(keepsoft.net)【水文模型】三水源新安江模型(qq.com)水文预报（第5版）_百度百科(baidu.com)2三水源新安江模型代码2.1导入相关库importos#导入os模块，用于文件操作fromopenpyxlimportload_workbook#导入openpyxl模块，用于读取Excel文件importnumpyasnp#导入numpy模块importmatplotlib.pyplotasplt#导入绘图模块frommatplotlibimportrcParams#设置matplotlib使用的字体fromscipy.interpolateimportmake_interp_spline#折线平滑处理的方法2.2读取Excel数据2.2.1数据展示2.2.2读取数据filename='test01_zhengsanfa_chanliu.xlsx'ifos.path.exists(filename):wb=load_workbook(filename)ws=wb.activerows=ws.rowsheader=[]fors_headerinnext(rows):header.append(s_header.value)date=[]forrowinrows:temprow_dict={}forheadercell,datecellinzip(header,row):temprow_dict[headercell]=datecell.valuedate.append(temprow_dict)defduqu(key):lst=[]foriinrange(len(date)):dic_1=date[i]#print(dic_1)lst.append(dic_1[key])returnlstp=duqu('P')em=duqu('EM')t=duqu('t')2.3蒸散发＋产流阶段##---------------蒸散发（三层蒸散发模型）+产流--------------------#以下为初始参数设定：WU=3#昨天的WL=15#昨天的WD=20#昨天的K=0.85WUM=20WLM=80WDM=20C=0.15#土壤深层蒸散发系数WM=120#流域平均蓄水容量IM=0.01B=0.4WMM=(1+B)*WM/(1-IM)W=WU+WL+WD#就是今天的初始土壤含水量（昨天结束的土壤含水量）lst1=[]#存储Rlst2=[]#EUlst3=[]#ELlst4=[]#EDlst5=[]#Elst6=[]#WUlst7=[]#WLlst8=[]#WDlst9=[]#Wforiinrange(len(p))=p[i]EM=em[i]a=WMM*(1-(1-W/WM)**(1/(B+1)))ifP-K*EM=K*EM:EU=K*EMEL=0ED=0ifP=C:EL=(K*EM-EU)*WL/WLMED=0WL-=ELelse:ifWL>=C*(K*EM-EU):EL=C*(K*EM-EU)ED=0WL-=ELelse:EL=WLED=C*(K*EM-EU)-WLWL=0ifWD>ED:WD-=EDelse:WD=0E=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进行对比，经验证结果一致，说明代码正确可行。《流域水文模拟——新安江模型与陕北模型》代码运算结果《流域水文模拟——新安江模型与陕北模型》P77表格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]#存储FRlst_Si=[]#Si=S0*FRO/FRlst_AU=[]#AUlst_RS=[]#RSlst_S=[]#Slst_RI=[]#RIlst_RG=[]#RGlst_S1=[S0]#存放下一阶段S0，即S1foriinrange(len(p))=p[i]EM=em[i]R=lst1[i]#产流：由上一步得到S0=lst_S1[i]FR0=lst_FR[i]ifP-K*EM