北方重点生态功能区生态系统服务权衡与协同

祝 萍,刘 鑫,郑瑜晗,王世豪,黄 麟

1 中国科学院地理科学与资源研究所生态系统网络观测与模拟重点实验室, 北京 100101

2 中国科学院大学, 北京 100049

3 生态环境部信息中心, 北京 100029

4 中国科学院地理科学与资源研究所区域可持续发展分析与模拟重点实验室, 北京 100101

5 中国科学院地理科学与资源研究所资源与环境信息系统国家重点实验室, 北京 100101

6 中国科学院地理科学与资源研究所陆地表层格局与模拟重点实验室， 北京 100101

北方防沙带是我国“两屏三带”生态安全战略格局的重要组成部分[1],在防风固沙、保持水土、调节气候、农牧产品供给、生物多样性维持等方面都具有重要生态意义[2-3]。然而,北方干旱半干旱区也是我国沙漠分布最广、荒漠化发展最快、受风沙危害最严重的生态极脆弱区,特别是半干旱地区向半湿润地区、高原向丘陵平原、草原牧区向农耕区过渡的农牧交错带[4-7],以及沿内陆河分布或位于内陆河下游的绿洲地区。长期以来,由于其自然生态的脆弱性、区域气候变化叠加超载过牧、草地开垦等人类不合理利用方式,该地区出现了退化、沙化、盐渍化等草地“三化”、沙尘暴灾害、河川断流、湖泊干涸、绿洲消亡等一系列严重生态问题[7-9]。

为了构建并提升国家生态安全屏障的服务功能,按照生态脆弱性和重要性,划定了25个承担涵养水源、保持水土、防风固沙、生物多样性维护等重要生态功能的区域作为国家重点生态功能区[10-11]。自2008年开始,通过实施中央对地方重点生态功能区转移支付政策[12],推进生态文明建设和高质量发展,引导地方政府加强生态环境保护,以提高生态功能重要地区所在地政府的基本公共服务保障能力。截至2020年,对重点生态功能区所涉及近700个县实行转移支付超过5000亿元。转移支付效果如何,实施转移支付区域的生态服务功能是否有所提升,是近年来管理部门和学术界的关注焦点。相关研究聚焦在转移支付政策、考核办法及县域生态环境质量评价[13- 15],生态功能区生态系统结构、核心服务功能变化与人类活动的影响[16-17]、生态承载力[18]、可持续发展能力[19]、乡村人居环境演变[20]等方面。

然而,由于强调重点生态功能区的关键生态功能,忽略了其他同样较为重要的服务功能,进而出现了调节与供给服务之间此消彼长的权衡现象。生态系统服务之间具有相互交织、复杂的非线性关系[21],叠加人类对生态系统服务选择的不同偏好[22],生态系统格局-过程-功能-服务的变化,都会导致生态系统服务之间的权衡/协同关系[23-24]。实际应用中,一般利用数学统计、空间制图、情景模拟、多目标决策和服务流动性分析等方法,开展不同时空尺度下不同服务类型之间的权衡与协同关系研究[25]。近年来,开展生态系统服务权衡分析,认识生态系统服务之间关系,已成为生态管理决策与调控的重要依据,应用于农业、渔业生产[26-27]、森林经营管理[28]、海洋空间规划[29]、能源管理[30]等方面。然而,生态系统服务权衡具有相对复杂的时空尺度,难点在于如何刻画不同尺度生态系统结构、过程、功能、服务变化之间的相互作用及其影响因素[25]。

以草地和荒漠为主的北方草原重要生态功能区,既要通过生态保护和修复保障生态安全,也要通过畜牧业可持续发展保障粮食安全,同步推进生态保护与生产利用是追求的目标,因此,分析过去一段时期生态系统调节与供给服务的时空变化,探究生态系统服务之间的权衡与协同关系,对于解决生态保护与生产利用之间的矛盾具有重要的科学意义。本文以北方农牧交错带和内陆河下游绿洲地区6个草原与荒漠化防治重点生态功能区为例,利用遥感监测、模型模拟、空间分析等手段,分析其2000—2015年主要生态系统类型与防风固沙、土壤保持、水源涵养、牧草供给等主要生态系统服务的时空变化状况,进而揭示重点生态功能区生态系统供给与调节服务之间的权衡与协同关系,为重点生态功能区生态保护与建设提出有针对性的措施建议。

1 研究区概况

本文以塔里木河荒漠化防治生态功能区(简称塔里木河)、阿尔金草原荒漠化防治生态功能区(简称阿尔金)、阴山北麓草原生态功能区(简称阴山北麓)、浑善达克沙漠化防治生态功能区(简称浑善达克)、科尔沁草原生态功能区(简称科尔沁)、呼伦贝尔草原草甸生态功能区(简称呼伦贝尔)作为研究区(图1)。塔里木河是“丝绸之路经济带”核心区和国家西北的重要生态屏障[31],阿尔金是我国西北荒漠区特有野生动植物的生物物种基因库[32],是我国西北的重要生态屏障。阴山北麓、浑善达克、科尔沁、呼伦贝尔地处第二阶梯和第三阶梯的过渡地带,是多种生态系统的交错地带,是京津乃至整个华北地区的重要生态屏障[33- 35]。

图1 研究区位置示意图Fig.1 Spatial distribution of study areas

2 数据与方法

本文首先分析了研究区内生态系统类型和植被覆盖度情况,然后分析了供给型服务(牧草供给)和调节型服务(水土保持、水源涵养、防风固沙)的协同权衡关系。

2.1 生态系统类型及变化数据

基于刘纪远等[36]完成的中国土地利用/覆被变化数据库,生成2000年和2015年全国100m空间分辨率的土地利用/覆被栅格数据集。基于该数据集将生态系统类型划分为森林、草地、湿地、农田、城镇、荒漠6个类型。森林生态系统即木本为主的植物群落,包括自然、半自然植被及集约化经营和管理的人工木本植被。草地生态系统是一年或多年生草本植被为主的植物群落。湿地生态系统包括沼泽、水域、永久性冰川、滩地等。农田生态系统是以收获为目的、有耕犁活动的人工种植草本植物覆盖表面。城镇生态系统即人工建造用于城乡居民点、工矿、交通等的陆地表面。荒漠包括年内最大植被覆盖度小于20%的地表。对各生态系统类型空间数据进行统计分析,综合评价6大生态系统类型。

2.2 植被覆盖度及变化数据

基于2000—2015年MODIS NDVI数据,利用像元二分模型理论计算植被覆盖度。该理论认为一个像元的NDVI值是由绿色植被贡献的信息与无植被覆盖贡献的信息组合而成,根据以下公式计算得到研究区2000—2015年的年最大植被覆盖度数据：

(1)

式中,FVC是植被覆盖度,本文采用5%的置信区间,NDVImax是纯植被像元的NDVI值,此处用95%处的NDVI值,NDVImin是完全无植被覆盖像元的NDVI值,此处用5%处的NDVI值。

采用最小二乘法分析植被覆盖度的年际变化趋势,计算公式为：

(2)

式中,i为2000年到2015年的年序号,slo为植被覆盖度的变化趋势是增加,如果值为正,说明植被覆盖度是增加的,反之则是减少。

2.3 生态系统服务估算方法2.3.1 防风固沙服务量

在充分考虑气候条件、地表土壤的粗糙度、植被状况、土壤可蚀性、土壤结皮的情况下,本文利用修正的土壤风蚀方程(RWEQ)估算区域土壤风蚀模数。

(3)

Qmax=109.8×WF×EF×SCF×K′×COG

(4)

式中,Mw表示土壤风蚀模数,x表示地块长度,Qx表示x处的沙通量(kg/m),Qmax表示风力的最大输沙能力(kg/m),s表示关键地块长度(m)。WF表示气象因子,EF表示土壤可蚀性成分因子,SCF表示土壤结皮因子,K′表示土壤糙度因子,COG表示植被因子。气候因子中的风和土壤湿度根据中国气象科学数据共享服务网(http://cdc.cma.gov.cn)下载的日均风速、降水、温度、日照时数等计算;雪盖因子利用从中国西部环境与生态科学数据中心(http://westdc.westgis.ac.cn)下载的中国雪深长时间序列数据集来计算。土壤可蚀性因子根据Fryear等[37]方程计算。土壤数据来源于西部环境与生态科学数据中心1∶100 万土壤类型图所附的土壤属性表和空间数据。根据文献收集地面测定的我国不同地区的风蚀模数结果[38-41]对本估算结果进行验证。

通过对生态系统防风固沙服务量进行定量分析来衡量生态系统防风固沙的能力。防风固沙服务量为裸土条件和地表覆盖植被条件下的土壤风蚀量的差值:

SLsv=SLs-SLv

(5)

式中,SLsv表示防风固沙服务量,SLs表示裸土条件下的潜在土壤风蚀量,SLv表示植被覆盖条件下的现实土壤风蚀量。

2.3.2水源涵养服务量

通过对比几种水源涵养量估算方法在量级与变化趋势上的差异,采用降水贮存量法[42-43]并改进了各项参数, 估算典型生态功能区草地、林地、湿地生态系统水源涵养量,表示为：

Q=M×J×R0

(6)

J=J0×K0

(7)

式中：Q为与裸地相比较,森林、草地和湿地生态系统涵养水分的增加量(m3),M为生态系统面积(hm2),J为产流降水量(mm),J0为年均降水量(mm),K0为产流降水量占降水总量的比例,R0为与裸地相比较,生态系统减少径流的效益系数。通过搜集已发表文献的实测降雨产流临界值,以临近国家气象台站实测日降水数据修正同时期的TRMM逐日3h降水量数据,累积单次降雨量大于降雨产流临界值的数值得到单点产流降雨量占降雨总量的比例,再与多年平均河川径流系数建立线性关系,得到区域产流降雨量占降雨总量比例的空间分布。森林减少径流的效益系数主要通过已有的文献资料收集得到,草地降雨径流率根据草地植被覆盖度计算得到,不同植被覆盖度[44]下高寒草甸的降水产流特征采用李元寿等[43]的研究结果。

2.3.3土壤保持服务量

采用修正通用水土流失方程(RUSLE)估算区域尺度土壤水蚀模数。通过对生态系统土壤保持量进行定量分析来衡量生态系统保育土壤的能力。土壤保持量为生态系统在极度退化状况下的土壤流失量与现实状况下土壤流失量的差值。

AC=R×K×L×S×CD-CR×P

(8)

式中,AC为土壤保持量,R是降雨侵蚀力因子,K为土壤可蚀性因子,L是坡长因子,S是坡度因子,C是覆盖和管理因子,CD和CR分别表示生态系统在极度退化状况下和现实状况下的土地覆盖和管理因子。P为水土保持措施因子。R采用基于日降雨量资料的半月降雨侵蚀力模型[45]来估算。K采用Nomo图法计算[46],L和S的估算基于McCool等[47]和刘宝元等[48]的方法,C根据蔡崇法[49]的方法计算确定。

2.3.4牧草供给服务量

利用MODIS的植被净初级生产力(NPP)数据产品与草地地下地上生物量比例估算草地产草量[50],得到2000—2015年研究区草地生态系统牧草供给量数据。公式如下：

(9)

式中,Ym为单位面积草地产草量(kg/hm2),r为草地植被地下/地上生物量比例,取值参考马安娜等[51]。t是生物量到生产力的转换系数,取值0.45[50]。

2.4 权衡与协同分析方法

基于逐像元空间相关分析方法,计算两组长时间序列生态系统服务之间的相关系数, 并利用其正负相关性关系判断两种生态系统服务之间的权衡与协同关系。相关系数计算公式如下：

(10)

根据相关系数的零假设检验T检验方法,判断生态系统服务之间相互关系的显著性。T检验公式如下：

(11)

式中,R是对应的偏相关系数,n是样本观测数,此处取值16。当T>T0.05,14时,即P<0.05,拒绝原假设,相关性结果显著;当T>T0.01,14时,即P<0.01,拒绝原假设,相关性结果极显著。T0.05,14和T0.01,14查阅T检验表确定。查阅T值表时,根据经验两种服务间可能存在正相关和负相关,所以用双侧界线。

3 结果与分析

3.1 2000—2015年北方重点生态功能区生态系统面积与植被覆盖度变化

2000—2015年,塔里木河以草地和农田之间的转换为主。草地面积净减少2425.6km2,占2000年草地面积的2.1%。其中,草地减少面积2818.6km2,主要转换为农田;草地增加面积393.0km2,主要来源于湿地。农田面积净增加3678.0km2,占2000年农田面积的33.4%。其中,农田减少面积108.4km2,主要转换为草地和城镇;农田增加面积3786.4km2,主要来源于草地(表1)。

表1 2000—2015年各重点生态功能区生态系统类型变化/km2Table 1 Area changes of ecosystem in key ecological function zones from 2000 to 2015

阿尔金以草地和荒漠之间的转换为主。草地面积净减少337.9km2,占2000年草地面积的0.5%。其中,草地减少面积438.1km2,主要转换为湿地和农田;草地增加面积100.1km2,主要来源于荒漠。荒漠面积净减少687.8km2,占2000年荒漠面积的0.3%。其中,荒漠减少面积为714.8km2,主要转换为湿地、农田和城镇。

阴山北麓以草地和荒漠之间的转换。草地面积净增加308.0km2,占2000年草地面积的0.5%。其中,草地减少面积423.8km2,主要转换为荒漠;草地增加面积863.0km2,主要来源于荒漠。荒漠面积净减少486.1km2,占2000年荒漠面积的2.4%。其中,荒漠减少面积696.6km2,主要转换为草地;荒漠增加面积主要来源于草地,面积为192.5km2。

浑善达克以草地和荒漠之间的转换为主。草地面积净减少384.4km2,占2000年草地面积的0.3%。其中,草地减少面积1541.3km2,主要转换为荒漠;草地增加面积1155.8km2,主要来源于荒漠。荒漠面积净增加197.0km2,占2000年荒漠面积的1.5%。其中,荒漠减少面积988.1km2,主要转换为草地;荒漠增加面积1184.6km2,主要来源于草地。

科尔沁以草地和农田、荒漠之间的转换为主。草地面积净减少1096.9km2,占2000年草地面积的2.1%。其中,草地减少面积1884.8km2,主要转换为农田和荒漠;草地增加面积787.9km2,主要来源于农田和荒漠。农田面积净增加1188.2km2,占2000年农田面积的4.5%。其中,农田减少面积427.6km2,主要转换为草地;农田增加面积1615.6km2,主要来源于草地。

呼伦贝尔以草地、荒漠、湿地之间的转换为主。草地面积净增加1.8km2,其中,草地增加面积90.3km2,主要来源于湿地和农田,草地减少面积88.5km2,主要转换为荒漠和湿地。湿地面积净减少108.4km2,占2000年湿地面积的2.4%。其中,湿地减少面积135.6km2,主要转换为荒漠和草地;湿地增加面积为27.23 km2,主要来源于草地。荒漠面积净增加155.8km2,占2000年荒漠面积的12.9%。其中,荒漠减少面积20.1km2,主要转换为草地;荒漠增加面积174.8km2,主要来源于湿地。

2000—2010年与2010—2015年的植被覆盖度变化对比来看,各区域多年平均植被覆盖度均有所提高。其中,覆盖度最高的呼伦贝尔和科尔沁增加最多,分别提升了7.0%和3.1%,覆盖度最低的阿尔金则提高了0.4%。

表2 2000—2015年各重点生态功能区植被覆盖度及其变化统计Table 2 Vegetation coverage of key ecological function zones from 2000 to 2015

3.2 2000—2015年北方重点生态功能区生态系统服务动态变化

3.2.1防风固沙量变化

2000—2015年,除呼伦贝尔的防风固沙量稳步增加(年增幅超过2.5 t/hm2)以外,其余各功能区防风固沙量整体呈减少趋势(图2),其中,阴山北麓、浑善达克和科尔沁的大部分区域年减幅超过2.5t/hm2;塔里木河、阿尔金的年减幅则低于2.5 t/hm2。

图2 2000—2015年重点生态功能区防风固沙服务量变化空间分布Fig.2 Spatial variations of sand fixation services in key ecological function zones from 2000 to 2015

3.2.2水源涵养量变化

2000—2015年,各生态功能区水源涵养量呈现总体稳定、局部减少态势(图3)。浑善达克与呼伦贝尔的水源涵养量年增幅约0.03m3/km2,水源涵养总量分别增加了0.5万m3和0.1万m3。然而,塔里木河东北部、阴山北麓西部、浑善达克东南部轻微下降,同时,科尔沁与呼伦贝尔的湿地水源涵养呈现下降趋势。

图3 2000—2015年重点生态功能区水源涵养服务量变化空间分布Fig.3 Spatial variations of water conservation services in key ecological function zones from 2000 to 2015

3.2.3土壤保持量变化

2000—2015年,各生态功能区的土壤保持服务呈现整体好转、局部恶化态势(图4)。其中,塔里木河南部、阿尔金中部、阴山北麓中部、浑善达克东南部、科尔沁和呼伦贝尔的土壤保持量提升明显。相反地,塔里木河东北部、阿尔金西南部、阴山北麓西部等区域呈现轻微下降,每年降幅约1—2 t/hm2。

图4 2000—2015年重点生态功能区水土保持服务量变化空间分布Fig.4 Spatial variations of soil conservation services in key ecological function zones from 2000 to 2015

3.2.4牧草供给量变化

2000—2015年,除呼伦贝尔西部牧草供给有所增长以外,各生态功能区的牧草供给量都呈现较为明显的下降趋势(图5),特别是科尔沁和浑善达克,单位面积产草量分别减少1.9 t/km2和1.2 t/km2,阿尔金由于单位面积产草量低,年减幅约0.4 t/km2。

图5 2000—2015年重点生态功能区牧草供给服务量变化空间分布Fig.5 Spatial variations of foliage supply in key ecological function zones from 2000 to 2015

3.3 北方重点生态功能区生态系统服务权衡与协同关系分析3.3.1 牧草供给与防风固沙的相关关系

2000—2015年,各生态功能区牧草供给与防风固沙之间的相关关系如图6所示,塔里木河大部分区域相关关系不显著;科尔沁大部分区域、阿尔金东北部、阴山北麓西部、浑善达克中部皆为显著协同关系;相反地,呼伦贝尔中部和东部以显著权衡关系为主。

图6 重点生态功能区牧草供给与防风固沙之间的相关关系空间分布Fig.6 Tradeoffs and synergies between foliage supply and sand fixation in key ecological function zones

3.3.2牧草供给与水源涵养的相关关系

2000—2015年,各生态功能区牧草供给与水源涵养之间的相关关系如图7所示,塔里木河中部、阿尔金中部呈现显著协同关系;阴山北麓东、西部差异明显,西部相关关系不显著,而东部为极显著协同关系;浑善达克中部和西部、呼伦贝尔中部和西部皆为极显著协同关系,而科尔沁北部和西南部、呼伦贝尔东端表现为为显著权衡关系。

图7 重点生态功能区牧草供给与水源涵养之间的相关关系空间分布Fig.7 Tradeoffs and synergies between foliage supply and water regulation in key ecological function zones

3.3.3牧草供给与土壤保持的相关关系

2000—2015年,各生态功能区牧草供给与土壤保持相关关系如图8所示,塔里木河和阿尔金大部分区域相关关系不显著;阴山北麓西部相关关系不显著,而东部为极显著协同关系;浑善达克东部和东南部为显著权衡关系,而中部和北部以显著协同关系为主;科尔沁北部和西部、呼伦贝尔东端为极显著权衡关系,而科尔沁东部和呼伦贝尔中、西部皆为极显著协同关系。

图8 重点生态功能区牧草供给与土壤保持之间的相关关系空间分布Fig.8 Tradeoffs and synergies between foliage supply and soil conservation in key ecological function zones

4 结论

本研究得到如下几点结论：2000—2015年,除阴山北麓以外,其他功能区开垦导致草地面积皆减少,东部三个功能区荒漠面积增加而西部三个功能区有所减少。草地植被覆盖度皆有所提高,增幅0.4%—7%。除呼伦贝尔以外,其他功能区的防风固沙服务皆下降;西部塔里木河与阿尔金的水源涵养服务先上升后下降,中部阴山北麓和科尔沁有所提高,而东部浑善达克与呼伦贝尔持续上升;各功能区土壤保持服务呈现不同程度提升,而牧草供给服务皆有所下降。除呼伦贝尔和浑善达克南端以外,牧草供给与防风固沙服务之间以协同关系为主; 除科尔沁西北部和呼伦贝尔东端以外,牧草供给与水源涵养、土壤保持服务之间以协同关系为主。通过生态系统类型变化及供给与调节服务相关关系分析可以看出,生态功能区开垦导致草地减少的现象需要遏制,出现权衡关系的东部区域也是生态保护与生产利用矛盾较大的农牧交错带,需要成为功能区后续保护和修复的重点区域。

5 讨论

本研究将长时间序列生态系统服务量化数据的两两正负相关性转化为其协同和权衡关系,分析生态系统供给与调节服务之间此消彼长或相互增益的变化趋势,对于如何全面提升重点生态功能区生态系统管理具有指导意义。初步结论表明,除呼伦贝尔和浑善达克南端以外,牧草供给与防风固沙服务之间以协同关系为主; 除科尔沁西北部和呼伦贝尔东端以外,牧草供给与水源涵养、土壤保持服务之间以协同关系为主。本研究的不确定性存在以下几个方面：1)有研究表明：人类活动对生态系统的选择性干预引起的生态系统服务变化是权衡与协同,而自然因素引起的生态系统服务之间的此消彼长,是一种竞争而非权衡关系[52],在权衡关系的研究中还需进一步区分人类活动和自然因素。2)生态系统提供多重服务,并且各种服务之间相互作用、相互联系、相互交织[53-54],本研究仅仅分析两两间的相互作用是不充分的。3)还应该考虑不同态系统服务对外界影响的响应时间和空间尺度不同,不同的生态系统服务之间的权衡协同关系,也可能发生于现在和未来的提供之间。有研究证明,在密西西比河谷大量施肥带来生产力增加,在农业增产后的20年,墨西哥湾才出现了死亡地带[55]。因此,需要深入量化产生协同和权衡关系的原因,在两两关系基础上开展多维生态系统服务相关关系分析以及通过权衡与协同关系的情景预测把握生态保护与管理措施的功能提升效果。