生态系统服务是指人类直接或间接从生态系统中获取的各种惠益, 涵盖人类赖以生存的物质基础和环境条件。多种生态系统服务在复杂的影响因素作用下往往表现出此消彼长的权衡关系(trade-off)和相互促进的协同关系(synergy)^[1-2]。识别生态系统服务间的相互关系, 厘清其时空变化特征及驱动机制, 可协同优化多种服务, 避免在提高一种服务的同时损害其他服务, 对促进区域生态系统可持续管理、指导自然资源的合理开发、提高人类福祉具有重要意义^[3], 也是生态系统服务从理论研究走向管理实践的关键环节^[4-5]。

近年来, 生态系统服务及相互关系研究已成地理学、生态学及生态经济学等多学科的前沿领域^{[1, 6-7]}。国内外学者普遍采用统计分析法^[8-9]、空间制图法^[10-11]、情景模拟法^[12-13]及模型量化法^[14-15]等方法, 借助InVEST、ARIES、Envisio、EcoAIM等软件^{[6, 16]}, 对生态系统服权衡/协同关系的时空表现特征、尺度效应、形成机制及不确定性等方面进行研究^{[8, 17]}。结果表明, 在气候、土地利用、社会偏好、激励政策等因素^[18-20]的影响下, 生态系统服务权衡与协同关系表现为空间异质性与时间动态性, 并且具有尺度依赖性^[21]和研究不确定性^[22]。然而, 当前研究尚处于起步阶段^[23], 主要集中于供给服务和调节服务之间的权衡关系, 且缺乏精确的数量模型, 急需通过增加数据样本、提高数据质量、优化研究方法等方式, 进行生态系统服务关系的度量。同时, 驱动机制的研究也较薄弱, 仅局限于一种或几种特定影响因素分析^{[12, 24]}, 缺少对驱动力的综合解析及主导驱动力识别, 从而限制生态系统服务关系研究在管理实践中的应用。

流域作为复合生态系统, 具有多样的服务类型, 在自然和人为因素的共同作用下, 生态系统服务之间表现出复杂的时空关系特征。同时流域也是人类活动干扰最强烈的地区之一, 面临水质退化、水量不足、水土流失等问题, 如何协调经济发展与环境保护的关系、缓解上下游利益相关者之间的矛盾成为流域可持续发展的关键。生态系统服务是衔接自然环境与人类需求的重要纽带, 通过研究流域多种生态系统服务的权衡与协同关系, 明确其驱动机制, 有助于制定流域发展与生态保护“双赢”的政策措施^{[1, 25]}。基于此, 本文以太湖流域江苏省为例, 考虑到农业非点源污染是该区域的首要水环境问题, 选择氮、磷净化作为主要服务, 其次是水量供给与土壤保持服务；通过空间分布式的生物物理模型与经验统计模型, 计算研究区每种生态系统服务的数量值；在此基础上, 借助GIS空间分析和数理统计法, 研究2种水体净化服务与其他服务之间的权衡/协同关系。最后, 结合驱动因子空间栅格化, 利用Logistic回归模型, 定量分析影响生态系统服务关系的主导驱动力, 以期为流域的水环境管理及生态环境保护奠定基础。

1 数据与方法 1.1 研究区概况

太湖流域地处长江三角洲城市群, 是中国经济最发达的地区之一。研究区选择太湖流域的江苏省(图 1), 因其位于相对完整的行政区和独立的水利分区^[26], 便于研究数据收集及生态环境模拟。该区北靠长江, 地势西高东低, 平均气温15—17 ℃, 多年平均降雨量约为1181 mm, 气候温暖湿润, 是中国重要的粮食产区；总面积为16984.48 km², 包括苏州市、无锡市、常州市、镇江市及南京市的部分地区, 城市化水平较高。近年来, 随着社会经济发展, 区域人地矛盾加剧, 生态环境日益退化, 非点源污染尤为严重。研究与非点源污染物产生、排放、传输、净化有关的生态系统服务类型, 明确其相互关系的时空特征及驱动机制, 从生态系统综合管理的角度为区域水环境管理, 尤其是非点源治理提供新思路。

1.2 数据来源

生态系统服务及驱动力空间量化所需的数据包括降雨、温度、太阳辐射等气象数据, 土壤类型、理化性质、含水量、深度等土壤数据, 用于植被覆盖度计算的遥感数据, 土地利用及道路的空间数据, 人口密度、GDP、农业等社会经济数据, 以及用于模型校正的水文、水质数据等, 具体来源如表 1所示。

1.3 研究方法 1.3.1 水体净化服务模型

生态系统的植被、土壤等要素通过储存、转换等方式移除径流中的营养盐污染物, 从而减少排放到外界的污染物, 达到净化水体的作用。生态系统净化功能越强, 营养盐输出量越小。因此, 利用氮输出量(Nitrogen Export, NE)和磷输出量(Phosphorus Export, PE)表征生态系统氮、磷净化服务, 这两个反向指标值越大表明水体净化服务越低, 反之亦然。本文利用InVEST软件的“营养盐净化”模型, 结合经验统计模型、专家咨询及文献调研等方法修正参数, 以实现对研究区2种水体净化服务的空间量化。由于该模型考虑了不同地类氮、磷营养盐参数的差异性, 在一定程度上降低了服务指标之间的相关性, 可用于后续的相互关系分析。模型原理如下：

(1)

式中, Exp_x为上游栅格x的氮、磷输出值(kg/hm²)；E_y表示每个下游栅格y过滤氮、磷营养盐的效率(%)；X代表氮、磷输出量从上游产生到下游水体的运输路径；ALV_x为校正的栅格x的氮、磷输出值(kg/hm²), 计算公式如下：

(2)

(3)

(4)

式中, pol_x表示栅格x的氮、磷营养盐输出系数, HSS_x是水文敏感分值, λ_x为栅格x的径流系数, λ_w则是研究区平均径流系数；的总产水量, 借助InVEST平台的“水产量”模型进行计算。

1.3.2 水量供给服务模型

生态系统通过水循环产生可供人类蓄水、发电等利用的水量供给服务(Water Supply, WS), 并且水循环过程会溶解、析出氮、磷营养盐, 进而间接影响区域非点源污染。本文采用InVEST软件的“水产量”模型来计算研究区的地表产水量, 实现水量供给服务的空间化。模型原理如下：

(5)

式中, WS(x)为栅格x的年产水量(mm)；Pre_x为栅格x的年降雨量(mm)；AET_x为栅格x的实际蒸散量(mm)。该公式可进一步表示为：

(6)

式中, AET_x/Pre_x为基于Budyko曲线的表达方式, 其计算公式如下：

(7)

式中, ω_x为非物理参数, R_x为栅格x的干燥指数, 各自计算公式如下：

(8)

(9)

式中, Z为降雨时间分布和水文生态特征的经验常数；AWC_x为栅格x的植被有效含水量(mm), 由公式(10)求得。k_j为不同地类的蒸散系数；ET₀为潜在蒸散量(mm/d), 通过Modified-Hargreaves法进行计算。

(10)

式中, soil.depth为限制根系生长的土壤最大深度(mm), root.depth为根系长度(mm), PAWC为植物有效持水量(%)。

1.3.3 土壤保持服务模型

生态系统在预防土壤侵蚀和保持土壤过滤污染物功能等方面发挥重要作用, 进而影响区域非点源污染。本文使用RUSLE模型估算土壤保持服务(Soil Retention, SR), 并基于小流域的实验数据对其进行校正, 以得到本地化的计算模型, 如下：

(11)

式中, R、K、LS、C分别为降雨侵蚀力因子、土壤可蚀性因子、坡长坡度因子、植被覆盖因子, 计算公式如下所示。P表示工程管理措施因子, 为0—1的无量纲数, 各地类赋值参考吕建树等的研究^[27]。A为本地化修正系数, 用来调整土地利用类型和地表光滑度对土壤侵蚀的影响, 根据研究区实测数据确定值为1.3。

(12)

式中, p为年平均降雨量(mm), p_i为月降雨量(mm)。

(13)

式中, S_a、S_i、S_n分别为砂粒、粉砂与粘粒的含量值(%)；C表示有机碳含量值(%)。

(14)

式中, F_a为汇水累积阈值, C_s为栅格大小, s为坡度, n为坡度相关参数, β为流向相关参数。

(15)

式中, C为0—1的无量纲数, c为植被覆盖度, 通过归一化植被指数(NDVI)求得。

1.3.4 驱动力分析方法

考虑研究数据的可得性, 在多因子相关分析的基础上, 结合专家咨询方法, 最终确定气候、土壤、经济、用地等10类31种驱动力, 具体如表 2所示。

在生态系统服务关系和驱动力空间量化的基础上, 将各自栅格值赋予GIS生成的10000个随机点, 形成生态系统服务关系为因变量, 驱动力为自变量的数据集, 再运用多元Logistic回归模型进行驱动机制研究。模型原理如下：

(16)

式中, p_i为生态系统服务权衡/协同关系发生的概率, x_i为驱动力；α为截距, 即模型的常数项；β_i为斜率, 即模型中每个驱动力的偏回归系数, 其值为正则表示因变量与自变量的变化方向相同, 负值则相反。

为更好的测量变量之间的关联性, 明确自变量对因变量发生概率的作用程度, 一般应用p/(1-p)来解释回归系数, 即发生比率(odds ratio, OR), 其计算公式如下：

(17)

本文选择Homsmer-Lemeshow(HL)指标进行模型拟合度检验, 当HL值不显著表示模型拟合好, 而显著则说明模型拟合不好。通过Wald统计量检验模型的回归系数, 取0.05为显著水平；当自变量的Wald所对应的P值小于0.05时, 保留在方程中, 而大于0.05则因未通过显著性检验而被去除。

2 结果分析 2.1 生态系统服务的时空变化特征

本文计算了2000、2005及2010年研究区4种生态系统服务指标的物理量, 如图 2。其中, 绿、黄、红色分别表示服务供给的高、中、低值区, 蓝色代表水体。结果显示, 氮净化服务高值区(氮输出量低值区)主要分布在研究区南部、西部及湖体周围, 低值区(氮输出量高值区)则零散分布于北部及中部, 中值区分布广泛。从2000到2010年, 氮输出量单位面积年均值呈现先增加后略有下降的趋势, 增幅由8.89%降低为8.08%。磷净化服务高值区(磷输出量低值区)主要分布于研究区周边的林地、草地、湿地等生态用地, 低值区(磷输出量高值区)在农村居民点及城镇建设用地附近突出, 并零散分布于部分农耕地。2000—2010年, 磷输出量年均值逐渐增加, 增幅为24.69%。过量的化肥施用、畜禽养殖等农业活动会增加氮、磷的排放, 而农村生活污染和城镇污水也会显著影响磷素富集^[28]。因此磷净化服务分布与农村居民点、城镇建设用地存在明显关系, 进而导致氮、磷指标时空变化的差异性。已有研究也表明磷负荷与城市用地、农村建设用地密切相关^[29-30]。水量供给服务的空间分布呈东西走向, 且变化显著, 高值区主要分布于研究区周边, 中值和低值区由研究区西部逐渐转向中、东部。2000—2010年, 水量供给服务年均值呈先下降后上升的趋势, 降幅由-27.88%变为-12.76%。研究区土壤保持服务普遍较低, 低值区分布广泛, 高值区零星分布。2000—2010年, 土壤保持服务年均值呈下降趋势, 降幅为-21.61%, 主要受降雨、日照等自然因子^[31]以及农田、林地、城镇等用地变化^[32-33]的影响。

2.2 生态系统服务权衡与协同关系分析

在ArcGIS软件中将2010年和2000年每种生态系统服务相减, 得到各自的空间变化图, 并将变化分为正向、负向及未变化3类。利用GIS的空间叠加分析功能得到两种服务变化的空间分布格局(图 3)。其中, 正向和负向变化区为两种服务均增加或均减少的区域, 表明两者具有相同变化方向的协同关系；相反变化区为两种服务变化方向相反, 表明两者具有权衡关系；未变化/正(负)向变化区为一种服务没变化而另一种服务增加(减少), 表明两者关系不显著；未变化区是两种服务均未变化, 其关系也不显著。统计每类变化区占区域总面积的百分比, 结果见表 3。本文以水体净化服务为主, 分析该类服务指标与其他服务的相互关系。

结果表明, 氮输出指标与水量供给关系中, 相同变化区(包括正、负向变化区, 表示氮净化与水量供给服务的权衡关系区)面积占统计面积的33.16%, 其中正向变化区占11.97%, 主要分布于研究区西北部、西南角及东南角地区, 负向变化区占21.19%, 与城镇建设用地分布较一致；相反变化区(表示氮净化与水量供给服务的协同关系区)占统计面积的29.98%, 主要分布于研究区中、东部。磷输出指标与水量供给关系中, 相同变化区和相反变化区的面积占比分别为27.21%和31.92%, 在研究区西北部、西南角及中部地区集中分布；其次是未变化/正向变化区, 占统计面积的26.43%, 主要位于研究区的中西部及北部一段沿江地区。氮、磷输出指标与土壤保持关系中, 未变化/正向变化区分布范围均最广, 分别占统计面积的48.06%和55.37%, 而正向、负向及相反变化区所占比例较小, 说明水体净化服务与土壤保持服务的权衡/协同关系不明显。氮与磷指标关系中, 正向变化区广泛分布, 占统计面积的58.86%, 说明研究区水体净化服务主要表现为协同关系。

2.3 生态系统服务关系的驱动力分析

选择显著的生态系统服务关系进行驱动力分析, 主要为氮净化服务与水量供给服务的权衡关系、两种水体净化服务的协同关系。通过多元Logistic回归模型进行定量分析, 模型中因变量是两种服务关系, 取值为0和1值, 0表示两者关系不显著, 即一种变化而另一种未变化或者两者均未变化；1表示两者具有协同关系(均正向变化或负向变化)或权衡关系(相反变化)。据此, 将生态系统服务关系的空间格局图进行重分类, 得到0和1的因变量图层。同时, 在多期驱动因子空间栅格化的基础上计算其平均值, 得到自变量图层。将因变量和自变量的图层值赋予随机点, 基于点的属性值确定Logistic回归模型的输入值。由于数据来源于点的随机抽样, 避免了空间自相关性；同时在随机点选取时, 尽量保证因变量的0和1值数量相等, 以避免对模型常数项的影响。模拟结果如表 4、5所示。

表 4模型的HL指标值为0.817, 统计不显著, 反映模型的拟合度较好；模型预测正确率为81.32%, 表明模型较稳定。根据Wald统计量可判断, 影响氮净化和水量供给服务关系的最重要因素是城镇建设用地密度和植被覆盖度, 其次是离农村距离、农村居民点密度、水网密度、离城镇距离及人口密度等。其中, 城镇建设用地密度和农村居民点密度的回归系数为正, 说明每平方公里增加1hm²的城镇或农村建设用地, 将使氮净化和水量供给服务关系权衡变化关系的发生概率分别增加29.100倍和10.282倍。而植被覆盖度和水网密度的回归系数为负, 表明此权衡关系将随着这些因素的增加而减小。此外, 年相对湿度、年均气温等气候因素的Wald统计变量也比较显著, 说明气候因素对两者权衡关系也具有一定的影响, 尤其是年平均气温对其贡献率达到28.979。土壤因素中, 总磷含量是最重要的解释因变量, 且其发生比率也最大, 达到37.490。社会经济因素中, 除人口密度的贡献率最大外, 农业GDP也是较为重要的影响因素。农业GDP增加可视为当地农业发展水平提高, 此时将带来两者权衡关系的下降, 其发生概率为0.997倍。因此, 当气候变化带来区域水量供给增加时, 由于农业的发展, 可采取适应性措施, 减弱对区域生态系统水体净化服务的影响。

表 5模型的HL指标值为1.017, 统计不显著, 反映模型拟合度较好, 同时模型也较稳定, 其预测正确率为89.49%。根据Wald统计量可判断, 影响两种水体净化服务协同关系的最重要驱动力是植被覆盖度, 其次是耕地比例、林地比例、离农村距离、水网密度, 并且这些解释变量的回归系数均为正, 表明这些因素的增加将带来两种服务协同关系的提高。从发生比率来看, 驱动因素的贡献率均不大, 普遍在1.000左右, 仅有农村居民点密度的概率达2.406。总体来看, 所有驱动力中植被用地的影响最明显, 通过增加耕地、林地及草地的面积或建立植被缓冲带将显著促进两种水体净化服务的协同变化；而增加林地、草地等生态用地可有效提高生态系统功能, 净化多种水质污染物, 达到协同减排的目的。此外, 土壤砂粒、粘粒等颗粒物组成也会促进两者的协同关系。因此, 通过调节土壤属性或改良土壤结构可同时改变两种水体净化服务, 为区域非点源污染物的协同治理提供有效途径。

3 结论与讨论

复杂的内外界环境因素作用于生态系统的结构、过程及功能, 进而影响生态系统服务供给, 表现为多种服务在时空尺度上的同增同减、一增一减或无明显变化关系。生态系统服务间相互关系研究有助于管理者识别多种服务的变化特征, 明确区域生态系统的变化程度, 以采取生态空间管制、环境污染防治等针对性政策措施进行调控, 促进区域生态系统的健康、持续发展。流域是自然环境和人类活动交互作用强烈地区, 面临诸多水环境问题, 结合流域的实际情况, 选择关键的生态系统服务, 进行多种服务之间的权衡与协同关系研究, 为实现流域生态系统可持续管理奠定基础。

本文以太湖流域江苏省为研究区, 运用空间显示的生物物理模型和经验统计模型计算区域4种生态系统服务指标, 在此基础上借助GIS空间分析手段进行服务权衡与协同关系研究, 并利用多元Logistic模型定量识别影响生态系统服务关系的主导驱动力。研究结果表明：

(1) 从2000到2010年, 研究区4种生态系统服务指标的单位面积年均值呈现不同程度的变化, 氮输出指标先增加后略有下降的趋势, 磷输出指标逐渐增加, 水量供给先下降后上升, 土壤保持则呈现逐渐下降趋势。研究区4种服务指标的空间分布及变化格局具有较明显的空间异质性。

(2) 氮、磷输出指标与水量供给关系的空间分布格局较为相似, 表现为广泛分布的相同变化区和相反变化区；氮、磷输出指标与土壤保持的正向、负向及相反变化区占比均较小, 表明水体净化服务与土壤保持服务的权衡/协同关系不明显；氮、磷输出指标的正向变化区广泛分布, 说明研究区水体净化服务主要表现为协同关系。

(3) 氮净化和水量供给服务权衡关系的正向驱动力为城镇建设用地密度和农村居民点密度, 且每平方公里增加1公顷的城镇或农村建设用地, 其权衡关系的发生概率将增加29.100倍8和10.282倍；而植被覆盖度和水网密度则具有显著的负作用。氮、磷服务协同关系的首要影响因素是植被覆盖度, 其次是耕地、林地比例等因素, 且这些因子均表现为正向促进作用。

尽管生态系统服务的研究已逐步深入, 然而从理论走向实践仍面临一系列挑战^{[1, 3, 5]}, 尤其是如何更精准的计算服务量, 如何科学识别多种服务之间相互关系, 如何定量分析服务关系的驱动机制等问题。本研究综合利用生态、地理、统计、计量等方法手段, 从时空多维度进行生态系统服务关系的表征及驱动力分析, 以期拓展相关研究, 为进一步理解生态系统服务之间复杂关系及其驱动机制提供借鉴, 并促进生态系统服务研究在流域水环境管理中的应用。生态系统服务具有时间积累性, 而不同类型服务的变化速度和阈值也存在差异, 影响服务关系的时间表现特征。由于数据的有限性, 本文仅进行了2000—2010年的研究, 未能充分考虑每种生态系统服务变化的周期。未来可结合生态系统变化的阈值分析进一步探究不同时间周期及更长时间范围上生态系统服务之间的相互关系, 并探索驱动力变化程度对生态系统服务关系的影响。