譚秀翠，莊會(huì)波，季 妤

(1.山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東 泰安 271018；2.山東省水文局，濟(jì)南 250002)

地表水與地下水的水量轉(zhuǎn)化是水文循環(huán)的重要過程，地表水與地下水的相互作用及轉(zhuǎn)化關(guān)系是水文地質(zhì)等領(lǐng)域研究的熱點(diǎn)和難點(diǎn)，分析和掌握其規(guī)律特征，對(duì)進(jìn)行水資源分析評(píng)價(jià)、物質(zhì)遷移和能量變化的相關(guān)研究具有重要意義[1]。目前，對(duì)地表水與地下水轉(zhuǎn)化的研究方法有：水文測(cè)量法[2]、基流分割法[3]、水量平衡法[4]、水文地球化學(xué)法[5,6]、溫度示蹤法[7,8]、數(shù)值模型法[9,10]等。

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

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

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

1 研究區(qū)域

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

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

流域內(nèi)分布有3個(gè)水文測(cè)站(1984-2016年逐日數(shù)據(jù))，19個(gè)雨量站(1984-2016年逐日數(shù)據(jù))，1個(gè)氣象站(1959-2018年逐日數(shù)據(jù))，各站點(diǎn)空間分布見圖2。

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

歸一化植被指數(shù)NDVI(normalized difference vegetation index)對(duì)植被的生長(zhǎng)勢(shì)和生長(zhǎng)量非常敏感，因此NDVI不僅可以反映植被、土地利用變化特征，也可以反映氣候?qū)χ脖挥绊?，通常被用于土地利?覆蓋變化、生態(tài)環(huán)境監(jiān)測(cè)以及植被動(dòng)態(tài)監(jiān)測(cè)等方面[17]。文中采用的NDVI數(shù)據(jù)來源于中國年度、月度植被指數(shù)(NDVI)空間分布數(shù)據(jù)集[18]，該數(shù)據(jù)集是基于連續(xù)時(shí)間序列的SPOT/VEGETATION NDVI衛(wèi)星遙感數(shù)據(jù)，采用最大值合成法生成的植被指數(shù)數(shù)據(jù)集。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突變檢驗(yàn)

Mann-Kendall突變檢驗(yàn)是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法[19]，其優(yōu)點(diǎn)是不需要樣本遵從一定的分布，也不受少數(shù)異常值的干擾，更適用于類型變量和順序變量，計(jì)算比較簡(jiǎn)便。

對(duì)于有n個(gè)樣本的時(shí)間序列x，構(gòu)造一個(gè)秩序列：

(1)

可見，秩序列Sk是第i時(shí)刻數(shù)值大于j時(shí)刻數(shù)值個(gè)數(shù)的累計(jì)數(shù)。在時(shí)間序列隨機(jī)獨(dú)立的假定下，定義統(tǒng)計(jì)量：

(2)

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

(3)

(4)

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

2.2 NDVI空間變化趨勢(shì)分析

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

(5)

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

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

2.3 SWAT數(shù)值模擬

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

(6)

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

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

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

表1 SWAT 模型所需數(shù)據(jù)類型Tab.1 Data types for SWAT model

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

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

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

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

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

KUSLE=fcsandfcl-siforgcfhisand

(7)

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

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

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

氣象資料，采用蒙陰站1984-2016年逐日溫度、風(fēng)速、相對(duì)濕度、蒸發(fā)等氣象數(shù)據(jù)及19個(gè)雨量站1984-2016年逐日降水?dāng)?shù)據(jù)。

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

對(duì)于模型適用性的驗(yàn)證，采用相關(guān)系數(shù)(R2)和納什效率系數(shù)(ENS)來判斷，但目前并沒有公認(rèn)的統(tǒng)一標(biāo)準(zhǔn)，綜合近年來的研究成果，Motovilov等[23]認(rèn)為模擬徑流量時(shí)ENS>0.75結(jié)果較好，在0.36～0.75之間結(jié)果較為滿意；郝芳華等[22]、Braemort[24]認(rèn)為R2>0.60、ENS>0.50時(shí)結(jié)果較為滿意。

3 結(jié)果分析

3.1 NDVI特征分析

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

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

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

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

3.2 SWAT模擬結(jié)果分析

采用SWAT-CUP對(duì)參數(shù)進(jìn)行調(diào)試、優(yōu)化，1984-1989年作為模型預(yù)熱期，1990-2016年作為校準(zhǔn)期和驗(yàn)證期。根據(jù)蒙陰、水明崖和岸堤水庫3個(gè)水文站的實(shí)測(cè)流量數(shù)據(jù)驗(yàn)證模擬結(jié)果，如圖9所示，3個(gè)水文站實(shí)測(cè)流量與模擬流量的過程線，變化趨勢(shì)基本一致，峰值相對(duì)應(yīng)，模擬效果良好，3個(gè)水文站徑流量模擬結(jié)果與實(shí)測(cè)結(jié)果的R2及ENS均>0.5，表明調(diào)試后的模型適用于研究流域。

圖9 實(shí)測(cè)流量與模擬流量過程線Fig.9 Process line of measured runoff and simulated runoff

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

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

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

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

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

表3 2006年前后水文要素統(tǒng)計(jì)Tab.3 Statistics of hydrological elements around 2006

4 結(jié) 論

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

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

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

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

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

□