基于SWAT模型的沂蒙山水土流失治理區地表水與地下水轉化規律研究

譚秀翠，莊會波，季 妤

(1.山東農業大學水利土木工程學院，山東 泰安 271018；2.山東省水文局，濟南 250002)

地表水與地下水的水量轉化是水文循環的重要過程，地表水與地下水的相互作用及轉化關系是水文地質等領域研究的熱點和難點，分析和掌握其規律特征，對進行水資源分析評價、物質遷移和能量變化的相關研究具有重要意義[1]。目前，對地表水與地下水轉化的研究方法有：水文測量法[2]、基流分割法[3]、水量平衡法[4]、水文地球化學法[5,6]、溫度示蹤法[7,8]、數值模型法[9,10]等。

植被作為陸地生態環境中的重要組成部分，對水分具有吸收、駐留、緩釋等作用，是水文循環的重要環節。氣候變化和土地利用/覆被變化引起的各種時間尺度下流域下墊面植被特性的改變將對流域水文循環產生顯著的影響[11]。Loch[12]對比分析澳大利亞昆士蘭州南部植被覆蓋度分別為0、23%、37%、47%以及100%區域的產流和地表入滲，結果顯示植被覆蓋度增加將大大增加地表入滲，從而減小地表產流。石峰[13]利用LCM模型對兩種不同情景下次降雨徑流進行模擬，表明林草植被增加都明顯減小了流域產流量。史曉亮等[14]研究發現諾敏河流域植被覆蓋變化對徑流影響較小。成向榮等[15]構建6種不同植被覆蓋情景，研究該植被覆蓋類型變化對徑流的影響。目前，關于植被覆蓋的水文效應主要集中在地表產流、入滲，對于地表水與地下水轉化規律的研究較少。

山東省沂蒙山區地處山東省中南部、淮河流域北部，是淮河流域重要的水源區和生態屏障，也是淮河流域水土流失最嚴重的地區之一。2006年，水利部公告了42個國家級水土流失重點防治區，山東省的沂蒙山區被列為重點治理區，其空間分布位置見圖1。經過多年治理，2016年，沂蒙山治理區的監測區域中蒙陰縣和沂水縣，水土流失面積1 577.61 km2，占土地總面積的39.34%，與2013年調查結果相比，水土流失面積減少11.29 km2(中國水土保持公報，2016)。本文以山東省沂蒙山區蒙陰縣岸堤水庫控制流域為研究對象，分析氣候與植被覆蓋時空變化對地表水與地下水轉化的影響規律，研究成果可為沂蒙山區水資源的合理開發利用、水土流失治理提供理論依據和技術支撐。

圖1 沂蒙山水土流失治理區位置分布圖Fig.1 The location map of Yimeng mongtain soil and water loss control area

1 研究區域

蒙陰縣岸堤水庫控制流域地處沂蒙山區腹地，蒙山北麓，東汶河上游。流域地貌受構造、巖性、氣候、河流等內外營力作用的控制和影響，整個地形南北高，中間低，由北西向東南傾斜。流域地處中緯度地帶，屬暖溫帶半濕潤大陸性季風型氣候，具有四季分明，光照充足，雨量充沛，無霜期長等特點，多年平均降水量781.8 mm，降水量多集中在6-8月，占全年的65%，多年平均溫度13.4 ℃(蒙陰氣象站，1959-2018年數據)。汛期徑流量占年徑流量的82.3%，易形成春旱夏洪，約有2/3徑流量形成洪水流失。蒙陰縣降水量和徑流量總的空間分布特征為：南多北少，東多西少，山區大，平地少，從東南向西北遞減[16]。

蒙陰縣城東部東汶河北部，屬層間巖溶裂隙水，水位埋深2 m左右。縣城南部東汶河南側，淺層屬松散堆積物孔隙水、沖積粗砂礫石層潛水，水位埋深3 m左右。山南莊村東部屬層間巖溶裂隙水。東汶河入云蒙湖口西南部屬裂隙巖溶水。淺層地下水主要靠大氣降水補給，其水位深度隨季節變化幅度較大。

流域內分布有3個水文測站(1984-2016年逐日數據)，19個雨量站(1984-2016年逐日數據)，1個氣象站(1959-2018年逐日數據)，各站點空間分布見圖2。

圖2 研究流域分布圖Fig.2 The location map of research watershed

歸一化植被指數NDVI(normalized difference vegetation index)對植被的生長勢和生長量非常敏感，因此NDVI不僅可以反映植被、土地利用變化特征，也可以反映氣候對植被影響，通常被用于土地利用/覆蓋變化、生態環境監測以及植被動態監測等方面[17]。文中采用的NDVI數據來源于中國年度、月度植被指數(NDVI)空間分布數據集[18]，該數據集是基于連續時間序列的SPOT/VEGETATION NDVI衛星遙感數據，采用最大值合成法生成的植被指數數據集。2000年流域平均NDVI為0.68，至2016年流域平均NDVI為0.76，NDVI空間分布圖見圖3。

圖3 流域NDVI空間分布Fig.3 Spatial distribution of NDVI in the watershed

2 研究方法

2.1 Mann-Kendall突變檢驗

Mann-Kendall突變檢驗是一種非參數統計檢驗方法[19]，其優點是不需要樣本遵從一定的分布，也不受少數異常值的干擾，更適用于類型變量和順序變量，計算比較簡便。

對于有n個樣本的時間序列x，構造一個秩序列：

(1)

可見，秩序列Sk是第i時刻數值大于j時刻數值個數的累計數。在時間序列隨機獨立的假定下，定義統計量：

(2)

式中：UF1=0，E(sk)和var(sk)是sk的均值和方差，在x1,x2,…,xn相互獨立，且有相同連續分布時，它們可由下式算出：

(3)

(4)

UFi為標準正態分布，它是按時間序列x順序x1,x2,…,xn計算出的統計量序列，給定顯著性水平α，若|UFi|>Uα，則表明序列存在明顯的趨勢變化。按時間序列的逆序xn,xn-1,…,x1，再重復上述過程，同時使UBk=-UFk(k=n,n-1,…,1)，UB1=0。

2.2 NDVI空間變化趨勢分析

基于一元線性回歸的方法，構建研究區域每個柵格單元NDVI值與時間的回歸方程，采用回歸方程的斜率表示植被指數的變化趨勢[20]，計算公式如下所示：

(5)

式中：k為柵格單元回歸方程的斜率；pi為第i年NDVI均值；n為時間序列，取值為17(2000-2016年)。

k>0，表示該柵格單元NDVI值呈增加趨勢，植被覆蓋情況改善；k<0，表明該柵格單元NDVI值呈減小趨勢，植被覆蓋情況退化。

2.3 SWAT數值模擬

SWAT模型是由美國農業部農業研究中心Jeff Arnold博士研發的基于物理過程的流域尺度模型，計算包括降水、地表徑流、蒸散發、壤中流、滲透、地下水回流和河道運移損失等水量的水平衡關系，水量平衡方程如下[21]：

(6)

式中：SWt為土壤最終含水量，mm；SW0為第i天土壤初始含水量，mm；t為時間步長，d；Rday為第i天降水量，mm；Qsurf為第i天的地表徑流量，mm；Ea為第i天的蒸散發量，mm；Wseep為第i天從土壤剖面進入包氣帶的水量，mm；Qgw為第i天回歸流的水量，mm。SWAT模型產流計算流程，如圖4所示。

圖4 SWAT模型計算流程圖Fig.4 The flow chart of SWAT model

SWAT模型需要的主要輸入數據包括：DEM數據、土地利用數據、土壤類型數據、氣象數據、水文數據及土壤屬性數據，各數據來源如表1所示。

表1 SWAT 模型所需數據類型Tab.1 Data types for SWAT model

借助ArcGIS10.2，建立SWAT2012模型數據庫，基于研究流域DEM數據，劃分38個子流域。將土地利用類型重分類，對土地利用類型進行編碼賦值以匹配土地利用類型數據庫及作物數據庫，土地利用類型分布見圖5。

圖5 土地利用類型分布圖Fig.5 Distribution of land use types

土壤數據包括：土壤類型分布圖(見圖6)與土壤屬性數據，土壤基礎屬性數據主要來源世界土壤數據庫(HWSD)，根據土壤的基礎參數，如黏土(Clay)、砂土(Sand)、有機質(Organic Matter)等含量運用SPAW軟件計算土壤有效持水量、濕密度、飽和導水率。

圖6 土壤類型分布圖Fig.6 Distribution of soil types

土壤侵蝕力因子(USLE_K)采用Williams提出的方程計算：

KUSLE=fcsandfcl-siforgcfhisand

(7)

式中：fcsand、fcl-si、fhisand分別為粗糙沙土地質、黏壤土、高沙質土壤侵蝕因子；forgc為土壤有機質因子，可通過土壤各粒徑組成和有機碳含量計算得到。一般情況下，USLE_K值在0.02～0.75之間。

影響土壤產流能力的屬性是指那些影響土壤在完全濕潤并且不凍的條件下的最小下滲率屬性[22]，根據最小下滲率確定土壤水文分組，見表2。

表2 土壤水文分組Tab.2 Soil hydrologic groups

氣象資料，采用蒙陰站1984-2016年逐日溫度、風速、相對濕度、蒸發等氣象數據及19個雨量站1984-2016年逐日降水數據。

水文資料，采用蒙陰、水明崖、岸堤水庫3個水文1984-2016年的逐日徑流數據。水庫數據，采用岸堤水庫1984-2016年逐月出流數據。

對于模型適用性的驗證，采用相關系數(R2)和納什效率系數(ENS)來判斷，但目前并沒有公認的統一標準，綜合近年來的研究成果，Motovilov等[23]認為模擬徑流量時ENS>0.75結果較好，在0.36～0.75之間結果較為滿意；郝芳華等[22]、Braemort[24]認為R2>0.60、ENS>0.50時結果較為滿意。

3 結果分析

3.1 NDVI特征分析

2006年，沂蒙山被劃為國家級水土流失重點治理區，經過多年的治理，植被覆蓋條件改善，由圖7(a)可以看出，2000-2016年流域NDVI呈增大趨勢，2006年前后，流域平均NDVI由0.7增大到0.74。由Mann-Kendall法突變分析結果可以看出，見圖7(b)，除2002年外，NDVI統計量UFk均大于0，表明NDVI呈增大趨勢，其中，2007-2009年，UFk明顯大于顯著性水平線，NDVI顯著增大。在0.05置信度線內，UFk與UBk存在一個交點，該交點為NDVI的突變點。

圖7 流域NDVI特征分析Fig.7 Characteristics analysis of NDVI in the watershed

由2000-2016年年度NDVI資料，計算每個柵格NDVI的一元回歸方程的斜率，計算結果如圖8所示。k取值范圍在-0.019～0.011之間，k值大于零的面積占流域面積84%，表明流域內大部分面積NDVI都呈現增加趨勢，植被覆蓋條件改善。k值小于零的區域主要集中在蒙陰縣城周邊，主要受城鎮發展，土地利用方式變化的影響。

圖8 2000-2016年NDVI變化趨勢空間分布圖Fig.8 Spatial distribution of NDVI trend from 2000 to 2016

3.2 SWAT模擬結果分析

采用SWAT-CUP對參數進行調試、優化，1984-1989年作為模型預熱期，1990-2016年作為校準期和驗證期。根據蒙陰、水明崖和岸堤水庫3個水文站的實測流量數據驗證模擬結果，如圖9所示，3個水文站實測流量與模擬流量的過程線，變化趨勢基本一致，峰值相對應，模擬效果良好，3個水文站徑流量模擬結果與實測結果的R2及ENS均>0.5，表明調試后的模型適用于研究流域。

圖9 實測流量與模擬流量過程線Fig.9 Process line of measured runoff and simulated runoff

SWAT計算結果顯示，1990-2016年，流域平均降水量611.96 mm，實際蒸散發量383.51 mm，徑流總量225.96 mm，地表徑流量126.07 mm，基流量72.07 mm，壤中流25.02 mm，潛水補給量119.52 mm。根據子流域輸出文件(output.sub)，繪制各子流域水文要素空間分布圖見圖10。

圖10 水文要素空間分布圖Fig.10 Spatial distribution of hydrological elements

SWAT模型的預熱期為1984-1989年，數據分析采用1990-2016年。由圖11可以看出，1990-2016年流域水文要素的時間變化，從水循環三要素角度考慮，降水、蒸散發、徑流呈周期性變化，其中徑流的變化趨勢及峰值位置與降水基本一致，兩者相關度較高(R2=0.92)，且線性趨勢變化平穩，蒸散發量變幅較小，與徑流量相關度較低，如圖11(a)。徑流量由3種水源匯聚而成，包括地表水、土壤水、地下水，其分布見圖10(b)，其中占比最大的是地表水，占55%，呈線性減小趨勢，在降水量變化平穩的情況下，主要受下墊面因素影響；其次是地下水，占27%，基流量年際變化比較明顯，呈線性增加趨勢，多年平均基流系數為0.27；土壤水對徑流量的貢獻最小，占18%。潛水的補給主要來源于降水和地表水，由圖10(c)可以看出，潛水補給呈下降趨勢，地表水對其貢獻率較小，僅占5%左右，主要補給水源為降水。

圖11 水文要素時間分布圖Fig.11 Temporal distribution of hydrological elements

為分析植被覆蓋變化對地表水與地下水轉換關系的影響，以2006年為分界點，采用SWAT數值模擬結果，分析各水文要素的變化量，統計結果見表3，可以看出，2006年前后，地表徑流、基流變幅分別為-22%、30%，除受降水影響外，主要受下墊面因素的影響。2006年前后，流域的NDVI值增大，表明植被覆蓋條件改善，促進了地表水流入滲，因而地表徑流減少，地表水補給地下水量減少，基流量增大。綜上，在流域降水量減少、植被覆蓋度增大的情況下，地表水向地下水轉化量減少，地下水向地表水轉化量增多。

表3 2006年前后水文要素統計Tab.3 Statistics of hydrological elements around 2006

4 結 論

本文對沂蒙山國家水土流失治理區典型小流域的NDVI進行特征分析，采用SWAT模型開展氣候與植被覆蓋影響下的地表水與地下水的轉化規律研究，得到如下結論。

(1)2000-2016年，岸堤水庫控制流域NDVI呈增大趨勢，Mann-Kendall突變分析表明，2007-2009年，NDVI顯著增大。流域內80%以上面積NDVI都呈現增加趨勢。

(2)通過驗證，SWAT模型適用于研究流域。1990-2016年，流域平均降水量611.96 mm，實際蒸散發量383.51 mm，徑流總量225.96 mm，地表徑流量126.07 mm，基流量72.07 mm，壤中流25.02 mm，潛水補給量119.52 mm。

(3)流域徑流與降水變化趨勢基本一致，兩者相關度較高(R2=0.92)。徑流量中，地表徑流占55%，基流量占27%，壤中流占18%。潛水的補給中，地表水對其貢獻率較小，僅占5%左右，主要補給水源為降水。

(4)2006年前后，岸堤水庫流域降水量減少、植被覆蓋度增大，地表水向地下水轉化量減少，地下水向地表水轉化量增多。

□