卧龙自然保护区典型生态系统服务时空变化研究

曹梦琪,蔡英楠,张 丽,徐建英

首都师范大学资源环境与旅游学院,三维信息获取与应用教育部重点实验室, 北京 100048

生态系统服务(Ecosystem Services, ESs)是指生态系统形成和提供给人们生存发展的直接产品或间接惠益。根据联合国千年生态系统评估,生态系统服务分为供应、调节、文化和支持服务4种类型[1]。一定时间和空间尺度上的生物与非生物特征及其生态过程构成了生态系统结构与功能的差异,成为了生态系统服务产生的基础[2]。不同区域的地形地貌特征、土地利用变化、气候变化以及人类活动等因素会改变生态系统的结构与空间特征,进而影响生态系统服务的时空分布差异和变化[3- 9]。全面了解特定时空下生态系统服务的时空差异,可为生态系统服务保护、管理和优化提供依据,促进自然资源的管理保护以及针对性的生态恢复[10- 14]。

生态系统服务的时空变化是当前的研究热点[15- 17],国内外学者已进行了大量的研究,业已形成相对成熟便捷的数量评估模型,如CASA、RUSLE和InVEST等,研究主题包括生态系统服务时空分布和变化趋势[18- 19]、时空格局和热点区域识别[20]以及生态系统服务变化的驱动因素分析[21-24]等方面。鉴于生态系统服务时空变化分析对于生态和社会经济系统的显著意义,这些研究往往关注生态环境比较脆弱的区域或人类干扰比较显著的区域,如王晓峰等研究了新疆地区3种生态系统服务的分区差异[18];潘竟虎等以干旱内陆河流域典型地区嘉峪关-酒泉地区为例,定量分析了2000和2010年4种生态系统服务的时空变化[19];许丁雪等研究了张家口-承德地区土地利用变化对生态系统服务的驱动作用[21],上述研究揭示我国生态环境脆弱区生态系统服务的时空变化并为区域生态恢复和建设提供了借鉴。而冉凤维[20]等、María Anaya-Romero[22]等、苏常红[23]等、潘梅[24]等学者分别研究了人口集中、人类干扰比较显著的鄱阳湖地区、地中海地区、汾河流域和京津冀地区,强调了人类活动对于生态系统服务变化的影响,分析了生态系统服务变化的驱动力。自然保护区是生态系统保护的关键地区和生态系统服务供给的重点区域,目前的相关研究主要集中于对各种保护项目和措施的成效分析[25- 26],少有针对保护区典型生态系统服务的定量研究[27]。同其他区域相比,保护区提供生态系统服务的能力尤为显著[28],是我国生态文明建设的重点,也是我国诸多生态保护/建设项目实施的重点区域。因此,对自然保护区生态系统服务进行定量研究,结合区域生态系统特征和生态保护/建设项目探讨分析其时空变化特征以及自然和人为驱动因素,可以为自然保护区生态系统保护、恢复和管理提供决策依据和建议。

本研究拟以卧龙自然保护区为例进行研究。卧龙自然保护区位于四川省汶川县,主要保护西南高山林区自然生态系统及大熊猫等珍稀动物,因大熊猫的保护在国内外著称,隶属我国主要的固碳贡献区——西南林地[29]。此外,保护区位于岷江上游,是岷江流域重要的水源涵养地,其水源涵养和土壤保持功能对维持区域生态安全有重要作用。因此,针对卧龙自然保护区保护对象及其承担的生态功能,本研究选择碳固定(NPP)、土壤保持、水源涵养和生境质量4种典型生态系统服务进行定量评估,分析主要生态建设工程实施以来保护区典型生态系统服务的时空分布和变化特征,确定生态系统服务供给的热点区域,以为生态保护/建设政策的实施和保护区的生态系统管理提供建议。

1 研究方法与数据来源

1.1 研究区概况

卧龙自然保护区(东经102°52′—103°25′,北纬30°45′—31°25′)位于四川省阿坝藏族羌族自治州汶川县境内,四川盆地的西缘,岷江的上游,是青藏高原向成都平原的过渡地带,地势西北高、东南低,南北长约62 km,东西宽约52 km,总面积约为2000 km2(图1)。保护区属典型的亚热内陆山地气候,年温差小,干湿季节分明,植被和土壤随海拔变化的垂直分布规律明显。区内动植物资源丰富,拥有近4000种植物以及大熊猫、金丝猴等多种珍稀濒危动物。水能资源丰富,生态系统类型复杂多样,是生态系统服务研究的典型区域。

图1 研究区位置及地形示意图Fig.1 Location and elevation of the study area

1.2 数据来源及预处理

本文使用的基础数据包括：土壤数据、植被数据、气象数据、高程数据和土地利用类型数据(表1)。研究时间为2000年和2015年两个年份的5—10月,该时段降水丰沛,雨热同期,是植被生长的最佳季节。

表1 数据来源及预处理

1.3 研究方法

1.3.1碳固定(NPP)

本研究碳固定服务用植被净初级生产力(Net Primary Productivity,NPP)予以指代。采用朱文泉[30]改进的CASA(Carnegie-Ames-Stanford Approach)模型估算研究区2000年和2015年的NPP,其计算公式如下：

NPP=APAR×ε

(1)

APAR=SOL×FPAR×0.5

(2)

ε=Tε1×Tε2×Wε×εmax

(3)

式中,APAR表示植被吸收光合有效辐射(MJ/m2);ε表示实际光能利用率(gC/MJ);SOL表示太阳总辐射量(MJ/m2),FPAR 表示植被层对入射的光合有效辐射(PAR)的吸收率,可通过归一化植被指数(NDVI)对其进行估算;常数0.5表示植被所利用的太阳有效辐射占太阳总辐射的比例;Tε1和Tε2表示温度对光能利用率的影响;Wε表示水分条件对光能利用率的影响;εmax表示在理想状态下植被的最大光能利用率(参考朱文泉等的研究成果)。

1.3.2土壤保持

修正通用土壤流失方程(Revised Universal Soil Loss Equation,RUSLE)经众多学者采用来评估土壤保持量,均取得较好效果[31- 33],计算公式如下：

Am=R×K×LS×C×P

(4)

Ap=R×K×LS

(5)

Ac=Ap-Am

(6)

式中,Am为实际土壤侵蚀量(t hm-2a-1),指当前地表覆盖和水土保持措施下的土壤侵蚀量;Ap为潜在土壤侵蚀量(t hm-2a-1),指在没有植被覆盖和水土保持措施情况下可能产生的土壤侵蚀量(C=1,P=1);Ac为土壤保持量(t hm-2a-1);R为降雨侵蚀因子(MJ mm hm-2h-1a-1);K为土壤侵蚀因子(t hm2h MJ-1hm-2mm-1);LS为坡长坡度因子;C为植被覆盖因子;P为水土保持措施因子。

(1)降雨侵蚀因子R采用 Wischmeier 等提出的基于月尺度降雨量数据的经验计算公式[34- 35]：

(7)

式中,pi为第i个月的降水量(mm),p为年降水量(mm)。

(2)土壤侵蚀因子K采用EPIC模型进行估算,计算公式如下[36]：

K=(-0.01383+0.51575×KEPIC)×0.1317

(8)

式中,ms为砂粒含量;msilt为粉粒含量;mc为粘粒含量;orgC为有机质含量。

(3)坡长坡度因子LS的计算公式如下[32]：

L=(λ/22.1)m

(9)

(10)

(11)

式中,λ为坡长(m);m为坡长效应指数;θ为坡度(°)。

(4)植被覆盖因子C的计算公式如下[37]：

(12)

(13)

式中,f为植被覆盖度;NDVI为植被归一化指数,NDVIsoil和NDVIveg分别为土壤和植被的NDVI。

(5)水土保持措施因子P指在采取保持措施情况下与顺坡耕作时土壤流失量的比值。P=0表示不发生侵蚀,P=1表示未采取任何水土保持措施。由于缺乏相关数据,本研究不考虑P因子影响,P值定义为1。

1.3.3生境质量

生境质量指生态系统为物种生存繁衍提供的条件或潜力[38]。本文中特指保护区生态系统可为大熊猫生存繁衍所利用的环境条件状况。InVEST模型生境质量模块以土地覆被为基础,利用外界威胁强度(影响的距离和权重)、生境适宜性、生境类型对威胁的敏感性进行生境质量评估。参照以往研究成果[39-41]和保护区实际情况,本研究选取耕地、主要道路、乡村小路和城镇居民用地作为威胁源(表2)并对各土地利用类型的生境适宜性进行打分以及对威胁源敏感性赋值(表3)。生境适宜性得分为0—1,其中1表示该生境具有最高适宜性,而非生境取值为0。计算公式如下：

表2 威胁源的最大影响距离及权重

表3 土地利用类型对威胁源的敏感性参数

(14)

(15)

式中,Qxj是土地利用类型j中栅格x的生境质量;Dxj是土地利用类型j中栅格x的生境胁迫水平;k为计算的半饱和参数;Hj为土地利用类型j的生境适宜性;R为胁迫因子;Yr为胁迫因子r所占的栅格数;Wr为胁迫因子权重,值为0—1;ry为栅格y的胁迫因子值;irxy为ry对栅格x的胁迫水平;βx为栅格x的可达性水平;Sjr为土地利用类型j对胁迫因子r的敏感性。

基于以上模型计算结果,根据欧阳志云[39]等提出的大熊猫生境自然环境因素评价准则中的适宜海拔和坡度对结果进行进一步修正,将研究区海拔和坡度分为适宜和不适宜两部分,与模型结果进行叠加。

1.3.4水源涵养

计算水源涵养,首先通过水量平衡原理得到保护区产水量,再根据土壤深度、土壤饱和导水率、地形指数和流速系数等对产水量进行修正得到水源涵养量,计算公式如下[42- 43]：

Y=P-ET

(16)

WR=min(1,249/V)×min(1,0.9×D/3)×min(1,Ksoil/300)×Y

(17)

(18)

式中,Y为产水量(mm);P为降水量(mm);ET为蒸发量(mm);WR为多年平均水源涵养量(mm);V为流速系数(采用模型参数表数据[44]);D为地形指数,无量纲;Ksoil为土壤饱和导水率(cm/d);Wpc(Watershed pixel count)为集水区栅格数,无量纲;Sd(Soil depth)为土壤深度(mm);Ps(Percent slope)为百分比坡度。由于缺乏2000年ET数据,本研究使用2001年数据代替。

1.3.5地形位指数

地形位指数是复合分析海拔和坡度属性信息的指标,能够综合反映空间内某点的地形条件。其计算公式如下[45]：

(19)

式中,T为地形位指数,E及E0分别为空间内任一栅格的高程(m)和平均高程(m),S和S0分别为空间内任一栅格的坡度值(°)和平均坡度值(°)。一般地,海拔低、坡度小的栅格地形位指数小,反之越大;海拔低坡度大或海拔高坡度小的栅格,其地形位指数居中。本研究根据计算结果采用自然断点法将研究区地形位指数划分为10个梯度(表4)。

表4 地形位指数分级标准

1.3.6热点分析方法

生态系统服务热点被定义为具有高服务多样性、高生物物理或货币服务价值或高服务能力的区域[18]。热点区域的识别,可进一步分析生态系统服务供给能力的强弱。本研究将保护区内各生态系统服务超过各自平均值的栅格定义为该类服务的热点区域,将这四类生态系统服务的热点区进行叠加分析,即可得到卧龙自然保护区多重生态系统服务热点区的分布状况。在本研究中,若某栅格单元内4种生态系统服务的值均未超过各自平均值,则将其定义为非热点区,若仅有一种生态系统服务的值超过平均值,则将其定义为I类热点区,以此类推,可得到II类、III类、Ⅳ类热点区[18-20]。

2 结果分析

2.1 生态系统服务的时空变化

从时间上看,2000—2015年,NPP、土壤保持、生境质量和水源涵养4种服务均呈现增加趋势,该增加趋势在不同类型生态系统服务的均值、最高值均表现明显,且4种服务的增益(即某生态系统服务2015年比2000年增加的区域)面积均显著大于减损(即某生态系统服务2015年比2000年减少的区域)面积(表5)。以NPP为例,其平均值由2000年的438.95 gC/m2增加到2015年的509.94 gC/m2,最高值由2000年的1125.32 gC/m2增加到2015年的1340.69 gC/m2,增益面积占保护区的86.3%而减损面积仅占13.7%。

表5 2000—2015年卧龙自然保护区生态系统服务变化情况

从空间上来看,4种典型生态系统服务从2000年到2015年均保持了相对稳定的空间格局以及增益面积超过减损面积的良好趋势,但具体的空间变化区域有所不同(图2)。NPP与土壤保持呈现近似的空间分布格局,即西北—东南方向上由低到高分布,增益面积主要分布在低海拔区,而减损面积主要分布在高海拔区。生境质量指数呈现东高西低的空间分布特点,高值区主要位于东部,北部和西部广大地区以及南部边缘地区生境质量指数为零,东南地区增益明显。水源涵养服务的空间分布无明显规律,其高值零星分布在北部、东部和南部边缘以及中部部分地区,低值区不同程度地分布在全域范围内,与生境质量相同的是,水源涵养的增益区也在东南区域表现显著。

图2 卧龙自然保护区2000—2015年生态系统服务及其变化的空间分布Fig.2 The spatial distribution of four ESs and their changes in Wolong Nature Reserve from 2000 to 2015NPP：植被净初级生产力 Net Primary Productivity;ESs：生态系统服务 Ecosystem Services

2.2 不同地形位梯度上生态系统服务的变化

从不同地形位梯度上各生态系统服务均值的分布来看,NPP、土壤保持、生境质量和水源涵养4种服务在2000年和2015年均各自保持相同的趋势,且2015年各服务在不同地形位梯度上的均值均大于2000年,但在具体地形位梯度上的分布有所不同(图3)。NPP和生境质量均呈现随地形位指数增大而递减的趋势,其高值多分布在较低地形位梯度即海拔较低且坡度较小的地区;土壤保持随地形位指数增大呈现先增高后降低的趋势,其高值主要分布在5—8级地形位梯度上;水源涵养呈现较为起伏波动的趋势,其高值主要分布在较低地形位梯度上,最高值分布在2级地形位上,最低值分布在9级地形位上。

从不同地形位梯度上生态系统服务变化类型所占的面积来看(图3),研究区内NPP、土壤保持和水源涵养3种服务的变化类型均以增益为主,随地形位指数的增大其增益面积比重均呈降低趋势,减损和不变的面积比重均呈递增趋势,其中水源涵养的减损类型面积比增加相对显著,在10级地形位梯度上甚至超过了50%。生境质量在区域内大部分地区表现为不变,其增益和减损面积比均在低地形位上相对显著,并随地形位指数增大而逐渐降低,因海拔高于3750 m、坡度45°以上的区域不是大熊猫生境适宜区,所以高地形位梯度上生境质量没有发生变化。由此来看,低地形位指数地区是生态系统服务增益优势区,海拔越高、坡度越大的区域减损面积相对增加。

图3 不同地形位梯度上的生态系统服务均值及其变化类型面积比Fig.3 The mean value of four ESs and their ratio of variation types at different terrain niche

2.3 不同土地利用类型上生态系统服务的差异

不同土地利用类型上生态系统服务供给量存在明显差异(图4)。其中林地和草地既是研究区主要的土地利用类型,也是生态系统服务供给的主体,其生态系统服务供给量占保护区总量的95%以上。从时间来看,2000年不同土地利用类型上4种生态系统服务的供给量均表现为：林地>草地>耕地>水域>建设用地>未利用地,其中NPP、土壤保持、水源涵养3种服务在林地上供给量分别为各自总量的49.9%、57.4%和50.5%,在草地上的供给量分别为其各自总量的48.7%、41.5%和48.9%。2015年不同土地利用类型上4种生态系统服务的供给量均表现为：林地>草地>水域>耕地>建设用地>未利用地。林地和草地依然是贡献较大的土地利用类型,但生态系统服务的供给量有所变化,其中林地的生态系统服务供给量有所提升,其提供的NPP、土壤保持、水源涵养3种服务均达到各自总量的59%以上,而草地生态系统服务的供给量有所减少,相应3种服务的供给量均低于各自总量的40%,低于2000年。此外水域提供的生态系统服务比重超过了耕地。

图4 2000年和2015年不同土地利用类型上生态系统服务供给量比Fig.4 The ratio of ecosystem service supply on different land use types in 2000 and 2015

2.4 热点区分析

研究区域不同类型热点区的面积变化分析表明(表6),尽管2000年和2015年四类热点区保持了相同的数量排序,即I类热点区>II类热点区>III类热点区>Ⅳ类热点区,但是与2000年相比,2015年Ⅳ类热点区面积增加显著,增加了134.59 km2,占保护区总面积的6.74%。非热点区的面积也有所增加,增加面积为34.97 km2。这两类热点类型增加的面积均来自其相邻类型的热点区,即Ⅳ类热点区增加的面积主要来自III类热点区,而非热点区增加的面积主要来自I类热点区(表7),II类热点区则主要转为III类热点区(104.66 km2),其次为I类热点区(62.58 km2)。生态系统服务热点区的变化表明,研究区域整体生态环境改善明显,但是局部区域生态环境有所退化,且生态环境改善和退化的区域具有渐变特征,因此需要加强I类和III类热点区的管理,预防I类热点区的生态环境退化,同时促进III类热点区向Ⅳ类热点区转化,实现区域生态系统服务供给的协同和最优化。

表6 研究区域不同类型热点区面积/km2

表7 研究区域生态系统服务热点区转移矩阵/km2

从空间分布来看,生态系统服务热点区及其变化具有明显的空间特征(图5)。研究区域的东部和中部地区Ⅳ类热点区面积明显高于西部,特别是西北部区域,呈现由东南向西北方向由强变弱的特征。东南部热点区增加明显,以Ⅳ类热点区最为显著,此外中部地区,特别是河流(皮条河)两侧Ⅳ类热点区的增加尤为显著。

图5 卧龙自然保护区2000和2015年生态系统服务热点区分布Fig.5 Hotspots of ESs in Wolong Nature Reserve in 2000 and 2015

3 讨论

生态系统服务的时空变化是自然和社会经济因素共同驱动的结果,总体来看,2000—2015年研究区域目标生态系统服务的供给能力明显加强,生态系统服务之间的协同程度有所提高,一定程度上表明研究区域生态环境在持续改善和提高,但是部分区域脆弱生态环境依旧存在,空间异质性特征明显。研究区域的东部、东南和中部地区目标生态系统服务供给量和生态改善特征明显,而南部边缘、西部和西北部生态脆弱性特征明显。研究区域目标生态系统服务的空间异质性及其变化与区域地形、土地利用、人类活动以及生态措施等密切相关。

研究区域目标生态系统服务的分布在不同地形位梯度上分异明显,除土壤保持外,其他生态系统服务的高值分布区和主要增益区均以低海拔低坡度区域为主,这可能与两个因素有关：一是研究区域为典型的高山峡谷区,气候和植被垂直分异特征明显,高海拔林线以上发育高山草甸、冰川雪地和裸地(未利用土地),而低海拔和低坡度区域的水热条件更有利于植被的生长和发育;其次高山峡谷区土地资源紧张,低海拔低坡度区域也是人类活动的主要分布区,居民点和耕地资源集中分布,但2000年前后研究区域大力推进退耕还林,部分耕地转化为林地,近十五年的生态恢复效果明显。作为我国生态恢复和建设的重大举措,退耕还林/还草政策的生态效应在各区域正在逐步显现[21,23],是近期我国生态改善的主要驱动力之一。

土地利用类型及变化作为影响生态系统服务时空变化的关键因子,通过直接或间接影响生态系统的格局与过程,改变生态系统服务的提供能力[46]。林地作为生态系统服务的供给量较高的土地利用类型,是研究区域目标生态系统服务增加和协同关系提高的主力,也是未来生态改善措施首选的土地利用类型。尽管草地的生态系统服务贡献较为突出,但是近年来其供给能力有所下降。研究区域的草地主要为林线以上的高山草甸,是在低温大风及瘠薄土壤等生态环境条件下发育的植被类型。草地生态系统服务供给能力下降,可能和研究区域人类活动有关,高山草甸是当地牧业(牦牛等)和众多中草药资源主要分布场所,人类利用程度的增加和干扰增强,可能是高山草甸生态系统服务供给力降低的主要原因。因此在生态保护和恢复工作中,除了目标物种(大熊猫)及其生境(林地)外,应该加强草地生态系统的保护和恢复,避免该脆弱生态系统的退化和过度利用。

研究区域生态系统的保护和恢复既需要考虑时空异质性,又要考虑自然和社会经济驱动力。从研究结果来看,以下三类区域应采取措施,加强保护、恢复或管理：一是以Ⅳ类热点区为代表的生态保护区,特别是东南部区域,应该加强生态环境的保护;其次是以I类热点区和III类热点区为代表的生态恢复区,预防生态环境的退化,促使其向多种生态系统服务协同的方向转化;三是以低海拔低梯度区域和高海拔草地为代表的生态协调区,人类活动对于该区域生态环境的恢复和退化有决定性作用,建议合理引导或控制人类活动,维护已经取得的生态恢复成果,防止脆弱生态系统的进一步退化。

4 结论

本研究以卧龙自然保护区为例,定量评估了2000年和2015年碳固定、土壤保持、生境质量和水源涵养4种典型生态系统服务的时空分布格局及变化特征,分析了生态系统服务变化的自然和社会经济驱动力,提出了分区管理的建议,对保护区的生态保护、恢复和管理有一定的借鉴作用。研究结果揭示了生态系统服务变化的时空异质性,表明在一定的区域内生态环境改善和退化可能并存,生态保护和恢复须因地制宜,甄别关键区域,制定针对性的政策。其次生态系统服务的变化既受自然因素的控制,又受人为因素的影响,生态恢复/建设工程不仅要关注生态系统服务供给能力较高的生态系统类型(如林地),也要关注人类活动干扰较大的相对脆弱的生态系统类型(如研究区域的草地),保护区域生态系统的多样性。

本文也存在如下不足：一是由于数据获取困难,本研究仅选取了4种典型生态系统服务,未来需增加生态系统服务类型,特别是与保护区生态保护和人类活动相关的生态系统服务类型;二是需要延长研究年份,将生态系统服务的连续变化与土地利用类型的连续变化相联系,定量评估生态建设工程对于生态系统服务变化的贡献;三是由于篇幅所限,仅考虑了4种典型生态系统服务的时空变化,没有考虑生态系统服务之间的权衡/协同关系,特别是生境质量与其他类型生态系统服务的权衡/协同关系。在以后工作中,须逐步增加数据的获取分析,拓展并深化研究主题,为提升生态系统服务供给水平,促进生态保护和恢复提供借鉴。