生态修复措施对干旱半干旱地区生态系统服务影响研究——以河北坝上地区为例

裴宏伟, 刘孟竹, 李雅丽, 张红娟, 肖雨霄, 杨国丽

(1.河北建筑工程学院 市政与环境工程系, 河北 张家口 075000;2.河北省水质工程与水资源综合利用重点实验室, 河北 张家口 075000)

生态系统服务是指人类从生态系统中获得的有形或无形的各种惠益，主要分为供给服务、调节服务、文化服务和支持服务4种类型[1]。根据生物多样性和生态系统服务政府间科学与政策平台(IPBES)2019年的评估，受人类活动的影响，全球生态系统正以“史无前例”的速度退化[2]。中国自改革开放以来，虽然经济上取得了显著增长却也损害了生态环境(如侵占自然用地的城镇化扩张[3])。对此，中国于2000年起便强力地推进一系列的生态修复工程以应对生态环境退化。Ouyang等[4]学者在对中国2000—2010年生态环境的评估中指出，生态修复工程的实施极大地增强了生态系统功能与服务供给，这一观点在区域尺度的研究中也得到了佐证[5]。在相关议题上，多数研究更加关注生态修复工程对过去时期生态系统服务的影响，如在特定区域(退耕还林区、生态保护区、湿地保护区)内分析生态修复措施前后的生态系统服务变化[6]，而对这种影响在未来的变化却很少给予关注。一般而言，生态修复工程的实施往往是逐阶段进行的，评估未来生态修复措施对生态系统服务的影响不仅有助于为生态恢复工程的规划提供决策依据，同时也能更好地理解土地利用/覆被变化对生态系统服务的影响机制。目前，几个大型的生态修复工程主要集中在中国北方的干旱或半干旱地区，采取的生态修复措施主要为大规模的人工种树。已有研究表明，在干旱半干旱这样的水分限制地区(降水量小于潜在蒸散发量)，大规模的植树(多为经济效益高的单一树种)反而会减少生态系统的供水服务以及生物多样性，进而造成生态环境退化[7]。因此，明确以植树为主的生态修复措施对未来生态系统服务的影响极为关键，对当下的生态管理及决策具有非常意义。

情景分析能够有效地评估不同生态恢复措施情景下的生态系统服务。不同情景的构建有多种选择，例如通过固定气候因子以过去时期的土地利用作为对比情景[8]；或者通过预测模型(CA-Markov，CLUE-S模型)来模拟特定发展方向(生态保护、城市扩张、耕地保留情景)下未来时期的土地情景[9]，以及设计未来的土地情景[10-11]。其中，通过设计土地情景，既能模拟出符合生态修复措施的未来土地情景，也能明确某一种或某几种土地利用类型变化对生态系统服务的影响。

对生态系统服务进行建模的模型常用的包括ARIES模型、EcoAIM模型、InVEST模型。其中，ARIES模型可以对生态系统服务流(生态系统到人的传递)中的关键要素(“源”、“汇”和“受益者”)的空间和数量进行建模，但只能应用于生物物理关系的解释[12]；EcoAIM模型可通过结合一系列的空间数据进行加权计算来绘制生态系统服务数值的空间分布，但该模型无法作为独立程序进行应用[13]；InVEST模型因其输入数据量少、过程相对简化而能够快速地量化结果，在区域乃至国家尺度上均受到相应的使用[14]。例如，Shoyama等[15]利用InVEST模型绘制了日本北部钏路农村流域的生态系统服务(产水量、碳储存、生境质量)；刘孟竹等[16]基于InVEST模型评估了中国北方农牧交错带碳储量变化；欧阳志云等[17]在全国生态环境10 a变化遥感调查评估中也使用了InVEST模型。InVEST模型能够对不同的生态系统服务进行多时空尺度的评估，为决策者提供生态管理相关的依据。

河北坝上地区地处我国干旱半干旱区过渡区，是典型的生态脆弱区。该区作为首都“生态屏障”“水源涵养区”肩负极其重要的生态意义和地理区位[18]。自退耕还林工程实施以来，坝上地区经历了一段长期的大规模蔬菜种植(蔬菜耗水量极大)，在该背景下，过量的地下水开采以及高强度的农业灌溉等导致坝上地区地下水位明显下降[19]，水资源稀缺加剧，生态可持续发展也受到了严重制约。为应对现有的生态危机，2019年国家发展与改革委员会联合河北省人民政府作出明确规划，旨在加强坝上地区发展为水源涵养功能区和生态环境支撑区的“两区”建设[20]。目前，针对坝上地区生态修复措施对生态系统服务影响的研究仍然缺乏，相关研究仅针对坝上地区的生态系统服务价值[21]和单一的服务[22]作出了评估。基于此，本文通过InVEST模型和土地利用数据，分析坝上地区2000—2018年土地利用以及3个重要的生态系统服务(水源涵养、碳储量和生境质量)的时空变化，并通过设计土地情景揭示未来生态修复措施对生态系统服务的影响，以期为当地的生态可持续发展提供理论支撑。

1 研究区概况

坝上地区位于河北省西北部(40°70′—42°60′N，113°70′—117°90′E)，涉及尚义县、沽源县、张北县、康保县、围场县和丰宁县共6个区域(图1)，总面积约1.9万km2，海拔831～2 215 m(填洼后)。坝上地区多为季节性河流，湖多河少，较大的湖有安固里淖、大青湖等，降水和径流是该区湖泊补给的主要方式。坝上地区属大陆性季风气候，年均气温1～2℃，多年平均降水量400 mm左右，70%集中在6—9月，年均蒸发强度高达1 800 mm[23]，属于典型的受水资源限制的干旱半干旱地区。坝上大部分区域属于半干旱草原地带，包括干草原和湿地草甸草原；植被类型主要有旱地作物、草地、原生林和人工防护林等；土壤以栗钙土、沙质栗钙土为主，有机质含量少且土质松软，易受侵蚀风化。根据《张家口经济年鉴2018》，坝上主要4县(尚义县、沽源县、张北县、康保县)2017年生产总值为240.96亿元，总人口为106.68万人。

图1 研究区海拔

2 数据与方法

2.1 数据源及其处理

本文中的数据介绍如下：(1) 土地利用数据来源于中国科学院资源环境科学数据中心，分辨率为30 m，时间为2000年和2018年，本研究参考坝上地区的实际情况将其分为6个一级分类(耕地、林地、草地、水域、建设用地、未利用地)；(2) DEM数据来源于地理空间数据云(http:∥www.gscloud.cn)GDEMV2 30 m分辨率数字高程产品，经裁剪后进行了填洼处理；(3) 土壤数据来自于国家冰川冻土沙漠科学数据中心(http:∥www.crensed.ac.cn/portal/)1∶100万世界土壤数据集；(4) 气象数据来源于中国气象数据网(http:∥data.cma.cn)，包括坝上及周边地区6个气象站2000—2018年降雨、气温等数据，年降雨量与年潜在蒸散量栅格图由克里金插值法所得；(5) 次级流域矢量图由ArcMap软件中水文分析工具提取得来；(6) 威胁源(耕地、城市、农村、交通)数据使用土地利用数据中的二级分类；(7) 碳密度数据来源于同类研究，详见研究方法。

2.2 研究方法

2.2.1 水源涵养量计算 由InVEST模型中的产水量模块求得水源涵养量，模型基于水量平衡原理，在考虑了地形因子和地表径流等因素对水源涵养量进行评估。具体公式[8]如下：

(1)

(2)

式中:WRij为栅格i中土地利用类型j的年水源涵养量(mm/a);Yij是第j类土地利用类型在栅格i的年产水量(mm/a);V是流速系数,为常数;TI是地形指数,无量纲;Ksat为土壤饱和导水率(cm/d),由Neuro Theta软件计算得到;Drainage_Area为集水区栅格数量;Soil_Depth为土层深度(mm);Percent_Slope为百分比坡度。

(3)

(4)

PETi=Kc,j·ET0,i

(5)

(6)

AWCi=min[max(layer_depthi)，Root_depthi]·PAWCi

(7)

式中:Pi为栅格i的年降雨量(mm/a);Yij为栅格i中土地利用类型j的年产水量(mm/a);AETij为栅格i中土地利用类型j的年实际蒸散量(mm/a);PETi为栅格i的年潜在蒸散量(mm/a);ET0，i为栅格i的年参考蒸散量(mm/a);Kc,j为土地利用类型j的蒸散系数,无单位;wi为气候—土壤的非物理参数;Z为季节性常数,参考同地区产水服务研究[24]取值4.02;max(layer_depthi)为栅格i的最大土层深度(mm);Root_depthi为栅格i的最大根系限制深度(mm);AWCi为土壤有效含水量(mm);PAWCi为栅格i的植物可利用水分含量,值为[0,1]。

2.2.2 碳储量计算 碳储量模型中需要的必要输入数据包括土地分类数据以及基本的四大碳库数据，即地上碳库、地下碳库、土壤碳库和死亡有机质碳库，输出可得到土地利用数据包含范围内的固碳总量和不同土地类型固碳量。由于所选研究区较大，除基本碳库以外的碳库数据获取上具有困难，故只考虑4大基本碳库，碳密度数据(表1)以及计算公式均参考针对河北省碳储量的研究[25]如下：

C=(Ca+Cb+Cc+Cd)×S

(8)

式中:C为总碳储量(t)；Ca为单位面积地上碳密度(t/hm2)；Cb为单位面积地下碳密度(t/hm2)；Cc为单位面积土壤碳密度(t/hm2)；Cd为死亡有机质碳密度(t/hm2)；S为各地类面积(hm2)。

表1 坝上地区土地利用类型碳密度参数 t/hm2

2.2.3 生境质量计算 利用InVEST模型生境模块来计算生境质量以及生境退化度。生境退化指数反映单元栅格受胁迫程度的大小，值越高代表受胁迫程度最大，生境退化程度越高，其公式[14]如下：

(9)

(10)

(11)

式中:Dxj为土地利用类型j中x栅格的生境退化度;R为威胁源个数;wr为危险源权重;Yr为威胁源的栅格数;ry为栅格y的胁迫值;irxy为栅格y的胁迫值ry对栅格x的胁迫水平;Bx为危险源对栅格x的可达性;Sjr为土地利用类型j对危险源r的敏感度;dxy为栅格x到栅格y的直线距离;drmax为危险源r的最大胁迫距离。

生境质量指数反映在胁迫条件下生境质量的优劣，其值在0～1。值越高代表区域生物多样性越丰富，生境质量越好；反之则生境质量越差，易受破坏，计算公式[14]如下：

(12)

式中:Qxj为土地利用类型j中x栅格的生境质量指数;Hj为土地利用类型j的生境适宜度;k为半饱和常数,取Dxj最大值的一半;z为归一化常量,通常取2.5。

本文将林、草地及水域定义为提供生境质量的地类，将耕地、村落、城镇用地、交通用地设为威胁源，各个威胁因子最大距离、权重、衰退类型参考文献[14]。各个地类生境适宜度、对威胁因子敏感度的设定详见表2—3。

表2 坝上地区威胁源参数

表3 坝上地区土地利用类型生境适宜度及威胁因子敏感系数

2.2.4 生态修复措施情景设置 为了分析相关的生态修复政策对坝上地区生态系统服务的潜在影响，本文以2018年土地利用情景为基础情景，在此基础上设计了4种未来生态修复情景，分别为水源缓冲带、植树造林、开垦荒地和综合发展等，具体见表4。

3 结果与分析

3.1 坝上地区2000-2018年土地利用变化分析

近20 a来，耕地是坝上地区最主要的土地利用类型(表5)，其面积占研究区总面积近49%，在研究期间呈轻微下降趋势(-1.07%)。草地和林地面积占比分别处于第二、第三，2000—2018年，草地缩减520.57 km2(-10.52%)，林地张602.61 km2(19.91%)。水域、建设用地和未利用地在整个研究期间占比均较小，但水域、建设用地相对变化较高；在2000—2018年，水域下降了31.75%，建设用地增长了43.16%。从土地数量变化来看，坝上地区整个研究阶段土地利用变化的主要特点是“两增四减”，即，林地面积显著增加，建设用地增幅明显，草地面积显著减少，水域面积降幅显著，耕地、未利用地面积减少不明显。

表4 生态修复措施情景设置

表5 坝上地区2000-2018年土地利用情况

由表6中土地利用转移可知，坝上地区2000—2018年退耕还林、还草面积分别达到230.97，191.19 km2，近201.59 km2的建设用地主要来自农田。同期，约295.63，534.96，130.25 km2的草地分别转变为耕地、林地和未利用地。此外，水域退化为未利用地的面积达到68.48 km2，这主要是由于安固里淖干涸造成的。未利用地转变为草地、耕地的面积分别达到140.56，76.05 km2，可以表明当地对未利用地的治理也取得了一定成效。综上而言，当地政府在植被恢复工作中取得了一些进展，但需要注意的是，林地的扩张几乎来源于草地的转入，未来的决策中应当提高对草地保护的关注。

表6 坝上地区2000-2018年土地利用转移矩阵 km2

3.2 生态系统服务评估

2000年、2018年坝上地区水源涵养量分别为10.71亿，9.96亿m3，水源涵养深度分别为56.79，52.80 mm，近20 a来下降约7%(表7)。分区来看，围场县植被覆盖水平较高，水源涵养功能较强(水源涵养深度最高为100.07 mm，下同)，丰宁县则其次(58.80 mm)；康保县和张北县耕地多、植被少，但常年雨水较为充沛，水源涵养功能较其余县处于中等水平，水源涵养深度分别为54.18，54.65 mm；尚义县(44.23 mm)和沽源县(38.25 mm)水源涵养量较低，主要受该区多年降雨量偏低的影响。另外，沽源县未利用地较多，导致该区域的水源涵养量整体偏低。

坝上地区碳储量在2000年、2018年分别为1.82亿t，1.87亿t，近18 a来增加了5.13×106t(2.82%)，单位面积碳储量由94.79 t/hm2增长到97.47 t/hm2。分区来看，丰宁县碳储量水平增长最显著(15.36 t/hm2)，主要是由于该区新增了较多植被，导致植被净初级生产力增强，从而提高了碳储量；整个研究区仅有康保县碳储量为减少趋势(-2.15 t/hm2)，这是由于该区未利用地较多以及城镇用地扩张，导致碳储量处于较低水平并呈减少趋势；此外，其余县区均处于增长趋势但程度不明显。

坝上地区生境质量在2000—2018年期间呈增长趋势，平均生境质量指数由0.49上升到0.51。近20 a来，丰宁县平均生境质量指数上升最显著(12%)，这主要是由于该区林地面积的增长较多，因此给当地生态系统提供了更多的动植物自然栖息地；沽源县和尚义县平均生境质量指数则处于下降趋势(-2%)，这是由于该区内广泛分布的未利用地以及扩张的建设用地增加了更多的威胁源，并缩减了生境适宜度高的自然栖息地，因此导致生态系统生物多样性减少，生境质量下降。

表7 坝上地区2000-2018年生态系统服务变化

3.3 生态修复措施对生态系统服务的影响

以2018年基础情景为基准水平，在林地缓冲带、植树造林、开垦荒地以及综合发展情景下，水源涵养深度分别提高了0.61%，1.38%，4.26%，6.12%(图2)，对应的总量增加分别为6.13×106，1.388×107，4.286×107，6.157×107m3。在林地缓冲带情景下，林地面积相对2018年水平增加了约4%，水源涵养量的增长非常轻微(0.61%)；植树造林情景下林地面积增加了约7%，水源涵养量增长1.38%。可以认为，造林为主的生态修复措施对该区水源涵养服务提高不够明显。在开垦荒地情景和综合情景下水源涵养量分别提升了4.24%，6.10%，这两种措施对水源涵养服务的提升作用较为明显。相比较而言，综合情景下的水源涵养量变化主要是受开垦荒地情景主导，在该情景下，未利用地减少而耕地增加，植被覆盖水平的增加提高了地表土壤的持水、固土能力，水源涵养服务也因此提高。开垦荒地情景对水源涵养量的提升作用明显高于林地缓冲带情景和植树造林情景，这是由于坝上地区土地利用结构以及地形的特殊性所致。坝上地区的未利用地占比约5%，且多分布在坡度平缓、海拔较低的区域。在开垦荒地情景下，大量未利用地(近94.47%)被整改为耕地，另外，当未利用地整改为耕地后，土壤的持水性、固土性会得到增强，这综合地导致了坝上地区以扩张耕地为主的生态修复措施比造林为主的措施在提高水源涵养服务上更明显。

注：S0为基础情景，S1为水源缓冲带，S2为植树造林，S3为开垦荒地，S4为综合发展。

碳储量在水源缓冲带、植树造林、开垦荒地以及综合发展情景下分别增长了3.78%，4.01%，5.81%，7.39%，单位面积碳储量分别达到了98.04，98.26，99.96，101.45 t/hm2。经分析，在水源缓冲带、植树造林情景下，新增林地仅占整个研究区面积的0.78%，1.38%，而碳储量增长均接近4%；开垦荒地情景下新增耕地占比为4.58%，相应的碳储量增长达到5.81%。可以认为，林地扩张的情景对碳储量提升更显著。一般而言，林地具有较高的植被覆盖水平，植物叶片经过光合作用固定的碳量远高于自身呼吸作用排放的碳量，林地净初级生产力较草地、耕地等更高，因此林地扩张情景对固碳功能的促进作用远高于耕地扩张情景。在综合情景下，碳储量有了较显著的上升，综合措施能够对坝上地区固碳、气候调节等服务均有明显的促进作用。

生境质量指数一定程度上反映的是生态系统的物种多样性以及自然栖息地的提供水平。在水源缓冲带、植树造林、开垦荒地以及综合发展情景下，坝上地区平均生境质量指数分别为0.510，0.514，0.519，0.534，相比于2018年基础情景增长了0.99%，1.78%，2.77%，5.74%。经分析，耕地、建设用地作为威胁源，当距离人类生活区越近时，对生态系统的干扰程度越大，生境质量指数相应地也就越低。在水源缓冲带、植树造林情景下，耕地、未利用地被林地代替，这既减少了威胁源数量，同时也提高了生态系统生境适宜度，进而促进了动植物栖息地面积和物种多样性的增加，因此生境质量指数得到了提高。开垦荒地情景中，未利用地到耕地的土地转变提高了区域内的生境适宜度，但也增加了威胁源数量，最终的生境质量水平较基准水平仍有所提高。

3.4 生态系统服务权衡与协同分析

生态系统服务的权衡关系是土地规划决策中必须重视的问题。本研究量化了由生态修复措施引起的土地利用类型改变对生态系统服务之间权衡与协同关系变化的影响。表8为不同生态修复措施情景下碳储量、生境质量和水源涵养服务在栅格尺度上两两之间的相关性，所有相关性系数均通过了p<0.01的显著性检验。总的来看，碳储量与生境质量之间、生境质量与水源涵养之间均存在较强的协同作用(r>0.42)，这表明随着生境质量提高，水源涵养与碳储量服务也随之增加；水源涵养与碳储量之间存在着较弱的协同作用(r<0.24)，这表明当水源涵养服务增加，碳储量尽管会增长但程度不明显。在水源缓冲带、开垦荒地措施情景下，生态系统服务之间协同作用相较于基准情景变化不大，在植树造林、综合发展措施情景下，这种协同作用均会变弱。相比于水源缓冲带与开垦荒地措施情景，植树造林与综合发展措施情景均采取了高坡度区域林地的扩张措施，可以认为，由生态修复措施所引起的土地利用转变对生态系统服务权衡与协同作用关系的影响主要取决于土地利用变化区域的地形，而并非变化的面积大小。

表8 不同情景下生态系统服务之间相关系数

4 讨论与结论

4.1 讨 论

从本文的结果来看，耕地为坝上地区主要的土地利用类型(其次为林地和草地)以及草地缩减、建设用地增加等结果与文献[21]研究结论一致，本文结果具有可信度。在生态修复的相关研究中，汪言在等[22]模拟了坝上地区在林地恢复情景下的土壤侵蚀程度，发现该情景有效地减少了土壤侵蚀(6%～14%)从而提高了坝上地区土壤保持服务。在坝上临近地区——潮河流域，吴一帆等[26]研究中发现该流域在河岸缓冲带、退耕还林下水源涵养量分别增长7%，20%，这远高于本研究中的结果。经分析，该研究将河岸缓冲带距离设置为1 km(本研究为0.1 km)，且在退耕还林情景中将未利用地、耕地一同转变为林地(本文仅考虑耕地)，不同的土地变化数量导致了与本研究水源涵养量增长的差异。另外，在同地区且更大的区域，相关研究均论证了生态修复措施对碳储量[16]、生境质量等[6]生态系统服务的提升作用，与本文结果基本一致。

从本文的误差分析，在数据处理上，由于二级分类中更细化的土地利用类型对应的径流系数等参数不易获取，土地利用数据仅采用了一级分类标准分为6个一级土地利用类型，该处理对结果的真实模拟产生了一定误差；此外，张北县的安固里淖于2004年已经干涸，然而土地利用数据中对应的区域仍然是水域。因此，本文对数据进行了人工纠正，将该区域的水域改变成未利用地，消除了数据来源误差。从情景分析来看，林地扩张措施(林地缓冲带、植树造林)并没有显著地提高水源涵养量。一般而言，林地具有较好的雨水下渗性和持水能力，也更适合于涵养水源。然而，林地较高的蒸散能力也会消耗区域内的水资源，导致涵养水源量偏低，相关研究已表明黄土高原地区由于人工造林导致的产水量平均减少幅度达到50 mm/a[27]；另外，就本研究区而言，坝上地区处于干旱半干旱地区，年降雨量约400 mm，但蒸发能力却高达1 800 mm，在这样的水分限制区开展造林活动是需要谨慎考虑的。一般而言，林地较长的根系在缺少雨水的补给下会吸收深层的地下水而导致浅层植物缺水而死亡，从而导致土壤失去保护更易受到侵蚀，使得土壤失去涵养水源的能力。需要强调的是，本研究仅以2018年为基准设计了生态修复措施在未来实施完成的情景，一般而言，林地在不同的生长阶段发挥的生态功能并非完全相同，因此未来的研究会更侧重于对时间尺度的考量。

未利用地在坝上地区占有一定的比例，由于其大部分分布在坡度平缓、海拔低这样适合农田耕作的区域，因此也非常适合于开垦荒地。在开垦荒地措施下，由于几近全部的未利用地转变为了耕地，农业生产将得到较大的提高，但农业的过度扩张也会导致水源的污染和物种多样性的丧失，而且在实际中将未利用地改造成耕地难度大、成本高等问题也有待解决。可见，坝上地区的农业开垦措施需要合理、慎重开展。总体而言，InVEST模型可以作为一个快速、多尺度评估生态系统服务的工具，但对于生态系统能流、物流、信息流的刻画仍缺乏详细的过程解释，未来的情景分析中应综合地考虑生态系统服务之间的权衡影响以及生态过程的其他影响因素。

4.2 结 论

(1) 2000年、2018年坝上地区土地利用结构上呈现“耕地为主，林、草地为辅”，耕地、林地、草地面积平均占比分别为49.54%，17.36%，24.46%；土地利用数量变化上表现为“两增四减”，即，林地面积明显增加(19.91%)，草地面积明显减少(-10.52%)，建设用地增幅显著(43.16%)，水域降幅显著(-31.75%)，耕地、未利用地减少相对不明显。

(2) 2000—2018年，林地的显著扩张整体提高了坝上生态系统服务的供给，尽管水源涵养服务有所减少，水源涵养量由10.71亿m3下降到9.96亿m3，但碳储量和生境质量增强，碳储量由1.82亿t增长到1.87亿t，平均生境质量指数由0.49增长到0.51。

(3) 生态修复措施中，当林地扩张4%～7%时，林地缓冲带、植树造林等措施能够提高碳储存(3.8%～4.0%)，但对水源涵养(0.6%～1.4%)和生境质量(1%～2%)提升不明显；开垦荒地和综合发展措施对水源涵养(4%～6%)、碳储量(4.6%～5.8%)和生境质量(2.8%～5.7%)的提升均相对较明显。以减少荒地扩张农业的生态修复措施对坝上地区生态系统服务的提升作用相对明显，以林地扩张为主的措施作用则相对有限。

(4) 碳储量、生境质量、水源涵养两两之间均存在正相关关系(协同)，陡坡植树措施会削弱生态系统服务之间的协同作用。