全球生态区临界转变的早期预警信号

《iScience》：Early-warning signals of critical transitions in global ecoregions

　　

摘要
全球约四分之一的生态区显示出结构性转变的早期信号
面临风险的生态区覆盖地球陆地面积的三分之一以上
植被变化与韧性下降在生态区内同时发生
许多受影响生态区位于人类改造程度较低的地区

陆地生态区是受共同环境条件和生物群落塑造的连贯生态单元。随着生态系统逼近临界转变，其韧性下降，表现为对扰动的恢复速度变慢、时间自相关性增强。我们分析了22年的全球卫星植被数据，以检测所有陆地生态区结构性生态变化的早期预警信号。通过将空间聚类、植被趋势分析和韧性指标相结合，我们识别出同时表现出植被变化和韧性下降的生态区。我们发现，约四分之一的生态区——覆盖地球陆地面积的三分之一以上——显示出结构性转变的持续信号。其中许多位于直接人类改造程度较低的地区，这表明大规模环境驱动因素发挥了重要作用。检测到的模式揭示了生态区内部碎片化程度的增加，并突出了可能正在逼近临界转变的地区。该框架提供了一种可扩展的方法，用于识别面临风险的生态系统，并为有针对性的监测和保护规划提供支持。

引言
陆地生态区是景观中具有独特物种组合的不同且基本的区域。由于生态区代表受共同环境条件和生态过程塑造的连贯生态单元，在此尺度上检测早期预警信号有助于了解生态系统层面的稳定性和结构性生态变化，而非孤立的局部变异性。研究表明，生态区通常在环境资源、生态系统和人为影响方面表现出相似性。生态区的划分主要用于全球和地方保护工作，因为研究表明生态区能够有意义地界定所有类群（尤其是植物类群）的生物群落分布。目前已划分了800多个陆地生态区。该分类的初衷有两方面：首先，以优于仅依赖遥感获取的光谱特征或基于温度、降水等环境属性模型的方式来准确代表物种分布；其次，将全球景观划分为适合保护规划的可操作尺度空间单元，生态区的中位面积为56,300平方公里（约240公里量级）。当然，生态区的划分并非完美，它更准确地描述小型物种、低扩散性物种和热带类群的群落。尽管如此，由于生态区代表具有特定生物多样性特征的区域，它们处于优先安排全球保护工作的启发式方法的核心地位，基于生态区的保护规划被证明可作为在物种、群落和生态系统层面保护生物多样性的有效指南。

无论是人为还是自然因素造成的生态压力源，都可能改变生态系统属性并引发扰动，若扰动足够严重，可能导致系统性崩溃。当定义生态系统的一个或多个变量超过一定阈值时，就会导致崩溃，从而引起生态系统属性的根本性转变。生态系统崩溃的可能原因包括气候变化、大规模土地利用变化、生物多样性丧失，以及上述因素的联合作用和反馈。生态系统崩溃可能涉及大规模森林枯死、草原荒漠化，或从有植被状态向稀疏植被状态的转变，导致生物多样性丧失、生态系统功能改变和生态系统服务减少。已记录的实例包括两河地区沼泽因水文改道而大面积退化，以及2000年代初北美西南部松类-杜松死亡事件在大范围旱热期间发生，说明直接人为影响和气候极端事件均可触发快速生态系统崩溃。

系统韧性定义为系统从扰动中恢复的能力。基于此定义，许多努力致力于检测正在逼近崩溃的生态系统。最近有研究表明，超过17,000个树种因人为威胁而面临风险；通过遥感信号分析，亚马逊雨林自2000年代初以来韧性下降；超过60%的陆地植被区域的韧性经历了从增长趋势到下降趋势的转化。此外，全球范围内森林树木从死亡事件中恢复的速率也在下降。由于生态系统崩溃很可能是不可逆的，早期识别对保护工作至关重要，尤其考虑到预测表明许多物种的扰动可能是突发的且几乎同时发生。

在理论层面，生态系统的建模工作是在生态系统在逼近分岔点时会经历临界减速过程的框架下开展的，这会导致信号自相关性增加，代表系统韧性下降。在临近临界点时，系统从扰动中恢复的速度更慢（即"临界减速"）。缓慢的恢复使系统对过去偏差的记忆更长，因此连续观测值变得更加相似（即自相关性增加）。

该方法的有效性最近已通过以下发现得到证明：大多数热带、温带和干旱森林因环境胁迫而出现韧性下降的早期预警信号，且这种下降与森林初级生产力的突然下降在统计上相关。此外，状态空间模型被证明有助于检测某些古湖泊学数据中的制度转变和临界转变。制度转变指生态系统结构和功能的突然且持久的变化，如有机间替代植被状态的转变或大规模生态重组。然而，利用自相关性进行韧性研究需要事先去除信号中的固有趋势，因为这些趋势会影响自相关性值。此外，虽然生态系统的韧性已被证明取决于其空间结构，但许多研究仍在单个格点层面评估韧性，忽略了空间模式和连通性。在这种情况下，大规模解释是从这些孤立点的统计分析中推导出来的，而未明确考虑空间背景或相互作用。

本研究探索生态区内的空间结构。空间结构指生态区内植被动态的空间一致模式，反映潜在的生态梯度、物种组成和扰动过程。识别这种空间一致的子区域，能够检测到分析单个像元无法捕捉的结构性生态碎片化和转变。分析生态区内空间一致的子区域，能够检测到将生态区视为同质单元或单独分析单个像元所会掩盖的结构性生态系统转变。

由于几乎所有环境信号都经历季节性变化，且在当前气候变化条件下许多信号存在趋势，韧性研究忽视了固有信息，而这些信息与生态系统对气候变化的响应密切相关。时间自相关性捕捉了生态系统状态变量在扰动后的恢复动态，提供了独立于系统状态方向性变化的韧性度量。相比之下，植被趋势的方向性反映了生态系统状态本身的结构性生态变化。结合这些互补指标可以区分韧性丧失和结构性生态系统转变，从而提供更全面的生态系统动力学视角。因此，我们建议一种综合方法，利用趋势分析结合韧性和信号以及时空分析来研究生态区单元，因为生态区代表具有独特属性的基本生态区域。我们指出，本研究并非直接检测生态系统崩溃本身，而是识别韧性丧失和结构性生态变化的早期预警信号，这些信号可能先于临界转变出现。

遥感已被证明在区域化方面相比传统生态学方法具有优势指标，特别是植被指数已被证明有助于生态区边界划定以及监测生态区对极端气候事件的响应。此外，遥感提供地球时间连续和空间连续覆盖的能力，以多种方式证明了其提供早期预警信号和检测复杂系统临界点的潜力。归一化差值植被指数是广泛使用的遥感产品，已被证明能有效估算和跟踪冠层叶绿素含量、叶面积覆盖和吸收的光合有效辐射比例，以及总初级生产力和光合通量。此外，NDVI被证明可作为生态区"状态"的代理变量，因为它在全球尺度上与多样性测量显著相关。在各类生物群落中，平均或最小NDVI被用于解释树木和鸟类物种丰富度空间变异的很大一部分，证实了其与支撑生态区概念的生物多样性模式的相关性。此外，基于NDVI的早期预警指标，如滞后1阶自相关性上升，已在亚马逊和其他森林中揭示了韧性下降，证明了其在检测即将到来的临界转变方面的效用。

最近引入了一种NDVI的修正版本，称为核植被指数（kNDVI）。研究表明，对于大多数实际应用，该修正指数可表示为kNDVI = tanh(NDVI2)。kNDVI的输出值限制在0至0.76的范围内，该修正指数在NDVI原本几乎达到饱和的范围内扩展了动态范围。研究表明：（1）在茂密冠层和常绿植被上，kNDVI在NDVI饱和时仍保持敏感性，改善了热带森林或针叶林的物候跟踪。（2）扩展的动态范围使kNDVI更好地捕捉干旱引起的褐化或细微的绿化趋势。（3）它减少了土壤、积雪或大气噪声的影响，产生更干净的时间序列，这对趋势/韧性分析是一个优势。

本研究使用全球、长期（22年）、高空间分辨率的月度kNDVI数据，检测生态区逼近临界转变的可能预警信号。

我们分析的基础在于这样一个理念：生态区是包含多个生态系统的生物地理区域，每个生态系统可能对生态压力的响应不同。这表明识别每个生态区内部的空间-时间聚类，对于检测生态行为的统计显著变化具有重要价值。我们的方法揭示此类聚类，并检查每个聚类内的时间自相关性（TAC），将其作为韧性的共同指标。TAC值增加的聚类（ΔTAC > 0）被解释为韧性下降，该生态区被怀疑正在逼近临界转变。在计算时间自相关性之前，对植被时间序列进行了去趋势处理，以避免与趋势或植被状态突变相关的虚假自相关性增加（详见STAR方法以及补充信息中的图S1–S5）。检测空间一致的聚类至关重要，因为结构性生态系统转变预计以有组织的空间模式出现，而非孤立统计异常。仅比较分布的方法即使变化在空间中随机发生也可能检测到统计差异，而聚类则识别出具有持续生态轨迹的连续子区域。

TAC被计算为时间序列中连续月度kNDVI值之间的Pearson相关系数。为量化韧性变化，分别对时间序列的前期和后期部分计算TAC，TAC的变化定义为ΔTAC = TAC_后期 - TAC_前期。正值ΔTAC表示自相关性增加，对应于从扰动中恢复速度变慢和韧性下降，与既有的早期预警信号理论一致。由于基于自相关性的韧性指标需要去趋势处理，韧性分析和趋势分析通常分别进行。在此，我们同时量化两者——趋势作为方向性状态变化的指标，ΔTAC作为韧性丧失的指标——并将它们的共现解释为转变风险升高。

本研究的目标是通过整合空间聚类、植被趋势分析和基于卫星观测的韧性指标，识别表现出结构性生态转变早期预警信号的陆地生态区。

结果

案例研究

我们通过两个详细案例展示该方法：非洲安哥拉湿润紫茎木稀树草原生态区（图1A–1C和图2）以及亚马逊盆地的马德拉-塔帕若斯湿润森林生态区（图1D–1F和图3）。安哥拉生态区因气候变率、土地利用变化和火灾活动而经历日益增加的生态胁迫，使其成为研究生态系统转变早期预警信号的合适案例。亚马逊生态区则受到大规模森林砍伐、火灾活动和气候驱动生态胁迫的影响，使其对韧性分析具有特殊相关性。生态区的kNDVI趋势通过多变量统计分析（K-Means）进行分析，以发现经历统计显著kNDVI趋势的空间聚类（详见STAR方法部分的详细描述）。

如下例所示，我们方法的创新之处在于检测并量化已划定生态区内部的特定区域，这些区域具有多年统计显著趋势，表现为韧性下降。这些生态区被认为正处于临界转变的早期预警信号之中。安哥拉生态区可能同时面临自然和人为的生态压力，而亚马逊生态区则至少部分受到严重人为压力的影响，其特征包括高频率的森林火灾和道路网络，非法砍伐和退化可能在此扩散。生态区的代表性kNDVI地图（由搭载于AQUA平台的MODIS传感器获取的NDVI值生成）见图1A、1B、1D和1E。相应的kNDVI差值图像（图1C和1F）突出了这些生态区内的时间动态，强调了对它们的生态压力。

图1 非洲西南部安哥拉湿润缅木林地生态区（上行）和亚马逊盆地马代拉-塔帕若斯湿润森林生态区（下行）的kNDVI图像

(A–F) 2003年11月（A）和2024年11月（B）、2003年8月（D）和2024年8月（E）的kNDVI图像及其对应的差值图像（C和F）突出了安哥拉生态区西部和森林生态区西部及南部的显著褐色化。A、B、D和E中的黄色标记标示了生态区边界，A和D中的插图展示了生态区的全球位置。

图2 安哥拉湿润缅木林地生态区的分析（详细说明见正文）

(A) 整个生态区（绿色）和两个聚类（红色和蓝色）的月度kNDVI值揭示了内部聚类之间的信号分离。(B–D) 具有双聚类分离的生态区空间范围（B）。两个聚类的每月十年期kNDVI趋势如图C所示；圆圈表示趋势具有统计显著性的月份。（D）显示了整个生态区的ΔTAC（数值使用高斯滤波器平滑以清晰展示：蓝色和橙红色比例尺分别用于两个聚类）。正的ΔTAC值表示时间自相关增加和韧性下降。

图3 马代拉-塔帕若斯湿润森林生态区，与图2类似

关于安哥拉生态区的主要发现见图2。整个生态区22年的月度平均kNDVI值见图2A（绿线），强调具有预期季节性的年循环。然而，对生态区内所有网格单元月度趋势的K-Means聚类分析（见STAR方法）揭示，生态区的kNDVI特征由两个具有不同时间行为的聚类组成。图B显示了生态区内两个不同聚类的空间范围，揭示蓝色区域占生态区面积的79%。两个聚类对应的时间kNDVI特征如图A所示。聚类代表具有空间一致性的生态子区域，表现出一致的植被动态。这些子区域反映了生态区内正在发生的结构性生态破碎化和重组。识别此类结构模式可为生态系统级别的转变提供见解，而这些转变无法仅通过像素级分析检测到。

值得注意的是，两个聚类的相互kNDVI值在时期开始（2003年）时相对相似，并随时间推移而出现分离。每个月和每个聚类的十年期kNDVI线性趋势如图C所示；圆圈表示趋势具有统计显著性的月份。蓝色聚类在4月至6月的月度kNDVI值出现小幅增加（"变绿"），而红色聚类全年的kNDVI值在下降（"褐色化"），特别是在10月至2月期间，这在11月图像（图1A–1C）中清晰可见。统计分析显示，红色聚类的月度kNDVI趋势在所有12个月中均具有显著性。此外，内部聚类的时间自相关差异（ΔTAC；见STAR方法）如图1D所示：蓝色区域使用蓝色色标，红色区域使用深红色色标。结果表明，总体而言，大部分生态区经历了TAC增加（以蓝色和深红色为主），这表明生态区韧性下降，特别是值得注意的是红色区域具有较高的ΔTAC值（以深红色为主）。具体而言，整个生态区、红色聚类和蓝色聚类的TAC值分别增加了3.3%、11.1%和0.3%。我们的分析表明，如果这些kNDVI趋势持续，预计生态区的植被群体将进一步分离，这可能对广泛的生态特性产生影响，最终可能导致生态区的临界转变。对亚马逊盆地马代拉-塔帕若斯湿润森林生态区的类似分析（图3）揭示，一个占生态区14%的聚类正在经历褐色化（kNDVI值下降），有8个月的kNDVI趋势具有统计显著性，伴随着ΔTAC增加（整个生态区、红色聚类和蓝色聚类分别增加了4.1%、4.9%和3.8%）。

另外四个生态区展示了分析所识别的模式范围（图S7–S14）。在马达加斯加次湿润森林中，36%的生态区在11个月内表现出显著的褐色化趋势；整个生态区的TAC增加了5.9%，红色和蓝色聚类分别增加了4.7%和7.7%（图S7和S8）。东德干湿润落叶林43%的区域在12个月内、德干多刺灌丛森林47%的区域在11个月内、西欧阔叶林35%的区域在7个月内出现显著的变绿趋势。在这些情况下，整个生态区和两个聚类的TAC均增加（图S9–S14）。

全球结果

所呈现的时空分析（图1、2、3和S7–S14）为每个生态区的分析奠定了基础。应用了几项标准以找到似乎正在经历潜在临界转变早期预警信号的生态区，包括仅选择至少在12个月中3个月具有统计显著kNDVI趋势的至少一个聚类的生态区（见STAR方法部分关于筛选生态区标准的描述）。在筛选程序的不同阶段，每个生态区识别出的聚类数量分布见图S6。

分析得出一幅全面的全球地图，显示27%的陆地生态区，覆盖全球陆地总面积的35.3%，经历了kNDVI信号的显著分离，并伴随韧性下降。图4展示了一幅全球地图，显示每个生态区内最小聚类的百分比（图4A和4B），以及在至少一个聚类内观察到统计显著kNDVI趋势的月份数（图4C和4D；从至少3个月开始，低于此值不考虑）。最小聚类的百分比面积代表生态区内结构性分化的程度，而非识别表现出趋势的具体聚类。无论是最小聚类本身还是其互补聚类表现出植被趋势，最小聚类的大小反映了遵循不同生态轨迹的子区域之间空间分离的幅度。该图指出了最易受生态临界转变影响的区域：南美洲（包括亚马逊盆地）、非洲和亚洲的热带地区；中亚和南亚；北美西北部；中欧；以及高纬度寒冷生态区。重要的是要注意，某些标记的生态区包含具有正kNDVI趋势（"变绿"）的聚类，而其他生态区包含具有负kNDVI趋势（"褐色化"）的聚类。当植被增加伴随着韧性下降时，无论是增加还是减少的植被趋势都可能代表结构性生态转变。植被增加可能反映生态重组、生物群系转变或木本植物入侵，而非生态系统稳定化。关键的早期预警信号是方向性变化和韧性丧失的结合，表明偏离了先前的稳定状态。

图4 经历潜在临界转变早期预警信号的生态区全球地图

(A和B) 每个生态区内较小聚类的相对面积（百分比），代表表现出不同植被动态的子区域之间的空间分化程度。(C和D) 在每个生态区内至少一个聚类中检测到的统计显著kNDVI月度趋势数量，表明方向性植被变化的持续性。(E和F) 经历早期预警信号的每个生态区的趋势推导20年kNDVI变化。彩色数值表示根据表现出显著植被变化的聚类的统计显著月度趋势计算得出的20年kNDVI变化。(G–I) 2023年北半球夏季（6月至8月）全球和生物群系kNDVI值。(G和H) kNDVI值全球地图。(I) 六种生物群系中0.05°像素的平均kNDVI，误差条显示像素间标准差。每个像素值是其2023年6月至8月的平均kNDVI：(1) 热带和亚热带湿润阔叶林；(2) 温带阔叶林和混合林；(3) 热带和亚热带草原、稀树草原和灌丛；(4) 温带草原、稀树草原和灌丛；(5) 地中海森林、林地和灌丛；(6) 沙漠和干旱灌丛。

图4E和4F显示了通过将统计显著的月度趋势乘以20年得出的kNDVI变化。因此，彩色数值代表表现出统计显著植被趋势的聚类的趋势推导20年kNDVI变化幅度。

值得注意的是，亚洲和欧洲的大部分生态区预计将出现变绿，而南美洲、热带非洲和北美西北部的生态区主要预计将出现褐色化。图2和3C中呈现的十年期kNDVI趋势以及图4E和4F中预期的kNDVI值变化可能显得很小，但检查全球kNDVI值提供了关于此类趋势可能影响的见解。

图4G和4H描绘了2023年北半球夏季的全球kNDVI值，以及六种生物群系的平均kNDVI值（图4GI）。值得注意的是，kNDVI信号中约0.2的平均差异代表从雨林到稀树草原或从地中海森林到沙漠的转变，强调小的累积趋势可能产生的影响。我们承认植被指数是生态系统功能的间接代理指标，可能无法捕捉生产力或生态过程的所有方面。然而，其全球空间一致性、长期可用性和对植被状态的敏感性使其非常适合检测大规模韧性变化和早期预警信号。

讨论

图4中显示的生态区涵盖了多样的气候和生物地理区域，每个区域都受到不同程度的人为压力。在本研究中，我们旨在区分直接人为影响，即人类活动对环境的局部影响，和由全球变暖引起的全球人为影响。这表明需要更详细的分析来识别对生态压力特别敏感的气候区域和生物群系。此外，了解人为活动对这些生态区的直接影响非常有价值，因为它为人类影响在生态系统中的作用提供了关键见解。图5展示了根据K?ppen-Geiger分类的陆地气候区域（图4A）、生物群系分类（图4C）和人类改造（HM）等级（图4E）。图4B、4D和4F分别展示了从每个气候区域、生物群系和HM等级中经历早期预警信号（如图4所示）的生态区面积百分比。K?ppen-Geiger极地气候分类包括没有植被的广阔区域，如南极洲和永久冰盖覆盖区域。由于早期预警信号只能在有植被存在的地方检测到，气候分类百分比仅从具有有效植被信号的像素计算得出。

因此，极地气候结果仅代表该气候类别中的植被覆盖区域，而非极地气候的全部范围。

图5 按气候区、植被带分类和人类改造程度分类统计显示早期预警信号的生态区

完整说明 | 图片查看器

（A）Koppen-Geiger分类的五个气候区。（B、D和F）还列出了每个气候区（B）、植被带（D）和人类改造程度（F）中经历早期预警信号的生态区面积（以百分比表示）。（C）植被带分类：热带和亚热带湿润阔叶林（1）、干旱阔叶林（2）、针叶林（3）、草原、稀树草原和灌木丛（7）；温带阔叶和混交林（4）、针叶林（5）和草原、稀树草原及灌木丛（8）；北方森林或泰加林（6）；洪水草原和稀树草原（9）；山地草原和灌木丛（10）；冻原（11）；地中海森林、林地和灌木丛（12）；沙漠和旱生灌木丛（13）；红树林（14）。（E）人类改造程度分类。

本分析表明，寒冷和极地气候生态区表现出受影响的植被面积比例最高（约45%）。大约40%的热带陆地面积和约20%的干旱和温带气候生态区面积出现了相同的现象。图4D揭示了受影响的植被带最多：热带和亚热带湿润阔叶林（1）面积的40%–65%；温带针叶林（5）；北方森林（6）；热带和亚热带草原、稀树草原和灌木丛（7）；洪水草原和稀树草原（9）；山地草原和灌木丛（10）；以及红树林（14）等植被带以生态临界转变的早期预警信号为特征。由于人类改造程度（HM）考虑了直接人为威胁（如居住、农业、交通、人类可达性等），我们将具有中等至较高人类改造程度（非常低至低）的生态区归类为在当地（全球范围内）受到人为活动影响。图4F显示，在当地（全球范围内）受影响的生态区中，26%（41%）的面积表现出潜在临界转变的早期预警信号。附加分析表明，本研究中突出显示的生态区总面积的74%属于全球范围内受到人为活动影响的区域。人类改造程度分析使得能够区分在具有直接人为影响的区域发生的韧性丧失和相对完整的生态系统中发生的韧性丧失。在低改造区域中早期预警信号的普遍存在表明，包括气候变化在内的大规模环境驱动因素对观察到的生态系统结构性变化有重要贡献。虽然我们的分析在全球尺度上识别了结构性生态变化的统计早期预警信号，但对特定生态驱动因素的解读需要超出本研究范围的具体区域生态、气候和土地利用分析。我们的结果在全球尺度上识别了统计早期预警信号，但对当地生态机制的解读需要进行详细的区域生态研究。

本研究所呈现的结果令人担忧。四分之一的生态区，覆盖地球陆地面积的35.3%，现在显示出韧性下降和植被方向性变化，而这些生态区覆盖面积的近四分之三属于直接人为威胁较低的地区，表明这些变化很可能由气候变化驱动。热点区域集中在潮湿和季节性干旱的热带地区，其中北方森林、湿润阔叶林、草原-稀树草原-灌木丛、温带针叶林、洪水和山地草原以及红树林中40%–65%的面积符合我们的风险标准。这些生态区显示出内部碎片化增加的特征，表明曾经连贯的生态单元可能面临解体，并发出如果当前趋势持续则可能出现突然生物多样性丧失的风险信号。这些发现提供了一种实用工具，用于识别处于风险中的生态系统，使人们能够在发生不可逆转的生态转变之前进行有针对性的监测、保护优先排序和主动管理。我们的全球早期预警信号图为有针对性的监测和主动保护工作提供了有用的候选清单，以防止不可逆转的临界转变。

研究局限性

本研究在全球生态区尺度上识别了结构性生态变化的统计早期预警信号，但需要考虑若干局限性。植被指数（如kNDVI）是生态系统状态的间接代理指标，可能受冠层结构、背景反射率和物候的影响。时间自相关被解释为韧性指标，但它也可能受外部驱动和随机变异的干扰，因此不能唯一指示临界点的临近。此外，我们的分析并未将变化归因于特定驱动因素，气候因素和人为因素都可能有所贡献。最后，全球和聚类框架强调大规模模式，可能无法捕捉细尺度的生态过程，后者需要区域性调查。

资源可用性

主要联系人

关于进一步信息和资源的请求应发给并将由主要联系人Ilan Koren（[e-mail]）满足。

材料可用性

本研究中使用的全部数据均可从在线来源公开获取。

数据和代码可用性

- 本研究所使用的卫星植被数据是MODIS/Aqua植被指数月度L3全球0.05°CMG版本061产品（MYD13C2.061），可从NASA陆地过程分布式主动存档中心获取（https://doi.org/10.5067/MODIS/MYD13C2.061）。在此分析的kNDVI值是根据该产品的NDVI值按照方法细节中描述的方法计算得出的。陆地生态区和植被带数据集可在https://ecoregions.appspot.com/获取；K?ppen-Geiger气候分类地图可在https://www.gloh2o.org/koppen/获取；全球人类改造数据集可在https://zenodo.org/records/14449495获取。
- 原始分析代码可在https://doi.org/10.34933/ca338033-388c-42b3-8120-bcdfe32ea903公开获取。重新分析本论文所报告数据所需的附加信息可向主要联系人索取。

致谢

本工作得到了以色列高等教育委员会规划和预算委员会（VATAT）的支持，作为"沙漠调谐"气候研究中心的资助项目的一部分。

作者贡献

E.H.和I.K.共同提出了主要想法，进行了分析，并撰写了论文。

利益声明

作者声明不存在竞争利益。

生成式AI和AI辅助技术在稿件准备过程中的声明

在准备本工作期间，作者使用了ChatGPT以提高英语语言、语法和可读性。使用该工具后，作者对内容进行了审查和编辑，并对出版物的内容承担全部责任。

STAR★方法

关键资源表

| 试剂或资源 | 来源 | 标识符 |
|---|---|---|
| 陆地生态区和植被带分类 | Dinerstein等人；世界自然基金会 | https://ecoregions.appspot.com/ |
| MODIS/Aqua植被指数月度L3全球0.05°CMG V061（MYD13C2），2002年7月–2024年8月 | Didan；NASA陆地过程分布式主动存档中心 | https://doi.org/10.5067/MODIS/MYD13C2.061 |
| K?ppen-Geiger气候分类图 | Beck等人 | https://www.gloh2o.org/koppen/ |
| 全球人类改造数据集，1990–2022 | Theobald等人 | https://zenodo.org/records/14449495 |
| 生态区聚类、植被趋势和韧性指标的分析代码 | 本论文 | https://doi.org/10.34933/ca338033-388c-42b3-8120-bcdfe32ea903 |

在新选项卡中打开表格

方法细节

材料

本研究中使用的全部数据均可从在线来源公开获取。

生态区和植被带

全球陆地生态区和植被带分类来自世界自然基金会（WWF）全球生态区数据集（https://ecoregions.appspot.com/）。该数据集根据生态相似性划分陆地面积，包括植被类型、气候和物种组成，广泛用于全球生态和保护分析。每个生态区被分配到一个代表广泛生态类别的植被带，如热带森林、草原和沙漠。在本研究中，植被带分类用于分析早期预警信号在主要生态系统类型中的分布（图5）。

植被数据

MODIS/Aqua植被指数月度3级全球0.05°气候模型网格版本061数据集从NASA Earthdata获取（https://search.earthdata.nasa.gov/search）。该数据集提供高空间分辨率的全球月度植被指数值。我们使用了由此数据衍生的kNDVI产品来分析植被动态和韧性。

气候分类

K?ppen-Geiger气候分类数据集从https://www.gloh2o.org/koppen/获取。该数据集根据长期温度和降水模式将全球陆地区域划分为气候区。气候分类用于检查表现出早期预警信号的生态区在不同气候制度中的分布（图5）。K?ppen-Geiger极地气候类别包括大面积永久冰盖和植被缺失的区域。由于在这些区域无法测量植被指数，在计算气候类别统计量时仅包含具有有效kNDVI信号的像元。

人类改造指数

人类改造（HM）数据集从https://zenodo.org/records/14449495获取。该数据集提供了对陆地人类改造的全球累积测量，纳入了土地利用、基础设施、人类定居点和可达性等因素。HM值范围从极低（人为影响最小）到极高（广泛的人类改造）。在本研究中，HM指数用于评估检测到的韧性丧失在不同程度直接人为影响区域中的范围（图5）。

方法

本研究采用的方法使用了2002年7月至2024年8月（共266个月）期间的"MODIS/Aqua Vegetation Indices Monthly L3 Global 0.05 Deg CMG"（MYD13C2）产品。该经过良好验证的NDVI产品以0.05°的高空间分辨率（对应3600 × 7200矩阵）提供全球月度的地球覆盖。基于NDVI地图，我们按照以下公式制作了月度kNDVI地图：

kNDVI = tanh(NDVI)

其中tanh为双曲正切函数，生成3600 × 7200 × 266的三维月度kNDVI值矩阵。

NDVI和kNDVI等植被指数是生态系统状态和植被结构的综合指标，其时间动态已被广泛用于检测生态系统的韧性丧失和临界点行为。虽然植被指数不是生产力的直接度量，但它们反映了生态系统状态变量，其恢复动态提供了关于系统韧性的关键信息。使用kNDVI具体改善了在高生物量区域的敏感性并减少了NDVI中存在的饱和效应，使其特别适用于在全球尺度上检测长期的结构性植被变化。

趋势分析

三维月度kNDVI矩阵中的每个像元都通过计算一年中每个月kNDVI值的线性趋势，被转换为一个12维向量。图S1显示了一个特定像元的266个kNDVI值（上行）。下方12个面板中的每一个显示了对应该特定像元的月度值以及这些月度值的线性拟合（红色）。红色趋势线的斜率被用于创建该特定检验像元的12值向量。

该过程产生了一个尺寸为3600 × 7200 × 12的三维矩阵。图S2显示了这个矩阵针对非洲安哥兰湿润金合欢林生态区的12个图层（对应图1A、1B、1C和2）。可以观察到显著的空间和时间差异，表明不同的空间区域具有不同的时间行为。

对于每个生态区，使用欧氏度量对12维趋势矩阵进行K-means分析。K-means是一种成熟的多元统计工具，能够找到将p个观测值（像元）按照n个变量（12个趋势）聚类为k个簇的最优解。该分析的第一阶段是估算生态区内的最优聚类数，使用Calinski-Harabasz准则，将最大可能聚类数限制为6。该方法假设生态区内存在簇。为了筛选可能的均匀生态区，我们定义了一个基于时间序列中各簇之间kNDVI差异变化的参数。

对于每个生态区，令Sk,m表示时间序列第一年月份m的第k个簇的平均kNDVI值，Ek,m表示最后一年相应的值。对于包含多个簇的生态区，计算时间序列开始和结束时所有簇对之间簇的平方差：

STi,j = Σ(m=1到12)(Si,m ? Sj,m)2
ENi,j = Σ(m=1到12)(Ei,m ? Ej,m)2

然后计算每对簇在开始和结束时期之间的相对变化：

Di,j = (ENi,j ? STi,j) / STi,j

生态区的筛选指标定义为所有簇对Di,j的平均值。

簇之间未表现出足够结构变化的生态区（定义为?0.5 ≤ Di,j ≤ 1）被归类为均匀型并从进一步分析中排除。对于两个簇的情况，该条件指的是E ≤ 2S或S ≤ 2E。该标准确保只有随时间表现出簇之间有意义的结构差异的生态区被纳入韧性分析。

K-means分析的结果是将生态区空间分解为具有不同时间特征的簇。图S3显示了安哥兰湿润金合欢林生态区分解为2个簇的空间结果。此外，K-means计算每个簇的中心（均值），代表我们分析中kNDVI值的12维月度线性趋势向量。

空间聚类识别出表现出一致植被动态的连贯生态子区域。生态系统结构性转变表现为空间上组织化的变化，而非孤立的像元级波动。因此，聚类使得能够在生态区尺度上检测结构性生态碎片化和重组。

根据均匀性准则未表现出足够类群差异性的生态区被归类为空间同质区域，并从进一步的韧性分析中排除。类群的空间结构也表明了分解结果的稳健性。可以预期，长期的生态过程会产生具有明确空间结构和高连通性的结果，而非类群的随机空间结构。我们的方法论利用了灰度共生矩阵，该矩阵广泛用于纹理分析和土地覆盖分类。我们通过测量标注的生态区类群图中所有像元对的联合概率出现情况，来计算整个生态区内一个像元与其相邻像元的关联程度（基于在4个代表角度：0°、45°、90°、135°处的单次偏移）。空间相关值低于0.15的生态区被过滤掉。一旦确立了空间分解，我们就计算生态区内每个类群的月度线性kNDVI趋势，并评估趋势的统计显著性。图S4的上方面板展示了安哥拉湿润金合欢林地生态区2个类群的月度kNDVI值，下方面板描绘了每个类群的月度值及其线性拟合。每个类群每个月度kNDVI趋势的统计显著性均被评估，采用95%的置信水平。对于每个生态区中的每个类群，计算具有统计显著性的月份数。仅考虑至少包含一个类群且该类群至少有3个月具有统计显著性的生态区。在图S4中，红色类群趋势的统计显著性在所有12个月中均超过95%置信水平，而蓝色类群趋势的统计显著性仅在5月超过95%置信水平。这种统计显著性上的差异并不令人意外，因为蓝色类群的趋势明显较小，需要更低的方差才能具有统计显著性。在趋势分析方法结束时，我们获得了详细列表，列出了具有以下特征的生态区：1. 类群具有稳健的空间结构。2. 初始和最终kNDVI值之间存在足够的变化。3. 至少一个类群中存在统计显著的趋势。

韧性分析

植被指数代表可观测的生态系统状态变量，其恢复动态提供了关于生态系统韧性的信息，并已被广泛用于检测临界点动态和韧性损失。如前所述，我们的数据库包含每个像元的266个月度kNDVI值。与之前的形式化方法类似，我们计算了坐标(i,j)处每个像元两个不同长期kNDVI值的滞后1阶时间自相关(TAC)。TACi,j,1在前133个月度kNDVI值上计算，TACi,j,2在后133个月度kNDVI值上计算。由于TAC计算的目的是量化信号对扰动的韧性，因此有必要从原始信号中减去平均值并去除信号中存在的任何趋势。在计算时间自相关之前，移除了长期确定性趋势，以确保自相关反映的是围绕局部系统状态的恢复动态，而不是被平均植被水平的方向性变化人为放大，这遵循标准的早期预警信号方法论。此去趋势程序还减少了植被状态突变式变化（如由森林砍伐、火灾或土地利用变化等扰动引起）的影响，这些变化可能会人为增加时间自相关。通过移除长期趋势和平均状态，所计算的时间自相关反映了植被信号围绕其局部平衡的恢复动态，而不仅仅是捕捉平均植被水平的变化。因此，本研究中识别出的时间自相关增加主要代表系统韧性的变化，而非突变扰动的直接伪影。

TAC的计算公式为：TAC = corr(X_t, X_{t?1})，其中X_t为去趋势后的月度kNDVI时间序列。为了量化韧性变化，分别计算了早期和晚期时期的TAC，并定义其差值为ΔTAC：ΔTACi,j = TACi,j,2 - TACi,j,1。正的ΔTAC值表示时间自相关增加，对应于系统对扰动的韧性降低。图S5展示了安哥拉湿润金合欢林地生态区的空间ΔTAC图。该图左方面板展示了ΔTAC参数有多大的噪声性。因此，对ΔTAC图应用了1°×1°（20个像元）的空间高斯滤波器（右方面板）。计算了每个生态区中每个类群的平均ΔTAC值，过滤掉了没有任何类群具有正ΔTAC的生态区。时间自相关的增加单独反映了韧性下降，但并不意味着生态系统正在经历结构性生态变化。结构性转变涉及生态系统状态的有方向性重组，并以空间一致的模式发生。通过将韧性损失与有方向性的植被变化及空间聚类相结合，我们的方法识别出正在经历活跃结构性生态转变的生态区，而不仅仅是表现出潜在脆弱性的生态区。

聚类算法使用Calinski-Harabasz准则识别1到6个类群，未施加预设的类群数量。然而，仅保留了满足空间一致性、差异性、植被趋势和韧性标准的类群配置。因此，大多数保留的生态区最优地由2个类群表示。这一结果反映了数据中观察到的生态系统差异化的主导结构模式，而非方法学限制。

图S6展示了在过滤程序各阶段（初始聚类、趋势准则、ΔTAC准则、空间一致性、类群差异性）每个生态区类群数量的分布。虽然聚类最初允许识别多达6个类群，但生态学筛选标准逐步减少了更复杂的类群配置，最终99.13%的保留生态区由2个空间一致的类群表示。这一模式表明，生态区内的结构性分化最常见地表现为两个空间一致子区域之间的分离。

NDVI与kNDVI分析

本研究中的分析基于kNDVI值。然而，由于NDVI的使用更为成熟，我们同样对NDVI重复了分析，并比较了使用NDVI和kNDVI时检索到的具有统计显著性的生态区数量（表S1）。该表显示，具有统计显著性的生态区数量对选择NDVI或kNDVI几乎无关，进一步证明了我们的方法的稳健性。额外一列显示了在每个显著月份阈值下使用NDVI和kNDVI识别出的生态区之间的空间重叠百分比。该百分比计算为两个集合交集的面积除以较小集合的面积，从而代表较小集合中也被另一种植被指数识别出的比例。此定义提供了保守的一致性度量，表明较小检测到的生态区集合中被另一种植被指数确认的比例。

量化与统计分析

所有分析均在MATLAB R2025a（MathWorks）中完成。本研究使用了2002年7月至2024年8月的266个月度MODIS/Aqua观测数据。对于每个有效像元，使用poly1模型分别对每个日历月的kNDVI观测值拟合单独的线性趋势。因此，每个月度拟合包含22个观测值，7月和8月除外，这两月各包含23个。由此产生的每个像元12个趋势斜率用于每个生态区内的空间聚类。

像元根据12个月度趋势斜率之间的欧氏距离通过k均值聚类进行分组。类群数量从1到6的候选值中通过Calinski-Harabasz准则选择。在此分析中，n为生态区内有效像元的数量；类群质心为其组成像元趋势斜率的算术平均值。保留空间一致且具有差异性的类群所采用的标准在方法细节中说明。各筛选阶段保留的类群数量见图S6。

对于每个类群，我们计算了每月的空间平均kNDVI，并使用MATLAB的fit函数及poly1模型对其每个日历月的时间序列分别拟合线性趋势。趋势显著性通过拟合斜率除以其估计标准误来评估。为了近似考虑滞后1阶时间自相关，有效样本量计算为neff = n×(1 ? |r1|)/(1 + |r1|)，其中r1为类群日历月时间序列的滞后1阶自相关。标准误使用neff-2代替普通残差自由度进行计算。当所得t统计量的绝对值超过2.086（20个自由度的双侧95%临界值）时，趋势被归类为显著。此固定阈值应用于所有日历月。当至少一个类群在3个或以上日历月中具有显著趋势时，该生态区满足趋势标准。示例时间序列和拟合趋势见图S1和S4。

为评估韧性变化，将每个像元的266个月kNDVI时间序列分为133个月的两个时期。在分别计算每个时期的滞后1阶时间自相关(TAC)之前，移除了均值和线性趋势。我们将ΔTAC定义为第二时期TAC减去第一时期TAC。ΔTAC图使用方法细节中描述的1°×1°高斯滤波器进行平滑处理。类群级别的ΔTAC值为每个类群内有效像元的算术平均值；每个均值的n为该类群内有效像元的数量。没有任何类群表现出正平均ΔTAC的生态区被排除。ΔTAC图的示例见图S5。

使用NDVI代替kNDVI重复了分析。表S1报告了在显著日历月趋势数量各阈值下识别到的生态区数量以及两个集合之间的空间重叠。重叠定义为交集面积除以较小集合的面积。除另有说明外，本研究报告的空间和类群级别摘要均为算术平均值。统计标准及筛选后剩余的生态区数量在结果、方法细节、图S6和表S1中报告。

补充信息（1）PDF（1.30 MB）
文档S1。图S1–S14和表S1