沙尘是影响中国北方干旱和半干旱地区的重要自然现象^［1-3］。为有效改善京津冀地区生态环境状态、减轻风沙危害、提升居民生活质量，2002年，国家紧急启动京津风沙源治理工程，目前治理范围已由第一阶段（2002—2012年）的45.8万km²扩大至第二阶段（2013—2022年）的70.6万km²，累计投资超过1 290亿元^［4-6］。项目启动后，区域通过封禁、飞播、退耕还草/林等多种手段进行生态环境的恢复和保护，当前京津风沙源治理工程区已经成为中国北方生态屏障的重要组成部分。因此，量化工程前后区域治理效果对于客观认识工程效益和治理措施的科学调控具有重要意义。

土地利用/覆盖变化是人类活动与自然环境相互作用最直接的表现形式，也是量化生态建设最有效的途径。其包括土地利用/覆盖类型的转换（如农田开垦、森林砍伐和城市扩展），也包括类型未变条件下土地覆盖特征的微变（如农业种植结构、退化草地封育）^［7］。无论哪类变化，均会带来相应的生物多样性、初级生产力、土壤质量、风蚀程度等生态系统功能和服务的改变^［8］。已有研究多关注土地利用/覆盖类型变化（例如，退耕还林）带来的生态环境效益，例如，盛晓雯等^［9］研究发现京津冀地区土地利用类型的变化对城市生态系统服务价值影响较小。然而，对于土地覆盖渐变（例如，禁牧育草）带来的效益，关注还不够。

生态系统服务价值是评估生态工程效益的有效途径。Costanza等^［10］率先采用价值系数法评估了全球生态系统服务价值，开启了生态系统服务价值评估的先河。随后，一大批学者针对中国自然生态系统特征开展了价值系数的修正，对中国森林、草地等生态系统开展了广泛的评估，掀起了生态系统服务价值评估的研究热潮^［11-15］。针对价值系数法的不足，一些研究尝试对评估方法进行改进，谢高地等^［16］基于单位面积当量因子运用时空调节因子法（初级生产力、降水、土壤保持）对其评估方法进行改进；郝林华等^［17］按照基准年价值水平运用其不同年份价格指数比值对区域调节服务的3项指标价值量进行修正。这些方法在客观评估生态系统服务价值方面为我们提供了丰富的参考。然而，对于如何反映土地覆盖渐变过程中生态系统服务价值的变化，仍然缺乏有效的处理方法。

京津风沙源治理工程区土地利用/覆盖变化及其生态环境效应一直是研究热点^{［5，18-22］}。吴丹等^［5］、迟文峰等^［20］基于区域土地利用/覆盖类型数据揭示了工程实施前后耕地、林地、草地和未利用地等类型之间的转换特征。Wu等^［23］、Yang等^［24］基于归一化植被指数（Normalized difference vegetation index，NDVI）数据检测了工程实施以来区域植被的动态变化。李屹峰等^［25］基于土地利用数据运用InVEST模型评估了密云水库流域的生态环境效益，江凌等^［26］基于土地利用数据采用空间统计分析评估了内蒙古土地利用变化对生态环境效益的影响。

总结现有研究，仍然存在以下两方面的不足：①现有研究均不能明确定量地回答不同土地覆盖类型上渐变的强度和时空格局；②土地利用/覆盖类型变化和渐变两种过程分别对生态系统服务价值变化有多大的贡献仍然不清楚。这些问题的解决对于有针对性地开展工程区治理措施的评估和土地的可持续管理是十分必要的。

因此，本研究基于土地利用/覆盖和NDVI等数据，从土地利用/覆盖类型转换和土地覆盖渐变两个角度揭示京津风沙源治理工程区1990—2018年土地利用/覆盖变化特征，并针对土地覆盖渐变程度修正了价值系数，评估了同期生态系统服务价值的变化。

京津风沙源治理工程区二期（2013—2022年）西起内蒙古的乌拉特后旗，东至内蒙古的阿鲁科尔沁旗，南起陕西省的定边县，北至内蒙古的东乌珠穆沁旗（36°49′—46°40′N，105°12′—120°59′E），东西横跨近1 300 km，南北纵跨近1 175 km^［19］。区域包括北京、天津、河北、山西、内蒙古和陕西等6个省（区、市）的138个县（旗、市、区），总国土面积为70.6万km^2［22］。研究区地貌复杂多变，由山地、平原、高原三大地貌类型组成。区域地处中纬度内陆地区，从东到西，年平均温度和年降水量均呈下降趋势，分别为-2—13 ℃和250—470 mm^{［19，22］}。受气候和地貌条件影响，区域西北部主要发育典型草原和荒漠草原，土壤以栗钙土为主，而东南部植被主要为温带落叶阔叶林和温带落叶灌丛，土壤以褐土和棕壤土为主^［27］。基于区域地貌、植被覆盖等自然条件，可以划分为11个亚区^{［4，19，22，27］}（图1）。

本研究使用的数据包括土地利用/覆盖、NDVI和气象数据等。1990、2000、2010、2018年1∶100 000的土地利用/覆盖数据来源于中国科学院全国土地利用/覆盖数据集，数据综合精度达91.2%^［28-29］。数据包括耕地、林地、草地、水域、建设用地和未利用地等6个一级类和25个二级类，为了阐明区域沙源地的变化特征，我们将未利用地中的沙地单独提出，构建了7个一级类。NDVI数据来源于美国国家航空航天局网站（https：//ecocast.arc.nasa.gov）提供的1990—2018年时间序列数据，其中1990年1月至2000年12月的数据来自全球库存监测与建模系统（Global Inventory Monitoring and Modeling System，GIMMS），空间分辨率为8 km，时间分辨率为15 d^［30］；2001年1月至2018年12月的数据为中分辨率成像光谱仪产品（Moderate-resolution Imaging Spectroradiometer，MODIS）空间分辨率为250 m，时间分辨率为16 d^［24］。将NDVI数据采用最大值合成方法进行处理和合成，得到年度NDVI数据^［31］。为了消除传感器和波距之间的绝对差异，利用像元二分模型进一步计算植被覆盖度^［32-33］。气象数据（年降水量、年均气温）来源于中国国家气象信息中心（http：//data.cma.cn/），包括来自分布在研究区内及周边的66个气象站，剔除了无关因素的干扰，时间序列为1990—2018年。

为保证计算中数据的匹配性，所有数据分辨率重采样为1 km×1 km。

采用地图代数方法分析土地利用/覆盖类型之间的转换^［34］。

利用该公式可以得出由k时期到k+1时期的土地利用变化图C_i×j，它表现了土地利用变化的类型及其空间分布。基于此，可获取1990—2000年、2000—2010年、2010—2018年3个阶段的转换图，并统计了土地利用转移矩阵表。

在1990—2018年土地利用/覆盖未变化区域，采用线性趋势法统计分析各土地利用/覆盖类型上植被覆盖度的变化特征，以揭示土地利用/覆盖的渐变过程^{［18，35］}。

式中：x为时间；y为对应年份某未变化地类上的覆盖度平均值；a为变化斜率；b为常数，x¯和y¯为x和y的平均。a值的正/负表示增加/减小趋势，在趋势显著的条件下，a的大小表示植被覆盖度变化的幅度。

我们提取了1990—2018年地类未发生变化的林地、草地、耕地和沙地区域，并计算区域1990—2018年逐年覆盖度平均值，进而获取各地类上土地利用/覆盖渐变特征。

根据土地利用/覆盖类型将区域划分为农田、森林、草地、水域和沙地5种生态系统，参照“千年生态系统评估”（Millennium Ecosystem Assessment， MA）的分类体系，将生态系统服务划分为供给服务、调节服务、支持服务和文化服务，具体包括食物生产、原材料生产、气体调节、气候调节、水文调节、废物处理、保持土壤、维持生物多样性和提供美学景观等9项^{［10，36］}。

采用当量因子法计算区域生态系统服务价值^{［10，16］}。

式中：ESV为生态系统服务总价值（元）；Ak为研究区k种生态系统的面积（hm²）；VCk为生态系统服务总价值系数（元·hm^-2·a^-1）；ESVf为单项生态系统服务价值（元）；VCkf为单项生态系统服务系数（元·hm^-2·a^-1）。

生态系统服务价值是量化和评价生态服务强弱的重要方法^［10］，是生态系统服务的货币表现，并不能够进行直接的实地验证。相关研究表明，生态系统服务价值与生态系统的生物量有密切关系，生物量越大，生态服务价值越高^{［10，14，37］}，同时，草地植被覆盖度是度量生物量的有效指标^［38］。因此，为了能够有效体现草地质量或盖度对生态系统服务价值的影响，我们在谢高地等^［16］提出的生态服务价值动态当量的基础上，根据草地和沙地的植被盖度对当量进行了修订。具体修订步骤为，首先，对于草地，按其覆盖度将其分为3类，包括高覆盖度草地（>50%）、中覆盖度草地（20%—50%）和低覆盖度草地（<20%）^［39］；对于沙地，由于植被盖度阈值20%对于局地土壤侵蚀和生态系统功能具有重要意义^［40］，因此，沙地按植被覆盖情况分为两类，包括高流动性沙地（>20%）、低流动性沙地（<20%）。然后，结合各类土地覆盖上植被平均盖度与全区域平均盖度的比值，来计算系数修正因子。草地和沙地价值当量因子的计算方法如下，

式中：Fnj为草地或沙地生态系统第j个亚类的生态系统服务n（包括食物生产、原材料生产、维持生物多样性等9项服务功能）的单位面积价值当量因子；Fn指谢高地等^［41］确定的该类生态系统生态服务n的当量因子；Gj指草地或沙地第j个亚类的修正系数；Bj为亚类j区域范围内的平均植被覆盖度；B¯为整个研究区对应生态系统类型的草地或沙地平均植被覆盖度。

在改进的修正当量的基础上，我们计算了区域各个生态系统单位面积生态系统服务当量（表1）。1个标准当量因子的生态系统服务价值量（元·hm^-2）参考谢高地等^［16］方法，为3 406.5元·hm^-2。基于此，我们计算了京津风沙源治理工程区1990—2018年区域生态系统服务价值的变化。

研究在全区和亚区两个尺度上开展，空间分析在ArcGIS10.6软件中完成，统计分析在SPSS26.0中完成。

京津风沙源治理工程区土地利用/覆盖以草地为主，总体上呈现草地和耕地集中分布，林地、沙地零散镶嵌的格局。2018年，区域草地面积达39.15万km²， 占区域总面积的55.25%，主要分布于典型草原东区（面积为8.15万km²，占草地总面积的20.82%）、鄂尔多斯高原（6.88万km²，17.58%）、浑善达克沙地（5.82万km²，14.87%）、典型草原中区（5.82万km²，14.86%）；耕地面积为12.45万km²，占区域总面积的17.57%，主要分布于农牧交错带区（面积为3.13万km²，占耕地总面积的25.14%）、燕山丘陵山地水源保护区（2.18万km²，17.47%）、鄂尔多斯高原（1.53万km²，12.26%）、晋北山地丘陵区（1.45万km²，11.61%）；林地、沙地面积相对较小，分别为7.64万、4.36万km²，但林地集中分布在燕山丘陵山地水源保护区和大兴安岭南部区；沙地集中分布在鄂尔多斯高原和浑善达克沙地（图2）。

1990—2018年，区域土地利用/覆盖类型的转换随时间段和亚区位置而具有明显的异质性（图3）。1990—2000年，土地利用/覆盖变化特征主要表现为耕地和沙地的增加、草地的减少。具体来看，有5 022.75 km²的草地被开垦为耕地，主要分布在科尔沁沙地（面积为1 703.58 km²，占变化耕地总面积的33.92%）、农牧交错带区（1 416.62 km²，28.20%）、大兴安岭南部区（760.23 km²，15.14%）；有2 841.04 km²的草地退为沙地，主要分布在鄂尔多斯高原（面积为1 394.00 km²，占变化沙地总面积的49.07%）、浑善达克沙地（815.45 km²，28.70%）、典型草原中区（272.55 km²，9.59%）；有2 776.43 km²的沙地恢复为草地，主要分布在鄂尔多斯高原（面积为1 827.49 km²，占变化草地总面积的65.82%）、科尔沁沙地（303.30 km²，10.92%）、浑善达克沙地（245.12 km²，8.83%）。

2000—2010年，主要表现为沙地、耕地和林地的增加，草地的减少。具体来看，有2 423.82 km²的草地退为沙地，主要分布在鄂尔多斯高原（面积为1 126.80 km²，占变化沙地总面积的46.49%）、浑善达克沙地（550.23 km²，22.70%）、典型草原中区（303.26 km²，12.51%）；有2 300.29 km²草地被开垦为耕地，主要分布在科尔沁沙地（面积为523.59 km²，占变化耕地总面积的22.76%）、农牧交错带区（471.51 km²，20.50%）、大兴安岭南部区（344.51 km²，14.98%）；有1 709.31、997.01 km²耕地分别恢复为草地和林地，其中，草地主要分布在农牧交错带区（面积为554.04 km²，占变化草地总面积的32.41%）、鄂尔多斯高原（379.28 km²，22.19%），林地主要分布在晋北山地丘陵区（面积为305.89 km²，占变化林地总面积的30.68%）、鄂尔多斯高原（265.94 km²，26.67%）、农牧交错带区（237.52 km²，23.82%）。

2010—2018年，土地利用/覆盖类型的转换强度有所减弱，总体表现为耕地、沙地的减少，草地和林地的增加。其中，有1 122.65 km²耕地变为建设用地，主要分布在农牧交错带区（面积为389.75 km²，占变化建设用地总面积的34.72%）、燕山丘陵山地水源保护区（244.59 km²，21.79%）；有816.45 km²沙地恢复为草地，主要分布在鄂尔多斯高原（面积为578.14 km²，占变化草地总面积的70.81%）、典型草原中区（101.10 km²，12.38%）、浑善达克沙地（68.78 km²，8.42%）；有230.89 km²沙地转变为林地，主要分布在鄂尔多斯高原（面积为215.29 km²，占变化林地总面积的93.24%）和典型草原西区（14.89 km²，6.45%）。

1990—2018年，京津风沙源治理工程区植被盖度表现为总体改善、局部退化的特征，在不同地类上表现也不同（表2）。区域沙地植被总体呈改善趋势，覆盖度每10年约增长4.22%。其中，鄂尔多斯高原的沙地植被改善尤为显著，覆盖度每10年增长达7.49%，而在农牧交错带区、典型草原东区、大兴安岭南部区和科尔沁沙地上植被盖度呈显著下降趋势，下降幅度分别为每10年14.78%、11.94%、11.45%和3.13%。

区域草地、林地和耕地植被总体呈波动变化，但在部分亚区体现出显著的趋势。从草地来看，黄河灌溉区、晋北山地丘陵区和鄂尔多斯高原等亚区植被改善趋势显著，覆盖度每10年分别增长11.34%、10.89%、9.34%；从林地来看，鄂尔多斯高原、晋北山地丘陵区和典型草原西区改善趋势显著，覆盖度每10年分别增长13.24%、7.91%、5.18%，而在大兴安岭南部区的呈下降趋势，下降幅度为每10年8.42%；从耕地来看，鄂尔多斯高原、晋北山地丘陵区、黄河灌溉区和典型草原西区等亚区改善趋势显著，覆盖度每10年分别增长17.77%、9.72%、8.73%、6.58%，而在大兴安岭南部区呈下降趋势，下降幅度为每10年6.66%。

区域草地对生态系统服务价值的贡献最大。2018年，京津风沙源治理工程区生态系统服务总价值为34 973.66亿元，其中，草地生态系统服务总价值为22 458.71亿元，占2018年区域总价值的64.22%；其次是森林和农田生态系统，其服务价值分别为7 315.39、3 349.75亿元，占比分别为20.92%、9.58%（表3）。

1990—2018年，区域生态系统服务价值总体呈增加趋势，各个时段和不同地类之间存在差异。1990—2000年，区域生态系统服务价值减少了1 221.14亿元，是1990年区域生态系统服务总价值的3.90%；高覆盖度草地生态系统服务价值减少量最大，为2 804.84亿元，占1990年区域生态系统服务总价值的8.96%；其次为高流动性沙地，减少了21.67亿元，占1990年区域生态系统服务总价值的0.07%；与此同时，其他生态系统服务价值均呈增加状态，其中，中覆盖度草地增加量最大，为1 339.20亿元，占1990年区域生态系统服务总价值的4.28%；其次是低覆盖度草地，为118.91亿元，占1990年区域生态系统服务总价值的0.38%。

2000—2010年，区域生态系统服务价值减少了599.22亿元，占2000年区域生态系统服务总价值的1.99%。中覆盖度草地生态系统服务价值减少量最大，为1 552.81亿元，占比为5.16%；其次为水域和高覆盖度草地，其服务价值共计减少了293.58亿元，占2000年区域生态系统服务总价值的0.98%；而森林、低覆盖度草地、高流动性沙地和农田生态系统服务价值呈增加状态，其中，森林生态系统服务价值的增加量最大，为830.57亿元，占2000年区域生态系统服务总价值的2.76%；其次是低覆盖度草地和高流动性沙地，其服务价值共计增加了416.12亿元，占2000年区域生态系统服务总价值的1.38%。

2010—2018年，区域生态系统服务价值增加了5 475.57亿元，占2010年区域生态系统服务总价值的18.56%。其中，高覆盖度草地生态系统服务价值增加量最大，为14 897.91亿元，占2010年区域生态系统服务总价值的50.5%；其次为森林、水域和低流动性沙地，其服务价值共计增加了158.56亿元，占2010年区域生态系统服务总价值的0.54%；其余各项生态系统服务价值均呈减少状态，其中，中覆盖度草地减少量最大，为8 636.01亿元，占2010年区域生态系统服务总价值的29.28%；其次是低覆盖度草地和高流动性沙地，其服务价值共计减少了919.66亿元，占2010年区域生态系统服务总价值的3.12%。

为了体现土地覆盖渐变对生态系统服务的影响，本研究对生态系统服务价值评估当量进行了修订，对比系数修订前后区域生态系统服务价值，我们发现，土地覆盖渐变的影响是不可忽视的。依据谢高地等^［41］提供的系数（原系数）计算，1990—2018年，草地、沙地生态系统服务价值减少了703.08亿、17.36亿元，而修订后的结果表明同期草地生态系统服务价值增加了2 790.83亿元、沙地生态系统服务价值减少了50.52亿元（表4）。主要原因是原系数对于区域土地覆盖不区分渐变，均为统一标准，而系数修订后，土地覆盖渐变导致的价值变化得以体现。具体来看，1990—2018年，草地总面积减少了2.24万km²；若从内部覆盖度亚类来看，高覆盖度草地面积增加了18.97万km²，中覆盖度草地和低覆盖度草地分别减少了22.26万、2.13万km²。同期，沙地总面积减少了0.7万km²；若从内部亚类来看，高流动性沙地面积减少了1.0万km²，而低流动性沙地增加了0.31万km²。因此，造成了两种方法价值评估的差异，而系数修订前，忽略了这种地类内部渐变的情况。

土地利用/覆盖类型变化和渐变均会使得生态系统服务发生变化，进而影响生态系统服务价值。土地利用类型的变化影响生态系统的能量交换、物质循环、土壤侵蚀与堆积等主要生态过程，从而改变生态系统服务的提供^{［36，42］}。例如，森林的调节与支持服务较高，而食物供给服务相对较低；农田的农产品供给服务较高，而调节、支持与文化服务相对较低^［43］。本研究结果表明，1990—2018年，生态系统服务价值总体增加了3 655.21亿元，其中，由类型变化导致生态系统服务价值增加了914.90亿元，主要由草地恢复为森林贡献；同时，土地覆盖渐变导致生态系统服务价值增加了2 740.31亿元，主要由高覆盖度草地面积增加贡献（表5）。总体来看，土地覆盖渐变导致的生态系统服务价值变化量更大。

各个时间段也表现出类似的特征。1990—2000年，生态系统服务价值总体减少了1 221.14亿元，其中，由土地利用类型变化导致生态系统服务价值增加了146.73亿元；同时，土地覆盖渐变导致的生态系统服务价值减少了1 367.87亿元，主要是高覆盖度草地退化使其服务价值减少了2 804.84亿元，而中覆盖度和低覆盖度草地服务价值分别增加了1 339.20亿、118.91亿元。2000—2010年，生态系统服务价值总体减少了599.22亿元，其中，由土地利用类型变化导致生态系统服务价值增加了647.64亿元；同时，土地覆盖渐变导致的生态系统服务价值减少了1 246.86亿元。2010—2018年，生态系统服务价值总体增加了5 475.57亿元，其中，由土地利用类型变化导致生态系统服务价值增加了120.53亿元；而土地覆盖渐变导致生态系统服务价值增加了5 355.04亿元（表5）。

土地利用/覆盖变化是自然和人类活动共同作用的结果。近30年京津风沙源治理工程区年均气温呈波动上升趋势，年降水量呈波动趋势^［44］。为探讨区域植被覆盖变化与气候因素的关系，我们采用区域降水和温度数据，计算了1990—2018年植被覆盖度与年降水量和年气温的相关性。结果显示，区域植被覆盖度受温度变化影响不大；而区域降水对沙地和草地植被覆盖度有积极影响（图4），体现在沙地植被覆盖度与年降水的相关系数为0.41（P<0.05），区域草地植被盖度与年降水的相关系数为0.56（P<0.01）。

我们进一步采用残差分析法厘定了自然和人类活动对土地覆盖渐变影响的主导特征。结果发现，人类活动对区域沙地和草地植被盖度较显著；2003年后，区域沙地、草地植被盖度残差值均为正值，表明人类活动对沙地、草地植被盖度变化均产生正面影响（图5）。这与国家和地方采取的一系列措施实现京津风沙源治理工程区的生态恢复密不可分（表6）。例如，区域在生态工程实施前（1990—2000年）是“沙进人退”的现象，实施中（2000—2010年）如内蒙古克什克腾旗封沙育林、农防林、荒山荒地造林等措施营造出工程实施后的“人沙和谐”局面^{［18，45］}；生态工程的实施促使鄂尔多斯高原和晋北山地丘陵区耕地、林地和草地的覆盖度呈上升趋势，黄河灌溉区草地、耕地覆盖度呈上升趋势^［46］。可见，人类活动在一定程度上改善了区域沙地、草地植被覆盖状态，成为区域土地利用/覆盖变化的主导驱动力。

1990—2018年，研究区域土地利用/覆盖类型变化随时间段具有异质性。2010年之前，耕地、草地、沙地面积减少，林地面积增加；而2010年后，沙地、耕地面积减少，草地和林地面积增加。

区域土地覆盖渐变总体改善、局部退化，不同地类存在差异。其中，沙地植被总体改善，覆盖度每10年约增长4.22%；而草地、林地和耕地植被总体呈波动变化，但在部分亚区增加显著（如，鄂尔多斯高原、黄河灌溉区等）。

1990—2018年，区域生态系统服务价值呈现先减少后增加、总体增加的特征。过去近30年，总体增加了3 655.21亿元，其中，土地覆盖渐变导致服务价值增加了2 740.31亿元，大于土地利用/覆盖类型变化的贡献。

我们建议京津风沙源治理工程区有针对性地开展工程治理措施的评估和土地资源的可持续管理，不仅要重视退耕还林、未利用地造林种草，更要关注草地和林地的修复和恢复，注重生态工程成果的维持和质量的提升。