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

1. 从零开始理解SenMK趋势分析第一次接触遥感数据分析时我被各种专业术语搞得晕头转向。直到真正用Sen斜率估计和Mann-Kendall检验分析植被变化才发现这套方法就像给地球做体检报告——不仅能看出植被变胖变瘦NDVI变化还能判断这种变化是不是真的生病统计显著性。在Google Earth EngineGEE平台上我们完全可以用代码实现这套分析流程。Sen斜率估计本质上是个老中医把脉的过程。它通过计算时间序列中所有数据点两两之间的斜率中位数来判断植被变化的整体趋势。比如2001-2023年间某个像素点的NDVI值从0.3增长到0.5Sen斜率会告诉我们这个增长是平缓还是剧烈。而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.05Z值的临界点是±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分省导出效率提升明显。