还林还草工程作为我国西部地区重要的生态保护措施，已经实施18年。取得明显的生态效益^[1]，植被迅速恢复^[2]。在大、中、小不同尺度下的研究表明：黄土高原地区土壤侵蚀得到有效的控制^[3-5]，但是植被恢复后，耗水量增加^[6]，近50年来该地区的降水呈现减少趋势^[7-8]，土壤侵蚀和植被恢复关系密切^[9]，水热条件对大规模的草木生长有一定的限制。未来植被的增长潜力与土壤侵蚀密切相关，目前黄土高原地区土壤侵蚀潜力和预测研究较少。在中国坡面水蚀预报模型计算2000—2017年土壤侵蚀的基础之上，笔者利用最小二乘法预测植被指数、神经网络方法预测在退耕还林和建设能源化工型城市的双重背景下的土地利用类型，预测榆林市土壤侵蚀潜力和变化情况，为该地区的水土保持和生态环境建设工作提供参考。

榆林市位于陕西省最北部，西邻甘肃环县、宁夏盐池县，北连内蒙古准格尔、伊金霍洛、乌审、鄂托克等4旗，东隔黄河与山西相望，南与陕西省延安市接壤。位于E 107°28′~111°15′、N 36°57′~39°35′。榆林市行政区划版图形似三角形。东西最大长度309 km，南北最大宽度295 km，总面积4万3 578 km²，约占陕西省的21%。地貌主要有风沙草滩区、黄土丘陵沟壑区、梁状低山丘陵区3大类。该区域气候属暖温带和温带半干旱大陆性季风气候，四季分明，无霜期短，年平均气温10 ℃，年平均降水量400 mm左右^[9]。

本研究使用的数据有：1)陕西省2000、2005、2008和2013年1∶10万土地利用数据, 空间分辨率30 m，来源于中国科学院资源环境科学数据中心(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年榆林榆阳区、绥德县、定边县、横山区、神木市等13个站点气象站。

目前RUSLE土壤侵蚀模型和RS、GIS已经深度融合。该模型在黄土高原地区不同尺度的土壤侵蚀研究中应用广泛，如羊圈沟小流域，面积约2.02 km^2[3]；榆林市，面积约4万3 578 km^2[10]，黄土高原整体，面积约64万km^2[5]。由于该地区细沟侵蚀比较严重，RUSLE估算值偏小，江忠善等^[11]将浅沟侵蚀因子加入到RUSLE模型中，提高估算精度^[10]。建立中国坡面水蚀预报模型：

式中：A为土壤流失量，t/(km²·a)；R为降雨侵蚀力因子，MJ·mm/(hm²·h)；L为坡长因子，量纲为1；S为坡度因子，量纲为1；K为土壤可蚀性因子，t·hm²·h/ (hm²·MJ·mm)；C为植被覆盖与管理因子，量纲为1；P为水土保持措施因子，量纲为1；G为浅沟侵蚀影响因子，量纲为1。

降雨侵蚀力R因子，依照谢云等^[12]和章文波等^[13]研究成果计算。LS因子是地形对土壤侵蚀的影响，用符素华等^[14]提供的基于DEM的LS计算工具计算。K因子利用几何平均粒径结合有机质模型计算^[15]。C因子利用蔡崇法等^[16]和林杰等^[17]的用NDVI推导的模型计算。我国2014—2020年完成退耕283万hm^2[9]，在2000—2017年NDVI数据基础上，用最小二乘法实现对未来NDVI的预测。对C因子分析采用RS时序分析和标准差分析方法, 用Hurst指数表示^[18]。浅沟侵蚀G因子的利用江忠善等^[11]的研究成果。GeoSOS-FLUS软件是根据FLUS模型原理开发的多类土地利用变化情景模拟软件。GeoSOS-FLUS利用神经网络算法^[19](artificial neuron net)获取各类用地的适宜性概率，再通过耦合系统动力学模型(system dynamics) 和元胞自动机(cellular automata) 模型以提高模型的适用性，在国内外土地利用变化模拟预测中得到了较好地应用^[20]。文中P因子的取值主要参照谢红霞^[21]的研究成果。文中所有因子栅格数据，重采样为30 m分辨率栅格数据再进行分析计算。

图 1为榆林市1988—2017年降雨侵蚀力因子分析，2015年最小为262.85 MJ·mm/(hm²·h)，2001年最大为2 394.59 MJ·mm/(hm²·h)，多年平均降雨侵蚀力均值为1 088.54 MJ·mm/(hm²·h)。图 2为1988—2017年榆林市多年平均降雨侵蚀力。丘陵沟壑区的佳县、米脂县、吴堡县和绥德县在1 150~1 350 MJ·mm/(hm²·h)之间，西部风沙区定边县、靖边县在800~1 050 MJ·mm/(hm²·h)之间，北部地区在1 000~1 150 MJ·mm/(hm²·h)之间。图 3为多年平均降雨侵蚀力标准差空间分布，标准差为682.44 MJ·mm/(hm²·h)，分布区间为(461.70，950.42)MJ·mm/(hm²·h)。1988—2000年R因子均值为998.67 MJ·mm/(hm²·h)，2001—2017年R因子均值为1 266.91 MJ·mm/(hm²·h)。21世纪的前13年和20世纪末的17年相比，降雨侵蚀力增加268.24 MJ·mm/(hm²·h)。

图 1(Fig. 1)

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

图 2(Fig. 2)

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

图 3(Fig. 3)

图 3 多年平均降雨侵蚀力标准差 Fig. 3 Standard deviation of multi-year average rainfall erosivity

LS因子分布在(0，80.29)之间，均值为3.98，榆阳区最低为1.57，清涧县最高为6.72。西部和北部地区比中东部地势平坦，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变为2013年的2.06；建设用地2000年LS均值由2.17变为2013年的1.91；未利用土地LS未发生变化。

图 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 km²的分辨率的GIMMS/NDVI数据在黄土高原地区，2000年NDVI和1984、1994年相近，都在0.250以下^[22]。2000年以后的MODIS NDVI，相比GIMMS/NDVI表现出更高的灵敏度，二者在黄土高原对植被的监测中表现出一致性。由MODIS NDVI计算的2000年C因子可以基本代表1984和1994年的C因子。图 5是利用最小二乘法预测的C因子空间分布，预测C因子为0.053，比2017年C因子的0.069，还有0.017的下降空间。

图 4(Fig. 4)

图 4 植被覆盖与管理因子 Fig. 4 Cover and managemt factor

图 5(Fig. 5)

图 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因子会增长。正向变化的和反向变化的区域分别达到51.53% 和48.74%。退耕还林(草)工程实施以后，丘陵沟壑植被覆盖增加显著，随之而来的是蒸发量增加，夏季蒸散量约占全年蒸散量一半^[23-24]，植树造林的人工林土壤含水量普遍低于农地，在降雨量>550 mm区域，土壤储水量损失最高。植被恢复越好，需水量越大，森林需水与土壤供水之间的矛盾突出，急需解决^[24]。未来极端干旱的年份，由于水资源的短缺，可能对植被生存产生威胁。

图 6(Fig. 6)

图 6 预测C因子的Hurst指数 Fig. 6 Predicting the Hurst index of C factor

结合榆林市“十四五”期间土地利用规划，重点考虑退耕还林和城市发展对土地利用的影响，把>25°坡耕地在土地利用类型转化成草地和灌木林地为高概率，对于平原地区的耕地在距离城市核心区>5 km的全部为禁止转换。利用神经网络网络模型在2013年土地类型基础上(图 7)预测未来到2022年土地利用类型(图 8)，2013年耕地面积为1万5 739.78 km²、林地面积为2 348.46 km²、草地面积为1万8 908.40 km²、水域面积为396.23 km²、建设用地面积为529.55 km²和未利用土地面积为4 390.12 km²。预测2022年耕地面积为1万4 468.05 km²、林地面积为3 036.88 km²、草地面积为1万9 487.55 km²、水域面积为396.10 km²、建设用地面积为662.50 km²和未利用土地面积为4 261.56 km²。

图 7(Fig. 7)

图 7 2013年土地利用类型 Fig. 7 Land use types in 2013

图 8(Fig. 8)

图 8 用神经网络网络模型预测土地利用类型 Fig. 8 Predicting land use types by artificial neuron net

利用预测的C、P因子和K、LS、K、R因子估算未来土壤侵蚀。将最终得到的各年土壤侵蚀强度划分为微度(0＜A≤500)、轻度(500＜A≤2 500)、中度(2 500＜A≤5 000)、强烈(5 000＜A≤8 000)、极强烈(8 000＜A≤15 000)和剧烈(A＞15 000)6级^[25]，单位为t/(km²·a)。2000—2017年土壤侵蚀模数如图 9a所示，相比年降雨侵蚀力下的土壤侵蚀计算，多年平均降雨侵蚀力下的土壤侵蚀估算实际意义更大，消除了年降雨侵蚀力差异引起的土壤侵蚀变化，评价更为客观。图 9b是多年平均降雨侵蚀力下土壤侵蚀模数序列。2000—2017年多年平均降雨侵蚀力因子R均值为1 220.95 MJ·mm/(hm²·h)。土壤侵蚀模数总体呈现波动下降的趋势。2000年土壤侵蚀模数最大，为3 559.99 t/(km²·a)；预测2022年最小，为793.27 t/(km²·a)。线性回归方程为y=-104.54x+3 103.1，R²=0.726 7，达到显著相关。2000—2017年多年平均降雨侵蚀力为1 101.02 MJ·mm/(hm²·h)。土壤侵蚀模数年际之间的差异，主要是植被覆盖与管理因子C和水土保持措施因子P的差异造成的。

图 9(Fig. 9)

图 9 2000—2017年土壤侵蚀模数 Fig. 9 Soil erosion modulus in 2000—2017

由图 9a可见，2000年的降雨侵蚀力最小，为439.54 MJ·mm/(hm²·h)；2001年的降雨侵蚀力最大，为2 394.59 MJ·mm/(hm²·h)。利用最大降雨、最小降雨和多年平均降雨侵蚀力条件下模拟2022年土壤各区县侵蚀情况(图 10)。最大降雨侵蚀力情况下的土壤侵蚀模数为1 752.68 t/(km²·a)；最小降雨侵蚀力情况下的土壤侵蚀模数为453.2 t/(km²·a)；多年平均降雨侵蚀力情况下的土壤侵蚀模数为793.27 t/(km²·a)。在这种预测模式下，除R因子之外，其余因子都是定值。3种不同情况下的图形线性分布呈现为增加或者减少趋势一致。在极大降雨侵蚀力情况下，定边县土壤侵蚀模数超过6 000 t/(km²·a)，子洲县超过3 000 t/(km²·a)，绥德县、米脂县和清涧县在2 000~3 000 t/(km²·a)之间。其余地区低于2 000 t/(km²·a)。在最小、最大和多年平均降雨侵蚀力下，榆林市总侵蚀量分别为1 900万t、7 370万t和3 370万t。

图 10(Fig. 10)

图 10 最小/最大和多年平均降雨侵蚀力下预测2022年的土壤侵蚀模数 Fig. 10 Predicted soil erosion moduli in 2022 under the min/max and multi-year average rainfall erosivity

在多年平均降雨侵蚀力条件下预测2022年各区县土壤侵蚀模数：神木市为306.32 t/(km²·a)、米脂县为738.37 t/(km²·a)、佳县为425.93 t/(km²·a)、绥德县为1 215.81 t/(km²·a)、吴堡县为793.23 t/(km²·a)、靖边县为1 115.91 t/(km²·a)、子洲县为1 556.88 t/(km²·a)、定边县为2 133.88 t/(km²·a)、府谷县为567.88 t/(km²·a)、榆阳区为240.09 t/(km²·a)、横山县407.22 t/(km²·a)、清涧县为1 599.75 t/(km²·a)。2022年榆林市大多数区县的土壤侵蚀模数将 < 2 000 t/(km²·a)，侵蚀以轻度和微度为主。

为更加直观地分析未来土壤侵蚀情况，绘制了多年平均降雨侵蚀力下各区县2000、2010、2016年和2022年各区县的土壤侵蚀模数(图 11)。2016年和2022年预测的土壤侵蚀模数相比较。神木市减少644.71 t/(km²·a)、米脂县减少174.86 t/(km²·a)、佳县减少812.1 t/(km²·a)、绥德县减少299.35 t/(km²·a)、吴堡县减少1 066.58 t/(km²·a)、靖边县减少659.18 t/(km²·a)、子洲县减少661.18 t/(km²·a)、府谷县减少510.36 t/(km²·a)、榆林市区减少469.7 t/(km²·a)、横山县减少592.76 t/(km²·a)、清涧县减少1 119.92 t/(km²·a)。定边县增加464.81 t/(km²·a)。到2022年榆林市大部分地区的土壤侵蚀量还有进一步下降的空间。

图 11(Fig. 11)

图 11 多年平均降雨侵蚀力下土壤侵蚀模数比较 Fig. 11 Comparison of soil erosion modulus under multi-year annual average rainfall erosivity

植被覆盖的大小，直接决定土壤侵蚀，NDVI越大，土壤侵蚀越小^[5]。40%~60% 植被覆盖度是植被对土壤侵蚀起明显干预作用的分界区间^[26]，笔者预测榆林市植被覆盖度可以达到62.14%，目前植被覆盖度为52.20%，NDVI指数还能增加约0.99，榆林市植被恢复已经接近潜力值。还林还草后水土保持措施对土壤侵蚀的减少起到主要作用，降雨对土壤侵蚀影响已处于次要的地位^[1]。大规模的植被生长使耗水增加，加剧了土壤干层，土壤净贮水量下降^[27]，植被恢复越好，需水量越大，Hurst指数表明未来，在丘陵沟壑区植被恢复较好的地区，森林需水与土壤供水之间的矛盾可能会更加突出，在极端干旱的年份，水资源的短缺，可能对植被生存产生威胁。模拟条件下，干化土壤中植被生长受当年降水量影响较大^[28]，植被迅速恢复过程中蒸发的增加是地表径流减少的主要原因之一^[29]，植被恢复和生长受控于降雨条件限制^[30]。下一阶段的还林还草应该“因水制宜”，避免加剧土壤干化。今后的还林还草工程要针对不同地区的降雨、土壤差异，合理选择草木品种。土壤侵蚀的进一步减少，还要进一步合理调配水资源的可持续利用，以保障还林还草工程取得生态和水土保持效益。

1) 1988—2017年多年平均降雨侵蚀力均值约为1 088.54 MJ·mm/(hm²·h)，丘陵沟壑区的降雨侵蚀力在1 150~1 350 MJ·mm/(hm²·h)之间，21世纪的前13年和20世纪末的后17年相比，降雨侵蚀力增加268.24 MJ·mm/(hm²·h)。

2) 还林还草工程实施18年来，植被覆盖持续改善，植被覆盖与管理因子C显著降低。2000—2017年多年平均降雨侵蚀力下，2000年的土壤侵蚀模数最大，为3 559.99 t/(km²·a)，2017年土壤侵蚀模数为1 369.19 t/(km²·a)，土壤侵蚀量减少958.65万t。预测2022年土壤侵蚀模数为793.27 t/(km²·a)，土壤侵蚀总体呈现波动下降的趋势。还林还草工程显著减少榆林市的水土流失。

3) 在最大、最小和和多年平均降雨侵蚀力条件下模拟计算未来土壤侵蚀，最大降雨侵蚀力条件下的土壤侵蚀模数为1 752.68 t/(km²·a)；最小降雨侵蚀力情况下的土壤侵蚀模数仅为453.20 t/(km²·a)。未来榆林市大多数区县的土壤侵蚀模数将 < 2 000 t/(km²·a)，土壤侵蚀以轻度和微度为主。土壤侵蚀模数约为793.27 t/(km²·a)。在最小、最大和多年平均降雨侵蚀力下，榆林市总侵蚀量分别为1 900万t、7 370万t和3 370万t。