退耕還林（草）工程對(duì)黃土高原地表產(chǎn)沙的影響

管亞兵 ，趙長(zhǎng)森 ，楊勝天 ，※，白 娟 ，彭睿文

（1. 北京師范大學(xué)水科學(xué)研究院，北京100875；2. 遙感科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京師范大學(xué)地理科學(xué)學(xué)部，北京100875；3. 貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院，貴陽(yáng)550001）

0 引言

作為陸地生態(tài)系統(tǒng)的主要組分，植被是影響土壤侵蝕與水土流失的重要因子。植被具有顯著的減水減沙效益，主要體現(xiàn)在：（1）植被冠層通過(guò)截留降雨減小降雨侵蝕力，同時(shí)增加地表糙度降低徑流對(duì)地表沖刷能力；（2）植被根系可以改善土壤結(jié)構(gòu)，增加土壤抗沖性能和入滲能力［1-2］。大量研究已經(jīng)證明植被對(duì)流域水沙變化有深刻影響，表征植被長(zhǎng)勢(shì)的植被覆蓋度指標(biāo)與產(chǎn)沙系數(shù)之間存在較好的相關(guān)性［3-7］。基于以上認(rèn)識(shí)，植被生長(zhǎng)狀況現(xiàn)已被廣泛用來(lái)作為土壤侵蝕估算和風(fēng)險(xiǎn)評(píng)價(jià)的重要指標(biāo)［7-10］。

黃土高原是世界上水土流失最嚴(yán)重的區(qū)域之一，約60%以上的土地已經(jīng)遭受水土流失［11］，土壤侵蝕模數(shù)可達(dá)15 000 t/km2/yr。 土壤侵蝕可能導(dǎo)致嚴(yán)重的土地問(wèn)題，如耕地蠶食、土壤有機(jī)質(zhì)和養(yǎng)分流失［12］、土壤物理性質(zhì)退化［13］。 進(jìn)入河道的泥沙又將在下游不斷淤積，抬高河床形成地上河，危害水壩和水庫(kù)安全。因此，黃土高原水土流失不僅限制了當(dāng)?shù)亟?jīng)濟(jì)社會(huì)的發(fā)展，而且威脅到黃河下游人民群眾的生命財(cái)產(chǎn)安全［14］。

自1999 年退耕還林草政策實(shí)施以來(lái)，黃土高原植被覆蓋率從1999 年的31.6%增加到2013 年的59.6%，增幅達(dá)88%，植被覆蓋呈明顯增加趨勢(shì)［15］。同期黃河水沙銳減現(xiàn)象日益明顯。黃河天然來(lái)沙量 16 億t，天然徑流量560 億m3，但在2003—2013 年降雨略偏豐情況下，潼關(guān)實(shí)測(cè)年均來(lái)沙量只有2.6 億t（其中河潼間偏少85%），花園口站天然徑流量只有488 億m3（其中河潼間偏少40%）［6］。在此背景下，很多學(xué)者研究了黃河流域林草植被變化特征［16］、林草植被變化對(duì)產(chǎn)流產(chǎn)沙的影響機(jī)制［6］、林草植被減沙方程構(gòu)建［4-5］、林草植被減蝕臨界效應(yīng)［17］。然而，以往黃河流域植被變化與地表產(chǎn)沙關(guān)系研究多在試驗(yàn)小區(qū)或流域尺度上開(kāi)展，以不同地貌分區(qū)為單元的研究還鮮有報(bào)道。該文基于多源衛(wèi)星遙感數(shù)據(jù)估算了1980 年、2000 年、2010 年、2016 年黃河流域植被覆蓋度，通過(guò)分析黃土高原重點(diǎn)研究區(qū)植被覆蓋度的時(shí)空變化特征識(shí)別植被變化背景下不同地貌分區(qū)侵蝕產(chǎn)沙風(fēng)險(xiǎn)變化，試圖回答對(duì)黃河減水減沙有重要貢獻(xiàn)的潛在區(qū)域。

1 研究區(qū)概況

中國(guó)黃土高原是地球上最集中最大的黃土地區(qū)，總面積達(dá) 64 萬(wàn)km2。 由于土壤松散、地形破碎、植被稀疏、溝壑縱橫，黃土高原是黃河流域主要的泥沙來(lái)源，也是世界上水土流失最為嚴(yán)重的和生態(tài)環(huán)境最脆弱的區(qū)域之一［18-19］。已有研究表明，黃土高原約 9.12 萬(wàn) km2土地土壤侵蝕速率超過(guò) 8 000 t/km2/yr［20］。為改善黃土高原水土流失現(xiàn)狀，自1999 年起，國(guó)家開(kāi)始提出并實(shí)施退耕還林（草）、荒山造林等一系列大規(guī)模生態(tài)修復(fù)工程，黃土高原地區(qū)成為政府首批“退耕還林還草”工程試點(diǎn)區(qū)域之一［4，21］。

該文以黃土高原中部區(qū)域?yàn)橹攸c(diǎn)研究區(qū)（圖1），該區(qū)域水土流失嚴(yán)重，是黃河流域的主要產(chǎn)沙區(qū)，面積約為20.4 萬(wàn)km2，氣候多為干旱半干旱類型。在地貌類型上，該重點(diǎn)研究區(qū)覆蓋了平原、黃土塬、黃土破碎塬、片沙黃土丘陵、峁?fàn)铧S土丘陵、梁狀黃土丘陵、寬谷長(zhǎng)梁丘陵、沙丘、風(fēng)蝕沙化丘陵、石質(zhì)山地、黃土低山、隴西丘陵12 種類型。

2 數(shù)據(jù)來(lái)源與方法

2.1 數(shù)據(jù)來(lái)源

研究所用的數(shù)據(jù)主要包括1 980 年前后MSS、2000 年前后TM、2010 年前后TM 和2016 年前后Landsat-8 OLI 多源遙感數(shù)據(jù)（表1）。這些數(shù)據(jù)均來(lái)源于美國(guó)地質(zhì)勘探局網(wǎng)站（https：//earthexplorer.usgs.gov/），獲取時(shí)間均為植被生長(zhǎng)季，較好代表了一年中植被生長(zhǎng)狀況。此外，還使用到1∶100 萬(wàn)土壤類型圖、四期（1980 年、2000 年、2010 年、2016 年）目視解譯的土地利用圖和黃土高原地貌分區(qū)圖，分別輔助植被蓋度信息提取和區(qū)域分析。

圖1 黃土高原重點(diǎn)研究區(qū)域空間位置Fig.1 Location of the key study area on the Loess Plateau

表1 遙感數(shù)據(jù)成像信息Table 1 Information of the data sources

2.2 方法

景可［22］認(rèn)為無(wú)論植被類型、降雨條件和其它下墊面條件如何，當(dāng)植被覆蓋度大于70%，地表侵蝕量都極其微弱。該文運(yùn)用≥70%植被蓋度面積變化指數(shù)揭示對(duì)地表侵蝕產(chǎn)沙有較強(qiáng)控制作用的林草植被變化情況，公式為：

式（1）中，VC≥70%表示≥70%植被蓋度變化指數(shù)；2016VC≥70%表示2016 年植被蓋度≥70%的面積百分比；1980VC≥70%表示1980 年植被蓋度≥70%的面積百分比。VC≥70%越大，表示對(duì)地表侵蝕產(chǎn)沙有較強(qiáng)控制作用的林草植被面積增加越多。

對(duì)于植被蓋度信息提取，由于歸一化植被指數(shù)（NDVI）與植被長(zhǎng)勢(shì)、生物量、蓋度和葉面積指數(shù)等有較強(qiáng)相關(guān)性，而且能部分消除輻照條件變化對(duì)反演參數(shù)的影響，該文采用式（2）所示像元二分模型反演植被蓋度［23-24］。

式（2）中，VC為像元植被蓋度；NDVI為像元?dú)w一化植被指數(shù)值；NDVIsoil為裸土或無(wú)植被覆蓋區(qū)NDVI 值；NDVIveg為完全被植被覆蓋區(qū)NDVI 值；NIR為近紅外波段反射值；R為紅光波段反射值。

為保證多時(shí)相遙感數(shù)據(jù)NDVI 有可比性，原始影像數(shù)據(jù)經(jīng)輻射定標(biāo)、簡(jiǎn)單大氣校正和幾何精校正后，生成的NDVI 用式3 所示相鄰圖幅遞進(jìn)回歸分析法消除時(shí)相差異［25］。

式（3）中，Image′i為校正后影像數(shù)據(jù)矩陣，Imagei為原始影像數(shù)據(jù)矩陣，ki-1和qi-1分別為影像數(shù)據(jù)矩陣Image′i-1⌒Imagei的回歸分析斜率與截距，n為影像圖幅數(shù)量。

由于受大氣效應(yīng)、地表水分和植被類型等影響，NDVIsoil和NDVIveg會(huì)隨時(shí)間和空間而改變。因此，使用土壤類型圖分圖斑統(tǒng)計(jì)NDVI 的累積頻率，選取土種單元內(nèi)累積頻率為5%的NDVI 值作為NDVIsoil；使用土地利用圖分類型統(tǒng)計(jì)NDVI 的累積頻率，選取林地和草地累積頻率為95%的NDVI 值作為NDVIveg。

3 結(jié)果與討論

3.1 植被覆蓋變化特征

3.1.1 植被蓋度均值變化特征

使用像元二分模型計(jì)算得到1980 年、2000 年、2010 年和2016 年黃土高原重點(diǎn)研究區(qū)植被蓋度。經(jīng)多次野外考察GPS 定位拍攝典型樣地照片，在植被蓋度無(wú)顯著變化區(qū)域內(nèi)，用照相法植被覆蓋度動(dòng)態(tài)監(jiān)測(cè)系統(tǒng)獲取樣地植被蓋度。其中，“十二五”期間獲取126 個(gè)樣地植被蓋度，2015—2017 年4 次野外考察獲取78 個(gè)樣地植被蓋度。檢驗(yàn)植被蓋度遙感反演精度達(dá)90.18%，表明遙感反演植被蓋度結(jié)果能滿足研究需求。統(tǒng)計(jì)得到黃土高原重點(diǎn)研究區(qū)1980 年全區(qū)植被蓋度均值為34.99%、2000 年為36.24%、2010 年為58.86%、2016 年為63.64%。表明黃土高原重點(diǎn)研究區(qū)植被蓋度均值呈增加趨勢(shì)，尤其是2000 年退耕還林（草）工程實(shí)施后增加的速度更為顯著，且2010—2016 年間增速相對(duì)于2000—2010 年間較緩（圖 2）。

3.1.2 植被蓋度分級(jí)變化特征

為了集中反映黃土高原重點(diǎn)研究區(qū)植被變化特征，將研究區(qū)1980 年、2000 年、2010 年和2016 年植被蓋度值從0～100%等距分為10 個(gè)等級(jí)，用GIS 空間分析方法分別統(tǒng)計(jì)4 個(gè)時(shí)期的植被蓋度等級(jí)變化，如圖3 和圖4 所示。低等級(jí)植被蓋度主要分布在研究區(qū)的北部和西北區(qū)域，高等級(jí)植被蓋度主要分布在研究區(qū)東南部，呈現(xiàn)出由西北至東南方向的植被狀況變好的空間趨勢(shì)。

圖2 黃土高原重點(diǎn)研究區(qū)全區(qū)植被蓋度均值統(tǒng)計(jì)Fig.2 Statistical characteristics of the mean value of vegetation coverage in the whole region

圖3 黃土高原重點(diǎn)研究區(qū)植被蓋度空間分布Fig.3 Spatial distribution of vegetation coverage in key study area on Loess Plateau

統(tǒng)計(jì)分析研究區(qū)各等級(jí)植被蓋度面積占比情況由圖4 所示，1980 年植被蓋度值主要集中在0～50%，共占全區(qū)面積比例為80.14%；2000 年植被蓋度值主要集中在10%～60%，共占全區(qū)面積比例為86.47%；2010 年植被蓋度值主要集中在30%～80%，共占全區(qū)面積比例為68.46%；2016 年植被蓋度值主要集中在40%～90%，共占全區(qū)面積比例為74.56%。可以看出從1980—2016 年，峰值植被蓋度等級(jí)由1980 年和2000 年的20%～30%過(guò)渡到2010年的50%～60%，再到2016 年的60%～70%，呈現(xiàn)出由低等級(jí)向高等級(jí)逐漸過(guò)渡的趨勢(shì)；低值植被蓋度區(qū)間0～20%占全區(qū)面積百分比呈現(xiàn)出持續(xù)減小的趨勢(shì)，高值植被蓋度區(qū)間50%～100%呈現(xiàn)出增加的趨勢(shì)，較低值植被蓋度區(qū)間20%～40%呈現(xiàn)出先增加后銳減的趨勢(shì)。總的來(lái)說(shuō)，黃土高原重點(diǎn)研究區(qū)植被蓋度等級(jí)向著植被持續(xù)恢復(fù)的趨勢(shì)發(fā)展。

圖4 黃土高原重點(diǎn)研究區(qū)全區(qū)植被蓋度分級(jí)統(tǒng)計(jì)Fig.4 Statistical characteristics of different vegetation coverage intervals in the whole region

3.2 產(chǎn)沙風(fēng)險(xiǎn)等級(jí)變化特征

該文根據(jù)野外實(shí)地調(diào)查確定產(chǎn)沙風(fēng)險(xiǎn)等級(jí)劃分標(biāo)準(zhǔn)如表2 所示。 黃土高原重點(diǎn)研究區(qū)1980 年、2000 年、2010 年和2016 年4 期產(chǎn)沙風(fēng)險(xiǎn)等級(jí)如圖5 所示。可以看出4 期產(chǎn)沙風(fēng)險(xiǎn)等級(jí)圖均是產(chǎn)沙高風(fēng)險(xiǎn)和易產(chǎn)沙區(qū)主要分布在研究區(qū)的北部和西北部，不產(chǎn)沙區(qū)和產(chǎn)沙微弱區(qū)主要分布在研究區(qū)東南部。

表2 產(chǎn)沙風(fēng)險(xiǎn)等級(jí)劃分標(biāo)準(zhǔn)Table 2 Criteria for classification of sediment yield risk

各產(chǎn)沙風(fēng)險(xiǎn)等級(jí)類型面積占比變化趨勢(shì)如圖6 所示。4 個(gè)時(shí)期均為易產(chǎn)沙區(qū)面積比例最高，1980 年和2000 年易產(chǎn)沙區(qū)面積比例均遠(yuǎn)高于產(chǎn)沙微弱區(qū)；2010 年和2016 年產(chǎn)沙微弱區(qū)面積比例與易產(chǎn)沙區(qū)面積比例差距逐漸縮小。1980—2016 年，產(chǎn)沙高風(fēng)險(xiǎn)區(qū)面積比例持續(xù)減小；易產(chǎn)沙區(qū)面積比例先小幅增加，后顯著減小；產(chǎn)沙微弱區(qū)面積比例先小幅減小，后顯著增加；不產(chǎn)沙區(qū)面積比例先小幅減小，后增大至穩(wěn)定。表明由于退耕還林（草）工程實(shí)施，植被狀況改善，研究區(qū)產(chǎn)沙高風(fēng)險(xiǎn)區(qū)面積比例顯著減小。

圖5 產(chǎn)沙風(fēng)險(xiǎn)等級(jí)分區(qū)圖Fig.5 Spatial distribution of different sediment yield risk zones

圖6 各產(chǎn)沙風(fēng)險(xiǎn)等級(jí)類型面積占比變化趨勢(shì)Fig.6 The trend of the proportion of the area of each sediment yield risk class

3.3 植被蓋度變化的地貌差異

3.3.1 植被蓋度變化特征

圖7 展示了不同地貌區(qū)的平均植被蓋度變化情況。4 個(gè)時(shí)期均為黃土低山地貌區(qū)植被蓋度均值最高，植被蓋度均值均在70%以上；2000 年、2010 年和2016 年均為風(fēng)蝕沙化丘陵植被蓋度均值最低；2010 年和2016 年黃土塬和黃土破碎塬植被蓋度均值均在70%以上；2016 年峁?fàn)铧S土丘陵地貌區(qū)植被蓋度也達(dá)到了70%以上；1980 年和2000 年除黃土低山外的所有地貌類型植被蓋度均值均在70%以下；各地貌區(qū)均呈現(xiàn)2016 年植被蓋度均值最高，2010 年植被蓋度均值次之。表明自1999 年實(shí)施退耕還林（草）工程以來(lái)，不同地貌分區(qū)的植被狀況均持續(xù)好轉(zhuǎn)。

圖7 不同地貌區(qū)平均植被蓋度變化Fig.7 Variations of the mean value of vegetation coverage in different landforms

3.3.2 ≥70%植被蓋度變化

圖8 和圖9 展示了不同產(chǎn)沙風(fēng)險(xiǎn)等級(jí)在不同地貌區(qū)的歷年絕對(duì)面積變化情況和不同地貌分區(qū)≥70%植被蓋度變化指數(shù)。對(duì)于不產(chǎn)沙區(qū)（VC＞90），黃土低山該類型面積比例較高，2010 年和2016 年面積占比均達(dá)到了70%；其余所有地貌區(qū)內(nèi)該類型面積比例均較低，在30%以下。對(duì)于產(chǎn)沙微弱區(qū)（VC ≥70），黃土低山、黃土塬、黃土破碎塬、峁?fàn)铧S土丘陵、梁狀黃土丘陵、隴西丘陵該類型面積比例顯著增加，至2016 年均達(dá)到了40%以上。對(duì)于易產(chǎn)沙區(qū)（VC＜70），除黃土低山外，其余地貌區(qū)1980 年和2000 年該類型面積比例均在80%以上，說(shuō)明這兩個(gè)時(shí)期研究區(qū)整體產(chǎn)沙風(fēng)險(xiǎn)很高；至2016 年，各地貌類型區(qū)該類型面積比例均呈現(xiàn)下降的趨勢(shì)，其中黃土塬、黃土破碎塬、峁?fàn)铧S土丘陵、隴西丘陵降幅比較大，均在40%以上。對(duì)于產(chǎn)沙高風(fēng)險(xiǎn)區(qū)（VC＜30），1980 年和2000 年平原區(qū)、片沙黃土丘陵、沙丘、風(fēng)蝕沙化丘陵該類型占比比較高，約在60%以上；各地貌類型區(qū)植被狀況均有所改善，至2016 年，該類型面積比例降幅比較大的有平原區(qū)、片沙黃土丘陵、寬谷長(zhǎng)梁丘陵、沙丘和風(fēng)蝕沙化丘陵，約在50%以上。

圖9 可知，1980—2016 年，≥70%植被蓋度變化指數(shù)最高的為黃土塬，其次為黃土破碎塬，最低的為沙丘。≥70%植被蓋度變化指數(shù)小于30%的地貌區(qū)排序?yàn)椋荷城穑硷L(fēng)蝕沙化丘陵＜平原＜黃土低山＜片沙黃土丘陵。由于黃土低山4 期植被蓋度均值均在70%以上，因此沙丘、風(fēng)蝕沙化丘陵、平原、片沙黃土丘陵地貌區(qū)是未來(lái)水土保持工作需要關(guān)注的區(qū)域，植被改善工作仍需加強(qiáng)。

圖8 各地貌類型區(qū)不同產(chǎn)沙風(fēng)險(xiǎn)等級(jí)面積變化趨勢(shì)Fig.8 Trends of variations of area of different sediment yield risk classes in each landform

圖9 1980—2016 年不同地貌區(qū)≥70%植被蓋度變化指數(shù)（VC≥70%）Fig.9 Variations of VC≥70% index in different landforms during the study period

4 結(jié)論與展望

基于多源遙感數(shù)據(jù)利用像元二分模型可實(shí)現(xiàn)大區(qū)域尺度長(zhǎng)時(shí)間序列的植被覆蓋動(dòng)態(tài)變化監(jiān)測(cè)。通過(guò)監(jiān)測(cè)分析黃土高原重點(diǎn)研究區(qū)1980 年、2000 年、2010 年及2016 年的植被覆蓋度狀況，得到結(jié)論如下。

（1）黃土高原重點(diǎn)研究區(qū)全區(qū)植被蓋度均值呈增加趨勢(shì)，由1980 年的34.99%提升至2016 年的63.64%，且2000 年退耕還林（草）工程實(shí)施后增加的速度更為顯著，2010—2016 年間增速相對(duì)于2000—2010 年間較緩。在空間上，植被狀況呈現(xiàn)出由西北至東南方向逐漸變好的趨勢(shì)。

（2）4 期產(chǎn)沙風(fēng)險(xiǎn)等級(jí)圖均是產(chǎn)沙高風(fēng)險(xiǎn)和易產(chǎn)沙區(qū)主要分布在研究區(qū)的北部和西北部，不產(chǎn)沙區(qū)和產(chǎn)沙微弱區(qū)主要分布在研究區(qū)東南部。由于退耕還林（草）工程實(shí)施，植被狀況改善，研究區(qū)產(chǎn)沙高風(fēng)險(xiǎn)區(qū)面積比例顯著減小。

（3）自1999 年實(shí)施退耕還林（草）工程以來(lái)，不同地貌分區(qū)的植被狀況均持續(xù)好轉(zhuǎn)。沙丘、風(fēng)蝕沙化丘陵、平原、片沙黃土丘陵地貌區(qū)是未來(lái)水土保持工作需要關(guān)注的區(qū)域，植被改善工作仍需加強(qiáng)。

相關(guān)研究表明黃土高原不同地貌類型區(qū)遏制侵蝕產(chǎn)沙的林草植被覆蓋率存在不同閾值［17］。此外，王萬(wàn)忠基于水文站年輸沙量數(shù)據(jù)，通過(guò)對(duì)比分析近50 年來(lái)黃土高原治理前、一期治理、二期治理和退耕后（ 2000—2009 年）4 個(gè)階段侵蝕產(chǎn)沙的時(shí)間變化特征，表明受坡改梯、退耕還林（草）、淤地壩建設(shè)等多種水土保持措施的綜合影響，侵蝕產(chǎn)沙量和減水減沙效益均表現(xiàn)出不同的區(qū)域分布特征［26］。該文采用了不同于上述研究的空間尺度，受研究觀測(cè)資料的限制，進(jìn)行了概化處理，主要包括：（1）在整個(gè)黃土高原重點(diǎn)研究區(qū)均按照景可等人提出的70%的植被覆蓋度標(biāo)準(zhǔn)來(lái)劃分黃土高原重點(diǎn)研究區(qū)侵蝕產(chǎn)沙風(fēng)險(xiǎn)等級(jí)；（2）僅考慮了植被覆蓋度這一指標(biāo)，分析了植被恢復(fù)背景下黃土高原重點(diǎn)研究區(qū)侵蝕產(chǎn)沙風(fēng)險(xiǎn)的時(shí)空變化。這些概化處理勢(shì)必會(huì)對(duì)黃土高原重點(diǎn)研究區(qū)侵蝕產(chǎn)沙風(fēng)險(xiǎn)的區(qū)域時(shí)空分布特征造成一定影響，在未來(lái)的研究中需要構(gòu)建考慮植被減沙閾值差異性、多水土保持措施綜合影響下的區(qū)域侵蝕產(chǎn)沙風(fēng)險(xiǎn)評(píng)價(jià)方法，全面評(píng)估今后黃土高原治理的重點(diǎn)區(qū)域。