生态系统服务是指生态系统与生态过程所形成及所维持的人类赖以生存的自然环境条件与效用, 是人类直接或间接从生态系统获得的所有收益^[1-2]。各项生态系统服务的供给在动态变化过程中存在着复杂的相互关系, 生态系统服务之间表现为此消彼长的状态, 称为权衡；两种或多种生态系统服务表现为同增同减的形式, 称之为协同^[3]。生态系统管理不能只追求单一的生态系统服务效益, 而要兼顾多种生态系统服务权衡的过程, 使其综合效益最大化, 促进整个区域的平衡发展^[4]。喀斯特地区脆弱的生态环境、破碎的地表形态以及水土流失导致的石漠化现象, 其实质是生态系统结构遭到破坏, 从而导致生态系统功能的下降与丧失, 更深层次地反映了协同关系的损害^[5-6]。厘清喀斯特地区生态系统服务之间权衡和协同的空间关系, 对喀斯特地区进行石漠化综合治理以及提升区域生态系统服务具有重要意义^[7-8]。

近年来, 喀斯特地区生态系统服务的研究较多关注于土地利用或气候变化背景下的单一生态系统服务^[8-9], 如Feng等^[10]运用RUSLE(the Revised Universal Soil Equation)模型和¹³⁷Cs方法进行广西峰丛洼地土壤侵蚀的对比模拟；侯文娟等^[11]基于SWAT(Soil and Water Assessment Tool)模型分析喀斯特山区的产流服务以及不同服务变量的空间变异。即使考虑多种生态系统服务, 也偏向于研究喀斯特地区生态系统服务价值的计算或各项生态系统服务的时空变化规律, 如张明阳等^[12]借助价值当量方法进行不同尺度的生态系统服务价值评估；尚二萍和许尔琪^[13]对黔桂喀斯特山地主要生态系统服务进行时空变化分析。在研究方法方面, 喀斯特生态系统服务权衡/协同关系的研究以定性分析或简单线性关系为主, 定量化的研究相对较薄弱^[9]。如何针对喀斯特地区强烈的空间异质性量化生态系统服务之间的空间权衡关系是生态系统服务研究中的关键问题。均方根偏差法是目前判断及量化任意两个或两个以上生态系统服务之间权衡关系的简单有效的方法^[14]。地理探测器通过分析地理现象的空间分层异质性来揭示自变量对因变量空间分异的解释程度^[15]。因此本文融合均方根偏差法权衡关系的方向性判断以及地理探测器空间分析的优势, 以空间权衡度的概念定量化表征喀斯特地区生态系统服务之间的空间权衡关系。

喀斯特石漠化问题在于水土流失引起的地表土壤损失、基岩裸露、土地生产力丧失和生态环境退化, 因而土壤保持与水源涵养服务成为喀斯特地区石漠化遏制与生态恢复的核心^[16], 同时土壤保持支撑着土壤中水分和营养物质的持续供给^[17], 是植被生长与固碳的必然需求。鉴于此, 本文聚焦于贵州典型峰丛洼地流域的水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系, 在海拔、坡度、降水和植被覆盖度等环境因子梯度下和地貌形态类型区内开展生态系统服务空间权衡度计算及其分异特征的研究, 以期推动喀斯特地区生态系统服务的优化以及促进区域可持续发展。

1 研究区概况

三岔河流域位于贵州省西北部, 发源于贵州省西部乌蒙山, 是乌江南源一级支流(图 1)。河流全长325.6km, 地理位置介于26°06′N—27°00′N, 104°54′E—106°24′E之间, 流域面积4860km²。三岔河流域地处亚热带季风气候区, 全年气候温和湿润, 降水主要集中在5—10月, 多年平均降水量为963—1460mm。研究区属于典型的喀斯特峰丛洼地区, 水源涵养能力弱、土层浅薄且破碎化、植被系统稳定性差致使区域生态环境敏感脆弱, 石漠化现象突出, 人地关系矛盾加剧^[18-19]。

2 数据来源与方法 2.1 数据来源 2.1.1 水文和气象数据

根据《中华人民共和国水文年鉴—长江流域水文资料乌江区》记载, 三岔河流域的上游和中游地区各设有一个水文站, 分别是阳长和龙场桥水文站。由于InVEST模型模拟的产流量需与流域出口处的实测径流数据进行验证, 因此本文选取鸭池河水文站(位于六冲河和三岔河汇流的乌江干流处)和洪家渡水文站(位于六冲河流域出口处), 将其逐月平均流量数据的差值作为三岔河流域的月径流量。气象数据来源于国家气候中心, 本文选取三岔河流域及周边区域29个气象站点的日值数据, 运用专业气象插值软件ANUSPLIN^[20]进行气温、降水、蒸散发等气象数据的插值。

2.1.2 地形和土壤数据

DEM(Digital Elevation Model)高程数据, 采用Google Earth影像下载, 空间分辨率为9m。地貌形态类型数据来源于中国科学院资源环境数据中心(http://www.resdc.cn)(图 2)。土壤机械组成数据通过基于世界土壤数据库(HWSD)的中国土壤数据集(v1.1)获取, 分辨率为1km, 来源网址为：http://westdc.westgis.ac.cn。土壤的根系深度数据通过国家地球系统科学数据共享服务平台(http://soil.geodata.cn/data/dataresource.html)获取。

2.1.3 植被覆盖和土地利用数据

NDVI(Normalized Differential Vegetation Index)(图 2)数据以HJ1A/B CCD(30m)、GF1 WFV(16m)(http://www.cresda.com/CN/)以及MODIS MOD09GQ(250m)(https://lpdaac.usgs.gov)产品数据为数据源, 经过辐射定标、大气校正、正射校正等预处理后, 通过近红外和红光波段反射率线性组合的方式获取。预处理数据经过异常值处理、数据镶嵌、投影变换和目标区域剪裁等工作后得到最终的NDVI数据, 空间分辨率为30m。西南地区降水充沛、空气对流旺盛导致当地云量较多, 对于云污染很严重的月份, 本文采用STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)时空适应性反射率融合模型^[21], 结合高时间分辨率的MODIS数据在HJ、GF-1的空间分辨率上进行NDVI的预测, 模拟缺失月份内与MODIS数据相应时间的高分辨率的NDVI数据。植被类型数据来源于中国科学院资源环境数据中心(http://www.resdc.cn), 空间分辨率为1km。土地利用数据以2015年两景Landsat TM遥感影像(美国地质调查局官网http://glovis.usgs.gov/)以及野外考察定点数据为数据源, 运用ENVI软件进行监督分类解译得到三岔河流域的土地利用类型。

2.2 方法 2.2.1 RUSLE模型

RUSLE模型由于其结构精简、参数合理、模拟效果较好, 是广泛应用的土壤侵蚀模型之一^[22]。但是, 由于喀斯特地区特殊的石漠化情况以及复杂的地形条件, 直接使用RUSLE模型会使模拟结果产生偏差^[13]。Dai等^[23]通过人工降雨模拟试验发现表层土壤侵蚀与基岩裸露率呈显著负相关关系, 相关系数(R)为-0.076。Gao和Wang^[24]在三岔河流域开展的土壤侵蚀研究表明, 确定性系数(R²)与基岩裸露率的乘积可以解释基岩裸露率对土壤侵蚀减少的贡献。因此本文根据研究区石漠化情况的不同, 采用不同程度基岩裸露率的校正系数α(表 1)对RUSLE模型进行修正。坡长坡度因子解释了地形对土壤侵蚀的影响, 喀斯特地区复杂的地上地下二元结构以及地表汇流的不连续性会增加坡长因子对地形数据的敏感性。Feng等^[10]研究发现由于岩溶地区地表径流的不连续性, RUSLE模型中的坡长因子与累积面积阈值和DEM的空间分辨率密切相关, 累积面积阈值越高, 均方根误差也越高。因此本文选择9m分辨率的高精度DEM数据作为基础数据计算坡长因子以提高模型模拟的精度。同时, 已有研究发现三岔河流域的土壤侵蚀量与气候因子之间存在以五年为周期的波动性变化^[24], 因此, 为消除气候年际差异对流域尺度空间权衡分析结果的影响, 本文基于2013年到2017年三岔河流域的土壤侵蚀、水源涵养和植被固碳的优化模拟, 以五年的滑动平均数值进行流域生态系统服务之间空间权衡关系的研究。RUSLE模型的计算公式如下：

(1)

式中, A为实际土壤侵蚀量(t hm^-2 a^-1)；α表示不同石漠化程度下基岩裸露率的校正系数；R为降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1), 反映由降水引起土壤侵蚀的潜在能力^[25]；K为土壤可蚀性因子(t hm² h hm^-2 MJ^-1 mm^-1)；LS为坡长坡度因子, 也称为地形因子；C是植被覆盖因子, 使用30米分辨率的高精度NDVI反演数据计算；P是水土保持措施因子, 与土地利用类型有关。本文参考前人在喀斯特地区的研究^[26-27], 根据研究区的土地利用类型对P因子赋值。

2.2.2 InVEST模型

InVEST模型中的Water Yield模块综合考虑研究区的气候、地形、土壤深度和土地利用类型等因素, 基于水量平衡原理以每栅格单元的降水量减去实际蒸散发后的水量计算栅格水源涵养量^[28]。公式如下：

(2)

式中, Y(x)表示研究区每个栅格单元x的年产流量(mm), AET(x)是每个栅格单元x的年实际蒸散发量(mm), P(x)是每个栅格单元x的年降水量。

在水量平衡公式(2)中利用Budyko假设^[29-30]计算水量平衡的蒸散部分AET(x)/P(x):

(3)

式中, 为Budyko干燥指数, 被定义为潜在蒸散PET(x)与降水P(x)的比; PET(x)表示每个栅格单元x的年潜在蒸散量(mm), 由标准Penman-Monteith(P-M)^[31]公式求得；ω为植物年需水量与降水量的比值, 它是描述自然气候与土壤性质的参数。

2.2.3 CASA模型

CASA模型是目前进行NPP研究最普遍和有效的方法之一, 适用于区域NPP的计算和全球尺度碳循环的研究^[32]。NDVI数据是CASA模型中多个因子计算的基础数据, 基础数据的精度及遥感影像反演的可靠性对NPP的计算结果至关重要。因此本文选择高时间分辨率和空间分辨率的遥感影像反演的高精度NDVI数据来计算NPP, 公式如下：

(4)

式中, NPP_t、APAR_t和ε_t分别表示月份t的植被净初级生产力(gC m^-2 month^-1)、植被吸收的光合有效辐射(gC m^-2 month^-1)和实际光能利用率(gC/MJ)。

植被吸收的光合有效辐射与太阳总辐射和植物自身的特征有关, 可用公式(5)计算。

(5)

式中, SOL_t表示t月的太阳总辐射量(MJ m^-2 month^-1)；FPAR_t表示植被对入射光合有效辐射的吸收分量；常数0.5为植被利用的太阳有效辐射占太阳总辐射的比例。

Potter等^[33]认为实际光能利用率主要受温度和水分的影响, 当温度和水分达到理想条件时, 植被具有最大光能利用率, 其公式为：

(6)

式中, ε_t表示实际光能利用率(gC/MJ)；T_{max, t}和T_{min, t}分别表示高温和低温的胁迫系数；W_t为水分胁迫因子；ε_max是理想条件下的最大光能利用率(gC/MJ), 不同植被类型的ε_max值不同。本文参考朱文泉等^[34]对中国典型植被类型ε_max的模拟结果, 同时结合董丹和倪健^[35]对西南地区植被类型ε_max的改进, 确定了三岔河流域各种植被类型的ε_max值。

2.2.4 空间权衡度

本文从数值变化与空间变异综合的角度, 构建了融合均方根偏差法与地理探测器的生态系统服务空间权衡度指标, 从传统意义上权衡关系的判断到融合了权衡方向、数值量化和空间分异性的探测。均方根偏差法定量化表征单个生态系统服务的标准化数值与平均生态系统服务标准化数值之间的差异, 通过某一对生态系统服务坐标点到1 : 1线的距离判断生态系统服务之间是否存在权衡关系^[14]。地理探测器是探测和利用要素的空间分层异质性并揭示其背后驱动力的统计学工具^[15]。地理探测器认为如果某个自变量对某个因变量有重要影响, 那么自变量和因变量的空间分布应该具有相似性^[36-37]。基于这样的假设, 该方法采用q值度量自变量对因变量空间分异的解释程度。具体地, 本文运用均方根偏差法判断水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间是否存在权衡关系, 进而分别以水源涵养量和植被固碳量为自变量(X), 土壤侵蚀量为因变量(Y), 以地理探测器中体现数值空间变异的解释力因子表征喀斯特生态系统服务空间权衡度。

3 结果与分析 3.1 喀斯特生态系统服务模拟验证及其空间特征

贵州省水利厅发布的2011—2015年水土保持公告中2015年全省喀斯特区域土壤侵蚀模数为279.47t km^-2 a^-1, 2009年贵州省毕节市鸭池高原山地区平均侵蚀量为4.72 t km^-2 a^-1[38]。本文中修正的RUSLE模型模拟的三岔河流域土壤侵蚀模数为3.19t hm^-2 a^-1, 与上述结果基本相近。根据水利部土壤侵蚀分类分级标准^[39], 侵蚀等级属于微度侵蚀(图 3)。通过与三岔河流域基岩裸露率情况的对比, 印证了Feng等^[10]提出的基岩露头和薄土层地区无土可流的情况可能会导致岩溶地区的土壤侵蚀量低于非岩溶地区。基于2015年的实测月径流资料, 本文采用Nash-Sutcliffe效率NES和确定性系数R²来验证InVEST模型模拟月、年产流量的可靠性。结果表明, InVEST模型模拟的水源涵养量与实测值之间具有很高的相关性, 即NSE=0.802, R²=0.938, 模拟效果良好, 结果较为可信。图 3显示三岔河流域五年平均的水源涵养量空间分布差异明显, 呈现出由北向南逐渐递减的分布规律, 水源涵养量范围为118.77—1230.14mm, 均值为917.58mm。张明阳等^[40]在桂西北西部喀斯特地区运用CASA模型模拟的NPP结果为422.73gC/m²；王冰等^[41]在贵州省运用光合作用与呼吸作用相分离的模型计算得出的喀斯特地区NPP值为407 gC/m²。本文得出的三岔河流域年NPP积累量范围为0—1035.68 gC/m², 平均值为459.13gC/m², 该值与上述喀斯特地区NPP的研究结果较为一致, 证明本文CASA模型模拟的NPP值较为可信。三岔河流域NPP的累积量呈现西北高、东南低的特征(图 3), 这一分布特征与研究区的植被类型密切相关。研究区的西北地区分布有大面积的草原, 较之于东南地区的灌丛, 最大光能利用率较高, 植被所能利用的水分条件较好。

3.2 环境因子梯度下的喀斯特生态系统服务空间权衡关系

生态系统服务之间的复杂关系受到自然因素(如海拔、坡度、气候等)和人为因素(包括政策、市场、偏好等)的作用。其中, 自然因素是生态系统服务时空分布的基础^[42]。王欢等^[43]基于地理探测器的喀斯特地区土壤侵蚀的定量归因研究发现土地利用和坡度是决定土壤侵蚀空间异质性的主导因子；Wang等^[44]认为大面积的陡坡垦荒是造成水土流失和岩溶石漠化的主要原因；熊康宁等^[38]在典型喀斯特石漠化治理区发现随着植被的生长和恢复, 水土流失量逐年降低, 石漠化治理区的植被覆盖度与保土作用存在明显的正相关关系。地质、地貌因素是脆弱生态系统得以存在、发展的载体与物质基础, 地貌类型从宏观上控制了自然生态环境的特征与区域水土流失强度, 直接决定生态系统服务的供给与维持^[45-46]。

3.2.1 环境因子梯度下水源涵养-土壤侵蚀空间权衡关系

水源涵养与土壤侵蚀之间的空间权衡度在不同环境因子梯度下表现出明显的空间分异性(图 4)。植被覆盖度对水源涵养-土壤侵蚀之间空间权衡度的影响最为显著。在植被覆盖度因子的梯度下, 水源涵养与土壤侵蚀之间空间权衡度的最大值出现在植被覆盖度0—0.2的区间内, 为0.55。随着植被覆盖度的增大, 二者之间的空间权衡度逐渐减小。不同梯度海拔因子的运行结果表明, 水源涵养与土壤侵蚀之间的空间权衡度在以1000m为界限的低海拔和中海拔地区表现出明显的空间差异性。具体表现为二者之间的权衡度在800—1000m的低海拔区域未通过显著性检验, 而在海拔大于1000m的中海拔区域权衡度较高。坡度和降水因子对水源涵养-土壤侵蚀之间空间权衡关系的影响较小, 但权衡度在不同梯度下仍表现出明显的差异性。

3.2.2 环境因子梯度下植被固碳-土壤侵蚀空间权衡关系

海拔、坡度、降水、植被覆盖度四项环境因子的不同梯度对植被固碳-土壤侵蚀之间的空间权衡关系表现出不同程度的影响, 其空间权衡度的大小随因子的梯度变化表现出一定的规律(图 5)。在海拔因子的梯度下, 植被固碳与土壤侵蚀之间的权衡度仍以1000m为分界线, 在800—1000m的低海拔区域内其权衡关系不显著；在海拔大于1000m的中海拔地区, 权衡度随着海拔的升高而逐渐增大, 最大值为0.451。在坡度梯度下, 植被固碳与土壤侵蚀之间的空间权衡度随着坡度的升高而逐渐增大, 最大值出现在坡度大于35°的地区, 其值为0.62。降水因子对植被固碳-土壤侵蚀之间的影响较为显著, 五个降水等级下其权衡度在0.195—0.297的范围内波动。植被固碳与土壤侵蚀之间的空间权衡度在植被覆盖度为0.6—0.8的区间内达到最大值0.1, 在植被覆盖度为0—0.4和0.8—1的区间内, 权衡度并不显著。

3.3 不同地貌形态类型区生态系统服务的空间权衡度

受地貌的宏观控制作用以及地貌形态特征的影响, 水源涵养-土壤侵蚀、植被固碳-土壤侵蚀的权衡度在不同地貌类型区表现出不同程度的敏感性, 其值随地形起伏度的变化表现出一定的规律(表 2)。具体表现为：水源涵养与土壤侵蚀之间的权衡度在相对平坦的中海拔平原地区达到最大值0.169, 该值是其余四类起伏度较大的地貌形态类型区权衡度的十倍及以上。植被固碳与土壤侵蚀之间权衡度的最大值出现在中起伏中山地区, 其值为0.334。随着不同地貌形态类型区地形起伏度的变化, 植被固碳与土壤侵蚀之间的权衡度表现出随地形起伏度升高而逐渐增大的趋势, 具体表现为权衡度在各地貌形态类型区的排序：中海拔平原 < 中海拔台地 < 中海拔丘陵 < 小起伏中山 < 中起伏中山。

4 讨论与结论 4.1 讨论

在当前喀斯特地区运用RUSLE模型进行土壤侵蚀的研究中, 模型中各个因子的计算大多采用全国性的普遍参数, 未能将岩溶地区和非岩溶地区区别对待, 并不能完全反映喀斯特地区特殊的石漠化现象和薄土层地区无土可流的情况^[47]。同时, Biesemans等^[48]的研究发现, LS因子和C因子对土壤侵蚀最为敏感, 其计算的精确性对RUSLE模型整体有效性的影响最为显著。因此本文考虑了喀斯特地区特殊的地质背景, 在RUSLE模型中引入不同程度基岩裸露率校正系数, 以9m高精度DEM数据计算L因子, 运用基于遥感影像反演的高分辨率NDVI数据(30m)计算C因子, 使得本文得到的土壤侵蚀量结果与Feng等^[10]在喀斯特小流域尺度上获得的土壤侵蚀模数较为相似, 与贵州省水利厅公布的2015年全省喀斯特区域土壤侵蚀模数较为一致。

本文融合均方根偏差法和地理探测器的优势, 将权衡关系的判断和空间信息的探测集为一体, 构建生态系统服务空间权衡度指标定量化研究水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系。结果表明, 在水源涵养与土壤侵蚀的空间权衡关系中, 植被覆盖度对水源涵养-土壤侵蚀之间空间权衡度的影响最为显著, 其权衡度随植被覆盖度的增大而大幅减小。原因在于植被根系及其凋落物具有较强的持水能力^[49], 高植被覆盖度可以增强土壤的抗蚀、渗透、蓄水能力, 因此植被覆盖度较高的地区土壤流失量较少, 水源涵养-土壤侵蚀之间的空间权衡度降低。同时, 这一结论从侧面印证了喀斯特地区石漠化治理等生态工程的实施对减小土壤侵蚀具有重要意义。海拔因子对水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间空间权衡度的影响在1000m以上的中海拔地区明显高于低海拔区域。此种现象可能是因为我国西南地区的海拔变幅较大, 随着海拔的升高, 气候、植被以及地形的综合性影响更为显著, 因此权衡关系的空间差异性较为明显。我国西南喀斯特地区地形复杂、山多坡陡, 森林一般分布在海拔较高, 坡度较陡的地区^[22]。这一结论与本文研究发现植被固碳-土壤侵蚀之间的空间权衡关系随海拔和坡度的升高而逐渐增大的情况相符合。地貌因素控制着中小尺度空间的水热资源分配, 直接决定生态系统服务的供给与维持^[42], 这一结论解释了本文得出的地貌形态类型对水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系具有宏观控制作用。

本文通过对模型中因子的修正提高了模型模拟的精度, 构建了生态系统服务空间权衡度指标, 定量研究水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系, 结果发现不同梯度下的环境因子与地貌形态对生态系统服务之间的空间权衡关系具有显著的影响。本文针对喀斯特地区特殊的地质条件对RUSLE模型进行了优化, 而目前陆面模式等其他模型在喀斯特地区应用过程中的不确定性仍需深入研究。此外, 基于多尺度的山地水土时空耦合及其效应在机理和驱动机制阐释方面仍是研究中的关键问题。今后应重视多个时间节点或长时间序列内生态系统服务权衡/协同关系的时空变化, 辨析不同尺度生态系统结构-过程-功能-服务的作用机制。

4.2 结论

本文基于对三岔河流域2013年到2017年五年平均土壤侵蚀、水源涵养以及植被固碳的优化模拟和空间分析, 在坡度、海拔、降水、植被覆盖度4项环境因子的梯度下运用空间权衡度指标阐释水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系及其随地貌形态类型变化的趋势。结论如下：

(1) 三岔河流域的土壤侵蚀模数为3.19t hm^-2 a^-1, 侵蚀等级属于微度, 占流域面积的79.95%；水源涵养量范围为118.77—1230.14mm, 其空间分布呈现出由北向南逐渐递减的趋势；年均植被NPP是459.13gC/m², NPP累积量呈现西北高、东南低的特征。

(2) 由于植被显著的保持水土能力, 植被覆盖度对水源涵养与土壤侵蚀之间空间权衡关系的影响最大, 其空间权衡度随植被覆盖度的增大而逐渐减小；水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡度在1000m以上的中海拔区域显著高于低海拔区域, 其中, 植被固碳与土壤侵蚀在中海拔区域内的空间权衡度是低海拔区域的4倍。在喀斯特山区, 由于地形因子对地表水土过程与植被功能的影响, 植被固碳和土壤侵蚀之间的空间权衡度随海拔和坡度的升高而逐渐增大, 其中, 陡坡地区的空间权衡度是缓坡地区的6倍。

(3) 地貌的宏观控制作用及其内部特征显著地影响水源涵养-土壤侵蚀、植被固碳-土壤侵蚀之间的空间权衡关系。随着不同地貌形态类型区地形起伏度的变化, 较平坦地区水源涵养与土壤侵蚀之间的空间权衡度相对山地丘陵区较大, 具体表现为较平坦地区的空间权衡度是起伏度较高的山地丘陵区的近两倍；植被固碳与土壤侵蚀之间的空间权衡度随地形起伏度的升高而逐渐增大。