引言

祁连山位于中国甘肃省与青海省交界处，是青藏高原东北部的重要山脉，地理位置重要，地貌类型复杂，具有显著的生态和环境意义。作为中国西部的生态屏障之一，祁连山在维持区域生态平衡、调节气候以及保护水资源方面起到了至关重要的作用。近年来，由于气候变化和人类活动的加剧，祁连山地区的生态环境面临着严重的威胁，尤其是土地利用的变化和景观格局的破碎化已成为显著的生态问题。研究祁连山的景观格局演变，能够为该地区的生态环境保护与资源合理利用提供科学依据。通过分析和预测未来景观格局的变化趋势，可以为决策者提供切实的管理建议，促进可持续发展，防止生态环境的进一步恶化。本研究的主要目标是利用多种遥感数据和景观分析方法，深入分析祁连山景观格局的演变过程及其未来趋势。本研究通过结合地理探测器模型、景观格局指数以及FLUS模型，不仅揭示了祁连山景观格局变化的规律，还对未来景观破碎化及多样性的趋势进行了预测。研究结果为今后的区域生态规划与管理提供了理论支持。

数据来源和研究方法

1.1数据来源和研究区概况

1.1.1 数据来源

本研究所用数据主要包括土地利用数据、地形与气候数据、社会经济数据、交通与居民点分布数据以及限制开发区域数据（表1）。其中，建设用地、居民点、铁路、公路和河流等基础地理数据来源于全国地理信息资源目录服务系统；GDP、人口、年平均降水量、高程和NDVI等数据来源于中国科学院资源环境科学与数据中心；2000年、2010年和2019年土地利用数据均由对应年份遥感影像解译获得。

为构建限制开发区域，本文引入了Geofabrik提供的OSM数据，并提取研究区内自然保护区、风景名胜区、历史遗址、宗教场所、冰川和湿地等空间要素，叠加祁连山生态保护红线范围，形成FLUS模型所需的限制发展区域数据。

为保证多源数据在后续分析中的一致性，所有数据均经过投影统一、空间裁剪、栅格化和重采样处理，最终统一为TIFF格式，空间分辨率设定为1000 m，栅格行列数为938×516，以满足FLUS模型及相关空间分析的输入要求。

1.1.2研究区域概况

祁连山国家公园位于中国甘肃省与青海省交界处，地理坐标为36°47'N-39°48'N，94°51'E-102°60'E（图1），处于青藏高原东北部，东南与永登县接壤，东北与肃南、民乐、永昌、凉州区和天祝县相连，西南与海西蒙古族藏族自治州德令哈市相接，西北与阿克塞、肃北县相邻。祁连山的总面积约为5.02×10^4 km²，常住人口约为54,665人。

祁连山地区以高山和盆地为主，海拔范围从1722米至5791米不等。该地区属于高原大陆性气候，年平均气温为4°C左右，年均降水量约为400毫米，气候干湿分明，昼夜温差较大。祁连山内水资源丰富，河流密布，具有显著的生态多样性。该区域的植被具有明显的垂直地带性，从低海拔的山地草原带、温带灌丛草原带、山地森林草原带，到亚高山灌丛草甸带和高山亚冰雪稀疏植被带，生态系统类型丰富。

1.2研究方法

1.2.1 土地利用转移矩阵分析

土地利用转移矩阵方法用于定量分析不同土地利用类型在不同时间阶段之间的转移方向和转移数量。转移矩阵的常见形式如下所示：

其中，表示从土地利用类型转移到类型的面积，代表面积的数量，为土地利用类型的总数，和（=1,2, ..., ）分别代表转移前和转移后的土地利用类型。该矩阵的构建可以揭示不同土地利用类型之间的时空转移过程，从而为后续的景观格局分析提供必要的基础数据。

在土地利用转移分析中，我们使用转移矩阵来揭示各类土地利用的扩展与减少趋势。通过计算不同时间段（如2000年和2019年）之间的转移概率，可以识别出关键的变化区域和变化趋势，为土地管理与生态保护提供数据支持。

1.2.2 景观格局指数计算方法

景观格局指数分析方法用于全面描述祁连山地区景观的空间特征，包括景观的破碎化程度、形态特征和空间异质性。在本研究中，我们选择了以下景观格局指数：斑块个数（NP）、景观类型面积比例（PLAND）、景观形状指数（LSI）、聚集度（AI）、破碎度（FN）、分维度（DIVISION）、分离度（SPLIT）、Shannon多样性指数（SHDI）以及集聚度指数（AI）。这些指数能够有效地刻画景观结构的变化，反映景观的空间特征。

为了避免信息冗余，本文采用了移动窗口法来分析景观格局的空间变化。具体方法是选择3000 m大小的窗口，每次窗口移动一个栅格，计算当前窗口内的景观格局指数。该方法可以帮助我们识别景观格局在研究区内的空间异质性，并揭示其破碎化和多样性的变化。移动窗口法的基本计算过程如下：

其中，为窗口(i,j)位置的景观格局指数，为窗口内栅格的数量，为窗口内第k个栅格的景观格局指数。通过这种方式，我们可以获得景观指数的空间分布特征图，从而分析景观格局的时空变化。

1.2.3 景观干扰度指数计算

景观干扰度指数用于量化景观类型对外界干扰的抵抗能力和自我恢复能力。景观的生态风险与人为及自然干扰的强度密切相关，干扰强度越大，生态系统的敏感性和景观的风险越高。景观干扰度指数通过结合景观破碎度、景观分离度和景观分维度等多个指标来构建，其计算公式为：

其中，, , 分别为景观破碎度、景观分离度和景观分维度的权重，赋值分别为0.5、0.3和0.2., , 分别为景观破碎度、景观分离度和景观分维度的数值。通过计算这些指数的加权和，我们可以得到每种景观类型的干扰度指数，从而量化景观的干扰程度。

在景观干扰度分析中，干扰度指数越高，表示该区域景观对外界干扰的抵抗能力较弱，生态风险较大。此分析为景观生态风险评估和保护区划提供了科学依据。

1.2.4 地理探测器分析

地理探测器（Geodetector）模型是一种用于分析多因子之间交互关系的空间分析工具。在土地利用变化研究中，地理探测器可以有效探测不同因素对土地利用变化的影响力和空间分异性。地理探测器主要包括四个子模型：因子探测器、风险探测器、生态探测器和交互作用探测器。因子探测器是用于分析不同因子对土地利用变化空间分异性的解释力；风险探测器是探测因子与土地利用情况之间是否存在风险关系；生态探测器是比较不同因子对不同区域的土地利用分布的影响；交互作用探测器则是分析单因子和双因子之间的交互作用，计算交互作用强度。若单因子作用的影响力大于双因子共同作用时的影响力，则表明该因子在变化过程中起主导作用。

交互作用探测器的计算公式为：

其中，为交互作用强度，为因子交互作用的总变异度，和分别为单因子作用时的变异度。该公式能够揭示不同因子对土地利用变化的独立和共同作用程度。

1.2.5 FLUS模型与土地利用预测

FLUS模型（基于元胞自动机的土地利用变化模拟模型）由两部分组成：适宜性概率计算模块（BP-ANN）和基于自适应惯性机制的元胞自动机模块（CA）。BP-ANN模块用于训练和预测土地利用变化的适宜性概率，CA模块则基于适宜性概率图进行土地利用的空间模拟。

BP-ANN模块的训练过程包括输入驱动因子（如地形、交通通达度、经济、人口等）并进行神经网络训练。训练样本的选择通过随机采样方式进行，训练过程的目标是最大化各类土地利用类型的适宜性预测精度。训练过程的优化目标为最小化以下误差函数：

其中，为误差函数，为实际值，为预测值，为样本总数。通过最小化误差函数，可以得到土地利用变化的适宜性概率图。

CA模块则通过元胞自动机的自适应惯性机制，模拟土地利用类型的变化，并考虑土地类型之间的相互转化难易度。模型的输入为初始土地利用数据，输出为未来土地利用格局。模拟过程基于以下转换规则：

其中，表示土地类型转换到土地类型的可能性，1表示允许转换，0表示不允许转换。通过此机制，FLUS模型能够进行未来土地利用变化的模拟和预测。

2 分析与结果

2.1祁连山土地利用变化特征

2.1.1 土地利用结构变化

2000-2019年，祁连山土地利用类型以草地、裸地和冰川—积雪地为主，三者构成研究区景观基质。其中，草地始终占主导地位，面积占比约为48%-53%；裸地次之，占比约为24%-35%；冰川—积雪地主要分布于高海拔区域。林地、耕地、水体、湿地和建设用地面积相对较小，整体呈局部集聚分布特征。

从时间变化看，研究期内土地利用结构总体保持稳定，但局部类型发生了持续调整。草地面积整体波动较小，说明其在区域景观结构中具有较强稳定性；裸地面积总体下降，表明部分裸地向草地、水体或冰川—积雪地转化；林地面积变化不大，但内部结构存在一定调整。耕地、水体和建设用地面积虽占比较低，但总体呈缓慢增加趋势。

总体而言，祁连山土地利用结构在研究期内未发生根本性改变，但草地、裸地、冰川—积雪地和水体之间存在较明显的相互转化，表明区域土地利用格局在相对稳定背景下仍表现出持续的动态调整特征。

时期	项目	冰川_积雪	草地	耕地	灌丛	建筑用地	林地	裸地	湿地	水体
2010	面积/km2	16742	98457	1362	2	56	3838	68090	335	5407
	比例/%	0.09	0.51	0.01	0.00	0.00	0.02	0.35	0.00	0.03
2010	面积/km2	11572	102664	1404	1	61	3856	68876	347	5508
	比例/%	0.06	0.53	0.01	0.00	0.00	0.02	0.35	0.00	0.03
2019	面积/km2	16038	101544	1439	2	65	3861	65254	340	5745
	比例/%	0.08	0.52	0.01	0.00	0.00	0.02	0.34	0.00	0.03

2.1.2 土地利用转移过程

土地利用转移矩阵表明，不同时段祁连山土地利用变化具有明显阶段差异。2000-2010年，土地利用转化较为活跃，主要转移方向为林地向草地转化、草地向裸地转化，以及裸地向草地转化。其中，草地与裸地之间的双向转换最为显著，表明这两类景观对区域环境变化具有较高敏感性。草地向裸地的转化主要集中在西部地区，而裸地向草地和林地的转化主要出现在中部和局部东部区域。

相比之下，2010-2019年土地利用变化幅度明显减弱，整体格局趋于稳定。除西部局部区域裸地变化相对明显外，耕地、林地、湿地和建设用地基本保持稳定。草地变化主要表现为少量草地向裸地转化，水体变化相对有限，裸地则主要向草地和冰川—积雪地发生局部转换。

总体来看，祁连山土地利用变化主要集中于草地、裸地和冰川—积雪地之间，自然景观类型之间的相互转化构成了区域景观演变的主体过程。2010 年之后祁连山土地利用类型变化较小，且变化主要集中在西部地区，裸地的变化

2.2 景观格局变化特征

2.2.1斑块尺度景观指数变化

斑块尺度分析结果表明，草地和裸地始终是祁连山地区的优势景观类型，其面积占比远高于其他地类，是维持区域景观结构稳定的主体。空间上，草地主要分布于中低海拔区域，裸地主要集中在西部及高寒干旱地区，二者均具有较高连片性。相比之下，自然过渡带和坡度较大区域的景观类型更加丰富，斑块分布更为破碎。

景观形状指数（LSI）结果显示，祁连山西北部裸地景观和东南部草地、林地景观均表现出较高的形状复杂度。2000-2010年，林地、湿地、水体和裸地的形状指数变化相对较小，而草地、灌丛和冰川—积雪地的形状指数有所下降，表明其斑块边界趋于规则。2010-2019年，草地、裸地、水体和建设用地的LSI普遍上升，表明这些景观类型的边界复杂性增强，景观分割和重组现象加剧。

总体上，斑块尺度指标表明，祁连山景观演变不仅表现为面积变化，还伴随斑块形态重塑和空间结构调整。

2.2.2 景观尺度景观指数变化

景观尺度分析结果表明，祁连山整体景观聚集度较高，说明研究区景观仍保持较好的连续性和完整性。聚集度高值区主要分布在青海湖周边及低海拔区域，反映出这些区域受单一或优势景观类型控制较强。2000-2010年，聚集度指数由77.2295上升至78.1150，表明景观聚合程度有所增强；2010-2019年，该指数下降至77.1607，说明后期景观聚集性有所减弱，但整体仍处于较高水平。

分析景观类型，其计算结果如下：

年份	香农多样性指数(SHDI)	香农均匀度指数(SHEI)	景观破碎度(FN)	聚集度(AI)
2000	1.1486	0.5227	0.0405	78.2295
2010	1.101	0.5011	0.04	77.115
2019	1.1437	0.5205	0.0445	77.1607

景观破碎度在空间上表现为山腰区域较高，并向山顶和山麓递减，说明中海拔过渡带景观镶嵌程度高，受自然条件变化和人类活动干扰更为明显。2000-2019年，景观破碎度由0.0405上升至0.0445，表明研究区景观连续性减弱，斑块切割程度增强。

香农多样性指数（SHDI）空间分布与破碎度具有较强一致性。高值区主要分布于山腰和生态过渡区，说明这些区域景观类型更丰富、异质性更强。时间上，2000-2010年SHDI略有下降，而2010-2019年有所回升，表明后期景观类型分布趋于均衡，景观异质性增强。

总体来看，祁连山景观格局在研究期内呈现出“聚集度总体较高、破碎度缓慢增加、多样性先降后升”的变化特征。

2.2.3 景观干扰度的时空变化

景观干扰度指数结果显示，祁连山景观以中度干扰和轻度干扰类型为主。2000年，中度干扰型景观面积最大，占研究区总面积的46.11%，主要分布于西北部地区；轻度干扰型景观占44.61%，主要分布于东南部和青海湖周边；重度干扰型景观面积最小，占9.28%，空间分布相对零散，主要位于西段与中段交界区域。

2000—2019年间，中度干扰型景观面积呈下降趋势，轻度干扰型景观面积持续增加，说明部分中度干扰区域向低干扰状态转化；与此同时，重度干扰型景观面积也有所增加，表明局部区域干扰强度仍在增强。总体而言，研究期内景观干扰度空间分布格局变化不大，但不同干扰等级之间发生了一定转换。

从空间分布看，高干扰区域主要对应建设用地和耕地等人为活动较强区域，低干扰区域主要对应冰雪覆盖地、水体和裸地，中度干扰区域则主要分布于草地、林地和灌丛集中区域。该结果表明，祁连山景观干扰度具有较明显的土地利用类型响应特征。

2.2.4 景观格局对景观干扰度的响应

相关分析结果表明，2000—2019年祁连山景观干扰度与SHDI和LSI呈正相关，与AI呈负相关。这说明随着干扰程度增加，景观多样性和形状复杂性增强，而景观聚集性下降。换言之，外界干扰会促进景观斑块分割，增加景观异质性，并削弱景观连通性。

其中，景观干扰度与SHDI和LSI的相关性相对较强，表明干扰对景观多样性和边界形态影响更为显著；与AI的负相关则进一步说明景观聚合程度会随干扰增强而下降。总体来看，景观干扰增强是推动祁连山景观破碎化和空间异质性上升的重要因素之一。

指标	香农多样性(SHDI)	景观形状指数(LSI)	聚集度(AI)
皮尔逊相关系数	0.967**	0.793**	-0.546**

2.3 地理探测器驱动因子分析

因子探测结果表明，不同驱动因子对祁连山土地利用空间分布的解释力存在明显差异，其排序为：NDVI >降水量> GDP >人口>气温。其中，NDVI的q值最高，为45.8%，说明植被状况是影响研究区土地利用格局空间分异的首要因子；降水量次之，表明水热条件对区域土地利用分布具有重要控制作用；GDP和人口则反映出人类活动在局部区域土地利用变化中的作用。

等级	气温/℃	人口/(人·KM²)	降水/mm	GDP/(万元·km⁻²)	NDVI
1	-12.1～ -5.5	0.1～17.7	62.6～188.3	0.1～4.6	0.0～0.1
2	-5.5～-3.5	17.7～58.3	188.3～282.6	4.6～19.2	0.1～0.3
3	-3.5～-1.5	58.3～110.2	282.6～355.2	19.2～42.1	0.3～0.4
4	-1.5～0.6	110.2～186.7	355.2～418.1	42.1～73.8	0.4～0.5
5	0.6～2.8	186.7～1009.9	418.1～490.6	73.8～639.1	0.5～0.7
6	2.8～8.9	1009.9～3302.1	490.6～679.3	639.1～2090.5	0.7～0.9

生态探测结果显示，部分因子之间存在显著差异，说明不同自然因子和社会经济因子在土地利用格局形成过程中发挥着不同作用。交互探测结果进一步表明，任意两个因子共同作用时的解释力均高于单因子作用结果，表现为双因子增强或非线性增强，说明祁连山土地利用格局演变是多因素耦合作用的结果。

变量	NDVI	TPOP	TEM	PER	GDP
q statistic	0.46	0.1	0.07	0.17	0.11
p value	0	0	0	0	0

	NDVI	TPOP	TEM	PER	GDP
NDVI		N
TPOP	N		N
TEM	N	N		Y
PER	N	Y	Y		N
GDP	N	Y	Y	N

2.4 基于FLUS模型的土地利用模拟与预测

以2000年土地利用数据为初始状态，在DEM、社会经济因子和生态环境因子的共同驱动下，对2019年土地利用格局进行了模拟。结果表明，当迭代次数设置为300次时，模拟结果与实际2019年土地利用分布具有较高一致性，Kappa系数达到0.85，总体精度达到0.91，说明FLUS模型能够较好地反映祁连山土地利用演变过程，可用于未来情景预测。

在模型验证基础上，以2019年土地利用格局为基准，对2038年土地利用进行了预测。结果显示，到2038年，建设用地增长最为明显，增长率达到32.14%；水体和耕地也呈增长趋势，增长率分别为12.11%和11.26%；裸地和冰雪覆盖地则表现为减少，变化率分别为−8.10%和−7.80%。这表明未来祁连山土地利用变化可能呈现出建设用地缓慢扩张、部分自然景观收缩调整的趋势。

从空间上看，建设用地扩张主要集中于已有居民点及交通可达性较高区域，裸地和冰雪覆盖地的变化则主要表现于高海拔及干旱区。总体而言，未来祁连山土地利用格局仍将受到自然环境条件和人类活动共同驱动。

2.5未来景观破碎化与多样性变化

根据2038年预测结果，祁连山未来景观格局将发生一定程度的重组。景观破碎度下降至0.034，聚集度指数上升至80.18，而多样性指数下降至1.11，表明未来景观整体将趋于更高聚合度和更低异质性。该结果说明部分景观斑块可能向更大尺度、更连续的空间格局整合，但景观类型组成的丰富度和均衡性有所减弱。

从空间分布看，未来景观破碎化高值区仍主要位于祁连山西段，并以东段为中心向两侧逐渐减弱，说明景观破碎化的空间分异格局在未来仍将持续。景观干扰度结果进一步表明，中段与西段过渡区域仍为高干扰区，而青海湖周边以轻度干扰为主。

总体来看，未来祁连山整体景观可能趋于聚合，但局部区域，尤其是中西段，仍存在较高生态脆弱性和干扰风险。

3讨论与结论

祁连山位于青藏高原东北缘，是我国西北地区重要的生态安全屏障，其景观格局演变受自然环境约束与人类活动扰动共同影响。研究表明，2000—2019年祁连山景观总体仍以草地、裸地和冰雪覆盖地为主，区域景观基底未发生根本改变，说明高寒山地景观格局具有较强的自然稳定性。这一特征与祁连山显著的垂直地带性、复杂地形地貌及高寒干旱气候密切相关，自然因子仍是区域景观格局形成与演变的基础力量。

尽管整体格局较稳定，不同景观类型之间仍发生了明显转化，主要集中于草地、裸地、冰雪覆盖地和水体之间，表明祁连山景观演变是在自然演替、气候波动及局地开发活动共同作用下持续调整的动态过程。景观指数分析显示，2000—2019年研究区景观破碎化程度整体上升，景观异质性增强，草地、灌丛、建设用地、裸地、冰雪覆盖地和水域等类型的形状复杂性增加，说明景观斑块边界和空间结构发生了明显重组。破碎化高值区主要分布于山腰地带和自然过渡区，反映出这些区域较高的生态敏感性。

景观干扰度分析表明，研究区以中度干扰和轻度干扰景观为主。2000年中度干扰景观面积最大，占46.11%，轻度干扰景观占44.61%，重度干扰景观占9.28%。研究期内，中度干扰景观减少，轻度干扰景观增加，重度干扰景观虽占比较低但呈上升趋势，说明祁连山整体生态格局未发生大范围恶化，但局部区域生态压力仍在增强。相关分析表明，景观干扰度与SHDI、LSI呈正相关，与AI呈负相关，说明外界扰动会增强景观异质性和边界复杂性，降低景观连通性和聚合程度，是推动景观破碎化的重要因素。

地理探测器结果表明，NDVI和降水量是影响祁连山土地利用空间分异的主导因子，GDP和人口等社会经济因子在局部区域也发挥了重要作用，各驱动因子之间普遍存在增强效应，说明祁连山景观演变是自然环境与人类活动多因子耦合驱动的结果。

FLUS模型预测显示，到2038年，祁连山建设用地、水体和耕地将增加，裸地和冰雪覆盖地将减少；景观破碎度下降至0.034，聚集度上升至80.18，多样性指数下降至1.11。未来祁连山景观整体可能趋于聚合，但中西段及其过渡区域仍是高干扰和高敏感区，应作为生态保护和景观修复的重点区域。未来土地利用优化应坚持生态保护优先原则，优先依托现有开发区域进行集约利用，严格控制新增建设活动对关键生态斑块的进一步切割与扰动，以维护祁连山景观格局稳定和区域生态安全。

参考文献：

1. [1] 靳芳,张振明,余新晓,等.甘肃祁连山森林生态系统服务功能及价值评估[J].中国水土保持科学,2005(01):53-57.
2. [2] 薛晓玉,王晓云,段含明,等.基于土地利用变化的祁连山地区生境质量时空演变分析[J].水土保持通报,2020,40(02):278-284+325.
3. [3] 于航,刘学录,赵天明,等.基于景观格局的祁连山国家公园景观生态风险评价[J].生态科学,2022,41(02):99-107.
4. [4] Liu X, Liang X, Li X, et al. A future land use simulation model (FLUS) for simulating multiple land use scenarios by coupling human and natural effects[J]. Landscape and urban planning,2017,168:94-116.
5. [5] Wang J, Zhang T, Fu B. A measure of spatial stratified heterogeneity[J]. Ecological indicators,2016,67:250-256.
6. [6] 汪孝贤,张秀霞,李旺平,等.基于遥感生态指数（RSEI）改进模型的祁连山国家级自然保护区生态环境质量评价[J].生态与农村环境学报,2023,39(07):853-863.
7. [7] Qiao B, Yang H, Cao X, et al. Driving mechanisms and threshold identification of landscape ecological risk: A nonlinear perspective from the Qilian Mountains, China[J]. Ecological indicators,2025,173:113342.
8. [8] Li Y, Qin Q, Wang D, et al. An analysis on the spatial heterogeneity characteristics of landscape ecological risk in Qilian Mountain National Park[J]. Frontiers in forests and global change,2024,7.
9. [9] 苏军德,赵晓冏,李国霞.基于InVEST-FLUS模型的祁连山国家公园碳储量演化分析及预测[J].环境工程,2024,42(07):190-199.
10. [10] Dong C, Yu H, Qian X, et al. Enhancing ecological connectivity in the Qilian Mountains: Integrating GCA and optimized MST models for ecological corridor construction[J].Ecological indicators,2024,166:112525.