使用 LandTrendr 绘制最大干扰/增长变化图：完整指南概述本博客介绍了一套基于 Google Earth Engine（GEE）的 LandTrendr 时间分割算法脚本，用于绘制森林干扰（主要是森林砍伐）和植被增长图。脚本能够自动检测每个像元在时间序列中变化幅度最大的植被损失事件，并支持基于参考数据（如国家森林覆盖图）进行参数校准和精度评估。作者：Justin Braaten, Zhiqiang Yang, Robert Kennedy原始代码：https://github.com/eMapR/LT-GEE脚本日期：2018-10-07（后续有更新）1. 准备工作1.1 添加 LT-GEE 公共 API使用本脚本前，需要将 LT-GEE 模块添加到你的 GEE 账户中：访问 https://code.earthengine.google.com/?accept_repo=users/emaprlab/public 并接受。1.2 参数辅助工具推荐使用在线应用辅助确定参数：https://emaprlab.users.earthengine.app/view/lt-gee-change-mapper2. 代码结构与核心功能本脚本整体分为六个部分：部分功能Part 1定义 ROI（研究区）、参考数据集、图像集合构建Part 2LandTrendr 参数设置、运行分割、提取分割段信息Part 3提取最大变化段（干扰）并过滤Part 4基于当地参考图进行参数校准（循环阈值）Part 5生成最终森林砍伐图并可视化Part 6使用 Hansen 全球森林变化数据验证结果3. 核心数据集说明脚本中使用了以下关键参考数据：数据集来源用途NFCGeobosques（秘鲁）2001–2022 年国家森林覆盖图，区分森林/非森林/水体GFCHansen et al. (UMD)2000–2023 年全球森林变化数据，loss 和 lossyear 波段GLC_FCS30D全球土地覆盖2000–2022 年逐年土地覆盖图，用于分析砍伐后的地类变化4. 主要参数说明4.1 时间和空间范围startYear,endYear：1984–2023年度合成时段：01-01至12-31研究区（ROI）：用户自定义矢量或矩形4.2 光谱指数主分割指数：TCW（Tasseled Cap Wetness），也可选NBR、NDMI、NDVI可选拟合指数列表：['NDVI']4.3 掩膜类别maskThese = ['cloud', 'shadow', 'snow', 'water']– 排除云、阴影、雪、水体像元5. LandTrendr 分割参数三组预设脚本提供了三组参数设置，分别对应不同检测目标：参数组maxSegmentsspikeThresholdrecoveryThreshold特点set180.91.0最大化干扰捕获set240.50.25减少误检（FP）set380.90.5减少误检同时提高检出率用户可以根据校准结果选择最佳参数组。6. 关键处理流程详解6.1 年度无云合成函数buildSRcollection：生成年度地表反射率 Medoid 合成函数transformSRcollection：计算光谱指数函数buildLTcollection：生成供 LandTrendr 输入的集合（第一波段为分割指数，其余为拟合指数）6.2 运行 LandTrendr函数runLT：执行时间分割，输出一个三波段数组影像：波段 0：LandTrendr 分段信息（4 行 × N 列）波段 1：拟合值（ftv）波段 2：RMSE6.3 提取分割段信息从LandTrendr数组中提取：每个分割段的起始年份、起始值、结束年份、结束值计算持续时间 (dur)、变化幅度 (mag)、变化速率 (rate)6.4 筛选最大损失段（森林砍伐事件）使用getChangeMap函数，按以下条件过滤：delta: 'loss'– 只关注植被损失sort: 'greatest'– 取变化幅度最大的段year: 2001–2021– 限定时间窗口dur: 4 年– 排除缓慢退化mag: 阈值– 幅度阈值（在校准中确定）7. 校准过程（Part 4 核心）7.1 参考数据准备将 NFC 图重分类为二值掩膜：1= 森林砍伐（2001–2022 年损失）0= 非砍伐（2022 年仍有森林）排除 2001 年前的非森林和水体（设为 null）7.2 幅度百分位数计算计算mag波段在 ROI 内的 0–95 百分位数（步长 5%）得到阈值序列：p0, p5, p10, …, p957.3 校准循环对每个幅度阈值：生成二值砍伐掩膜应用形态学滤波（开运算、闭运算、众数滤波）与参考数据（NFC）构建误差矩阵输出 OA（总体精度）、Kappa、PA（生产者精度）、CA（用户精度）7.4 选择最佳阈值根据校准结果选择平衡误检与漏检的幅度阈值。示例中最终采用mag_thres = 175.1。8. 最终成果输出8.1 砍伐图层次产品产品描述DF_changeImg包含 yod（砍伐年份）、mag、dur、rate 等波段的影像DF_mask_mode二值掩膜（1=砍伐，0=非砍伐，null=掩膜）DF_mask_class分类影像：0=掩膜，1=非砍伐，2=砍伐8.2 砍伐后土地覆盖变化分析使用 GLC_FCS30D 逐年土地覆盖数据提取每个砍伐像元在砍伐年份的土地覆盖类型生成逐年像素计数直方图，展示砍伐后地类组成变化8.3 导出选项所有产品均可导出到 Google Drive 或 GEE Assets，支持：GeoTIFF 格式，30 米分辨率可指定导出区域和金字塔策略9. 验证：对照 Hansen 全球森林变化数据脚本最后部分使用 Hansen 等 (2023) 的全球森林变化数据作为独立参考：提取lossyear波段中 1–21（对应 2001–2021 年）构建与参考样本点的误差矩阵可以导出带有参考值标记的样本集合，用于后续精度评估10. 使用建议先在小区域测试参数，使用辅助 App 可视化不同参数的效果。充分利用校准循环，根据本地参考数据（如高分辨率影像或国家森林图）确定最优幅度阈值。注意计算限制：GEE 对最大像素数有限制，可适当调整maxPixels和tileScale。形态学滤波：对二值掩膜应用focal_mode（众数滤波）可以有效去除椒盐噪声。11. 引用与致谢如果使用本脚本或其中部分功能，请引用：Braaten, J., Yang, Z., Kennedy, R. (2018).LandTrendr Greatest Disturbance/Growth Mapping. GitHub repository: https://github.com/eMapR/LT-GEE以及相关的 LT-GEE API 论文（参见项目主页）。12. 扩展到自己的研究区只需修改以下变量即可应用到其他区域：roi：导入自己的研究区矢量或绘制多边形startYear/endYear：调整时间范围index：选择合适的光谱指数changeParams：根据当地干扰特征调整最大持续时间（dur）和幅度（mag）过滤条件全部代码//########################################################################################################//# #\\//# 使用 LandTrendr 绘制最大干扰/增长变化图 #\\//# #\\//########################################################################################################// 日期: 2018-10-07// 作者: Justin Braaten | jstnbraaten@gmail.com// Zhiqiang Yang | zhiqiang.yang@oregonstate.edu// Robert Kennedy | rkennedy@coas.oregonstate.edu// 参数说明: https://emapr.github.io/LT-GEE/api.html#getchangemap// 项目网站: https://github.com/eMapR/LT-GEE// 注意:// - 运行本脚本前，必须将 LT-GEE API 添加到您的 GEE 账户中。// 访问以下 URL 添加：// https://code.earthengine.google.com/?accept_repo=users/emaprlab/public// - 可使用以下 App 辅助参数化：// https://emaprlab.users.earthengine.app/view/lt-gee-change-mapper// 导入研究区矢量、参考数据和显示参数varROI_1=ee.FeatureCollection("users/mariarisconardizrs/ROI_1"),ROI_2=ee.FeatureCollection("users/mariarisconardizrs/ROI_2"),ROI2_cal=/* color: #081dd6 *//* shown: false */ee.Geometry.Polygon([[[-74.48218398702033,-7.901392279180413],[-74.48218398702033,-8.250827203424928],[-74.08530288350471,-8.250827203424928],[-74.08530288350471,-7.901392279180413]]],null,false),ROI1_val=/* color: #d63000 *//* shown: false */ee.Geometry.Polygon([[[-74.92071088403225,-8.389587067577216],[-74.92071088403225,-8.416758021241236],[-74.90011151879787,-8.416758021241236],[-74.90011151879787,-8.389587067577216]]],null,false),NFC=ee.Image("users/mariarisconardizrs/GT_NFC_Bosque_No_Bosque_2022"),GFC=ee.Image("UMD/hansen/global_forest_change_2023_v1_11"),GT_NFC_DF_VisParam={"opacity":0.51,"bands":["b1"],"min":0,"max":1,"palette":["fdffef","fe3300"]},GT_GFC_LossBand_VisParam={"opacity":1,"bands":["loss"],"min":1,"max":2,"palette":["18b114","ff1602"]},DF_mask_VisParam={"opacity":0.79,"bands":["mag"],"palette":["fffaf8","ff8d0b"]},NPixels_VisParam={"opacity":0.43,"bands":["B1"],"min":0,"max":2,"palette":["ff1e01","ffd005","049900"]},MAG_VisParam={"opacity":1,"bands":["mag"],"min":16.210510261489162,"max":1621.8816873116205,"palette":["5eff07","fffe09","ffad05","ff6501","ff2003"]},DF_class_VisParam={"opacity":1,"bands":["remapped"],"palette":["121212","1b7e02","ff3204"]},DF_yod_VisParam={"opacity":1,"bands":["yod"],"min":2001,"max":2021,"palette":["08ffff","0273ff","f00eff","12ff87","cdff06","f0ff0c","ffde3b","ff8c3b","ffd17e","ff8f76","ff8966","ff4a41","ff6745","ff2929","ab2721","ff2b2b","ff3127","ff1800","ff2704","ff5616"]};// 加载 LandTrendr.js 模块varltgee=require('users/emaprlab/public:Modules/LandTrendr.js');//##########################################################################################// 第一部分：研究区、参考图及影像集合构建//##########################################################################################// 选择研究区（ROI）//var punto = ee.Geometry.Point([-74.51200192757253,-8.084441268366168]);varroi=ROI_1varroiName='ROI1'//** 根据研究区选择更新名称Map.centerObject(roi,10)// 在地图上显示研究区Map.addLayer(roi,{color:'#000000'},roiName);// 参考图：varNFC=NFC.clip(roi);// 国家森林覆盖图：2001-2022 年毁林区域 - Geobosques 数据集//print('National Forest Cover:', NFC); //波段1类别//Map.addLayer(NFC, {}, 'National Forest Cover');varGFC=GFC.clip(roi);// 全球森林变化数据集（GFC），Hansen 等// 波段 'loss'：研究期（2000-2023）森林损失，定义为林地向非林地的替代性干扰。0: 无损失；1: 损失// 波段 'loss year'：森林损失年份，编码为 0（无损失）或 1-23，分别对应 2001-2023 年//print('Global Forest Cover:', GFC);varGLC=ee.ImageCollection("projects/sat-io/open-datasets/GLC-FCS30D/annual").filterBounds(roi);// 全球土地覆盖数据集（GLC_FCS30D）print('Global Land Cover in '+roiName,GLC)// 23 个波段：2000年(b1), 2001年(b2)...2021年,2022年varGLC_1=GLC.first()// 提取影像集合中的第一张影像varGLC_2=ee.Image(GLC.toList(2,0).get(1))// 提取影像集合中的第二张影像// // 在地图上显示 ROI 的全球土地覆盖// Map.addLayer(GLC_1.clip(roi), {}, 'Global Land Cover 1'); //第一张影像// if (roi === ROI_1) {// Map.addLayer(GLC_2.clip(roi), {}, 'Global Land Cover 2')}; //仅对 ROI1 显示第二张影像// 定义影像集合参数：varstartYear=1984;varendYear=2023;varstartDay='01-01';varendDay='12-31';varaoi=roi;varindex='TCW';// 可选：'TCW'、'NBR'、'NDMI'varftvList=['NDVI']// 用于拟合的指数列表varmaskThese=['cloud','shadow','snow','water'];// 需要掩膜的类别// 输出1：每年未掩膜像元计数（用于年度合成）varnClearCollection=ltgee.buildClearPixelCountCollection(startYear,endYear,startDay,endDay,aoi,maskThese);// 每年未掩膜像元的平均值：varMean=function(image){varmean=image.reduceRegion({// 计算每幅影像的平均值reducer:ee.Reducer.mean(),geometry:aoi,scale:30,maxPixels:1e9});returnimage.set('promedio',mean);// 将平均值添加为影像属性};// 对每幅影像应用函数并打印varnClearCollection=nClearCollection.map(Mean);print('Number of clean pixel/year_'+roiName,nClearCollection);// // 可视化每年未掩膜像元数量// var image00 = nClearCollection.filter(ee.Filter.eq('system:index', '16')).first();// print('Image 00 Unmasked Pixels:', image00);// Map.addLayer(image00, NPixels_VisParam, 'Image 2000');// 输出2：年度无云/阴影 Medoid 合成的光谱指数影像集合varannualSRcollection=ltgee.buildSRcollection(startYear,endYear,startDay,endDay,aoi