黑 哲, 王 飛, 韓劍橋, 史尚渝, 戈文艷

(1.中國科學(xué)院 水利部 水土保持研究所, 陜西 楊凌 712100; 2.西北農(nóng)林科技大學(xué)水土保持研究所, 陜西 楊凌 712100; 3.中國科學(xué)院大學(xué), 北京 100049)

植被是陸地生態(tài)系統(tǒng)中的重要組成部分，在全球物質(zhì)循環(huán)和能量流動中有重要作用[1]。研究地表植被蓋度及其變化，對于揭示植被生長狀況，空間分布規(guī)律，評價區(qū)域生態(tài)環(huán)境質(zhì)量具有重要意義[2-3]。同時植被在防治土壤侵蝕過程中有重要作用[4]。一方面，植物冠層能夠攔截降雨，植物枯枝落葉層能夠增加入滲來減少徑流[5]。另一方面，植被能夠阻擋泥沙運移，減少水土流失[6-8]。在植被與土壤侵蝕的關(guān)系研究中，植被覆蓋是影響土壤侵蝕的重要指標。傳統(tǒng)的植被覆蓋測量雖然在不斷改進，但仍然有野外操作不便、成本較高、難以快速大范圍動態(tài)監(jiān)測等缺點，而遙感技術(shù)可以有效解決這一問題[9-11]。歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)能夠相當精確地反映植被蓋度、生長狀況、生物量以及光合作用強度等，是目前最常用的植被指數(shù)[12-14]。為了進一步認識植被與土壤侵蝕的關(guān)系，有必要研究植被的時空變化對土壤侵蝕風(fēng)險性的影響。

土壤侵蝕風(fēng)險性研究主要依據(jù)土壤侵蝕預(yù)報模型，國際上先后開發(fā)了通用土壤流失方程(Universal Soil Loss Equation,USLE)、修正的通用土壤流失方程(Revised Universal Soil Loss Equation,RUSLE)、水蝕預(yù)報模型(Water Erosion Prediction Project,WEPP)、歐洲水蝕預(yù)報模型(European Soil Erosion Model,EUROSEM)等。此外一些學(xué)者基于USLE建立了若干侵蝕預(yù)報模型[15]，但是以上土壤侵蝕預(yù)報模型都需要獲取降雨侵蝕力、土壤可蝕性、坡度坡長、作物覆蓋和水土保持措施因子基礎(chǔ)數(shù)據(jù)，且需要多年長期監(jiān)測數(shù)據(jù)，在基礎(chǔ)設(shè)施薄弱的地區(qū)并不適用。有學(xué)者基于植被蓋度、地形和土地利用進行水土流失風(fēng)險評價，使得大范圍長時間尺度的土壤侵蝕研究成為可能[16-18]。遙感技術(shù)能夠?qū)Φ匦?、植被覆蓋、土地利用方式和地表裸露程度等土壤侵蝕的重要標志進行規(guī)則重復(fù)的觀測，是進行動態(tài)環(huán)境監(jiān)測的重要技術(shù)手段，目前被廣泛應(yīng)用于土壤侵蝕監(jiān)測[19-20]。彭文甫等[21]利用Landsat5/8影像研究了岷江汶川—都江堰段植被覆蓋動態(tài)變化；彭雙云等[22]利用RUSLE模型對滇池流域1999—2014年土壤侵蝕時空演變進行了分析；陳學(xué)兄[23]利用坡度和植被蓋度的指標對中國水土流失進行評價。

白于山區(qū)是陜西省水土流失最為嚴重的地區(qū)之一，該地區(qū)氣候干旱，降雨集中，地形破碎，加之退耕實施不徹底，經(jīng)濟發(fā)展落后等因素的影響，導(dǎo)致區(qū)域水土流失現(xiàn)象并未明顯緩解，土壤侵蝕風(fēng)險一直較高。但目前對其土壤侵蝕的評估還較少涉及，因此快速評估白于山區(qū)植被狀況和土壤侵蝕風(fēng)險，對水土保持措施決策提供依據(jù)具有非常重要的意義。為此，本文利用MODIS-NDVI和Landsat數(shù)據(jù)，結(jié)合GIS手段，分析白于山區(qū)退耕還林(草)以來植被覆蓋和土壤侵蝕強度變化特征，為該地區(qū)生態(tài)環(huán)境建設(shè)提供科技支持。

1 材料與方法

1.1 研究區(qū)概況

白于山區(qū)(107°16′17″—109°1′55″E，36°49′11″—37°32′55″N)位于陜北榆林市與延安市的交界處，北鄰毛烏素沙地，南為黃土丘陵溝壑區(qū)。全區(qū)包括定邊縣、靖邊縣的南部，吳起縣的大部，總面積約8 706.48 km2，最高海拔1 913 m，最低海拔1 178 m。白于山梁長溝深，溝坡陡峻，由主梁分出的次一級長梁分別向東、南、北方向延伸，呈波狀階梯式降低，為無定河、洛河、延河、清澗河等陜北主要河流的發(fā)源地。該區(qū)屬于溫帶半干旱大陸性季風(fēng)氣候，年均氣溫為7.8～10℃，年降雨量320～490 mm，主要集中在5—9月。土壤類型多為細沙黃綿土，質(zhì)地均一，土質(zhì)疏松，抗侵蝕能力差，加之坡面流水侵蝕、重力侵蝕活躍，退耕還林實施面積小，坡耕地較多，水土流失嚴重。

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

本研究使用的遙感數(shù)據(jù)來源于美國國家航天航空局(National Aeronautics and Space Administration,NASA)EOS/MODIS衛(wèi)星的MOD13Q1 v006數(shù)據(jù)產(chǎn)品數(shù)據(jù)集(https:∥lpdaac.usgs.gov/)，空間分辨率為250 m，時間分辨率為16 d，2000—2018年內(nèi)共有434期影像。該產(chǎn)品屬于MODIS三級數(shù)據(jù)產(chǎn)品，進行過大氣校正和幾何糾正。研究中使用MRT(MODIS Reprojection Tool)工具對數(shù)據(jù)進行波段提取、投影轉(zhuǎn)換、格式變換等處理并用最大值合成法(maximum Value Compositess,MVC)獲得月柵格圖像和年柵格圖像。為了更精確地研究植被蓋度的地形分布特征和土壤侵蝕特征，本研究利用30 m分辨率的Landsat數(shù)據(jù)與相同分辨率的DEM數(shù)據(jù)和土地利用數(shù)據(jù)進行疊加。Landsat數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http:∥www.gscloud.cn/)，空間分辨率為30 m，選取了2010年、2017年和2018年7月份無云覆蓋的影像。利用ENVI 5.1進行輻射定標和大氣校正并用像元二分模型計算植被蓋度。數(shù)字高程模型(DEM)是空間分辨率為30 m的ASTER GDEM。土地利用數(shù)據(jù)來源于清華大學(xué)Finer Resolution Observation and Monitoring of Global Land Cover數(shù)據(jù)庫(http:∥data.ess.tsinghua.edu.cn/)，空間分辨率為30 m。

1.3 研究方法

1.3.1 最大值(MVC)合成法 最大值合成法是將每個月的最大值假設(shè)為天氣晴朗時獲取的影像，通過最大值合成方法可以較好地去除云和大氣的干擾?；颈磉_式如下：

NDVIi=(NDVIj)max

(1)

式中：NDVIi為第i月NDVI影像(i=1，2，…，12)，NDVIj為當月兩景NDVI影像(j=1，2)。

1.3.2 NDVI趨勢分析 為了對NDVI變化逐像元分析，采用一元線性回歸分析法可計算每個柵格的NDVI值在研究時段中的變化趨勢，并估計變化幅度。本文使用該方法模擬2000—2018年19 a平均的變化趨勢，其計算公式如下：

(2)

表1 植被覆蓋變化趨勢分類標準

1.3.3 土壤侵蝕分級 研究中使用ArcGIS將研究區(qū)的土地利用類型、坡度和植被蓋度數(shù)據(jù)疊加分析，并參照《土壤侵蝕分類分級標準》(SL190—2007)，將土壤侵蝕強度分為5個等級[24](表2)。

1.3.4 像元二分法 為了評估研究區(qū)土壤侵蝕強度，需要利用像元二分模型將植被指數(shù)轉(zhuǎn)換為植被蓋度(Fractional vegetation coverage，F(xiàn)VC)，再對植被指數(shù)進行等級劃分。植被覆蓋度計算公式如下：

(3)

式中：FVC為植被蓋度值；NDVI為各像元的NDVI值；NDVIsoil，NDVIveg分別為純土壤覆蓋像元和純植被覆蓋像元的NDVI值。根據(jù)李苗苗等[25]提出的估算NDVIsoil和NDVIveg的方法，將研究區(qū)內(nèi)NDVI累計頻率為0.5%，99.5%時的值作為全土壤覆蓋和全植被覆蓋下的NDVI值，然后利用上式計算每像元植被覆蓋度。

表2 土壤侵蝕分級指標

2 結(jié)果與分析

2.1 植被蓋度空間變化特征

2.1.1 NDVI隨高程變化特征 白于山區(qū)高程分布如圖1所示，大部分區(qū)域高程范圍在1 300～1 700 m，高程最高的區(qū)域位于定邊縣南部白灣子鎮(zhèn)，從中心到四周高程環(huán)形遞減。隨著山區(qū)內(nèi)高程的增加(圖2)，植被蓋度均呈現(xiàn)出“增加—減少”的趨勢，不同年份均在1 250～1 400 m范圍內(nèi)達到峰值。

圖1 白于山區(qū)高程分布

圖2 2000年、2010年、2018年白于山區(qū)NDVI隨高程的變化

2.1.2 NDVI隨坡度變化特征 以2°為間距來估算研究區(qū)NDVI與坡度關(guān)系表明(圖3)，2000年、2010年、2018年NDVI隨坡度的變化趨勢均為“增加—減少”，2000年植被蓋度整體偏低，在坡度≤50°范圍內(nèi)NDVI呈現(xiàn)緩慢增長趨勢，>50°陡坡處達到最高值0.34。2010年NDVI均值為2000年的1.38倍，2018年NDVI相比2010年增加了45%，其中在坡度0°～25°范圍內(nèi)快速增長，在坡度為56°時達到最高值0.65，在56°以上，各年份植被蓋度均顯著減少。

圖3 2000年、2010年、2018年白于山區(qū)NDVI隨坡度的變化

2.1.3 NDVI隨坡向變化特征 將研究區(qū)坡向按照平地(-1°)、陰坡(0°—45°,315°—360°)、半陽坡(45°—135°)、陽坡(135°—225°)、半陰坡(225°—315°)分為5個坡向，按照年份分別統(tǒng)計研究區(qū)NDVI隨坡向的變化特征研究，結(jié)果表明植被覆蓋在坡向上的差異較小，2000年、2010年總體植被蓋度隨坡向變化趨勢相同，從小到大依次為：陽坡<半陽坡<平地<半陰坡<陰坡，2018年平均植被蓋度隨坡向變化由小到大依次為：平地<陽坡<陰坡<半陽坡<半陰坡，表明坡面植被得以較快恢復(fù)。

2.1.4 NDVI空間分布特征 白于山區(qū)2000—2018年平均NDVI空間分異比較明顯(圖4)，其中植被覆蓋最好的地區(qū)位于吳起縣，NDVI最高值可達0.72，白于山區(qū)的NDVI多年平均值主要分布在0.2～0.3，占整個區(qū)域73.63%，整個研究區(qū)NDVI平均值僅有0.25，整體上植被覆蓋偏低(表3 )。

圖4 2000-2018年白于山區(qū)平均NDVI空間分布

2.2 植被蓋度時間變化特征

2.2.1 NDVI年際變化特征 通過最大值合成法(MVC)處理得到2000—2018年研究區(qū)各像元年最大NDVI值并在整個研究區(qū)進行平均，以年平均值代表研究區(qū)植被覆蓋最佳狀況。圖5是2000—2018年研究區(qū)植被覆蓋年變化趨勢，植被蓋度在2000—2018年呈極顯著上升趨勢(p<0.001)，其增速為0.065/10 a?？傮w來說研究區(qū)植被蓋度的年際變化大致可以分為3個階段，第1階段為2000—2003年，是大規(guī)模生態(tài)建設(shè)初期，植被覆蓋隨年份呈穩(wěn)定增加趨勢，其增速達到0.36/10 a(p<0.001)。第2階段為2004—2015年，該階段植被覆蓋處于相對穩(wěn)定狀態(tài)，NDVI最低值(0.41)出現(xiàn)在2015年。第3階段為2016—2018年，植被呈快速增長趨勢，增速達到0.59/10 a(p<0.001)。

表3 2000-2018年白于山區(qū)平均NDVI空間分布

注：研究區(qū)平均NDVI值為0.248。

圖5 2000-2018年白于山區(qū)NDVI年際變化趨勢

為了進一步分析2000—2018年白于山區(qū)植被覆蓋的年際空間變化趨勢，經(jīng)過Slope分析并通過Mann-Kendall檢驗計算出置信度水平。根據(jù)計算結(jié)果將趨勢變化劃分為嚴重退化到明顯改善，將M-K檢驗計算結(jié)果劃分為無顯著變化、弱顯著變化(p<0.1)、顯著變化(p<0.05)和極顯著變化(p<0.01),結(jié)果見圖6—7和表4，空間格局的主要變化特征有：(1) 呈退化趨勢的區(qū)域比重占3.18%，分布在定邊縣白于山南部丘陵溝壑區(qū)和靖邊縣中部的零星區(qū)域；其中定邊縣東部零星區(qū)域呈現(xiàn)極顯著(p<0.01)退化趨勢，主要是因為工程建設(shè)導(dǎo)致植被退化。(2) 基本不變區(qū)域比重為32.62%，連片分布在定邊縣中部地帶即白于山南緣和靖邊縣中南部區(qū)域；(3) 輕微改善區(qū)域占研究區(qū)面積的45.66%，連片分布于吳起縣東部、靖邊縣北部和南部區(qū)域、定邊縣北部區(qū)域。而且這些區(qū)域的改善趨勢都通過了0.001的信度檢驗，為極顯著改善；(4) 中度改善的區(qū)域比重為16.38%，主要是吳起縣的東部區(qū)域、定邊縣白于山北部和靖邊縣的部分區(qū)域，改善趨勢為極顯著(p<0.01)；(5) 明顯改善的區(qū)域較小，比重為2.16%，零星分布在吳起縣的東部和靖邊縣的南部。

圖6 白于山區(qū)多年平均NDVI時空演變

圖7 白于山區(qū)NDVI變化趨勢置信度水平

表4 2000-2015年白于山區(qū)植被覆蓋演變趨勢

2.2.2 NDVI年內(nèi)變化特征 為了研究白于山區(qū)季節(jié)尺度上NDVI變化特征，采用氣象學(xué)標準進行季節(jié)劃分：春季3—5月，夏季6—8月，秋季9—11月，冬季12月—次年2月。圖8為白于山區(qū)4個季節(jié)平均NDVI的年際變化，可以看到冬季的NDVI值最低，其次為春季、秋季，夏季最高，這反映了植被與物候的關(guān)系。統(tǒng)計得出，4個季節(jié)均呈現(xiàn)顯著(p<0.05)上升趨勢，但是增速有所不同，植被增速由快向慢依次為夏季(0.069/10 a)>秋季(0.044/10 a)>春季(0.039/10 a)>冬季(0.017/10 a)。

圖8 白于山區(qū)季度NDVI變化趨勢

2.2.3 NDVI月變化特征 由圖9可知，多年NDVI的月變化均隨年份呈現(xiàn)不同程度的增長趨勢，其中1月、2月呈顯著(p<0.05)上升趨勢，其余10個月均呈極顯著(p<0.001)上升趨勢。增速最快的為7月(0.082/10 a)，最慢為2月(0.015/10 a)。研究區(qū)植被覆蓋隨月份變化顯著，NDVI為0.15～0.46，其中2月最低，8月最高，平均值為0.25，其中6—9月都在平均值以上,植被蓋度在年內(nèi)的變化曲線整體呈單峰型波動。

2.3 土壤侵蝕變化特征

2.3.1 土壤侵蝕等級結(jié)構(gòu)變化特征 根據(jù)水利部標準，利用植被蓋度、坡度和耕地數(shù)據(jù)對研究區(qū)2010年和2017年土壤侵蝕狀況進行評價，并統(tǒng)計出土壤侵蝕面積(圖10，表5)。2010年和2017年土壤侵蝕等級分布面積最大的等級均為中度侵蝕，2010年土壤侵蝕面積為中度侵蝕>輕度侵蝕>強烈侵蝕>極強烈侵蝕>劇烈侵蝕，2017年土壤侵蝕面積為中度侵蝕>強烈侵蝕>輕度侵蝕>極強烈侵蝕>劇烈侵蝕。2010—2017年期間，輕度侵蝕面積顯著減少，由2010年的3 437.87 km2減少至2017年的1 728.22 km2，輕度侵蝕面積占比減少了19.64%。中度侵蝕的面積變化較小，減少了5.23%。強烈侵蝕的面積變化較大，增加面積為15.24%。極強烈侵蝕和劇烈侵蝕面積均為小幅度增加，分別增加了7.64%，1.98%。

對比分析2010年和2018年研究區(qū)土壤侵蝕等級結(jié)構(gòu)變化，計算得到兩年土壤侵蝕強度轉(zhuǎn)移矩陣(表6)。2010—2017年期間，輕度侵蝕面積減少19.64%，其中45.64%轉(zhuǎn)化為中度侵蝕，5.68%轉(zhuǎn)化為強烈侵蝕。中度侵蝕面積減少5.23%，其中1.43%的面積好轉(zhuǎn)為輕度侵蝕，44.14%的面積轉(zhuǎn)化為強烈侵蝕，11.80%的面積轉(zhuǎn)化為劇烈侵蝕。強烈侵蝕的面積增加了15.24%，其中31.43%來源于極強烈侵蝕，7.92%來源于劇烈侵蝕，說明土壤侵蝕狀況略有好轉(zhuǎn)。極強烈侵蝕的面積增加7.64%，其中有20%來自劇烈侵蝕，11.04%來源于中度侵蝕。劇烈侵蝕的面積增加較少，僅有1.98%，其中4.69%來自于強烈侵蝕，15.56%來自于極強烈侵蝕，說明極強烈侵蝕的土壤在2010—2017年改善不明顯。

圖9 白于山區(qū)月均NDVI變化趨勢

圖10 2010年和2017年土壤侵蝕強度

表5 2010年、2017年白于山區(qū)土壤侵蝕面積

2.3.2 土壤侵蝕空間變化特征 為了探究2010—2017年土壤侵蝕等級變化的空間異質(zhì)性，將白于山區(qū)2010年和2017年土壤侵蝕強度圖疊加、相減，并分區(qū)得到土壤侵蝕等級變化分布圖(圖11)。其中4～-4代表侵蝕狀況由顯著改善到顯著退化。對比2010年和2017年土壤侵蝕等級，其中，51.30%的面積土壤侵蝕有所改善，42.21%的面積土壤侵蝕無明顯變化，6.49%的面積土壤侵蝕狀況加劇，即研究區(qū)內(nèi)大部分面積的土壤侵蝕有所改善。其中定邊縣土壤侵蝕加劇的面積占比為8.06%，稍大于吳起縣的6.08%和靖邊縣的5.97%，整體退化趨勢不顯著。定邊縣和靖邊縣內(nèi)侵蝕無變化的區(qū)域最大，面積占比分別為48.33%，41.47%，說明定邊縣和靖邊縣從2010—2017年土壤侵蝕狀況未有明顯好轉(zhuǎn)，吳起縣內(nèi)土壤侵蝕輕微改善的面積最大，為41.08%。土壤侵蝕中度改善和明顯改善的面積為吳起縣>靖邊縣>定邊縣，即吳起縣土壤侵蝕狀況有較好改善。

3 討 論

植被覆蓋和土地利用方式是影響土壤侵蝕的重要指標，目前已經(jīng)受到廣大學(xué)者的關(guān)注[26-27]。本研究基于植被蓋度、坡度和土地利用3個因子評價區(qū)域土壤侵蝕風(fēng)險性，與傳統(tǒng)的土壤流失方程相比，本文重點在于快速便捷地評價區(qū)域內(nèi)土壤侵蝕風(fēng)險性[28-29]。同時利用長時間序列的遙感數(shù)據(jù)對研究區(qū)平均NDVI時空變化進行了分析，關(guān)于植被與土壤侵蝕風(fēng)險性的研究表明，隨著植被蓋度的增加，土壤侵蝕風(fēng)險逐漸降低。這是由于植被具有較高的水土保持功能，其冠層能夠防止雨滴直接擊打地面，根系能夠吸收降水，共同阻止土壤侵蝕的發(fā)生。如張雪峰等[30]對內(nèi)蒙古錫林河流域草地生態(tài)系統(tǒng)的研究表明，提高植被覆蓋可以顯著提升草地生態(tài)系統(tǒng)的土壤保持能力。孫文義等[31]對黃土高原不同生態(tài)系統(tǒng)水土保持服務(wù)功能評價研究中表明植被抗土壤侵蝕能力隨植被覆蓋度的增加而增加，認為植被蓋度大于70%時，地表土壤侵蝕很微弱。本研究與上述結(jié)論基本一致。除植被蓋度之外，植被類型對土壤侵蝕的影響也有所不同[32]。例如實施“退耕還林(草)”等生態(tài)修復(fù)工程時植樹造林的植被與農(nóng)耕地分布有較大差異[33]。之后可以進一步分析白于山區(qū)不同植被類型對土壤侵蝕風(fēng)險的影響，降雨侵蝕力與不同類型植被NDVI的一致性。

表6 2010年和2017年白于山區(qū)土壤侵蝕強度轉(zhuǎn)移矩陣 km2

圖11 白于山區(qū)土壤侵蝕等級變化分布

4 結(jié) 論

(1) 研究時段內(nèi)，年均NDVI呈現(xiàn)顯著上升趨勢，增長速率為0.065/10 a。

(2) NDVI空間分布上，吳起縣NDVI較高，而定邊縣西部和靖邊縣中部NDVI較低。隨著海拔和坡度的增加,NDVI變化趨勢均為先增大后減小，其中在海拔1 250 ～1 400 m達到峰值,在坡度50°～56°達到峰值；NDVI的坡向差異性較小，基本趨勢為陰坡植被覆蓋大于陽坡。

(3) 研究區(qū)植被退化面積為276.88 km2，主要集中在白于山南部丘陵溝壑區(qū)，植被保持穩(wěn)定的面積為2 839.96 km2，主要集中在定邊縣中部，植被改善的面積為5 589.67 km2，主要集中于白于山北部和靖邊縣北部區(qū)域。

(4) 2010年和2017年二期土壤侵蝕狀況表明，輕度侵蝕和中度侵蝕的面積減少，強烈侵蝕、極強烈侵蝕和劇烈侵蝕的面積增加。

(5) 2010年和2017年研究區(qū)土壤侵蝕加重的面積占比僅為6.49%，主要集中在定邊縣的北部和靖邊縣的中西部區(qū)域；土壤侵蝕改善的面積占比為51.30%，主要集中在吳起縣。