matlab液体湿润模拟,【水文模型】01 三水源新安江模型
目錄
一、概述
1.1 流域概況
1.2 產流方式論證
1.3 設計基本任務
1.4 設計原始資料
二、新安江模型
2.1 流域劃分
2.2 產匯流計算
2.2.1 蒸散發計算
2.2.2 產流計算
2.2.3 水源劃分
2.2.4 匯流計算
2.3 模型參數
2.3.1 參數物理意義
2.3.2 參數率定結果
2.4 模擬結果
2.4.1 精度統計表
2.4.2 計算與實測流量過程線比較
三、精度統計與誤差分析
3.1 精度統計
3.2 誤差分析
3.2.1 模型參數的影響
3.2.2 人類活動的影響
3.2.3 資料代表性的影響
四、參考文獻
五、程序代碼
數據文件
一、概述
本博客基于河海大學水文水資源學院水文預報課程設計編寫。采用MATLAB語言編寫三水源新安江模型，并用人機交互率定模型參數。模型運行所需文件請從新安江模型數據下載。
1.1 流域概況
呈村流域控制面積為290 km2，地勢南高北低，相對高差較大，平均海拔高程為583m，流域河道平均坡度為0.95%，最大匯流路徑長度為36 km。流域內植被良好，雨量充沛，多年平均降雨量約為2100 mm，流域降水在年內年際分配極不均勻，為典型的濕潤流域。該流域植被類型主要包括常綠針葉林、落葉闊葉林、混合林、森林地、林地草原、牧草地與作物地，土壤類型主要為黏壤土。
1.2 產流方式論證
根據水文氣象條件、下墊面條件及流量過程線分析：
多年平均降雨量約為2100 mm>1000 mm；
流域降水在年內年際分配極不均勻，為典型的濕潤流域；
土壤類型主要為黏壤土，土質疏松，不易超滲；
綜上，該流域符合蓄滿產流模式。
1.3 設計基本任務
編制預報方案，根據流域概況和原始資料編制預報方案；
編寫新安江模型，時間尺度：日模型(24h)，水源劃分：兩水源(可選)和三水源(可選)；
根據已給的呈村流域資料，利用編制的新安江模型進行日徑流模擬，率定新安江模型參數；
分析日模型模擬結果，精度評定時，日模型采用徑流深相對誤差與確定性系數。
1.4 設計原始資料
逐日資料，包括呈村站流量資料、蒸發資料，呈村、汪村、樟源口、棣甸、董坑塢、用功城、左龍、馮村、田里與大連站雨量資料。資料起止時間見表1-1。
參數建議取值見表1-2。
各雨量站控制面積及比例，見表1-3。
二、新安江模型
2.1 流域劃分
根據流域地形、地貌條件及布設的雨量站網，用泰森多邊形法將呈村流域劃分為10塊單元面積，各單元面積的權重見表1-3。
2.2 產匯流計算
對每塊單元面積采用三水源新安江模型分別進行蒸散發計算、產流計算、水源劃分和匯流計算，得到單元面積出流過程；由于本流域面積較小，單元面積的出流過程不采用馬斯京根分段連續演算法進行出口以下的河道洪水演算，而是通過調整滯時 L 求得單元面積在流域出口的流量過程線；將每個單元面積在流域出口的流量過程線線性疊加，即為呈村流域的出口流量過程。
2.2.1 蒸散發計算
(1)當W U + P ≥ E P WU+P≥E_PWU+P≥EP?時
E U = E P ， E L = 0 ， E D = 0 E_U=E_P，E_L=0，E_D=0EU?=EP?，EL?=0，ED?=0
(2)當W U + P < E P ， W L ≥ C ? W L M WU+PWU+P
E U = W U + P ， E L = ( E P ? E U ) W L / W L M ， E D = 0 E_U=WU+P，E_L=(E_P-E_U)WL/WLM，E_D=0EU?=WU+P，EL?=(EP??EU?)WL/WLM，ED?=0
(3)當W U + P < E P ， C ( E P ? E U ) ≤ W L < C ? W L M WU+PWU+P
E U = W U + P ， E L = C ( E P ? E U ) ， E D = 0 E_U=WU+P，E_L=C(E_P-E_U)，E_D=0EU?=WU+P，EL?=C(EP??EU?)，ED?=0
(4)當W U + P < E P ， W L < C ( E P ? E U ) WU+PWU+P
E U = W U + P ， E L = W L ， E D = C ( E P ? E U ) ? E L E_U=WU+P，E_L=WL，E_D=C(E_P-E_U )-E_LEU?=WU+P，EL?=WL，ED?=C(EP??EU?)?EL?
式中：
W U WUWU——上層土壤含水量，mm；
W L WLWL——下層土壤含水量，mm；
P PP——降雨量，mm；
E P E_PEP?——流域蒸發能力，mm；
E U E_UEU?——上層蒸發量，mm；
E L E_LEL?——下層蒸發量，mm；
E D E_DED?——深層蒸發量，mm；
C CC——蒸發擴散系數；
W L M WLMWLM——下層土壤含水容量，mm。
2.2.2 產流計算
流域平均土壤水蓄積容量關系
W M = W M M 1 + b WM=\frac{WMM}{1+b}WM=1+bWMM?
a = W M M [ 1 ? ( 1 ? W W M ) 1 1 + b ] a=WMM[1-(1-\frac{W}{WM})^{\frac{1}{1+b}} ]a=WMM[1?(1?WMW?)1+b1?]
(1)當P E ≤ 0 PE≤0PE≤0時
R = 0 R=0R=0
(2)當P E > 0 ， a + P E ≤ W M M PE>0，a+PE≤WMMPE>0，a+PE≤WMM時
R = P E + W ? W M + W M ( 1 ? P E + a W M M ) b + 1 R=PE+W-WM+WM(1-\frac{PE+a}{WMM})^{b+1}R=PE+W?WM+WM(1?WMMPE+a?)b+1
(3)當P E > 0 ， a + P E > W M M PE>0，a+PE>WMMPE>0，a+PE>WMM時
R = P E + W ? W M R=PE+W-WMR=PE+W?WM
式中
W M WMWM——流域平均蓄水容量，mm；
W M M WMMWMM——流域最大蓄水容量，mm；
a aa——流域最大初始土壤含水量，mm；
b bb——常數，反映流域包氣帶蓄水容量分布的不均勻性
W WW——流域平均初始土壤含水量，mm；
P E PEPE——扣除雨期蒸發后的降雨量，mm；
R RR——產流量，mm。
2.2.3 水源劃分
(1)二水源劃分
α = R P E α=\frac{R}{PE}α=PER?
當PE>FC時
R G = α ? F C RG=α?FCRG=α?FC
R S = α ? ( P E ? F C ) RS=α?(PE-FC)RS=α?(PE?FC)
當P E ≤ F C PE≤FCPE≤FC時
R G = R RG=RRG=R
R S = 0 RS=0RS=0
式中
α \alphaα——產流面積比例；
P E PEPE——扣除雨期蒸發后的降雨量，mm；
R RR——產流量，mm；
F C FCFC——流域的穩定下滲量，mm；
R G RGRG——地下徑流，mm；
R S RSRS——地面徑流，mm。
(2)三水源劃分
流域平均自由水蓄積容量關系
S m = S m m 1 + E X S_m=\frac{S_{mm}}{1+EX}Sm?=1+EXSmm??
A U = S m m [ 1 ? ( 1 ? S 1 F R 1 F R S m ) 1 1 + E X ] AU=S_{mm} [1-(1-\frac{S_1\frac{FR_1}{FR}}{S_m})^{\frac{1}{1+EX}}]AU=Smm?[1?(1?Sm?S1?FRFR1???)1+EX1?]
當P E + A U < S m m PE+AUPE+AU
R S = F R [ P E + S 1 ? F R 1 F R ? S m + S m ( 1 ? P E + A U S m m ) E X + 1 ] RS=FR[PE+\frac{S_1?FR_1}{FR}-S_m+S_m (1-\frac{PE+AU}{S_{mm}} )^{EX+1} ]RS=FR[PE+FRS1??FR1???Sm?+Sm?(1?Smm?PE+AU?)EX+1]
當P E + A U ≥ S m m PE+AU≥S_{mm}PE+AU≥Smm?時
R S = F R ( P E + S 1 ? F R 1 F R ? S m ) RS=FR(PE+\frac{S_1?FR_1}{FR}-S_m )RS=FR(PE+FRS1??FR1???Sm?)
本時段的自由水蓄量為
S = S 1 ? F R 1 F R + R ? R S F R S=\frac{S_1?FR_1}{FR}+\frac{R-RS}{FR}S=FRS1??FR1??+FRR?RS?
相應的壤中流和地下徑流為
R I = K I ? S ? F R RI=KI?S?FRRI=KI?S?FR
R G = K G ? S ? F R RG=KG?S?FRRG=KG?S?FR
本時段末的自由水蓄量變為
S 1 = S ( 1 ? K I ? K G ) S_1=S(1-KI-KG)S1?=S(1?KI?KG)
式中
F R 1 、 F R FR_1 、FRFR1?、FR——上一時段和本時段的產流面積比例；
S m S_mSm?——流域平均自由水蓄量，mm；
S m m S_{mm}Smm?——流域最大自由水蓄量，mm；
A U AUAU——流域最大初始自由水蓄量，mm；
E X EXEX——自由水蓄量分布曲線指數
P E PEPE——扣除雨期蒸發后的降雨量，mm；
R G RGRG——地下徑流，mm；
R S RSRS——地面徑流，mm；
R I RIRI——壤中流，mm；
R RR——產流量，mm；
S 1 S_1S1?——時段初始自由水蓄量，mm；
S SS——本時段的自由水蓄量，mm。
2.2.4 匯流計算
(1)地面徑流的坡地匯流
地面徑流的坡地匯流時間不計，直接進入河網，計算公式為：
Q S ( I ) = R S ( I ) × U QS(I)=RS(I)×UQS(I)=RS(I)×U
式中：
Q S QSQS——地面徑流，m3/s；
U UU——單位轉換系數；
R S RSRS——地面徑流量，mm。
(2)壤中流匯流
表層自由水以K I KIKI側向出流后成為表層壤中流，進入河網。但如土層較厚，表層自由水尚可滲入深層土，經過深層土的調蓄作用，才進入河網。深層自由水用線性水庫模擬，其消退系數為C I CICI，計算公式為：
Q I ( I ) = C I × Q I ( I ? 1 ) + ( 1 ? C I ) × R I ( I ) × U QI(I)=CI×QI(I-1)+(1-CI)×RI(I)×UQI(I)=CI×QI(I?1)+(1?CI)×RI(I)×U
式中：
Q I QIQI——壤中流，m3/s；
C I CICI——壤中流消退系數；
R I RIRI——壤中流徑流量，mm。
(3)地下徑流匯流
地下徑流匯流用線性水庫模擬,其消退系數為C G CGCG，出流進入河網。表層自由水以K G KGKG向下出流后，再向地下水庫匯流的時間不另計，包括在C G CGCG之內，計算公式為：
Q G ( I ) = C G × Q G ( I ? 1 ) + ( 1 ? C G ) × R G ( I ) × U QG(I)=CG×QG(I-1)+(1-CG)×RG(I)×UQG(I)=CG×QG(I?1)+(1?CG)×RG(I)×U
式中：
Q G QGQG——地下徑流，m3/s；
C G CGCG——地下水消退系數；
R G RGRG——地下徑流徑流量，mm。
(4)單元面積河網匯流
單元面積的河網匯流用滯后演算法，參數有滯后量L LL與消退系數C S CSCS，計算公式為：
Q T ( I ) = Q S ( I ) + Q I ( I ) + Q G ( I ) QT(I)=QS(I)+QI(I)+QG(I)QT(I)=QS(I)+QI(I)+QG(I)
Q ( I ) = C S × Q ( I ? 1 ) + ( 1 ? C S ) × Q T ( I ? L ) Q(I)=CS×Q(I-1)+(1-CS)×QT(I-L)Q(I)=CS×Q(I?1)+(1?CS)×QT(I?L)
式中：
Q T QTQT——總徑流，m3/s；
C S CSCS——河網水流消退系數；
Q QQ——單元面積河網總入流，m3/s；
L LL——河網匯流滯時，d。
2.3 模型參數
新安江模型是一個通過長期實踐和對水文規律認識基礎上建立起來的一個概念性水文模型。模型大多數參數都具有明確的物理意義，它們在一定程度上反映了流域的基本水文特征和降雨徑流形成的物理過程。因此，原則上可以按其物理意義通過實測、實驗、比擬等方法確定。本課程設計采用參數的概念分析方法，即先按實測值或參數的物理意義初定參數初值范圍；然后根據輸入，通過模型計算輸出；再將輸出過程與實測過程進行比較，作優化調試；根據特定的目標函數——年徑流相對誤差和確定性系數確定參數的最優值。
2.3.1 參數物理意義
(1)蒸散發能力折算系數K KK
K是影響產流量計算最為重要和敏感的參數，產流計算中K控制著水量平衡。K主要反映流域平均高程和蒸發站高程之間差別的影響和蒸發皿蒸散發與陸面蒸散發間差別的影響。蒸散發能力的地區分布大體上反映了氣候和自然地理條件的影響，具有較為明顯的區域性規律。
(2)自由水蓄水容量S M SMSM
SM反映表土蓄水能力，其值受降雨資料時段均化的影響明顯。當以日作為時段長時，在土層很薄的山區，其值為10mm或更小一些；而在土深林茂透水性很強的流域，其值可取50mm或更大一些；一般流域在10~20mm之間。SM對地面徑流和地下徑流的比重起著決定性作用；SM大，則地下徑流所占比重相對大，地面徑流所占比重相對小，洪峰流量相對小；反之，SM小，則地下徑流所占比重相對小，地面徑流所占比重相對大，洪峰流量相對大。
(3)自由水蓄水庫對地下水和壤中流的日出流系數K G 、 K I KG、KIKG、KI
KG的大小反映基巖和深層土壤的滲透性，KI的大小反映表層土的滲透性。KG+KI代表出流的快慢，KG/KI代表地下徑流與壤中流的比，對于一個特定流域它們都是常數，1-(KG+KI)為消退系數，它決定了直接徑流的退水快慢。
(4)地下水消退系數C G CGCG
CG可根據枯季地下徑流的退水規律來推求，C G = Q t + ? t Q t CG=\frac{Q_{t+?t}}{Q_t}CG=Qt?Qt+?t??。若以日作為計算時段長，則CG=0.950-0.998，大致相當于消退歷時為20~500d。
(5)壤中流消退系數C I CICI
若無壤中流則CI→0，若壤中流豐富則CI→0.9，相當于匯流時間為10d。
(6)河網水流消退系數C S CSCS
CS代表坦化作用，其值取決于河網的地貌條件，可通過河網地貌推求。因與時段長短有關，其值應視河道特征和洪水特性而定。
(7)河網匯流滯時L LL
L代表平移作用，其值取決于河網的地貌條件，可通過河網地貌推求。
2.3.2 參數率定結果
原始資料中給定的參數值如表2-1所示。
利用人機交互率定方法得到的參數率定結果如下表2-2所示。
2.4 模擬結果
2.4.1 精度統計表
呈村流域日模型模擬結果及精度統計見表2-3。
2.4.2 計算與實測流量過程線比較
其余年份比較圖見附圖。另因計算時間尺度長，僅附1989年逐日計算數據，見附表。
三、精度統計與誤差分析
3.1 精度統計
從表2-3中可見，年產流量絕對誤差小于100mm的有7次，占總數的87.5%，絕對誤差小于45mm的有4次，占總數的50%；所有年份產流量的相對誤差均小于10%；確定性系數最大為0.9006，最小為0.6872，平均為0.8317。
精度統計表明，率定的模型參數基本上是合理的。
3.2 誤差分析
影響流域降雨徑流過程的因素很多，三水源新安江模型的結構和參數能夠反映濕潤地區降雨徑流過程的主要規律和特點，因而能獲得較好的精度。但是，模型本身以及模型計算中有許多的概化，會造成誤差。造成呈村流域日模型預報方案誤差來源主要有以下三個方面：
3.2.1 模型參數的影響
模型參數是根據輸入，通過模型計算輸出，再將輸出過程與實測過程進行比較，用人工調試的方法進行優化的，所率定出的參數可能不是最優的。
3.2.2 人類活動的影響
隨著經濟建設的發展，人類活動的影響加劇，流域內可能建有塘壩、橡膠壩等蓄水水利工程。這些水利工程沒有固定的調度規則和調洪方式，干旱季節或汛初時，流域內的水利工程均會蓄水，會使實測徑流量偏小。
3.2.3 資料代表性的影響
本次課程設計采用的樣本僅有8年實測水文系列資料，由于水文自然規律的隨時間變化，使觀測到的水文氣象資料代表性不夠，難免造成預報方案的誤差。
四、參考文獻
[1] 包為民.水文預報.北京：中國水利水電出版社，2017.
[2] 李巧玲.水文預報課程設計任務書.南京：河海大學水文水資源學院，2019．
[3] 翟家瑞.常用水文預報算法和計算程序. 河南：黃河水利出版社，1995
五、程序代碼
clear;
% data=xlsread('呈村流域資料.xls');
data=xlsread('1995參數檢驗.xls');
data1=xlsread('呈村流域資料.xls',2); %載入原始資料數據
QZ=data(:,2); %QZ表示實測流量
%%%%%%%%%%%%%需率定的參數%%%%%%%%%%%%%
par=load('input.txt');
K=par(1); %蒸散發折算系數 范圍：0.8～1.2
SM=par(2); %自由水容量 范圍：10～30
KG=par(3); %地下水出流系數 范圍: 0.01～0.69
CG=par(4); %地下水消退系數 范圍：0.98～0.998
CI=par(5); %壤中流消退系數 范圍：0.01～0.9
CS=par(6); %河網水流消退系數 范圍：0.1～0.9
L=1; %河網匯流滯時 范圍：0～2
%%%%%%%%%%%%%需率定的參數%%%%%%%%%%%%%
E0=data(:,3); %E0表示呈村流域器皿蒸發能力過程,位于原始數據表第三列
EP=K*E0; %計算流域蒸發能力
F=data1(:,3); %讀取各子流域面積
Fr=data1(:,4); %讀取各子流域面積比例
U=F/(3.6*24); %U為單位轉換系數
WM=120;WUM=20;WLM=60;WDM=WM-WUM-WLM;C=0.18;B=0.3;FC=2.5;FR0=0;PTT=0;%PTT表示流域年降雨量 %常數
EX=1.5;IM=0.01;KI=0.7-KG;WMM=(1+B)*WM;SMM=(1+EX)*SM;SZ=length(E0); %常數
S1=zeros(SZ,1);S1(1)=SM; %為S1預分配內存；S1表示本時段初產流面積上的平均自由水深
FR=zeros(SZ,1); %為FR預分配內存；FR表示本時段產流面積比例
EU=zeros(SZ,1); %為EU預分配內存
EL=zeros(SZ,1); %為EL預分配內存
ED=zeros(SZ,1); %為ED預分配內存
E=zeros(SZ,1); %為E預分配內存
WU=zeros(SZ,1); %為WU預分配內存
WL=zeros(SZ,1); %為WL預分配內存
WD=zeros(SZ,1); %為WD預分配內存
W=zeros(SZ,1); %為W預分配內存
R=zeros(SZ,1); %為R預分配內存
RG2=zeros(SZ,1); %為RG2預分配內存;RG2表示二水源地下徑流
RD2=zeros(SZ,1); %為RD2預分配內存；RD2表示二水源地面徑流
RS3=zeros(SZ,1); %為RS3預分配內存；RS3表示三水源地面徑流
RI3=zeros(SZ,1); %為RI3預分配內存；RI3表示三水源壤中水徑流
RG3=zeros(SZ,1); %為RG3預分配內存；RG3表示三水源地下徑流
PE=zeros(SZ,1); %為PE預分配內存；PE表示凈雨量
QI3=zeros(SZ,1); %為QI3預分配內存；QI3表示三分水源壤中流匯流
QG3=zeros(SZ,1); %為QG3預分配內存；QG3表示三分水源地下徑流匯流
Q3=zeros(SZ,1); %為Q3預分配內存；Q3表示單元面積河網匯流
Q=zeros(SZ,1); %為Q預分配內存；Q表示呈村流域出口斷面流量過程
TTT=zeros(SZ,18); %為TTT預分配內存；
WU(1)=0; %流域上層初始土壤含水量 范圍：0～20
WL(1)=50;WD(1)=WDM; %WL、WD都表示本時段初相應土層土壤含水量
W(1)=WU(1)+WL(1)+WD(1); %W表示本時段初土壤含水量
%%%%%%%%%%%%%計算10個子流域%%%%%%%%%%%%%
for J=1:10
P=(1-IM)*data(:,J+3); %P表示某子流域降雨過程，data數組中第4列到第13列分別為
% 呈村、汪村、樟源口、棣甸、董坑塢、用功城、左龍、馮村、田里、大連的降雨資料
RB=IM*data(:,J+3); %RB表示降落在不透水面上的降雨量
PTT=PTT+sum(data(:,J+3))*Fr(J);
for I=1:SZ
%%%%%%%%%%%%%三層蒸散發模型%%%%%%%%%%%%%
if WU(I)+P(I)>=EP(I)
EU(I)=EP(I);EL(I)=0;ED(I)=0;
else if WL(I)>=C*WLM
EU(I)=WU(I)+P(I);EL(I)=(EP(I)-EU(I))*WL(I)/WLM;ED(I)=0;
else if WL(I)>=C*(EP(I)-(WU(I)+P(I)))
EU(I)=WU(I)+P(I);EL(I)=C*(EP(I)-EU(I));ED(I)=0;
else
EU(I)=WU(I)+P(I);EL(I)=WL(I);ED(I)=C*(EP(I)-EU(I))-EL(I);
end
end
end
E(I)=EU(I)+EL(I)+ED(I);PE(I)=P(I)-E(I);
%%%%%%%%%%%%%三層蒸散發模型%%%%%%%%%%%%%
%%%%%%%%%%%%%計算總徑流深%%%%%%%%%%%%%
a=WMM*(1-(1-W(I)/WM)^(1/(B+1))); %a表示流域初始土壤含水量最大值
if PE(I)<=0
R(I)=0;
else if a+PE(I)<=WMM
R(I)=PE(I)+W(I)-WM+WM*(1-(PE(I)+a)/WMM)^(B+1);
else
R(I)=PE(I)+W(I)-WM;
end
end
%%%%%%%%%%%%%計算總徑流深%%%%%%%%%%%%%
%%%%%%%%%%%%%計算下一時段土壤含水量%%%%%%%%%%%%%
WU(I+1)=WU(I)+P(I)-EU(I)-R(I);
WL(I+1)=WL(I)-EL(I);
WD(I+1)=WD(I)-ED(I);
if WU(I+1)>WUM
WL(I+1)=WL(I)-EL(I)+WU(I+1)-WUM;
WU(I+1)=WUM;
end
if WL(I+1)>WLM
WD(I+1)=WD(I)-ED(I)+WL(I+1)-WLM;
WL(I+1)=WLM;
end
if WD(I+1)>WDM
WD(I+1)=WDM;
end
W(I+1)=WU(I+1)+WL(I+1)+WD(I+1);
%%%%%%%%%%%%%計算下一時段土壤含水量%%%%%%%%%%%%%
%%%%%%%%%%%%%計算產流面積FR=R/PE%%%%%%%%%%%%%
if PE(I)==0
FR(I)=0;
else if R(I)/PE(I)>1
FR(I)=1;
else
FR(I)=R(I)/PE(I);
end
end
%%%%%%%%%%%%%計算產流面積FR=R/PE%%%%%%%%%%%%%
%%%%%%%%%%%%%二分水源計算%%%%%%%%%%%%%
if PE(I)<=FC
RG2(I)=R(I);RD2(I)=0;
else
RG2(I)=FC*FR(I);RD2(I)=(PE(I)-FC)*FR(I);
end
%%%%%%%%%%%%%二分水源計算%%%%%%%%%%%%%
%%%%%%%%%%%%%三分水源計算%%%%%%%%%%%%%
%%%%%%%%%%%%%第一時段數據計算%%%%%%%%%%%%%
if I==1
AU=SMM*(1-(1-(S1(I)*FR0/FR(I))/SM)^(1/(1+EX)));
if PE(I)<=0
RS3(I)=0;RI3(I)=0;RG3(I)=0;S1(I+1)=S1(I)*(1-KI-KG);
else if PE(I)+AU
數據文件
模型運行所需文件請從新安江模型數據下載。
總結
以上是生活随笔為你收集整理的matlab液体湿润模拟,【水文模型】01 三水源新安江模型的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： VS 2010 OpenGL 配置与实例
- 下一篇： 数字图像处理（MATLAB）（第三版）