[發(fā)明專利]一種基于自然子流域的通用流域水循環(huán)模擬計(jì)算方法

有效

• 申請(qǐng)?zhí)枺?/td>: 202011463721.8
• 申請(qǐng)日: 2020-12-11
• 公開（公告）號(hào): CN112651189B
• 公開（公告）日: 2023-03-10
• 發(fā)明（設(shè)計(jì)）人: 楊開斌;夏建榮;韓兵;劉陽(yáng)容;盧鵬;周鵬程
• 申請(qǐng)（專利權(quán)）人: 中國(guó)電建集團(tuán)昆明勘測(cè)設(shè)計(jì)研究院有限公司
• 主分類號(hào): G06F30/28
• 分類號(hào): G06F30/28;G06F113/08;G06F119/14
• 代理公司: 昆明盛鼎宏圖知識(shí)產(chǎn)權(quán)代理事務(wù)所(特殊普通合伙) 53203
• 代理人: 王輝

1.一種基于自然子流域的通用流域水循環(huán)模擬計(jì)算方法，其特征在于，包括以下步驟：

步驟一：流域水文分析及子流域河道匯流計(jì)算順序生成

基于DEM地形資料對(duì)研究流域進(jìn)行水文分析，將研究流域劃分為若干個(gè)子流域，明確各子流域相應(yīng)下游子流域編號(hào)，采用拓?fù)潢P(guān)系模型分析出子流域河道匯流計(jì)算順序；

步驟二：計(jì)算子流域面雨量過程、選擇流域代表性蒸發(fā)過程

基于研究流域內(nèi)及周邊雨量站降雨資料，采用高容錯(cuò)面雨量計(jì)算方法計(jì)算得到各子流域面雨量過程；選取流域內(nèi)代表站蒸發(fā)過程作為流域代表性蒸發(fā)過程：

對(duì)各自然子流域劃分正方形網(wǎng)格，各網(wǎng)格以有效面積占自然子流域面積的比例為權(quán)重，自然子流域面雨量采用下式計(jì)算得到；

式中：n為子流域內(nèi)網(wǎng)格個(gè)數(shù)；A_{網(wǎng)格i}為第i個(gè)網(wǎng)格有效面積；A_子流域為子流域面積；P_{網(wǎng)格i}為第i個(gè)網(wǎng)格雨量；

對(duì)于單一網(wǎng)格的雨量推求，通過排序算法搜索最接近網(wǎng)格中心點(diǎn)的N個(gè)雨量站點(diǎn)，根據(jù)站點(diǎn)雨量值判斷、挑選出資料完備、準(zhǔn)確的m個(gè)站點(diǎn)，m小于N，基于m個(gè)站點(diǎn)的雨量資料采用距離平方倒數(shù)法插值得，距離平方倒數(shù)法公式如下所示：

式中：為m個(gè)站點(diǎn)中第j個(gè)站點(diǎn)的雨量；d_ij為第j個(gè)站點(diǎn)距第i個(gè)網(wǎng)格中心點(diǎn)的距離；lat為緯度；lon為經(jīng)度；

步驟三：根據(jù)氣候及下墊面特性，選擇各子流域產(chǎn)匯流計(jì)算方法

根據(jù)地形資料識(shí)別研究流域所屬地形地貌；根據(jù)柯本氣候分類識(shí)別研究流域內(nèi)氣候類別分布；繪制研究流域下墊面特性：土壤類別分布、土地覆被類別分布及坡度類別分布圖；結(jié)合氣候類別分布及下墊面特性分布，設(shè)定各子流域選取的產(chǎn)匯流計(jì)算方法；

步驟三：產(chǎn)匯流計(jì)算方法針對(duì)的子流域不同水文氣候條件為：1)氣候干旱，土層厚度大，超滲產(chǎn)流顯著；或2)氣候濕潤(rùn)，地表易蓄滿產(chǎn)流；或3)子流域缺乏代表性的蒸發(fā)資料輸入；或4)子流域高寒高海拔，存在融雪徑流；

針對(duì)子流域不同水文氣候條件為：1)氣候干旱，土層厚度大，超滲產(chǎn)流顯著；算法采用霍頓下滲曲線進(jìn)行超滲產(chǎn)流計(jì)算；

霍頓下滲曲線的相關(guān)公式為：

f＝f_c+(f₀-f_c)e^-kt (4)

聯(lián)立式(4)和式(5)，有：

式中：f為子流域時(shí)段內(nèi)平均下滲率，mm/h；f₀、f_c分別為子流域平均最大、最小下滲能力，mm/h；t為歷時(shí)，h；k為下滲能力衰減系數(shù)，h^-1；W為土壤含水量，mm；

需用迭代求出f～W的關(guān)系；迭代過程為：以T＝W/f₀作為t的第一次近似值，由式(5)計(jì)算得到W的第一個(gè)近似值ST，如果|ST-W|允許誤差e，則由式(4)計(jì)算出f的第一個(gè)近似值U，然后T＝T+(W-ST)/U，迭代多次后，直到|ST-W|≤e，求得所需的f值；

針對(duì)子流域不同水文氣候條件為：1)氣候干旱，土層厚度大，超滲產(chǎn)流顯著；算法采用菲利普下滲曲線進(jìn)行超滲產(chǎn)流計(jì)算；

菲利普下滲曲線的相關(guān)公式為：

式中：B和A為兩個(gè)待定參數(shù)；其他參數(shù)含義如上；

菲利普下滲曲線，在給定一組系數(shù)A、B參數(shù)值，求得f～W的關(guān)系；

在得到f～W關(guān)系后，根據(jù)下式計(jì)算超滲產(chǎn)流量：

F＝f△t (11)

PE＝P-E (12)

f_mm＝f(1+BX) (13)

式中：R_IE為超滲產(chǎn)流量，mm；F為時(shí)段下滲量，mm；f_mm為流域平均下滲能力為f時(shí)流域最大的點(diǎn)下滲能力；BX為指數(shù)系數(shù)；

針對(duì)子流域不同水文氣候條件為：2)氣候濕潤(rùn)，地表易蓄滿產(chǎn)流；算法采用基于地形指數(shù)的蓄滿產(chǎn)流計(jì)算方法；

①蒸發(fā)量計(jì)算

式中：E_a,i為點(diǎn)i處實(shí)際蒸發(fā)量，m；E_P為蒸散發(fā)能力，m；S_rz,i為植被根系區(qū)缺水量，m；S_rmax,i為植被根系區(qū)最大蓄水容量，m；

②產(chǎn)流量計(jì)算

式中：a_i為點(diǎn)i處單寬集水面積，m²；tanβ_i為點(diǎn)i處的地表坡度；z_i為點(diǎn)i處地下水距地表深度，m；z為飽和地下水水面平均深度，m；S_zm為非飽和區(qū)最大蓄水深度，m；

若zi為負(fù)值，飽和地下水將漫出地面，形成地表徑流；

點(diǎn)i處下滲率計(jì)算公式為：

式中：S_uz,i為點(diǎn)i處非飽和區(qū)土壤含水量，m；SD_i為非飽和區(qū)土壤蓄水能力，m；t_d為時(shí)間參數(shù)，h；

整個(gè)流域的下滲率為：

式中：A_i為地形指數(shù)數(shù)值相同的各處面積之和，m²；

Q_b＝AT₀exp(-λ)exp(-z/S_zm) (19)

式中：T₀為飽和導(dǎo)水率，m²/h；

飽和地下水水面平均深度z的計(jì)算公式為：

針對(duì)子流域不同水文氣候條件為：2)氣候濕潤(rùn)，地表易蓄滿產(chǎn)流；算法采用蓄水容量～面積曲線的蓄滿產(chǎn)流計(jì)算方法；

①蒸發(fā)量計(jì)算

當(dāng)WU+P≥E_P時(shí)

E_U＝E_P E_L＝0 E_D＝0 (21)

當(dāng)WU+PE_P，WL≥C·WLM時(shí)

E_U＝WU+P E_L＝(E_P-E_U)WL/WLM E_D＝0 (22)

當(dāng)WU+PE_P，C(E_P-E_U)≤WLC·WLM時(shí)

E_U＝WU+P E_L＝C(E_P-E_U) E_D＝0 (23)

當(dāng)WU+PE_P，WLC(E_P-E_U)時(shí)

E_U＝WU+P E_L＝WL E_D＝C(E_P-E_U)-E_L (24)

E＝E_U+E_L+E_D (25)

式中：E_P為蒸散發(fā)能力；P為降雨量；WL為下層土壤含水量；WU為上層土壤含水量；WLM為下層土壤含水容量；C為蒸發(fā)擴(kuò)散系數(shù)；E_U為上層土壤蒸發(fā)量；E_L為下層土壤蒸發(fā)量；E_D為深層土壤蒸發(fā)量；E為總蒸發(fā)量；

②蓄滿產(chǎn)流量計(jì)算

當(dāng)a+PE≤WMM流域部分面積產(chǎn)流：

當(dāng)a+PEWMM全流域產(chǎn)流：

R＝PE-(WM-W) (27)

其中：

式中：R為蓄滿產(chǎn)流量，mm；a為與流域初始平均蓄水量W相應(yīng)的流域最大蓄水量，mm；b為拋物線指數(shù)；WM為流域平均蓄水容量；WMM為流域最大蓄水容量；PE為扣除蒸發(fā)后的凈雨；

③產(chǎn)流分配計(jì)算

當(dāng)PE+AUS_mm：

當(dāng)PE+AU≥S_mm：

其中：

FR＝R/PE (32)

時(shí)段自由水蓄水量為：

壤中產(chǎn)流量和地下產(chǎn)流量分別為：

RI＝KI·S·FR (34)

RG＝KG·S·FR (35)

式中：S_mm為流域最大自由水蓄水容量，mm；S_m為流域平均自由水蓄水容量，mm；EX為拋物線指數(shù)；S₁為時(shí)段初始流域平均自由水蓄水量，mm；AU為與S₁相應(yīng)的流域最大自由水蓄水容量；FR₁和、FR分別為上一時(shí)段及本時(shí)段的產(chǎn)流面積比例；PE為扣除蒸發(fā)后的凈雨；R為蓄滿產(chǎn)流量；RS為地表產(chǎn)流量；RI為壤中產(chǎn)流量；RG為地下產(chǎn)流量；

針對(duì)子流域不同水文氣候條件為：3)子流域缺乏代表性的蒸發(fā)資料輸入；算法推薦采用基于增益因子的產(chǎn)流計(jì)算方法：

在基于增益因子的產(chǎn)流計(jì)算方法中，產(chǎn)流R表現(xiàn)為降雨P(guān)和增益G的乘積：

R(t)＝G(t)P(t) (36)

增益G與前期土壤含水量W有關(guān)，為：

經(jīng)泰勒展開后：

G(t)＝g₁+g₂W(t) (38)

R(t)＝g₁P(t)+g₂W(t)P(t) (39)；

針對(duì)子流域不同水文氣候條件為：4)子流域高寒高海拔，存在融雪徑流；算法推薦采用基于增益因子的產(chǎn)流計(jì)算方法：

算法采用度日因子計(jì)算方法進(jìn)行積融雪計(jì)算；

JR_snow＝D(T_t-T_c) (40)

式中：JR_snow為正時(shí)表示為融雪量，為負(fù)值時(shí)表示為積雪量，mm；D為度日因子，mm/d；T_t為日平均氣溫；T_c為臨界氣溫，設(shè)為0℃；

若時(shí)段降雨為P，前期積雪深為S，則計(jì)算時(shí)段末積雪深為S-JR_snow限制不小于0，時(shí)段凈雨為P+JR_snow(JR_snowS)或P+S(JR_snow≥S)；

此外，對(duì)于地表徑流匯流計(jì)算，采用線性水庫(kù)或地貌單位線方法；地下徑流匯流計(jì)算，采用線性水庫(kù)計(jì)算方法；

1)線性水庫(kù)法

線性水庫(kù)法公式為：

Q_t+1＝R_t+1(1-C)U+Q_t*C (41)

U＝AREA/(△t*3.6) (42)

式中：Q_t+1、Q_t分別為t+1、t時(shí)刻的流量，m³/s；R_t+1為t+1時(shí)刻產(chǎn)流量，mm；C為消退系數(shù)；U為單位轉(zhuǎn)換系數(shù)；AREA為流域面積，km²；△t為時(shí)段長(zhǎng)，h；

2)地貌單位線

地貌單位線形式如下：

式中：N為反映流域調(diào)蓄能力的參數(shù)，K為線性水庫(kù)的蓄泄系數(shù)；Γ(N)為Γ函數(shù)，即

在參數(shù)N及參數(shù)K推求上，根據(jù)霍頓地貌幾何率：面積比、河長(zhǎng)比、分叉比推求：

式中：R_B,R_L,R_A—流域水系的分叉比、河長(zhǎng)比和面積比，基于斯特拉勒級(jí)別通過DEM資料求得；

使用如下關(guān)系式：

τ＝1-(1-λ)(1-ρ) (45)

其中：

得到如下關(guān)系式：

τ＝λ^1-mλ (47)

由式(45)和式(47)推導(dǎo)出：

用霍頓河長(zhǎng)定律，推導(dǎo)出：

以上式中：τ為凈雨質(zhì)點(diǎn)自河源至下游某一斷面的平均匯流時(shí)間與河源至流域出口斷面的平均匯流時(shí)間的比值；ρ為與河流長(zhǎng)度及河底比降有關(guān)的參數(shù)；n為自河源至下游某斷面的子河段數(shù)；N為自河源至流域出口斷面的子河段數(shù)；△l_j為自河源開始劃分的第j子河段長(zhǎng)；p_j為第j子河段的平均坡度；m為反映河道縱剖面特性的綜合參數(shù)；Ω為河系最高級(jí)別河流的級(jí)數(shù)；V_Ω為流域出口斷面的流速，一般由出口斷面洪水過程線的漲洪段的平均流速給定；α為流域形心至流域出口斷面的距離與流域長(zhǎng)度的比值；

m參數(shù)認(rèn)為是反映河道縱剖面特性的綜合參數(shù)；

根據(jù)霍頓河長(zhǎng)定律及比降定律，針對(duì)ρ參數(shù)計(jì)算，構(gòu)造河長(zhǎng)比降比R_LS這一概念，其為各級(jí)河流lp^-0.5值的平均比值，則有：

在此基礎(chǔ)上，聯(lián)立式(45)及式(47)，通過迭代求解，計(jì)算參數(shù)m；

步驟四：耦合拓?fù)潢P(guān)系模型、水文模型及水力模型，優(yōu)化產(chǎn)匯流參數(shù)

根據(jù)各子流域河道流經(jīng)區(qū)域的土地利用情況，確定河道糙率值；耦合拓?fù)潢P(guān)系模型、水文產(chǎn)匯流計(jì)算模型及水力學(xué)模型，逐時(shí)段進(jìn)行聯(lián)合計(jì)算，以模擬得到的研究流域出口斷面流量過程與實(shí)測(cè)流量過程的Nash效率系數(shù)作為參數(shù)優(yōu)化對(duì)象，優(yōu)化產(chǎn)匯流計(jì)算的參數(shù)，得到各子流域出流過程；

拓?fù)潢P(guān)系模型為自然子流域拓?fù)潢P(guān)系模型；算法為：步驟1)統(tǒng)計(jì)所有子流域的入度，入度即為匯流至某一子流域的相鄰子流域個(gè)數(shù)；步驟2)分離出對(duì)于入度為0的子流域，將入度為0的子流域匯流至的相鄰子流域的入度減1；步驟3)重復(fù)步驟2直至所有子流域被分離出來，完成子流域河道匯流計(jì)算順序的排序計(jì)算；

水力學(xué)模型采用一維非恒定流模型，對(duì)于有n個(gè)子流域的研究流域，需建立(n-1)/2個(gè)河段的單河道一維非恒定流模型；模型建立包括連續(xù)方程和動(dòng)量守恒方程；

1)連續(xù)方程

2)動(dòng)量守恒方程

在有限差分中采用四點(diǎn)隱式差分格式，其中n和j分別表示時(shí)間和空間的離散，n時(shí)段的水位和流量已知，待求的是n+1時(shí)段的水位和流量；差分形式為：

按此形式寫出各個(gè)項(xiàng)的差分并代入連續(xù)方程和動(dòng)量守恒方程；對(duì)于方程中的非線性項(xiàng)進(jìn)行適當(dāng)線性化，有：

經(jīng)過整理得：

A_1j△Q_j+B_1j△Z_j+C_1j△Q_j+1+D_1j△Z_j+1＝E_1j (57)

A_2j△Q_j+B_2j△Z_j+C_2j△Q_j+1+D_2j△Z_j+1＝E_2j (58)

式中A_1j、B_1j、C_1j、D_1j、E_1j、A_2j、B_2j、C_2j、D_2j、E_2j均為系數(shù)，如

上邊界采用上游子流域入流，下邊界采用相應(yīng)單河道最下游斷面的水位流量關(guān)系，由曼寧公式計(jì)算；

式中：Q為流量；n為糙率；A為斷面過水面積；R為斷面水力半徑；S為斷面位置水面比降，以底坡代替；

步驟五：若存在引水/調(diào)水過程，則分析對(duì)流域水循環(huán)的影響若特定子流域內(nèi)存在引水/調(diào)水過程，通過概化為該子流域出流過程的增加/減少，根據(jù)步驟四優(yōu)化得到的參數(shù)及其他已確定參數(shù)，計(jì)算各子流域出流過程，分析引水/調(diào)水過程對(duì)流域水循環(huán)的影響。