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

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

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

1 研究區(qū)概況

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

2 數據與方法

2.1 數據來源

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

2.2 中國坡面水蝕預報模型因子計算

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

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因子利用幾何平均粒徑結合有機質模型計算[15]。C因子利用蔡崇法等[16]和林杰等[17]的用NDVI推導的模型計算。我國2014—2020 年完成退耕283萬hm2[9]，在2000—2017年NDVI數據基礎上，用最小二乘法實現對未來NDVI的預測。對C因子分析采用RS時序分析和標準差分析方法,用Hurst指數表示[18]。淺溝侵蝕G因子的利用江忠善等[11]的研究成果。GeoSOS-FLUS軟件是根據FLUS模型原理開發(fā)的多類土地利用變化情景模擬軟件。GeoSOS-FLUS利用神經網絡算法[19](artificial neuron net)獲取各類用地的適宜性概率，再通過耦合系統(tǒng)動力學模型(system dynamics) 和元胞自動機(cellular automata) 模型以提高模型的適用性，在國內外土地利用變化模擬預測中得到了較好地應用[20]。文中P因子的取值主要參照謝紅霞[21]的研究成果。文中所有因子柵格數據，重采樣為30 m分辨率柵格數據再進行分析計算。

3 結果與分析

3.1 主要因子的結果

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)之間，西部風沙區(qū)定邊縣、靖邊縣在800～1 050 MJ·mm/(hm2·h)之間，北部地區(qū)在1 000～1 150 MJ·mm/(hm2·h)之間。圖3為多年平均降雨侵蝕力標準差空間分布，標準差為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世紀的前13年和20世紀末的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 多年平均降雨侵蝕力標準差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分布在總體也表現出西部、北部大于中部、東部的特征。2000年耕地LS均值由4.55減少為2013年的4.43，還林還草工程后，坡耕地退耕植樹種草，LS因子呈現下降趨勢，與現實相符合。2000 年林地LS均值由4.04增加到2013年的4.31；2000年草地LS均值由3.99增加到2013年的4.09。林地和草地的LS因子增大，也和實際情況相符。水域2000年LS均值由2.11變?yōu)?013年的2.06；建設用地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因子總體呈現出波動下降的遞減趨勢。2000年最大為0.163，2016年最小為0.067，2015年的C因子異常偏大，和2015年的降水偏少有關。C值隨時間的減少達到非常顯著的水平。相關研究認為1 km2的分辨率的GIMMS/NDVI數據在黃土高原地區(qū)，2000年NDVI和1984、1994年相近，都在0.250以下[22]。2000年以后的MODIS NDVI，相比GIMMS/NDVI表現出更高的靈敏度，二者在黃土高原對植被的監(jiān)測中表現出一致性。由MODIS NDVI計算的2000年C因子可以基本代表1984和1994年的C因子。圖5是利用最小二乘法預測的C因子空間分布，預測C因子為0.053，比2017年C因子的0.069，還有0.017的下降空間。

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

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

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

圖6 預測 C因子的Hurst指數Fig.6 Predicting the Hurst index of C factor

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

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

圖8 用神經網絡網絡模型預測土地利用類型Fig.8 Predicting land use types by artificial neuron net

3.2 土壤侵蝕預測分析

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

圖9 2000—2017年土壤侵蝕模數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)。最大降雨侵蝕力情況下的土壤侵蝕模數為1 752.68 t/(km2·a)；最小降雨侵蝕力情況下的土壤侵蝕模數為453.2 t/(km2·a)；多年平均降雨侵蝕力情況下的土壤侵蝕模數為793.27 t/(km2·a)。在這種預測模式下，除R因子之外，其余因子都是定值。3種不同情況下的圖形線性分布呈現為增加或者減少趨勢一致。在極大降雨侵蝕力情況下，定邊縣土壤侵蝕模數超過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 最小/最大和多年平均降雨侵蝕力下預測2022年的土壤侵蝕模數Fig.10 Predicted soil erosion moduli in 2022 under the min/max and multi-year average rainfall erosivity

在多年平均降雨侵蝕力條件下預測2022年各區(qū)縣土壤侵蝕模數：神木市為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年榆林市大多數區(qū)縣的土壤侵蝕模數將<2 000 t/(km2·a)，侵蝕以輕度和微度為主。

為更加直觀地分析未來土壤侵蝕情況，繪制了多年平均降雨侵蝕力下各區(qū)縣2000、2010、2016年和2022年各區(qū)縣的土壤侵蝕模數(圖11)。2016年和2022年預測的土壤侵蝕模數相比較。神木市減少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 多年平均降雨侵蝕力下土壤侵蝕模數比較Fig.11 Comparison of soil erosion modulus under multi-year annual average rainfall erosivity

4 討論

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

5 結論

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

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

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