生态系统向人类社会与环境提供的各种直接和间接的资源供给、环境调节、文化休憩和永续支持的功能，可被理解为生态系统服务功能^［1-2］。因此，在全球气候变化、生物多样性保护和可持续发展受环境制约日益严峻的背景下，开展生态系统服务功能评价，成为判断区域环境现状的主要途径。在干旱半干旱地区，随着气候变化和人类活动加剧，土壤侵蚀成为区域生态系统服务功能变化的主要表征，其发生发展变化过程直接影响到水土资源开发利用和社会经济发展^［3］。土壤侵蚀包括了水力、风力、冻融和重力侵蚀，其中水力侵蚀和风力侵蚀涉及范围最广，成为影响生态系统服务功能的主要过程。

在生态系统服务功能研究中，通过赋值量化生态系统服务功能，不仅可以揭示生态系统生态和经济价值，同时可为制定区域可持续发展政策提供基础数据支持^［3-5］。在生态系统服务功能研究方法方面，模型应用为主要研究手段。其中InVEST（Integrated Valuation of Ecosystem Services and Tradeoffs）模型是目前被广泛使用的评价工具^［4-6］。该模型内置许多模块并支持用户自行输入新建模块，能定量评估水源涵养、水土保持、防风固沙、碳固定和生物多样性维持等多种生态系统服务功能，输出结果具有二维可视性，可结合遥感和地理信息分析技术揭示生态系统服务功能的空间格局时空变化特征^［6］。在黄河和长江流域，InVEST模型典型研究案例解释了以流域（或行政区）为单元的生态系统服务功能时空变化以及驱动机制，研究人员借助模型不仅评估了流域的水源涵养功能，还揭示了土地利用和气候变化等因素对流域生态系统服务功能的影响^［7］，指出研究区生态系统服务功能变化受气候变化和人类活动共同驱动。生态系统服务功能表现出显著的时空差异，尤其是在水源涵养、土壤保持和碳固定等方面^{［6，8-9］}。位于黄河“几字弯”的内蒙古十大孔兑地区，地貌类型差异明显，自南向北呈带状分布着砒砂岩黄土区、库布齐沙漠和洪积冲积平原区^［10］。其中西柳沟流域是十大孔兑中土壤侵蚀最严重地区。1998年7月12日发生的洪水携带泥沙淤积堵塞黄河干流，形成一座长10 km、宽1.5 km、厚6.27 m的巨型沙坝，是历史上十大孔兑区域最大的水土流失灾害事件。同时，西柳沟流域黄土、沙漠和平原三类地貌面积比例在十大孔兑地区具有较好的代表性，而且也是近30多年来国家和地方生态建设的重点区域^［11］ 。

因此，本文选取1990年和2020年为研究时段，采用InVEST模型与遥感技术相结合的方法，对西柳沟流域生态系统服务功能从水源涵养、土壤保持、防风固沙和固碳4个方面分别进行单项和综合评估。主要目的是解释地貌分异对不同生态系统服务功能空间格局的影响，揭示近30年来生态系统服务功能变化特点，为十大孔兑区域今后的生态建设规划与布局提供基础。

西柳沟流域呈北低南高，海拔900~1 600 m，总面积为2 078.6 km²。地貌南部为砒砂岩黄土区、中部为库布齐沙漠区、北部为洪积冲积平原区，在十大孔兑中西柳沟流域在地貌类型分布方面具有最好的代表性（图1），三类地貌区面积分别为916.41 km²、596.88 km²和565.06 km²。黄土区以草地为主，面积为561.6 km²，占该类型区面积61.3%；沙漠区以流动沙地为主，面积为287.3 km²，占该类型区面积48.1%；平原区以农田为主，面积为332.0 km²，占该类型区面积58.8%。研究区植被类型以南部砒砂岩黄土区的典型草原植被、中部库布齐沙漠的沙生草本+灌丛植被和北部洪积冲积平原的人工乔木林+农业植被为显著特征。

十大孔兑流域边界矢量数据：源于内蒙古自治区“十三五”水利重大专项（NSK2017-M2）基础数据集。

DEM数据：源于地理空间数据云（https：//www.gscloud.cn/），采用空间分辨率为30 m的SRTM1 V3.0的中国区域数据。

土壤数据：来源于国家科技资源共享服务平台——国家地球系统科学数据中心-土壤分中心（http：//soil.geodata.cn）。

降水和蒸散发数据：源于国家地球系统科学数据中心（https：// www.geodata.cn/main/）。

NDVI数据：源于中国科学院地理科学与资源研究所土地利用与全球变化遥感团队发布的中国30 m年最大NDVI数据集。

土地利用数据：源于中国科学院资源环境科学数据中心基于Landsat 5和Landsat 8数据解译的全国土地利用数据集（分辨率 30 m）。

所有数据在输入InVEST模型前，在ArcGIS中采用重采样功能将分辨率都调整为30 m。

InVEST模型是用于对陆地、海洋和淡水生态系统服务功能进行评估的主要工具（更新至第7版），该模型内含的各类模块运行成熟，适用于不同生态系统、不同时空尺度的生态系统服务功能评估。其中，与陆地生态系统有关的计算模块主要有土壤侵蚀、碳储量和产水量等。InVEST模型通过输入较少的数据就可以获得较为可信的结果，打破了传统方法的局限性，使复杂问题简易化，运行结果可以在ArcGIS进行二次编辑，较纯文字的表达方式更加直观^{［3，5-6］}。

本文主要从水源涵养、水土保持、固碳和防风固沙4个指标进行单项评估。然后，依据各评估因子权重赋值对归一化后的单项指标进行综合评估。

式中：R是水源涵养量（mm）；P是年降水量（mm）；AET是实际蒸散量（mm），通过潜在蒸散量和土壤湿度条件估算；K是土地利用类型的植被涵养系数（无量纲），反映不同植被类型的水源涵养能力；C是土壤可渗透性系数（无量纲），代表土壤对水的渗透能力。其中， 植被涵养系数K依据研究区土地利用类型分别赋值为草地0.6、林地0.8、耕地0.4、沙地0.1。土壤可渗透性系数C依据研究区土壤类型分别赋值为沙土0.6、壤土0.4、黏土0.2。

式中：RKLS_i 为第i个栅格单元的年潜在土壤侵蚀量（t·hm^-2·a^-1）；USLE_i 为第i个栅格单元的年实际土壤侵蚀量（t·hm^-2·a^-1）；SDR_i 为第i个栅格单元的土壤保持量（t·hm^-2·a^-1）。年潜在土壤侵蚀量与年实际土壤侵蚀量之差为年土壤保持量（模型细节见文献［12］）。

式中：S是防风固沙量（t·hm^-2·a^-1）；V是植被覆盖度，表示植被对风沙侵蚀的阻止能力（百分比，无量纲）；W是风力指数，由风速数据计算；C_v是植被防风系数，反映不同植被类型在防护风沙中的效果，分别赋值为草地0.7、林地0.9、荒漠0.2^［13］；风力指数W通过风速的平方函数表示，是基于区域内和周边气象站风速数据采用公式W=v²计算获得（v是站点风速，1990、2020年的平均风速分别为5.6、5.8 m·s^-1）。植被覆盖度V通过NDVI数据计算并结合土地利用数据修正获得。

采用InVEST模型中的碳储存与封存模块计算固碳量：

式中：C为固碳量（t·hm^-2·a^-1）；V为植被覆盖面积（hm²）；Cv是植被固碳系数（t·hm^-2·a^-1）；S为土壤面积（hm²）；Cs为土壤固碳系数（t·hm^-2·a^-1）。各参数赋值为： Cv：草地0.5 t·hm^-2·a^-1，林地1.0 t·hm^-2·a^-1，耕地0.3 t·hm^-2·a^-1，沙漠0.1 t·hm^-2·a^-1；Cs：沙土0.4 t·hm^-2·a^-1，壤土0.6 t·hm^-2·a^-1，黏土0.8 t·hm^-2·a^{-1［14-15］}。

在单项服务功能评价的基础上，采用因子累加法进行综合服务功能的评价。不同生态系统服务功能权重用等值赋值法，即：水源涵养、水土保持、防风固沙和固碳能力的权重都为0.2。由于采用的4个评价因子单位不同，所以在累加评价之前对各因子进行归一化处理，获得无量纲结果，然后采用自然断点法做5级分类进行变化分析。

式中：X_j 为第j个单项服务功能归一化后的栅格值；X_i 为第j个单项服务功能归一化前的栅格值；X_minj 为第j个单项服务功能归一化前所有栅格中的最小值，X_maxj 为第j个单项服务功能归一化前所有栅格中最大值。综合生态服务功能大小由1级到5级逐渐增高，分别定义为1（低）级0.119~0.335、2（较低）级0.335~0.419、3（中等）级0.419~0.492、4（较高）级0.492~0.562和5（高）级0.562~0.762。综合评价的格局分析主要选取衡量不同级别斑块类型大小、密度、形状和空间关系的4个指标进行，计算在Fragstats中完成。

在InVEST模型中，水源涵养等4个模块的输出结果见图2，不同生态系统服务功能对应的空间属性统计结果见表1。在流域尺度上，1990年和2020年水源涵养量分别为3.78 mm和8.40 mm。从水源涵养量分布格局来看，不同年份都以研究区南部黄土区和北部平原区较高，而中部沙漠区较低为特征（图2），水源涵养量由大到小依次为黄土区、平原区和沙漠区。数量化评估细节表明，1990年黄土区和平原区水源涵养量分别高于全流域0.86 mm和0.52 mm，而沙漠区则低于全流域1.8 mm；2020年只有黄土区水源涵养量高于全流域2.36 mm，而平原区和沙漠区都低于全流域，分别低了0.31 mm和3.35 mm。从土壤保持量分布格局来看，研究区内不同年份均是北部平原区较高，而中部沙漠区和南部砒砂岩黄土区较低，土壤保持量依次为平原区、沙漠区和黄土区。数量化评估细节表明，1990年平原区和沙漠区土壤保持量分别高于全流域43.10 t·hm^-2·a^-1和21.44 t·hm^-2·a^-1，而黄土区则低于全流域40.54 t·hm^-2·a^-1；在2020年平原区和沙漠区土壤保持量分别高于全流域58.96 t·hm^-2·a^-1和18.44 t·hm^-2·a^-1，而黄土区则低于全流域48.36 t·hm^-2·a^-1。从防风固沙量分布格局来看，不同年份差异明显，1990年高值区主要分布在不同地貌类型交错区东侧，呈分离状态；2020年已经连片分布，而且沙漠区防风固沙量明显提高；1990年不同地貌类型区防风固沙量相差不大，为2.02~2.27 t·hm^-2·a^-1，而2020年防风固沙量明显提高，为2.86~3.72 t·hm^-2·a^-1。数量化评估细节表明，1990年黄土区防风固沙量高于全流域0.11 t·hm^-2·a^-1，而平原区和沙漠区则分别低于全流域0.04 t·hm^-2·a^-1和0.14 t·hm^-2·a^-1；2020年沙漠区和黄土区防风固沙量高于全流域0.26 t·hm^-2·a^-1和0.19 t·hm^-2·a^-1，而平原区则低于全流域0.60 t·hm^-2·a^-1。从固碳量分布格局来看，不同年份特征一致，都是中部沙漠区和北部平原区较高，南部黄土区偏低，不同地貌类型区固碳量在不同年份相差不大，1990年为6.06~6.73 t·hm^-2·a^-1，2020年为6.28~6.66 t·hm^-2·a^-1。数量化评估细节表明，1990年黄土区固碳量低于全流域0.34 t·hm^-2·a^-1，而平原区和沙漠区则分别高于全流域0.21 t·hm^-2·a^-1和0.33 t·hm^-2·a^-1；2020年黄土区固碳量仍低于全流域0.20 t·hm^-2·a^-1，而平原区和沙漠区则分别高于全流域0.14 t·hm^-2·a^-1和0.18 t·hm^-2·a^-1。

1990年西柳沟生态系统服务功能较高级别（4级）分布面积最大，为891.3 km²，占研究区面积44.8%（表2）；较低（2级）和中等（3级）级别分布面积相差不大分别为454.5 km²和393.4 km²，占比分别为22.8%和19.8%；高（5级）级别分布面积为233.6 km²，占比为11.7%；低（1级）级别分布面积为16.4 km²，占比不足1%。2020年西柳沟生态系统服务功能仍以较高级别分布面积最大，与1990年相比增加了213.6 km²，达到了1104.9 km²，占研究区面积比例提高到了55.5%；分布面积居第二位是高级别，分布面积较1990年增加了138.1 km²，达到了371.7 km²，占比也提高到了18.7%；其余3个级别分布面积较1990年都呈减少状态，中等、较低和低三个级别分别减少了109.1、235.0、7.5 km²。

从分级格局来看（图3），1990年较低以下级别集中分布在中部沙漠区，较高级别主要在平原区，而高级别主要集中在平原与沙漠交错区和黄土区最南部。到2020年，较低以下级别面积明显减少（表2，图3），主要分布在沙漠区西部，其他部分主要转换为中等和较高级别。高级别区域在1990年基础上呈扩大状态，其中黄土区东北部变化最明显。

从1990年综合评价分级格局数量化特点来看（表3），斑块破碎化方面以较低级别最高，斑块密度达到了3.01 个·km^-2；较高、中等和低级别破碎化程度差别不大，斑块密度变化在1.27~1.89 个·km^-2；高级别斑块破碎化最低，斑块密度为0.78 个·km^-2。从图3可以明显看出黄土区是各种级别斑块破碎化最为严重区，沙漠区次之，平原区最轻。从最大斑块指数来看，1990年较低和较高级别拥有连片分布的大斑块，最大斑块指数占流域面积的16.0%和14.1%；中等和高级别最大斑块指数分别为7.0%和5.9%，连片分布的最大斑块面积不到上述二类的一半；低级别斑块最大斑块指数相对甚微。从图3中可以明显看出，高和较高级别最大斑块分布在沙漠和黄土交错区，低级别最大斑块分布在沙漠区，中等级别分布在黄土区。从不同斑块级别平均分维数来看，由于本研究采用的分辨率为30 m，不同级别拥有大量独立栅格（斑块）分维数为1，使斑块分维数在不同级别相差无几，为1.039~1.046。从同级别斑块之间的距离来看，1990年低级别类型斑块平均距离最大，为163.5 m，说明该级别斑块之间的离散性最大，而高级别类型斑块间平均距离最小，为98.3 m，说明该级别斑块间聚集度较强，呈集中分布。

从2020年综合评价分级格局数量化特点来看（表3），斑块破碎化仍以较低级别最高，斑块密度为2.80个·km^-2；中等和较高级别破碎化程度相近，斑块密度分别为1.57 个·km^-2和1.50 个·km^-2；高级别斑块密度为0.92 个·km^-2，低级别为0.35 个·km^-2。从最大斑块指数来看，较高级别拥有连片分布的大斑块，最大斑块指数占流域面积的29.1%；高和中等级别最大斑块指数相差不大，分别为8.5%和7.7%，最大斑块面积不到较高级别的1/3；与1990年相比，较低级别最大斑块指数显著减少，为4.8%，低级别最大斑块指数保持稳定，仍为0.05%。不同斑块级别平均分维数变化幅度与1990年相似，也为1.039~1.046，略有不同的是除低级别类型分维数保持不变外，其他级别各有微量增减。从不同级别斑块间距离来看，2020年仍以低级别类型斑块平均距离最大，为222.0 m，保持最大的斑块间离散性，而高和中等级别类型斑块间平均距离相差不大，分别为137.1 m和134.0 m；与1990年相比较，较低和较高级别斑块间距都为减少状态，斑块间距分别为122.0 m和110.2 m，分别减少了0.4 m和4.9 m。

十大孔兑生态系统服务功能在不同地貌区空间差异比较明显，这些差异不仅表现在不同服务功能单项评价中，而且在时间尺度上也具有较显著的特点^{［6-7，10］}。其中，水源涵养量以黄土区最高、平原区次之和沙漠区最低为特点，而且这一特点在此次两个研究时段中具有一致性。从变化特点来看，1990—2020年全流域和不同地貌分区水源涵养量都呈增加趋势，全流域水源涵养量增加了4.62 mm，是1990年的2.2倍；各地貌分区以沙漠区增加最明显，增加了3.07 mm，是1990年的2.55倍。这与西柳沟流域在过去30多年中因植被恢复与生态建设工程极大改善了地表植被覆盖状态密切相关，而植被因子在InVEST模型水源涵养功能模块计算中具有较高的赋值。

土壤保持量以平原区最高、沙漠区次之和黄土区最低为特点，这一特点同样在此次两个研究时段中具有一致性。从变化特点来看，1990—2020年全流域和不同地貌分区土壤保持量都呈减少趋势，全流域土壤保持量减少了22.62 t·hm^-2·a^-1，是1990年的99.2%；各地貌分区减少不明显，其中黄土区减少最多，为30.44 t·hm^-2·a^-1，是1990年的98.9%。产生这一结果与当年降水量和降水分布格局有关，虽然从1990—2020年，研究区的植被状况和水土保持与荒漠化治理得到逐步改善与加强，但是土壤保持量计算中最大影响变量是降雨侵蚀力^［8］。而离西柳沟最近的达拉特旗气象站观测数据表明2020年的降雨量（359.1 mm）高于1990年（347.0 mm），差别不大，造成2020年土壤保持量略有下降可能与年内大于12 mm降水量的差异有关。防风固沙量在不同时段存在差异，1990年以黄土区最高，2020年以沙漠区最高。从变化特点来看，1990年到2020年全流域和不同地貌分区防风固沙量都呈增加趋势，全流域防风固沙量增加了1.30 t·hm^-2·a^-1，是1990年的1.6倍；各地貌分区以沙漠区增加最明显，增加了1.70 t·hm^-2·a^-1，是1990年的1.8倍。西柳沟流域防风固沙功能改善在沙漠区最大，这与该区域近30多年坚持荒漠化治理与利用密切相关，多年累积的“三北”防护林工程、环京津风沙源治理工程的成果，是该区域防风固沙功能改善的根本原因^［10-11］。固碳量以沙漠区最高、平原区次之和黄土区最低为特点，这一特点在本次研究的两个时段中也具有一致性（表1）。从变化特点来看，1990—2020年除沙漠区为弱减少趋势（减少了0.07 t·hm^-2·a^-1）外，全流域和其他地貌分区固碳量都呈弱增加趋势，全流域仅增加了0.08 t·hm^-2·a^-1；各地貌分区以黄土区增加相对明显，增加了0.22 t·hm^-2·a^-1，是全域增量的2.75倍。固碳量变化不大，与2020年遥感数据选择受限、无法获取最佳（云量为0）生长旺期数据有关。

考虑到现状生态系统服务功能具有多维度的特点，采用综合分级格局评价对生态系统服务功能具有更好的概括性和客观性。在1990年，综合评价分布面积最大的较高级别，斑块密度、最大斑块指数和斑块间平均距离都小于分布面积居次位的较低级别。这与西柳沟流域地貌分异明显有关，较低级别在沙漠区集中连片分布拉高了最大斑块指数，而在黄土区零散分布导致了斑块密度增加，在平原区的离散分布又导致了斑块间距加大。较高级别最大斑块分布在平原区，而黄土区同样也是该级别零散分布密集区。所以，从地貌分区差异对生态系统服务功能影响来看，黄土区最复杂导致分级格局空间异质性最高；沙漠区和平原区空间异质性相对较低。在2020年，综合评价分布面积最大仍是较高级别（表2），斑块密度和最大斑块指数都大于分布面积居次位的高级别，而斑块间距相反。这种结果表现在地貌分异上呈现出以下两个特点，一是在黄土区和平原区较高级别呈连片分布，二是高级别主要分布在平原-沙漠交错区和沙漠-黄土交错区。综合评价分级格局的变化（1990—2020年）主要发生在黄土区和沙漠区。其中，高和较高级别分布面积增加，分别增加了138.1 km²和213.6 km²，低、较低和中等级别面积减少，分别减少了7.5、235.0、109.1 km²。导致变化的原因是近30年的砒砂岩黄土区水土保持工程实施、草原科学禁牧和库布齐沙漠化治理与开发利用并举等生态工程建设^［12-15］。

西柳沟地貌类型空间分异格局对生态系统服务功能有明显影响。1990—2020年，水源涵养量增加了4.62 mm，其中沙漠区增加幅度最大，达到2.55倍；土壤保持量减少了22.62 t·hm^-2·a^-1，其中黄土区减少最多，为30.44 t·hm^-2·a^-1；防风固沙量增加了1.30 t·hm^-2·a^-1，其中沙漠区增加最大，为1.7 t·hm^-2·a^-1；固碳量增加了0.08 t·hm^-2·a^-1，其中黄土区增加最大，为0.22 t·hm^-2·a^-1。

在时间尺度上，生态服务功能综合评价较高级的部分都是分布面积最大类型。1990年较低级别在沙漠区集中连片分布，2020年较高级别在黄土区和平原区呈连片分布，都是最大斑块指数居高的根本原因。综合评价分级格局变化主要发生在黄土区和沙漠区，在研究时段内高和较高级别面积分别增加了138.1 t·hm^-2·a^-1和213.6 t·hm^-2·a^-1，低、较低和中等级别面积分别减少了7.5、235.0、109.1 t·hm^-2·a^-1。