GEE中线性回归分析、Sen‘s slope分析、Mann-Kendall检验、FORMA趋势分析一NDVI为例

GEE中线性回归分析、Sen’s slope分析、Mann-Kendall检验、FORMA趋势分析一NDVI为例

1. 线性回归斜率（Linear Regression Slope）

• 含义：通过最小二乘法计算NDVI随时间变化的线性趋势，每个像元的斜率反映植被指数年际变化速率
• 代码实现：ee.Reducer.linearFit()计算线性拟合模型，"scale"分量即为斜率
• 结果解读：
• 正值表示NDVI随时间增长（植被改善）
• 负值表示NDVI随时间减少（植被退化）
• 绝对值大小反映变化速率
• 特点：假设数据呈正态分布，对异常值敏感，适用于线性变化明显的数据

2. Sen’s斜率（Sen’s Slope）

• 含义：基于Theil-Sen Median的非参数趋势估计方法，通过计算所有时间点对的斜率中位数，反映NDVI变化趋势
• 代码实现：ee.Reducer.sensSlope()输出像元级中位数斜率
• 优势：
• 不受异常值影响（如云污染像元）
• 不要求数据服从特定分布
• 适合长时序、噪声较多的遥感数据
• 应用场景：常用于植被变化、气候趋势分析，与Mann-Kendall检验结合使用

3. Kendall相关系数（Kendall’s Tau）

• 含义：通过Mann-Kendall非参数检验评估NDVI时序的单调趋势显著性
• 代码实现：ee.Reducer.kendallsCorrelation()计算tau值
• 结果解读：
• Tau ∈ [-1,1]，绝对值>0.5表示强相关
• 结合p值判断显著性（常用置信度1.65/1.96/2.58对应90%/95%/99%）
• 特点：
• 不要求线性关系
• 能检测非线性单调趋势
• 适合小样本数据分析

4. Forma趋势分析

• 含义：基于时间序列突变检测的算法，识别NDVI序列中的持续性趋势和转折点
• 代码实现：.formaTrend()方法（GEE内置函数）
• 应用特性：
• 可检测渐变和突变过程
• 擅长识别森林砍伐/恢复等突变事件
• 输出趋势强度等级分类结果
• 优势：综合时间序列特征，比单一统计量更直观展示植被动态

---

方法对比表

方法	参数类型	异常值敏感性	显著性检验	适用场景
线性回归	参数	高	需要	线性趋势明显的数据
Sen’s斜率	非参数	低	需结合MK	长时序/噪声数据
Kendall相关系数	非参数	低	自带	单调趋势检验
Forma趋势	非参数	中	内置	突变检测/复杂趋势分解

这些方法在植被动态研究中常组合使用：Sen’s斜率量化变化幅度，Kendall检验判断显著性，线性回归辅助验证线性假设，Forma分析突变特征。建议重点关注Sen’s斜率与Kendall检验的组合结果，这是当前遥感趋势分析的主流方法。

代码实现

利用Modis

/*
本代码演示了如何利用 MODIS 数据计算植被趋势和变化指标，
包括线性回归斜率（linear regression slope）、Sen's slope、Kendall相关系数以及 Forma 趋势等。
*/

// ==================================================================
// 1. 定义研究区和时间范围
// ==================================================================

// 定义县级行政区划数据源
var county = ee.FeatureCollection("projects/ee-xyt556/assets/Data/county");
// 设置筛选条件：地级市和县级行政区名称
var cityName = "徐州市";
var countyName = "云龙区";
// 创建复合过滤器，同时满足地级和县级条件
var filter = ee.Filter.and(
  ee.Filter.eq("地级", cityName),
  ee.Filter.eq("县级", countyName)
);
// 应用过滤器，获取感兴趣区域（ROI）
var roi = county.filter(filter);

var time_start = "2001-01-01",
    time_end = "2025-01-01";

// 可视化研究区 ROI
Map.centerObject(roi, 10); // 将地图中心设置为 ROI，缩放到适当级别
Map.addLayer(roi, {color: 'red'}, '研究区边界', false); // 添加 ROI 边界图层，显示为红色

// 将地图中心定位到指定区域
Map.centerObject(roi);



// ==================================================================
// 2. 加载 MODIS NDVI 数据集
// ==================================================================

// 从 MODIS/061/MOD13A2 数据集中选择 NDVI 波段，时间范围为2001年至2021年
var modis = ee.ImageCollection("MODIS/061/MOD13A2")
             .select('NDVI')
             .filterDate('2001', '2021');

// ==================================================================
// 3. 为每幅影像增加时间信息
// ==================================================================

/**
 * ndvi_time 函数
 * 提取影像的 system:time_start 属性，转换为秒并与影像合并，生成新的波段
 * 同时保留原始所有属性
 *
 * @param {ee.Image} img - 输入影像
 * @return {ee.Image} 带有时间信息的新影像
 */
function ndvi_time(img) {
  // 提取 system:time_start 属性并除以 1e9 转为秒数
  var time = img.metadata('system:time_start').divide(1e9);
  // 将时间信息作为一个波段添加到影像中，同时复制所有原始属性
  return img.addBands(time)
            .copyProperties(img, img.propertyNames());
}

// 对 MODIS 影像集合应用 ndvi_time 函数，构建带时间信息的时序集合
var modis_time = modis.map(ndvi_time);

// ==================================================================
// 4. 计算线性回归斜率（Linear Regression Slope）
// ==================================================================

// 选择时间（system:time_start）和 NDVI 波段，利用线性回归求解
var linear_reg = modis_time.select('system:time_start', 'NDVI')
                   .reduce(ee.Reducer.linearFit())
                   .select('scale');  // 'scale' 对应斜率

// 将线性回归斜率图层裁剪到研究区，并添加至地图（调色板： red, black, green）
Map.addLayer(linear_reg.clip(roi),
             {palette: ['red', 'black', 'green']},
             'linear_reg',
             false);

// ==================================================================
// 5. 计算 Sen's 斜率（Sen's Slope）
// ==================================================================

// 同样选择时间和 NDVI 波段，利用 ee.Reducer.sensSlope() 计算 Sen's 斜率
var sen_slope = modis_time.select('system:time_start', 'NDVI')
                   .reduce(ee.Reducer.sensSlope())
                   .select('slope');

// 将 Sen's 斜率图层裁剪到研究区，并添加至地图
Map.addLayer(sen_slope.clip(roi),
             {palette: ['red', 'black', 'green']},
             'sens_slope',
             false);

// 同时将 Sen's 斜率结果导出到 Google Drive
Export.image.toDrive({
  image: sen_slope.clip(roi),
  description: 'sen_slope',
  scale: 1000,
  region: roi,
  maxPixels: 1e13,
  crs: 'EPSG:4326'
});

// ==================================================================
// 6. 计算 Kendall's 相关系数（Kendall's Tau）
// ==================================================================

// 对 NDVI 波段应用 kendallsCorrelation reducer，计算 Kendall's Tau 值
var mannkendall = modis_time.select('NDVI')
                   .reduce(ee.Reducer.kendallsCorrelation());

// 将结果中的 'NDVI_tau' 波段裁剪到研究区，并添加至地图以展示（调色板同样设置为 red, black, green）
Map.addLayer(mannkendall.select('NDVI_tau').clip(roi),
             {palette: ['red', 'black', 'green']},
             'mannkendall',
             false);

// ==================================================================
// 7. 计算 Forma 趋势
// ==================================================================

// 利用 formaTrend 方法（可能是自定义扩展）计算 NDVI 的趋势
var forma = modis_time.select('NDVI')
            .formaTrend();

// 将 Forma 趋势图层裁剪到研究区，并添加至地图
Map.addLayer(forma.clip(roi), [], 'forma', false);

/*
研究区定义
MODIS 数据加载 从数据集 "MODIS/061/MOD13A2" 中提取 NDVI 波段，设置时间范围为 2001 至 2021 年。 MODIS 数据通常经过预处理并标度处理，此处后续计算中也需要注意数值范围（例如在时间函数中未另外转换）。
增加时间信息 使用 ndvi_time 函数从每幅影像中提取 system:time_start 属性，并将其转换为秒后加入影像作为新波段（便于后续时序分析），同时复制所有原始影像属性。
线性回归斜率计算 选择 system:time_start 与 NDVI 波段，使用 ee.Reducer.linearFit() 求解线性回归，返回的结果中 scale 波段即为斜率；添加至地图时使用了红色、黑色和绿色调色板表示不同斜率值范围。
Sen's 斜率计算 利用 ee.Reducer.sensSlope() 计算 Sen’s 斜率，它常被用于非参数趋势分析，对异常值不敏感。 结果同样通过调色板展示，并通过 Export.image.toDrive 导出。
Kendall’s 相关系数计算 使用 ee.Reducer.kendallsCorrelation() 对 NDVI 数据进行计算，获得 Kendall Tau 值，即衡量趋势的统计显著性。 输出的图层中选择 'NDVI_tau' 波段进行展示。
Forma 趋势 调用 formaTrend() 方法计算 NDVI 趋势，由于此方法不属于 Earth Engine 内置标准方法，一般为外部扩展或自定义方法，此处假设其已定义并可使用；展示在地图上。
*/

####利用Landsat

// ==================================================================
// 1. 定义研究区和时间范围
// ==================================================================

var county = ee.FeatureCollection("projects/ee-xyt556/assets/Data/county");
var cityName = "徐州市";
var countyName = "云龙区";
var filter = ee.Filter.and(
  ee.Filter.eq("地级", cityName),
  ee.Filter.eq("县级", countyName)
);
var roi = county.filter(filter);

var time_start = "2015-01-01",
    time_end = "2024-12-31";

Map.centerObject(roi, 10);
Map.addLayer(roi, {color: 'red'}, '研究区边界', false);

// ==================================================================
// 2. 加载 Landsat 8 数据集
// ==================================================================

var landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
                .filterBounds(roi)
                .filterDate(time_start, time_end)
                .filter(ee.Filter.lte('CLOUD_COVER', 20));

// 定义 NDVI 计算函数
function calculateNDVI(image) {
  var red = image.select('SR_B4').multiply(0.0000275).add(-0.2);
  var nir = image.select('SR_B5').multiply(0.0000275).add(-0.2);
  var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
  return ndvi.copyProperties(image, image.propertyNames());
}

var landsat_ndvi = landsat8.map(calculateNDVI);
print("NDVI:\n",landsat_ndvi);

// ==================================================================
// 3. 为每幅影像增加时间信息
// ==================================================================

function ndvi_time(img) {
  // 将 system:time_start 转换为秒数（ee.Number）
  var time = ee.Number(img.get('system:time_start')).divide(1000);
  // 设置属性供图表使用
  img = img.set('system_time', time);
  // 获取 NDVI 波段的投影信息，用于重新投影时间波段
  var ndviProj = img.select('NDVI').projection();
  // 创建常数时间波段，并强制转换为 Float 类型，再使用 NDVI 的投影
  var timeBand = ee.Image.constant(time)
                    .rename('system_time')
                    .toFloat()
                    .reproject(ndviProj);
  return img.addBands(timeBand);
}

var landsat_time = landsat_ndvi.map(ndvi_time);

// ==================================================================
// 4. 计算线性回归斜率（Linear Regression Slope）
// ==================================================================

var linear_reg = landsat_time.select('system_time', 'NDVI')
                    .reduce(ee.Reducer.linearFit())
                    .select('scale');
Map.addLayer(linear_reg.clip(roi),
             {palette: ['red', 'black', 'green']},
             'linear_reg',
             false);

// ==================================================================
// 5. 计算 Sen's 斜率（Sen's Slope）
// ==================================================================

var sen_slope = landsat_time.select('system_time', 'NDVI')
                    .reduce(ee.Reducer.sensSlope())
                    .select('slope');
Map.addLayer(sen_slope.clip(roi),
             {palette: ['red', 'black', 'green']},
             'sens_slope',
             false);

Export.image.toDrive({
  image: sen_slope.clip(roi),
  description: 'sen_slope',
  scale: 30,
  region: roi,
  maxPixels: 1e13,
  crs: 'EPSG:4326'
});

// ==================================================================
// 6. 计算 Kendall's 相关系数（Kendall's Tau）
// ==================================================================

var mannkendall = landsat_time.select('NDVI')
                    .reduce(ee.Reducer.kendallsCorrelation());
Map.addLayer(mannkendall.select('NDVI_tau').clip(roi),
             {palette: ['red', 'black', 'green']},
             'mannkendall',
             false);

// ==================================================================
// 7. 计算 Forma 趋势
// ==================================================================

var forma = landsat_time.select('NDVI').formaTrend();
print(forma);
Map.addLayer(forma.clip(roi), [], 'forma', false);

// ==================================================================
// 8. 输出时序图表
// ==================================================================

print(ui.Chart.image.series(landsat_ndvi, roi, ee.Reducer.mean(), 1000, 'system:time_start')
      .setOptions({title: 'Landsat NDVI 时序图'}));

研究方法概述

NDVI（归一化植被指数）时间序列趋势分析是一种用于揭示植被生长状态随时间变化规律的重要方法。整个过程包括数据预处理、时间序列构建、趋势计算（包括各种统计方法）、显著性检验以及结果可视化。下面将详细介绍这一过程与方法。

---

一、数据预处理

1. 数据获取

• 数据源选择：常用的数据集有 Landsat、MODIS、Sentinel 等。这些数据集经过辐射和大气校正处理，适合直接进行 NDVI 计算。
• 研究区提取：首先确定一个兴趣区域（ROI）或行政区划边界，从数据集中只筛选出与该区域有关的影像。

2. NDVI计算

• 计算公式：
  $$\text{NDVI}=\frac{(\text{NIR}-\text{Red})}{(\text{NIR}+\text{Red})}$$
  Landsat 中通常用 SR_B5（近红外）和 SR_B4（红光）波段计算 NDVI。
• 比例因子和偏移量：有些数据（例如 Landsat Collection 2 Level-2）需要应用比例系数和偏移值（例如程序中的 SCALE_FACTOR 和 OFFSET）。
• 云遮挡与噪音去除：进行云和阴影掩膜处理，保证每个 NDVI 值的可靠性。

---

二、构建NDVI时间序列

1. 影像集合构建

• 将研究区内所有满足时间与质量要求的 NDVI 影像构建成一个影像集合（ImageCollection）。

2. 时间聚合（选取月均或年均）

• 目的：由于单幅影像可能受云、季节性波动等因素影响，直接分析可能存在噪声。
• 方法：
  • 按照时间（如月份或年份）对影像进行分组
  • 对每组影像计算统计平均值（例如均值），生成月均或年均 NDVI 影像
• 在程序中，构建年均 NDVI 集合就是一个典型示例：
  • 为每图像标注采集年份
  • 对同一年中的所有影像进行平均，从而得到一幅代表该年度的 NDVI 影像。
  • 同时将每幅影像的时间属性设为当年1月1日（或直接用年份），使后续趋势分析时的时间变量（通常叫做 decimal_year）清晰表达时间信息。

---

三、趋势分析方法

在得到 NDVI 时间序列后，我们希望了解植被状态是否在持续变化、变化的方向和速率是多少。常用的趋势分析方法包括：

1. 线性回归

• 方法：把时间作为自变量（如 (t)），NDVI均值作为因变量，拟合线性模型：

  [
  \text{NDVI} = a + b \times t
  ]

  回归系数 (b) 即表示 NDVI 单位时间内的变化率。

• 优点：简单直接，易于理解。

• 缺点：对异常值（噪声）较为敏感，要求数据满足线性关系。

2. Sen’s 斜率估计

• 基本原理：对所有可能的观测对计算斜率，然后取中位数作为整体斜率。
• 优点：是一种非参数方法，对异常值鲁棒，适合于噪声较多的遥感数据。
• 解释：该斜率表示 NDVI 随时间变化的稳健估计值。

3. Mann-Kendall检验

• 基本原理：一种非参数统计检验方法，根据各观测值之间的秩次关系（即成对比较一致或不一致的次数）来计算 Kendall tau 相关系数。
• 含义：
  • tau 值范围在 -1 到 1 之间
  • 正值表明 NDVI 随时间整体上升，负值则表明下降
  • 此方法还能检验趋势的统计显著性，帮助判断观察到的变化是否为随机波动或确实存在长周期变化。

4. 其他方法

• 如 LOESS 局部加权回归、FORMA趋势分析等，适用于捕捉非线性变化特征。

---

四、统计检验与结果解释

在趋势分析中：

• 斜率 (b)：表示 NDVI 每年（或每月）变化量，正斜率表示植被在改善，负斜率表示植被在退化。
• Mann-Kendall 检验：通过计算 Kendall tau 值和对应的 p-value，我们可以判断时间序列中的趋势是否显著。
• 多方法结果对比：通常综合使用线性回归、Sen’s 斜率与 Mann-Kendall 检验，可以相互验证，增强结论的可信度。

---

五、结果可视化与应用

1. 时序图

• 绘制 NDVI 均值随时间变化的图表，横轴为时间（年或月），纵轴为 NDVI 平均值。
• 可以直观展示植被变化趋势及季节性/年际波动。

2. 趋势地图

• 根据回归或 Sen’s 斜率结果，生成地图图层，利用颜色表达 NDVI 趋势的方向和幅度：
  • 正趋势区域显示为绿色，负趋势为红色，变化幅度由色调深浅表示。

3. 应用领域

• 植被监测和生态系统健康评估
• 土地退化和荒漠化监控
• 农业与水资源管理

---

总结

NDVI时间序列趋势分析的核心步骤如下：

1. 数据预处理：收集多时相遥感影像、校正和计算NDVI、筛选研究区影像。
2. 构建时间序列：通过空间统计方法（如年均或月均）生成连续的NDVI时间序列，并附加时间属性。
3. 趋势计算：用线性回归、Sen’s 斜率及 Mann-Kendall 检验等统计方法，定量描述 NDVI 随时间的趋势和变化速率。
4. 显著性检验与解释：评估趋势的统计显著性，进而判断实际植被变化是否可靠。
5. 结果可视化：通过时序图和趋势图等直观图形展示 NDVI 的长期变化，为环境监测与管理提供数据支持。

这种趋势分析方法能帮助研究者捕捉到长期、稳定的植被动态变化信息，即使数据中存在噪声或异常值，也能较为准确地反映环境变化的真实趋势。