当前位置：首页 > news >正文

【GEE实战】Sen+MK趋势分析：从代码到地图，解锁植被变化时空密码

news 2026/4/24 22:09:42

1. 从零开始理解Sen+MK趋势分析

第一次接触遥感数据分析时，我被各种专业术语搞得晕头转向。直到真正用Sen斜率估计和Mann-Kendall检验分析植被变化，才发现这套方法就像给地球做"体检报告"——不仅能看出植被"变胖变瘦"（NDVI变化），还能判断这种变化是不是"真的生病"（统计显著性）。在Google Earth Engine（GEE）平台上，我们完全可以用代码实现这套分析流程。

Sen斜率估计本质上是个"老中医把脉"的过程。它通过计算时间序列中所有数据点两两之间的斜率中位数，来判断植被变化的整体趋势。比如2001-2023年间，某个像素点的NDVI值从0.3增长到0.5，Sen斜率会告诉我们这个增长是平缓还是剧烈。而Mann-Kendall检验则是"化验科医生"，它不关心变化幅度，只判断这种变化是否具有统计学意义——就像化验单上的"↑""↓"符号，告诉我们指标异常是否超出正常波动范围。

为什么要用这两个方法搭配？我在分析黄河流域植被变化时就吃过亏。单看Sen斜率，某些区域显示植被明显改善，但MK检验却发现这种变化可能只是随机波动。后来发现，2008年该地区有次异常降雨导致NDVI骤增，单一方法很容易被这种异常值欺骗。二者结合就像"中医+西医"会诊，既看趋势强弱，又看统计可信度。

2. GEE环境准备与数据加载

2.1 初始化GEE工作环境

在GEE代码编辑器（https://code.earthengine.google.com/）新建脚本时，我习惯先做好三件事：

// 1. 定义研究区（以河南省为例） var roi = ee.FeatureCollection("users/your_account/henan_boundary"); // 2. 设置地图中心点和缩放级别 Map.centerObject(roi, 6); // 3. 清除可能存在的旧图层 Map.layers().reset();

遇到过新手常踩的坑是直接使用行政边界名称加载区域，比如ee.FeatureCollection("FAO/GAUL/2015/level2")，这会导致后续计算效率低下。建议提前下载边界文件并上传到个人Assets，速度能快3-5倍。

2.2 加载MODIS NDVI数据

MOD13A1数据集是植被分析的"面包黄油"，但原始数据需要预处理：

// 定义时间范围（2001-2023年） var startYear = 2001; var endYear = 2023; // 创建年度NDVI合成图像集合 var ndviCollection = ee.ImageCollection( ee.List.sequence(startYear, endYear).map(function(year) { var startDate = ee.Date.fromYMD(year, 1, 1); var endDate = ee.Date.fromYMD(year, 12, 31); return ee.ImageCollection('MODIS/006/MOD13A1') .filterDate(startDate, endDate) .select('NDVI') .max() .multiply(0.0001) // 缩放因子转换 .set('year', year) .addBands(ee.Image.constant(year).toFloat().rename('year')); }) ); print('NDVI Collection', ndviCollection);

这里有个关键细节：原始NDVI值需要乘以0.0001转换为实际范围（-1到1）。我曾在西藏植被分析中漏掉这步，结果得出荒谬的Sen斜率值（有的像素点斜率高达300多，明显不合理）。

3. Sen斜率计算实战

3.1 核心算法实现

GEE内置的ee.Reducer.sensSlope()让计算变得简单：

// 计算Sen斜率 var senSlope = ndviCollection.select(['NDVI', 'year']) .reduce(ee.Reducer.sensSlope()) .clip(roi); // 可视化参数 var visParams = { bands: ['slope'], min: -0.005, // 经验值：年际变化通常在这个量级 max: 0.005, palette: ['red', 'white', 'green'] }; Map.addLayer(senSlope.select('slope'), visParams, 'Sen Slope');

参数设置经验：

min/max值需要根据研究区特点调整。在干旱区建议用[-0.01,0.01]，湿润区用[-0.005,0.005]
颜色方案中，红色通常表示退化（负斜率），绿色表示改善（正斜率）

3.2 结果解读技巧

Sen斜率结果包含两个关键波段：

slope：变化趋势率（单位：NDVI/年）
intercept：趋势线截距

我曾帮某林业局分析退耕还林效果，发现虽然整体斜率是正的，但截距差异很大——有些区域初始NDVI就高，同样斜率下实际改善效果不如低值区明显。这时候就需要结合两个指标综合判断。

4. Mann-Kendall显著性检验

4.1 统计原理与代码实现

MK检验的核心是计算Z统计量：

// 创建辅助图像 var ones = ee.Image(1); var zeros = ee.Image(0); var epsilon = 0.001; // 防止除零错误 // 将ImageCollection转为List处理 var imgList = ndviCollection.toList(ndviCollection.size()); // 计算S统计量 var mkS = ee.ImageCollection( ee.List.sequence(0, imgList.size().subtract(2)).map(function(i) { return ee.ImageCollection.fromImages( ee.List.sequence(ee.Number(i).add(1), imgList.size().subtract(1)).map(function(j) { var diff = ee.Image(imgList.get(j)).select('NDVI') .subtract(ee.Image(imgList.get(i)).select('NDVI')); return diff .where(diff.abs().lt(epsilon), zeros) .where(diff.gt(0), ones) .where(diff.lt(0), ones.multiply(-1)) .rename('S'); }) ).sum(); }) ).sum(); // 计算方差 var n = ndviCollection.select('NDVI').count().rename('n'); var varS = n.multiply(n.subtract(1)).multiply(n.multiply(2).add(5)).divide(18); // 计算Z值 var mkZ = mkS .where(mkS.abs().lt(epsilon), 0) .where(mkS.gt(epsilon), mkS.subtract(1).divide(varS.sqrt())) .where(mkS.lt(-epsilon), mkS.add(1).divide(varS.sqrt())) .clip(roi) .rename('Z_score'); Map.addLayer(mkZ, {min: -2.5, max: 2.5, palette: ['blue', 'white', 'red']}, 'MK Z Score');

4.2 显著性判断

在95%置信水平下（α=0.05），Z值的临界点是±1.96：

// 创建显著性掩膜 var significant = mkZ.abs().gte(1.96).rename('significant'); Map.addLayer(significant, {min: 0, max: 1, palette: ['gray', 'yellow']}, 'Significant Areas');

有个实用技巧：将Sen斜率和MK检验结果组合显示：

// 组合显示显著的变化趋势 var trendResult = senSlope.select('slope') .updateMask(significant) .rename('significant_slope'); Map.addLayer(trendResult, { min: -0.005, max: 0.005, palette: ['red', 'white', 'green'] }, 'Significant Trend');

5. 高级可视化与结果导出

5.1 交互式可视化技巧

GEE的ui.Chart接口可以生成时间序列曲线：

// 选择典型像素点分析 var point = ee.Geometry.Point([113.5, 34.5]); var chart = ui.Chart.image.series({ imageCollection: ndviCollection.select('NDVI'), region: point, reducer: ee.Reducer.mean(), scale: 500 }).setOptions({ title: 'NDVI Time Series at Sample Point', vAxis: {title: 'NDVI'}, hAxis: {title: 'Year'}, lineWidth: 2, pointSize: 4 }); print(chart);

5.2 结果导出策略

导出GeoTIFF到Google Drive：

// 导出Sen斜率结果 Export.image.toDrive({ image: senSlope.select('slope'), description: 'Sen_Slope_Export', fileNamePrefix: 'sen_slope', region: roi, scale: 500, maxPixels: 1e13 });

对于大区域分析，建议分块导出。我曾导出全国数据时遇到内存溢出，后来改用Export.map.toCloudStorage分省导出，效率提升明显。

查看全文

http://www.jsqmd.com/news/694678/