基于多场景模型的沙漠-绿洲交错带林草生态网络模拟

苏 凯 于 强 YANG Di 张启斌 杨 斓 孙小婷

(1.北京林业大学精准林业北京市重点实验室， 北京 100083； 2.佛罗里达大学地理系， 盖恩斯维尔 FL 32611)

0 引言

生态系统是人类生存和发展的物质基础，也是人类的生命支持系统[1]。人类的生存和发展取决于生态系统的稳定性和健康。在干旱地区，沙漠-绿洲生态交错带作为缺水沙漠和富水农田系统之间的过渡区，在确保区域生态安全和维护区域内部稳定方面发挥着关键作用。然而，沙漠-绿洲生态交错带也是脆弱的，其状态极易受发展策略的影响[2]。当发展策略以经济优先时，景观实现了经济的发展，但导致了地下水位下降，森林遭受破坏，荒漠化加速，湿地退化等生态问题日益严重[2-3]，使区域生态系统不平衡且难以恢复[4]。如何有效管理荒漠化和遏制沙漠扩张一直是一个世界性的问题。

生态网络通过带状植被、河流等形成的线性廊道或者是潜在线性生态廊道将荒漠绿洲交错区相互孤立的生态斑块连接，形成“点、线、面”相结合的空间生态安全格局，提升景观的自我调节能力，维持区域生态环境稳定，抵御沙漠化风险[5]。从生态网络的视角来看，土地沙漠化及其造成的一系列影响的过程实际是生态网络中生态节点与生态廊道不断被破坏与消失的过程[6]。在这一过程中，网络结构的完整性不断降低，生态功能逐步下降，网络中生态、物质与能量的交换逐渐受阻。

生态网络的构建已经形成基本的研究范式：生态源地的识别与潜在生态廊道的提取。识别生态源地，主要是通过评估生态适宜性、生态重要性或生态连通性。其中基于生态系统服务的生态重要性评估最为常见[7]。提取生态廊道的方法通常基于土地覆盖的价值分配构建阻力面，阻力面描述了物种穿越不同生境斑块的难度，反映了物种对生态过程的水平阻力。阻力面表征了景观异质性对生态过程流量的影响。最小成本分析通常用于提取生态走廊。但这种方法忽略了生态源地之间的引力，没有考虑生态源地之间的相互作用。引力模型早期被应用于城市间的空间结构和零售市场研究[8]，后来也被广泛用于区域经济联系[9]、城市间相互作用结构，以及城市间贸易、物流联系等研究[10]。目前引力模型广泛应用于度量区域经济联系量或空间相互作用量。

然而，沙漠-绿洲生态交错带的森林-草地生态网络(FG生态网络)是一个动态网络，而且是一个脆弱的系统，受自然条件和经济活动的影响。因此，本文以中国西北地区典型的沙漠-绿洲交错带为研究对象，提出林草生态网络(FG生态网络)构建的综合模型，该模型包括两个子模型，利用改进的生态阻力模型获取生态流流动的最小累计阻力面，利用改进的生态引力模型计算生态源地之间的引力。通过逐一比较各生态源地之间的合力，利用最小成本分析提取生态廊道，构建连接生态源地的生态网络。同时，制定多种发展战略，探讨发展战略中生态保护与经济发展的比例对生态网络演变的影响。采用复杂网络分析方法，对研究区林草生态网络进行拓扑分析和空间结构分析。旨在为干旱区经济发展和生态环境保护研究提供新的视角，为干旱区可持续发展提供理论支持。

1 材料与方法

1.1 研究区概况

研究区域(106°7′～107°23′E，39°56′～40°56′N)位于乌兰布和沙漠东北缘，面积5 923.9 km2，海拔968～1 338 m。研究区东部是中国的母亲河——黄河；东北部的平原是中国重要的粮食产区——河套平原；西侧是狼山；西南部是乌兰布和沙漠的腹地。根据1983—2016年磴口沙漠生态站的气象资料，该地区年平均降水量为144.5 mm，年平均蒸发量为2 398 mm。研究区域的位置如图1所示。

图1 研究区位置Fig.1 Location of research area

1.2 数据来源与处理

使用2018年夏季的Sentinel 2遥感影像，分辨率为20 m。数据集由欧洲航天局提供(http:∥www.esa.int/ESA)。通过辐射定标、大气校正、图像增强和几何校正等进行图像的预处理。利用最大似然监督分类方法对遥感图像进行分类，并提取研究区域内的土地利用类型信息。从土地利用数据中提取水网和道路网，通过多缓冲分析获取距离因子数据。从土地利用数据中提取住宅区、道路网和水网，通过密度分析获取密度因子数据。从处理后的图像中提取归一化差异植被指数(NDVI)、修正归一化差异水体指数(MNDWI)和温度植被干旱指数(TVDI)。DEM由中国科学院计算机网络信息中心地理空间数据云平台(http:∥www.gscloud.cn)提供。地下水埋深的空间分布数据来源文献[11]。

1.3 方法

1.3.1生态网络多场景仿真模型

FG生态网络多场景仿真模型由两个子模型组成。其中利用生态阻力模型模拟不同发展策略中生态源地之间的最小累计耗费生态阻力，利用生态引力模型模拟生态源地之间的相互作用。通过比较不同发展方案时生态源的综合受力情况，分析FG生态网络的变化。

在干旱地区，尽管水资源稀缺，但仍然可以自然形成相对稳定的生态系统[12]。一些研究认为，生态资源地之间存在相互依赖的关系，通过生态网络，在干旱地区建立了较为稳定的生态格局，维护了区域生态安全。FG生态网络由生态源和连接生态源的生态廊道组成[13]。在生态网络中，生态源是提供各种物质和能量的空间单元或生态系统，生态廊道是能量和物质流动的重要通道。在FG生态网络中，生态源地之间的相互作用会产生生态引力，但同时也会受到来自基质的阻力，生态源地之间互作用的示意图如图2所示，图中F12=F21，为生态引力，fmin表示生态源地之间的最小累积生态阻力。

图2 生态源地间作用力示意图Fig.2 Schematic of interaction forces between ecological sources

生态网络多场景仿真模型的详细算法步骤如下：

(1)计算生态源地之间的最大生态引力和最小累积阻力。当生态引力大于生态阻力时，源地之间将形成生态廊道，生态流将通过廊道在源地之间传递。当生态引力小于生态阻力时，源地之间就没有生态廊道。当生态引力等于生态阻力时，源地之间的生态廊道将非常脆弱，基质的微小变化可能导致生态廊道断裂。

(2)重复计算，直至遍历完成整个FG生态网络的源地，最终获得FG生态网络空间分布网格模拟结果。然后模拟下一个开发计划中FG生态网络的空间分布，直到所有发展策略模拟结束。

由于政策的实施是统一的，因此往往不适合当地条件。本文通过调整发展策略中经济发展和生态保护的比例，设定了11个方案。其中2个是极端发展计划：没有经济发展的纯生态保护模式和没有生态保护的纯经济发展模式。实际上，研究区的发展既是生态保护，也是经济发展。但是，该地区的生态环境极其脆弱，是控制荒漠化的重要区域。我国政府颁布了更为严格的生态保护政策，严格限制了该地区的经济发展，重点是改善和保护生态。因此，研究区目前的发展模式可以看作是一种纯粹的生态保护模式。不同情景的发展模式将导致生态环境条件的差异。每个发展计划的发展策略定义为1，生态保护的比例用m表示，而n代表经济发展的比例，则m+n=1。当m=1.0时，n=0，(1.0，0)代表纯生态保护模式，即以生态保护为重点的研究区发展策略。当m=0，n=1.0时，(0,1.0)代表纯经济发展模型，即研究区域优先考虑经济发展。最后，获得了11个开发方案：(0,1.0)、(0.1,0.9)、(0.2,0.8)、(0.3,0.7)、(0.4,0.6)、(0.5,0.5)、(0.6,0.4)、(0.7,0.3)、(0.8,0.2)、(0.9,0.1)、(1.0,0)。

1.3.2生态阻力模型

生态资源具有向外扩展的能力[14]。在扩张过程中，由于土地利用和景观的差异，会遇到生态障碍，受到生态阻力的影响。荷兰生态学家KNAAPEN提出了最小累积阻力(MCR)[15]，本研究根据研究区性质对MCR模型进行改进。本文认为不同源地具有不同的生态能量，它们辐射和传递能量的能力不仅与其面积有关，而且与源的类型和形状有关。目前，一些学者对此进行了研究和讨论[16]。生态源的大小将对物质交换或能量流产生强烈影响，其边缘也会影响物质、能量和未来的物种分布。因此本文增加了形状指数(Li)和NDVI、MNDWI来描述生态资源的特征。由于研究区位于干旱地区，蓝源(河流、湖泊)水平高于绿源(森林、草地)水平，本文提出了生态源地质量的概念，并改进了MCR模型。生态源地质量的定义是任何生态源地都具有向外辐射，传递能量或接收能量的特性，其生态源地质量与生态源地的面积、形状、属性有关。

(1)

式中Ei——生态源地i的周长

Ai——生态源地i的面积

(2)

式中mi——生态源地i的生态源地质量

本研究将基质阻力分为两类：①建设用地集中区和农田等生态屏障。②地形坡度、植被覆盖度、土地利用类型、地下水埋深、土壤干旱程度等生态障碍。一般来说，当地表变得陡峭时，植被覆盖度和土壤干旱指数会降低，地下水埋藏得更深，阻力系数会越大，这将阻止生态资源的扩大。建设用地集中区和耕地是人类活动的主要集中区，土地利用方式难以改变。通过对研究区域的分析，改进的MCR模型更适合于研究区域的生态阻力模拟，改进的公式为

(3)

式中VMCR-m——生态源地扩展的最小累积生态阻力

Dij——生态源地i和j的空间距离

Ri——生态源地i运动过程中的阻力系数

1.3.3生态引力模型

生态引力模型是研究两个物体或空间相互作用的模型，源于牛顿万有引力定律[17]。19世纪，学者研究城市系统的相互作用时将引力模型引入到地理空间的研究中。后来，地理学家对引力模型进行了广泛的研究，并用于分析城市之间贸易的相互作用，以及国家之间贸易交流形式的研究。2010年，KEUM[18]证明了引力模型在国际贸易和旅游业的可用性。

本文尝试将引力模型引入景观生态领域，计算生态源地之间的生态引力。生态引力与最小累积耗费距离的平方成反比，并且与生态源地质量成正比。生态引力模型为

(4)

式中Fij——生态源地i和生态源地j之间的生态引力

αi、αj——生态源地i、j的类型

G——常数，通常为1

mj——生态源地j的生态源地质量

1.3.4生态网络性能分析

本文选择复杂网络分析中衡量网络拓扑结构和节点重要性的6个图论指标揭示各发展策略下FG生态网络的连接特征，为FG生态网络的优化提供重要参考[19]。从网络整体性与节点的重要性角度对生态网络的性能进行分析，其中连通性c(G)[20]、核数k(G)[21]、网络密度d(G)[22]、平均节点连接数a(G)4个指标评价网络的整体特征，而介数、PageRank值则评价节点的的重要性程度[23]。

表1 生态网络性能分析指标Tab.1 Indicators for ecological network performance analysis

1.3.5骨架廊道识别

克鲁斯卡尔(Kruskal)算法是骨架廊道识别的常见算法[24]，其主要思想是将初始最小生成树边为0，在迭代中，每次选择最小代价边(即本文中的生态廊道)加入最小生成树集合，最终提取出骨架廊道[25]。记原网络Graph中有v个节点，e条边，将Graph中的顶点提取，形成新的网络Graph(new)，此时该网络中没有边，将原图中的e条边按照边权从小到大排序，以边权最小的边开始，遍历原网络中的所有边，直到v个节点都在同一个连通分量中。循环的条件为：若某条边两侧的节点在Graph(new)中不在同一个连通分量中，将该边增加到Graph(new)中。Kruskal算法的图例描述如图3所示，图3a为原始的加权网络，图3b～3i为骨架廊道的提取过程。

图3 Kruskal算法图例描述Fig.3 Legend description of Kruskal algorithm

2 结果与分析

2.1 生态源地识别与提取

生态源地提取是生态网络提取的第一步。根据景观生态学中的“源-汇”理论:生态源地是多种生态服务的源头，具备景观连续性与完整性，对于维持生态系统稳定具有重要作用。选取对乌兰布和沙漠东北缘具有重要生态意义的林地、草地、水体3种景观类型斑块，通过各斑块面积、NDVI平均值、MNDWI平均值及斑块形状指数筛选生态源地。利用ArcMap软件对各斑块面积进行统计，利用Zonal Statics工具提取NDVI平均值、MNDWI平均值，采用Fragstats软件计算斑块形状指数。通过熵值法确定斑块面积、斑块NDVI平均值、MNDWI平均值及斑块形状指数的权重，综合评价各生态景观斑块的生态功能重要性程度，根据各斑块生态功能重要性程度的统计分布特征，选取重要性排序中从高到低占据重要性总和前60%的生态用地斑块作为生态源地。最终获得研究区生态源地共481个，细碎及孤立的斑块在筛选过程中被大量去除，景观单元中的核心斑块得到了保留，提取结果如图4所示。

图4 研究区生态源地分布Fig.4 Distribution of ecological sources in study area

其中水体型源地共37个，面积占所有生态用地斑块的25.47%，水体型源地主要分布在研究区东部，西部区域也有零星水体分布。黄河是研究区最大面积的水域，黄河、奈伦湖等面积较大的水域均由黄河水引流形成。林草地型源地与水体型源地相比更加破碎，生态源地总数为439个，占全部生态源地的91.26%，而面积仅占全部源地的74.52%，林草地分散分布在研究区中部农田集中区、沙漠交错地带、狼山南麓以及黄河东岸，其他区域也有少量分布。

2.2 生态阻力模型构建

景观生态流从源向外流动需要克服阻力，受多种因素的影响。研究区的生态环境特征由12个因素评估，包括地形因子(海拔、坡度)、植被覆盖因子、水文因子(MNDWI、TVDI、地下水埋深)、景观因子和密度因子(居民点密度、路网密度、水网密度)等。每个因子的阻力根据评估系统设定，如表2所示。

此外，发展策略将影响某些因素的阻力，并最终影响研究区域中基质的最小累积消耗阻力。例如，植被覆盖度越高，水网越近，生态保护模式中的生态阻力越小，而经济发展模式则相反。纯生态保护模式((1.0,0)模式)和纯经济发展模式((0,1.0)模式)的生态阻力如图5所示。在(1.0,0)模式下，基质的生态阻力为17～32，从东北向西南逐渐增大。阻力较低的地区主要位于东北部，由于黄河附近水资源丰富；而阻力较大的地区位于乌兰布和沙漠、狼山和黄河的东岸，这里是库布其沙漠的西缘。在这种发展模式下，严格限制经济发展，优先考虑生态保护的发展策略，使得水资源和林草地得到较好的保护。研究区的生态源地及其周边区域的生态阻力都较低。在(0,1.0)模式中，水资源主要用于农业生产和工业制造，而很少用于生态建设和保护，这使得在东北地区引起较高的基质阻力。在中部的乌兰布和沙漠地区以及狼山的前冲积扇中，基质阻力相对较低。以上情况表明，在现状的基础上大力发展经济将严重阻碍生态流在源地之间的流动，影响生态的发展。

最小累积生态阻力是生态流从一个生态源地到另一个生态源地，克服生态阻力面的最小累积成本。在改进的MCR模型的基础上，采用成本距离模型计算不同生态源地之间的累积生态阻力。(1.0,0)模式和(0,1.0)模式的最小累积生态阻力如图6所示。由图6a可知，在(0,1.0)模式中，累计生态阻力最大值是380 727，最小值是28 459。累计阻力较大的区域主要位于南部的乌兰布和沙漠。一方面是因为该区域生态环境条件较为恶劣，地下水埋深、水体密度、NDVI等多因子的阻力均较高。另一方面，该区域生态景观斑块较少，且与现存的生态用地距离较远，生态阻力的累积作用较为明显。研究区北部同样出现了累积阻力较大的区域，尽管该位置生态条件较好，但这个区域是耕地，为生态屏障，阻力较大。且在该发展模式中，农业生产中大量从地下水、地表水中使用水资源，导致生态源地的供水水量较少，这也导致了在(0,1.0)发展模式中基质有比较大的生态累积阻力。这表明，在研究区内大力发展经济将严重阻碍生态流在源地之间的流动，严重影响生态源的发展。由图6b可知，在(1.0,0)发展模式中，累积生态阻力最大值是346 774，最小值是17 467。由于生态条件恶劣，累积阻力较大的区域同样主要位于南部的乌兰布和沙漠。不同于(0，1.0)模式，在该模式中生态阻力较小的区域显著增加。

表2 生态阻力评价指标体系Tab.2 Ecological resistance evaluation index system

图5 (0,1.0)模式和(1.0,0)模式的生态阻力Fig.5 Ecological resistance in (0, 1.0) mode and (1.0, 0) mode

图6 (0,1.0)模式和(1.0,0)模式的最小累积生态阻力Fig.6 Minimum cumulative ecological resistance in (0, 1.0) mode and (1.0, 0) mode

2.3 生态引力模型

通过将生态源地的质量和其间的最小累积耗费距离代入生态引力模型，可以计算出任何生态源地之间的生态引力。虽然生态源地会受到基质的阻碍作用，但它们之间仍然可能通过生态廊道形成联系。生态引力越大，生态源地之间的联系越紧密，它们之间的生态流传递能量、物质就越频繁。而且当生态引力足够大时，生态源地会通过不断吸收和合并外部资源的方式，扩大其自身规模。但如果生态引力很小，在生态源地之间形成生态廊道的可能性不大，除非基质发生变化，否则它们之间将不会形成联系。

以1号生态源地为例，表3显示了1号生态源地与其他生态源地之间的生态引力。由于生态引力在1号生态源地和342号生态源地或233号生态源地之间引力最大，因此在它们之间形成生态廊道的可能性最高。相比之下，由于在1号生态源地和35号生态源地或66号生态源地之间生态引力最小，形成生态廊道最困难。

表3 生态源地之间的生态引力(源地1)Tab.3 Ecological gravitation between ecological sources (source 1)

2.4 生态网络多场景模拟

图7 (0.1,0.9)模式到(0.9,0.1)模式的最小累积生态阻力Fig.7 Minimum cumulative ecological resistance from (0.1, 0.9) mode to (0.9, 0.1) mode

生态引力由生态源地的性质和其间的引力半径决定，并且不随着基质的变化而波动。相比之下，源地之间的最小累积生态阻力与基质密切相关。通过使用栅格计算器模块获得本研究中其他9个发展策略的累积生态阻力。如图7所示，9种模式的累积生态阻力范围为27 743～379 056、27 064～377 915、26 132～375 922、25 875～373 169、23 509～270 376、22 205～368 743、21 413～363 205、19 976～359 492和18 603～352 301。累积阻力大的区域位于南部的乌兰布和沙漠以及北部河套平原附近的农业开发区。随着经济发展的重要性增加，生态保护的重要性逐渐降低，生态积累阻力较小的灰色地带逐渐减少，这将使生态资源地之间能源和物质交换的可能性变小。

通过模拟不同开发方案中基质的最小累积阻力，比较生态源地之间的生态阻力和引力，确定其是否能形成生态廊道，然后利用Python脚本语言模拟生态引力与生态源阻力之间的相互作用，构建FG生态网络。模拟结果如图8所示。

图8中FG生态网络分布不均匀。例如，北方生态源地之间的生态廊道数量较多，生态网络较为复杂，而南方的FG生态网络则具有相反的特征。(1.0,0)模式的FG生态网络比(0,1.0)模式更复杂，(1.0,0)模式下生态源地之间的联系更紧密。在(1.0,0)模式向(0,1.0)模式的转变中，发展策略中经济发展的比例增加，一些地区的生态阻力逐渐增加。生态源地之间的最小累积生态阻力大于生态引力，这导致生态廊道破裂并导致FG生态网络的破坏。

2.5 生态网络性能分析

采用Matlab计算不同情景下FG生态网络的连通性，核心数和密度、平均节点连接数以及介数和PageRank值。计算结果如表4所示。(1.0,0)模式是目前研究区域的发展模式，其网络连通性最大，为329，占网络生态资源总量的68.24%。结果表明，去除329个生态源地可以将(1.0,0)模型的FG生态网络转化为平凡图，从而完全丧失网络的生态功能。

图8 (0,1.0)模式和(1.0,0)模式中的FG生态网络Fig.8 Forest-grass ecological network in (0, 1.0) mode and (1.0, 0) mode

表4 11个模拟场景下FG生态网络的结构特征Tab.4 Structure characteristics of forest-grass ecological network in 11 simulated scenarios

从网络的整体特征来看，随着发展策略中经济发展的比重增加、生态保护的比重减小，FG生态网络的连通性逐渐降低，并且FG生态网络的稳定性随着连通性的降低而恶化。在纯粹经济发展的(0,1.0)模式下，FG生态网络的连通性最低。仅需要去除265个生态资源可以将(0,1.0)模式的FG生态网络转化为平凡图。核数的变化趋势与连通性的相同。随着经济发展的比重逐渐增大，FG生态网络的核数减少，同时网络的连通性和稳健性也在逐渐减弱。网络平均节点连接数的变化趋势与核数和连通性相同。就网络密度而言，网络密度在从(1.0,0)到(0,1.0)模式的演变过程中先增加后减小。FG生态网络的密度在(0.9,0.1)模式中最大。从节点的重要性程度来看，节点的最大介数在(0.9,0.1)模式中最大，随着发展策略中经济发展的比重增加而持续减小。随着发展策略中经济发展比重的增加、生态保护的比重减小，节点的最大PageRank值逐渐降低。

2.6 关键节点及骨架廊道识别

提取研究区生态节点的介数、PageRank值，统计分布中前10%的生态节点作为关键生态节点，利用最小生成树算法，提取研究区潜在生态廊道中的骨架廊道，结果如图9所示。

图9 研究区关键生态节点及骨架廊道Fig.9 Key ecological nodes and skeleton corridors of study area

研究区识别出36个关键生态节点。关键生态节点是研究区生态网络的核心组成部分，其生态功能的完好，有利于区域生态物质、能量和信息的交换，对于维持区域物种多样性、保持生态环境稳定、防止沙漠扩张具有至关重要的作用，应着重加强对关键生态节点的生态建设与保护，实行较严格的生态管控措施，防止关键生态节点被破坏。绝大部分关键生态节点位于研究区中部，如扎木温多尔、嘎咕塔拉、阿茨嘎查以及纳林湖附近，少量位于东部如奈伦湖国家湿地公园及黄河附近，南部仅有个别点位分布，位于哈拉特尔瑞辉附近。

经提取的骨架廊道基本保留了研究区生态廊道的骨架部分，细碎部分在提取过程中被舍弃。对研究区骨架廊道的长度分布进行统计可知，骨架廊道总长度为1 236.89 km。因此应对骨架廊道进行重点建设与保护，尤其是廊道110、160、235、297、304、311、457、487。

3 结论

(1)考虑了生态网络中的相互作用，并使用引力模型来量化生态源地之间的相互作用。基于引力模型和MCR模型，构建了沙漠-绿洲交错带FG生态网络多场景模拟模型，用于提取不同情景下的生态网络。

(2)从(1.0,0)模式到(0,1.0)模式，随着经济发展比例的增加和生态保护比例的下降，FG生态网络逐渐被破坏。在(0,1.0)模型中，FG生态网络的结构被破坏，FG生态网络的连通性、稳定性和核数达到最低水平。这表明，在目前的自然条件下，干旱地区的大规模经济发展将对FG生态网络造成巨大破坏。研究发现(0.9,0.1)模式的FG生态网络更复杂，网络的稳定性和连通性指标更大。然而，随着经济的发展，FG生态网络的结构遭到越来越严重的破坏。这表明研究区可能有经济发展空间，但空间不大。在干旱地区，经济发展与生态破坏之间存在平衡。应根据当地条件发展经济，避免过度开采造成的生态恶化，加剧荒漠化。

(3)根据潜在生态网络拓扑结构特征，识别出36个关键生态节点；利用最小生成树算法，识别出研究区骨架廊道总长度为1 236.89 km。对于关键生态节点的生态建设，应实行较严格的生态管控措施，防止因关键生态节点被破坏导致的区域生态恶化。骨架廊道是研究区生态廊道的主体部分，是网络其他部分连通的基本保障，因此应对骨架廊道进行重点建设与保护。