黄土高原是中国"两屏三带"生态安全格局的重要组成部分^[1-2], 该地区沟壑纵横, 水土流失严重。水土流失过程会破坏原有土壤结构, 降低土地生产力, 造成人民财产损失, 同时在流域范围内淤积河道, 造成面源污染, 加重洪涝和干旱灾害等潜在威胁, 已成为全球普遍关注的生态环境问题之一^[3-4]。为控制黄土高原水土流失, 该地区自20世纪70年代起开始实施一系列水土保持治理措施^[5]。退耕还林(草)工程(1999)实施后, 黄土高原地区植被覆盖度大幅增加, 有效降低了该地区水土流失的危害。因此, 研究该地区植被恢复与土壤侵蚀之间的相互关系引起众多学者的关注。植被覆盖对土壤侵蚀的影响研究主要集中于植被的防蚀功效、减水减沙效应。其中, 植被冠层消减了雨滴动能, 地表覆被物分散了径流动能, 根系提高了土壤抵抗径流侵蚀能力^[6], 并且不同的植被类型能在不同程度上减轻土壤侵蚀的危害。植被恢复具有延长坡面产流时间、增加入渗、减少产流产沙的作用^[7]。程复^[8]研究指出退耕还林(草)工程实施后, 2008年黄土丘陵沟壑区所有流域单元的输沙模数整体平均降低54%, 土壤侵蚀程度明显减弱；赵巧巧^[9]研究发现黄土高原地区流域年均归一化植被指数(Normalized Difference Vegetation Index, NDVI)和年均输沙量呈现出负相关的关系；王怡菲^[10]研究表明当退耕还林(草)面积增加1%, 流域含沙量每立方米减少1.894kg；陈浩^[11]研究提出以植被恢复为主的退耕还林(草)工程对黄土高原地区北洛河流域土壤侵蚀减少的贡献率为85%。另一方面, 部分学者^[12-14]研究土壤侵蚀对植被的影响后, 发现土壤侵蚀过程会降低土壤的持水能力和养分积累, 影响植被群落的发育和演替, 并通过侵蚀过程的机械力对植物产生干扰和破坏。

退耕还林(草)工程距今已实施二十余年, 植被恢复能够抑制土壤侵蚀, 改善生态环境, 促进生态系统服务的提升, 同时, 生态系统服务功能的提升有利于植被恢复。由此可见, 植被恢复的生态环境效应及其带来的一系列人类福祉是维持区域可持续发展的关键因素。目前植被与土壤侵蚀的关系方面已有大量研究, 但关于黄土高原植被格局及其对土壤侵蚀的影响和植被恢复带来的生态效益尚未明确。此外, 生态系统服务的提升也是研究的热点。鉴于此, 本文以黄土高原延河流域为研究对象, 探讨2000—2015年该区的植被恢复状况, 并借助SWAT(Soil and Water Assessment Tool)模型估算土壤侵蚀量和土壤保持量, 揭示植被格局对土壤侵蚀的影响, 定量研究植被恢复带来的生态系统土壤保持服务价值增益, 对今后植被建设和水土流失治理工作以及生态补偿具有一定的理论和实际指导意义。

1 研究区概况与数据来源 1.1 研究区概况

延河流域(36°27'—37°58'N, 108°41'—110°29'E)作为黄土高原丘陵沟壑区的典型流域, 地表破碎, 沟间地以丘陵为主, 梁、峁状丘陵大约占流域沟间地的80%, 为流水侵蚀和重力侵蚀提供了"有利"的地貌条件。流域处于东部季风湿润区与内陆干旱区的过渡地带, 年降水量偏少, 多年平均降水量为520mm, 降水季节分配不均, 集中于夏季, 秋季次之, 春季较少, 冬季有少量降雪。该流域水土流失严重, 河流含沙量大, 泥沙多为悬移质, 水土流失面积接近流域总面积的80%, 是黄河中游水土流失最严重的区域之一。延河流域的植被特性也具有过渡性, 随着植被类型、覆盖度等生态特点的不同, 植被抗蚀作用也有所差异。因此本文选取延河流域作为研究区域, 其位置分布如图 1所示。

1.2 数据来源与处理

本文使用的数据包括NDVI数据及构建SWAT模型所需的数字高程模型(Digital Elevation Model, DEM)、土地利用数据、土壤数据、水文数据和气象数据。(1)NDVI数据采用美国地质勘探局(http://glovis.usgs.gov/)提供的MODIS植被指数产品—MOD13Q1, 研究选用的时间范围为2000—2015年, 每年23期影像, 共368期, 时空分辨率为16d及250m。利用MRT(MODIS Reprojection Tools)软件对该数据产品进行格式转换和投影转换, 将HDF格式转换成TIFF格式, 并将SIN投影转换成WGS_1984_UTM_Zone_49N投影。同时, 利用MVC(Maximum Value Composites)合成月NDVI最大值, 来进一步消除云、大气和太阳高度角等因素对遥感影像产生的影响^[15]。由于温带半干旱区域存在明显的季节性, 生长季与非生长季植被覆盖存在较大差异^[16], 因此实验采用延河流域生长季5—10月NDVI的平均值表征植被覆盖的年际变化状况^[17]；(2)DEM数据来自中国科学数据云地理空间数据云(http://www.gscloud.cn/), 空间分辨率为30m；(3)土地利用数据来源于中国科学院资源科学环境数据中心(http://www.resdc.cn/), 空间分辨率为30m, 本文选用2005年和2010年土地利用数据, 主要包括八类土地利用类型：耕地、林地、园地、高覆盖草地、低覆盖草地、建设用地、水域和未利用地, 按照SWAT模型自带的土地利用类型编码方式建立相对应的属性索引表, 构建土地利用数据库；(4)土壤数据采用寒区旱区科学数据中心(http://westdc.westgis.ac.cn/)提供的世界土壤数据库(Harmonized World Soil Database, HWSD), 其中的中国土壤数据是第二次全国土地调查南京土壤所提供的1 ∶ 100万土壤数据, 利用SPAW(Soil-Plant-Atmosphere-Water)软件计算所需土壤参数, 建立土壤数据库；(5)水文数据来自《中华人民共和国水文统计年鉴》的黄河流域资料, 包括甘谷驿水文站2003—2015年逐日径流量和产沙量；(6)气象数据来自中国气象数据共享网(https://data.cma.cn/), 包括延河流域内及周边站点的位置信息和1999—2015年逐日降水量、最高和最低温度、相对湿度、风速数据。由于目前国家基本气象站点无太阳辐射数据, 所以该数据由天气发生器模拟提供。为了减小空间数据带来的偏差, 本文将所有的空间数据投影变换为统一的坐标系(WGS_1984_UTM_Zone_49N), 且栅格分辨率重采样为30m。

2 研究方法 2.1 NDVI变化特征

NDVI可从整体上表征植被的生长状况和覆盖变化^[18-19], 被众多学者广泛应用^[20-22]。本文运用趋势线分析法逐像元分析2000—2015年延河流域的NDVI变化趋势, 反映植被的时间演变和空间差异特征。依据趋势线的斜率判断变量在研究时段的变化趋势, 若斜率大于0, 则植被覆盖呈增加趋势；若斜率等于0, 则植被生长趋于稳定；若斜率小于0, 则植被覆盖呈减少趋势。斜率绝对值的大小表示研究时段植被覆盖增加或减少的速率。计算公式如下：

(1)

式中, θ为趋势线的斜率；n为研究时段年数；x_j为研究时段内第j年年际或生长季NDVI值。

2.2 SWAT模型

美国农业部农业研究中心开发了具有很强物理机制、开源的SWAT分布式水文模型, 其可以综合利用土壤类型、土地利用类型及管理措施预测和模拟流域长时期内的径流、泥沙、化学元素和有机物质的迁移运动, 评估整个流域范围内的水分平衡^[23]。该模型广泛应用于黄土高原地区的产流产沙模拟^[24-25], 因此, 本文采用SWAT模型中的产沙模块, 即改进的通用水土流失方程(Modified Universal Soil Loss Equation, MUSLE)模拟延河流域的土壤侵蚀量, 包括实际土壤侵蚀量和潜在土壤侵蚀量。实际土壤侵蚀量指流域实际上发生的土壤侵蚀量, 潜在土壤侵蚀量指假设流域没有地表植被覆盖与管理措施下发生的土壤侵蚀量, 二者之间的差值为土壤保持量。计算公式如下：

(2)

式中, S_sj为模拟的实际土壤侵蚀量(t)；S_qz为模拟的潜在土壤侵蚀量(t)；ΔS为研究区生态系统土壤保持量(t)；Q_surf为地表径流(mm/hm²)；q_peak为径流峰值(m³/s)；area_hru为水文响应单元(Hydrologic Research Unit, HRU)的面积(hm²)；K_usle为土壤可蚀性因子；C_usle为地表植被覆盖与管理因子, 当无植被覆盖与管理因子的影响时, 设C_usle=1；P_usle为水土保持措施因子；LS_usle为地形因子, 包括坡长因子L和坡度因子S；CFRG为直径大于2mm的粗碎块因子。

采用德克萨斯农业大学提供的SWAT-Calibration and Uncertainty Programs(SWAT-CUP)软件中的Sequential Uncertainty Fitting(SUFI-2)方法对SWAT模型进行敏感性参数分析、率定和验证。使用决定系数(Coefficient of determination, R²)和纳什效率系数(Nash-Sutcliffe efficiency coefficient, NSE)作为评价SWAT模型模拟精度的指标。R²可评价模拟值与实测值变化趋势的一致性, NSE可反映模拟值与实测值之间的拟合程度。一般来说, 当R²>0.6, 且NSE>0.5时, 即可认为SWAT模型模拟取得显著性效果^[26]。评价指标的表达式如下：

(3)

(4)

式中, Q_mi为模拟值；Q_oi为实测值；Q_o为实测均值；Q_m为模拟均值。

通过基础数据的准备, 构建数据库, 借助ArcSWAT2012版本软件建立延河流域的SWAT模型。首先是子流域的划分：基于DEM数据生成河网, 通过设定阈值(形成子流域的最小给水面积)为10000hm²重新划分河网, 并设置流域总出水口, 流域划分为47个子流域；其次划分HRU：通过叠加土地利用数据、土壤数据和坡度数据, 获得各个子流域中具有相同土地利用类型、土壤类型和坡度的1987个HRUs；最后将气象数据等导入模型, 重新读取基础数据库, 运行SWAT模型, 实现对延河流域1999—2015年径流量和产沙量的模拟。以1999—2002年作为预热期, 再分别以2003—2008年和2009—2015年作为率定期和验证期, 并依据甘谷驿水文站实测数据进行参数率定和验证。由于延河流域的地形、土壤等基础数据长时间内几乎不变, 而流域内的土地利用状况不断变化。为提高SWAT模型的精度, 本文将1999—2008年(预热期+率定期)作为整体, 运用2005年的土地利用数据驱动SWAT模型；2009—2015年(验证期)选用2010年的土地利用数据驱动SWAT模型。经验证, 本文构建的延河流域SWAT模型率定和验证结果均达到模型要求(如图 2), 可用于研究区的水文模拟。

2.3 生态系统土壤保持服务价值评估

采用市场价值法、机会成本法、恢复费用法和影子工程法等生态系统服务价值核算方法, 依据每个HRU上的生态系统土壤保持量, 从保持土壤肥力、减少土地废弃和减轻泥沙淤积(3个指标)方面定量估算土壤保持服务价值, 最后汇总计算流域总价值量, 阐明植被恢复的生态效益。

(5)

(6)

(7)

(8)

式中, E为生态系统土壤保持服务价值(元)；E₁为保持土壤肥力价值(元)；E₂为减少土地废弃价值(元)；E₃为减轻泥沙淤积价值(元)；ΔS为研究区生态系统土壤保持量(t)；B_i为土壤中氮、磷、钾及有机质的平均含量, 分别为0.17%、0.06%、1.40%、0.90%^[27]；C_i为土壤中的氮、磷、钾及有机质折算成相应肥料(尿素、过磷酸钙和氯化钾)及碳比率的系数, 分别为2.164、4.065、1.923、0.5^[28]；D_i为化肥的市场价格(元/t), 依据全国肥料的平均售价, 分别为1825、522、1948、550元/t；D₁为单位面积土地的机会成本, 取值为2.35×10⁵元km^-2 a^-1[27]；p为单位土壤体积质量分数, 取值为1.20t/m³；h为土壤厚度, 江忠善、郑粉莉等^[29-30]提出的黄土高原地区平均土壤厚度为0.2m；依据中国主要河流泥沙运动规律, 土壤流失泥沙淤积在河道、湖泊、水库中的系数取0.24^[31]；D₂为建设水库工程的费用, 取5.714元/m^3[27]。

3 结果与分析 3.1 NDVI时空变化特征 3.1.1 年际NDVI时空变化

从图 3可以看出, 2000—2015年NDVI均值表现出递增的趋势, 增速达到0.0245/10a。其中, NDVI平均值由2000年的0.514增加到2014年的0.553, 其增长率达到7.05%, 主要是由于退耕还林(草)工程的实施, 在政策因素和人为因素的影响下, 延河流域大量耕地实现还林还草, 植被覆盖状况逐渐改善。为更好评价延河流域植被恢复状况, 分析该流域2000—2015年NDVI时空变化趋势。经统计相关系数显著性检验, 参考岳本江等^[32]对延河流域NDVI变化趋势分级标准理论并结合本文研究对象实际情况, 将该流域的NDVI变化情况划分为明显退化(θ≤－0.007)、退化(－0.007＜θ≤－0.001)、基本不变(－0.001＜θ≤0.001)、改善(0.001＜θ＜0.004)和明显改善(θ≥0.004)5个等级。

由图 4可知, 绿色区域的NDVI呈现出增加的趋势, 其中84%的区域改善, 9.6%的区域明显改善, 主要位于流域上游和下游。由于流域上游地势陡峭, 坡度大于25°的陡坡占70%以上^[6], 而退耕还林(草)工程对于坡度大于25°的坡耕地实现了全部还林, 导致植被覆盖变化大, 植被恢复好；而下游地区地势平缓, 土壤肥沃, 农业活动频繁, 土地缺乏植被保护, 对该地区平缓的坡耕地进行退耕还草, 植被得到保护, 使得植被覆盖逐渐增大。还有6.1%的区域NDVI变化趋势基本不变, 主要位于流域中游, 这是因为延安市以南地区梢林广布, 植被生长状况好造成的结果。流域仅有0.3%的区域NDVI出现减少的趋势, 主要分布于延安市新区, 由于城市化使得植被退化显著。

3.1.2 生长季NDVI时空变化

生长季NDVI时间变化分析。本文把一年中植被开始返青到落叶的时期作为生长季^[33-34], 根据延河流域的物候特征, 其生长季从5月初至10月底。由于MODIS-NDVI的影像为半月合成产品, 把NDVI均值分为上、下半月, 可以在时间上更细致的表达植被覆盖的变化情况。如图 5所示, 以时间为横轴(刻度为半月), 以NDVI均值为纵轴, 用回归法生成拟合曲线, 即NDVI变化趋势线, 其斜率反映变化率。经统计2000—2015年生长季各月(5—10月), 时间与NDVI的相关系数大于0.6(P＜0.01), 2000—2015年生长季NDVI变化呈现出逐年增加的总趋势, 但不同月份的增加程度不同。其中每年8月的NDVI值最大, 每年6月的NDVI增长速率最快(拟合曲线斜率最大), 主要原因是流域气候条件和植被物候的季节差异^[35]。而一年中NDVI随物候变化, 5—8月NDVI值逐渐增大, 9—10月NDVI值逐渐减少。

生长季NDVI空间变化分析。如图 6所示, 2000—2015年各月NDVI变化的空间差异显著, 主要由于植被在一年内经历开始返青(生长季始期)到落叶(生长季末期)的过程, 植被进入生长季始期和生长季末期的面积比例均为先增加后减少^[33-34]。依据上文NDVI变化趋势分级标准, 其中生长季始期(5月)NDVI增加的区域主要分布在流域中游, 该地区位于延安市城区附近, 受人类活动影响较大, 植被生长较上游和下游更好；6月NDVI增加的区域主要位于下游, 由于下游地区退耕还林(草)工程实施效果显著, 使得下游植被状况逐渐改善；7月NDVI增加的区域占97.8%, 全流域植被长势较好；8月NDVI增加的区域主要位于上游和下游, 中游植被生长状况趋于稳定, 这与图 5中8月NDVI值最大形成呼应, 全流域植被生长状况最好；9月全流域NDVI增加趋缓, 植被进入生长季末期；10月NDVI变化趋势呈增加的区域主要位于上游, 退耕还林(草)工程实施中流域上游的坡耕地还林的效果显著。

总体上看, 该流域NDVI时空变化特征为2000—2015年生长季各月NDVI值及NDVI增长速率不同, 表现在空间上, 流域上中下游植被变化呈现出不同的趋势, 这与植被类型、地势条件、气候条件、人类活动等因素相关。

3.2 植被覆盖对土壤保持服务的作用

为探究植被覆盖措施对土壤侵蚀的影响, 本文借助SWAT水文模型, 以植被覆盖措施作为模型的变量, 以HRU为研究单元, 定量模拟土壤侵蚀量, 以此说明植被覆盖变化对流域土壤保持的积极作用。从图 7可以看出, 在HRUs尺度上, 2000—2015年年均实际和潜在土壤侵蚀模数的拟合曲线整体上呈减少趋势, 其中多年平均实际和潜在的土壤侵蚀模数为31783t km^-2 a^-1和624015t km^-2 a^-1, 主要是因为流域的土壤类型多为黄土母质上发育的黄绵土, 土壤质地均一, 土质疏松, 抗侵蚀能力差^[36], 在没有植被覆盖措施下, 土壤极易发生侵蚀, 且侵蚀模数大。由此定性说明植被覆盖对土壤侵蚀具有强大的抑制作用。

如何定量探究植被覆盖发挥的效果是一个值得研究的问题。基于此, 本文以2000年土壤保持量为基准, 计算2001—2014年逐年土壤保持增量, 从定量的角度深究植被恢复带来土壤保持量的变化特征。从表 1可知, 2001—2014年土壤保持量均存在不同程度的增加。植被覆盖抑制土壤侵蚀, 土壤得到有效保持, 土壤保持又促进植被生长, 从而形成良性循环, 不断改善生态环境。因此, 在今后的生态恢复中要注重植被覆盖措施。

植被覆盖措施是减少土壤侵蚀的重要因素之一。从不同植被覆盖度的角度分析土壤保持量的变化是极具价值的。参照宋敏敏等^[37]在延河流域已有的研究经验, 采用自然间断点分级法, 将NDVI分为低覆盖(0—0.3)、中低覆盖(0.3—0.45)、中覆盖(0.45—0.6)、中高覆盖(0.6—0.75)、高覆盖(0.75—1)5个等级。经统计, 2000—2015年延河流域的植被覆盖情况为：低覆盖占0.26%, 中低覆盖占15.8%, 中覆盖占66.43%, 中高覆盖占17.05%, 高覆盖占0.46%, 说明延河流域植被主要以中覆盖为主。为更明确土壤保持量的变化, 标准化多年(2000—2015)平均土壤保持量, 参照邓伟等^[38]的研究, 采用等间隔分级法, 将标准化后的土壤保持量划分为低(0—0.2)、较低(0.2—0.4)、中(0.4—0.6)、较高(0.6—0.8)、高(0.8—1)5个等级。延河流域各等级土壤保持量如下：低保持占75.77%, 较低保持占13.64%, 中保持占6.72%, 较高保持占3.17, 高保持占0.7%, 表明该流域土壤保持量以低保持为主。

运用叠加分析法, 研究植被覆盖对土壤保持量的影响。从图 8可以看出, 流域植被覆盖变化引起土壤保持量的变化特征明显。低覆盖下土壤保持量几乎为0, 由于植被覆盖度小于0.3, 地表景观为裸岩、裸土和稀疏植被^[39], 土壤侵蚀严重；中低覆盖下, 存在低、较低、中和较高土壤保持, 几乎没有高保持, 反映出中低植被覆盖条件下, 生态系统土壤保持功能仍有限；而中覆盖条件下, 高保持面积占比最大, 为85.35%, 反映出生态系统土壤保持功能强, 主要原因是对先前坡度大于25°、农业活动频繁和矿区等地区还草还林力度大, 植被恢复显著、覆盖度高, 所以土壤侵蚀量大幅减少；流域高覆盖区面积占比最小, 仅不到1%, 因此土壤保持功能强但总量小。因此, 今后水土保持工作要注意对各个等级植被覆盖采取不同的措施：低覆盖区重点治理, 中低、中和中高覆盖区预防保护, 着重扩大高覆盖区。

3.3 生态系统土壤保持服务价值时空分布 3.3.1 时间尺度变化

土壤保持服务作为一项重要的生态系统调节服务, 在防止土壤侵蚀, 减少泥沙等方面发挥重要的作用^[40-41]。为定量评估植被恢复所带来的生态系统土壤保持服务价值, 本文通过公式(5—8)计算土壤保持服务价值并分析其时空分布特征, 探究植被恢复对土壤保持服务的增益效果。

从图 9可以看出, 时间尺度上, HRUs年均土壤保持服务价值随时间波动, 主要原因是2000—2015年的降水量和降水强度不同, 导致土壤侵蚀程度不同。其中2013年土壤保持服务价值达到多年最大值, 因为2013年延河流域遭遇有气象记录以来的强降水^[42], 该年降水量约为790mm, 是2000—2015年降水量的最大值。总体看, 生长季HRUs月均土壤保持服务价值在年内呈单峰型变化, 峰值在7月, 主要是因为流域多年以来年内降水量也呈单峰型变化, 7月降水量约为111mm, 达到年内最大值, 且以阵雨或暴雨为主, 局部地区出现暴雨的情况较多, 暴雨强度也大。本文土壤保持服务价值随时间的变化特征与张乐涛等^[43]研究产沙量变化情况相吻合。黄土高原地区几乎所有的土壤侵蚀量都发生在雨季, 且随着降水量的增加土壤保持量增加, 由此表明土壤保持量在时间尺度上受降水因素的影响, 产生的生态系统土壤保持服务价值也随之波动。在多雨季节, 适宜的温度和降水促进了植被的生长, 当NDVI达到峰值时, 产沙量迅速下降, 土壤保持量快速增加。因此, 植被覆盖是减少黄土高原地区水土流失的主要因素。

2000—2015年HRUs年均土壤保持服务价值为3.64×10⁸元, 生长季月均土壤保持服务价值为6.06×10⁷元。从表 2和表 3可以看出, 年际和年内保持土壤肥力价值E₁、减少土地废弃价值E₂和减轻泥沙淤积价值E₃对土壤保持服务价值的贡献不同, 其中E₁占96.7%, E₂占1.5%, E₃占1.8%, 表明2000—2015年延河流域生态系统土壤保持服务价值以保持土壤肥力价值为主。

总体上看, 2000—2015年延河流域土壤保持服务价值呈波动增长的态势, 生态系统土壤保持服务增益明显, 从2000年仅为2.63×10⁶元, 至2014年已增长为5.48×10⁸元, 翻了一番。

3.3.2 空间尺度变化

本文重点分析退耕还林(草)工程实施以后, 延河流域生态系统土壤保持服务价值的空间差异, 希望能为该地区今后退耕还林(草)工程实施的方向提供参考依据。从图 10可以看出, 2000年和2014年流域的土壤保持服务价值空间差异大。2000年流域的土壤保持服务价值分布多集中在上游地区, 主要是因为延河流域地势西北高、东南低, 西北部的坡度大于25°的坡耕地多, 退耕还林(草)工程实施初期对该地区的还林力度大；2014年土壤保持服务价值在流域内趋于均匀分布的态势, 因为退耕还林(草)工程实施以来, 全流域植被恢复显著, 全面提升了生态系统土壤保持服务, 进而带来显著生态效益。流域年均土壤保持服务价值在下游地区最高, 上游次之, 中游最低, 其原因为：下游残塬平梁沟壑区地势平坦, 土壤肥沃, 农业活动频繁, 先前土地缺乏植被保护导致水土流失严重, 通过退耕还草, 土壤保持量大幅提升；上游低山地区植被生长环境差, 通过退耕还林、园林绿化、水保工程措施等, 使得土壤侵蚀减弱, 相较于中游, 土壤保持服务价值更高；中游梁峁丘陵沟壑区侵蚀严重, 并且延安市宝塔区人口较多, 人类活动频繁, 城市建设导致不透水面不断增加, 因此土壤保持服务价值低。今后水保工作要注意对流域上中下游采取不同措施：上游限制土地开发, 以生态恢复为重；中游生态与建设并重, 加强土地规划与保护；下游因地制宜, 合理开展高附加值农业活动提高土地生产效率, 主要种植果树、牧草等植被覆盖度较高的作物, 尝试免耕法等保护性耕作方式。

4 结论

本文以黄土高原延河流域为研究对象, 运用MODIS-NDVI数据探究2000—2015年研究区植被恢复的时空变化特征, 同时基于SWAT模型揭示植被恢复对抑制土壤侵蚀、提高土壤保持的作用, 并定量评价生态系统土壤保持服务价值增益效果。主要得出以下结论：

(1) 通过分析2000—2015年延河流域年际和生长季NDVI的时空变化特征, 得出该流域植被覆盖整体上处于增加态势, 退耕还林(草)工程实施效果显著, 持续向好。但是上中下游存在空间差异, 流域上游植被覆盖变化大、植被恢复好；下游地区植被覆盖逐渐转好；而流域中游, 出现城市化快速增长区植被退化问题。因此, 在黄土高原生态脆弱区, 城市化发展要立足于生态可持续, 节约高效利用土地资源, 合理规划建设用地, 为生态建设留足余地。

(2) 通过定量分析, 探明植被恢复对抑制土壤侵蚀的作用十分显著。从退耕还林(草)工程实施以来, 随着植被覆盖度增加, 在HRUs尺度上年均实际土壤侵蚀模数处于逐年递减的态势, 说明植被恢复措施能够有效抑制土壤侵蚀, 促进生态系统土壤保持量增益明显。而总体上, 延河流域植被主要以中覆盖为主, 流域土壤保持量以低保持为主。通过运用GIS叠加分析法, 研究植被覆盖对土壤保持量的影响, 表明植被低覆盖条件下生态系统土壤保持量几乎为0；植被中低覆盖条件下, 生态系统土壤保持功能仍有限；而中覆盖条件下, 土壤保持功能较强；植被高覆盖区所占面积最小, 仅不到1%, 因此植被高覆盖区的生态系统土壤保持功能虽然强, 但是对流域土壤保持总量的贡献很小。综上所述, 今后的水土保持工作应该注重能够有效提高植被覆盖度的生态恢复措施, 要注意对低覆盖区重点治理, 中低、中和中高覆盖区预防保护, 着重扩大高覆盖区。

(3) 实施退耕还林(草)工程带来了巨大的生态系统土壤保持服务价值增益, 从2000年至2014年翻了一番, 流域年均(2000—2015)生态系统土壤保持服务价值为3.64×10⁸元，生长季月均土壤保持服务价值为6.06×10⁷元，而且土壤保持服务价值以保持土壤肥力价值为主。流域生态系统土壤保持服务价值时空分布特征也有规律可循。时间尺度上, 土壤保持服务价值受降水因素影响而波动变化；月均土壤保持服务价值在年内一般呈单峰型变化, 峰值在7月, 与降水的单峰型变化相一致。空间尺度上, 上中下游土壤保持服务价值的空间差异大, 今后防治水土流失工作要注意对流域上中下游因地制宜, 采取不同措施：上游改善生态, 中游加强保护, 下游科学耕作。

5 讨论

(1) 当前土壤侵蚀的定量化模拟存在一定困难。土壤侵蚀是影响植被发育并受植被反作用的一种生态应力^[44-45]。大多数研究^[46-49]土壤侵蚀采用的模型一般是美国通用水土流失方程(Universal Soil Loss Equation, USLE), 但该模型是经验模型缺乏土壤侵蚀机理, 在实际应用中不能详细描述土壤侵蚀的物理过程。为此基于物理过程的水蚀预报模型应运而生, 例如美国农业部研发了WEPP(Water Erosion Prediction Project)模型, 该模型几乎涉及与土壤侵蚀相关的所有过程^[50-51], 然而这类模型存在参数众多、所需数据难以获取等缺点, 推广应用很困难。因此, 本文选择目前最流行的SWAT分布式水文模型。SWAT模型采用改进的通用水土流失方程MUSLE模拟基于HRU的土壤侵蚀量, MUSLE中地表径流、径流峰值、土壤可蚀性等因子的计算方法较为科学, 使模型具有了一定的物理机制。

(2) 本文基于两期土地利用数据构建SWAT模型, 并将HRU作为研究单元细化了水文模拟过程。为提高模型模拟精度, 本文在构建SWAT模型时基于土地利用数据、土壤数据和坡度数据将子流域划分成HRUs的离散集合, 把HRU作为研究单元, 每一个分布式地理参数(如土壤导水率、土壤可蚀性等)都可以针对每个HRU进行定义, 从而细化了水文过程。另一方面, 目前多数研究在借助SWAT模型模拟长时间序列的水文状况时, 将一期土地利用数据作为输入数据源。由于土地利用数据代表流域下垫面的状况, 而下垫面对模拟水文过程具有重要的意义, 尤其是退耕还林(草)工程实施后, 土地利用变化大, 长时间序列的研究则更需要考虑土地利用数据的真实性。因此, 本文采用2005年和2010年两期土地利用数据作为输入数据源, 有利于提高模型模拟的准确性。

(3) 植被动态对土壤侵蚀的影响是一个不可忽视的问题, 尚需要进一步研究。本文定量研究植被恢复的生态系统服务价值增益问题, 但植被恢复是一个长期动态的过程, 如何实现基于植被动态模拟水文过程尚需要进一步研究。MODIS遥感数据产品MOD13Q1具有免费性、连续性、长时间序列等优势, 而SWAT模型具有开源性的特点, 为植被动态变化的生态水文效应研究提供了必要条件。将NDVI数据引入SWAT模型实现植被动态变化条件下的水文过程模拟, 是今后研究的重难点。

(4) 本文的研究结论认为退耕还林(草)工程因地制宜, 将不良耕地转换为林地或草地, 实现植被恢复显著, 使生态系统土壤保持服务价值增益显著。这与众多学者研究提出的实施退耕还林(草)工程以来黄土高原水土流失减少, 生态环境改善, 生态系统服务价值增强的观点一致。

本文存在的不足之处：(1)假设研究时段地形未发生变化的问题。本文将DEM数据和土壤数据作为长期不变的数据源输入SWAT模型, 但实际上地形一直处于微小的变化, 这种变化宏观遥感数据难以探测发现, 对实验结果可靠性存在一定影响, 因此微地形变化也是一个重要的研究问题；(2)假设研究时段土壤未发生变化的问题。随着退耕还林(草)工程的实施, 植被不断恢复的同时, 伴随着土壤肥力不断恢复, 但本文中假设土壤性质也不变, 将土壤数据作为常量, 使得研究结果具有一定偏差。因此, 在今后研究中应该进一步考虑微地形、土壤变化对土壤侵蚀的影响, 使研究更加科学合理。

致谢: 感谢西安科技大学测绘科学与技术学院周自翔副教授对本研究的帮助。