簡介這份資源是一份面向水文水資源專業(yè)學(xué)生與從業(yè)者的三層蒸發(fā)蓄滿產(chǎn)流模型新安江模型Python計算程序用于完成基于降水與蒸發(fā)能力數(shù)據(jù)的產(chǎn)流模擬計算。程序核心實現(xiàn)新安江模型中的三層蒸發(fā)與蓄滿產(chǎn)流機(jī)制適合學(xué)習(xí)流域水文模型、開展課程設(shè)計或進(jìn)行科研驗證的初級至中級Python使用者。資源包為1個py文件體積約1KB內(nèi)容緊湊主程序需配合作者提供的PEdata.xlsx降水蒸發(fā)數(shù)據(jù)運行可直接復(fù)制使用或在此基礎(chǔ)上修改參數(shù)。截至目前已有2642人瀏覽學(xué)習(xí)說明該程序在同類資源中具有較高實用參考價值。通過運行這份代碼讀者可以直觀理解新安江模型的計算流程快速上手實現(xiàn)產(chǎn)流模擬為后續(xù)開展流域水文分析提供可復(fù)用的腳本基礎(chǔ)。1. 三層蒸發(fā)蓄滿產(chǎn)流模型新安江模型python計算程序.py從降雨到徑流的確定性推演做水文預(yù)報的人遲早會碰到這個場景手里有流域的逐日降雨和蒸發(fā)數(shù)據(jù)想算出對應(yīng)的徑流過程但商業(yè)化軟件黑匣子太多參數(shù)一改就不知道內(nèi)部發(fā)生了什么。三層蒸發(fā)蓄滿產(chǎn)流模型常說的新安江模型就是干這個事的——它把流域降水轉(zhuǎn)化成蒸散發(fā)、土壤蓄水、地表徑流和壤中流是一個完整閉合的水量平衡計算框架。標(biāo)題里這個python計算程序.py本質(zhì)是一份可以直接運行的腳本解決的是沒有現(xiàn)成水文模型庫時如何用幾百行代碼實現(xiàn)經(jīng)典蓄滿產(chǎn)流計算的現(xiàn)實問題。適合水文專業(yè)學(xué)生做課程設(shè)計、科研人員做對比實驗以及工程師做中小流域的洪水預(yù)報方案初篩。只要你有逐日降雨、蒸發(fā)和流量數(shù)據(jù)這套代碼就能從零開始推演出模擬徑流過程還能人工調(diào)參逼近實測值。核心并不神秘蓄滿產(chǎn)流的核心就是先算土壤缺水量再決定產(chǎn)不產(chǎn)流。2. 蓄滿產(chǎn)流與三層蒸發(fā)模型運行的底層邏輯2.1 張力水和自由水蓄滿產(chǎn)流的核心區(qū)分蓄滿產(chǎn)流模型能成立是因為它區(qū)分了兩種性質(zhì)完全不同的土壤水。張力水是被土粒吸附、能被植物根系吸收的水分它不參與重力排水自由水則是填充在土壤孔隙中、受重力支配可以橫向或縱向流動的水分。降雨先補(bǔ)充張力水張力水蓄滿后才開始產(chǎn)生自由水自由水超過一定容量就形成徑流。這個先蓄滿再產(chǎn)流的機(jī)制是模型最關(guān)鍵的判斷條件。實際計算中每個計算單元都有一個張力水蓄水容量WM它等于上層張力水容量WUM、下層張力水容量WLM和深層張力水容量WDM之和。降雨先補(bǔ)上層上層滿了溢到下層下層滿了溢到深層。只有當(dāng)整個土柱的張力水都蓄滿降雨才轉(zhuǎn)為自由水蓄量進(jìn)而劃分產(chǎn)流。這個順序保證了水量平衡的閉合性也讓參數(shù)物理意義明確——WUM和WLM可以根據(jù)土壤類型和經(jīng)驗值初定再通過率定微調(diào)。2.2 三層蒸發(fā)為什么不能只算一層很多初學(xué)模型的人會問一個流域的總蒸發(fā)量乘以一個系數(shù)不就完了為什么要分三層答案在于干旱期和濕潤期的蒸發(fā)規(guī)律截然不同。濕潤期土壤含水充足蒸發(fā)接近蒸散發(fā)能力干旱期表層先干蒸發(fā)衰減很快但深層仍有水分緩慢供給。三層蒸發(fā)結(jié)構(gòu)正是為了模擬這種表層快、深層慢的衰減過程。模型里每層有自己的蒸發(fā)能力分配系數(shù)。上層蒸發(fā)EU取蒸散發(fā)能力EM乘以系數(shù)K如果上層有足夠水量就直接消耗不夠時下層蒸發(fā)EL取剩余蒸散發(fā)能力乘以一個折減系數(shù)但最多不能超過下層張力水蓄量的某個比例再不夠才動用深層ED深層蒸發(fā)速率接近常數(shù)。這一套規(guī)則寫進(jìn)代碼就三個條件分支卻決定了連續(xù)干旱期的模擬精度。2.3 從降雨到徑流水源劃分和匯流張力水全部蓄滿后超出的降雨進(jìn)入自由水蓄量SM。自由水蓄量有一個最大的容量超過了就產(chǎn)地面徑流RS同時自由水還按一個固定的消退系數(shù)KG向下補(bǔ)給深層地下水形成地下徑流RG。壤中流RI則和自由水蓄量成正比計算為KI乘以自由水蓄量其中KI和KG是消退系數(shù)。降雨產(chǎn)流過程結(jié)束后各水源按不同路徑進(jìn)入河道地面徑流幾乎即時匯入壤中流衰減稍慢地下徑流衰減最慢。模型通常用一個線性水庫或單位線把三種水源分別演算到流域出口。這個先分水源、再分路徑的做法讓模型既能模擬洪水過程的陡漲陡落又能模擬退水過程的緩降比單一水庫模型靈活得多。3. 編寫新安江模型Python計算程序蓄滿產(chǎn)流演算模塊的完整實現(xiàn)3.1 數(shù)據(jù)輸入和初始化直接讀取Excel更方便我一般用pandas讀取Excel里的三列數(shù)據(jù)日期、降雨量、實測流量。蒸散發(fā)能力可以單獨一列也可以由程序根據(jù)月份自動調(diào)整。初始化時先定義土壤含水量初值、蓄水容量參數(shù)、蒸發(fā)系數(shù)以及各水源的消退系數(shù)。代碼定義得越清晰后面調(diào)參越容易。import pandas as pd import numpy as np # 讀取數(shù)據(jù) df pd.read_excel(流域日數(shù)據(jù).xlsx, sheet_nameSheet1) dates df[日期].values P df[降雨].values.astype(float) # 日降雨單位mm EM_input df[蒸發(fā)].values.astype(float) # 日蒸散發(fā)能力單位mm Q_obs df[流量].values.astype(float) # 實測流量單位m3/s用于對比 # 模型參數(shù)初始值 WUM, WLM, WDM 20, 90, 40 # 上、下、深層張力水容量(mm) WM WUM WLM WDM # 總張力水容量(mm) K 1.0 # 蒸散發(fā)能力折算系數(shù) SM 15.0 # 自由水蓄水容量(mm) KG, KI 0.3, 0.35 # 地下水和壤中流出流系數(shù) IM 0.1 # 不透水面積比例 FC 10.0 # 穩(wěn)定下滲率(mm/day) # 注原文此處為幫助記憶的占位符 # 初值 WU, WL, WD 10.0, 50.0, 10.0 # 三層初始張力水蓄量(mm) S_free 5.0 # 初始自由水蓄量(mm) Q_sim []這段代碼解決的問題是模型跑起來之前所有狀態(tài)變量從哪來。參數(shù)初值不要求精確后續(xù)會率定。比較關(guān)鍵的是前三層初始蓄水量的設(shè)置——如果給得太小模擬初期會有一段填洼過程模擬流量偏低如果給得太大初期徑流會虛高。3.2 三層蒸發(fā)計算模塊核心的日水量平衡每天的演算遵循固定順序先算蒸發(fā)再算產(chǎn)流最后更新各層蓄量。蒸發(fā)計算的先后不能顛倒因為蒸發(fā)消耗的是當(dāng)天的土壤蓄水而產(chǎn)流用的是當(dāng)天的降雨增量。# 三層蒸發(fā)計算邏輯 def calc_evap(P_day, EM_day, WU, WL, WD): # 1. 上層蒸發(fā)優(yōu)先消耗上層水量 EU EM_day * K if WU EU: WU WU - EU else: EU WU WU 0 # 2. 下層蒸發(fā)剩余蒸散發(fā)能力按折減系數(shù)分配 EL (EM_day * K - EU) * 0.6 # 0.6是下層蒸發(fā)折減系數(shù) if WL EL: WL WL - EL else: EL WL WL 0 # 3. 深層蒸發(fā)只占剩余蒸發(fā)的小部分 ED (EM_day * K - EU - EL) * 0.1 if WD ED: WD WD - ED else: ED WD WD 0 return EU, EL, ED, WU, WL, WD函數(shù)內(nèi)部的三個if判斷體現(xiàn)的是能供多少就蒸發(fā)多少的規(guī)則。比如上層含水量只有5mm但當(dāng)天蒸發(fā)能力是8mm那就只能蒸發(fā)5mm余下的蒸發(fā)需求依次往下層、深層傳遞。這也解釋了為什么連續(xù)干旱天模擬蒸發(fā)衰減快——上層很快耗盡下層折減系數(shù)0.6進(jìn)一步限制了蒸發(fā)速率。3.3 蓄滿產(chǎn)流和自由水蓄量降雨轉(zhuǎn)化過程蒸發(fā)算完后再算產(chǎn)流。核心操作是判斷加雨后的總張力水蓄量是否超過WM。如果超過超出部分進(jìn)入自由水蓄量如果沒超過降雨全部用于補(bǔ)充土壤缺水量不產(chǎn)流。這部分的順序、條件、更新邏輯是模型的重中之重。# 產(chǎn)流計算模塊 def calc_runoff(P_day, WU, WL, WD, S_free): # 凈雨扣除不透水面積直接產(chǎn)流 Pe P_day * (1 - IM) P_imp P_day * IM # 不透水面積上的直接徑流 # 1. 降雨先補(bǔ)充張力水 remaining P_day * (1 - IM) # 上層補(bǔ)充 if WU WUM: dWU min(WUM - WU, remaining) WU dWU remaining - dWU # 下層補(bǔ)充 if WL WLM and remaining 0: dWL min(WLM - WL, remaining) WL dWL remaining - dWL # 深層補(bǔ)充 if WD WDM and remaining 0: dWD min(WDM - WD, remaining) WD dWD remaining - dWD # 2. 張力水蓄滿后剩余水量進(jìn)入自由水蓄量 dS remaining # 剩余水量全部進(jìn)入自由水蓄量 # 3. 自由水蓄量超過容量時產(chǎn)生地面徑流 if S_free dS SM: RS (S_free dS - SM) S_free SM else: RS 0 S_free S_free dS # 4. 自由水的壤中流和地下徑流分配 RI KI * S_free RG KG * S_free # 更新自由水蓄量扣除流出部分 S_free S_free - RI - RG # 總徑流 地面徑流 壤中流 地下徑流 不透水面積直接徑流 R_total RS RI RG P_imp return R_total, WU, WL, WD, S_free, RS, RI, RG產(chǎn)流劃分的關(guān)鍵在于自由水蓄量的庫容有限。降雨水超過SM的部分直接成為地面徑流而壤中流和地下徑流的量由KI、KG和當(dāng)前S_free共同決定。如果把SM調(diào)大地面徑流會減少壤中流增多洪水過程會變得坦化如果把SM調(diào)小地面徑流占比大洪峰更尖更陡。這部分參數(shù)直接控制模擬的洪峰形態(tài)。3.4 匯流計算三種水源分別演算產(chǎn)流計算結(jié)束后地面徑流RS直接進(jìn)入出口斷面壤中流RI經(jīng)過一個線性水庫衰減地下徑流RG經(jīng)過另一個衰減更慢的線性水庫。線性水庫的退水公式是今天出流量等于昨天的剩余蓄量乘以消退系數(shù)。# 匯流模塊線性水庫演算 def route_flow(RS_day, RI_day, RG_day, S_RI_prev, S_RG_prev, CS, CG): # 地面徑流直接匯入但做單位線或線性水庫處理 # 這里使用一個簡單線性水庫模擬地面匯流 Q_surface CS * S_RI_prev (1 - CS) * RS_day S_surface S_RI_prev * (1 - CS) (1 - CS) * RS_day # 簡化示意 # 壤中流線性水庫 Q_inter CS * S_RI_prev (1 - CS) * RI_day S_RI_prev_new S_RI_prev * (1 - CS) RI_day * (1 - CS) # 地下徑流線性水庫 Q_ground CG * S_RG_prev (1 - CG) * RG_day S_RG_prev_new S_RG_prev * (1 - CG) RG_day * (1 - CG) # 總出流量轉(zhuǎn)換為m3/s Q_total (Q_surface Q_inter Q_ground) * 流域面積_km2 * 1000 / 86400 return Q_total, S_RI_prev_new, S_RG_prev_new匯流部分在實際工程中經(jīng)常被簡化。如果只有一個日尺度的流量序列用線性水庫就夠了不需要復(fù)雜的單位線。參數(shù)CS反映地面匯流速度取值0.3~0.8越大表示匯流越快CG一般取0.9以上表示地下水消退非常緩慢。這里單位換算要注意每天的徑流深(mm)乘以流域面積再除以86400秒才能轉(zhuǎn)換為m3/s。4. 參數(shù)率定和敏感性分析讓模擬結(jié)果靠近實測值4.1 人工試錯法和納什效率系數(shù)模型參數(shù)有七八個直接上來就全自動尋優(yōu)容易陷入局部最優(yōu)而且你不知道哪個參數(shù)起了主要作用。常見做法是先手動調(diào)一遍讓模擬過程線大致貼合實測再用自動優(yōu)化微調(diào)。評價模擬好壞的指標(biāo)最常用納什效率系數(shù)NSE它衡量模擬和實測的吻合程度計算公式是1減去誤差方差除以實測方差。def calc_nse(Q_obs, Q_sim): obs_mean np.mean(Q_obs) numerator np.sum((Q_obs - Q_sim) ** 2) denominator np.sum((Q_obs - obs_mean) ** 2) return 1 - numerator / denominatorNSE等于1表示完美擬合0表示模擬效果等同用實測均值預(yù)測負(fù)數(shù)表示模型比均值預(yù)測還差。實際工程中NSE達(dá)到0.7以上就算可用0.8以上效果良好。調(diào)參順序也有講究先調(diào)WM和K讓年水量平衡閉合再調(diào)SM控制洪峰形態(tài)最后調(diào)KI、KG控制退水過程。4.2 參數(shù)敏感性哪個參數(shù)最影響結(jié)果我做過一組實驗把各個參數(shù)單獨拉偏20%觀察NSE變化幅度。結(jié)果最敏感的是蒸發(fā)折算系數(shù)K它對年徑流總量影響極大——K從1.0調(diào)到1.2年徑流可能會減少20%以上因為多算的蒸散發(fā)直接減少了產(chǎn)流。其次是WM它控制著流域的總體蓄水能力WM太小模擬洪峰偏高WM太大則洪峰偏低。SM和KI、KG的敏感性相對中等它們主要影響過程的形態(tài)而不是總量。參數(shù)率定建議先跑一次基準(zhǔn)模擬輸出年徑流量和實測年徑流量對比。如果模擬偏大大概率是K偏小或者WM偏小優(yōu)先調(diào)這兩個。如果年總量對上了但洪峰偏高再調(diào)SM。這種分步調(diào)整比一次性盲目搜索效率高得多也更容易攢出調(diào)參經(jīng)驗。4.3 率定腳本批量搜索加人工判斷# 簡單的網(wǎng)格搜索率定 best_nse -999 best_params {} for test_WM in [100, 120, 150, 180]: for test_SM in [10, 15, 20, 25]: for test_K in [0.8, 1.0, 1.2]: # 運行模型得到Q_sim nse calc_nse(Q_obs, Q_sim) if nse best_nse: best_nse nse best_params {WM: test_WM, SM: test_SM, K: test_K}網(wǎng)格搜索看起來粗暴但勝在直觀可靠。三層嵌套的循環(huán)把所有組合跑一遍每組參數(shù)都調(diào)用一次完整的模擬計算最后選出NSE最高的一組。注意每組參數(shù)跑模擬之前要讓模型有一段預(yù)熱期比如前30天的模擬結(jié)果不參與NSE計算用來消除初始條件的影響。5. 模型調(diào)試避坑運行報錯和結(jié)果失真排查5.1 數(shù)值溢出Python報numpy.inf或nan現(xiàn)象運行一段時間后模擬流量變成inf或nan過程線斷掉。原因通常是某個狀態(tài)變量被算成了負(fù)值或無窮大。比如蒸發(fā)計算時上層含水量已經(jīng)是0但代碼里仍然嘗試減去蒸發(fā)量導(dǎo)致負(fù)值產(chǎn)生。累計下去自由水蓄量或張力水蓄量變成無意義的大數(shù)最終溢出。解決在每個計算函數(shù)入口加一個保護(hù)判斷任何蓄量變量如果出現(xiàn)負(fù)值就強(qiáng)制歸零另一個常見原因是蒸發(fā)輸入數(shù)據(jù)有缺測或零值導(dǎo)致EM_day * K為0但后續(xù)計算中分母出現(xiàn)了0。檢查原始Excel數(shù)據(jù)里的空值用pandas的dropna或者在讀取后直接用fillna(0)填充。WU max(WU, 0) WL max(WL, 0) WD max(WD, 0) S_free max(S_free, 0)5.2 模擬徑流一直為0檢查蓄滿條件是否永遠(yuǎn)不滿足現(xiàn)象模型跑完模擬流量全是0一條直線。原因大概率是WM過大或降雨太小所有降雨都被土壤吸收從未蓄滿。另一種情況是蒸散發(fā)輸入值異常大導(dǎo)致土壤始終處于虧水狀態(tài)。解決第一步在產(chǎn)流計算模塊里加一行print把每天的WUWLWD和WM對比輸出觀察蓄水量是否在逼近WM。如果連續(xù)幾十天蓄量都在WM以下說明降雨和蒸發(fā)不匹配。第二步把SM調(diào)小讓它更容易產(chǎn)流。第三步檢查蒸發(fā)數(shù)據(jù)是否使用了總的潛在蒸散發(fā)而不是實際蒸散發(fā)能力。如果確實用了潛在的考慮把K調(diào)小到0.5~0.7來折減。5.3 洪水過程模擬不錯退水過程偏高或偏低現(xiàn)象洪峰對上了但退水尾巴太慢或者太快。原因在于KI和KG的比例不對。KG太大地下水補(bǔ)充太快退水過程降不下來KG太小退水過程又太平緩。另一個因素時線性水庫的消退系數(shù)CS和CG設(shè)置不合理CS取值太小造成地面匯流延遲過大導(dǎo)致洪峰后還有一個不小的拖尾流量。解決單獨拿出退水段的數(shù)據(jù)畫出半對數(shù)圖。如果實測退水段在半對數(shù)坐標(biāo)下是一條直線說明退水主要受單一線性水庫控制這時候就只需要調(diào)CG。如果實測退水是折線前段陡后段緩說明地面-壤中流-地下水三分水源的消退系數(shù)都需要調(diào)整。具體到代碼先固定CS0.5、CG0.95再單獨微調(diào)KG從0.2到0.4看NSE變化。5.4 年徑流總量偏差大引入水量平衡誤差檢查現(xiàn)象模擬年徑流總量比實測偏大30%以上。原因可能是不透水面積比例IM設(shè)置過高或者蒸散發(fā)能力數(shù)據(jù)偏小。IM直接決定有多少降雨不參與蓄滿產(chǎn)流而直接成為徑流IM每增加0.1年徑流就會顯著增加。解決在代碼末尾加一個年度水量平衡的統(tǒng)計輸出打印降雨總量、蒸散發(fā)總量、徑流總量三項。如果降雨減去蒸發(fā)減去徑流后剩余的水量不為零說明存在計算bug或蓄量變化不閉合。逐項檢查后優(yōu)先把IM調(diào)回0.05以下同時把K往1.1~1.2方向調(diào)整讓蒸散發(fā)消耗掉更多水分。5.5 初值影響前幾天的模擬結(jié)果預(yù)熱期不夠現(xiàn)象模擬前10天結(jié)果和實測差距很大后面逐漸正常。原因是初始土壤含水量設(shè)置不合理。如果初始WU、WL、WD設(shè)得太小模型需要一段時間填蓄這段時間模擬徑流幾乎為0。解決常見的處理方式有兩種。第一種是把模擬起點提前一年最后使用后半段的結(jié)果。第二種是在代碼中增加一個預(yù)熱循環(huán)使用實測降雨的前90天數(shù)據(jù)連續(xù)計算讓狀態(tài)變量收斂到合理范圍然后再開始記錄模擬結(jié)果。第二種更簡單不用額外拼接數(shù)據(jù)只需在計算循環(huán)開始前多跑幾十次同樣的模擬邏輯。6. 參數(shù)邊界和率定順序的實操心得用好模型的關(guān)鍵技巧調(diào)參順序比參數(shù)值本身更重要。這條建議反復(fù)說因為很多同學(xué)上來就把全部參數(shù)扔給自動尋優(yōu)算法結(jié)果幾千次迭代后NSE還沒到0.6回頭問為什么。原因是參數(shù)之間存在明顯的相關(guān)性比如WM和K都會影響年徑流總量如果不固定其中一個直接一起搜算法會陷入斜向的山脊里怎么跳都找不到最優(yōu)。我現(xiàn)在的習(xí)慣是先固定各層張力水容量比例為2:9:4K取1.0SM取15KG取0.3也就是用一套普通參數(shù)把模型跑通看過程線的基本形態(tài)再進(jìn)入調(diào)參階段。第一步只調(diào)K緊盯年徑流總量讓水量平衡誤差在5%以內(nèi)。第二步調(diào)WM確保濕潤期產(chǎn)流不偏早不偏晚。第三步調(diào)SM修正洪峰的坦化和銳化。第四步調(diào)KG和KI修正退水過程。第五步才允許微調(diào)各層蒸發(fā)系數(shù)。中間還有幾個細(xì)節(jié)容易被忽略自由水蓄量的初始值S_free如果設(shè)置過大前幾天的壤中流會異常偏大這個問題在北方干旱流域尤其嚴(yán)重另一個常見問題是蒸發(fā)數(shù)據(jù)的量綱日蒸散發(fā)能力一般小于8mm如果某人給的數(shù)據(jù)單位是月蒸發(fā)量而你直接當(dāng)成日值用了模型會一直蓄不滿水還有流域面積的單位換算調(diào)用匯流模塊前確認(rèn)面積單位是km2而不是m2或者畝這個錯了輸出流量會差好幾個數(shù)量級。我的做法是在代碼開頭統(tǒng)一加一個assert斷言檢查關(guān)鍵參數(shù)的范圍比如0KG1、WM50、SM0超出范圍直接報錯避免帶著錯誤參數(shù)跑幾百行代碼最后才發(fā)現(xiàn)問題。跑完模型之后還有一件值得做的事把模擬和實測的雙累積曲線畫出來對比。雙累積曲線的斜率如果在中后段明顯改變大概率是參數(shù)K需要分段取值或者流域內(nèi)有水利工程調(diào)蓄導(dǎo)致下墊面條件變化。我自己做過一個案例某河流域前半段模擬效果很好后半段模擬流量整體偏低查了資料發(fā)現(xiàn)后段有多個小型水庫開始蓄水蒸發(fā)和調(diào)度都變了再調(diào)K也救不回來只能分段率定。水文模型不是萬能鑰匙參數(shù)只能反映下墊面的平均狀態(tài)子流域的調(diào)蓄改變屬于結(jié)構(gòu)性誤差這時候往參數(shù)上找補(bǔ)只會讓模型過擬合。經(jīng)驗上一篇完整的新安江模型Python程序應(yīng)該包含參數(shù)配置文件建議用yaml或json分開寫、核心計算模塊蒸發(fā)、產(chǎn)流、匯流三個函數(shù)、率定模塊網(wǎng)格搜索或遺傳算法、結(jié)果輸出模塊生成過程線圖和NSE、水量平衡誤差報告。這些代碼加起來大概500行不算長但結(jié)構(gòu)清晰后續(xù)無論是換流域還是改時間尺度都能快速適配。希望這篇筆記能把你的新安江模型程序從能跑帶到跑得明白、調(diào)得有方向祝你好運。本文還有配套的精品資源點擊獲取