內(nèi)蒙古草原流域降水-徑流模擬試驗(yàn)及水文演變驅(qū)動因素辨析

張阿龍，于慶峰，2，于 嬋，房麗晶，高瑞忠

（1.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院，呼和浩特010018；2.內(nèi)蒙古農(nóng)業(yè)大學(xué)職業(yè)技術(shù)學(xué)院，內(nèi)蒙古 包頭014109；3.內(nèi)蒙古自治區(qū)水文總局，呼和浩特010020）

近年來，氣候波動和人類活動導(dǎo)致極端降水事件頻繁、徑流銳減［1-3］，流域產(chǎn)匯流過程發(fā)生變化［4-6］，進(jìn)而草原流域草地退化、土地沙化和沙塵暴等生態(tài)環(huán)境問題日益突出。國內(nèi)外學(xué)者對于流域產(chǎn)匯流過程從物理模擬、數(shù)學(xué)及數(shù)值模擬進(jìn)行研究分析，管曉祥［7］對VIC 模型、新安江模型、WBM 模型和GR4J 模型在黃河典型子流域徑流過程模擬效果進(jìn)行對比；Hamed［8］和Martinez［9］發(fā)現(xiàn)，在相同的降水歷時(shí)內(nèi)，徑流量和產(chǎn)沙量對降水強(qiáng)度存在正向響應(yīng)；姚沖［10］通過模擬降水對南方典型黏土坡面進(jìn)行研究，查明降水量、降水強(qiáng)度和坡度對坡面產(chǎn)流產(chǎn)沙量和土壤侵蝕過程的影響，并建立坡面土壤侵蝕模型；SWAT模型在域水文研究中應(yīng)用廣泛，Aronld［11］、Van Liew［12］運(yùn)用SWAT 模型分別對美國不同流域徑流和小沃希托河流域進(jìn)行水文模擬，均取得了滿意的模擬效果，Demirel［13］同時(shí)建立人工神經(jīng)網(wǎng)絡(luò)模型和SWAT模型，得出SWAT 模型顯示的均方誤差值更小，模擬效果更好。

綜上所述，人工物理模擬降水-徑流試驗(yàn)對地表產(chǎn)匯流模擬具有科學(xué)性，SWAT 模型對水文模擬效果具有優(yōu)越性，因此，本文將人工物理模擬降水-徑流試驗(yàn)和SWAT 模型共同引入徑流模擬及變化歸因分析中，研究草原流域產(chǎn)匯流特征，定量分析驅(qū)動徑流變化的氣候波動與人類活動的貢獻(xiàn)，旨在探討環(huán)境變化對草原天然徑流演變的影響，進(jìn)而為草原流域植被建設(shè)、生態(tài)環(huán)境保護(hù)和生態(tài)水文研究提供參考。

1 材料與方法

1.1 研究區(qū)概況

巴拉格爾流域（117°32'～118°30'E，43°26'～44°39'N）位于內(nèi)蒙古錫林郭勒盟（圖1），是以針茅和羊草為主的典型草原內(nèi)陸河流域，其流域面積2 703 km2，冬季寒冷、夏季干熱，平均氣溫2.39 ℃、年均降水量330.3 mm；流域主要經(jīng)濟(jì)產(chǎn)業(yè)為畜牧業(yè)和工業(yè)［2，5］。

1.2 數(shù)據(jù)來源與處理

土壤植被數(shù)據(jù)通過野外調(diào)查獲取。結(jié)合流域多年降水資料、植被特征、地形地貌和草地管理方式等設(shè)計(jì)試驗(yàn)方案，共布設(shè)取樣點(diǎn)15個(gè)（圖1），分別測得0～30 cm 土壤干容重、飽和滲透系數(shù)、飽和含水率、毛管上升含水率和田間持水率；選取50 cm×50 cm 的植被樣方，將鮮草剪去裝入保鮮袋，并在該區(qū)挖取深30 cm 土層，篩洗出根系，裝入保鮮袋，在實(shí)驗(yàn)室進(jìn)行烘干稱重測取地表和地下植被干物質(zhì)量。

圖1 研究區(qū)位置及土壤植被采樣點(diǎn)分布Fig.1 Location of balaguer river basin and distribution of sampling points

降水-徑流模型由供水、降水和徑流裝置3 部分構(gòu)成（圖2），供水裝置由水箱，水泵，供水槽組成，由雨強(qiáng)調(diào)節(jié)器控制水量大小，人工模擬雨強(qiáng)均勻的降水，降水器支架為角鋼焊接成的矩形框架，降水量大小利用雨量筒測定，經(jīng)過多次重復(fù)試驗(yàn)后在雨強(qiáng)調(diào)節(jié)器上標(biāo)定。降水-徑流模擬試驗(yàn)方案設(shè)計(jì)為在典型區(qū)隨機(jī)采集100（長）cm×50（寬）cm×30（深）cm 的低擾動土樣塊，進(jìn)行不同降水強(qiáng)度下草原產(chǎn)匯流及水土侵蝕過程的試驗(yàn)?zāi)M，分別設(shè)置不同降水強(qiáng)度、持續(xù)時(shí)間和土地坡度等多組合模擬場景，即降水強(qiáng)度為60、90 和120 mm/h，降水持續(xù)時(shí)間為1、2、3、4、5、6、7、8、9、10、11、12、13、14、15 min，坡度為10°、16°、28°等。

圖2 降水-徑流模擬裝置示意圖（單位：cm）Fig.2 Schematic diagram of simulated precipitation-runoff device

研究區(qū)氣候數(shù)據(jù)包括降水、氣溫、風(fēng)速、輻射和相對濕度等，來源于研究區(qū)內(nèi)氣象站實(shí)時(shí)監(jiān)測和中國氣象數(shù)據(jù)網(wǎng)（http：//data.cma.cn），DEM 數(shù)據(jù)來源于地理空間數(shù)據(jù)云（http：//www.gscloud.cn），土壤數(shù)據(jù)來源于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心（http：//wastdc.westgis.ac.cn）（圖3），土地利用數(shù)據(jù)來源于美國地質(zhì)勘查局（United States Geological Survey，USGS）數(shù)據(jù)共享網(wǎng)站（http：//glovis.usgs.gov/）提供的Landsat TM/OLI 影像（時(shí)間分辨率16 d，空間分辨率為30 m）解譯獲?。▓D4）。

圖3 土壤類型分布Fig.3 Distribution of the soil type

圖4 土地利用和植被覆蓋分布Fig.4 Distribution of the land use and vegetation cover

統(tǒng)計(jì)分析采用Excel@2010軟件完成，土壤-植被-根系的空間分布圖通過ArcGIS@10.3軟件制作。

1.3 研究方法

Mann-Kendall（M-K）檢驗(yàn)是一種廣泛用于時(shí)間序列趨勢和突變分析，該法計(jì)算簡便，可以識別突變開始時(shí)間及指出突變區(qū)域［14］。

克里金插值法是常用的地質(zhì)統(tǒng)計(jì)格網(wǎng)化方法，其考慮了變量空間相關(guān)特性，廣泛應(yīng)用于地質(zhì)、水文、氣象等領(lǐng)域的空間分析［15］。

利用降水-徑流模擬試驗(yàn)建立不同坡度降水強(qiáng)度與產(chǎn)流系數(shù)的數(shù)理回歸模型，由降水強(qiáng)度經(jīng)驗(yàn)頻率與降水強(qiáng)度、不同坡度空間分布面積占比的變化范圍分別劃分區(qū)間，計(jì)算各區(qū)間經(jīng)驗(yàn)頻率，公式如下：

式中：Ri為坡度在i區(qū)間的徑流系數(shù)；Qi為坡度在i區(qū)間的面積占比；Pj為降水強(qiáng)度在j區(qū)間的經(jīng)驗(yàn)頻率；rij為坡度在i區(qū)間，經(jīng)驗(yàn)頻率在j區(qū)間的徑流系數(shù)；i為坡度，0°～10°、10°～16°、16°～28°；j為降水強(qiáng)度，0～10、10～60、60～90、90～120 mm/h；R為模擬徑流系數(shù)。

SWAT 模型進(jìn)行水文演變驅(qū)動因素分析，以Nash-Suttcliffe效率系數(shù)（Ens）和決定系數(shù)（R2）進(jìn)行模型參數(shù)率定和模型驗(yàn)證評價(jià)及檢驗(yàn)?zāi)M擬合過程是否合理［17］。

以流域徑流突變前的多年平均徑流量為天然時(shí)期流域徑流的基準(zhǔn)值，則天然徑流量與徑流突變后期的多年實(shí)測平均徑流量之間的差值就是受氣候變化與人類活動的影響部分，因此，氣候變化與人類活動的影響分割公式如下［18］：

式中：QH為突變后的實(shí)測徑流；QN為突變前徑實(shí)測流均值；QA為突變前實(shí)測徑流均值；QB為突變前模擬徑流均值；QC為突變后模擬徑流均值；QCI為修正突變后模擬徑流均值；I為修正系數(shù)；ΔQC為氣候變化引起的均值變化度為；ΔQH為人類活動引起的均值變化度；WC為氣候變化引起的貢獻(xiàn)率；WH為人類活動引起的貢獻(xiàn)率。

2 結(jié)果與分析

2.1 氣象水文突變檢驗(yàn)

通過M-K 檢驗(yàn)法得到氣象水文要素的突變年份（表1），突變主要集中在20 世紀(jì)90年代到20 世紀(jì)末，1982年最低氣溫突變最早，2001年降水量突變最晚，1988年徑流深突變，突變前后變化率為-35.97%，而徑流系數(shù)突變年份為1994年，突變前后變化率為-38.57%，徑流深與徑流系數(shù)突變前后變化率較為接近。

表1 氣象水文要素突變前后變化Tab.1 Distinction of the hydrometeorological factors before and after the shift year

2.2 土壤物理及水力特性

土壤物理參數(shù)（表2）中飽和滲透系數(shù)介于0.19～14.23 m/d，但均值僅為2.87 m/d，表明土壤物理性質(zhì)空間異質(zhì)性顯著，干容重最大值、最小值與均值接近，毛管上升含水率極差最大，飽和含水率次之。

表2 土壤物理特征統(tǒng)計(jì)Tab.2 Statistical parameters of the soil physical characteristics

2.3 土壤類型及植被分布

流域下墊面特征與土壤質(zhì)地、植被分布和植被根系密切相關(guān)。研究區(qū)土壤類型、草干重、草根干重空間分布特征見圖5，土壤以砂質(zhì)壤土為主，粉砂質(zhì)壤土相間的多礫沙質(zhì)壤土和少礫砂質(zhì)壤土呈島狀分布的分布特征，草干重最大值80～90 g 主要分布在上游山谷區(qū)，分布面積較小，呈現(xiàn)出中部最少，由南向北，由東向西遞減的趨勢，草根重以200～300 g 分布區(qū)域?yàn)橹鳎畲笾到橛?00～500 g，分布于上游東北部，100～200 g呈星點(diǎn)狀分布，草根重于草種分布存在較大差異。

圖5 土壤類型、草干重、草根干重空間分布圖Fig.5 Spatial Distribution of soil type，grass dry weight and grass root dry weight

2.4 降水-徑流模擬與分析

研究區(qū)多年降水頻率分布曲線表明流域降水強(qiáng)度主要集中在0～10 mm/h，約占80%，地面坡度空間分布圖（圖6）中流域坡度主要集中在0～10°，約占流域面積的86%。

圖6 流域坡度空間分布圖Fig.6 Distribution of the basin slope

根據(jù)降水經(jīng)驗(yàn)頻率特征和不同坡度空間分布面積占比，分別將降水-徑流模擬試驗(yàn)分為4組和3組（表3）進(jìn)行。

表3 流域坡度面積與降水強(qiáng)度經(jīng)驗(yàn)頻率占比Tab.3 Proportion of different slope area and empirical frequency of precipitation Intensity

降水徑流模擬試驗(yàn)中每場降水前土壤平均含水率保持在10%～15%之間，產(chǎn)匯流特征參數(shù)見表4，隨著降水強(qiáng)度和坡度的增加，徑流總量和徑流系數(shù)均增大，地面坡度為10°時(shí)，降水強(qiáng)度為120 mm/h 的產(chǎn)流總量是降水強(qiáng)度為90 mm/h 的2.19 倍，是降水強(qiáng)度為60 mm/h 的4.69 倍；而地面坡度為16°時(shí)，降水強(qiáng)度為120 mm/h 的產(chǎn)流總量是90 mm/h 的2.07 倍，分別是60 mm/h 的4.71 倍；而坡度為28°時(shí)，降水強(qiáng)度120 mm/h 的產(chǎn)流總量是90 mm/h的1.92倍，分別是60 mm/h的4.66倍，綜上，總產(chǎn)流量不僅與降水強(qiáng)度有關(guān)，而且受坡度影響較大。當(dāng)坡度為10°時(shí)，3種降水強(qiáng)度產(chǎn)流機(jī)制主要以蓄滿產(chǎn)流為主，隨著坡度的增加，產(chǎn)流量不斷加大，產(chǎn)流機(jī)制由蓄滿產(chǎn)流為主逐漸向超滲產(chǎn)流過度，徑流系數(shù)加大。

表4 不同坡度和降水強(qiáng)度下產(chǎn)流特征參數(shù)Tab.4 Statistical parameters of runoff characteristics

根據(jù)模擬降水試驗(yàn)結(jié)果，建立不同坡度降水強(qiáng)度與產(chǎn)流系數(shù)的多項(xiàng)式回歸方程（圖7），各回歸方程R2值均大于99%，擬合效果較好。

圖7 不同坡度、降水強(qiáng)度與徑流系數(shù)關(guān)系Fig.7 Relationship between slope，rainfall intensity and runoff coefficient

不同區(qū)間降水強(qiáng)度值以中值代替，不同區(qū)間坡度值以最大值代替，計(jì)算不同坡度與降水強(qiáng)度對應(yīng)的徑流系數(shù)，結(jié)果如表5。

表5 不同坡度和降水強(qiáng)度下的徑流系數(shù)Tab.5 Runoff coefficient of the different slope and rainfall intensity

根據(jù)公式（1）、（2）計(jì)算得到天然徑流系數(shù)R值為0.424（表6），可看出該區(qū)坡度為16°～28°區(qū)間徑流系數(shù)相對較大，而坡度為0°～10°區(qū)間徑流系數(shù)較小，表明流域坡度對徑流系數(shù)影響較大。

表6 流域天然徑流系數(shù)計(jì)算Tab.6 Calculation of the basin runoff coefficient

2.5 流域徑流變化驅(qū)動因素

基于人工降水模擬試驗(yàn)建立降水強(qiáng)度-坡度-徑流系數(shù)的數(shù)學(xué)模型結(jié)果，可得突變前人類活動對巴拉格爾流域徑流系數(shù)變化的貢獻(xiàn)率為24%，氣候變化貢獻(xiàn)率為76%，突變后人類活動對巴拉格爾流域徑流系數(shù)變化的貢獻(xiàn)率為53%，而氣候變化貢獻(xiàn)率為47%（表7）。

表7 流域徑流系數(shù)變化歸因分析Tab.7 Attribution analysis of runoff coefficient

基于水文氣象要素實(shí)測數(shù)據(jù)，利用SWAT 模型研究流域徑流變化驅(qū)動因素，率定期模擬值和實(shí)測值的納什系數(shù)Ens為0.88，決定系數(shù)R2為0.87，驗(yàn)證期模擬值和實(shí)測值的Ens為0.78，R2為0.81，由SWAT 的評價(jià)標(biāo)準(zhǔn)可知，SWAT 模型在該流域模擬效果較好，驗(yàn)證期的模擬值總體略大于實(shí)測值，但徑流極大值實(shí)測值均大于模擬值（圖8）。

圖8 流域年徑流模擬值與實(shí)際值對比Fig.8 Comparison between the observed and simulated runoff

由公式（3）～（9）分析流域氣候和人類活動驅(qū)動影響徑流的貢獻(xiàn)，氣候變化貢獻(xiàn)率為79%，人類活動貢獻(xiàn)率為21%，徑流變化主要受氣候變化影響，而受到的人類活動影響較?。ū?）。

表8 流域徑流變化歸因分析Tab.8 Attribution analysis of basin runoff

3 討 論

巴拉格爾流域水文氣象要素突變主要發(fā)生在20 世紀(jì)90年代到20 世紀(jì)末，突變前后水文氣象要素差異較大，主要表現(xiàn)為氣溫升高，降水量、徑流量減少，生態(tài)環(huán)境惡化加劇，導(dǎo)致生態(tài)環(huán)境與氣象水文要素互饋陷入惡性循環(huán)，人類活動（對地下水吸奪，過度放牧，工業(yè)興起，城市擴(kuò)張等）加重了惡化趨勢，造成巴拉格爾河豐水期斷流、草場退化等生態(tài)水文現(xiàn)象凸顯。Zhang等［2］通過對巴拉格爾流域氣象水文要素突變時(shí)序研究，發(fā)現(xiàn)平均最低氣溫、平均氣溫、平均最高氣溫、徑流量、降水對環(huán)境敏感性由強(qiáng)到弱，由氣溫的突然升高導(dǎo)致了巴拉格爾河流域春季徑流量的急劇下降，因春季徑流占全年的40%以上，春季徑流的突變減少和溫度的突變增加，導(dǎo)致降水的突然減少，突變時(shí)序與白勇［19］、高瑞忠等［20］對巴拉格爾流域水文氣象突變年份一致；土壤類型和植被的空間分布對產(chǎn)匯流過程影響較大，不同的土壤類型滲透性差異較大，不同的植被對降水的截留和產(chǎn)匯流速率的影響較大，由于研究區(qū)土壤類型的質(zhì)地和土壤水力參數(shù)相近，植被分布主要以草地為主，其他地物類型星點(diǎn)分布，雖然對徑流系數(shù)有一定的影響，但相對流域整體而言可忽略不計(jì)。魯克新等［21］對黃土區(qū)降水產(chǎn)流產(chǎn)沙過程研究認(rèn)為裸坡坡面下，認(rèn)為雨強(qiáng)是影響徑流的主要因素，坡度是影響產(chǎn)沙的主要因素，忽略了模擬降水過程中蒸散發(fā)作用的影響和流域下墊面的空間異質(zhì)性等因素的干擾，徑流系數(shù)可以間接反映流域內(nèi)自然地理要素對降水-徑流關(guān)系的影響，由于人工模擬降水-徑流關(guān)系可排除人類活動造成的干擾，且各要素均可人為控制，根據(jù)研究區(qū)地形地貌多年水文氣象資料計(jì)算出各要素的分布頻率，再通過人工降水模擬得理想條件較接近天然狀態(tài)的徑流系數(shù)，雖然該方法對天然徑流系數(shù)的刻畫存在一定的誤差，但從模擬結(jié)果看，通過人工降水模擬試驗(yàn)建立降水強(qiáng)度-坡度-徑流系數(shù)的數(shù)學(xué)模型得出天然徑流系數(shù)值為0.424，模擬值與突變前的徑流系數(shù)接近，較突變前多年平均實(shí)測徑流系數(shù)大，因巴拉格爾河流域產(chǎn)流機(jī)制主要為超滲產(chǎn)流，多年平均徑流系數(shù)的改變主要受下墊面條件控制，下墊面的變化主要受人類活動的影響；試驗(yàn)結(jié)果表明突變前人類活動對巴拉格爾流域徑流系數(shù)變化的貢獻(xiàn)率為20%左右，突變后為50%左右。利用分布式水文模型SWAT對徑流變化歸因分析發(fā)現(xiàn)巴拉格爾流域上游氣候變化貢獻(xiàn)率為79%，人類活動貢獻(xiàn)率為21%，徑流變化主要受氣候變化影響，而受到的人類活動影響較小，與人工降水模擬試驗(yàn)建立降水強(qiáng)度-坡度-徑流系數(shù)的數(shù)學(xué)模型得出的貢獻(xiàn)率存在較大差異，由于SWAT 模型假設(shè)突變前的多年平均徑流量為天然時(shí)期流域徑流的基準(zhǔn)值，縮小了人類活動貢獻(xiàn)率，進(jìn)而導(dǎo)致突變后人類活動對徑流的影響的貢獻(xiàn)率偏小，但SWAT 模型模擬結(jié)果與蘇輝東等［22］對黃河、黑河、長江、雅魯藏布江流域四個(gè)寒區(qū)流域得出氣候變化對徑流的影響貢獻(xiàn)率達(dá)到78%以上是徑流演變的主導(dǎo)作用的結(jié)論相近，而與高瑞忠［20］對巴拉格爾流域和王威娜［23］、于嬋［24］等對錫林河流域的研究相反，可能是使用的水文數(shù)據(jù)序列和徑流歸因分析的方法不同造成。

變化環(huán)境下徑流量的減少逐年加劇，由于突變后經(jīng)濟(jì)飛速發(fā)展，人類活動范圍的不斷擴(kuò)大（過度放牧對草場的破壞、地下水開采對徑流量的吸奪、礦產(chǎn)資源開發(fā)對地下水的輸干和對下墊面條件的破壞等）和全球持續(xù)增溫等因素的影響，導(dǎo)致草原生態(tài)系統(tǒng)退化，地下水位遞減，草場面積萎縮等環(huán)境問題突出。氣溫突變增高導(dǎo)致徑流減少，徑流突變減少和氣溫的突變增高共同導(dǎo)致了降水的突變減少，促使巴拉格爾流域水文氣象極端事件增多，干旱加劇，水文氣象要素與生態(tài)環(huán)境互饋機(jī)制改變，但氣候變化也受人類活動的干擾。張樹磊等［25］認(rèn)為在氣候較為干燥的地區(qū)，徑流對氣候和下墊面變化都更為敏感，且區(qū)域差異性明顯；劉昌明等［26］和張曉明等［27］認(rèn)為森林植被增加不僅能顯著減少豐水、平水及枯水期年的徑流，也改變了地下水與地表水的分配。

學(xué)界對高原內(nèi)陸河草原流域降水強(qiáng)度-坡度-徑流系數(shù)之間關(guān)系的深入研究相對較少，需進(jìn)一步探討。對不同降水模式下流域產(chǎn)匯流試驗(yàn)研究中由于缺乏不同植被蓋度的對比研究，降水-徑流模擬設(shè)備還需進(jìn)一步改進(jìn)完善，需要加入測量土壤實(shí)時(shí)含水率、徑流無間斷測定等設(shè)備，考慮壤中流與地下徑流的水文過程，以及各徑流中的化學(xué)成分變化等，徑流系數(shù)模型建立未考慮蒸散發(fā)，壤中流和地下徑流等，由于降水強(qiáng)度取值較少，對降水強(qiáng)度與徑流系數(shù)的擬合曲線可能存在一定的影響，對徑流系數(shù)數(shù)學(xué)模型還需進(jìn)一步研究，對徑流歸因分析方法間的差異仍需要進(jìn)一步探究。

4 結(jié) 論

（1）流域水文氣象要素突變主要集中在20 世紀(jì)90年代到20世紀(jì)末，最低氣溫突變最早，降水量突變最晚。

（2）流域產(chǎn)流量不僅與降水強(qiáng)度有關(guān)，而且受坡度影響較大，當(dāng)坡度為10°時(shí)，3 種降水強(qiáng)度產(chǎn)流機(jī)制主要以蓄滿產(chǎn)流為主，隨著坡度的增加，產(chǎn)流量不斷加大，產(chǎn)流機(jī)制由蓄滿產(chǎn)流為主逐漸向超滲產(chǎn)流過度，徑流系數(shù)加大。

（3）根據(jù)人工降水模擬試驗(yàn)發(fā)現(xiàn)天然徑流系數(shù)R值為0.424，突變前人類活動對巴拉格爾流域徑流系數(shù)變化的貢獻(xiàn)率為20%左右，突變后為50%左右。

（4）基于SWAT 模型研究流域徑流變化驅(qū)動因素表明巴拉格爾流域突變后徑流變化的人類活動貢獻(xiàn)率為21%。