2011—2020年全球表面温度较1850—1900年平均水平上升约1.1 ℃，全球变暖趋势日益加剧^［1］。气温每升高1 ℃，饱和水汽压增加约7%，大气水汽含量随之上升，从而导致全球极端降水事件的频率、强度和持续时间普遍增强^［2］。在中国西北干旱区，1961—2020年气温以约0.33 ℃/10a的速率持续升高；自2000年以来，降水以29 mm/10a的速率显著增加，呈现明显的暖湿化特征^［3］。持续的气候扰动显著抑制植被的光合效率^［4-5］，削弱生态系统生产力并改变生态系统结构，致使生态系统稳定性的不确定性提高^［1］，从而对生态系统的服务功能、区域的生态环境和社会经济的可持续发展构成系统性挑战。因此，探究气候变化对生态系统稳定性的影响及其驱动机制，对提升生态系统应对气候变化的能力，保障生态环境和社会经济可持续发展具有重要意义。

生态系统稳定性是指生态系统在长期发展与演化过程中维持或恢复其结构与功能相对稳定的能力^［6-8］。近年来，生态系统稳定性主要集中在研究方法和影响机制的研究上。在研究方法上，建立了以抵抗力、恢复力与时间稳定性为核心的多维度量化框架^［9］，并提出了综合评价法、试验观测法和遥感分析法等评估生态系统稳定性方法。综合评价法通过构建多维度指标体系，并结合层次分析法^［10-12］、模糊综合评价^［13］、主成分分析^［14］等方法进行评估，能全面评估生态系统的多维状况，但指标选取和权重设定较为主观；试验观测法依托野外调查和控制试验，利用物种组成^［15］、生物多样性^［16-17］、种群动态^［18］等数据解释生态系统的抵抗力与恢复机制，适合小时空尺度机理研究；遥感分析法则利用归一化植被指数（Normalized Difference Vegetation Index，NDVI）、增强型植被指数（Enhanced Vegetation Index，EVI）、太阳诱导叶绿素荧光（Solar-Induced Chlorophyll Fluorescence，SIF）、总初级生产力（Gross Primary Productivity，GPP）、叶面积指数（Leaf Area Index，LAI）等指标表征生态系统稳定性^［19-20］，覆盖范围广、时间序列完整，适合揭示大尺度生态系统稳定性的格局和趋势。另外，探讨了气候、生物、人类活动对生态系统稳定性的影响。其中气候因子是最核心的外部驱动力^［8，21］。气候变暖会改变物种组成与群落结构，进而影响种间作用与生态系统稳定性^［22］；长期的水分胁迫则会导致植被代谢受限、光合效率下降，从而降低生态系统的抵抗力和恢复力^{［21，23-24］}。生物多样性则通过调节群落功能波动与养分循环增强生态系统稳定性^［25-26］，但物种同步性和种群变异性的增加可能导致生态系统稳定性下降^［27］。人类活动（如大规模的土地利用变化和森林砍伐）则导致生境丧失、生态功能退化，从而削弱生态系统的稳定性^［28-29］，而合理放牧等适度干扰又能缓解种间竞争，提升和恢复群落的稳定性^［30-31］。尽管相关研究取得了很大的进展，但是生态系统稳定性的概念和度量方法仍不统一，多重扰动的交互效应与稳定性下降区的精准识别仍缺乏系统性的探索，区域间的差异比较研究相对匮乏^［32］。典型区域尤其是干旱区高寒山地生态系统，由于自然环境复杂、观测条件受限，现有研究对其生态系统稳定性的变化特征及其驱动机制仍缺乏系统认识，亟需对相关问题进行系统分析。

祁连山是中国西北干旱区生态安全的重要屏障，特别是其乔木林生态系统起着核心作用。祁连山的乔木林生态系统主要分布在黑河、石羊河和疏勒河的河川径流形成区，承担着涵养水源、调节径流的关键功能^［33］，在水土保持、区域气候调节^［34］和固碳增汇^［35］中也发挥着重要作用，素有“绿色水库”之称。气候变化将改变乔木林生态系统的结构、组成和功能^［36-37］，进而可能降低其生态系统稳定性，削弱其涵养水源的功能，加剧水资源短缺^［38］，从而威胁河西走廊生态环境的安全和社会经济的可持续发展。因此，选择祁连山国家公园为研究对象，开展气候变化对乔木林生态系统稳定性影响的研究，以揭示祁连山乔木林生态系统稳定性的时空变化规律，探究气候变化对其生态系统稳定性的影响机制，将对建立科学的干旱区山地乔木林生态系统应对气候变化的管理策略有重要意义。

祁连山国家公园，地处青藏、蒙新、黄土三大高原交会地带的祁连山北麓，处于甘肃、青海交界处，总面积5.02×10⁴ km²（图1）。研究区由西北-东南走向的高山、沟谷和山间盆地组成，平均海拔超过4 000 m，最高点为疏勒南山团结峰（5 808 m），属高原大陆性气候区。研究区东部气候湿润，多年平均降水量为400 mm，西部气候干燥，多年平均降水量为200 mm，降水集中在6—9月，年均气温-2~6 ℃。地带性植被有山地草原、温带灌丛、山地乔木林、亚高山灌丛、高山亚冰雪稀疏植被等类型，有维管植物95科451属1 311种，有野生脊椎动物28目63科294种。针叶林建群种有青海云杉（Picea crassifolia）、祁连圆柏（Juniperus przewalskii），阔叶林树种有山杨（Populus davidiana）、白桦（Betula platyphylla）等。

研究所用遥感数据为美国地质调查局（USGS，https：//www.usgs.gov）提供的30 m空间分辨率的Landsat系列影像，时间为1987—2021年，包括Landsat 5 TM和Landsat 8 OLI/TIRS的Collection 2， Tier 1 Level-2表面反射率产品。数据通过辐射校正、大气校正、几何校正和云掩膜处理后，用于植被指数计算与时间序列分析。对于云层遮盖严重、影像缺失的年份，使用相邻年份影像通过插值法获得。

气温与降水数据来源于国家冰川冻土沙漠科学数据中心（http：//www.ncdc.ac.cn）发布的已处理高分辨率气候栅格产品，经过重采样后用于后续气候因子与生态系统稳定性分析。该数据基于英国东英吉利大学发布的 CRU TS v4.02 月尺度气候数据（1941—2021年，空间分辨率为0.5°），结合30 m 分辨率的数字高程模型（Digital Elevation Model，DEM）与祁连山及周边地区19个气象站的观测资料，通过地理加权回归模型（Geographically Weighted Regression Model，GWR模型）进行空间插值与精度提升，最终生成空间分辨率为1 km 的月均气温与降水数据。

首先，基于遥感解译的土地利用数据^［39］，采用8邻接原则将乔木林像元聚合为斑块。然后，采用空间距离阈值法将相距小于1 000 m的乔木林斑块合并，并进一步以流域范围为生态边界，将合并后的斑块整合为初始分析单元。最后，结合地形特征和土地利用背景信息，对整合后的初始分析单元的边界进行优化，剔除过小和孤立斑块，确保分析单元在空间上相对完整、功能上相对一致，形成了空间尺度较大、生态功能明确的乔木林景观分析单元Ⅰ~Ⅷ。其次，利用R语言工具计算香农多样性指数（Shannon's Diversity Index， SHDI）、辛普森多样性指数（Simpson's Diversity Index， SIDI）、香农均匀度指数（Shannon's Evenness Index， SHEI）、景观形状指数（Landscape Shape Index， LSI）、平均斑块周长（Mean Patch Perimeter， PERIM_MN）、边界密度（Edge Density， ED）、斑块数（Number of Patches， NP）、斑块密度（Patch Density， PD）、平均斑块面积（Mean Patch Area， AREA_MN）、最大斑块指数（Largest Patch Index， LPI）、平均最近邻距离（Mean Euclidean Nearest Neighbor Distance，ENN_MN）、平均连续度指数（Mean Contiguity Index，CONTIG_MN）等景观格局指数^［40］，以定量分析分析单元的景观格局特征。

采用连续小波分析（Continuous Wavelet Transform， CWT）对1991—2021年的气温与降水时间序列进行多尺度多时频分析，以揭示气候变化的动态过程。

式中：x（t）为气温或降水序列；ψ为Morlet母小波；s为尺度参数，控制提取周期的长短；τ为时间平移参数；ψ^* 为母小波的复共轭函数。

采用线性回归模型刻画气候长期趋势。

式中：X_t 为第t年的气温或降水；b为回归斜率，表征气候因子的变化速率（℃·a⁻¹或 mm·a⁻¹）。

气候弹性指数（Climate Elasticity， CE）是指 NDVI 变化中未能由气温和降水的线性关系解释的波动，相对于气候因子总体变异性的比值，该指数可用于反映植被变化中气候因子之外的综合扰动效应^［41-43］。基于以上定义，在给定的滑动时间窗口（5 a）内建立 NDVI 与气温、降水之间的线性回归方程，并基于回归残差构建气候弹性指数。

式中：T_i，j，t 与P_i，j，t 分别为第t年时间滑动窗口内的气温与降水的Z-score标准化值；α、β₁ 、β₂ 为回归系数；ε_i，j，t 为残差；σ（ε_i，j，t ）为残差的标准差，反映NDVI中非气候驱动的波动性；σ（T_i，j，t ）和σ（P_i，j，t ）为气温和降水的标准差，表示气候因子的总体变异性。CE值越大，说明NDVI的变化中非气候驱动因素贡献越高。

计算归一化植被指数（NDVI）在每个滑动时间窗口（5 a）内的变异系数（Coefficient of Variation，CV），以CV作为表征乔木林生态系统的时间稳定性的指标。

式中：CV_i，j，t 表示像元（i，j）在第t个滑动时间窗口内的NDVI变异系数；σ_i，j，t 和μ_i，j，t 分别为该时间窗口内的NDVI的标准差与均值。CV值越小，说明植被波动性越低，生态系统越稳定^［44-45］。

式中：NIR和RED 分别为近红外和红光波段反射率。

研究区内98.2%的乔木林生态系统多年NDVI变异系数小于0.223，该阈值可作为稳定性分级的重要参考^［46-47］。为增强CV值在低值区的区分度，对乔木林生态系统的时间稳定性进行负对数转换，从而定义乔木林生态系统的森林稳定性指数（Forest Stability Index， FSI）为：

式中：FSI_i，j，t 值越大，代表在t时间段内像元（i，j）的稳定性越高。

结合CV值的分布特征，将FSI划分为五个等级：极高稳定性（≥2.0）、高稳定性（1.5~2.0）、中等稳定性（1.0~1.5）、较低稳定性（0.5~1.0）及极低稳定性（FSI<0.5，表1）。

在乔木林生态系统的约470万个像元中，基于FSI分值区间进行分层抽样，总体抽样比例为 0.5%，以1991年为基准抽取23 476个样点。随后在1991—2021年提取相同位置的逐年数据，在排除空值后共得到727 687条样本记录。基于样本数据进行线性回归拟合，以判定乔木林生态系统稳定性变化趋势。

式中：FSI_t 表示年份t的稳定性值；α为回归截距；β为趋势斜率；ε_t 为残差。在显著性水平α=0.05下采用双侧t检验评估拟合模型的显著性。

采用全局莫兰指数（Global Moran's I）进行空间自相关分析，以量化乔木林生态系统稳定性在空间上的分布特征。

式中：N为空间分析单元总数；x_i 为第i个分析单元的属性值；w_ij 为空间权重矩阵元素。采用标准化Z值检验莫兰指数的显著性，当∣Z∣>1.96时，表明空间自相关在α=0.05水平显著。

以乔木林生态系统的森林稳定性指数（FSI）为因变量，气温、降水和气候弹性指数为自变量，逐像元建立基于时间序列的多元线性回归模型（Pixel-wise Multiple Linear Regression），用于定量评估各因子对稳定性变化的影响。

式中：FSI_i，j，t、Z_{CE，i，j，t​}、Z_{T，i，j，t​}、Z_{P，i，j，t} 分别表示像元（i，j）在年份t的乔木林生态系统稳定性指数、气温、降水和气候弹性指数，经逐像元的Z-score标准化后的值；γ₁、γ₂、γ₃为回归系数，分别表征气温、降水和气候弹性指数对乔木林生态系统稳定性的边际作用； ϵi,j,t为残差项，表示模型未能解释的稳定性波动。此外，通过计算气温（T）、降水（P）和气候弹性指数（CE）的方差膨胀因子（Variance Inflation Factor，VIF）对多元共线性进行检验，各自变量均满足VIF<5，不存在显著多重共线性。

首先，在ArcGIS平台中对所有栅格数据进行空间预处理，主要包括数据裁剪、投影转换、空间重采样及掩膜处理，以统一不同数据源的空间分辨率与坐标系统。

然后，在 Python环境中开展核心的数据计算与统计分析工作。利用Rasterio库完成栅格数据的读取与写出，借助NumPy和SciPy等库实现指标计算、变量标准化及逐像元多元线性回归建模等，并对回归结果进行显著性检验（t检验，P<0.05）。

此外，景观格局指数的计算在R语言环境中完成。通过调用landscapemetrics等成熟的景观生态学分析包，批量计算香农多样性指数（SHDI）、辛普森多样性指数（SIDI）、景观形状指数（LSI）、边缘密度（ED）和斑块数量（NP）等景观指标，并对计算结果进行汇总与输出。

研究区乔木林在8个分析单元间差异显著，总面积为9 509.93 km²，其中乔木林面积1 652.25 km²，平均覆盖率17.37%。分析单元Ⅶ乔木林面积最大，为795.94 km²，占全部乔木林面积的48.17%，覆盖率21.95%，是主要分布区；分析单元Ⅰ次之，面积287.41 km²，占全部乔木林面积的17.40%，覆盖率16.20%；分析单元Ⅳ和Ⅱ面积分别为190.34 km²和172.01 km²，占全部乔木林面积的比值分别为11.52%和10.41%，覆盖率16.68%和21.78%；分析单元Ⅷ 面积88.27 km²，占全部乔木林面积的5.34%，覆盖率6.23%；分析单元Ⅲ、Ⅴ和Ⅵ的乔木林面积较少，分别为52.82、41.79 km²和23.68 km²，覆盖率分别为29.57%、12.92%和9.18%（表2）。

分析单元Ⅱ、Ⅶ、Ⅰ、Ⅴ、Ⅲ、Ⅷ、Ⅳ和Ⅵ的香农多样性指数依次降低，分别为1.37、1.34、1.34、1.26、1.25、1.22、1.19、0.90。分析单元Ⅶ、Ⅱ、Ⅰ、Ⅲ、Ⅷ、Ⅴ、Ⅳ和Ⅵ的辛普森多样性指数依次降低，分别为0.70、0.69、0.68、0.66、0.65、0.64、0.60、0.44。分析单元Ⅱ、Ⅲ、Ⅰ、Ⅶ、Ⅴ、Ⅳ、Ⅷ和Ⅵ的香农均匀度分别为0.66、0.64、0.61、0.61、0.60、0.57、0.56、0.43。结果表明，分析单元Ⅰ、Ⅱ和Ⅶ的景观多样性水平最高；分析单元Ⅲ、Ⅳ、Ⅴ和Ⅷ的景观多样性水平适中；分析单元Ⅵ最低，其群落结构相对简单。

分析单元Ⅲ、Ⅱ、Ⅴ、Ⅳ、Ⅶ、Ⅰ、Ⅷ和Ⅵ的平均形状指数分别为1.84、1.81、1.80、1.77、1.77、1.75、1.69和1.64。分析单元Ⅲ、Ⅶ、Ⅱ、Ⅰ、Ⅳ、Ⅴ、Ⅵ和Ⅷ的平均斑块周长分别为3.35、2.55、2.28、2.22、2.13、2.06、1.51、1.39 km。分析单元Ⅲ、Ⅱ、Ⅳ、Ⅶ、Ⅰ、Ⅴ、Ⅵ和Ⅷ的边界密度分别为4.68、4.07、3.41、3.30、2.85、2.77、1.73、1.55 km·km^-²。结果表明，分析单元Ⅱ、Ⅲ的平均形状指数、平均斑块周长和边界密度均较高，说明其斑块形态复杂、边界效应突出，是空间异质性最突出的区域；分析单元Ⅰ、Ⅳ、Ⅴ和Ⅶ处于中等水平，其斑块形态和边界特征相对适中；分析单元Ⅵ、Ⅷ的最低，其斑块形态规则、边界效应弱，是空间异质性最小的区域。

分析单元Ⅶ、Ⅰ、Ⅳ、Ⅷ、Ⅱ、Ⅴ、Ⅵ和Ⅲ的乔木林斑块数量依次为4 865、2 341、1 808、1 651、1 453、435、301、256个。分析单元Ⅱ、Ⅳ、Ⅲ、Ⅴ、Ⅶ、Ⅰ、Ⅵ和Ⅷ的斑块密度依次降低，分别为1.84、1.58、1.43、1.35、1.34、1.32、1.17、1.16 个·km^-2。分析单元Ⅲ、Ⅶ、Ⅰ、Ⅱ、Ⅳ、Ⅴ、Ⅵ和Ⅷ的平均斑块面积依次为0.20、0.16、0.12、0.12、0.10、0.09、0.08、0.05 km²。分析单元Ⅲ、Ⅰ、Ⅶ、Ⅵ、Ⅳ、Ⅱ、Ⅴ和Ⅷ的最大斑块指数分别为35.03%、21.23%、15.63%、12.76%、11.59%、9.79%、9.27%、3.53%。分析单元Ⅵ的平均最近邻距离为0.18 km，分析单元Ⅴ、Ⅷ的均为0.17 km，分析单元Ⅰ的为0.15 km，分析单元Ⅲ、Ⅶ、Ⅳ、Ⅱ的平均最近邻距离均为0.13 km。分析单元Ⅴ、Ⅵ、Ⅲ、Ⅰ、Ⅳ、Ⅶ、Ⅷ、Ⅱ的平均连续度指数分别为0.60、0.59、0.59、0.58、0.58、0.57、0.57、0.57。结果表明，分析单元Ⅵ和Ⅷ的斑块密度最低，平均斑块面积最小，且斑块孤立，破碎化程度最高；分析单元Ⅱ和Ⅳ的斑块密度较高，斑块数量也较多，但平均斑块面积较小，最大斑块指数偏低，说明该分析单元中斑块普遍较小，破碎化程度较高；分析单元Ⅴ的斑块数量少但斑块密度高，同时斑块面积普遍较小且孤立，因此其破碎化程度亦偏高；分析单元Ⅰ和Ⅶ的斑块数量最多，但规模适中或存在较大斑块，其破碎化程度属中等水平；分析单元Ⅲ以少数大斑块主导，最大斑块指数最高，破碎化程度最低。

小波重构趋势结果表明1991—2021年乔木林生态系统的气温总体呈波段式上升（图2）。1991—2002年为快速升温阶段，气温由6.91 ℃增至7.54 ℃，年均升温速率为0.057 ℃·a⁻¹；2003—2010年气温相对平稳，维持在7.47~7.55 ℃；2011-2021年气温再次加速上升，由7.49 ℃增至7.85 ℃，升温速率为0.035 ℃·a⁻¹。2006年以后，气温以3~5 a周期的高频振荡为主导，在此周期范围内的平均功率为0.49，3.6 a周期时达到最大平均功率（0.63）。在2008—2015年形成强能量聚集区，表明此阶段气温波动最为强烈。乔木林生态系统在1991—2021年的气温增幅为0.024 ℃·a⁻¹（表3）。单元Ⅲ、Ⅱ和Ⅳ的升温最为强烈，升温速率分别为0.036、0.033、0.035 ℃·a⁻¹，分析单元Ⅰ的升温速率仅为0.016 ℃·a⁻¹；分析单元Ⅵ、Ⅶ的升温速率分别为0.028、0.016 ℃·a⁻¹，比分析单元Ⅲ的升温速率分别低0.0080、0.020 ℃·a⁻¹；分析单元Ⅷ的升温速率最低，仅为0.015 ℃·a⁻¹，比分析单元Ⅲ低0.021 ℃·a⁻¹。各分析单元的气温年际波动性也存在差异，气温标准差从大到小依次为分析单元Ⅲ、Ⅳ、Ⅱ、Ⅴ、Ⅵ、Ⅷ、Ⅰ、Ⅶ，其大小分别为0.51、0.50、0.47、0.47、0.42、0.42、0.40、0.35 ℃，表明分析单元Ⅲ的气温波动性最大，而分析单元Ⅶ的气温波动性最小。

1991—2021年，乔木林生态系统的年降水量呈现下降-上升-下降的波动格局（图3）。小波重构趋势显示：1991—1999年的年降水量由306.70 mm降至260.94 mm，年均速率为-5.72 mm·a⁻¹，为下降期；2000—2011年的年降水量由265.33 mm增至峰值388.19 mm，年均速率为11.17 mm·a^-1，为持续增湿阶段，是整个时段内降水量增加速率的最高水平；2011年后年降水量再度减少，从388.19 mm降至2021年的223.70 mm，年均速率为-16.45 mm·a⁻¹。在1991—1995年和2014—2021年年降水存在高强度能量聚集，分别集中在2~4 a和4~8 a周期，对应时段和周期内的平均功率分别为0.25和0.28。这表明降水在中短周期尺度上波动显著。这说明区域降水变化在长期趋势之外，受多尺度气候系统调控，叠加了明显的周期性扰动。

尽管降水存在明显的周期性波动，但在1991—2021年总体呈上升趋势（表3）。乔木林生态系统降水量的年均增幅为7.32 mm·a⁻¹，较国家公园的平均降水量增幅（1.43 mm·a⁻¹）高5.89 mm·a⁻¹。分析单元间的降水量的年均增幅也存在差异，从高到低依次为分析单元Ⅱ、Ⅷ、Ⅲ、Ⅳ、Ⅰ、Ⅵ、Ⅶ、Ⅴ，其年降水量的增长速率分别为9.76、8.62、7.70、7.48、6.88、6.47、5.53、2.91 mm·a⁻¹。各分析单元的降水年际波动性也存在差异，降水标准差从大到小依次为分析单元Ⅱ、Ⅷ、Ⅲ、Ⅳ、Ⅰ、Ⅶ、Ⅵ和Ⅴ，其大小分别为187.85、152.68、145.90、131.52、116.46、103.75、96.07、76.45 mm，表明分析单元Ⅱ的降水波动性最大，分析单元Ⅴ的降水波动性最小（表3）。

图4显示，乔木林生态系统稳定性在研究期内呈显著上升趋势（∣t∣=335.92>1.96，n=727 687，α=0.05），年变化率为0.027 a⁻¹。因模型的决定系数（R²）为0.134，时间因素能够解释约13.4%的长期变化，其余波动主要来源于气候因子、人类活动以及生态过程等复杂驱动。

1 276.36 km²的乔木林生态系统处于极高稳定性，占乔木林总面积的77.25%；347.14 km²处于高稳定性，占乔木林总面积的21.01%；中等及以下稳定性等级总计仅28.91 km²，占比1.75%（表4，图5）。

不同分析单元的乔木林生态系统稳定性水平存在明显差异。在极高稳定性占比超过90%的分析单元中，分析单元Ⅰ、Ⅱ和Ⅲ分别有278.90、157.89、50.34 km²面积的乔木林生态系统处于极高稳定性，分别占其分析单元内乔木林生态系统面积的97.04%、91.79%和95.31%；分别有8.16、14.02、2.47 km²面积处于高稳定性，占比为2.84%、8.15%和4.67%；分别有0.35、0.10、0.01 km²面积处于中、低稳定性，占比为0.12%、0.06% 和0.02%。这些分析单元整体表现出高度稳定的格局。在极高稳定性占比70%~90%的分析单元中，分析单元Ⅳ和Ⅴ分别有168.11 km²和35.51 km²面积的乔木林生态系统处于极高稳定性，占其分析单元内乔木林生态系统面积的88.32%和84.97%；分别有22.10 km²和6.24 km²面积处于高稳定性，占比为11.61%和14.93%；有0.11 km²和0.04 km²面积处于中、低稳定性，占比为0.06%和0.10%。这类分析单元极高稳定性虽略低于前一类，但总体仍以极高和高稳定性为主。在极高稳定性占比不足70%的分析单元中，分析单元Ⅶ、Ⅵ和Ⅷ分别有529.46 km²、14.29 km²和34.01 km²面积的乔木林生态系统处于极高稳定性，占其分析单元内乔木林生态系统面积的66.52%、60.36%和38.53%；分别有253.83、8.85、42.89 km²面积处于高稳定性，占比为31.89%、37.37%和48.59%；分别有12.66、0.54、11.37 km²面积处于中、低稳定性，占比为1.59%、2.27%和12.88%。这些分析单元的稳定性水平整体偏低，尤其是分析单元Ⅷ，稳定性分布最为分散，中、低稳定性面积占比较高。

空间自相关分析结果表明，分析单元Ⅵ、Ⅷ、Ⅶ、Ⅰ、Ⅴ、Ⅱ、Ⅳ和Ⅲ的全局莫兰指数依次降低，分别为0.89、0.87、0.82、0.74、0.72、0.72、0.68、0.61（∣Z∣>1.96），说明FSI在研究区内各分析单元中空间聚集显著（图6）。单元Ⅵ、Ⅶ和Ⅷ的聚集性最强，但以低值-低值集聚为主；单元Ⅰ、Ⅱ、Ⅲ、Ⅳ和Ⅴ的聚集性相对较弱，但多表现为高值-高值聚集。

经显著性检验表明，气温和降水的变化对不同乔木林生态系统分析单元稳定性的影响存在差异。分析单元Ⅲ、Ⅰ、Ⅱ、Ⅳ、Ⅴ、Ⅵ、Ⅷ和Ⅶ的气温标准化值与FSI的回归系数依次降低，分别为0.12、0.12、0.11、0.09、0.07、0.07、0.04、0.03，因此，气温每升高1 ℃时分析单元Ⅲ、Ⅰ、Ⅱ、Ⅳ、Ⅴ、Ⅵ、Ⅷ和Ⅶ的FSI的变化量分别为0.24、0.30、0.23、0.18、0.15、0.17、0.10、0.09（表5）。气温标准化序列与FSI的相关性显著的像元比例较高，在α=0.05下，分析单元Ⅰ、Ⅱ和Ⅲ的显著像元的比例分别占到35.65%、31.43%和42.50%，而分析单元Ⅳ、Ⅴ、Ⅵ的显著像元的比例相对较低，分别为27.20%、27.13%和23.98%。分析单元Ⅶ和Ⅷ的显著像元的比例最低，在α=0.05下仅为7.08%和11.50%（图7）。

降水标准化值与FSI的回归系数从大到小依次为分析单元Ⅶ、Ⅲ、Ⅱ、Ⅷ、Ⅳ、Ⅴ、Ⅵ、Ⅰ，分别为0.07、0.04、0.04、0.03、0.02、0.02、0.00、-0.03，因此，降水量每增加1 000 mm时分析单元Ⅶ、Ⅲ、Ⅱ、Ⅷ、Ⅳ、Ⅴ、Ⅵ、Ⅰ的FSI变化量分别为0.68、0.27、0.21、0.2、0.15、0.26、0.00、-0.26（表5）。降水标准化序列与FSI的相关性显著的像元比例低于温度的，且空间差异明显，其在α=0.05下的相关性显著的像元比例从大到小依次为分析单元Ⅰ、Ⅱ、Ⅶ、Ⅲ、Ⅳ、Ⅵ、Ⅷ和Ⅴ，分别占到28.13%、23.20%、20.56%、19.25%、18.82%、16.30%、10.51%、7.50%（图7）。

分析单元Ⅷ、Ⅱ、Ⅴ、Ⅲ、Ⅵ、Ⅰ、Ⅳ和Ⅶ的气候弹性指数与FSI的回归系数分别为-0.25、-0.27、-0.28、-0.31、-0.31、-0.31、-0.32、-0.34，因此，CE值每增加0.10时分析单元Ⅷ、Ⅱ、Ⅴ、Ⅲ、Ⅵ、Ⅰ、Ⅳ和Ⅶ的FSI的变化量分别为-1.40、-1.15、-1.35、-1.00、-1.15、-1.25、-1.10、-1.75（表5）。气候弹性指数标准化序列与FSI的相关性显著（α=0.05）的像元比例最高，其在α=0.05下的相关性显著的像元比例从大到小依次为分析单元Ⅳ、Ⅶ、Ⅲ、Ⅴ、Ⅰ、Ⅵ、Ⅱ、Ⅷ，分别占到97.36%、96.73%、96.56%、95.03%、93.88%、91.12%、88.80%、87.65%（图7）。

经过多元统计回归分析发现气温标准化序列与FSI的偏R²在分析单元Ⅲ、Ⅰ、Ⅱ、Ⅳ、Ⅴ、Ⅵ、Ⅶ和Ⅷ的值依次为0.13、0.11、0.10、0.09、0.09、0.08、0.05和0.05，说明温度对FSI的贡献率仅占5%~13%。降水标准化序列与FSI的偏R²在不同分析单元间差异显著，从大到小依次为分析单元Ⅰ、Ⅱ、Ⅶ、Ⅵ、Ⅳ、Ⅲ、Ⅷ和Ⅴ，分别为0.10、0.09、0.08、0.07、0.07、0.07、0.05和0.04，说明降水对FSI的贡献率仅占4%~10%。气候弹性指数标准化序列与FSI的偏R²在所有分析单元中均为最高，从大到小依次为Ⅶ、Ⅳ、Ⅴ、Ⅲ、Ⅰ、Ⅵ、Ⅱ和Ⅷ，分别为0.47、0.44、0.43、0.43、0.39、0.38、0.33和0.33，说明气候弹性指数对FSI的贡献率达到33%~47%。但是各分析单元的多元统计模型的R²在0.41~0.56，说明气温、降水和气候弹性指数可以解释FSI的41%~51%的变化（表5）。

本研究结果表明，1991—2021年祁连山国家公园乔木林生态系统中极高稳定性面积占比77.25%，且乔木林生态系统整体的FSI以0.027 a⁻¹（α=0.05，t_α=±1.96<335.92，n=727 687）的速率显著提升，这一趋势与区域内其他生态系统功能的研究一致。李玉辰等^［48］指出，祁连山国家公园的NDVI和NPP整体上升，其中87.29%的区域NPP显著提高，生态系统服务功能（如生境质量、土壤保持、水源涵养和碳储存等）亦显著增强^［49-50］。Zhang等^［51］基于NDVI计算青藏高原植被覆盖度的研究也表明，2000—2015年青藏高原60.30%的面积上植被覆盖度呈增加趋势，且显著增加区域主要分布在青藏高原北部。

在研究时段内，回归结果显示气候变暖促进了乔木林生态系统稳定性的提升。而气温升高对生态系统稳定性的正效应可能与光合作用增强和生长季延长有关。温度上升在一定范围内可以提高光合作用效率和生产力，进而提升森林的抵抗力和恢复力^［52］。Wang等^［53］利用北美洲的森林永久样地的研究同样证明，气温升高能够促进树木的生长，提高森林的生产力，从而增强生态系统稳定性，但白云杉（Picea glauca）、山杨（Populus tremuloides）、杰克松（Pinus banksiana）的生产力在夏季最高温分别到达23.4、27.0、26.2 ℃时开始下降，说明升温对生态系统稳定性的影响温度存在阈值。过度升温可能导致干旱胁迫和呼吸消耗增加，从而削弱稳定性^［4］。Ma等^［54］在青藏高原高寒草地的试验中发现，气候变暖显著降低了该地区生态系统的时间稳定性。Zhao等^［55］通过评估温度升高背景下的19个浮游生物食物网的稳定性，发现温度升高会降低其生态稳定性，这意味着生态系统稳定性在不同群落中对气温变化的响应可能存在差异。本研究中气温变化对生态系统稳定性的影响为“正效应”，这可能是气温变化仍处于乔木林生态系统可承受温度的适宜区间，而未来持续升温是否对生态系统稳定性造成负面影响还需进一步研究。

本研究结果显示，在分析单元Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅶ、Ⅷ中降水量的增加促进了生态系统稳定性的提升。这与现有研究结论一致，即在干旱与半干旱地区，降水量增加通过提高植被的平均生产力、物种异步性等，进而增强群落的时间稳定性^［56-58］。Liu等^［59］通过对青藏高原高寒草地生态系统的79个样地开展实地调查获得的植物物种多样性、年净初级生产力和土壤理化性质的研究，指出降水量的增加通过提高物种多样性和增强地下生态系统多功能性，间接促进群落生产力的提升，从而提高生态系统稳定性。Smith等^［60］的全球尺度研究也表明，植被稳定性随着水分可利用性的提高而增强。然而在分析单元Ⅵ中，降水量对生态系统稳定性的影响相对微弱；在分析单元Ⅰ中，降水量增加对生态系统稳定性反而表现为负效应（-0.26）。这一现象表明，降水对稳定性的作用并非单向度的正效应，而是受到区域背景水热条件的调控^［61］。在本研究中，单元Ⅰ位于祁连山国家公园的最东端，其年降水量高于中西部地区，在这一相对湿润的本底条件，额外的降水可能带来多重负面影响。其一，持续增加的降水可能造成土壤过湿甚至渍涝，导致根系缺氧、病原滋生与个体抗逆性的下降^［62-63］；其二，降水量的持续增加还会加剧土壤侵蚀和养分流失，从而削弱群落维持稳定性的能力^［64］；其三，降水量增加可能伴随云量增加、辐射减弱，从而导致区域光照条件恶化，限制光合作用，降低植被生产力、削弱生态系统稳定性^［65］。

本研究发现，气候弹性指数在所有分析单元中对乔木林生态系统稳定性的影响均为负效应，其CE值每增加0.10，FSI值为-1.00~-1.75。这一结果表明，FSI除受到气温和降水的影响外，还受到非气候因子的制约。分析单元Ⅶ和Ⅷ中CE值增加对FSI影响的负效应最强，这些区域景观破碎化程度较高、空间异质性较低，说明景观条件差可能放大了非气候扰动对FSI的影响^［66］。相比之下，在景观多样性高、异质性强的分析单元（Ⅰ、Ⅱ、Ⅲ），CE 虽仍表现为负效应，但数值较小（-1.00~-1.25），表明景观结构较好时可以在一定程度上缓冲非气候扰动的不利影响^［67］。

本研究结果显示，气温与降水对乔木林生态系统稳定性的独立解释力相对有限，其与FSI显著相关的像元比例分别为7.08%~42.50%和7.50%~28.13%，对FSI的贡献率仅为5.00%~13.00%和4.00%~10.00%，与已有研究结果一致，即单一气候因子往往难以充分解释生态系统过程的主要变异^［68-70］。Li等^［71］的研究进一步显示，全球中高纬地区春季复合暖干（CWD）、温度和降水对夏季植被变异的显著相关面积分别为27.00%、15.00%和24.00%，也表明了复合指标的解释力通常优于单因子。本研究中的气候弹性指数（CE）正是利用残差法构建的用于代理非气候因素的综合指标，其与FSI的显著相关的像元比例为87.65%~97.36%，对FSI的贡献率为33.00%~47.00%，明显高于气温的或降水的。多元回归模型的R²为0.41~0.56，表明气温、降水和CE等3类因子能够在一定程度上解释乔木林生态系统稳定性的变化，但仍有44%~59%的未解释部分。这是因为CE仅是对非气候因素的近似度量，不足以完全涵盖植被结构、土壤环境、地形格局和人类活动等因素^［72-73］。因此，未来研究应纳入更多直接表征非气候因素的指标，以进一步提升模型对生态系统稳定性变异机制的解释能力。

尽管本研究分析了祁连山国家公园乔木林生态系统稳定性对气候变化的响应，但仍存在一定局限性。首先，乔木生态系统稳定性主要基于NDVI时间序列构建，虽能有效表征植被功能状态的时间稳定性，但对乔木林生态系统的结构差异和生物量变化的刻画仍较为有限。其次，本文重点关注年尺度气温和降水的影响，未进一步区分极端气候事件及季节尺度气候因子的作用，可能低估其对生态系统稳定性的扰动效应。此外，气候弹性指数（CE）作为非气候因素的综合代理指标，难以区分不同非气候驱动因子的具体作用机制。未来研究可结合森林结构、土壤和人类活动等多源数据，以进一步提升对乔木林生态系统稳定性变化机制的解释能力。

近30年来祁连山乔木林生态系统稳定性显著增强，说明干旱区森林在这30年表现出一定的恢复力与适应性。

气温升高对祁连山乔木林生态系统稳定性的影响呈现正效应，但降水对乔木林生态系统稳定性的作用在空间上存在差异，其在大部分区域对乔木林生态系统的稳定性提升有正向作用，而对部分区域的影响呈现负效应。

气候弹性指数在祁连山所有分析单元均与乔木林生态系统的稳定性呈负相关，表明非气候因素对乔木林生态系统稳定性的影响具有抑制作用。

乔木林生态系统稳定性对气候变化的响应具有显著的空间差异，该特征表明在不同区域背景下，稳定性的维持过程可能存在差异性，这为后续开展分区尺度的生态管理研究提供了科学依据。