自2005年千年生态系统评估发表以来, 生态系统服务(Ecosystem Services, ES)研究已逐渐受到学术界的关注^[1—2]。其中, 探讨和分析不同生态系统服务之间的权衡与协同作用已成为重要方向之一^[3—5]；且如何辨识生态系统服务之间多重非线性关系及其影响因素的研究已成为该方向的研究热点之一, 它对实现生态系统管理可持续发展和改善人类福祉, 具有重要的理论和现实意义。

目前, 有关生态系统服务之间的相关关系, 即权衡与协同关系及其驱动因子分析已取得较大进展。研究大多采用相关分析^[6—9]、空间制图^[10—13]、情景分析^[14—16]、最小二乘回归和地理加权回归^[17]等方法, 虽可直观地揭示生态系统服务间的相关性, 但在反映生态系统服务间的非线性关系方面略有不足, 在准确表达生态系统服务间权衡协同关系的内在机理及其驱动因素方面还有较大提升空间。基于贝叶斯网络模型在表达多要素间非线性关系的优势, 已有学者采用该模型探讨单一生态系统服务与多因素间的相关性研究。比较有代表性的如李晶等^[18]和Dang等^[19], 前者利用贝叶斯网络模拟水源涵养服务过程, 根据设置的不同土地利用情景预测2050年的水源涵养服务状态概率, 对关天经济区水源涵养服务空间格局进行优化；后者基于贝叶斯网络模型探讨环境和人为因素对水稻供给服务产生的影响, 优化水稻供给能力并为农业的可持续发展提供决策支持。上述研究主要是基于贝叶斯网络探讨多种因素对单一服务的影响, 对于多种ES权衡协同及其驱动因素效应差异的研究较少。即在一定区域范围内, 如何构建基于多种生态系统服务的贝叶斯网络模型、辨识生态系统服务间的多重非线性关系及其驱动因子, 正是本文的研究重点和要解决的关键问题。

基于此, 本文以陕北黄土高原为研究对象, 在定量评估其2018年土壤保持、产水服务、粮食供给和固碳服务的基础上, 构建基于贝叶斯网络的生态系统服务模型, 利用联合概率分布分析4种ES的权衡协同关系, 最后探讨影响ES权衡协同关系的驱动因子及其效应差异。研究对于提升对多种生态系统服务复杂关系的认识、提高研究区生态系统的整体效益均具有重要意义。

1 研究区概况和数据来源 1.1 研究区概况

陕北黄土高原位于中国黄土高原地区的中心地带, 陕西省北部, 地理坐标35°16′—39°34′ N, 107°15′—111°15′ E。地势西北高、东南低, 基本地貌类型为黄土塬、梁、峁、沟、塬。地域南北长、东西窄, 总面积为92521.4 km², 占全部黄土高原总面积的12.6%(图 1)。气候类型为温带大陆性半干旱气候, 年均温为8—12 ℃, 年降水为350—600 mm。区域土地利用类型以耕地、林地和草地为主, 呈西北耕地向东南林地和草地过渡的空间格局。该区行政区划包括榆林市和延安市, 下辖府谷、神木和榆阳等25个县(市/区), 2018年总人口达767.72万人, 人口密度为83人/km²。

作为高度敏感的生态脆弱区之一, 陕北黄土高原一直是政府和学者关注的焦点区域^[20—21], 连续15年的退耕, 使得该区各类ES变化明显。加之陕西乡村振兴战略的逐步落实, 势必对该区ES的供给、需求和调节产生多样性影响。本文以陕北黄土高原为案例地, 为分析区域生态系统服务权衡协同关系及其驱动因子提供了良好的研究平台。

1.2 数据来源与处理

研究数据主要包括：(1)气象数据：气温、降水空间插值数据, 分辨率1 km, 来源于国家地球系统科学数据中心(http://www.geodata.cn)；(2)2018年土地利用数据：分辨率30 m, 来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)；(3)NDVI、NPP数据：分辨率为1 km、500 m, 来源于美国国家航天航空局(NASA)提供的数据(https://search.earthdata.nasa.gov/search)；(4)DEM数据：分辨率90 m, 来源于地理空间数据云平台(http://www.gscloud.cn)；(5)人口密度数据：分辨率1 km, 来源于世界人口数据集(http://www.worldpop.org)；(6)土壤数据：采用基于世界土壤数据库(HWSD)的中国土壤数据集, 比例尺为1∶100万, 来源于寒区旱区科学数据中心(http://westdc.westgis.ac.cn)。以上数据均在ArcGIS 10.6平台上进行数据预处理, 将所有栅格数据重新采样, 空间分辨率为1 km。此外, 数据均采用Albers投影。

2 研究方法 2.1 生态系统服务量化

结合研究区水资源短缺、水土流失严重、生态环境脆弱等地理特征和资源环境状况, 本文选择土壤保持、产水服务、粮食供给和固碳服务等评估其生态系统功能, 并采用InVEST模型、ArcGIS等方法对上述服务进行量化。

2.1.1 土壤保持

利用InVEST模型的泥沙输移比模块对土壤保持服务进行评估, 该模块以通用的水土流失方程为基础, 考虑上游水源截留沉积的泥沙量。计算公式：

(1)

式中, SR为土壤侵蚀总量(t)；R为降雨侵蚀性因子(MJ mm hm^-2h^-1a^-1), 采用Wischmeier和Mannering^[22]提出的公式计算；K为土壤可蚀性因子(t hm²h MJ hm^-2mm^-1), LS为坡度坡长因子, 无量纲；C为植被覆盖和作物管理因子, 无量纲, 参照蔡崇法等^[23]的研究成果；P为水土保持措施因子, 无量纲。

2.1.2 产水服务

利用InVEST模型中的产水模块计算产水量, 计算公式：

(2)

式中, Y(x)、AET(x)、P(x)分别为栅格单元x的年产水量(mm)、年实际蒸发量(mm)、年降水量(mm), 具体计算步骤详见包玉斌等^[24]的研究。

2.1.3 固碳服务

利用InVEST模型的碳储量模块对固碳服务总量进行计算, 模块包括地上物质、地下物质、土壤和死亡有机物碳储量, 求其和获得总碳储量。计算公式：

(3)

式中, C_tot为总碳储量(t), C_above、C_below分别为地上、地下生物量中的碳量(t/hm²), C_soil为土壤中的碳量(t/hm²), C_dead为死亡有机物中的碳量(t/hm²), 各部分碳储量参数根据文献资料获得^[25]。

2.1.4 粮食供给

粮食供给服务作为生态系统服务中最为基础的供给服务之一, 对于人类生存与发展有着不可或缺的作用。研究表明粮食供给与NPP之间具有很强的相关性^[26], 本文将粮食的总产量按照栅格NPP值与耕地总NPP值的比值来分配, 确定各栅格的粮食供给能力, 计算公式：

(4)

式中, Grain_i和Grain_j分别是网格i和县j的粮食产量(kg)；NPP_i是网格i耕地的净初级生产力, NPP_j是县j耕地的净初级生产力之和(kgC/m²)。

2.2 基于贝叶斯网络的生态系统服务模型 2.2.1 模型构建

贝叶斯网络(Bayesian belief network, BBN)由Pearl^[27]于1988年提出, 是一种基于概率推理的图形网络。BBN的构建包括两个步骤：一是网络结构学习, 确定贝叶斯概念网络的节点和连接节点的箭头。其中, 节点表示变量的概率分布, 箭头表示变量之间的相互依赖关系；二是网络参数学习, 确定各节点的条件概率表, 反映节点之间因果关系的强弱^[28—30]。

本文使用Netica软件构建基于贝叶斯网络的生态系统服务模型(简写为BBN-ES模型)。首先, 根据生态系统服务的过程机理, 选择各生态过程中的相关变量作为节点, 对贝叶斯网络进行结构学习。构建的贝叶斯概念网络共包含14个节点和17个箭头。其中, 土壤保持、产水服务、粮食供给和固碳服务等4种生态系统服务为目标节点, 人口密度、坡度、降水、土地利用、气温、NDVI、土壤可蚀性、降雨侵蚀、蒸散发、NPP等10个变量为影响因子节点。其次, 利用ArcGIS软件对数据图层进行离散化处理, 具体变量离散化分级标准(表 1), 土地利用节点为离散变量, 根据中科院土地利用覆盖分类体系划分为6类；其它节点为连续变量, 利用自然断点法划分为低、较低、中、较高、高5种状态。并创建1 km×1 km渔网提取各图层栅格值, 对贝叶斯网络进行参数学习。最后, 将提取的78366个样本数据带入网络进行样本训练, 获得贝叶斯网络所有节点的条件概率表, 即BBN-ES模型构建完成(图 2)。

2.2.2 模型验证

本文利用误差矩阵评估模型精度。误差矩阵是通过比较分类结果的预测值与实际值, 即将正确分类的样本数除以样本总数得到分类结果的精度^[31]。误差矩阵的每一列代表预测类别, 列总和表示预测类别数目；每一行代表实际类别, 行总和表示实际类别数目。通过在研究区创建3 km×3 km渔网生成采样点作为验证数据, 共8962个。利用Netica软件对4个目标节点进行准确性测试, 计算其误差矩阵评估模型精度。

2.3 ES权衡协同关系的驱动因子分析 2.3.1 ES的权衡协同关系

利用贝叶斯网络的联合概率分布分析ES间的相关关系。其中, 联合概率由贝叶斯网络所有节点的条件概率分布确定, 计算公式如下：

(5)

式中, P(X₁, X₂, …, X_n)表示一个联合概率分布, P(X_i|parents(X_i)表示条件概率分布。在联合概率分布图中, ES间的相互作用可视化为成对的联合概率分布, 每个像素代表服务1和服务2发生的联合概率P(ES₁=x, ES₂=y), 由ES₁和ES₂在各自状态下的条件概率x和y计算得到。像素颜色的深浅程度代表两服务间的联合概率大小, 颜色越深, 表明概率越高, 一个联合概率分布内的像素值总和为1。若两种服务的联合概率分布呈“I”型分布, 则表示这两种服务为协同关系；若两种服务的联合概率分布呈“L”型分布, 则为权衡关系^[32]。

2.3.2 影响因子对生态系统服务的重要性分析

基于Netica软件敏感性分析, 通过计算方差缩减探讨BBN-ES模型中影响因子节点对生态系统服务节点的相对重要性大小^[33]。计算公式：

(6)

式中, VR为方差缩减, 表示相对重要性大小；V(Q)、E(Q)分别为生态系统服务Q的方差和期望；V(Q|F)、E(Q|F)为变量F条件下生态系统服务Q的方差和期望；X_q是与状态q对应的真实数值。VR值越大表明该节点对目标节点的相对重要性大, 即其对ES节点的影响越大；反之表明对ES节点的影响越小。本文选取VR>0.5%的影响因子作为影响生态系统服务的关键节点。

2.3.3 驱动因子分析

本文借鉴Feng等^[34]的研究, 根据生态系统服务的权衡协同关系和节点的重要性分析, 通过设置BBN-ES模型中ES节点的不同状态概率值作为情景, 利用贝叶斯网络的概率推理, 计算4种情景下影响生态系统服务关键节点状态的后验概率相对先验概率的变化, 并依据研究区的实际情况, 将概率变化最大值大于35%的关键节点确定为影响ES权衡协同关系的驱动因子。

依据ES的相关关系, 设置如下4种情景：

情景Ⅰ：将具有协同关系的A服务和B服务最大化, 即设置其“高=100%”；

情景Ⅱ：将具有协同关系的A服务和B服务最小化, 即设置其“低=100%”；

情景Ⅲ：将具有权衡关系的A服务最大化、B服务最小化, 即设置A服务“高=100%”、B服务“低=100%”；

情景Ⅳ：将具有权衡关系的A服务最小化、B服务最大化, 即设置A服务“低=100%”、B服务“高=100%”。

3 结果与分析 3.1 生态系统服务的空间分布

图 3为2018年研究区的土壤保持、产水服务、粮食供给和固碳服务的空间分布图。由图 3可知, 4种生态系统服务均呈现出较强的空间异质性。其中, 土壤保持、产水服务和固碳服务的空间分布较为一致, 呈现出“南高北低”的特征；而粮食供给则呈现出“中部高、南北低”的空间分布。

结合表 1和图 3可知：土壤保持服务的高值主要分布在研究区的东南部、中部和东北部区域, 呈带状分布且与河流的走向保持一致, 所占面积2183.5 km²；低值则在北部和西北部聚集成片, 面积为66245.3 km²。土壤保持服务区域差异的原因在于高值区所在的延长、宜川和黄龙等县植被覆盖度高, 而低值区所在的神木、榆阳和横山等市(区、县)建设用地占比较大, 易发生水土流失。

产水服务的高值集中在研究区的南部区域, 且北部呈现高值、较高值交错分布, 高值区所占面积2359.3 km²；低值分布在研究区的西部和中东部, 中南部有零星分布, 面积为48851.3 km²。高值区所在的黄陵、宜川和黄龙等县区域林草地分布广泛, 有利于产水服务；而低值区所在的定边、靖边和横山等县植被覆盖度较低。

粮食供给的高值集中在研究区中西和中东部区域；低值分布范围较为广泛, 分布在研究区东北部和南部, 面积为62451.9 km²。这是由于高值区所在的定边、横山和子洲等县域的耕地面积大, 粮食供给能力较强；而低值区所在的富县和宜川县等区域林草地分布较广, 耕地面积相对较小, 粮食供给能力相对较弱。

固碳服务的高值区和低值区的分布与产水服务基本一致, 低值集中在西北部外, 在研究区中部有零星分布, 且固碳服务的中值集中在研究区东北部、中西部和中南部, 所占面积49961.6 km²；高值和低值分布的原因与产水服务类似。

3.2 BBN-ES模型及验证

构建的2018年研究区4种生态系统服务的贝叶斯网络如图 4所示。经过参数学习后的贝叶斯网络反映了研究区2018年真实情况的先验概率。其中, 土壤保持服务由低到高的5种状态概率分别为71.6%、18.4%、5.01%、2.67%和2.36%；产水服务5种状态概率分别为52.8%、15.6%、9.04%、20%和2.55%；粮食供给服务各状态概率为67.5%、2.33%、25.2%、3.21%和1.81%；固碳服务各状态概率分别是22.5%、6%、54%、10.9%和6.59%。土壤保持、产水服务和粮食供给低状态概率均超过50%, 表明研究区控制侵蚀和拦截泥沙的能力较弱、总体产水服务水平较低、粮食供给服务明显不足；而固碳服务22.5%(< 50%)的地区处于低水平, 其在中等及以上的地区占71.49%, 表明研究区总体固碳服务水平较高, 主要是退耕还林(草)工程实施后, 研究区的林草地面积增加, 植被覆盖率的增强有利于固碳服务。

模型验证以土壤保持服务为例, 预测土壤保持状态与实际土壤保持状态误差矩阵见表 2, 土壤保持精度为74.4%；表中第一行第一列表示有5970个实际为低状态的样本被正确预测, 第二列的100表示有100个低状态样本被错误预测为较低状态。同理, 计算产水服务、粮食供给、固碳服务精度分别为70.0%、75.9%、73.5%, BBN-ES模型总体精度为73.5%, 表明模型具有较高的精度, 模拟效果好, 对上述4种ES节点概率预测具有较强的可靠性。

3.3 ES权衡协同关系的驱动因子分析 3.3.1 ES权衡协同关系分析

基于BBN-ES模型, 利用联合概率分布和相关系数得到土壤保持、产水服务、粮食供给和固碳服务等4种服务两两间的相关关系(图 5)。据图可知, 2018年研究区的土壤保持、产水服务与固碳服务两两间呈“I”型概率分布, 互为协同关系；其中, 产水服务与固碳服务的协同关系最强(r=0.52)。究其原因, 研究区自1999年实施退耕还林(草)工程后, 截止2018年其林草地面积迅速增加, 林地增加12.2%, 草地增加3.5%, 林地将吸收的土壤水分凝结成靠近林草地的降雨, 有利于产水服务^[35]；同时, 1999—2018年研究区年均降水量为183.2—432.6 mm, 而2018年降水量为325.7—475.1 mm, 有明显的增加趋势, 这也会提高其产水服务。且植被覆盖度的增加也会提高固碳能力, 同时也会减少雨水对土壤的侵蚀增加土壤保持能力。而粮食供给分别与土壤保持、产水服务、固碳服务呈“L”型概率分布, 即存在权衡关系；其中, 产水服务与粮食供给的权衡关系最强(r=-0.62)。结合图 1和图 3可知, 2018年研究区土地利用类型以耕地、林地、草地为主, 粮食供给能力较高的耕地区域, 其土壤保持、产水服务和固碳服务相对较弱；而土壤保持、产水服务和固碳服务较高的林草地区域, 粮食供给能力较弱。

3.3.2 影响因子对生态系统服务的重要性分析

影响因子节点对生态系统服务节点的重要性分析结果见表 3。依据2.3.2中设置判断关键节点的规则(VR>0.5%), 结合表 3可知, 土壤保持服务的关键节点包括坡度、土地利用、土壤可蚀性、降雨侵蚀、降水和NPP, 其中坡度的VR值最大(7.91%), 表明其对土壤保持的变化影响最大, 原因在于研究区实施退耕还林(草)工程后, 大量耕地转化为林地和草地, 植被覆盖度的提升增强了植被的固土持沙能力, 有利于土壤保持^[36]。影响产水服务的关键节点包括降水、降雨侵蚀、NPP和蒸散发, 且4个节点的VR值均在2%以上, 表明其对产水服务增加影响较大。对于粮食供给服务而言, 主要受到NPP、降水和降雨侵蚀等3个因子的影响, 其VR值均>10%。影响固碳服务的依次为NPP、降水、NDVI和降雨侵蚀, 其中NPP的VR值最大(4.83%), 表明其对固碳服务的变化影响最为显著, 原因在于NPP与区域的植被覆盖度密切相关, 植被覆盖度的增加有利于固碳服务。

因此, 确定坡度、土地利用、土壤可蚀性、降雨侵蚀、降水、NPP、蒸散发和NDVI为4种服务的8个关键节点, 为下文分析ES权衡协同关系的驱动因子奠定基础。

3.3.3 驱动因子分析

依据3.3.1权衡协同关系结果设置4种情景(表 4)。利用贝叶斯网络的概率推理, 计算4种情景下影响生态系统服务8个关键节点状态的后验概率相对先验概率的变化(表 5), 确定土地利用、降水和NPP是影响ES权衡协同关系的主要驱动因子。

由表 5可知：在协同-情景Ⅰ和情景Ⅱ中, 土地利用、降水、NPP等3个节点的概率变化最大值均大于35%；在权衡-情景Ⅲ和情景Ⅳ中, 土地利用、降水和NPP结点的概率变化最大值均大于35%, 由此确定土地利用、降水、NPP是影响ES权衡协同的驱动因子。其中, 土地利用主要影响生态系统服务的协同关系, 而降水和NPP在影响协同关系的同时对权衡关系也有一定地制约作用。

(1) 在协同情景中, 不同驱动因子在不同状态下产生相同的协同方式。如情景Ⅰ, 降水较高值、NPP高值状态下, 均可促使3种ES的高值协同；同时, 相同驱动因子在不同状态则会造成协同方式的差异。譬如, 降水在低值状态下促使3种ES的低值协同(情景Ⅱ), 而较高值状态下表现为高值协同。

(2) 在权衡情景中, 不同驱动因子在不同状态下产生相同的权衡方式。如情景Ⅲ, NPP高值、降水较高值状态下, 均可促使3种高值ES(土壤保持、产水服务和固碳服务)与低值粮食供给权衡；同时, 相同驱动因子在不同状态下亦会产生权衡方式的差异。譬如, NPP在较低值状态下促使3种低值ES与高值粮食供给的权衡(情景Ⅳ), 而高值状态下则呈现3种高值ES与低值粮食供给的权衡关系。

4 讨论 4.1 BBN-ES模型

本文基于构建的贝叶斯网络, 探讨多种ES权衡协同及其驱动因素效应差异, 对于认识多种ES间复杂关系, 为后续提出改善生态系统管理和提高人类福祉的对策与措施, 具有重要的理论和实践价值。本文构建的BBN-ES模型总体精度为73.5%, 4种生态系统服务中粮食供给精度最高(75.9%)、产水服务精度最低(70.0%), 表明贝叶斯网络在模拟生态系统服务供给过程准确, 对生态系统服务节点概率预测具有较强的可靠性。模型虽具有较高的精度, 但简单通过自然断点法对连续变量分级有可能导致模型参数的设定具有潜在的不确定性^[37]。因此, 依据研究区的实际情况对变量进行分级, 可进一步提升模型的精度。

4.2 ES权衡协同关系

本文基于BBN-ES模型, 利用联合概率分布和相关系数分析ES间的权衡协同关系。2018年研究区的土壤保持、产水服务与固碳服务之间互为协同关系, 粮食供给分别与土壤保持、产水服务和固碳服务为权衡关系, 该结论与孙艺杰等^[38]、韩磊等^[39]的研究成果一致, 可见贝叶斯网络方法在揭示ES间复杂的非线性关系方面具有一定的优势。需要注意的是, 本研究只是对研究区2018年的ES间的权衡协同关系进行探讨, 并未涉及随时间变化对ES间关系的影响, 因此, 未来有必要采取切实有效的模型与方法, 基于长时间序列分析ES权衡协同关系的时空动态变化趋势。

4.3 ES权衡协同关系的驱动因子

本文在情景分析和相对概率变化的基础上, 确定影响ES权衡协同关系的驱动因子, 并进一步探讨了不同情景下主要驱动因子的效应, 即驱动因子及其在不同状态下对权衡与协同的影响。通过分析发现：不同驱动因子在不同状态下产生相同的协同(或权衡)方式；相同驱动因子在不同状态则会造成协同(或权衡)方式的差异。

研究区自实施退耕还林(草)工程后, 林草地面积的迅速增加提升了植被覆盖度, NPP与降水、植被覆盖度等密切相关, 降水量和植被覆盖度的增加有利于NPP^[40—42]。因此, 为了优化研究区生态系统的整体效益, 应尽可能加强生态系统服务之间的协同关系, 减弱权衡关系。政府部门在进行相关决策工作时, 应重点关注土地利用、降水和NPP等驱动因子, 特别是土地利用的的结构和格局, 如耕地、森林和草地的空间分布。需要说明的是, 不同驱动因子判定标准的设定, 可导致确定主要驱动因子及其效应的差异。本文将概率变化最大值高于35%的关键节点确定为驱动因子, 具有一定的主观性, 未来仍需进一步讨论。因此, 确定适宜于研究区的驱动因子判断标准, 是后续提出针对性对策与措施的前提和基础。

5 结论

本文以陕北黄土高原为研究区, 定量评估其2018年的土壤保持、产水服务、粮食供给和固碳服务, 将生态系统服务与贝叶斯网络结合构建BBN-ES模型；在ES的权衡协同关系与情景分析基础上, 探讨影响多种生态系统服务权衡协同关系的驱动因子。主要结论如下：

(1) 从空间分布特征来看, 陕北黄土高原的土壤保持、产水服务和固碳服务呈“南高北低”的分布特征；粮食供给服务呈现出“中部高、南北低”的分布特征。

(2) 从权衡协同关系来看, 陕北黄土高原土壤保持、产水服务和固碳服务之间互为协同关系, 粮食供给分别与土壤保持、产水服务、固碳服务为权衡关系。其中, 产水服务与固碳服务的协同关系、产水服务与粮食供给的权衡关系较强。

(3) 从情景分析结果来看, 土地利用、降水和NPP是影响ES权衡协同关系的主要驱动因子, 其中土地利用主要影响生态系统服务的协同关系, 而降水和NPP在影响协同关系的同时对权衡关系也有一定地制约作用。研究发现：不同驱动因子在不同状态下产生相同的协同(或权衡)方式, 相同驱动因子在不同状态则会造成协同(或权衡)方式的差异。