延河流域生態系統土壤保持服務時空變化

劉 婷, 周自翔, 朱 青, 白繼洲

(西安科技大學, 西安 710054)

土壤侵蝕是全球性環境問題之一，不僅是造成土地退化、農業生產力下降的主要原因，也是阻礙流域可持續發展的重要生態環境問題，一直以來都是大家的關注熱點[1-3]。延河流域地處黃土高原，是中國水土流失嚴重的區域之一，為了遏制嚴重的水土流失，中國政府在黃土高原區域采取了一系列的水保措施，包括梯田、造林、淤地壩等[4]。隨著水保措施的實施，延河流域植被狀況逐漸好轉，土壤侵蝕狀況不斷改善，土壤保持服務所受的關注度也越來越高。基于此，諸多學者對黃土高原的土壤保持服務進行了研究，如韓永偉等[5]、楊波等[6]、王森等[7]利用多種方法評估了土壤保持服務功能及其產生的經濟價值，也有眾多學者研究土壤保持與地形因子之間的關系，例如陸傳豪等[8]基于RUSLE模型研究萬州區土壤保持服務功能空間分布特征與坡度和土地利用類型之間的關系;劉睿等[9]研究三峽庫區重慶段土壤保持功能與地形因素之間的分布關系;饒恩明等[10]以四川省為研究區域，探究土壤保持功能空間特征及地形、植被等對其的影響。以上研究在探討地形因子對土壤保持的影響時，并未充分考慮地形因子的精細度，目前以高分辨率數字高程模型(Digital Elevation Model,DEM)為基礎數據，探討地形因子對土壤保持影響的研究還比較欠缺，故可將高分辨率DEM獲取的地形因子納入到土壤保持的研究中。

本文以延河流域為研究區，利用合成孔徑雷達(Interferometric Synthetic Aperture Radar,INSAR)技術生成的高分辨率DEM，并基于SWAT(Soil and Water Assessment Tool)水文模型模擬1999—2014年產沙量，獲得流域實際土壤侵蝕量，通過修改模型源代碼，獲得潛在土壤侵蝕量，進而獲得土壤保持量，從而分析土壤保持服務的時空動態變化特征及其約束條件，為區域水土流失研究和治理提供科學依據。

1 研究區概況及數據來源

1.1 研究區概況

延河流域位于東經108°45′—110°28′，北緯36°23′—37°17′，面積為7 725 km2，平均坡度4.4%，河網長度4.7 km/km2。該區域屬于暖溫帶大陸性半干旱季風氣候，年均降水量為500 mm，年均氣溫為8.8～10.2℃。土壤類型主要以黃綿土為主，占流域面積80%以上，該土體疏松、軟綿，抗侵蝕能力差，因此，延河流域是我國水土流失的重要區域之一。研究區地勢西北高東南低，地貌類型主要包括，殘原平梁溝壑區、黃土丘陵溝壑、黃土臺地、以及黃土覆蓋的山地，其中分布最為廣泛的是丘陵溝壑區，約占研究區面積的90%。由于流域內水土流失嚴重，近年來采取了一系列的水土保持措施及植被覆蓋措施，使得流域生態狀況逐漸好轉，圖1為2003年、2014年的年NDVI空間分布特征。

1.2 數據來源

本研究所使用的數據:(1) 研究區DEM是通過選用Sentinel-1A數據為數據源，利用INSAR技術獲取的，分辨率為16.87 m，通過與水經注獲取的19.11 m的DEM進行對比驗證，發現該DEM精度較高。(2) 土地利用數據來源于中國科學院資源科學環境數據中心，分辨率為30 m，研究中主要選用2005年、2010年土地利用數據。(3) 土壤數據來源于世界土壤數據庫(Harmonized World Soil Database,HWSD)，通過土壤參數計算軟件(Soil Parameter Calculation Software,SPCS)獲取土壤參數數據。(4) 氣象數據來源于中國氣象數據共享網，包括延河流域內及周邊站點的1999—2014年的逐日降雨、溫度、相對濕度和風速數據。(5) 水文數據來源于《中華人民共和國水文年鑒》，包括甘谷驛水文站1999—2014年的逐日徑流數據和1999—2014年逐日產沙數據。(5) 梯田、水庫數據來自于Sentinel-2A數據，通過解譯影像得到2015年研究區梯田、水庫分布情況。(6) 歸一化植被指數(Normalized Difference Vegetation Index,NDVI)數據來自MOD13Q1產品，該產品具有250 m的空間分辨率，每隔16 d提供一次，通過合成處理，得到2003—2014年NDVI數據。

2 試驗材料與方法

2.1 延河流域SWAT模型的構建

2.1.1 SWAT模型土壤侵蝕子模塊以及土壤保持量的獲取 目前已有大量研究利用SWAT水文模型對黃土高原區域的產流產沙進行模擬[11-13]，因此本研究選用SWAT模型，對延河流域土壤保持量進行模擬。該模型中產沙量的計算采用改進版通用水土流失方程(Modified Universal Soil Loss Equation,MUSLE)，如式(1)所示。

Sed=11.8·(Qsurf·qpeak·areahru)0.56·Kusle·
Cusle·Pusle·LSusle·CFRG

(1)

式中:Sed為模擬產沙量;Qsurf為地表徑流(mm H2O/hm2);qpeak為徑流峰值(m3/s);areahru為HRU面積(hm2);Cusle為地表植被覆蓋因子;Pusle為水土保持措施因子;Kusle為土壤可蝕性因子;LSusle為地形因子,包括坡長因子L和坡度因子S;CFRG為直徑大于2 mm的粗碎塊因子。

(1) 土壤可蝕性因子(Kusle)可以表示土壤對侵蝕的敏感性，SWAT中計算Kusle因子的方法，具體公式如下式所示。

Kusle=fcsand·fcl-si·forgc·fhisand

(2)

式中:fcsand為沙土質地土壤侵蝕因子;fcl-si為黏壤土土壤侵蝕因子;forgc為土壤有機質因子;fhisand為高沙質土壤侵蝕因子。

(3)

(4)

(5)

(6)

式中:ms,msilt,mc分別為砂粒、粉粒、黏粒的含量(%);orgC為有機碳含量(%)。

(2) 植被覆蓋因子(Cusle)是影響土壤侵蝕的最敏感因子，主要受到覆蓋度、植被類型等的影響，SWAT模型中計算公式如下。

Cusle=exp([ln(0.8)-ln(Cusle,mn)]·

exp[-0.00115·rsdsurf]+ln[Cusle,mn])

(7)

式中:Cusle,mn為土地利用的最小植被覆蓋度因子;rsdsurf為土壤表面殘留物(kg/hm2)。

Cusle,mn因子可由已知的Cusle,aa因子用下式估計:

Cusle,mn=1.463ln[Cusle,aa]+0.1034

(8)

式中:Cusle,mn為土地利用的最小植被覆蓋度因子;Cusle,aa為土地利用年均植被覆蓋因子。

(3) 水土保持措施因子(Pusle)是指采取水保措施后，土壤流失量相對于未采取任何措施時的比例，取值范圍為0～1[14]。0表示不存在侵蝕現象的區域，1表示未采取水土保持措施的區域。SWAT模型中Pusle因子是采用Wischmeier和Smith給出的信息計算，為了使模擬結果更適合研究區，因此本文對Pusle因子進行了修改。SWAT模型中可對每一個HRU的Pusle因子進行修改，HRU中包含土地利用、坡度等信息，因此本文根據土地利用類型及坡度范圍，并考慮前人在黃土高原的研究以及黃土高原上水保措施的實施，包括梯田、水平溝以及魚鱗坑等的影響[15-16]，最終確定表1中Pusle因子的取值。

(4) 坡度坡長因子(LS)可以表示地形因子對土壤侵蝕的影響，SWAT模型中采用如下公式進行計算。

(9)

式中:Lhill為坡長(m)；m為指數項；αhill為坡度。

m的計算公式如下:

m=0.6(1-exp[-35.835·slp])

(10)

式中:slp為HRU的坡度。

slp與αhill的關系為:

slp=tanαhill

(11)

研究中涉及潛在土壤侵蝕量(Sedq)、實際土壤侵蝕量(Seda)、土壤保持量(Sedc)3種數據量的計算。潛在土壤侵蝕量是指在假定沒有地表植被覆蓋和水土保持措施下的侵蝕量，即通過將SWAT源代碼和數據庫中的Cusle，Pusle因子修改為1得到的產沙量。實際土壤侵蝕量是指流域實際發生的侵蝕量，即由SWAT模型直接模擬得到的產沙量。土壤保持量是二者的差值，計算公式如下式所示。

Sedc=Sedq-Seda=11.8(Qsurf·qpeak·areahru)0.56·
Kusle·LSusle·CFRG·(1-Cusle·Pusle)

(12)

2.1.2 數據庫的構建

(1) 土地利用數據庫。研究中主要使用2005年、2010年的土地利用數據，將延河流域土地利用數據分為耕地、林地、園地、高覆蓋草地、低覆蓋草地、建設用地、水域和未利用地8類，模型在輸入土地利用數據時，同時還需要輸入土地利用索引表，方便將流域的土地利用類型與SWAT模型自帶的植被生長模型庫和農業管理數據庫進行連接。因此，按照SWAT模型中土地利用的編碼對各種土地利用重新進行編碼，見表2，從而建立土地利用索引表。

表2 土地利用類型SWAT編碼

(2) 土壤數據庫。本文所使用的土壤數據來源于HWSD，流域的土壤類型主要包括黃綿土、沖積土、粗骨土、紅黏土、潛育始成土等，見圖2。通過在HWSD中獲取各種土壤的部分參數以及利用SPCS軟件計算部分土壤參數，建立土壤參數數據庫，并根據各類土壤在SWAT中的編碼，建立土壤參數索引表。

圖2 土壤類型

(3) 氣象數據的準備。本文選用流域內及距離流域較近的11個傳統氣象站點1999—2014年的逐日氣象資料輸入到SWAT模型中，其中降水、最高溫度、最低溫度、相對濕度、風速可直接在中國氣象數據共享網站獲取，太陽輻射數據通過天氣發生器模擬得到。

2.1.3 模型的建立 本研究使用ArcSWAT 2012版本，首先是子流域的劃分，基于DEM生成河網，通過設定閾值(形成子流域的最小給水面積)為10 000 hm2，設置流域出水口點，將延河流域劃分為41個子流域;其次是水文響應單元(HRU)的劃分，模型在劃分子流域的基礎上，通過疊加土地利用數據、土壤數據和坡度數據，獲得各個子流域中具有相同土地利用類型、土壤類型和坡度的HRU。最后將氣象數據等導入模型，對延河流域1999—2014年產沙量進行模擬。模型以1999—2002年作為預熱期，利用2003—2008年甘谷驛水文實測值進行率定，利用2009—2014年的水文實測值進行驗證。由于研究時間較長，土地利用會發生一定程度地變化，因此將預熱期和率定期視為整體，選用2005年土地利用進行模擬;在驗證期時選用2010年土地利用進行模擬。

2.1.4 模型的參數敏感性分析、率定、驗證以及不確定性分析 本文通過選用SWAT-CUP(SWAT-Calibration and Uncertainty Programs)中的SUFI-2算法對SWAT模型進行敏感性參數分析、率定、驗證以及不確定分析。在該算法中，采用判定系數(Cofficient of determination,R2)和納什效率(Nash-Sutcliffe efficiency coefficient,NSE)對模型進行適用性評價，其中R2可評價模擬值與實測值變化趨勢的一致性，NSE可反映模擬值與實測值之間的擬合程度[17]。

(13)

(14)

式中:Qmi為模擬值;Qoi為實測值;Qm為模擬均值;Qo為實測均值。

2.2 地形因子的獲取

本文選用的地形因子包括:坡度、坡向和高程，其中坡度由DEM直接獲取，將其分為0°～15°，15°～25°，25°～35°以及>35°共4個等級，獲得坡度分級圖;坡向由DEM直接獲取，將其分為平地、陽坡、半陽坡、半陰坡及陰坡5類，獲得坡向分級圖;高程由DEM直接獲取，通過設置50 m的高程間隔，分為26個等級，獲得高程分級圖。

3 結果與分析

3.1 模型率定與驗證

本文選用甘谷驛水文站1999—2014年的逐月徑流數據和泥沙數據作為實測值，由于未能獲取枯水期時泥沙數據，因此選用0值代替。在SWAT-CUP中進行敏感性參數分析后，完成對模擬結果的率定和驗證。圖3為研究區徑流和泥沙實測值與模擬值的比較結果。一般研究認為，R2>0.6,NSE>0.5時，模擬結果較可靠，校準期和驗證期模擬結果都較好時，參數的最佳校準值適用于研究區域[18]。

從圖3(A1—A2)可以看出，徑流量在率定期(2003—2008年)和驗證期(2009—2014年)模擬結果的R2和NSE均達到模型要求;從圖3(B1—B2)可以看出，泥沙量在率定期(2003—2008年)和驗證期(2009—2014年)模擬結果的R2和NSE也均達到模型要求。因此，SWAT模型可以對延河流域的水文狀況進行較好地模擬。從而可得到延河流域基于HRU的2003—2014年逐月產沙量。

圖3 延河流域甘谷驛水文站點月徑流量和月產沙量在率定期、驗證期的模擬結果

3.2 土壤保持量分布及變化

3.2.1 時間尺度變化 從圖4可以看出，在時間尺度上，潛在土壤侵蝕量、實際土壤侵蝕量和土壤保持量變化趨勢相同，均會產生不同幅度的波動，主要原因是受降水因子的影響，各個月份的降水量、降水強度均不同，其中在2013年7月份3種數據量均產生了大幅度增加且均達到多年最大值，主要原因是該月降雨總量為435 mm，遠遠超過其他月份的降水量。因此，可以發現延河流域土壤保持在時間尺度上的變化主要受降水因子的影響。

圖4 潛在土壤侵蝕量、實際土壤侵蝕量和土壤保持量隨時間變化趨勢

延河流域土壤保持量在2003—2014年，整體上呈現出波動的趨勢，但土壤保持量仍有所增加，其中2003年土壤保持總量為677.1萬t,2014年土壤保持總量為1 345.02萬t。二者的空間分布格局也產生了很大的變化。由圖5可以看出，2003年土壤保持量主要集中在流域下游，主要原因是該下游區域，地勢平坦，侵蝕輕微，因此產生的土壤保持量較高，而在2014年土壤保持量在整個流域均有分布，主要原因是該區域近年來實行了一系列的水保措施及植被覆蓋措施產生了良好的生態效應。

圖5 延河流域2003年和2014年土壤保持量

為了更加明確土壤保持服務在不同年份空間分布差異性的原因，本文嘗試分析研究區水土保持措施及植被覆蓋措施的變化及其對水土保持服務的影響。因此，本文通過解譯高分辨率影像得到研究區梯田、水庫分布情況，并用此來表示研究區水土保持措施實施的狀況。由于梯田數據的提取需要高分辨率影像的支持，因此選用10 m分辨率的Sentinel-2A數據，但該數據是最早發行時間為2015年，因此通過解譯影像得到研究區2015年梯田、水庫數據，通過觀察可知，經過多年來對坡耕地的改造，到2015年為止流域上中游梯田及水庫分布較廣泛，下游分布相對較少。本文通過處理MOD13Q1數據得到的研究區2003年、2014的年NDVI數據，并用此來表示研究區植被覆蓋狀況的變化。可以看出，經過多年來實施的退耕還林(草)工程，整個流域的植被覆蓋狀況發生了明顯地好轉，尤其是流域的上中游區域。由此可知，梯田、水庫和NDVI數據在空間上的變化可以直接證明2003年和2014年土壤保持量在空間上的差異，也進一步說明水土保持措施及植被覆蓋措施在流域治理過程中發揮了至關重要的作用。

3.2.2 空間尺度變化 從圖6可以看出，在空間尺度上，土壤保持量主要受地形因子的影響。2003—2014年延河流域年均土壤保持量為1 002.02萬t，在流域上中下游區域表現不同，其中在下游年平均土壤保持量最大，上游次之，中游最低，造成這種現象的主要原因是，流域下游的區域主要為殘原平梁溝壑區，該區域黃土侵蝕方式以雨滴擊濺侵蝕為主，侵蝕輕微，且地面平坦、土層肥厚[19]，因此有著最高的土壤保持量;流域上游主要為黃土覆蓋的山地區域，該區域25°以上的陡坡占總面積的70%以上[19]，雖易受土壤侵蝕的影響，導致該區域的潛在土壤侵蝕量大于其他區域，但由于近年來這個區域實施退耕還林力度最大，使得25°以上的區域全部退耕，進而植被不斷恢復，植被覆蓋度不斷增大，該區域實際土壤侵蝕量隨之不斷減少，所以該地區的土壤保持量相對中游較高;流域中游主要為梁峁丘陵溝壑區，該地區人類活動頻繁，墾殖嚴重，因此該區域的土壤保持量最低。

圖6 延河流域多年(2003-2014)年平均土壤保持量

為了更清楚延河流域土壤保持量受地形因子的影響，將多年平均土壤保持量進行標準化，參照蔣春麗等人的研究[20]，將標準化后的土壤保持量分為低(0～0.2)、較低(0.2～0.4)、中(0.4～0.6)、較高(0.6～0.8)、高(0.8～1.0)5個等級，見表3。延河流域主要以低保持為主，面積為5 432.94 km2，占總面積的73%。

表3 土壤保持量分級統計

3.3 不同地形因子土壤保持分布特征

3.3.1 坡度與土壤保持分布特征 經統計得到各土壤保持量在不同坡度等級所占面積比例(表4)。由表可知，各級土壤保持量隨坡度的增大均呈先增加后減小的趨勢，其中在<25°區域，低保持占比較多。主要是因為耕地在這一區域分布范圍較廣[21]，人類活動比較頻繁，易發生土壤侵蝕。>25°區域，低保持占比較小，較低、中和高保持量占比較多。這主要是因為延河流域自退耕還林以來，坡度是執行退耕的首要標準，>25°以上的坡耕地必須退耕[21]，這一措施導致>25°的區域耕地面積逐漸下降，草地、林地面積相應增加，中度和強度以下侵蝕的面積占比顯著增加，劇烈侵蝕面積顯著下降。而>35°的區域各級保持量均低于25°～35°區域，主要是因為大于35°屬于急陡坡，其耕作類型雖發生改變，但該區域地形陡峭，易受土壤侵蝕的影響。

表4 不同坡度土壤保持面積分布情況

總之，流域土壤保持隨坡度變化規律明顯，25°的坡面是土壤保持量的關鍵帶，也是今后水土保持工作應重點關注的區域。

3.3.2 坡向與土壤保持分布特征 通過將土壤保持量分級圖與坡向進行疊加分析，得到土壤保持量空間分布的坡向特征圖(圖7)。由此可知，在平地上，基本只存在低、較低和中保持量，其他各級保持量幾乎為0。經統計可得，土壤保持量在不同坡向上的規律表現為半陽坡>陰坡>陽坡>半陰坡>平地。這主要是因為對于黃土丘陵溝壑區，平地上人類活動頻繁，因此基本只存在低、較低和中保持量，且占比都較小;陽坡所獲得的太陽輻射總量高于陰坡，但陽坡土壤濕度低于陰坡，且陽坡上人類耕作現象較多[22-24]，并參照湯巧英等[24]的研究成果，可知延河流域2000年、2010年植被覆蓋度總體表現為陰坡>陽坡>平地，因此陰坡上的土壤保持量高于陽坡。又根據李勉等人[25]的研究成果可知，黃土高原不同坡向土壤侵蝕速率大小表現為半陰坡>陽坡>半陽坡，從這一規律可以看出半陽坡的土壤保持量應相對較高。由此可知，土壤保持量在不同坡向的分布與坡向的植被覆蓋、土壤濕度、侵蝕速率等有著緊密的聯系。

圖7 土壤保持量坡向分布特征

總之，流域土壤保持隨坡向變化規律明顯，主要是與坡向上的植被覆蓋、濕度、光照、侵蝕速率等因子關系密切。半陽坡為土壤保持量最大的區域，半陰坡為土壤保持量較小的區域，二者均為日后開展水保工作的重點對象。

3.3.3 高程與土壤保持分布特征 將流域高程分級圖與土壤保持量分級圖進行疊加分析，得到土壤保持量空間分布的高程特征圖(圖8)。由此可知，隨著高程的增加，土壤保持量大致呈先增大后減小的現象，土壤保持量在低、較低、中和較高保持的高程分布特征大致相同，其峰值均在1 250～1 300 m高程范圍內，土壤高保持峰值主要分布在950～1 000 m的高程范圍內，主要原因是，在1 250～1 300 m高程范圍內，主要為黃土丘陵溝壑區，是延河流域典型的地貌類型，該種地形上雖易發生土壤侵蝕[26]，但近年來在丘陵溝壑區開展了淤地壩工程，給該地區帶來了良好的生態和經濟效益，使得該地區的土壤保持量不斷增加;在950～1 000 m的高程范圍內，一般為黃土臺地，土壤肥厚，地面平坦，不易發生土壤侵蝕現象，因此存在高土壤保持量。

圖8 土壤保持量高程分布特征

總之，流域土壤保持隨高程變化分異特征明顯，950～1 000 m和1 250～1 300 m是土壤保持量最多的地帶，也應是今后水土保持的關注重點。

4 結 論

(1) 基于高分辨率的DEM，應用SWAT模型開展流域生態系統土壤保持服務的模擬，不僅提高了產沙的精度，而且拓展了SWAT模型的應用領域。

(2) 延河流域在2003—2014年期間，年均土壤保持量為1 002.02萬t，土壤保持功能有所增強，在時間上的變化趨勢主要受降水量的影響，空間上的變化趨勢主要受地形因子的影響。其中流域2003年和2014年土壤保持量的空間分布特征表明，近年來流域內實施的一系列水土保持及植被覆蓋措施使水土流失現象得到了一定程度的改善，對流域生態建設發揮了極其重要的作用。

(3) 土壤保持量受坡度、高程、坡向影響規律明顯。具體表現為，土壤保持量隨坡度、高程的增大呈先增大后減小的趨勢，土壤保持量在不同坡向上的規律表現為半陽坡>陰坡>陽坡>半陰坡>平地。在>25°的坡面、950～1 000 m和1 250～1 300 m的高程范圍、以及半陽坡上均存在較高的土壤保持量，應是今后水保工作關注的重點。

本文的研究方法和研究結果能比較真實地模擬延河流域土壤保持服務的真實狀況，但在研究過程中仍然存在一些不足之處，如在利用SWAT模型計算土壤保持量的過程中，土壤參數對產沙的模擬有很大影響，但SWAT土壤數據庫中輸入的一些土壤參數是直接根據經驗方程計算得到的，雖然具有一定的科學性，但仍然會存在一定誤差;在進行模擬結果的率定和驗證時，未能獲取枯水期時泥沙數據，因此選用0值代替;在研究地形因子對土壤保持量的影響時，發現土壤保持量同時還會受到人類活動、植被覆蓋、氣候等因子的綜合影響，因此，后續研究將進行多種因子對土壤保持量的共同影響，從而更好地分析土壤保持服務變化的原因。