退耕還林(草)背景下榆林市土壤侵蝕變化與未來趨勢預(yù)測

楊 波,王全九,周 佩，許曉婷,黨江茹

(1.咸陽師范學(xué)院資源環(huán)境與歷史文化學(xué)院，712000,陜西咸陽；2.西安理工大學(xué)水利水電學(xué)院,710048,西安;3.福建師范大學(xué)地理科學(xué)學(xué)院,350007,福州)

還林還草工程作為我國西部地區(qū)重要的生態(tài)保護措施，已經(jīng)實施18年。取得明顯的生態(tài)效益[1]，植被迅速恢復(fù)[2]。在大、中、小不同尺度下的研究表明：黃土高原地區(qū)土壤侵蝕得到有效的控制[3-5]，但是植被恢復(fù)后，耗水量增加[6]，近50年來該地區(qū)的降水呈現(xiàn)減少趨勢[7-8]，土壤侵蝕和植被恢復(fù)關(guān)系密切[9]，水熱條件對大規(guī)模的草木生長有一定的限制。未來植被的增長潛力與土壤侵蝕密切相關(guān)，目前黃土高原地區(qū)土壤侵蝕潛力和預(yù)測研究較少。在中國坡面水蝕預(yù)報模型計算2000—2017年土壤侵蝕的基礎(chǔ)之上，筆者利用最小二乘法預(yù)測植被指數(shù)、神經(jīng)網(wǎng)絡(luò)方法預(yù)測在退耕還林和建設(shè)能源化工型城市的雙重背景下的土地利用類型，預(yù)測榆林市土壤侵蝕潛力和變化情況，為該地區(qū)的水土保持和生態(tài)環(huán)境建設(shè)工作提供參考。

1 研究區(qū)概況

榆林市位于陜西省最北部，西鄰甘肅環(huán)縣、寧夏鹽池縣，北連內(nèi)蒙古準(zhǔn)格爾、伊金霍洛、烏審、鄂托克等4旗，東隔黃河與山西相望，南與陜西省延安市接壤。位于E 107°28′～111°15′、N 36°57′～39°35′。榆林市行政區(qū)劃版圖形似三角形。東西最大長度309 km，南北最大寬度295 km，總面積4萬3 578 km2，約占陜西省的21%。地貌主要有風(fēng)沙草灘區(qū)、黃土丘陵溝壑區(qū)、梁狀低山丘陵區(qū)3大類。該區(qū)域氣候?qū)倥瘻貛Ш蜏貛О敫珊荡箨懶约撅L(fēng)氣候，四季分明，無霜期短，年平均氣溫10 ℃，年平均降水量400 mm左右[9]。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

本研究使用的數(shù)據(jù)有： 1)陜西省2000、2005、2008和2013年1∶10萬土地利用數(shù)據(jù),空間分辨率30 m，來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn/)；2)土壤數(shù)據(jù)來源于《陜西省第2次土壤普查數(shù)據(jù)集》《陜西土壤》和中國土壤數(shù)據(jù)集(V1.1)，其中中國土壤數(shù)據(jù)集(V1.1)來源于(http:∥westdc.westgis.ac.cn/)；3)DEM數(shù)據(jù)分辨率為30 m，2000—2017年NDVI(normalized difference vegetation index)分辨率為250 m，數(shù)據(jù)來源于(http:∥earthexplorer.usgs.gov/)； 4)2000—2017年榆林榆陽區(qū)、綏德縣、定邊縣、橫山區(qū)、神木市等13個站點氣象站。

2.2 中國坡面水蝕預(yù)報模型因子計算

目前RUSLE土壤侵蝕模型和RS、GIS已經(jīng)深度融合。該模型在黃土高原地區(qū)不同尺度的土壤侵蝕研究中應(yīng)用廣泛，如羊圈溝小流域，面積約2.02 km2[3]； 榆林市，面積約4萬3 578 km2[10]，黃土高原整體，面積約64萬km2[5]。由于該地區(qū)細溝侵蝕比較嚴(yán)重，RUSLE估算值偏小，江忠善等[11]將淺溝侵蝕因子加入到RUSLE模型中，提高估算精度[10]。建立中國坡面水蝕預(yù)報模型：

A=RLSKCPG。

(1)

式中：A為土壤流失量，t/(km2·a)；R為降雨侵蝕力因子，MJ·mm/(hm2·h)；L為坡長因子，量綱為1；S為坡度因子，量綱為1；K為土壤可蝕性因子，t·hm2·h/ (hm2·MJ·mm)；C為植被覆蓋與管理因子，量綱為1；P為水土保持措施因子，量綱為1；G為淺溝侵蝕影響因子，量綱為1。

降雨侵蝕力R因子，依照謝云等[12]和章文波等[13]研究成果計算。LS因子是地形對土壤侵蝕的影響，用符素華等[14]提供的基于DEM的LS計算工具計算。K因子利用幾何平均粒徑結(jié)合有機質(zhì)模型計算[15]。C因子利用蔡崇法等[16]和林杰等[17]的用NDVI推導(dǎo)的模型計算。我國2014—2020 年完成退耕283萬hm2[9]，在2000—2017年NDVI數(shù)據(jù)基礎(chǔ)上，用最小二乘法實現(xiàn)對未來NDVI的預(yù)測。對C因子分析采用RS時序分析和標(biāo)準(zhǔn)差分析方法,用Hurst指數(shù)表示[18]。淺溝侵蝕G因子的利用江忠善等[11]的研究成果。GeoSOS-FLUS軟件是根據(jù)FLUS模型原理開發(fā)的多類土地利用變化情景模擬軟件。GeoSOS-FLUS利用神經(jīng)網(wǎng)絡(luò)算法[19](artificial neuron net)獲取各類用地的適宜性概率，再通過耦合系統(tǒng)動力學(xué)模型(system dynamics) 和元胞自動機(cellular automata) 模型以提高模型的適用性，在國內(nèi)外土地利用變化模擬預(yù)測中得到了較好地應(yīng)用[20]。文中P因子的取值主要參照謝紅霞[21]的研究成果。文中所有因子?xùn)鸥駭?shù)據(jù)，重采樣為30 m分辨率柵格數(shù)據(jù)再進行分析計算。

3 結(jié)果與分析

3.1 主要因子的結(jié)果

3.1.1 降雨侵蝕力因子 圖1為榆林市1988—2017年降雨侵蝕力因子分析，2015年最小為262.85 MJ·mm/(hm2·h)，2001年最大為2 394.59 MJ·mm/(hm2·h)，多年平均降雨侵蝕力均值為1 088.54 MJ·mm/(hm2·h)。圖2為1988—2017年榆林市多年平均降雨侵蝕力。丘陵溝壑區(qū)的佳縣、米脂縣、吳堡縣和綏德縣在1 150～1 350 MJ·mm/(hm2·h)之間，西部風(fēng)沙區(qū)定邊縣、靖邊縣在800～1 050 MJ·mm/(hm2·h)之間，北部地區(qū)在1 000～1 150 MJ·mm/(hm2·h)之間。圖3為多年平均降雨侵蝕力標(biāo)準(zhǔn)差空間分布，標(biāo)準(zhǔn)差為682.44 MJ·mm/(hm2·h)，分布區(qū)間為(461.70，950.42)MJ·mm/(hm2·h)。1988—2000年R因子均值為998.67 MJ·mm/(hm2·h)，2001—2017年R因子均值為1 266.91 MJ·mm/(hm2·h)。21世紀(jì)的前13年和20世紀(jì)末的17年相比，降雨侵蝕力增加268.24 MJ·mm/(hm2·h)。

圖1 1988—2017年降雨侵蝕力因子Fig.1 Rainfall erosivity factors from 1988 to 2017

圖2 多年平均降雨侵蝕力Fig.2 Multi-year average rainfall erosivity

圖3 多年平均降雨侵蝕力標(biāo)準(zhǔn)差Fig.3 Standard deviation of multi-year average rainfall erosivity

3.1.2 坡長坡度因子LS因子分布在(0，80.29)之間，均值為3.98，榆陽區(qū)最低為1.57，清澗縣最高為6.72。西部和北部地區(qū)比中東部地勢平坦，LS分布在總體也表現(xiàn)出西部、北部大于中部、東部的特征。2000年耕地LS均值由4.55減少為2013年的4.43，還林還草工程后，坡耕地退耕植樹種草，LS因子呈現(xiàn)下降趨勢，與現(xiàn)實相符合。2000 年林地LS均值由4.04增加到2013年的4.31；2000年草地LS均值由3.99增加到2013年的4.09。林地和草地的LS因子增大，也和實際情況相符。水域2000年LS均值由2.11變?yōu)?013年的2.06；建設(shè)用地2000年LS均值由2.17變?yōu)?013年的1.91；未利用土地LS未發(fā)生變化。

3.1.3 植被覆蓋與管理因子 圖4是2000—2017年植被覆蓋與管理因子分析，2000年和2017年的C因子均值分別為0.164和0.069，降低136.75%。18年來C因子總體呈現(xiàn)出波動下降的遞減趨勢。2000年最大為0.163，2016年最小為0.067，2015年的C因子異常偏大，和2015年的降水偏少有關(guān)。C值隨時間的減少達到非常顯著的水平。相關(guān)研究認(rèn)為1 km2的分辨率的GIMMS/NDVI數(shù)據(jù)在黃土高原地區(qū)，2000年NDVI和1984、1994年相近，都在0.250以下[22]。2000年以后的MODIS NDVI，相比GIMMS/NDVI表現(xiàn)出更高的靈敏度，二者在黃土高原對植被的監(jiān)測中表現(xiàn)出一致性。由MODIS NDVI計算的2000年C因子可以基本代表1984和1994年的C因子。圖5是利用最小二乘法預(yù)測的C因子空間分布，預(yù)測C因子為0.053，比2017年C因子的0.069，還有0.017的下降空間。

圖4 植被覆蓋與管理因子Fig.4 Cover and managemt factor

圖5 最小二乘法預(yù)測C因子Fig.5 Predicting C factor via least square method

Hurst指數(shù)預(yù)測未來C因子的變化情況，Hurst指數(shù)介于(0.11, 0.98)之間(圖6)，均值為0.50。Hurst指數(shù)在0～0.500，表示未來C因子是反向變化，會降低；0.50～1.00之間是表示正向變化，C因子會增長。正向變化的和反向變化的區(qū)域分別達到51.53% 和48.74%。退耕還林(草)工程實施以后，丘陵溝壑植被覆蓋增加顯著，隨之而來的是蒸發(fā)量增加，夏季蒸散量約占全年蒸散量一半[23-24]，植樹造林的人工林土壤含水量普遍低于農(nóng)地，在降雨量>550 mm區(qū)域，土壤儲水量損失最高。植被恢復(fù)越好，需水量越大，森林需水與土壤供水之間的矛盾突出，急需解決[24]。未來極端干旱的年份，由于水資源的短缺，可能對植被生存產(chǎn)生威脅。

圖6 預(yù)測 C因子的Hurst指數(shù)Fig.6 Predicting the Hurst index of C factor

3.1.4 土地利用 結(jié)合榆林市“十四五”期間土地利用規(guī)劃，重點考慮退耕還林和城市發(fā)展對土地利用的影響，把>25°坡耕地在土地利用類型轉(zhuǎn)化成草地和灌木林地為高概率，對于平原地區(qū)的耕地在距離城市核心區(qū)>5 km的全部為禁止轉(zhuǎn)換。利用神經(jīng)網(wǎng)絡(luò)網(wǎng)絡(luò)模型在2013年土地類型基礎(chǔ)上(圖7)預(yù)測未來到2022年土地利用類型(圖8)，2013年耕地面積為1萬5 739.78 km2、林地面積為2 348.46 km2、草地面積為1萬8 908.40 km2、水域面積為396.23 km2、建設(shè)用地面積為529.55 km2和未利用土地面積為4 390.12 km2。預(yù)測2022年耕地面積為1萬4 468.05 km2、林地面積為3 036.88 km2、草地面積為1萬9 487.55 km2、水域面積為396.10 km2、建設(shè)用地面積為662.50 km2和未利用土地面積為4 261.56 km2。

圖7 2013年土地利用類型Fig.7 Land use types in 2013

圖8 用神經(jīng)網(wǎng)絡(luò)網(wǎng)絡(luò)模型預(yù)測土地利用類型Fig.8 Predicting land use types by artificial neuron net

3.2 土壤侵蝕預(yù)測分析

利用預(yù)測的C、P因子和K、LS、K、R因子估算未來土壤侵蝕。將最終得到的各年土壤侵蝕強度劃分為微度(015 000)6級[25]，單位為t/(km2·a)。2000—2017年土壤侵蝕模數(shù)如圖9a所示，相比年降雨侵蝕力下的土壤侵蝕計算，多年平均降雨侵蝕力下的土壤侵蝕估算實際意義更大，消除了年降雨侵蝕力差異引起的土壤侵蝕變化，評價更為客觀。圖9b是多年平均降雨侵蝕力下土壤侵蝕模數(shù)序列。2000—2017年多年平均降雨侵蝕力因子R均值為1 220.95 MJ·mm/(hm2·h)。土壤侵蝕模數(shù)總體呈現(xiàn)波動下降的趨勢。2000年土壤侵蝕模數(shù)最大，為3 559.99 t/(km2·a)；預(yù)測2022年最小，為793.27 t/(km2·a)。線性回歸方程為y=-104.54x+3 103.1，R2=0.726 7，達到顯著相關(guān)。2000—2017年多年平均降雨侵蝕力為1 101.02 MJ·mm/(hm2·h)。土壤侵蝕模數(shù)年際之間的差異，主要是植被覆蓋與管理因子C和水土保持措施因子P的差異造成的。

圖9 2000—2017年土壤侵蝕模數(shù)Fig.9 Soil erosion modulus in 2000—2017

由圖9a可見，2000年的降雨侵蝕力最小，為439.54 MJ·mm/(hm2·h)；2001年的降雨侵蝕力最大，為2 394.59 MJ·mm/(hm2·h)。利用最大降雨、最小降雨和多年平均降雨侵蝕力條件下模擬2022年土壤各區(qū)縣侵蝕情況(圖10)。最大降雨侵蝕力情況下的土壤侵蝕模數(shù)為1 752.68 t/(km2·a)；最小降雨侵蝕力情況下的土壤侵蝕模數(shù)為453.2 t/(km2·a)；多年平均降雨侵蝕力情況下的土壤侵蝕模數(shù)為793.27 t/(km2·a)。在這種預(yù)測模式下，除R因子之外，其余因子都是定值。3種不同情況下的圖形線性分布呈現(xiàn)為增加或者減少趨勢一致。在極大降雨侵蝕力情況下，定邊縣土壤侵蝕模數(shù)超過6 000 t/(km2·a)，子洲縣超過3 000 t/(km2·a)，綏德縣、米脂縣和清澗縣在2 000～3 000 t/(km2·a)之間。其余地區(qū)低于2 000 t/(km2·a)。在最小、最大和多年平均降雨侵蝕力下，榆林市總侵蝕量分別為1 900萬t、7 370萬t和3 370萬t。

圖10 最小/最大和多年平均降雨侵蝕力下預(yù)測2022年的土壤侵蝕模數(shù)Fig.10 Predicted soil erosion moduli in 2022 under the min/max and multi-year average rainfall erosivity

在多年平均降雨侵蝕力條件下預(yù)測2022年各區(qū)縣土壤侵蝕模數(shù)：神木市為306.32 t/(km2·a)、米脂縣為738.37 t/(km2·a)、佳縣為425.93 t/(km2·a)、綏德縣為1 215.81 t/(km2·a)、吳堡縣為793.23 t/(km2·a)、靖邊縣為1 115.91 t/(km2·a)、子洲縣為1 556.88 t/(km2·a)、定邊縣為2 133.88 t/(km2·a)、府谷縣為567.88 t/(km2·a)、榆陽區(qū)為240.09 t/(km2·a)、橫山縣407.22 t/(km2·a)、清澗縣為1 599.75 t/(km2·a)。2022年榆林市大多數(shù)區(qū)縣的土壤侵蝕模數(shù)將<2 000 t/(km2·a)，侵蝕以輕度和微度為主。

為更加直觀地分析未來土壤侵蝕情況，繪制了多年平均降雨侵蝕力下各區(qū)縣2000、2010、2016年和2022年各區(qū)縣的土壤侵蝕模數(shù)(圖11)。2016年和2022年預(yù)測的土壤侵蝕模數(shù)相比較。神木市減少644.71 t/(km2·a)、米脂縣減少174.86 t/(km2·a)、佳縣減少812.1 t/(km2·a)、綏德縣減少299.35 t/(km2·a)、吳堡縣減少1 066.58 t/(km2·a)、靖邊縣減少659.18 t/(km2·a)、子洲縣減少661.18 t/(km2·a)、府谷縣減少510.36 t/(km2·a)、榆林市區(qū)減少469.7 t/(km2·a)、橫山縣減少592.76 t/(km2·a)、清澗縣減少1 119.92 t/(km2·a)。定邊縣增加464.81 t/(km2·a)。到2022年榆林市大部分地區(qū)的土壤侵蝕量還有進一步下降的空間。

圖11 多年平均降雨侵蝕力下土壤侵蝕模數(shù)比較Fig.11 Comparison of soil erosion modulus under multi-year annual average rainfall erosivity

4 討論

植被覆蓋的大小，直接決定土壤侵蝕，NDVI越大，土壤侵蝕越小[5]。40%～60% 植被覆蓋度是植被對土壤侵蝕起明顯干預(yù)作用的分界區(qū)間[26]，筆者預(yù)測榆林市植被覆蓋度可以達到62.14%，目前植被覆蓋度為52.20%，NDVI指數(shù)還能增加約0.99，榆林市植被恢復(fù)已經(jīng)接近潛力值。還林還草后水土保持措施對土壤侵蝕的減少起到主要作用，降雨對土壤侵蝕影響已處于次要的地位[1]。大規(guī)模的植被生長使耗水增加，加劇了土壤干層，土壤凈貯水量下降[27]，植被恢復(fù)越好，需水量越大，Hurst指數(shù)表明未來，在丘陵溝壑區(qū)植被恢復(fù)較好的地區(qū)，森林需水與土壤供水之間的矛盾可能會更加突出，在極端干旱的年份，水資源的短缺，可能對植被生存產(chǎn)生威脅。模擬條件下，干化土壤中植被生長受當(dāng)年降水量影響較大[28]，植被迅速恢復(fù)過程中蒸發(fā)的增加是地表徑流減少的主要原因之一[29]，植被恢復(fù)和生長受控于降雨條件限制[30]。下一階段的還林還草應(yīng)該“因水制宜”，避免加劇土壤干化。今后的還林還草工程要針對不同地區(qū)的降雨、土壤差異，合理選擇草木品種。土壤侵蝕的進一步減少，還要進一步合理調(diào)配水資源的可持續(xù)利用，以保障還林還草工程取得生態(tài)和水土保持效益。

5 結(jié)論

1)1988—2017年多年平均降雨侵蝕力均值約為1 088.54 MJ·mm/(hm2·h)，丘陵溝壑區(qū)的降雨侵蝕力在1 150～1 350 MJ·mm/(hm2·h)之間，21世紀(jì)的前13年和20世紀(jì)末的后17年相比，降雨侵蝕力增加268.24 MJ·mm/(hm2·h)。

2)還林還草工程實施18年來，植被覆蓋持續(xù)改善，植被覆蓋與管理因子C顯著降低。2000—2017年多年平均降雨侵蝕力下，2000年的土壤侵蝕模數(shù)最大，為3 559.99 t/(km2·a)，2017年土壤侵蝕模數(shù)為1 369.19 t/(km2·a)，土壤侵蝕量減少958.65萬t。預(yù)測2022年土壤侵蝕模數(shù)為793.27 t/(km2·a)，土壤侵蝕總體呈現(xiàn)波動下降的趨勢。還林還草工程顯著減少榆林市的水土流失。

3)在最大、最小和和多年平均降雨侵蝕力條件下模擬計算未來土壤侵蝕，最大降雨侵蝕力條件下的土壤侵蝕模數(shù)為1 752.68 t/(km2·a)；最小降雨侵蝕力情況下的土壤侵蝕模數(shù)僅為453.20 t/(km2·a)。未來榆林市大多數(shù)區(qū)縣的土壤侵蝕模數(shù)將<2 000 t/(km2·a)，土壤侵蝕以輕度和微度為主。土壤侵蝕模數(shù)約為793.27 t/(km2·a)。在最小、最大和多年平均降雨侵蝕力下，榆林市總侵蝕量分別為1 900萬t、7 370萬t和3 370萬t。