基于贝叶斯网络的生态系统服务空间格局优化——以泾河流域为例

余玉洋， 李 晶， 周自翔， 唐承延

（1.陕西师范大学地理科学与旅游学院，陕西 西安 710119；2.河南师范大学旅游学院，河南 新乡 453007；3.西安科技大学测绘科学与技术学院，陕西 西安 710054）

生态系统服务是指人们从生态系统中获得的所有受益，包括直接和间接受益［1-2］。由于生态系统服务与人类的生命和生存密切相关，因此研究生态系统服务对于人类的生存和地球发展都非常重要［3-6］。随着科学技术水平的提升，越来越多的学者对生态系统服务研究集中于供需关系、空间权衡与优化、生态福祉、服务流等方面［7-12］。其中生态系统服务空间格局优化受到越来越多的学者关注，学者们开始尝试运用不同模型与方法进行深入探究，但是由于研究时限不长，所以促使该方面研究处于萌芽时期［13-15］。目前多数空间优化研究是以土地利用空间格局为支撑，其目的是实现一定规模的需求和规划的目标［16］。刘慧灵等［17］通过单位面积上不同土地利用的生态系统服务价值来作为变量系数，建立生态效益函数，进而对土地利用进行空间优化；王观湧等［18］、向芸芸等［19］通过土壤有机碳的最大储量为目标，确定生态效益函数，最终达到优化土地利用的目的。李秀霞等［20］通过设置经济和生态的双重目标，利用系统动力学模型和多目标优化模型，对吉林省的土地利用进行仿真模拟和优化。上述研究大多数是对土地利用进行空间优化，然后在以优化的土地利用为基础，进行生态系统服务的空间制图和分析，不仅忽略影响生态系统服务的其他影响因素，而且没有真正意义上实现生态系统服务的空间优化，即在研究区范围内，如何合理和优化配置自然因素和土地利用等因素，使该区域的生态系统服务能力最强，如何选取不同生态系统服务的影响因素以及量化关键因子对生态系统服务的影响程度，从而达到优化区域生态系统服务的目的，正是该研究所解决的问题。

贝叶斯网络是一种大规模的数据分析工具，具有强大的逻辑分析能力。该网络采用图形化建模的方法，将变量之间的相互作用以网络的形式表达出来，并构建一个系统概率模型［21］，将贝叶斯网络和生态系统服务相结合，可以很好地为生态系统服务的优化管理提供决策支持［22］。Fox等［23］研究了贝叶斯网络在水质和自然资源管理中的价值；Landuyt等［24］利用贝叶斯网络评估池塘综合管理的机会和风险；Dang等［25］将生态系统服务和自然环境因子纳入贝叶斯网络模型，通过不同指标设置不同的研究情景，为农业发展提供概率支持。综上看出贝叶斯网络可以很好地和生态系统服务相结合，从而来实现研究目标。

泾河流域位于黄土高原中部，在六盘山和子午岭之间，对于此区域的研究大多集中在生态和自然环境方面。Chen 等［26］使用主成分分析和相关分析来分析泾河流域的归一化植被指数（Normalized difference vegetation index,NDVI）和降水量，发现植被的动态变化与气候变化、洪水和人类活动密切相关；Xie 等［27］分析了降水-景观-径流变化和相互作用；Guo 等［28］使用滑动偏相关系数法分析降雨和径流关系的变化，并使用双累积曲线法验证其分析；甄霖等［29］通过参与式社区评估法对泾河流域景观管理中存在的相关问题进行深入剖析。从现有研究中看出对于泾河流域生态系统服务的研究还处在初级阶段，并且随着社会经济的快速发展，土地利用格局发生很显著的变化，最终也会导致生态系统服务的格局发生变化，如何通过现有的土地利用格局来优化生态系统服务也显得至关重要，因此，深入研究泾河流域生态系统服务以及空间格局优化会对区域的生态-经济和谐发展奠定基础。本文在生态系统服务空间评价的基础上，引入贝叶斯网络模型，通过筛选关键因子和最优状态子集的方式，来获取4种生态系统服务的空间优化区域，为区域生态健康持续发展奠定基础。

1 数据与方法

1.1 研究区概况

泾河发源于泾源县泾源乡，源区地处六盘山东麓，自西北向东南流经固原市、平凉市、庆阳市、咸阳市等38 个县市，于咸阳市高陵县陈家滩村汇入渭河。泾河流域地理位置介于106°14′～108°42′E、34°46′～37°19′N 之间，地跨陕西省、甘肃省、宁夏回族自治区（图1），流域面积达45421 km2。区域多年平均气温8 ℃，年降水量350～650 mm，降水主要集中在夏季。泾河流域地势较为平坦，耕地面积约占1/3，农业发达；矿藏量较为丰富，经济开发潜力大。泾河流域内黄土高原区域和泾惠渠灌区域，是陕西省主要商品粮油基地。除此之外，还有煤炭、石油、天然气、建筑材料等矿产，如长庆油田。但长期过度开发利用和气候变迁，导致流域内生态环境十分脆弱，其中水资源缺乏已成为研究区社会经济以及生态发展的主要制约因素。

图1 研究区概况Fig.1 Overview of the study area

1.2 数据来源

本文主要数据包括遥感数据（土地利用数据）和非遥感数据（地理基础数据、高程数据、气象数据、土壤类型数据和水文数据等），具体信息见表1所示。

表1 数据来源Tab.1 Data sources

1.3 生态系统服务模型

1.3.1 农业生产力服务模拟首先，从2000—2020年地表覆盖数据中获取耕地资源数据作为计算区域农业生产潜力的基础；其次是基于气象站点数据，通过ANUSPLIN 软件获取2000—2020 平均气温和年降水量数据；最后，采用桑斯维特纪念模型［30］，计算区域2000—2020年农业生产力，公式如下：

式中：CP为以耕地为基础的农作物生产力（t·hm-2）；Pv 为以年平均蒸散量为自变量计算出的植物干物质总量（g·m-2·a-1）；Izrd为土地利用水平，是由中国农用地分级规定决定的［31］。Pv的计算公式如下：

式中：v为年平均实际蒸散量（mm）；p为年降水量（mm）；l为年平均最大蒸散量（mm）；t为年平均气温（℃）。

1.3.2 产水服务模拟水文评价模型（Soil and water assessment tool, SWAT）可以对流域内一系列复杂的物理过程进行模拟，模型中的水量平衡始终是流域内所有过程的驱动力［32］。并且此模型可以基于水量平衡模拟每个水文响应单元的地表径流量和洪峰流量，水量平衡方程如下［33］：

式中：SWt为最终土壤含水量（mm）；t为时间天数；SW0为第i天的初始土壤含水量（mm）；Rday为第i天的降水量（mm）；Qsurf为第i天的地表径流量（mm）；Ea为第i天的蒸散发（mm）；Wseep为第i天通过土壤剖面进入包气带的水量（mm）；Qgy为第i天的回归流量（mm）。

1.3.3 净初级生产力（NPP）服务模拟参考朱文泉等［34］的遥感估算模型（Carnegie-ames-stanford approach,CASA）计算植被NPP，由植物的光合有效辐射（Absorbed photosynthetic active radiation, APAR）和实际光能利用率（ε）2个因子表示［35］。计算公式如下：

式中：NPP(x,t)为像元x在t月内所得到的植被净初级生产力（g C·m-2）；APAR(x,t)为像元x在t月吸收的光合有效辐射（MJ·m-2）；ε(x,t)为像元x在t月的实际光能利用率（g C·MJ-1）。

1.3.4 土壤保持服务模拟采用通用土壤流失方程估算土壤保持量，20世纪50年代由Wischmeier等［36］提出，基于理论以及大量实地观测数据统计分析的基础上，对各因子都有具体解释。通用土壤流失方程的表达如下：

式中：A为土壤保持量（t·hm-2）；R为降雨侵蚀力因子（MJ·mm·hm-2·h-1·a-1）；K为土壤可蚀性因子（t·h·MJ-1·mm-1）；L和S分别为坡长和坡度因子；C为地表植被覆盖及经营管理因子；P为水土保持措施因子。

1.4 基于BBNs的农田生态系统服务空间格局优化

1.4.1 贝叶斯网络模型贝叶斯网络又称信度网络，是Bayes方法的扩展，目前不确定知识表达和推理领域最有效的理论模型之一［37］。1988 年Pearl 提出后，已经成为近几年来研究的热点。一个贝叶斯网络是一个有向无环图（DAG），由代表变量节点及连接这些节点有向边构成［38］。

其中条件概率表（CPT）表达了父节点X与子节点Y之间的关系强度，表中每一行由父节点的状态组合和条件概率P(y|x)组成，即父节点X在离散状态x下，子节点对应状态y发生的概率。子节点的概率分布由父节点的概率分布及条件概率表决定，对于没有父节点的节点，其概率分布为先验概率分布P(X)。所有节点的条件概率分布相乘得到贝叶斯网络中所有变量的联合概率［式（8）］。贝叶斯网络能够有效计算网络中任意节点概率分布的基础，为本文优化生态系统服务格局提供了理论［39-40］。

式中：P（X1，X2，…，Xn）为一个联合离散概率分布，这种分布的随机变量具有1，2，…，n个离散值；n为随机变量的最大值；（X1，X2，…，Xn）为随机变量。

以4 种生态系统服务量化模型为依据，选择与其相关的影响因子作为节点，构建生态系统服务的贝叶斯概念网络（图2）。借助ArcGIS 软件，离散节点对应的栅格图层，根据变量的实际分布情况，参考学者对贝叶斯网络栅格分层的方法［9,22］，运用Arc-GIS软件的自然间断法将上述每个变量栅格数据分为4 个等级：最高、高、中、低，具体离散化分级标准如表2所示。提取离散数据值到属性表并将概率条件输入贝叶斯网络模型，便可得到泾河流域生态系统服务的贝叶斯网络。其他节点的条件概率表可由条件概率公式［式（9）］计算获得。

表2 因子状态分级Tab.2 State classification of factors

图2 生态系统服务贝叶斯网络构建原理Fig.2 Schematic diagram of Bayesian belief network construction for farmland ecosystem service

式中：P(AB)为事件A 与B 的联合概率，即2 个事件共同发生的概率；P(A)为A的先验概率，即事件A发生的概率；同理，P(B)为事件B发生的概率。

1.4.2 农田生态系统服务格局优化构建出泾河流域生态系统服务贝叶斯网络模型，计算出条件概率表和敏感性因子，确定出关键变量关键状态子集。以ArcGIS 软件为工具，可视化子集分布空间，最终绘制不同等级所对应的关键变量关键状态子集同时发生的区域。筛选关键变量关键状态子集，包括以下步骤：

（1）确定各变量的关键状态

根据节点各状态的概率及两两节点状态组合的联合概率，计算每个影响因子与各服务间的条件概率，选取不同等级对应各个影响因子的条件概率最大状态作为变量的关键状态。当服务处在最高、高、中、低4 种不同级别时，将对应4 种变量状态集合。

（2）确定关键变量

为了评价贝叶斯网络节点的相对重要性，Netica 提供了敏感性分析［23,41］，主要是用于评估生态系统服务节点是否会敏感地感知其他影响因子节点的变化［9,22］。据此确定生态系统服务关键变量关键状态子集，然后选取变量关键状态子集的面积与研究区总面积比值最大的情景作为生态系统服务最优的情景，将该情景下的关键变量关键状态子集选定为最优状态子集，然后结合县区的发展和管理等划分出生态系统服务的优化区。

2 结果与分析

2.1 SWAT模型率定与验证

SWAT 模型有26 个与径流模拟相关的敏感性因子，不同因子影响各异。本文根据SWAT-CUP 软件进行模拟率定，通过T检验确定水文参数的敏感性相对显著性，P值确定其显著性特征［42-43］。进而选取相关参数（表3），通过调整参数，误差能够控制在一定范围之内。

表3 参数敏感性分析Tab.3 Sensitivity analysis result of parameters

由于张家山水文站观测数据限制，本文以2000—2010年实测数据和模拟数据进行率定，根据不断调整得到的参数，获得模拟效果最佳的一组数据。如图3a所示，其中决定性系数（R2）为0.74，纳什效率系数（NSE）为0.73，都大于0.5，表示模拟值与实测值具有较高的相关性，模型的模拟结果是可靠的。模拟2011—2020年泾河流域的月径流量，与张家山水文站数据进行比对进行验证，可以得出R2为0.78（图3b），径流实测与模拟值具有较高的相关性，模拟结果是可信的，据此模拟2000—2020年泾河流域产水服务。

图3 率定期和验证期逐月径流量模拟值与实测值Fig.3 Simulated and measured values of monthly runoff during the regular period and the verification period

2.2 生态系统服务时空格局分析

生态系统服务功能之间相互作用与影响，对流域生态格局发展与优化具有显著影响。本研究通过模型模拟了2000—2020 年泾河流域产水量、NPP、农业生产力和土壤保持4 种生态系统服务量的空间分布（图4）。区域产水量在时间上呈增长趋势，2000 年平均产水量为489.6 m3，而2020 年平均产水量为729.6 m3。在空间上呈现北部和西南部高，中间低的特征，高值区域主要集中在环县、千阳县和凤翔县等区域。此外，低值区主要集中在西部的泾源县以及北部的盐池县。但是2020 年产水量高值区域较为分散，主要呈现中间低，四周高的分布格局，主要是因为2020年相较于其他年份降水重心有所转移，并且受季风天气的影响，导致区域降水量年际差异明显，所以在一定程度上影响产水量的空间分布。NPP 在年际变化中呈现增长的趋势，2000 年NPP 的平均值为185.1 g C·m-2，总量为7.7×106t，而2020年NPP的平均值为461.9 g C·m-2，总量为2.4×107t，总量相较于2000 年大约增加了3 倍。NPP 在空间上呈现西北向东南逐渐递减的趋势，高值区域集中在流域两侧子午岭山系及西部六盘山地区。低值区主要原因是植被以典型草原为主，并且逐渐东南方向退化，沙漠化问题较严重。农业生产力在时间上从由2000 年的280.9×107t，增长到至2020 年的368.3×107t；在空间上自西向东呈增加趋势，南北向呈“马鞍形”分布，空间分布差异较小。高值分布面积有所增加，集中分布在陕西省长武县、淳化县以及宁夏回族自治区泾源县。土壤保持从时间尺度上看，总体呈波动变化趋势，但年际差异较大，土壤保持从2000年的3.0×108t增加到2020年的5.3×108t，土壤保持呈现显著增长趋势。从空间分布特征看出，土壤保持强度空间分布存在一定差异，整体呈现出北部、中部低，西南部和东南部高的分布格局。土壤保持强度低值区主要分布在流域北部和中部，北部为地形起伏变化的黄土丘陵区，加之年降水量低于区域平均值，植被生长水热条件很差，导致该地区植被覆盖度较低。高值区主要分布在流域西南部和东南部，部分区域大面积为林地，海拔较高，并且规划了森林保护区，人类活动对地表的干扰小。

图4 泾河流域生态系统服务量空间变化与分布Fig.4 Change and distribution of ecosystem services in Jing River Basin

2.3 生态系统服务空间格局优化

贝叶斯网络模型以概率相关理论为基础，得到研究区不同生态系统服务关键变量最优状态子集，进而筛选出需要优化的区域。对于泾河流域生态系统来说，人为因素和自然环境的双重作用，对生态系统服务的影响不容小觑。不同的生态系统服务，筛选不同的关键变量，得到不同生态系统服务的条件概率表，如表4～7所示。根据表4所示，土地利用水平和蒸散发是影响农业生产的关键变量，当农业生产力出现最高（0.353）、高（0.422）、中等（0.620）和低水平（0.491）概率时，所对应的土地利用水平和蒸散发量分别处于不同程度的水平。对比得到土地利用水平处在高级水平即农业用地水平，蒸散发量处于高值区间5362～9827 mm 时，农业生产力达到最大可能概率（0.620），处于中值范围（57.7～60.0 t·hm-2），对其进行空间可视化展示（图5），优化区域主要分布在环县、华池县、彭阳县、庆阳县以及镇原县等北部黄土区，对这些地区加大农业生产的投资和管理力度，能够较为迅速地提高泾河流域整体的农业生产力。当土地利用类型为草地，土地利用水平处在中级水平，降水量处于最高值584～693 mm 时，气温小于8 ℃时，产水量达到最大可能概率（0.805），处于高值范围（408.61～1001.75 m3），优化区域主要集中在环县的中南部以及彭阳县和镇原县附近，对区域水资源进行高效利用，能在一定程度上增加地表径流量，最终提升区域的产水服务能力。当土地利用类型为未利用地，太阳辐射、气温和植被覆盖度都为最高水平时，NPP服务达到最大可能概率（0.878），处于中值范围（11.97～135.89 g C·m-2），优化区域主要集中在彭阳县和华池县、镇原县、环县中部和南部，吴起县南部和合水县北部，提升区域的植被覆盖率对提高NPP服务有深远的影响。通过关键因子的筛选，影响土壤保持的关键因素是降水量、坡度、水土保持措施因素、植被和土壤类型，降水量概率是最重要的因素。当土壤类型、降水量、坡度以及水土保持措施因素位于中级水平，植被覆盖度为最高水平时，所对应的土壤保持出现最高水平的概率最高，达到0.824，该子集主要分布在彭阳县和华池县、庆城县北部和镇原县以及环县南部，对该区域增加植被覆盖，抚育和保护森林资源，加强生态环境的建设和投入，有效控制水土流失带来的损失，会对区域的水土保持服务起到正向的作用。

表4 农业生产力节点条件概率Tab.4 Conditional probability of the node crop production

表5 产水节点条件概率Tab.5 Conditional probability of the node water yield

表6 NPP节点条件概率Tab.6 Conditional probability of the node NPP

表7 土壤保持节点条件概率Tab.7 Conditional probability of the node soil conservation

最后对4种生态系统服务的优化区域进行空间叠加分析，得到综合的优化像元（图5）。优化区域主要集中在彭阳县的中部和西南部以及环县的零星区域，部分区域应该采取相应的生态环境保护措施，加快生态、民生和多功能林业的发展，对于草原生态要坚持保护与建设并举，保护优先的原则，设置一定的惩奖措施，在一定程度上提升居民对于生态环境保护意识，促进区域多种生态系统服务良性发展。

图5 综合生态系统服务优化区Fig.5 Integrated ecosystem service optimization zone

3 讨论

3.1 生态系统服务空间差异性分析

生态系统过程是自然和人类活动共同作用的结果，如何构建合理的模型来模拟区域生态系统服务是深刻了解生态过程的关键，本文通过不同的生态系统服务评估模型得到4种生态系统服务的空间格局。为了探究研究结果的准确性，将本文的研究结果和其他学者进行对比。刘宇等［44］通过用Penman-Monteith 公式计算泾河流域近70 a 的径流变化，以及黄晨璐等［45］在统计数据的基础上，运用相关的模型来模拟泾河流域径流状况，在时空尺度上和本研究结果具有相似性；在NPP和土壤保持服务的时空分布方面与卞鸿雁等［46］以及Zheng等［47］的结果具有相似性，都在年际间趋于稳定上升的趋势，并且时空差异显著。此外，根据本文研究结果，4种生态系统服务的优化区域都集中在黄土高原和泾河流域北部的山区，制约区域发展的关键因子主要是土壤、植被和降水，而关键因子对于评估区域4种生态系统服务有着至关重要的影响。泾河流域产水服务的空间变化特征和土地利用、蒸散发和植被水分利用系数相关，2000—2020年区域耕地面积减少，森林和草地面积增加，进而增加区域植被水分利用系数，植被类型从禾本科植物向森林植被转化所引起。2020 年相较于2000 年，降水量增加，而部分区域径流量减少，暗示区域内潜在蒸散发则逐渐增大，表明区域径流量对于降水和潜在蒸散发较为敏感［44］。产水和降水的空间分布较为一致，主要取决于气候和地形因素相关，西北部属于黄土高原区，受副热带高气压控制，干旱少雨，产水量也较少。泾河流域大多区域经济相对落后，居民的生活需求对区域生态系统服务依赖性很强，农业生产力较强的区域，居民能够满足生活，少数区域粮食供给能力不足［48］。北部区域以草地为主，牧业长期集中与此，长期的粗放经营，会造成土地严重的沙化，进而影响区域的农业生产力以及NPP服务，造成两种生态系统服务在区域上呈现南北差异化明显。而随着退耕还林、封山育林、飞播、人工植树等生态措施的实施，促使区域内林地面积增加，提升区域的NPP 和水土保持能力。本研究建议对生态系统服务优化区实施提升肉类的产量以及薪柴供给的手段等生态措施和生态补偿政策［49-50］，不仅可以提升居民的收入和生活质量，而且还可以保障原有的土壤保持和NPP等生态系统服务的增强，促进区域生态经济和谐发展。为了深入研究生态系统服务的区域差异性，不仅要定性分析，还要利用模型进行定量的研究，未来应该通过多视角、多指标对生态系统服务驱动因素进行定量分析，促使研究结果更加精准。

3.2 生态系统服务空间格局优化

贝叶斯网络可以将多元知识进行图解可视化的概率模型，它不仅包括各个网络节点之间的因果相关性，而且结构还比较灵活，既可以优化现在的空间格局，还可以预测未来的空间格局分布。李婷等［9］通过贝叶斯网络和固碳服务相结合，来探究固碳服务的格局优化问题；曾莉等［22］通过不同发展情景水源涵养的状态概率分布，最终给出水源涵养服务的空间格局优化策略。但是贝叶斯网络是具有不确定性的，主要包括结构不确定性、输入数据不确定性和参数不确定性［51］。本文将产水量、NPP、土壤保持和农业生产4种生态系统服务和贝叶斯网络相结合，探讨了不同生态系统服务以及综合优化区域，为优化区的保护和生态经济和谐发展奠定基础。目前，对生态系统服务空间格局优化的研究主要集中在围绕土地利用模式的变化设计方案，预测各种方案下的服务变化，并提供优化建议［16-17］。然而，土地利用并不是唯一影响生态系统服务的因素，而各因素对不同生态系统服务的影响程度也是不同的。本文在贝叶斯网络的基础上，提出了关键子集法，对泾河流域生态系统服务的空间格局进行了优化，研究结果为研究区生态环境保护区划提供决策支持，丰富生态系统服务空间格局优化的研究方法。但是自然和人类活动共同影响着区域的生态环境，本文仅考虑降水量、气温、土地利用等因素对生态系统的影响是不够的，未来应该结合更多的影响指标，比如人类活动、生态政策以及社会经济等因素来构建多情景的模式，将情景分析和贝叶斯网络相结合，更加全面的为区域空间格局优化提供依据。

4 结论

（1）SWAT 模型在泾河流域有较好的适用性。校准后的SWAT模型对径流量的模拟结果与实测值相比，具有较高的决定性系数（R2＞0.6）、纳什效率系数（NSE＞0.5），表明该模型能较为准确地模拟流域生态水文过程，为进一步评估产水服务提供保障。

（2）2000—2020年泾河流域4种生态系统服务的时空差异性较为显著。在时间尺度上，随着退耕还林还草以及居民生态环保意识的提升，促使4 种生态系统服务都呈现波动中上升的趋势；而在空间尺度上，NPP、农业生产力以及土壤保持服务变化较小，呈现较为稳定的趋势，而产水量和区域气候有着直接的关系，导致产水服务在区域内呈现显著变化。

（3）贝叶斯网络模型是一种强大的神经网络工具，可以很好地与生态系统服务相结合，在生态系统服务空间优化过程中发挥重要作用。本研究通过对4 种生态系统服务优化区域进行叠加分析，发现综合优化区域集中在彭阳县的中部和西南部以及环县的零星区域，对区域应该采取相应的生态环境保护措施，要坚持保护与建设并举，保护优先的原则，进而促进区域多种生态系统服务良性发展。