生态系统服务不仅为人类提供了生产生活必需品, 还维持了地球生命支持系统, 如保持物种多样性, 净化环境, 保持大气平衡与稳定等。长期以来由于人们对生态系统服务重要性缺乏认识, 对生态系统的过度使用导致了生态系统服务的退化, 对人类的安全与健康造成影响, 并影响到区域和全球生态安全^[1]。随着“3S”空间分析技术的发展, 基于模型的区域生态系统服务研究得以实现^[2]; 特别是随着多重生态系统服务综合研究的进展, 生态系统服务评估模型集成技术也逐渐发展起来(如GUMBO、MIME和InVEST等); 斯坦福大学等机构联合开发了InVEST模型, 通过预测不同土地利用情景下生态系统服务的变化, 实现了对多重生态系统服务集成和权衡分析。

生态系统服务及其变化是一个由多种因素所驱动的复杂过程; 这些驱动因素涉及自然、社会经济、以及人为干扰等多个方面^[3-4]。探讨生态系统服务变化及其驱动因素, 有助于认识生态系统服务形成机理, 实现生态系统服务的可控化, 在制定生态保护政策时发挥其驱动与反馈机制, 实现社会经济可持续发展。生态系统服务驱动因素中, 自然因素相对稳定, 具有累积效应; 而人文驱动因素尤其是社会经济驱动力相对更活跃, 其包含的因素也更为广泛复杂, 存在着更多的不确定性且难以定量化。生态系统服务驱动机制分析大致分定性与定量两类方法; de Lima等通过历史文献挖掘, 定性分析了巴西大西洋沿岸森林景观格局变化及对生态系统服务提供能力的影响, 表明甘蔗/咖啡轮作、道路建设、贵重金属采掘是驱动当地生态系统服务的主要因素, 从历史来看, 经济活动、乡村建设、城市化、旅游共同导致了区域生态系统服务退化^[5]。Wang等系统分析了退耕还林/草政策下自然环境(降雨、温度、坡度)、社会经济(人口、耕地、收入、一二产比例、粮食产量、园艺作物产量)、政策(政策监管区域、道路、粮食补贴、节水渠)等要素对生态系统服务的驱动作用, 并采用主成分方法(PCA)构建了生态系统服务驱动模型^[6]; 唐秀美等采用STIRPAT模型和地理加权回归模型, 分析了北京市生态系统服务驱动因素, 表明人口、GDP、绿化率、三产、GDP能耗、城市化率是影响北京生态系统服务的主要因素^[7]。此外, 还有研究分析不同因素对生态系统服务影响进程中“自然因素”和“人文因素”的转折点, 对自然与人文影响按时序进行“剥离”^[8]。

在生态系统服务诸多驱动因素中, 人类活动定量化研究存在着较多的不确定性, 尤其是政策、经济、文化等因素难以定量化。人类活动包含了直接因子(如土地利用变化、污染等)和间接因子(如人口、政治、经济、文化、科技等)。早期的人类活动强度研究多用单一的因素如人口、农业生产、资源开采、放牧等来表征人类活动强度; 联合国开发署(UNDP)于1990年提出一个综合性的人类发展指数(HDI)涵盖了人口密度、农田比率、自然资源禀赋、经济发展、城市化、政策等多个因素^[9]。此基础上, 综合自然、经济、社会等多方面因子构建人类活动强度指标体系的研究逐渐开展起来^[10-11], 但这些人类活动评价多缺乏空间定量化分析。胡志斌等^[12]在ArcGIS软件支持下, 基于地理学自相关理论, 选择居民点、道路、地形等因子, 对岷江上游人类活动强度实现了定量化空间表达。

本研究基于前人研究, 采用InVEST、CASA等生态模型对汾河上游流域生态系统服务进行评估, 搜集区域自然和社会经济因素对生态系统服务驱动机制进行定量分析; 其中人文驱动因素的空间分析参照胡志斌等^[12]的方法, 综合居民点、道路、人口等信息对区域人文驱动因素进行空间刻画。通过研究, 为汾河上游流域生态保护规划和社会经济可持续发展提供决策依据。

1 研究区概况

研究区位于管涔山至太原上兰村以北的汾河上游段, 行政区域涉及宁武、静乐、岚县、娄烦、古交、阳曲及太原市尖草坪区、万柏林区部分区域(图 1), 共计89个乡镇, 约70万人口; 地理位置为111°21′—112°27′E, 37°51′—38°59′N之间, 流域长126 km, 宽93 km, 面积5253.56 km²。其中耕地面积2207 km²、林地1822 km²。地貌单元由丘陵沟壑、土石山区和河川阶地构成, 地质环境复杂, 水土流失严重, 流域年径流量6.87亿m³; 丘陵沟壑区年均侵蚀模数6080 t/km², 属强度侵蚀, 也是汾河的主要沙源地。流域气候特征属温带大陆性季风气候, 干旱少雨, 日照充足, 湿热同季。本区工业基础薄弱, 支柱产业为种植业与养殖业。

2 研究方法 2.1 生态系统服务类型的选择

针对汾河上游流域水土流失严重、区域工农业及生活用水的严峻形势, 资源开采与毁林开荒造成的植被破坏、黄土高原气候暖干化大背景、以及退耕还林等生态保护政策对粮食安全的影响等, 本研究筛选了泥沙截持、产水量、NPP生产、固碳释氧、粮食生产五类生态系统服务进行研究。

2.2 生态系统服务评价方法 2.2.1 泥沙截持

采用InVEST模型下的RUSLE模块对土壤保持服务进行评估:

式中, A为单位面积减少土壤流失量(t/hm²), R为降雨因子, 单位(MJ mm hm^-2 h^-1), K为土壤可蚀性因子(t h MJ^-1 mm^-1), LS为地形因子, 包括坡长因子和地面坡度因子, C为植被覆盖因子, P为土壤保持措施因子。本研究参照Fu等^[13]在黄土高原的研究成果, 制订如下C值表(表 1)。

2.2.2 产水服务

采用InVEST模型中water yield模块, 基于budyko曲线和年均降雨量计算产水量^[14], 计算公式：

式中, Y_xj指j类土地利用类型的x栅格的产水量, AET_xj指j类土地利用类型的x栅格的实际蒸散发, P_x指x栅格每年实际降雨量。

2.2.3 NPP

NPP是基于CASA(Carnegie Ames Stanford Approach)过程模型来计算的; 该模型综合考虑植物吸收的光合有效辐射与光能转换率, 同时考虑到光能利用率、土壤水分、降水量、平均温度等因素的影响^[15]。

式中, t表示时间, x表示空间位置, NPP_{(x, t)}表示像元x在t时间的植被净初级生产力(g C m^-2月^-1), APAR_{(x, t)}表示像元x在t时间吸收的光合有效辐射(MJ m^-2月^-1), ε_{(x, t)}表示像元x在t时间的实际光能利用率(单位：gC/MJ)。

2.2.4 固碳释氧

根据光合作用方程式6CO₂+6H₂O→6O₂+C₆H₁₂O₆, 植被生产的有机物、固定的CO₂以及释放O₂的比例为1:1.47:1.07, 植被固碳释氧的量可以基于NPP的计算结果换算得到^[16]。

2.2.5 粮食生产

将县区统计部门获取的以乡镇为单元的粮食生产数据, 除以乡镇面积换算为单位面积粮食产量, 赋值到汾河上游行政区划矢量图并转换为栅格格式, 生成汾河上游流域粮食生产服务。

2.3 生态系统服务驱动机制分析

本研究筛选降雨量、相对湿度、日照时数、气温等气象因子以及坡长和坡向地形因子; 同时, 选取道路、居民点、人口等人文因素整合成综合性的人类活动强度指数。在ArcGIS软件支持下对降雨量、相对湿度、日照时数、气温等气象因子进行克里格插值; 基于DEM数字高程图, 采用ArcGIS软件中栅格表面数字分析模块中的坡度坡向分析工具进行提取; 所有自然因素栅格图通过ArcGIS中zonal statistic工具分配到各乡镇供统计分析。基于SPSS软件, 采用典型相关分析方法(Canonical correlation analysis, CCA)分析生态系统服务驱动因子。

2.4 人类活动定量化研究方法

整合人口密度、居民点、道路相关数据, 制订综合性的人类活动强度指数(HAI, human activity intensity), 采用层次分析法和专家打分法对3个因素进行权重分配居民点(0.3)、道路(0.3)、人口密度(0.4)。首先采用ArcGIS软件, 对纸质的县乡行政区划图、居民点、道路图进行矢量化。根据国家基础地理信息数据编码标准, 结合《山西省土地利用现状数据集》对居民点和道路影响力进行赋值(表 2, 表 3)。在ArcGIS软件支持下, 居民点采用克里格方法进行空间插值; 道路的插值采用线源距离衰减方法首先将道路按级别分级并计算其欧氏距离, 采用公式FEN=道路分级影响力×(1-实际距离/欧氏距离)计算各级道路影响力栅格图并进行加和, 生成研究区道路影响力栅格图。

2.5 数据来源

土地利用图由Landsat TM和Cbers-2B多光谱影像从网上下载并进行目视解译, 研究区DEM从1:50000地形图提取, 土壤类型从当地农业部门获取, 围绕汾河上游流域周围选取21个气象站点, 从国家气象科学数据共享服务平台下载降雨量、相对湿度、日照时数、气温、蒸散量、太阳辐射等气象因子。人口数据从研究区各县(市、区)统计部门获取, 乡镇行政区划图及居民点分布从各县(市、区)民政部门区获取。道路图由山西科学技术出版社出版的《山西省地图册》及山西省地图出版社出版的《山西省分县地图册》进行提取(图 2)。

3 结果与分析 3.1 生态系统服务时空格局

基于DEM栅格图将整个汾河上游流域分成47个子流域, 泥沙截持服务结果以子流域为单元输出。2000年研究区的北部和西部泥沙截持服务较高, 而中部区域较低。2008年研究区西部泥沙截持服务除北部区域泥沙截持服务较高外, 其他区域均较低(图 3)。2000—2008年间研究区泥沙截持服务有所降低, 从2000年的280.28 t/hm²降低到2008年的236.98 t/hm²; 除区域北部局部区域有所增强外, 占研究区87.5%的区域呈现减少态势; 泥沙截持服务减少幅度从东北到西南递增。

产水服务结果以子流域为基本单元输出。2000年流域产水服务大致呈现西北高东南低的格局, 产水量最高的区域集中于流域西部的岚县境内。2008年产水量服务与2000年空间格局有较大的差异, 大致呈现为北高南低的格局(图 3)。产水服务最高值集中于宁武和静乐区域内。总体上产水服务从2000年的146.33 t/hm²增加到2008年的151.07 t/hm²。栅格图层相减显示2000—2008年间产水服务增加的区域占整个汾河上游流域的89.6%, 尤其是流域东北部的宁武境内的子流域产水量增加幅度最高。

NPP生产服务结果以栅格为基本单元输出。2000年和2008年两年间汾河上游流域NPP生产服务均呈现出周边高、中间低的格局, 尤其以东南部区域NPP生产服务最高(图 3)。2000—2008年间, NPP生产服务平均值由4.07 t/hm²增加到4.93 t/hm²。NPP生产服务增高的区域占整个研究区总面积的93.7%;仅流域周边零星区域NPP呈减少态势。

固碳释氧与NPP时空格局一致, 总量从2000年的10.34 t/hm²增加到2008年的12.52 t/hm²。固碳释氧两年的变化量呈现出中间高边缘低的格局, 尤其是流域南部区域固碳释氧增加量最为明显, 流域边缘地方零星分布固碳释氧减少的斑块。

粮食生产服务以乡镇为单元输出。总体上研究区粮食生产服务有所降低, 从2000年的175.11kg/hm²减少到2008年的159.46 kg/hm²。2000—2008年间除宁武静乐大部、岚县西部粮食增加外, 占流域65%的区域粮食生产服务减少。从空间分布尺度来看, 2000年研究区粮食生产呈现西北高东南低的格局, 粮食生产较高的区域主要位于流域西部岚县、静乐西部、娄烦北部以及尖草坪区; 2008年, 粮食生产格局呈现北高南低的格局, 除岚县和宁武静乐西部区域外, 古交市及阳曲县的北小店及西凌井乡粮食生产减少明显(图 3)。

3.2 人类活动强度(HAI)

2000年人类活动强度(HAI)呈现北-南-东南逐渐增加的空间态势(图 4), HAI最低区域主要位于流域北部宁武县的舍庄(0.059)、涔山(0.061)、迭合寺(0.09)、东马坊(0.10)、东寨(0.11)、圪廖(0.11)、堂儿上(0.11)、怀道(0.11)等乡镇, HAI最高的区域主要位于流域东南部太原市区万柏林(0.51)、尖草坪(0.436), 以及古交市西曲(0.33)等乡镇。2008年HAI空间格局略有变化, 呈现北-西南-中部-东南逐渐增加的空间态势, HAI较低区域主要分布在流域北部宁武县余庄乡(0.071)、涔山乡(0.074), 迭合寺乡(0.11);岚县的王狮乡(0.11)、梁家庄乡(0.14);娄烦县的米峪镇(0.13)、盖家庄乡(0.14)、马家庄乡(0.14)、古交市的岔口乡(0.14);HAI最高值位于太原市区的万柏林区(0.50)和尖草坪区(0.44)。从时间上来看, HAI总体平均值由2000年的0.18增加到2008年的0.21;HAI随时间增加的区域主要分布在宁武县和静乐县搭界的堂儿上和娑婆乡等乡镇, 而流域西南部岚县、娄烦县部分区域HAI明显减少, 太原市区万柏林和尖草坪区HAI也略有降低。总体上, 城市化程度高的地方人类活动强度相对较高。

3.3 生态系统服务驱动因素分析

气象因素中只有气温在2000年和2008年呈现出相似的北低南高的空间格局, 其他因素如年降雨量、相对湿度、日照时数2000年和2008年间有较大的变化(图 5)。年降雨量和相对湿度由2000年的西高东低变为2008年的北高南低, 日照时数由2000年的东北高西南低变为2008年的南高北低。坡向和坡度2000年和2008年变化不大; 道路影响力、居民点影响力、人口密度在2000年和2008年均呈现相似的东南高西北低的格局, 但区域内变化幅度有所增大。

典型相关分析结果表明(图 6), 2000年第一个典型变量将泥沙截持服务从其他服务中区分出来(0.961);自变量组中对应的解释变量是坡度(0.832)。第二个典型变量将产水量分离出来(0.955), 其对应的解释变量是降雨量(0.849)。2008年分析结果与2000年类似, 第一个典型变量将泥沙截持服务分离出来(0.959), 其对应解释变量为坡度(0.960);第二个典型变量将产水服务分离出来(0.977), 对应解释变量为降雨量(0.837)与温度(0.798);与2000年相比, 2008年温度对产水量影响凸显出来。2000—2008两年差值的典型相关分析结果显示, 泥沙截持服务(0.964)与产水量服务(0.893)受日照时数(0.856)、相对湿度(0.840)和降雨(0.735)等自然因子的影响较大; 第二个典型变量表明NPP生产/固碳释氧(0.970)对应的解释变量为温度(0.685)和坡向(0.652)。

3.4 生态系统服务之间关系分析

生态系统服务之间存在着复杂的关系, 主要体现为权衡与协同关系。本研究基于SPSS统计分析软件, 采用皮尔逊相关系数显著性双尾检验对NPP生产、泥沙截持、产水量三类生态系统服务与人类活动强度2000—2008年间的变化量进行相关性分析(表 4)。从生态系统服务总量看2000—2008年间泥沙截持量降低, 而产水服务则呈现增加, 二者呈现一定的权衡关系, 但在空间耦合性上, 二者又呈现一定的正相关; 粮食生产和NPP与二者的关系不显著。人类活动强度与泥沙截持和产水量呈正相关关系。生态系统服务的权衡协同有复杂的地域性特征, 岚县、娄烦、市区(尖草坪区和晋源区)与流域全域呈现相似的关联特征, 即泥沙截持下降, NPP生产和固碳释氧增加; 古交市与阳曲县呈现出泥沙截持与粮食生产下降的同时NPP生产与固碳释氧呈现增加态势; 而宁武则呈现出四类生态系统服务均增加的协同格局(图 7)。

4 结论与讨论 4.1 结论

2000—2008年间, 泥沙截持服务减少了43.30 t/hm², 除北部宁武县部分区域增加外, 占流域87.5%的区域截持服务减小, 减小幅度从北到南依次加大; 同期产水服务增加了4.74 t/hm², 除西南局部地区有所减少外, 占流域89.6%的区域产水服务增加; NPP生产服务增加了0.86 t/hm², 仅流域边缘零星区域NPP生产服务减小, 占流域93.7%的区域NPP生产服务增加; 固碳释氧服务与NPP生产呈现相同的时空格局, 两年间增加了2.18 t/hm²。粮食生产服务总体上减少了15.65 kg/hm², 减少区域占到整个流域的65%。

相关性分析表明, 泥沙截持与产水服务在空间上呈现一定的协同关系, 这可能与二者在形成机理(生态水文过程)上的重叠有一定的关系; 泥沙截持服务和产水服务共同受到土壤-植被-大气系统, 以及径流、入渗、蒸发等水文过程等的影响。生态系统服务间消长关系具有复杂地域特征, 除流域北端宁武呈现出各服务协同关系外, 其他区域泥沙截持服务与NPP生产服务及固碳释氧呈现出权衡的关系。

人类活动强度空间格局与城市化程度呈一定的正相关, 由北向南依次升高; 2000—2008年间, 人类活动由0.18增至0.21;增加幅度最高的区域位于北部宁武、静乐交界处, 流域西南部岚县、娄烦以及万柏林区和尖草坪区人类活动强度有所降低。相关性分析表明, 人类活动强度与泥沙截持和产水量二者呈正相关, 表明由于退耕还林还草及城市生态建设等生态保护工程的实施, 人类活动趋于合理, 对泥沙截持和产水量有积极的影响。

典型相关分析表明地形在影响泥沙截持服务中起着决定性作用, 降雨和温度对产水服务起着重要作用; 年际间变化典型相关分析表明, 日照时数、相对湿度、降雨等共同驱动了泥沙截持服务和产水服务, 而温度和地形对NPP生产及固碳释氧服务有重要影响。

4.2 讨论

生态系统服务作为联系自然环境与人类福祉的的桥梁, 对其进行驱动机制分析是实现区域生态系统科学管理、协调人地关系, 实现社会经济可持续发展的重要内容^[17-18]。非生物环境(地理特征、地形地貌、气候、土壤)、生物特征(动植物微生物区系), 生态过程(光合、呼吸、蒸散发、有机质分解、生物地化循环等)是生态系统存在的基础, 也是生态系统服务形成的条件^[19]。生态系统服务驱动因素与其依存的生态过程密不可分, 本研究生态系统服务驱动分析结果与各类生态服务生态过程有很强的一致性, 如产水量作为降雨量扣除蒸散发后剩余的分量, 降雨量是其中重要的输入参数, 同时蒸散发本身又与降雨量和温度密切相关; 泥沙截持服务评估USLE模型中, 地形因子(坡长坡度)是重要的输入参数, 对泥沙的输出具有决定性影响; NPP生产服务评估所应用的CASA模型本身就是基于遥感数据、温度、降水、太阳辐射、植被、土壤数据共同的光能利用率模型, 坡向和温度对植被吸收光能有效辐射具有重要影响。另一方面也说明了典型相关分析方法在驱动机制中较强的适用性;

作为间接影响因素, 人类活动通过影响生境、生态系统结构、生物地球化学循环来影响生态系统过程进而影响生态系统服务。汾河上游流域作为传统的资源开采区, 矿区和城市区并存, 影响生态系统服务的因素繁多, 选择合适的人类活动表征指标困难, 简单的土地利用变化不能客观反应人类活动的内涵, 而资源开采的数据尤其是私挖乱采的数据无法获取, 对人类活动定量刻画的准确性造成了影响。

生态系统服务的尺度依赖性及影响因素空间单元的不一致也给生态系统服务驱动机制分析带来了很多困难。研究表明, 随着研究尺度的加大, 生态系统服务驱动因素有从自然因素为主向社会经济因素为主转换的特点^[20]; 这种尺度依赖性一方面可能与生态系统服务本身尺度特征有关, 另一方面, 也可能与尺度分析中的粒度有关^[21]。驱动因素中自然因素数据多基于自然地理单元, 而人文因素数据多基于行政单元, 空间单元上的不匹配势必影响结果的确定性; 自然和社会经济数据的有机整合也是生态系统服务驱动机制及人地关系等复合生态系统研究的难点。

致谢: 兰州大学张金茜、山西大学张娇、山西环保厅荣月静在统计数据收集与分析方面提供帮助, 特此致谢。