基于SoilGrids的柵格新安江模型參數(shù)空間分布估算

童冰星，李致家，姚 成

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

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

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

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

1 計(jì)算方法

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

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

(1)

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

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

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

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

在數(shù)字高程模型(DEM)基礎(chǔ)上，GXM模型采用D8法[18]識(shí)別柵格單元流向，確定柵格匯流演算次序，提取水系網(wǎng)絡(luò)，區(qū)分河道和坡地柵格單元。模型基于蓄滿產(chǎn)流理論計(jì)算柵格單元產(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分別為壤中和地下徑流進(jìn)入坡面溝道延遲時(shí)間，h，采用式(5)計(jì)算。

(5)

2 實(shí)例分析

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

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

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

圖2 陳河流域高程與站點(diǎn)分布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é)果

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

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

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

2.3 土壤水空間分布

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

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

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

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

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

圖6 2003090319號(hào)洪水過程的土壤水飽和度空間動(dòng)態(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ù)值，對(duì)2003090319號(hào)洪水進(jìn)行模擬，分別統(tǒng)計(jì)漲洪、洪峰、落洪和洪尾段的ΔR和CD，量化參數(shù)SM對(duì)不同洪水階段的影響(圖8)。

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

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

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

大量研究和野外觀測實(shí)驗(yàn)表明，SM空間分布的合理與否對(duì)土壤含水量等水文要素時(shí)空分布動(dòng)態(tài)變化模擬有重要影響[20]。為進(jìn)一步驗(yàn)證參數(shù)SM空間分布的合理性，根據(jù)陳河流域支流多分布在南側(cè)的地理特點(diǎn)，本文以經(jīng)緯度為(107°48′ E，33°48′ N)和(108°18′ E，33°52′ N)的2個(gè)端點(diǎn)，作1條跨多支流的取樣線段(圖9)。在取樣線段上以250 m為間隔設(shè)置樣本點(diǎn)，統(tǒng)計(jì)樣本點(diǎn)處SM值與高程、距離河道遠(yuǎn)近(G)和土壤有機(jī)質(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值減小，而山脊附近有機(jī)質(zhì)含量較高，SM值較大。

4 結(jié) 論

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

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

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

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