童冰星，李致家，姚 成

(河海大學水文水資源學院，江蘇 南京 210098)

水文模型參數(shù)是對流域下墊面特征的量化，其空間分布估算需考慮下墊面特征空間異質(zhì)性，是分布式水文模型研究和發(fā)展的關(guān)鍵環(huán)節(jié)和難點問題[1- 2]。中國集水面積200～3 000 km2的中小流域約9 000個[3]，廣泛分布在下墊面地理特征觀測薄弱的偏遠山區(qū)，不利于參數(shù)空間分布的合理估算，制約分布式水文模型在中小流域的構(gòu)建與應用。因此，基于高分辨率地球物理數(shù)據(jù)，開發(fā)分布式水文模型參數(shù)空間分布估算方法，對于中小流域降雨徑流過程精細化模擬，開展洪水風險預警與影響預報、防旱減災和水土保持等工作具有重要意義[4]。

國內(nèi)外已有研究學者提出一些分布式水文模型參數(shù)空間分布估算方法[5- 7]。Todini等研發(fā)的TOPKAPI模型[8]根據(jù)世界土壤數(shù)據(jù)庫(Harmonized World Soil Database,HWSD)劃分土壤類型，估算土壤飽和滲透系數(shù)等參數(shù)；Arnold等開發(fā)的SWAT模型[9]采用土地利用、植被覆蓋數(shù)據(jù)對土壤濕容重等參數(shù)進行估計；Yao等[10]提出的柵格新安江模型(Grid- Xin′anjiang Model，GXM)基于土壤分類推求張力水和自由水蓄水容量。位于秦嶺北麓的陳河流域?qū)儆谥行×饔?，土壤類型空間異質(zhì)性較小，然而砂粒、粉粒含量等土壤質(zhì)地差異可能導致土壤物性參數(shù)不同，影響蓄水容量等參數(shù)的空間分布[11- 12]。國際土壤信息中心(ISRIC)研發(fā)的新版數(shù)字土壤制圖系統(tǒng)(SoilGrids250mTM，V2.0，以下簡稱SoilGrids)在2020年5月開始運行[13]。SoilGrids基于WoSIS數(shù)據(jù)庫23萬個土壤剖面觀測數(shù)據(jù)和一系列環(huán)境協(xié)變量擬合土壤脊線預測模型，采用機器學習方法繪制全球范圍內(nèi)250 m空間分辨率、6個標準深度間隔(5 cm、15 cm、30 cm、60 cm、100 cm和200 cm)的土壤特性圖(重度、黏粒含量和砂粒含量等)。相較原始版本，新版SoilGrids更新補充土壤剖面觀測數(shù)據(jù)，是目前全球范圍內(nèi)較為完善的數(shù)字土壤制圖系統(tǒng)。因此，結(jié)合新版SoilGrids提供的高分辨率數(shù)據(jù)，量化土壤質(zhì)地等下墊面特征的空間異質(zhì)性，合理估算GXM等分布式水文模型的參數(shù)空間分布，是中小流域降雨徑流過程精細化模擬研究的重要突破方向。

本文基于新版SoilGrids構(gòu)建GXM模型參數(shù)空間分布估算方案，對陜西省陳河流域16場洪水過程進行模擬，得到降雨徑流過程自由水含量動態(tài)空間分布，與實測徑流過程和新安江模型[14]的計算結(jié)果進行對比，量化GXM模型關(guān)鍵參數(shù)自由水蓄水容量空間分布特征，開展基于洪水過程劃分的參數(shù)敏感性分析。

1 計算方法

1.1 土壤物理性質(zhì)參數(shù)空間分布估算

土壤含水量隨深度增加變化幅度逐漸減小，自上而下大致可分為活躍層和相對穩(wěn)定層，通常水分活躍層有機質(zhì)含量較高[15- 16]。本文基于SoilGrids提供的垂向剖面有機質(zhì)含量變化趨勢，分層概化土壤，估算活躍層厚度空間分布。

(1)

式中：Lh為土壤水分活躍層厚度，mm；Lmin為水分活躍層最小可能厚度，mm；Lmax為水分活躍層最大可能厚度，mm；La為土壤包氣帶厚度，mm；LM為流域土壤包氣帶最大厚度，mm，由SoilGrids提供。Lmin和Lmax分別可設(shè)置為有機質(zhì)含量占總量比例為α和β處的土壤深度(圖1)，α和β取值參考流域下墊面實地調(diào)查和洪水模擬經(jīng)驗。

圖1 土壤垂向剖面有機質(zhì)變化示意Fig.1 Change of the organic matter on the soil vertical profile

土壤物理性質(zhì)參數(shù)主要包括凋萎含水量(W)、田間持水量(F)、飽和含水量(B)和飽和滲透系數(shù)(K)等，參數(shù)估算依據(jù)文獻[17]。

1.2 GXM模型參數(shù)估算

在數(shù)字高程模型(DEM)基礎(chǔ)上，GXM模型采用D8法[18]識別柵格單元流向，確定柵格匯流演算次序，提取水系網(wǎng)絡，區(qū)分河道和坡地柵格單元。模型基于蓄滿產(chǎn)流理論計算柵格單元產(chǎn)流量，采用自由水蓄水庫結(jié)構(gòu)劃分地表水、壤中流與地下水3種徑流成分，根據(jù)柵格間匯流演算次序，將徑流演算至流域出口。GXM模型主要參數(shù)如表1所示。

表1 GXM模型主要參數(shù)

SM、WM(=WUM+WLM+WDM)參數(shù)與土壤厚度和土壤物理性質(zhì)參數(shù)相關(guān)：

SM=Lh(B-F)WM=La(F-W)

(2)

KI、KG參數(shù)反映土壤排水能力：

(3)

CI、CG參數(shù)量化徑流消退快慢，主要與坡段長度和坡度相關(guān)：

(4)

式中：S為地形坡度；Ku和Km分別為水分活躍層和包氣帶的飽和滲透系數(shù)，m/h；H為坡段長度，m；Ti和Td分別為壤中和地下徑流進入坡面溝道延遲時間，h，采用式(5)計算。

(5)

2 實例分析

2.1 流域簡介與參數(shù)空間分布

陳河流域位于秦嶺北麓，流域面積約1 350 km2，流域內(nèi)主河道自西南向東北匯入渭河，其余支流呈扇形匯入主河道。該區(qū)域大部分為山區(qū)，海拔高程630～3 747 m，受大陸性季風氣候影響，多年平均降水量700～900 mm。觀測站網(wǎng)包括流域出口處的陳河水文站以及流域內(nèi)麥場、板房子等9個雨量觀測站(圖2)，該觀測站網(wǎng)收集的2003—2012年降雨、徑流觀測數(shù)據(jù)用于洪水模擬。

根據(jù)SoilGrids可得到土壤包氣帶厚度分布。陳河流域包氣帶最大厚度、分層比例α和β取值分別為1 419 mm、25%和60%，估算得到GXM模型參數(shù)空間分布如圖3。

圖2 陳河流域高程與站點分布Fig.2 Observation stations and Digital Elevation Model of the Chenhe watershed

圖3 GXM模型水文參數(shù)在陳河流域的空間分布Fig.3 Estimated spatial parameters of the GXM model in the Chenhe watershed

2.2 洪水模擬結(jié)果

對陳河流域2003—2012年的16場洪水過程進行模擬，統(tǒng)計洪量相對誤差(ΔR)、洪峰相對誤差(ΔP)、峰現(xiàn)時間誤差(ΔT)和確定性系數(shù)(CD)等指標，并與新安江模型(XAJ)的模擬結(jié)果對比分析(圖4)。

GXM模擬的洪量和洪峰相對誤差水平為15.3%和14.9%。新安江模型模擬的洪量和洪峰相對誤差水平為16.8%和15.4%。新安江模型和GXM模擬結(jié)果的確定性系數(shù)均值分別為0.79和0.76。相較于新安江模型，雖然GXM模擬結(jié)果的確定性系數(shù)略小，但洪峰和洪量模擬精度更高，且GXM模擬的陳河流域峰現(xiàn)時間誤差水平降低約0.31 h。GXM模擬的部分洪水過程如圖5所示。

圖4 洪水模擬的誤差評價指標結(jié)果Fig.4 Result of error evaluation indexes for floods simulation

2.3 土壤水空間分布

與新安江模型相比，GXM能夠?qū)ν寥浪柡投鹊人囊氐膭討B(tài)空間分布進行較合理地模擬。以陳河流域2003090319號洪水為例(圖5(a))，分別在洪水開始后25 h、45 h、65 h、85 h、105 h和125 h輸出GXM模擬的土壤水飽和度空間分布結(jié)果(圖6)?？梢钥闯觯跏纪寥浪柡投容^低(圖6(a))，隨著降雨持續(xù)，土壤水飽和度逐漸增加(圖6(b))，在洪峰出現(xiàn)時接近飽和(圖6(c))；降雨結(jié)束后土壤水飽和度逐漸減少(圖6(d)、圖6(e))，但仍高于初始時刻的土壤水飽和度(圖6(f))。

圖5 GXM模型洪水模擬結(jié)果Fig.5 Floods simulation results using the GXM model

3 參數(shù)SM的量化分析

3.1 基于洪水過程劃分的SM敏感性分析

研究不同洪水階段參數(shù)SM敏感性，有助于進一步開展參數(shù)不確定性分析，為實時洪水預報參數(shù)動態(tài)調(diào)整提供參考，促進降雨徑流過程精細化模擬[19]。本文提出五段法，將洪水劃分為初始、漲洪、洪峰、落洪和洪尾階段，研究參數(shù)SM在不同階段對模擬結(jié)果的影響。以陳河流域2003090319號洪水為例，將時間t作為自變量計算洪水過程一階導數(shù)(Q′=?f(t)/?t)，采用3 h為周期加權(quán)移動平均后依次確定：一階導數(shù)值增大位置(A)、最大值(B)、最小值(C)和恢復穩(wěn)定位置(D)，將洪水劃分為初始、漲洪、洪峰、落洪和洪尾階段(圖7)。

圖6 2003090319號洪水過程的土壤水飽和度空間動態(tài)變化Fig.6 Spatial dynamics of soil moisture during No.2003090319 flood

圖7 洪水階段劃分Fig.7 Phases division of flood

在1～30 mm內(nèi)，采用0.1 mm為步長設(shè)置SM參數(shù)值，對2003090319號洪水進行模擬，分別統(tǒng)計漲洪、洪峰、落洪和洪尾段的ΔR和CD，量化參數(shù)SM對不同洪水階段的影響(圖8)。

圖8 參數(shù)SM的敏感性Fig.8 Sensitivity of the parameter SM

對因參數(shù)SM變化導致的ΔR和CD在不同洪水階段的方差進行統(tǒng)計，作為衡量參數(shù)敏感性的指標。漲洪段洪量相對誤差受SM參數(shù)影響較大，敏感性為0.37，明顯高于洪峰、落洪和洪尾段(敏感性分別為0.11、0.04和0.003)；SM參數(shù)對洪峰段和漲洪段的確定性系數(shù)影響較大，敏感性分別為3.76和3.05，對落洪段和洪尾段影響較小,敏感性分別為0.03和0.05。說明SM參數(shù)對洪峰和漲洪過程的確定性系數(shù)及漲洪段的洪量相對誤差影響較大，對退水過程影響小。

3.2 SM空間分布特征的量化分析

大量研究和野外觀測實驗表明，SM空間分布的合理與否對土壤含水量等水文要素時空分布動態(tài)變化模擬有重要影響[20]。為進一步驗證參數(shù)SM空間分布的合理性，根據(jù)陳河流域支流多分布在南側(cè)的地理特點，本文以經(jīng)緯度為(107°48′ E，33°48′ N)和(108°18′ E，33°52′ N)的2個端點，作1條跨多支流的取樣線段(圖9)。在取樣線段上以250 m為間隔設(shè)置樣本點，統(tǒng)計樣本點處SM值與高程、距離河道遠近(G)和土壤有機質(zhì)含量(M)等下墊面因子的定量關(guān)系(圖10)。

圖9 參數(shù)SM與下墊面因子(高程、G和M)的空間分布Fig.9 Spatial distribution of the parameter SM and the geographical characters (DEM,G and M)

圖10 取樣線段上參數(shù)SM與下墊面因子(高程、G和M)的相關(guān)性Fig.10 Correlation between the parameter SM and the geographical characters (DEM,G and M) in sampling line

由圖9和圖10可知，SM值隨高程、G和M的增加，先減小后增大，相關(guān)性分別為0.48、0.44和0.29。水流攜帶的泥沙在陳河流域谷底沉積，河道附近SM值較大。坡段中部水流流速較快，土壤侵蝕頻發(fā)，使土壤厚度和SM值減小，而山脊附近有機質(zhì)含量較高，SM值較大。

4 結(jié) 論

本文基于新版全球數(shù)字土壤制圖系統(tǒng)(SoilGrids)，構(gòu)建柵格新安江模型的參數(shù)估算方案。對陜西省陳河流域2003—2012年16場洪水過程進行模擬，并與新安江模型的模擬結(jié)果和實測徑流比較，得到以下主要結(jié)論：

(1) 結(jié)合參數(shù)估算方案，柵格新安江模型模擬的峰現(xiàn)時間誤差水平降低約0.31 h，洪峰和洪量模擬精度較高，模型能夠?qū)ν寥浪柡投鹊人囊氐膭討B(tài)空間分布進行較合理地模擬。

(2) 自由水蓄水容量對洪峰和漲洪過程的確定性系數(shù)以及漲洪段的洪量相對誤差影響較大，對退水過程影響小。

(3) 自由水蓄水容量在陳河流域河谷和山脊附近較大，坡段中部較小。