基于InVEST模型的门头沟区生态系统土壤保持功能研究

刘晓娜, 裴 厦, 陈 龙, 刘春兰

(北京市环境保护科学研究院 国家城市环境污染控制工程技术研究中心, 北京 100037)

生态系统服务是生态学和地理学的研究前沿和热点，是人类赖以生存和发展的资源与环境基础，是将自然过程与人文过程联系起来的桥梁和纽带[1-3]。土地利用与土地覆被变化(LUCC)通过改变生态系统类型、格局以及生态过程直接影响生态系统服务功能，是生态系统服务变化的重要驱动力[4-5]，定量评估LUCC与生态系统服务变化之间的关系成为当前生态学研究的热点[6-8]。

土壤保持功能是指森林、草地等地表植被具有防止与减少土壤侵蚀的功能，是生态系统服务功能的重要组分，对于土壤形成、涵养水源、防土固沙以及减少水土流失具有重要意义[1]。目前对土壤保持功能的研究主要是采用通用土壤流失方程(USLE)[9-10]，但USLE模型并未考虑地块自身拦截上游沉积物的能力，通过USLE模型计算的土壤保持量存在一定的问题[11]。随着"3S"技术的快速发展，基于USLE和GIS的模型方法也相应的产生[12-14]，生态系统服务价值化和权衡得失综合评价模型(InVEST)，在USLE模型基础上考虑了地块自身拦截上游沉积物能力，并加入了水库数据，使土壤保持功能评估的合理性和准确性均得到提升[11,15]。目前，InVEST模型在土壤保持功能方面相对成熟，不仅促进了区域尺度土壤保持功能的量化研究，也使生态系统服务功能综合管理和决策有了坚实的保障。

门头沟区是距离北京市城区最近、植被覆盖度最高的区域，对于减少下游地区河流淤积、防止风沙危害等具有重要作用[16]。历史上以煤矿工业为支柱产业的发展方向对其生态环境打击沉重，作为首都西部生态安全屏障和生态涵养发展区，自2005年以来实施生态林补偿、矿山修复等生态环境保护和生态修复工程，LUCC、生态系统结构与类型变化以及生态系统土壤保持功能受到诸多学者的关注[17-18]。韩永伟等[16]利用USLE和风力侵蚀模型得出门头沟生态系统每年可以减少土壤的水蚀量为14.93万t/a，风蚀量为2 254.38万t/a，植被覆盖度高且坡度小的地区土壤保持能力最强。冯朝阳等[19]采用USLE模型发现门头沟区植被生态系统的土壤保持量为15.7万t/a。周彬等[20]应用InVEST模型对北京山区不同森林类型土壤侵蚀进行模拟研究，发现北京山区林地下总的土壤侵蚀量为176万t，平均保持土壤能力为220 t/hm2。黄从红[18]采用InVEST模型计算2002年和2011年门头沟区生态系统的土壤保持量分别为1.64亿t和1.65亿t，平均土壤保持能力为1 135.19 t/hm2，1 141.44 t/hm2。已有研究鲜少关注LUCC与土壤保持功能时空变化的关系，模型中关于植被与经营管理因子的获取，大多是通过查阅相关研究资料和文献对不同地类进行直接赋值，无法体现植被覆盖度的时空变化分异特征，在土壤保持量的研究中均未提及对地块自身泥沙持留量的计算，且大多研究仅为单一时间尺度，缺少较高空间尺度的土壤保持功能时空变化研究。

开展门头沟区的土壤保持功能的研究，不仅有益于门头沟区防土固沙而且也对整个北京的生态效益的研究至关重要。本研究基于Landsat TM/OLI和HJ-1 CCD遥感影像，分析北京市门头沟区自实施生态修复工程以来，2005—2013年的土地利用与土地覆被的结构、类型与空间变化特征；在此基础上，采用InVEST模型对研究区内生态系统的土壤保持功能进行定量的分析和动态评估，探析LUCC对土壤保持功能的影响，发现并揭示土壤保持功能的变化规律及驱动因素，为门头沟区水土保持研究、综合生态管理以及生态恢复提供数据支持与决策依据。

1 材料与方法

1.1 研究区概况

门头沟区位于北京城区正西偏南，东经115°25′00″—116°10′07″，北纬39°48′34″—40°10′37″，总面积为1 450 km2，海拔为50～2 303 m，山地面积占98.5%，是北京市唯一的纯山区，下辖4个街道、3个地区、9个镇。门头沟区属中纬度大陆性季风气候，降水量自东向西逐渐减少，受中纬度大气环流的不稳定和季风影响，降水量年际变化大，2013年门头沟区平均降雨量为469.3 mm，多年平均为528.70 mm。

北京市2010—2012年开展的第一次水务普查显示[21]，门头沟区的土壤侵蚀面积为396.44 km2，占北京市土壤侵蚀面积的12.38%，侵蚀类型主要是水力侵蚀，较2000年全国第三次土壤侵蚀遥感调查的减少了362.59 km2，说明近年来开展的水土保持与生态环境建设卓有成效。根据2005—2013年《北京市水土保持公报》(2005年为《北京市水土流失监测报告》)公布数据显示，北京山区坡地土壤流失量由2005年的248.3万t增加到2013年的483.07万t，门头沟区的土壤流失量也增加了近35万t，土壤流失量总体呈增加趋势，土壤侵蚀问题仍较为严重(图1)。门头沟区作为北京市西部安全屏障和北京市唯一生态修复科技试验区，自2005年以来实施了一系列封山育林、矿山修复、生态小流域综合治理、京津风沙源治理等生态环境保护和修复工程，截止到2013年，区内林木绿化率达到了62.5%，森林覆盖率达到了41%[22]，植被覆盖度和植被质量显著提高，区内的生态环境获得了极大改善。

1.2 数据基础

本研究主要基于Landsat TM/OLI和HJ-1 CCD遥感影像，通过采用决策树分类方法，获得2005年和2013年两期的土地利用与土地覆被分类数据。考虑到云盖量和植被覆盖度的要求，本研究主要采用秋季(9—10月)的影像。此外，为了获得不同年份逐月的植被覆盖度数据，2005年为Landsat TM遥感影像，2013年遥感影像主要为HJ-1 CCD遥感影像，辅助Landsat OLI遥感影像，保证每月至少一景影像覆盖。因Landsat系列数据重访周期和影像质量问题，缺失月份影像可用上下年份相同月份的影像替代(表1)。

图1 门头沟区和北京山区土壤流失量年际变化表1 门头沟区遥感影像数据

借助门头沟区2000年TM精校正遥感影像和2013年9月的野外采样数据(图2)，对2005年的Landsat TM和2013年Landsat OLI以及HJ-1 CCD遥感影像，采用立方卷积插值和多项式变换的校正方法，控制影像校正误差控制在0.5个像元内，并将UTM投影转换成Albers等面积割圆锥投影。基于ENVI软件平台，采用FLAASH模型对2005年和2013年逐月的遥感影像进行大气校正，为反演归一化植被指数和植被覆盖度参数奠定基础。

DEM数据来源于从中国科学院计算机网络信息中心国际科学数据服务平台，空间分辨率30 m。气象数据来源于中国国家气象局，为2005年和2013年的月均降雨量数据。对数据进行精度验证，剔除不可替代的错误数据，采用ArcGIS软件提供的反距离插值方法，获得空间分辨率为30 m的逐月降雨量和年均降雨量空间分布数据。土壤类型数据来源于北京市农业局。

图2 2013年野外采样点

1.3 土地利用与土地覆被分类

根据门头沟区的土地利用与土地覆被特征，同时考虑生态系统服务功能评价内容以及Landsat TM/OLI和HJ-1 CCD的遥感影像的解译能力，本研究所采用的土地利用与土地覆被分类体系主要参考《全国土地利用现状分类》(2007年)国家标准的一级分类体系，包括耕地、园地、草地、林地、建设用地、水域和未利用地7个地类。基于ENVI软件平台，结合2013年野外采样数据和Google Earth高清影像，建立7种地类的训练样本，每种地类至少50个样本区。采用决策树分类方法，获得2005年和2013年的土地利用与土地覆被分类图，分类精度分别为87.56%和89.23%，均可以满足空间分析的要求(图3)。

图3 门头沟区土地利用与土地覆被分类

1.4 模型原理与因子计算

沉积物保留模型(Sediment Delivery Ratio model)的土壤保持量包括土壤侵蚀减少量和泥沙持留量两部分，前者反映各地块对自身潜在侵蚀的减少，以潜在侵蚀与实际侵蚀的差表示；后者表示该地块对进入它的上坡来沙的持留，以来沙量与泥沙持留率的乘积表示[23]。模型计算公式如下：

SEDRETx=Rx×Kx×LSx×(1-Cx×Px)+SEDRx

(1)

(2)

USLEx=Rx×Kx×LSx×Cx×Px

(3)

式中：SEDRETx和SEDRx分别为栅格x的土壤保持量和泥沙持留量；USLEx和USLEy分别为栅格x及其上坡栅格y的实际土壤侵蚀量；Rx代表栅格x的降雨侵蚀力因子；Kx代表栅格x的土壤可蚀性因子；LSx代表栅格x的坡度坡长因子；Cx代表栅格x的植被与经营管理因子；Px代表栅格x的水土保持措施因子；SEx代表栅格x的泥沙持留效率。

本研究采用Wischmeier等[9]提出的月尺度计算降雨侵蚀力(R)计算公式，获得门头沟区2005年和2013年降雨侵蚀力空间分布图(图4)。

(4)

式中：Yi为月平均降雨量(mm)；Y为年平均降水量(mm)。

考虑到本研究区的土壤可蚀性(K)已有专门的研究成果[20,24]，根据土壤数据中包含的土壤类别、名称及其分布图，将土壤类型与已有的K值关联起来，获得门头沟区土壤可蚀性分布图(图5)。

由于InVEST模型原始的计算公式是在美国坡度较小的试验区上进行测定，因此需要对坡度坡长因子(LS)因子进行相应的调整，参照已有研究成果，本研究将边坡阈值设为25°[20]。

图4 门头沟区降雨侵蚀力因子

图5 门头沟区土壤可蚀性因子

当坡度小于边坡域值时，坡度坡长的计算公式为：

(5)

(6)

当坡度大于边坡域值时，坡度坡长的计算公式为：

LS=0.08β0.35s0.6

(7)

式中：Fa为汇水累积量阈值；Cs为栅格大小；s为坡度；n为与坡度有关的参数；β为根据流向确定的栅格大小值。

根据对已有的门头沟区河流水系的现状分析，基于干流全部提取、支流最多提取，同时考虑到各支流的汇水量情况，本文通过反复试验，确定汇水累积量阈值为10 000，经过InVEST模型计算的门头沟区坡度坡长空间分布图(图6)。

本文采用蔡崇法等[25]的植被与经营管理因子(C)研究方法，反映植被或作物管理措施对土壤侵蚀量的影响。计算公式为：

(8)

fc=(NDVI-NDVIs)/(NDVIv-NDVIs)

(9)

NDVI=(ρnir-ρred)/(ρnir+ρred)

(10)

式中：C为植被与经营管理因子，其值介于0～1；当C=0时表示土地表面的植被覆盖度良好，土壤几乎没有受到侵蚀；当C=1时表示土地表面是没有植被覆盖即完全裸露的。fc为植被覆盖度；NDVI为植被归一化指数；NDVIs为完全是裸土或无植被覆盖区域的NDVI值；NDVIv为完全被植被所覆盖的像元的NDVI值，在实际计算中，取NDVI最大值代替NDVIv和最小值代替NDVIs；ρred为遥感影像数据红波段反射率，ρnir为遥感影像数据近红外波段反射率。

图6 门头沟区坡度坡长因子

本研究为获得高精度的2005年和2013年的年植被覆盖度数据，在获得逐月植被覆盖度的基础上，采用最大值合成方法，获得两个年份的植被覆盖度数据(图7)。根据公式(8)计算获得两个年份的植被与经营管理因子(P)空间分布图(图8)，结合土地利用与土地覆被数据，统计不同地类的植被与经营管理因子(C)的平均值。

图7 门头沟区植被覆盖度空间分布

土壤保持措施是指土地采取相应的水土保持措施后的土壤侵蚀量与顺坡种植时的土壤侵蚀量的比值(P)，其值介于0～1。当P值为0时表示该区域土地表面采取了较好的水土保持措施，土壤几乎不会受到侵蚀；当P值为1时表示该区域土地表面没有采取任何水土保持措施。本研究通过实地状况以及查阅相关文献资料得到相应P值[26]，耕地、园地、有林地、灌木林地、草地、水域、建设用地、未利用地分别为0.4，0.7，1，1，1，0，0，1。

泥沙持留效率反映了侵蚀产生的泥沙在输移过程中因植被过滤、拦截等作用而发生沉积的过程，被拦截泥沙比例越大，持留效率越高[23]。本研究参照InVEST模型数据库获得不同植被类型的泥沙持留效率。

图8 门头沟区植被与经营管理因子空间分布

2 结果与分析

2.1 土地利用与土地覆被变化分析

基于2005年和2013年两期的土地利用与土地覆被分类数据，从数量变化、类型变化和空间变化3个尺度，分析门头沟区自实施生态修复工程以来(2005年)，区域的土地利用与土地覆被变化特征，并分析其驱动因素与驱动机制。

根据对门头沟区土地利用与土地覆被结构特征的分析(表2)，发现2005年和2013年，门头沟区林地占据绝对主导的地位，是土地利用的基质与最主要的土地覆被类型，面积比例占到89%以上，有林地与灌木林地的比重接近1.35∶1；建设用地、耕地和园地为主导的土地利用类型。基于ArcGIS的空间叠加分析，2005—2013年的土地利用与土地覆被类型和空间变化(表3)，发现门头沟区耕地、灌木林地、水域和未利用地面积呈减少趋势，园地、有林地、草地、建设用地呈增加趋势，其中耕地为减速最快的地类，园地为增速最快的地类；门头沟区的土地利用动态变化的过程中，存在着频繁的地类间的转化现象，耕地、园地、草地、水域转化较为显著，主要转为建设用地、园地和耕地；发生转化的地类主要位于地势平坦、交通便利、土壤质量较高的地区。

具体来看，耕地面积急剧减少，是减少最明显和减速最快的地类，结构比重下降了1.1%，年均减少率为3.73%，主要转为建设用地、草地和园地，三者的转移比重占总转移比例的78.03%，减少的耕地主要位于清水镇、斋堂镇和雁翅镇，主要转为草地和园地，主要受退耕还林还草等政策影响；此外，永定镇的耕地主要转为草地，受经济发展驱动影响显著。园地面积增加了679.06 hm2，是增速最快的地类，年均增长率为5.97%，主要由耕地转化而来，贡献率达到65.66%，新增园地主要位于雁翅镇和清水镇，形成了以龙泉务为主的香白杏产业带、以东山、孟悟为主的京白梨产业带、以潭柘寺和陇驾庄为主的盖柿产业带、以妙峰山樱桃沟为主的樱桃产业带、边远山区仁用杏核和薄皮核桃种植带、以太子墓为主的苹果种植带，主要受沟域经济发展以及农业结构调整等政策的影响。林地面积减少了249.06 hm2，是减速最慢的地类，主要是由于灌木林地大面积的减少，有林地的增加主要来源于灌木林地和耕地，减少的灌木林地主要位于东部和南部地区，主要转为有林地和建设用地，这主要是受生态公益林建设与改造、退耕还林等政策的影响。草地面积增加了400.21 hm2，增速仅次于园地，主要来源于耕地和建设用地，主要转为耕地和林地，增加的草地主要位于斋堂镇北部，主要是因为退耕还林还草以及生态搬迁等政策影响。建设用地增加了837.70 hm2，遵从"一主一副两重点"的城镇发展模式，门头沟区的建设用地主要集中在东部地区，主要由耕地转变而来。

2.2 土壤保持功能评价与驱动因子分析

采用InVEST模型中的土壤保持模块，通过导入相关因子，获得门头沟区的2005年和2013年的现实土壤侵蚀量、潜在土壤侵蚀量及土壤保持量的数据(图9)。在此基础上，从土壤保持空间分布、土壤保持功能变化、不同土地利用与土地覆被类型的土壤保持能力3个尺度对比研究，分析门头沟区自实施生态修复工程2005年以来，生态系统土壤保持能力的变化规律与特征。

2.2.1 土壤保持量及其空间分布 2013年，门头沟区生态系统平均现实土壤侵蚀量为7.46 t/(hm2·a)，平均潜在土壤侵蚀量为121.00 t/(hm2·a)，平均土壤保持量为113.56 t/(hm2·a)，土壤保持总量达到1 638.30万t。2005—2013年，门头沟区土壤保持能力呈增加趋势，8年间平均土壤保持能力提升了26.15 t/(hm2·a)，土壤保持总量增加了377.38万t，虽然受降雨侵蚀力增大影响，由2005年的2 395.25 (MJ·mm)/(hm2·h·a)增加到2013年的3 144.80 (MJ·mm)/(hm2·h·a)，潜在土壤侵蚀总量增加了400.80万t，但由于植被覆盖度和植被与经营管理措施的提高和改善，现实土壤侵蚀量较2005年的84.17万t，减少了23.42万t，水土保持措施初显成效。

表2 门头沟区土地利用与覆被类型的面积变化情况

表3 门头沟区土地利用与土地覆被类型转移概率矩阵 %

门头沟区北部、南部、西部和东北地区多为海拔比较高的山地，其土地利用类型主要是林地和草地，而土壤保持量大小与林地的覆盖率密切相关，林地覆盖率高的区域的土壤保持量明显高于其他地区的土壤保持量，所以这些地区土壤现实侵蚀量较小，土壤保持情况较好。中部、东南部地区地势较为平坦，土地利用类型主要是耕地、园地和建设用地，这些地区土壤保持量较低，大力发展特色种植业以及“一主一副两重点”的城镇发展模式，因采取适当的植被经营与管理工程措施，土壤保持能力相对较高。2005—2013年，西北部和东北部地区受植被覆盖度显著提高影响，土壤保持呈现增长趋势；西南地区在降雨侵蚀力减弱主导下，土壤保持呈现减少趋势；东部地区虽然降雨侵蚀力有所增强，但地类主要转为建设用地、耕地和园地等类型，因此总体土壤保持能力变化不大。

2.2.2 不同土地利用与土地覆被类型土壤保持能力 门头沟区自2005年以来实施了多项生态修复工程，包括退耕还林、草等，促使大量耕地面积减少，并转化为林地、草地；受北京山区发展沟域经济驱动，门头沟区大力发展特色种植业，导致园地面积大量增加，从而使得其植被覆盖度增加。不同的土地利用与土地覆被类型，受植被覆盖度和面积变化的影响，土壤保持能力差异显著。2013年，林地的单位面积土壤保持能力最高，达到121.36 t/(hm2·a)，其次是草地；水域的最低，仅为25.04 t/(hm2·a)；耕地和园地具有相似的土壤保持能力，介于42.20～43.20 t/(hm2·a)；单位面积土壤保持能力排序依然为：林地>草地>未利用地>建设用地>耕地>园地>水域(图10)。2005—2013年，不同土地利用/覆被类型的单位面积土壤保持能力均呈增加趋势，其中未利用地增速最快，水域增速最慢；耕地受面积减少显著影响，土壤保持总量减少了2.57万t，其他地类的土壤保持总量均呈增加特征；虽然林地、水域、未利用地面积呈减少趋势，但因3种地类的单位面积的土壤保持能力增长显著，因此其土壤保持总量均呈增加趋势。

图9 门头沟区土壤保持量空间分布

图10 门头沟区不同土地利用类型土壤保持能力

3 讨论与结论

3.1 讨 论

由于模型中相关因子并不适合门头沟的实际情况，因此本研究对InVEST模型相关参数与计算方法进行了校正、检验，使之适应于门头沟区的生态系统服务功能评估。本研究仅针对生态系统服务的土壤保持功能进行了动态评估，但自然生态系统供给的服务多种多样，其他的服务功能如碳储存、水源涵养、生物多样性功能也同样重要，未来有必要针对这些服务功能进行动态评估，从而实现对区域生态系统服务功能的全面了解。本研究只利用物质量法计算了门头沟区生态系统的土壤保持功能，并没有进一步计算它们的价值量。价值化后生态系统服务功能可以与生态系统提供的直接经济价值进行对比，促使人们加大对生态系统的保护力度，因此门头沟区生态系统土壤保持和价值量估算应该成为后期研究中的重点。相对于过去，应该更关注生态系统服务功能在未来的变化趋势，有必要对未来门头沟区生态系统服务功能的变化趋势进行评估。

3.2 结 论

(1) 门头沟区景观格局特征是林地为主导的土地覆被类型，建设用地、耕地和园地为主导的土地利用类型。2005—2013年，耕地、灌木林地、水域和未利用地面积呈减少趋势，园地、有林地、草地、建设用地呈增加趋势，其中耕地为减速最快的地类，园地为增速最快的地类；耕地、园地、草地、水域地类间的转化较为显著，主要转为建设用地、园地和耕地；发生转化的地类主要位于地势平坦、交通便利、土壤质量较高的地区。生态公益林建设与改造、退耕还林还草、农业结构调整、生态搬迁以及山区绿色生态经济发展等政策，是驱动门头沟区土地利用与土地覆被变化的主要因素。

(2) 2005—2013年，门头沟区生态系统土壤保持能力呈增加趋势，8年间平均土壤保持能力提升了26.15 t/(hm2·a)，土壤保持总量增加了377.38万t，由于植被覆盖度和植被与经营管理措施的提高和显著改善，虽然降雨侵蚀力增加了749.55 MJ·mm/(hm2·h·a)，但现实土壤侵蚀量减少了23.42万t，自2005年以来实施生态修复与水土措施初显成效。西北、东北部地区主要是林地和草地类型，植被覆盖率提高显著，土壤保持呈现增长趋势，西南地区主要受降雨侵蚀力减弱影响，土壤保持呈现减少趋势，东部地区则在降雨侵蚀力增强和土地利用与土地覆被类型快速转变的综合影响下土壤保持能力变化不大。

(3) 土壤保持功能与土地利用与土地覆被类型密切相关。不同土地利用与土地覆被类型单位面积土壤保持能力排序依然为：林地>草地>未利用地>建设用地>耕地>园地>水域；2005—2013年，不同土地利用/覆被类型的单位面积土壤保持能力均呈增加趋势，其中未利用地增速最快，水域增速最慢；除耕地受面积减少显著影响外，其他地类的土壤保持总量均呈增加趋势。因此，可通过调整土地利用与土地覆被的空间结构与人为干预土地利用与土地覆被变化，为门头沟区及京津冀山区的土地利用科学管理及决策提供参考。