基于复杂网络理论的西藏巴宜区森林景观空间结构研究

卢 杰 王 戈 马 骏 于 强 高 郯 牛 腾

(1.西藏农牧学院高原生态研究所, 林芝 860000； 2.西藏高原森林生态教育部重点实验室， 林芝 860000；3.北京林业大学精准林业北京市重点实验室， 北京 100083)

0 引言

西藏巴宜区森林资源极为丰富，森林覆盖率达到46%，是中国第三大林区，拥有藏西80%的森林蓄积量，丰富的森林资源构成了西藏巴宜区独有的森林景观。西藏巴宜区雨量丰富，森林景观作为重要的水源涵养系统，其削洪补枯能力以及防止水土流失、降低山体滑坡风险的作用日益凸显。在较大的空间尺度上，森林景观的空间结构即空间分布格局决定其生态功能能否正常发挥。对区域森林景观结构进行准确评估是目前的研究热点与重点，目前，许多学者利用GIS空间分析方法、地学信息图谱法、景观生态分析法和空间格局分析法等对区域植被格局、景观生态格局和森林景观格局等进行了研究，但研究方法主要侧重于空间景观结构的分析，对景观斑块之间的关系缺乏定量化的表达，尤其是森林景观斑块之间的关联[1-3]。虽有研究者利用耦合网络模型研究森林小班之间的联系，但利用相邻小班建立联系主观性过强，并没有从数据本身反映相邻小班是否有信息的传递和能量的交换[4]。复杂网络可以将复杂系统简化成节点和边的集合，节点代表复杂系统的基本单元，边代表单元之间的关系。利用复杂网络理论的分析方法可以很好地反映空间节点与空间节点之间的联系，可以定量描述区域森林景观斑块之间的内在联系[5-7]。

本文以西藏巴宜区为研究区，针对该地区多发的水土流失和滑坡泥石流现象，从森林景观网络的连通性和抗破坏能力角度出发，利用复杂网络理论的分析方法，在区域整体尺度研究其森林景观空间网络的复杂系统特性，旨在反映西藏巴宜区森林景观空间的网络特性，为后续的森林景观结构空间格局研究提供指导。

1 研究区概况

巴宜区位于西藏东南部(35°2′～95°17′E，29°21′～30°15′N)，行政区域总面积为1.02×104km2。南邻米林县，西部和西北部与工布江达县相连，东北部与波密县接壤，东临墨脱县(图1)。全区平均海拔3 300 m，巴宜区南部地区属于冈底斯山脉的余脉，北部属于念青唐古拉山支脉的高山地段。巴宜区的气候类型为温带湿润季风气候，雨水充沛，日照时间较长。年平均温差差异较小，但是日均温差差异则较大。雅鲁藏布江自西向东流经林芝市东南部，雅鲁藏布江支流帕隆藏布江和尼泽河由西北往东南汇入雅鲁藏布江。巴宜区森林植被受海拔和气候的影响，在分布上呈明显的规律性变化特点。

2 数据及研究方法

2.1 数据来源与处理

选取巴宜区夏季且少云的Landsat OLI影像(2018年)为研究素材，空间分辨率30 m的地面高程数据(DEM)作为辅助数据。利用ENVI 软件对影像进行波段合成、图像增强和几何校正处理，选择最大似然监督分类法对遥感影像进行目视解译，提取巴宜区的景观类型信息，使用ArcMap软件进行细碎斑块处理，运用叠加分析工具进行空间数据分析，最终在Arc Info环境下完成拓扑和改错处理，针对影像分类精度进行实地验证，分类误差在3个像元以内。

2.2 森林景观空间网络识别与构建

区域内森林景观斑块之间存在物质和能量的流动，即各个森林景观斑块是相互联系的，其联系的紧密程度决定了区域森林景观系统发挥生态系统服务功能[8]。因此通过基于复杂网络理论和生态阻力面模型识别西藏巴宜区森林景观斑块之间的联系，进而构建西藏巴宜区森林景观空间网络。

2.2.1生态阻力面模型

生态阻力面模型能够从宏观上反映生态过程的趋势，可以从地形坡度、植被覆盖、水文分布、土地覆盖类型4方面构建森林景观生态阻力评价体系。坡度越平缓、植被覆盖度越高、距水体越近生态能量流动越畅通，农田与建筑建设用地对生态能量流动具有刚性阻碍[9]。按照表1(表中NDVI为归一化植被指数，MNDWI为归一化差异水体指数)将各项生态阻力划分为5个等级，分别用1、3、5、7、9来表示阻力值，利用ArcGIS软件制作各个因子的评价结果，进行叠加栅格计算得到生态阻力基面的综合评价结果，利用ArcGIS软件中cost-distance模块完成累积阻力面计算。

2.2.2生态源地识别

生态源地一般是指物种栖息地，是生态安全格局中的核心区域，也是各种生态流的源头，它的提取是生态安全格局构建的第一步[10]。通过ArcMap软件统计森林景观斑块面积，利用Zonal Statics工具进行平均NDVI/MNDWI提取以及采用Fragstats软件计算斑块的形状指数，通过熵值法确定以上4个要素的权重，确定计算各景观斑块重要程度，将重要性统计结果中前15%的森林景观斑块识别为生态源地。

表1 阻力评价体系Tab.1 Resistance evaluation system

斑块形状指数用于衡量斑块形状的复杂程度[11]。将相同面积的规则图形(圆形或正方形)作为参照物，通过计算斑块形状与参照之间的差异程度来判断斑块形状的复杂程度。一般来说，景观斑块形状越复杂，其物质能量信息交换等生态过程越频繁[12]。斑块形状指数计算公式为

(1)

式中Sshape——斑块形状指数

P——斑块周长A——斑块面积

2.2.3森林景观生态廊道识别

森林景观生态廊道是森林景观斑块之间的物质传递和能量流动的通道，利用最小累积阻力模型提取西藏巴宜区森林景观生态廊道[13-14]。森林景观生态廊道具有提高生境连接度、保护生物多样性、防止水土流失等作用。根据西藏巴宜区森林景观生态阻力评价体系，构建阻力面，提取最小阻力路径，生成潜在生态廊道[15-16]。模型一般公式为

(2)

式中VMCR——最小累积阻力值

fmin——土地单元的累积阻力最小值

n——生态源地数目

m——土地单元数目

Dij——从生态源地j到土地单元i的空间距离

Ri——土地单元i对物种运动的阻力

2.2.4生态脆弱节点识别

生态节点是相邻源地间连通的重要支撑点，是生态廊道中对于物质、能量流动和促进生态功能的健康循环具有重要作用的景观组分[17-18]。脆弱节点是生态廊道累积阻力最大的区域，也就是最薄弱的区域，其生态功能相对较弱，但对于保持景观格局的连通度至关重要。本文将生态累积阻力面作为基础地形数据，利用ArcGIS水文分析工具提取累积生态阻力栅格面的山脊线，累积生态阻力值最大的点就位于这条山脊线上，山脊线和生态廊道的交点就是生态廊道最薄弱的区域，即为生态脆弱节点[19]。

2.2.5空间网络构建

基于复杂网络理论，将西藏巴宜区森林景观生态源地和生态脆弱节点作为森林景观空间的节点，将森林景观生态廊道作为森林景观空间的连接边，构成了“点-线-面”交织的森林生态网络。

2.3 森林景观空间网络特性分析

2.3.1度及度分布

生态节点度是与该生态节点相连廊道数量，节点度越大表明该生态节点重要性越高。网络平均度是网络中每个生态节点度平均值[20]。

在复杂生态空间网络中，节点的度在统计学上服从一定分布函数，生态节点度分布p(k)代表在一个生态网络中，度为k的节点占比，度分布函数p(k)是节点度为k的节点被抽到的概率。

说起从省城来到农村，每天和老百姓同吃同住同劳动，“吃苦受累倒不怕，就怕做不出成绩，3年时间一晃就过去了，我们更应该齐心协力、只争朝夕。找到合适的工作方法，为乡村振兴做一些实事，无愧于青春、无愧于事业，也就不枉驻村一场了。”林燕玲说。

2.3.2聚类系数

复杂生态网络平均聚类系数为0～1，C=0表明所有节点没有任何边连接；C=1表明网络中任意两节点均直接相连；C越大，表明生态网络中节点联系越紧密，聚集系数越大的网络小世界特性越强，反之越弱[22]。

2.3.3空间网络生态节点介数

网络生态节点介数是生态网络中任意两个生态节点之间最短生态廊道路径所通过的生态节点数量，生态节点介数可以反映生态网络中生态节点的重要性[23]。生态节点介数Bi计算公式为

(3)

式中njl——生态节点vj与vi之间的最短生态廊道数量

njl(i)——生态节点vj与vi之间最短生态廊道路径经过生态节点vi数量

N——节点个数

森林景观空间网络是维持区域生态环境稳定的重要保障，完整的空间结构能保证其发挥正常功能[24]。森林景观空间网络结构鲁棒性是用来衡量森林景观空间网络在遭受到外界干扰破坏时，维持其正常结构和功能的能力，即潜在森林景观空间网络抵抗能力，与之相对应的就是潜在森林景观空间网络恢复能力。在其空间结构遭到破坏后，潜在森林景观空间网络能够恢复的能力即鲁棒性。结构鲁棒性指数公式为

R=c/(N-Nr)

(4)

式中Nr——从森林景观空间网络中去除的节点个数

c——当节点被去除后森林景观空间网络最大连通子网络中的生态节点数量

针对生态节点和生态廊道，恢复鲁棒性指数计算公式分别为

D=1-(Nr-Nd)/N

(5)

E=1-(Mr-Me)/M

(6)

式中D——节点恢复鲁棒性指数

Nd——通过某种策略恢复节点个数

E——边恢复鲁棒性指数

Mr——从网络中去除的边个数

Me——通过某种策略恢复的边数量

M——初始网络中边的数量

通过构建实际森林景观空间网络表示节点关系的邻接矩阵来评价森林景观空间网络的结构鲁棒性，制作无向无权网络拓扑图，并对网络进行随机攻击和恶意攻击[25]。随机攻击即从网络中随机去除若干个节点，恶意攻击即从潜在森林景观空间网络中，同时去除度最大的Nr个节点及其对应边。

3 结果与分析

3.1 森林景观空间网络构建

3.1.1生态源地提取

基于2018年遥感影像分类结果共提取出森林景观斑块共453个，面积为3 085 km2，并根据熵值法计算得到影响森林生态源地的斑块面积指数、斑块形状指数和NDVI的熵值权重分别为0.5、0.3和0.2。依据熵值法统计各景观斑块重要程度，确定出巴宜区森林生态源地68个(图2)，巴宜区森林生态源地分布很不平衡，主要分布在尼洋河附近以及东北方向的山谷，海拔为1 000～2 600 m。

3.1.2累积阻力面构建

基于研究区基本情况，选取高程、土地利用类型、归一化差异水体指数(MNDWI)、归一化植被指数(NDVI)来构建巴宜区生态阻力评价体系(图3)。巴宜区高程为1 032～7 326 m，中部和南部海拔较高，东部海拔较低。巴宜区森林景观主要分布在河流两岸，尼洋河横穿巴宜区在东南部汇入雅鲁藏布江，未利用地主要为常年积雪的高海拔山地，积雪融化导致在未利用地周围形成大量水体景观。研究区内高山积雪区MNDWI最高，尼洋河附近MNDWI较高，草地景观MNDWI较低。巴宜区内森林景观NDVI较高，水体景观与积雪区NDVI最低。

考虑4方面阻力因子进行森林生态阻力面模拟，并叠加生成巴宜区生态阻力面(图4)，森林景观生态阻力值较低，东部生态阻力值最低，为10，积雪区生态阻力值较高，最高为32。基于已识别的森林生态源地斑块和生态阻力面，运用cost-distance模型进行生态能量最小累积阻力模拟。

3.1.3生态廊道提取

基于最小累积阻力面与森林生态源地，运用cost-path模型模拟生态能量从源地出发，克服生态阻力到达另一个生态源地的最大概率流动路径。研究区内共提取66条森林生态廊道(图5)，在南部与东部源地斑块距离较近、生态廊道较短，中部与西南部廊道较长。

3.1.4生态节点提取

生态廊道与最小累积阻力面的阻力值山脊线的交点即确定为生态节点，是生态网络的最薄弱、最至关重要的点。研究区内共提取出潜在生态节点22个(图6)，主要分布在研究区中部与东南部附近。

3.2 森林生态网络结构评价

3.2.1度及度分布

由图7可知，该森林生态网络的度为1的节点最多，有38个，即该38个森林生态源地只与一条生态廊道连接。度最大值为19，仅有一个森林生态源地。对森林生态网络度分布散点图进行曲线拟合，拟合结果表明幂律分布特征更为明显，Poisson分布特征不明显，表明该森林生态网络无标度性特征要强于均匀性特征。

3.2.2森林生态网络节点重要性分析

由图8可知，该森林生态网络中大部分森林生态源地的介数较小，有38个森林生态源地介数为0，介数为0的森林生态源地度为1，当其中任意一个遭到破坏时不会对整个森林生态网络的连通性造成影响。介数在0～400之间的森林生态源地有29个，介数大于400的森林生态源地有1个，其度为19，在整个森林生态网络中重要性最高。

3.2.3聚类系数

由图9可知，该森林生态网络聚类系数为0的森林生态源地斑块占95%，达到65个，这些森林生态源地与其他森林生态源地不具有集群特点。该网络中聚类系数最大值为0.35，具有集群特点，表明该网络具有非均匀特性，网络结构稳定性较差。

3.2.4森林生态网络结构鲁棒性分析

恶意攻击和随机攻击下的边恢复鲁棒性均呈现下降趋势(图10)，随机攻击下，当去除的边小于18时，其网络结构能够恢复。而随机攻击下随着边去除数目的增加，边恢复鲁棒性的下降呈现增加的趋势。在恶意破坏下，当除去的边小于3时，网络完全可以恢复，对森林生态网络边进行恶意攻击对其网络结构的破坏非常明显。

当去除节点的数目较小时，遭到破坏的节点可以完全恢复。随着破坏的节点数目的增加，恶意攻击和随机攻击的节点恢复鲁棒性均呈下降趋势，且均为凸曲线。随着网络中随机攻击去除的节点数目超过58个，越来越多的丢失节点得不到恢复。而恶意攻击下去除超过60个节点后丢失节点得不到恢复，随机攻击下的节点恢复鲁棒性要优于恶意攻击。

4 结论

(1)基于西藏巴宜区具体情况，选取斑块面积指数、斑块形状指数、NDVI和MNDWI 4方面指标，运用熵值法共筛出68个森林生态斑块作为生态源地。

(2)选取高程、NDVI、MNDWI与土地利用类型作为阻力因子进行空间叠加分析，基于最小累积阻力模型提取66条森林生态廊道，识别森林生态廊道上脆弱生态节点22个，在县域范围内构成了“点-线-面”交织的森林生态网络。

(3)基于复杂网络原理构建西藏巴宜区森林生态复杂网络，对森林生态网络结构进行评价。结果表明，该地区森林生态网络度分布呈现幂律分布，具有无标度性特征；该网络具有非均匀特性，网络结构稳定性较差；在随机攻击下，边与节点的恢复鲁棒性优于恶意攻击下的恢复鲁棒性。