基于GEE的祁连山自然保护区景观格局变化分析

刘宽梅

（兰州大学 资源环境学院，甘肃 兰州730000）

景观动态变化研究具体表现为景观结构、景观功能、景观空间格局随时间的变化，是景观生态学的核心问题[1-2]。景观格局是指在空间上大小、形状不同、属性有别的景观斑块在景观空间上的排列[3]。研究景观格局的动态过程，能掌握其动态变化规律，进一步准确认识研究区面临的各种生态问题。随着国内外景观格局研究成果不断丰富，大致可以将景观格局分析方法分为三大类：（1）景观格局指数法，是由文献[4]中Turner和Gardner、Li和Reynolds等研究发展的反映景观结构与空间配置特征的量化指标，能高度浓缩景观格局信息。Fragstats等软件的出现让指数计算更便利，如Yang C Y、Wang J J及Zhuang G S等均通过Fragstats计算相应的景观格局指数，对不同地区的景观格局进行分析[5-7]；Song H H等通过景观指数分析发现景观格局的不足从而优化景观格局[8]。（2）空间统计学方法，包含了地统计学方法、空间自相关及空间插值法等[9]。如黄玲[10]、肖飞艳[11]、潘竟虎[12]、Pan D Y[9]等利用转移矩阵法分别对不同地区景观格局变化进行了探讨。（3）景观格局模型与模拟法，便于研究者对数据隐含的内在规律进行探索[13-14]。张小松[15]、廖芳均[16]、罗攀[17]、陈颐[18]、罗紫薇[19]等基于CA-Markov模型分别对城市、自然保护区等景观格局进行预测。从研究内容来看，主要有景观格局的时空变化及其驱动力因素[20]、景观空间格局的粒度效应[21-22]、景观格局的演化过程与特征及景观格局变化的尺度效应等方面的内容[23-24]。

祁连山作为水源涵养林和生物多样性保护热点区，其生态环境问题受到广泛关注。目前，关于景观格局的研究多在祁连山东段开展。草地为该区的景观基质，其优势度高，连接性好[25-28]。斑块水平格局特征为：草地形状结构趋于松散且简单化，森林和冰雪形状结构趋于紧密但简单化，灌木形状结构趋于松散但复杂化，裸地和水域形状结构趋于紧密但简单化[27，29]。景观水平格局特征为：景观整体多样性水平增加，破碎化程度增大[28-29]，景观均匀度增大，景观优势度、分离度减小[27]。对于景观整体异质性的研究，不同学者研究时段及选取指标不同，会导致结果有所差异，刘晶[27]和刘宇[30]等对祁连山东段1994-2008年的景观格局进行探讨，刘晶等认为该区景观异质性增大，景观向着均匀化方向发展；而刘宇等认为该区景观异质性呈下降趋势；徐亚男等[29]发现2000-2016年间景观异质性是增大的。人类活动是影响该区景观格局变化、生态敏感性变化和景观脆弱性变化的决定性因素[26-29]。而对于祁连山其他区域，不同地区景观格局变化不同。祁连山地区最主要的土地利用类型是草地和未利用地[31]；祁连山南坡基质由裸地向草地转变，景观破碎化水平较高，放牧活动导致草地和灌丛破碎化，水域形状差异大，沼泽分布集中成片，居民区和工矿用地分布破碎[32-33]；中祁连山浅山区草地面积最大，受海拔和坡向的影响，各景观组分空间分布不均衡[34]，乔木林、灌木林和农田面积在增加，疏林、草地和裸地面积减少[35]；祁连山三河源地区景观格局向多样化、破碎化趋势演进，人类活动对景观格局的干扰在增强[36]；石羊河上游祁连山区草地、工矿居民用地、裸土和戈壁的破碎度增大，耕地和林地的破碎度减小，景观多样性略有上升，单一性略微下降[37]；祁连山甘肃区域景观破碎化程度较高，Ⅰ级生态景观格局比较稳定，Ⅱ级生态景观格局受自然和社会环境等因素影响变化较大[38]。

祁连山地区云量分布有明显的西部向东部逐渐减少的趋势，各季节总云量空间分布也有明显的变化，遥感影像数据质量会影响分类结果的好坏，已有研究多基于RS和GIS平台获取景观分类结果，其会受计算机空间存储和处理环境配置限制，还会受单一数据质量限制，而谷歌地球引擎（Google Earth Engine，简称GEE）云平台利用多个CPU同时对遥感数据进行运算，可实现数据快速处理，而且可以提取研究时段内多景遥感影像合成云量最小的影像，克服受云量影响难以获得高质量数据的缺陷，提高数据可用性。因此，基于GEE平台解译保护区近30年的景观类型图，并根据斑块水平和景观水平选取景观格局指数对该区域内景观格局变化进行分析，以期为保护区管护、生态规划提供参考。

1 研究区概况

祁连山自然保护区（简称保护区）位于北纬36°29′-39°43′，东经96°46′-103°45′之间，位于祁连山中、东段区域范围内，面积约282.19万公顷（包括甘肃祁连山国家级自然保护区198.72万公顷和青海祁连山省级自然保护区83.47万公顷）。祁连山东西长达1200km，南北宽50-120km，海拔3000-5000m。祁连山深居内陆，具有典型的高原大陆性气候特征，全年降水量主要集中在5-9月。祁连山气候要素的垂直地域分异明显，一般山前低山是荒漠气候，年降水量约150mm，年均温度6℃左右；山体中下部属于半干旱草原气候，年降水量250-300mm，年均温2-5℃；山体中上部为半湿润森林草原气候，年降水量400-500mm，年均温0-1℃；山体上部为亚高山和高山冷湿气候，年降水量约800mm，年均温-5℃左右。研究区位置示意图如图1所示。

图1 研究区位置示意图

2 数据来源与方法

2.1 数据来源与处理

（1）研究区边界：青海祁连山省级自然保护区边界是基于《青海祁连山自然保护区科学考察集》[39]中的边界图，对其进行地理配准矢量化得到，甘肃祁连山国家级自然保护区边界来源于国家地球系统科学数据中心（http：//www.geodata.cn）。

（2）景观分类标准体系：由于景观是具有特定结构、功能和动态特征的不同生态系统镶嵌组成的异质性区域[40]，因此分类体系以《全国生态系统分类体系》为基础，参考欧阳志云等人[41]的《基于遥感技术的全国生态系统分类体系》，结合保护区景观特征，将研究区的景观类型划分为八类，分别为湿地、高寒荒漠、冰川/永久积雪、城镇、农田、森林、草地、灌丛。

（3）遥感数据及预处理：基于GEE平台调用研究区1989、1995、2000、2005、2010、2015、2019年每年4-10月已经过预处理的Landsat地表反射率产品（Surface Reflectance，SR）（https：//code.earthengine.google.com）；利用多景影像重构及合成目标年份最小云量影像集；计算归一化建筑指数、归一化水体指数及归一化植被指数作为进行随机森林分类时的特征变量以提高分类精度；结合GEE与Google Earth平台，采用人工目视判读得到8类景观类型各年的总样本数分别为1193、1092、1100、1023、1101、1110、1161个，其中70%用于分类器训练，30%用于结果验证。

2.2 研究方法

2.2.1 基于GEE的随机森林景观类型分类

随机森林分类方法是由多棵决策树构成的机器学习算法[42]，该方法已被证明比其他方法具有更高的分类精度[43]。设置决策树数量为100对训练样本集进行训练，并提取每一目标年份的景观类型信息，以总体精度、Kappa系数作为衡量景观分类精度的指标，具体各指标计算公式参见相关文献[44-45]，七个时期的总体精度分别为：90.8%、91.21%、91.78、85.81%、88.90%、83.71%及89.66%，其Kappa系数分别为0.88、0.88、0.87、0.82、0.85、0.79及0.87。

2.2.2 景观格局指数

景观格局指数能够高度浓缩景观格局信息[4]，而且其计算简单，现被广泛使用。根据以往的研究来看，指数间有很高的相关性，若是引入同一类型指数太多，容易造成信息冗余[46]。因此，采用基于统计的景观指数选择方法，该方法能解决个人选择指数的主观性，降低选取大量景观指数进行计算的工作量，具有一定的合理性[17]。以“自然保护区景观格局”为主题检索词，在知网上检索2015－2020年的论文数据（共计76篇）；然后统计论文中出现的景观指数，最后提取到使频次为10次及以上的景观指数17个，对入选指数进行皮尔逊相关分析，最后综合考虑指数生态学意义、不重复性、常用性及研究需求，斑块类型水平选择4类特征指标共8个：（1）优势度指标：最大斑块面积指数（LPI）、平均斑块面积（MPS）；（2）破碎化指标：斑块密度（PD）、边缘密度（ED）；（3）形状特征指标：景观形状指数（LSI）、分维数（PAFRAC）；（4）聚集度指标：散布与并列指数（IJI）、聚合度指数（AI）。景观水平选择4类特征指标共6个：（1）聚集度指标：聚合度指数（AI）、散布与并列指数（IJI）；（2）优势度指标：最大斑块面积指数（LPI）；（3）形状特征指标：景观形状指数（LSI）；（4）景观多样性指标：香农多样性指数（SHDI）、香农均匀度指数（SHEI）。指标含义见表1，指数计算采用Fragstats4.2软件完成，计算方法参考Fragstats4.2软件帮助文件（http：//www.umass.edu/landeco/research/fragstats/fragstats.html）。

表1 景观格局指数意义

3 结果与分析

3.1 景观类型结构变化

由图2可知，保护区景观类型以高寒荒漠和草地为主。1989-2019年间，高寒荒漠面积总体上呈减少趋势，净减少面积为3130.45km2，其优势地位慢慢被草地代替，说明荒漠化环境有所改善；草地面积在各研究时段内波动变化，但总体是上升的，净增加面积为4223.62km2，说明草地景观比较脆弱，极易受放牧活动等人类活动的影响；森林面积经历先减少后增加再减少的趋势，总体上面积是减少的，净减少面积为261.09km2，灌丛面积变化呈现先增加后减少的趋势，总体上面积是减少的，净减少面积为242.47km2，可以看出植被结构遭到破坏，出现退化现象；冰川/永久积雪面积呈现缓慢减少趋势，净减少面积为662.18km2；湿地、城镇与农田面积占比均较小，面积变化不明显。

图2 1989-2019年祁连山自然保护区景观类型面积

3.2 基于景观指数的景观格局分析

3.2.1 类型水平的景观指数分析

（1）优势度分析

由表2可知，最大斑块面积和平均斑块面积最大值由高寒荒漠和草地主导，2019年之前，高寒荒漠是保护区的基质景观，但在1989-2019年间，由于封沙育林育草、黑土滩治理、荒漠化治理等环境保护政策和手段的实施，使得高寒荒漠被草地等景观分割的程度逐渐增加，其最大斑块面积指数和平均斑块面积值不断减小，破碎化程度增大，优势度降低，对景观的控制作用正在削弱，而草地的最大斑块面积指数和平均斑块面积不断增加，优势度则不断增强，逐渐代替高寒荒漠成为新的优势景观，说明荒漠化环境得到改善。其他景观类型的最大斑块面积指数值变化相对较小，结合其平均斑块面积变化，可以发现冰川/永久积雪和森林的优势度次于高寒荒漠和草地，灌丛和湿地的优势度较低，城镇和农田的优势度最低。

表2 最大斑块面积指数和平均斑块面积

（2）破碎度分析

由表3可知，高寒荒漠和草地的斑块密度值和边缘密度值比其他景观类型的值高，说明要素斑块与其相邻异质斑块间的接触多，边缘复杂，破碎化程度较其他景观类型高。1989-2019年间，高寒荒漠和草地的斑块密度值波动幅度小，但高寒荒漠的边缘密度值呈下降趋势，说明其边缘由复杂向简单化演变，景观的破碎化一定程度上在减轻，草地的边缘密度值呈上升趋势，其边缘复杂化，破碎度增加，道路交通及放牧等人类活动导致其破碎度增加；灌丛的边缘密度值在不断下降，斑块密度值先上升后下降，说明该景观的破碎化程度在由高向低转变；冰川/永久积雪、农田的边缘密度值在不断下降，但冰川/永久积雪的斑块密度值变化不大，农田的斑块密度呈下降趋势，这说明这两种景观要素的破碎化程度好转；森林的斑块密度和边缘密度值的波动趋势一致，均呈“倒N型”分布，破碎化程度好转；湿地和城镇的斑块密度和边缘密度值均较小且波动程度也小，其破碎化程度较低，这与城镇呈团聚式发展有关。

表3 斑块密度和边缘密度指数

（3）形状特征分析

由表4可知，1989-2019年间所有景观类型的形状指数都较大，说明研究区内各景观斑块形状相对复杂。结合表5可看出，各类景观类型的形状指数和周长面积分维数变化规律基本一致。具体来看：草地的形状指数值和周长面积分维数指数值均较大且变化平稳，说明草地斑块形状复杂度较大且较稳定；灌丛、农田、冰川/永久积雪的形状指数和周长面积分维数值均呈下降趋势，说明这三类景观类型的形状在趋于简单化；森林和高寒荒漠的形状指数在增加，周长面积分维数基本不变，说明其形状在趋于复杂；湿地和城镇的斑块形状指数和周长面积分维数均呈波动变化状态，但形状指数是除冰川/永久积雪外最小的，说明这两类景观在所有景观中的形状最简单。

表4 景观形状指数与周长面积分维数

（4）聚集度分析

从表5可以看出，聚合度指数最大的是高寒荒漠，其次为冰川/永久积雪和草地，城镇和农田最低。说明高寒荒漠、草地和冰川/永久积雪聚集程度高，趋向连片分布，前两者的聚集度高与其面积占优势有关，而冰川/永久积雪则往往在高海拔地区成片分布，而城镇和农田趋于分散分布，这与保护区内城镇分布实际情况相符。在1989-2019年间，高寒荒漠、冰川/永久积雪、草地和森林的聚合度指数基本没发生变化，说明这四类景观各自聚集度高且稳定；湿地和城镇的斑块聚集度呈增加趋势；而灌丛和农田的聚合度指数呈波动变化，灌丛的斑块分布由聚集走向分散再走向聚集，2010年以后趋于平缓，而农田的聚集度则是波动下降的。高寒荒漠和草地的散布与并列指数值比较大，说明这两者与其他景观斑块间的连通度较大，这与高寒荒漠和草地在景观中的优势地位有关，但在1989-2019年间，它们与其他斑块的连通度都呈先上升后下降的趋势；灌丛和农田与其他的景观斑块类型的连通度也是下降的；湿地和城镇与其他斑块的镶嵌程度不断增加；而冰川/永久积雪的散布与并列指数尽管有上升趋势，但仍然是最低的，说明它与其他斑块的连通度最小，这跟该景观类型分布在高海拔地区有关，除高寒荒漠外，其余景观类型极少分布在冰川/永久积雪所在的高海拔位置。

表5 聚合度指数和散布与并列指数

3.2.2 景观水平的景观指数分析

（1）聚集度分析

由图3可知，散布与并列指数表现出逐渐下降的趋势，说明各种景观类型之间空间邻接分布的复杂程度有所减小，彼此之间相邻的景观要素类型由多变少，空间关系趋于简单化，景观的连通性变弱。通过聚合度指数变化趋势可以看出，该值在研究时段内聚合度指数值均较大（均大于88），且有缓慢增长趋势，说明同类斑块趋于集中分布，景观破碎化减小，景观间的空间异质性逐渐缩小。

图3 景观水平上的散布与并列指数和聚合度指数

（2）景观优势度及形状特征分析

由图4可知，最大斑块面积指数呈现出先增加再减小再增加的变化趋势，进一步验证高寒荒漠景观优势度由强变弱，草地景观优势度由弱变强的发展过程。1989-1995年的景观形状复杂度减小，1995-2005年景观复杂度呈增加趋势，2005年以后景观斑块形状在趋于简单化发展，说明景观形状变化较频繁，人类活动对景观格局的影响呈现出有目的、有规划改造的趋势。

图4 景观水平上的最大斑块面积指数和景观形状指数

（3）景观多样性分析

由图5可知，香农多样性指数和香农均匀度指数值变化趋势一致，1995-2010年间香农多样性指数和香农均匀度指数值呈缓慢增长趋势，表明该时段内景观多样性呈增长趋势，且各景观类型分布趋向均匀化发展，1989-1995年和2010-2019年间两个指数为下降趋势。但从总体来看，香农多样性指数值在1.14-1.24之间，反映出研究区景观类型多样性低，所有景观类型在研究区内呈现不均衡化分布趋势，香农均匀度指数值为0.55-0.60之间，说明各景观类型所占比例差异较大，景观受优势景观支配，进一步验证了本保护区内高寒荒漠、草地景观占主导地位的现状。

图5 景观水平上的香农多样性指数和香农均匀度指数

4 结论

（1）保护区景观类型以高寒荒漠和草地为主。1989-2019年间，高寒荒漠面积总体上呈减少趋势，净减少面积为3130.45km2，其优势地位慢慢被草地代替，荒漠化环境有所改善；草地面积总体上是波动上升的，净增加面积为4223.62km2，说明草地景观比较脆弱，极易受人类活动尤其是放牧活动的影响；森林面积总体上是减少的，净减少面积为261.09km2，灌丛面积总体上是减少的，净减少面积为242.47km2，可以看出植被结构遭到破坏，出现退化现象；冰川/永久积雪面积呈现缓慢减少趋势，净减少面积为662.18km2；湿地、城镇与农田面积占比均较小，面积变化不明显。

（2）斑块类型水平景观格局变化：30年间，优势景观由高寒荒漠向草地过渡；草地的斑块密度值和边缘密度值在各年都是最大的，其破碎度最大且呈增加趋势；形状指数值和周长面积分维数值最大值由草地和灌丛主导，两者的斑块形状复杂度较大，草地的斑块形状复杂度稳定，而灌丛的斑块形状由复杂趋向简单化发展；高寒荒漠和草地在景观中的优势地位使得高寒荒漠和草地的聚合度和连通度均比其他景观高，两者的聚合度变化不大，而连通度都呈先上升后下降的趋势，冰川/永久积雪因为其分布海拔高且成片分布，致使其聚合度跟草地相当但连通度低。

（3）景观水平格局变化：景观整体间的连通度在减小，同类景观趋于集中分布，景观破碎度在减小；景观斑块形状总体在趋向简单化发展；景观类型多样性低，景观类型呈现不均衡化分布趋势，景观受优势景观所支配。

5 讨论

有研究表明，草地是祁连山大部分区域内的基质[25-27]，本文中高寒荒漠占比最大的原因是由于本研究区域青海境内大片高寒荒漠分布，这在汤萃文[32]等对祁连山南坡植被景观格局的研究中得到印证，而农田、城镇及水体小面积零散分布与的王清忠[35]的研究结果相同。刘晶[27]，蒋刚[33]等学者对祁连山不同地区的研究都表明，人类活动是景观频繁变化的动因，甚至使某些区域的植被类型出现严重退化现象，不同地区受到的人类活动影响不同，各研究学者选取的景观格局指数及研究时段也不同，所以对景观格局变化结果也会有差异。如辛雨琼[26]等发现东祁连山农牧交错区1988-2018年景观呈复杂化和破碎化趋势，刘晶[28]等发现祁连山东段1994-2008年间景观多样性、均匀度、破碎度指数增大，景观优势度、分离度减小，但是徐亚男[29]等发现祁连山东段2000-2008年景观整体多样性水平降低，异质性和破碎化程度减小，各景观组分分配不均匀且朝着聚集的趋势发展，2008-2016年，景观整体多样性水平增加，异质性和破碎化程度增大，各景观组分向着均匀化且分散的趋势发展。综上可以看出，即使是同一个大区域，也可能会因为指标选取及小区域内人类活动影响程度不同等因素导致景观格局结果有差异，因此祁连山这样的山地景观格局指标选取还有待进一步探讨。此外，研究区内植被存在抵御能力弱及退化现象，仍需对其进行加强管护。