替代稳态下阜康北部荒漠生态弹性的时空格局

许毓哲，林 涛，李 君,4

（1.中国科学院新疆生态与地理研究所，荒漠与绿洲生态国家重点实验室，新疆 乌鲁木齐 830011；2.中国科学院大学，北京 100049；3.自然资源部荒漠-绿洲生态监测与修复工程技术创新中心，新疆 乌鲁木齐 830022；4.中国科学院新疆生态与地理研究所，阿克苏绿洲农田生态系统国家野外科学观测研究站，新疆 乌鲁木齐 830011）

生态弹性是衡量生态系统稳定性的重要指标，其定义为在转换到一个由不同结构和功能的替代稳态之前，系统所能承受干扰的大小［1-3］。以系统变量变化趋势为基准，系统变量围绕变化趋势上下波动，其实际数值偏离变化趋势的时间越长，生态弹性越弱。而系统变量存在多种偏离变化趋势的情况，并会表现出稳定性，这就是替代稳态。系统变量受正面环境扰动时高于变化趋势变化的替代稳态称为稳定状态I，反之为稳定状态II［3-4］。不同替代稳态下生态系统表现出的弹性不同，其系统变量所处的波动区间差异明显，这种现象普遍存在于湖泊、森林、荒漠等各类生态系统［5］。

古尔班通古特沙漠是发育典型的温带荒漠系统，沙漠南缘植被受一定的绿洲排水补给，沙漠腹地植被主要依赖少量降水生存。阜康北部的荒漠植被，其空间格局在整个古尔班通古特沙漠具有代表意义。研究表明，该地区植被对气候变化反应十分敏感［6］，过去数十年呈现恢复-退化-恢复交替变化趋势［7］，意味着该区域在不同替代稳态间反复转化。在变化环境中，生态弹性是生态系统维持功能稳定的重要机制。以植被功能的时间序列为切入点，探讨阜康北部荒漠生态系统不同替代稳态下生态弹性的时空属性特征以及生态弹性内在的作用机制，在荒漠生态系统过程认识、区域生态保护和修复方面兼具重要的理论和实践意义。

地球大数据相关技术的研究和构建，为利用遥感数据时间序列量化生态系统弹性提供了可能性。在区分固有的季节变化和环境干扰方面最重要的方法就是BFAST(Breaks for Additive Season and Trend)方法［8-10］。该方法在检测植被功能变化趋势方面已得到广泛应用，并可用于量化生态系统弹性［11-12］。通过时间序列量化生态系统弹性的代表性方法分为两类，一类是自相关系数法，另一类是自回归系数法。前者通过计算系统变量的时间自相关系数量化弹性［12-15］；后者则是计算系统变量自回归系数量化弹性［16］。两种方法均建立在生态系统只存在单一稳态的前提下，所以无法量化替代稳态下的生态弹性。最近，生态系统从不同稳定状态退出所用的时间解决了替代稳态条件下生态弹性的量化问题，并应用于湖泊生态系统藻类生物量的生态弹性量化中［17-19］，但是尚未应用于空间尺度上的生态弹性研究。Meng等［20］利用Landsat时间序列量化出不同时期湘江流域生态弹性，并对弹性的时空变化做了初步探索，但是缺少对替代稳态下弹性的差异性和其变化驱动因子的深入探讨。

为了深入认识阜康北部荒漠生态系统弹性的时空演变，本文选择了一条从沙漠南缘至沙漠腹地的样带，运用2001 年1 月至2020 年12 月MODIS 的全球植被指数遥感数据产品，将该时间序列产品分成四段不同时期，并使用BFAST和状态空间建模对数据进行处理和提取，通过退出时间法计算了替代稳态条件下的生态弹性，并进一步分析降水变化与生态弹性时空演变的内在机制。研究能够为量化荒漠生态系统弹性及保护生态系统稳定性提供参考。

1 研究区概况与方法

1.1 研究区概况及样带选择

研究样带南起阜康北部古尔班通古特沙漠南缘，北至沙漠腹地，地理位置位于87°42′～88°12′E、44°21′～45°21′N（图1），南北长111.3 km，东西宽40.1 km，总面积约4.47×103km2。样带内发育南北走向的固定半固定沙丘，沙丘相对高度10～30 m，样带内沙丘与丘间低地相间而列，地表起伏。样带内降水由南至北梯度递减，区域年平均降水量在100～250 mm［21］。春季积雪融水是土壤水的重要补给［22］，也是荒漠植物的重要水源［23］。区域内发育以梭梭（Haloxylon ammodendron）、白梭梭（H.persicum）等为建群种的典型温带荒漠植被，常见灌木包括淡枝沙拐枣（Calligonum leucocladum）等是该区沙丘呈固定半固定状态的主要因素［24］。灌木以外，区域内草本层植物稀疏分布，常见一年生草本植物如猪毛菜（Solsolaspp.）、角果藜（Ceratocarpus arenarius）［25］。样带南部、中部、北部生物多样性指数呈递减趋势［26］，且荒漠植物物种组成不同，因此沿南北方向对样带进行三等分（图1a），并进一步分析降水与生态弹性时空演变。

图1 研究区位置及概况Fig.1 Location and overview of the study area

1.2 研究方法

1.2.1 生态弹性量化 生态弹性的量化流程如图2所示。文中数据来源于GEE（Google Earth Engine），调用的遥感数据为MODIS MOD13Q1 级NDVI 产品数据，时间分辨率为16 d，空间分辨率为250 m，地理坐标系为WGS84，由于产品数据始于2000 年2月，为保证时间序列数据完整性，故研究时段为2001 年1 月—2020 年12 月。本研究利用最大值合成法对NDVI 数据进行校正，通过提取同一月份不同时间节点上最大的NDVI数值，作为该时期的ND⁃VI 值，从而进一步消除环境干扰［27］，最后得到共20 a NDVI的月时间序列数据。

图2 数据处理流程及残余项与替代稳态间的关系Fig.2 Methods for data processing and relationship between remainder component and alternative stable states

生态弹性作为长期观测得到的稳定性指标，需要较长的时间序列进行量化。该区域植被指数最高值出现在2016年［28］，图2显示2015年前后到2020年NDVI 趋势变化较前期具有较大差异，故将2016—2020年数据划分为一个时间段，并利用该时段数据量化同时期弹性。为比较不同时期弹性的变化，研究中量化所需的时间序列数据需等长，结合地球表层系统分析方法［29］，故以5 a为一个周期进行划分，将数据以5 a 为一个周期，分为t1、t2、t3、t44个时间段，分别代表2001—2005年、2006—2010年、2011—2015 年、2016—2020 年。其中，样带内每一个栅格具体的时间序列为Yt，其中Yt由下列公式中多项构成：

式中：St为季节项，表示NDVI 时间序列随季节变化的函数项；Tt为趋势项，表示NDVI 因季节变动外的压力扰动所产生变化的函数项；Et为残余项，表示NDVI 因脉冲扰动所产生变化的函数项，本研究使用BFAST方法提取以上各项。

将获取的NDVI 数据转化成时间序列对象后，数据处理的主要步骤如图2 所示。首先，提取不同时间段NDVI时间序列，利用BFAST检测变化趋势，其中的季节项使用默认的调和模式来拟合。之后，利用趋势项和残余项进行状态空间建模。该建模有效地反映自然界的真实状态［30］，具体建模方法调用R 语言相应的MARSS 包［31］。完成状态空间建模后，再次使用BFAST 分离趋势项和残余项，并通过残余项划分出稳定状态I（Et>0）和稳定状态II（Et<0）。最终，利用残余项计算退出时间，从而得到生态弹性。

退出时间计算的是环境扰动使生态系统状态退出当前稳定状态所需的平均时间［18-19］。由于状态空间建模还原了数据变化的连续性，这使研究可以准确地计算出系统在不同替代稳态之间转换所用的“退出时间（EXITTIME）”。这里借助统计方法将不同时段的EXITTIMET作为样本数据，然后计算各阶段平均值，方法如下：

1.2.2 生态弹性变化趋势及降水因素相关分析方法 为深入认识弹性与降水波动变化的关系，本研究选取了不同时期多年平均年降水量与t1、t2、t3、t4相应时间段内弹性及其变化趋势做相关分析。降水数据使用的是CHIRPS（Climate Hazards Group Infra⁃Red Precipitation with Station data）日降水量数据，其空间分辨率为0.05°，并通过GEE获取。经分区统计提取生态弹性不同分区平均值，然后利用规则网格（图1a）提取多年平均年降水量数据以及生态弹性数据并计算Spearman 相关系数［32］。生态弹性变化趋势主要通过后期数据减去前期数据得到，如t1到t2时间段弹性变化趋势等于t2时间段弹性减去t1时间段弹性，降水量的变化趋势同理。

2 结果与分析

2.1 生态弹性时空分布

弹性数值及其变化趋势的时空格局如图3。弹性+取值范围在0.075～0.500。t1时段弹性+整体相对较高，t3时段弹性+整体相对较低，可能在该区域不同时间段降水量变化趋势影响了区域弹性。弹性+高值在样带南部分布较为集中，这主要受降水量由南至北梯度递减影响。弹性-取值范围在0.085～0.417。t4时段弹性-数值高于其他时间段，t1到t3时间段高值呈斑块状混杂分布，说明弹性的变化受局部生境条件影响较大。弹性变化趋势显示，t2到t3时间段内样带南部弹性+呈大幅下降趋势，t3到t4时间段内样带整体弹性呈上升趋势，其中样带南部的弹性+存在大幅度上升趋势。不同时期弹性+与弹性-变化趋势均存在较大差异，可能与区域植被对不同类型环境扰动存在适应性差异有关。

图3 不同时间段弹性数值（左）及变化趋势（右）的时空格局Fig.3 Spatial distribution of resilience value and trends in different periods

综合来看，弹性对不同类型环境扰动的响应在空间分布和时间变化上均存在差异，这种差异可能受降水量变化趋势、局部生境条件和植被对环境扰动的适应性共同影响。

2.2 生态弹性与降水量关系的空间差异

样带整体多年平均年降水量在t1到t4时间段内出现大幅度上升。将样带从南到北等分成3个部分（图1a），图4 为样带不同区域各时段的平均年降水量变化，样带南部和中部呈现先升高后下降的趋势，而北部则保持上升趋势。

图4 样带年平均降水量变化Fig.4 Variations in averaged annual precipitation

整体上，弹性与降水量之间具备显著的相关性（表1）。弹性除与本时段内降水相关外，也与前期的降水显著相关。如样带整体在t2时间段的生态弹性+与t1时间段降水量具备较高的相关性，相关系数为0.592。但这种相关性存在局部差异，研究基于降水数据分辨率大小构建矢量网格（图1a），经规则采样各部提取60 个样点的不同时期降水及生态弹性数据平均值，最终计算得出降水与弹性相关系的时空变化。样带内南部区域与降水量保持显著正相关，但样带中部及北部t1到t3时间段弹性与降水量存在显著负相关。此外，不同替代稳态条件下的弹性与降水量相关性差异较大。如样带南部t2时间段的弹性+与不同时期降水量之间的相关性为0.701～0.741，而弹性-与降水量之间的相关性则低于0.34。

表1 弹性与降水量相关性Tab.1 Correlation between resilience and precipitation

样带整体弹性变化趋势与降水量变化趋势的相关性较低（表2），而样带整体弹性变化趋势相比局部与降水量变化趋势之间存在更显著的相关性。综合来看，弹性+和弹性-与降水量相关性随时间和空间变化复杂多样，可能与区域植被对降水变化趋势存在适应性差异有关，并且这可能还与降水季节分布的时空差异有关。

表2 弹性与降水量变化趋势的相关性Tab.2 Correlation between the trends of resilience and precipitation

2.3 生态弹性与降水季节分布的关系

图5a为稳定状态I和稳定状态II在各月份出现的频数，二者的频数差反映了整个样带在不同月份下稳定状态变化的规律。整个生态系统有更大概率在6 月和9—11 月处于稳定状态I，在3—5 月及7月处于稳定状态II（图5b）。

图5 稳定状态频数与降水量变化的关系Fig.5 Stable state frequency in relation to precipitation variability

当降水量变化促进植被生长时（图5c），生态系统在5—6 月的稳定状态I 所停留的时间增加，7 月在稳定状态II 停留的时间减少，甚至存在转换为替代稳态的可能。因此，生态系统在稳定状态I 的平均退出时间增加，弹性+数值降低。反之，弹性-数值降低。

当春夏两季降水量变化对稳定状态的影响有差异时（图5e，图5f），稳定状态I、II的平均退出时间将更加复杂，与不同时期降水量的相关性降低。因此，这可能是替代稳态下的弹性与部分时段的降水量之间相关性不高的主要原因。

3 讨论

本文以阜康北部荒漠植被样带为例，通过结合BFAST、状态空间模型和退出时间法，初步量化并探讨了不同替代稳态下阜康北部荒漠生态系统南缘至腹地生态弹性的时空变化。研究表明，该区域生态弹性总体呈先下降后上升的趋势，不同替代稳态条件下弹性数值从古尔班通古特沙漠边缘至沙漠腹地的空间分布差异显著，且弹性数值对降水波动存在滞后响应，降水波动的空间差异、植物群落的异质性分布及其对降水变化的差异响应可能是造成这种差异的主要成因。

通过NDVI 计算得出的替代稳态下的生态弹性，表明了生态系统功能在承受正面或负面的环境扰动后恢复的时间长度。生态系统恢复时间越长，表明其弹性越小，承受扰动和恢复的能力越低。对荒漠生态系统而言，降水是生态系统生产力的最主要限制因素。根据Bai等［33］对整个中亚地区生态系统弹性研究发现，阜康北部的古尔班通古特荒漠生态系统南缘生态弹性大于沙漠腹地，与样带北部弹性+在t2和t4时间段低于南部的研究结果一致。

本研究结果显示沙漠南缘至腹地降水变化趋势差异明显，样带南部及中部降水量呈现先升高后下降的趋势（图4b，图4c），北部降水量则一直保持上升趋势（图4d），与新疆大部分地区一致［34］。而降水量变化趋势的空间分布差异是不同替代稳态条件下的弹性在t1、t2、t3时间段内空间分布差异的主要原因，结果显示弹性变化趋势与降水量变化趋势间显著相关。如样带整体弹性+在t2到t3时间段的变化趋势与同时间段降水量的相关性为0.418，这说明弹性+会随着同期降水量的上升而上升。

弹性指标反映了系统在历经变化时保持其原有功能、结构、特性和抗干扰及重组的能力。因而，生态系统内物种的丰富度将有助于增加生态弹性［3］。已有研究表明，荒漠植被群落的弹性受物种丰富度的影响显著［35］，而降水量高的区域植物物种丰富度明显大于降水量低的区域［21］。降水由沙漠南缘向腹地的复杂变化导致的植物物种丰度变化［36］，这是阜康北部古尔班通古特沙漠南缘生态弹性高于腹地的自然基础。此外，研究表明立地条件如坡向和坡度也会对生态弹性产生影响［11］。阜康北部荒漠生态系统，沙丘纵横起伏，植物稀疏分布，且主要分布于沙丘底部和丘间地［37］。这是因为地形的起伏导致丘间地有更好的水分条件，使丘间低地具备更高的物种丰富度，进而影响弹性的空间分布，而不同类型沙丘（如固定沙丘、半固定沙丘与流动沙丘）由于水分条件差异也会使得荒漠植被种群的稳定状态在不同地貌类型条件下存在差异［38］。因此，在降水格局和局地立地条件综合作用下，阜康北部荒漠生态系统弹性的空间分布既有南北方向上的趋势性差异，在局部又呈现弹性斑块混杂分布的特征（图3）。

以往研究表明古尔班通古特沙漠降水量与ND⁃VI之间存在滞后响应［39］，本文研究进一步证明这种滞后响应也会反映到基于NDVI计算的生态弹性变化。如样带南部存在弹性+与前期降水量之间的相关性大于0.7。这是因为降水波动对同期植被物种丰富度的影响过程较为缓慢，通过改变植被群落结构从而反映在生态弹性上的时间较长，因此生态弹性对降水量变化存在滞后响应。此外，并非所有的荒漠植物盖度与降水量呈显著正相关［21］，如白梭梭、淡枝沙拐枣等的密度和盖度随降水的增加而减少，而梭梭与降水则没有相关性，这可能与植被在适应温带荒漠气候环境过程中形成了不同生存策略有关［40］。本文研究结果显示该区域存在弹性与降水量呈负相关的现象（表1）。如样带北部在t1到t2时间段弹性与同期降水量变化趋势之间的相关系数在-0.409～-0.511。结合弹性与降水季节分布的关系可知，当样带中部及北部植被的盖度与降水量呈负相关时，弹性随降水量增加而下降。由于降水的季节分布变化多样，且对弹性的影响过程复杂，研究区荒漠生态系统弹性与降水量的相关性也呈复杂多样。综合而言，研究区生态弹性与降水量的相关性取决于植被群落构成、降水的变化趋势及降水的季节分布。

为了解决不同替代稳态条件下的生态弹性量化问题，本研究采用了BFAST变化趋势检测和状态空间建模的方法划分出不同替代稳定状态（稳定状态I 和稳定状态II），并利用退出时间法量化得出生态弹性，通过阜康北部古尔班通古特荒漠生态系统研究的应用，证实了荒漠生态系统普遍存在替代稳态［18］，以及不同替代稳态条件下的生态弹性的时空差异。与已有研究工作［33,35］相比，本文研究通过划分不同替代稳定状态和不同时段研究荒漠生态系统弹性，更好地解释了荒漠生态系统弹性的时空变化规律，为认识荒漠生态稳定性长期动态与环境干扰的关系提供了新的思路。阜康北部的古尔班通古特沙漠稳定性是保障区域生态环境安全、促进区域可持续发展的重要因素［41］，其生态弹性的研究结果有助于深入认识环境干扰下荒漠生态系统功能稳定性的维持机制，对指导荒漠生态系统保护和修复也具有重要的理论指导意义。

4 结论

本研究运用BFAST和状态空间建模的方法，通过对遥感时间序列数据的分期处理提取出了不同时段替代稳态条件下的弹性，并量化不同时段阜康北部的荒漠生态系统弹性的空间分布，然后分区探讨了不同时段生态弹性与降水变化之间的关系。弹性分布总体受降水格局控制，但立地条件导致的植被空间异质性增加了弹性空间分布的复杂性；在时间上，生态弹性在研究时段内呈先下降后上升的趋势，对降水变化的响应存在滞后性，生态弹性与降水量的相关性取决于植被群落构成、植物对降水变化的响应以及降水量变化趋势和季节分布。本文为量化替代稳态条件下的生态弹性提供了方法与技术支撑，在深入认识荒漠生态系统功能稳定性维持机制以及荒漠生态系统保护和修复方面具有理论与实践意义。