生态系统服务是指生态系统所提供的能够满足和维持人类生活需要的条件和过程, 从而为人类直接或间接提供的所有收益^[1-2]。它主要包括自然生态系统为人类提供的供给服务、调节服务、文化服务和支持服务^[3]。由于不同区域的生态系统结构与空间差异, 使得各项服务功能也具明显差异^[4]。而地形常常是影响山区生态结构和空间格局分异的重要因子^[5], 人类活动频繁、强烈的区域, 地形特征影响着景观格局和生态系统服务的空间分布^[6-7]。因此, 生态系统服务的空间异质性研究成为当前生态评估的热点课题。近年来, 大多学者主要集中在地形因子与土地利用类型分布、格局变化及生态系统服务价值的空间差异的研究中, 而缺少对生态系统服务功能在地形梯度上的定量表达及空间分异的定量化研究。如邱扬等^[8]分析了黄土丘陵小流域土地利用的时空变化及其与地形因子的关系；幸瑞燊等^[9]开展了基于地形梯度的三峡库区土地利用时空变化分析；李博等^[10]测算了石羊河流域各生态服务功能空间异质性；Yao等^[11]分析了浑河上游生态系统服务价值及其时空差异；杨锁华等^[12]分析了长江中游地区县域生态系统服务价值及其时空特征与地形梯度效应；徐煖银等^[13]利用网格法分析赣南地区土地利用结构及生态系统服务价值的时空分异。总的来说, 地形与生态系统服务存在着紧密联系, 影响着生态系统的分布、结构和功能^[14-15]。但当前涉及地形因素对生态系统服务空间异质性影响的研究报道较少。认知地形要素对生态系统服务空间差异的影响, 对深化认识生态系统服务分异, 合理调控人类活动, 因地制宜、开展有效生态系统管理具有重要意义。

甘肃白龙江流域是长江上游重要的水源地和生态屏障, 也是中国滑坡、泥石流灾害四大高发区之一^[16]。流域内沟壑纵横、地形起伏, 生态系统空间差异明显。然而, 白龙江流域虽然为区域人类活动提供了各种不同的生态系统服务, 但这些服务的空间差异及地形对其影响如何还有待进一步研究。基于此, 本研究以甘肃白龙江流域为例, 从高程、坡度、地形起伏度、地形位指数四方面分析其典型生态系统服务地形梯度效应, 旨在明确流域生态系统服务垂直截面的空间异质性, 为流域生态系统可持续管理和生态文明建设提供科学依据。

1 研究区概况

甘肃白龙江流域(32°36′—34°24′N, 103°00′—106°30′E)位于甘肃省东南部(图 1)。白龙江是长江支流嘉陵江的最大支流, 发源于甘南碌曲县与川北若尔盖县交界的郎木寺, 主要流经甘肃省的迭部县、舟曲县、宕昌县、武都区和文县, 入川汇入嘉陵江。流域面积约为1.84万km²。流域地势西北高、东南低, 地形起伏差异明显, 地貌多样, 以山地、高原、河谷、丘陵及黄土为主。气候属亚热带向北温带过渡区, 上中下游依次为温带湿润-暖温带湿润-北亚热带湿润气候, 年降水量约500—900mm, 年均温约2℃—15℃。研究区森林广布, 具有明显的水平和垂直分异规律。近年来, 区内社会经济的发展、城镇化的扩张, 加上特殊的地理条件及降水特征, 水土流失严重, 滑坡、泥石流等自然地质灾害频发, 对流域生态环境产生一定压力^[17], 另外, 退耕还林、天然林保护、野生动物及自然保护区建设等工程使得该流域生态环境质量得到一定的改善。

2 数据处理与方法 2.1 数据来源

所用数据主要有：来源于中国科学数据库, 空间分辨率为30m的DEM数据；参考《土地利用现状分类标准》(GB/T21010―2007), 并结合研究区实际情况, 对2016年Landsat ETM+遥感影像进行目视解译获得的土地利用数据；2016年流域及周边18个气象站点的降雨、气温等气象数据；甘肃省第二次土壤普查资料、《甘肃土壤》和《陇南土壤志》等土壤理化性质数据；2016年流域各县区的乡镇尺度人口、农作物生产、收入与消费状况等社会经济统计数据等。

2.2 InVEST模型

InVEST模型是由斯坦福大学、大自然保护协会(TNC)和世界自然基金会(WWF)实施的“自然资源计划”(Natural Capital Project)中开发, 量化多种生态系统服务功能的生态系统服务功能综合评价模型^[18-19]。主要包括碳储存、生物多样性、土壤保持、水体净化、木材收获管理等模块^[18]。根据研究区流域生态环境特点及数据的可获取性, 选取碳储量、水源涵养、生境质量和土壤保持4种生态系统服务进行评估。

2.2.1 碳储量

该模型主要利用各土地利用与覆被类型空间分布格局及基本碳库来评估区域生态系统碳储量, 生成碳储量空间分布图。其计算公式^{[18, 20-22]}为：

(1)

式中, C_tot表示总碳储量(t), C_above表示地上生物碳储量(t), C_below表示地下生物碳储量(t), C_soil表示土壤碳储量(t), C_dead表示枯落物碳储量(t)。碳密度数据来源于收集和整理李克让^[20]、解宪丽^[21]、方精云^[22]等成果和参考模型数据库^[18]。

2.2.2 水源涵养

该模型主要利用降雨量、蒸散量、土壤深度、植物可利用水等数据和水量平衡的方法计算流域产水量, 在此基础上结合土壤渗透性、地表径流差异和地形等因素计算得到水源涵养。其计算公式^{[18, 23]}为：

(2)

式中, Y_{x, j}、AET_{x, j}分别是土地类型j中栅格单元x的年产水量(mm)和年实际蒸散量(mm), 其中AET_{x, j}为Budyko曲线的近似值, P_x是栅格单元x的年降水量(mm)。

(3)

式中, Watershed pixel count为集水区汇流累积栅格数量(无量纲), Soil depth为土壤深度(mm), Percent slope为百分比坡度。

(4)

式中, WR为水源涵养量(mm), Velocity为径流系数, TI为地形指数(无量纲), K_soil为土壤饱和导水率(cm/d), Y为产水量(mm)。

2.2.3 生境质量

该模型主要利用土地利用、生态胁迫因子及敏感性属性表和社会经济数据等数据来评估区域的生物栖息地质量。参数设置来源于InVEST模型使用手册^[18]、前人相关研究^[24]以及专家打分获得, 其计算公式^[24]为：

(5)

式中, Q_{x, j}是土地类型j中栅格单元x的栖息地质量(无量纲), H_j土地类型j的栖息地属性/生境适宜度, D_{x, j}是土地类型j中栅格单元x的栖息地退化程度, k是半饱和系数, 其值通常占栅格单元分辨率大小的一半, z为模型默认的归一化常数(z=2.5)。

2.2.4 土壤保持

该模型主要利用DEM、土地利用、气象数据、土壤属性等数据, 在基于通用土壤流失方程(USLE)上, 同时考虑植被在缓减侵蚀和截留上坡来沙的作用计算得到土壤保持量。其计算公式^[25]为：

(6)

(7)

(8)

(9)

式中, PKLS_x是栅格单元x的土壤潜在侵蚀量(t), R_x是降水侵蚀力(MJ·mm/(hm²·h·a))、K_x是土壤可蚀性(t·hm²·h/(hm²·MJ·mm))、LS_x是坡长坡度因子, USLE_x是栅格单元x的土壤实际侵蚀量(t), SEDR_x是栅格单元x的土壤实际持留量(t), USLE_y为上坡栅格y的实际侵蚀量(t), SE_x为栅格单元x的持留率；SE_z为上坡栅格的泥沙持留量(t), SEDRET_x是最终求得的土壤保持量(t/hm²)。

2.3 地形数据处理 2.3.1 地形起伏度

地形起伏度是指某一分析窗口内最高与最低点间的高差, 是定量描述地形形态、划分地貌类型的重要指标。该指标主要利用ArcGIS中栅格邻域计算工具获取, 其计算公式^[26]为：

(10)

式中, C_max为分析窗口内的最大高程(m)、C_min为分析窗口内的最小高程(m)。

2.3.2 地形位指数

地形位指数是复合分析空间任意一点高程和坡度属性信息的指标, 可综合反映出地形条件的空间分异。其计算公式^[27]为：

(11)

式中, T为地形位指数, E及E₀分别为空间内任一栅格的高程(m)和平均高程(m), S及S₀分别为空间内任一栅格的坡度值(°)和平均坡度值(°)。一般地, 高程低、坡度小的栅格地形位指数小, 反之越大。

2.3.3 地形数据分级标准

为分析不同地形梯度下生态服务空间分布, 根据甘肃白龙江流域地形分布特点, 采用分位数法^[28]将高程、坡度、地形起伏度和地形位指数分别分为5级。按照其数值大小分别命名为Ⅰ、Ⅱ、Ⅲ、Ⅳ和Ⅴ级(表 1)。

3 结果与分析 3.1 甘肃白龙江流域地形特征

利用ArcGIS软件, 基于DEM数据得到甘肃白龙江流域地形特征图(图 2)。由图 2可见, 白龙江流域地形差异较大, 水系所处位置各地形因子均处于较低值, 呈现出较明显的河谷、黄土地貌及高原山地等地貌特征。其中白龙江流域高程由西北向东南递减：海拔最高达4860m, 最低为568m。流域东北部坡度较缓, 西部最陡处坡度达88.01°。地形起伏度由北、东两个方向逐渐向西、南升高。地形位指数最大值为2.45, 说明局部地区坡度大、高程高等地形特征。

3.2 甘肃白龙江流域生态系统服务分布特征

基于土地利用、气象、土壤、DEM等基础数据, 利用InVEST模型估算出白龙江流域2016年碳储量、水源涵养、生境质量、土壤保持服务功能的空间分布图(图 3)。甘肃省白龙江流域碳储存服务在0—116.43t之间。碳储量高值区(72.46—116.43t)主要分布在林地区域, 如宕昌南部林区、白水江南岸高山区、拦坝河中上游等；中值区(54.94—72.46 t)主要集中在迭部北部、宕昌东部和武都区；低值区(0—54.94t)主要分布在农业种植区、城镇、迭部北部高海拔高寒草甸区, 重大滑坡地带及高山裸岩地区。水源涵养服务在0—321.10mm之间, 其空间分布格局基本呈东北低西南高, 高值区(71.78—321.10mm)集中于武都东南部和文县南部；中值区(46.59—71.78mm)分布在迭部、宕昌西部和舟曲中西部；低值区(0—46.59 mm)连片分布于宕昌-舟曲-武都段白龙江两岸及其以北的人口密集区和文县北部农作区。生境质量分布范围为0—0.92之间, 高值区(0.61—0.92)主要是自然保护区和林业管护区, 如白水江国家自然保护区、迭部的达拉、阿夏、多尔自然保护区及林业分布区。中值区(0.32—0.61)集中于林区与人类活动干扰区的过渡地带；低值区(0—0.32)则主要分布在人类活动相对频繁的河谷地带、城乡聚落区、耕作区、裸地和低覆盖草地等区域。土壤保持服务范围在0—2326.19t/hm²之间, 其高值区(706.30—2326.19t/hm²)主要分布在武都东南部和文县南部, 这些区域多属于地表扰动较少的石质性山区或自然保护区；中值区(479.05—706.30t/hm²)主要分布在迭部南部、宕昌局部和文县西北部等人类干扰强度较弱的过渡地区；低值区(< 281.30t/hm²)主要分布在人类活动相对频繁、工农业相对发达的舟曲—武都—文县段白龙江河谷地带。

3.3 生态系统服务的地形梯度分布特征

利用ArcGIS软件, 将碳储量、水源涵养、生境质量、土壤保持4种典型生态系统服务分别与不同地形梯度叠加分析, 获取每个地形梯度级上的分布趋势图(图 4)。其中, 碳储量、生境质量和水源涵养在3000m以下区段, 随高程增加而不断增加, 而3000m以上的海拔区域, 这三项服务功能则逐渐减少；随着坡度、地形起伏度的增加呈递增趋势；随着地形位指数的增加, 碳储量一直呈递增趋势, 其他服务则均呈先增加后降低的特点；而土壤保持随着高程的增加呈先降低后增加的趋势, 随着坡度、地形起伏度的增加则呈先增加后降低的趋势, 而随着地形位指数的增加则呈起伏波动的特点。

在高程上, 由于流域内高程呈现西北高东南低, 且地处秦岭地质构造褶皱带, 沟深坡陡, 各地类分布较为分散。森林主要分布在海拔范围1716—3049m之间(高程Ⅱ-Ⅳ级), 如武都南部、文县西南部、舟曲大部及迭部局部沟壑处。而高程第Ⅴ级主要以高寒草甸、沙地裸地、高寒荒漠及高山积雪等地类为主, 因此, 各生态系统服务随高程增加而增加, 到第Ⅴ级则下降的特点。另外, 由于高程中低值区也是人类活动相对频繁、滑坡泥石流高发地区, 从而使得土壤保持服务在高程上呈现出先降低后增加的特点。

由于不同坡度带上人类活动强度的差异, 相关生态系统服务功能也会随之变化。如坡度较小处主要集中在白龙江流域东北部和中南部, 以耕地和建设用地为主, 受人类活动干扰大。随着坡度的增加, 水分、热量等自然条件适中, 人为活动影响逐渐减少, 山地林区、草地等植被长势渐好, 因此, 碳储量、生境质量等也随之呈上升趋势。而坡度升高至第Ⅴ级, 植被覆盖较少, 不利于水土的保持, 使得土壤保持呈现略微下降的趋势。

由于白龙江流域内地形起伏度和坡度的空间格局较为相似, 因此地形起伏度上各生态系统服务变化特点与在坡度梯度上具有类似性。相较于其他各地形因子, 土壤保持服务在地形起伏度上变化最大, 说明其对水土流失具有突出影响。

流域中, 地形位指数高值区主要为迭部南部、舟曲大部、文县西北部等地区；低值区在文县大部、宕昌东北部及武都大部等地区, 主要以农业种植区及城镇为主；中值区人类活动干扰较少, 以林区为主；从而各生态系统服务随着地形位指数增加而呈较为波动趋势。

3.4 多重服务重要区的地形梯度效应

为了解各项生态系统服务在垂直空间上供给能力大小, 本研究将各项服务功能值分别大于各自平均值的地区视为该项服务的重要区域。通过空间分析功能, 将碳储存、水源涵养、生境质量和土壤保持4项生态系统服务的栅格图层进行叠加, 统计每个栅格单元上提供各项生态系统服务能力, 即多重服务供给的数目多少以及是否为该类服务的重要区域, 并根据各项服务重要区域的频次划分为五级：极度重要、高度重要、中度重要、一般重要、一般(表 2), 得到白龙江流域多重生态系统服务分布图(图 5)。

如图 4所示, 甘肃白龙江流域多重生态系统等级为极度重要的区域主要集中白水江自然保护区、迭部县中部地区、舟曲县武平—插岗梁—博峪等山林区、宕昌县迭山林区、武都东南部地区等。高度重要的区域主要分布在自然保护区和林区。中度重要分布相对比较均匀、分散, 在流域农耕区、林区、农牧交错区均有分布。等级为一般重要主要分布在宕昌县、文县、舟曲—武都段的白龙江两岸。等级为一般的区域主要集中在文县至碧口镇段白水江两岸、宕昌东北部和武都北部。

利用编码分级法和空间叠置法分别提取各地形梯度的白龙江流域多种生态系统服务面积占比, 随各梯度的升高其综合生态系统服务的面积所占的比例也随之变化(图 6)。各地形因子第Ⅰ、Ⅱ梯度中占主导地位的类型均为一般等级, 第Ⅲ、Ⅳ和Ⅴ梯度中则以高度重要为主, 但各类型面积所占比例随地形梯度变化较大, 面积比例的变化大致表现为：①在高程上, 等级为一般的区域在第Ⅰ、Ⅱ和Ⅴ梯度中占据优势, 随梯度增加先降低后略增；而多种生态系统服务重要区在第Ⅲ、Ⅳ梯度中面积占比较大, 一般重要的区域面积占比呈先持续降低到第Ⅴ梯度略增趋势, 其他等级重要区面积比重则呈现先增加后减少的趋势。②坡度与地形起伏度层面, 各类型区在5个梯度上的占比变化趋势相较一致。等级为一般的区域在第Ⅰ梯度占比最大, 随梯度增加逐渐降低；中度重要和高度重要的面积占比随梯度增加呈先增后减趋势；极度重要区域则呈持续增加趋势。③地形位指数等级为一般与一般重要的面积占比呈持续减少的趋势；中度重要与高度重要的面积比重则呈先增后减的态势；极度重要面积比重较小, 随梯度呈轻微波动。综合来看, 白龙江流域内生态系统服务重要区逐渐向两极分化, 多种生态系统供给服务的重要区域在中高梯度上占据主导, 随着地形梯度的上升, 人类活动变小, 其供给的服务类型越来越多、也越来越明显；而等级为一般的区内各项生态系统服务都比较弱或逐渐越偏向于某一单一类型的生态系统服务, 其比重在低梯度占据优势。

4 讨论与结论 4.1 讨论

地形是影响生态系统服务分布、结构及功能的重要基础, 本文选取高程、坡度、地形起伏度和地形位指数, 探讨了甘肃白龙江流域重要生态系统服务分别在这4种地形要素上的变化, 旨在对生态系统服务进行垂直空间上的截面剖析, 为流域社会生态协调发展和生态系统服务优化提供科学参考。这4种地形指数也被广泛应用于山地生态系统服务在地形梯度上的分布特征的研究中^{[6, 12-13, 29-31]}。根据白龙江流域“山大沟深”、地形差异突出等特点, 相较于其他划分方法, 为避免分级时某梯度内面积过大或过小, 而导致在后续操作过程中易造成较大损失, 最终结果偏差较大, 本研究中选用分位数法作为等级划分方法, 相关领域的一些学者也采用了这种分级方法^[32-33]。

本文通过选取碳储量、水源涵养、生境质量、土壤保持4种典型服务功能进行评估, 探讨甘肃白龙江流域综合服务的不同地形梯度效应。相较于以往只基于单一服务的地形梯度效应研究, 本文则采用多重服务综合效应较全面反映其地形梯度特征, 结果表明：各生态系统服务功能值随梯度的增加大致呈先增后减的趋势, 可能原因是由于中低梯度处适宜人类生产生活, 随梯度的升高, 人为扰动影响减小, 植被覆盖更高, 各项服务也随之增加。而至第Ⅴ级, 主要以裸岩、荒漠及未利用地等类型为主, 土壤保持能力较低, 各项服务呈现下降的趋势。这与谢余初^[17]和郜红娟^[29]的研究结果相似。

明确生态系统格局和功能服务最终目的是为流域生态系统服务管理和优化提供科学依据和理论支持^[34]。基于以上研究结果提出相关生态保护建议与措施：(1)在多重服务低值(Ⅰ、Ⅱ梯度)区段, 即宕昌-舟曲-武都段白龙江两岸及其以北的人口密集区和文县北部农作区。应严格管控人类行为, 加强土地的集约化利用, 缓解人地矛盾, 同时注重污水处理、城区绿化和生态宜居等环保工程的建设；(2)在多重服务中高值(Ⅲ、Ⅳ梯度)区段, 即迭部大部、宕昌南部、武都东南部。应建立和扩大生态系统服务保护区, 加强森林的管护工作, 重视林草恢复, 强化生物多样性保护, 提升环境承载力。(3)在多重服务低值(Ⅴ梯度)区段, 即各县区裸地、高寒荒漠等低植被覆盖区域, 应加强高寒草甸保护、合理放牧, 减少人类活动干扰, 降低生态环境脆弱性, 提高抗风险能力。

4.2 结论

(1) 甘肃白龙江流域碳储量、水源涵养、生境质量、土壤保持服务功能的空间分布格局上, 碳储量范围在0—116.43t之间；水源涵养量范围在0—321.10mm之间；生境质量在0—0.92之间；土壤保持在0—2326.19t/hm²之间。碳储量、生境质量及土壤保持高值区主要集中在自然保护区、林区等人类干扰较少的地区, 低值区则集中于人类活动较频繁、自然灾害较多的河谷地区。水源涵养高值区主要分布在林地分布区, 低值区则在河谷、城镇、农耕地区。

(2) 各生态系统服务功能在梯度上呈现较为明显的地理层级性分布的特点。碳储量、水源涵养和生境质量随着坡度和地形起伏度的增加而增加, 土壤保持则呈先增加后降低的趋势。随地形位指数增加, 碳储量一直呈递增趋势, 其他服务则均呈先增加后降低的特点。在高程上各服务呈现波动特点, 在3000m以下的区段内, 碳储量、水源涵养和生境质量随高程增加而不断增加, 而3000m以上的海拔区域, 则逐渐减少。

(3) 流域内综合生态系统服务重要区逐渐向两极分化, 多种生态系统供给服务的重要区域在中高梯度上占据主导, 其供给的服务类型越来越多、也越来越明显；而供给服务能力一般的区内各项生态系统服务都比较弱或逐渐越偏向于某一单一类型的生态系统服务, 其比重在低梯度占据优势。