
简介这是一份面向遥感与地理信息分析学习者的CCDC连续变化检测与分类实用教程聚焦如何在谷歌地球引擎GEE中完成土地覆盖变化监测全流程。内容从土地覆盖变化背景与监测需求切入系统讲解CCDC算法原理、关键参数、Landsat影像预处理、算法调用及结果制图并附可直接运行的JavaScript代码只需替换研究区即可复现柬埔寨、莫桑比克、哥伦比亚等典型场景。资源包仅1个PDF文件大小约904KB便于快速阅读与本地存档。目前已有269人学习。掌握该教程后读者能独立在GEE中完成时间序列变化检测输出光谱与土地覆盖变化图适用于科研人员、研究生及希望进阶GEE的遥感初学者。1. CCDC把Landsat时间序列拆成可解释的“变化事件”土地覆盖变化影响自然和人为环境全球气候观测系统已将其列为核心气候变量。荒漠化、森林砍伐、城市扩展这些过程在国家到区域尺度上需要用量化方法捕捉而遥感时间序列是唯一可行的数据源。CCDCContinuous Change Detection and Classification正是为这个场景设计的算法它不依赖单一波段也不预设变化方向而是把全部可用Landsat影像逐像元拟合为季节趋势回归模型再用统计检验识别断点。GEE上的实现由Arevalo等人在2020年封装成API支持脚本和图形界面两种方式。配套脚本可以直接作为课程设计或科研预实验的起点。读完这篇你能在GEE里跑通从数据获取、参数设置、变化检测到分类制图的完整流程适合做森林、农业、城市变化监测的遥感工程师和科研人员。2. CCDC算法内核季节回归模型、LASSO拟合与卡方检验2.1 为什么不用年度合成影像年度合成是传统变化检测的主流做法但合成过程会抹掉季节内变化两个年份之间的NDVI差值无法区分物候差异和真实扰动。CCDC换了一条路径把所有云掩膜后的Landsat地表反射率逐像元按时间排序对每个波段拟合一个包含截距、斜率、三组正弦余弦对的回归模型。模型中的季节项能解释旱季和雨季反射率的周期波动趋势项则捕捉多年缓慢变化。当一个像素后续观测持续偏离模型预测才被标记为潜在变化。这类设计的优势在多云地区尤其明显。热带地区每年可用无云影像数量并不固定早期Landsat 5时代一年可能只有两三景有效数据固定合成窗口容易落在不同季节。CCDC用全部观测拟合模型不受固定合成窗口限制时间跨度越长模型对季节基线的估计越稳定。2.2 模型系数的生态学含义回归参数拟合时使用LASSO正则化让不重要的系数收缩为0于是不同土地覆盖类型天然呈现不同的系数组合。拟合完成后三组正弦余弦对被转换为振幅和相位更便于解释。表2-1列出模型参数对应的生态学含义。表2-1 CCDC模型参数与生态学含义模型参数可解释含义Intercept平均反射率水平Slope长期趋势方向与速率Amplitude季节性波动强度Phase物候峰值出现时间RMSE非季节性的观测噪声举个例子GLANCE项目8类训练数据的统计显示森林在NIR波段截距高且振幅大因为树冠物候导致反射率季节性明显城市建成区在SWIR2波段截距高但振幅小因为人工地表反射率年波动弱草本和裸地的RMSE在多数波段高于其他类别因为地表状态本身不稳定。这些系数差异为后续Random Forest分类提供了区分依据。2.3 变化检测的统计机制模型拟合完成后进入监测阶段。GEE实现中监测窗口是一个移动窗口长度等于minObservations参数。窗口内每个新观测与模型预测比较残差计算出的检验统计量服从卡方分布。只有窗口内所有观测的检验统计量都超过chiSquareProbability阈值才判定发生一次光谱变化否则窗口继续向后滑动。// 在GEE中调用CCDC算法的最小脚本 var ccdcResults ee.Algorithms.TemporalSegmentation.Ccdc({ collection: filteredLandsat, breakpointBands: [GREEN, RED, NIR, SWIR1, SWIR2], tmaskBands: [GREEN, SWIR2], minObservations: 4, chiSquareProbability: 0.99, minNumOfYearsScaler: 1.33, dateFormat: 1, lambda: 0.002, maxIterations: 25000 }); print(ccdcResults.bandNames());这里collection接收经过质量筛选和反射率换算的Landsat集合bandNames()会输出结果中包含的波段其中tStart、tEnd、tBreak、numObs分别表示分段起始时间、结束时间、断点时间和观测数。dateFormat: 1表示日期以小数年存储例如2015.62换算成日期时比儒略年更直观。注意minObservations和chiSquareProbability都调大时误检率降低但真实变化被检出也会更滞后。在多云地区建议保持minObservations4因为在云掩膜后连续4个有效观测可能跨越半年时间再增加会显著推迟断点识别。表2-2列出了完整的参数推荐值其中breakpointBands建议使用五个波段tmaskBands只需Green和SWIR2两个波段就能获得稳定的云和云阴影掩膜效果。表2-2 CCDC主要参数推荐值参数名作用推荐值breakpointBands参与变化检测的波段GREEN, RED, NIR, SWIR1, SWIR2tmaskBands云和云阴影掩膜波段GREEN, SWIR2minObservations触发变化所需连续超阈值观测数4chiSquareProbability卡方检验概率阈值0.99minNumOfYearsScaler训练期重新拟合模型所需年数1.33dateFormat日期存储格式1小数年lambdaLASSO正则化系数0.002maxIterationsLASSO最大迭代次数25000lambda取0.002是Zhu和Woodcock原始实验的经验值如果研究区季节信号较弱可以尝试0.01让模型更稀疏减少过拟合。maxIterations在高分辨率大区域上偶尔需要提高到50000才能收敛。2.4 光谱变化不等于土地覆盖变化CCDC分割出来的片段只是光谱轨迹的分段光谱断点可能是干旱导致的植被枯黄、洪水导致的地表覆盖变化也可能是真实的土地覆盖转变。这些变化是否符合项目目标必须通过分类组件来判断。这也是CCDC名称中Classification的由来每个片段用模型系数作为特征配合训练数据训练Random Forest分类器为片段赋予土地覆盖标签再把相邻片段的标签差异映射为变化类型。3. GEE中跑通CCDC柬埔寨国家尺度案例3.1 加载CCDC API并定义参数在GEE JavaScript编辑器中第一步是加载CCDC API。该API由GLANCE项目托管通过require方式引入。// 加载CCDC API var utils require(projects/GLANCE:ccdcUtilities/api); // 定义变化检测参数 var changeDetectionParameters { breakpointBands: [GREEN, RED, NIR, SWIR1, SWIR2], tmaskBands: [GREEN, SWIR2], minObservations: 4, chiSquareProbability: 0.99, minNumOfYearsScaler: 1.33, dateFormat: 2, lambda: 0.002, maxIterations: 25000 }; // 定义研究区柬埔寨全境 var studyRegion ee.FeatureCollection(USDOS/LSIB_SIMPLE/2017) .filterMetadata(country_na, equals, Cambodia) .union();filterMetadata按国家名字段过滤出柬埔寨边界union将可能存在的多个多边形合并为单一几何对象方便后续统一导出。dateFormat: 2在这里用Unix时间存储日期如果你的后续分析在Python或R环境可以改用1以避免时间转换的额外步骤。3.2 获取Landsat影像并执行变化检测var inputParams { start: 2000-01-01, end: 2020-01-01 }; var filteredLandsat utils.Inputs.getLandsat() .filterBounds(studyRegion) .filterDate(inputParams.start, inputParams.end); print(filteredLandsat.size()); changeDetectionParameters[collection] filteredLandsat; var results ee.Algorithms.TemporalSegmentation.Ccdc(changeDetectionParameters); print(results);utils.Inputs.getLandsat()返回经过预处理的Landsat集合内部已经完成Landsat 5/7/8/9串联、pixel_qafMask云掩膜和反射率换算。在2000年到2020年之间柬埔寨大约有7889景有效影像。Ccdc算法输出的是一个Array图像每个像素上存储了不定数量的分段每个分段包含模型系数、起始时间、结束时间和断点标志。表3-1 Landsat输入集合的关键参数参数取值作用start2000-01-01时间序列起点end2020-01-01时间序列终点scale30米输出像元分辨率collectionfilteredLandsat传入Ccdc的影像集合实际导出时scale取30对应Landsat原始分辨率如果最终产品需要5公里格网可以在后续步骤用reduceResolution聚合不必在算法运行阶段改动。3.3 导出结果为AssetArray图像每个像素的段数不同不能直接存成普通栅格必须以Array图像形式导出Asset。导出时pyramidingPolicy必须设置为sample否则金字塔重采样会混合不同分段的系数导致后续分类出错。var paramsCombined ee.Dictionary(changeDetectionParameters) .combine(inputParams) .remove([collection]); Export.image.toAsset({ image: results.setMulti(paramsCombined), scale: 30, description: ccdc_change_results, maxPixels: 1e13, region: studyRegion, assetId: path/to/asset, pyramidingPolicy: { .default: sample } });setMulti把参数写入图像元数据后续分析时可以直接从图像属性中读取避免另存一份参数记录。maxPixels: 1e13适用于整个国家范围如果GEE提示超出像素限制需要把区域拆分成多个格子分别导出。3.4 大区域导出网格拆分与任务循环国家尺度导出很难一次跑完。常见做法是先用makeAutoGrid生成覆盖研究区的规则网格然后用循环提交多个导出任务。var grid utils.Inputs.makeAutoGrid(studyRegion.geometry().bounds().buffer(150000), 2) .filterBounds(studyRegion.geometry()) .toList(100); grid.size().evaluate(function(s) { print(## of grids: , s); for (var i 0; i s; i) { var outGeo ee.Feature(grid.get(i)).geometry() .intersection(studyRegion.geometry()); Export.image.toAsset({ image: results.setMulti(paramsCombined), scale: 30, description: ccdc_change_results, maxPixels: 1e13, region: outGeo, assetId: Cambodia_Change_Results_Grid_ i, pyramidingPolicy: { .default: sample } }); } });makeAutoGrid的第二个参数是网格密度2表示粗略切分格子数量少但每个范围大如果你的研究区跨多个生态区建议用4或6让每个任务更快完成也方便单独排查失败任务。每个格子导出前用intersection裁剪到研究区边界避免导出国界外的无效区域。4. 分类组件与图形界面从莫桑比克到哥伦比亚4.1 分类输入特征与Random Forest训练CCDC的分类过程不是直接对原始光谱分类而是对每个分段保存的模型系数分类。每个分段在每个波段上都有四个特征截距、振幅、相位和RMSE。如表4-1所示森林在NIR波段截距高、振幅明显建筑区在SWIR2波段截距高但振幅小这些差异使Random Forest能够区分不同土地覆盖类型。表4-1 CCDC分段分类特征特征计算方式对分类的作用Intercept回归模型截距区分高反射率与低反射率地表Amplitude正弦余弦幅度区分季节性强的自然地表与季节性弱的人工地表Phase正弦余弦相位区分物候周期不同的植被类型RMSE模型残差标准差辅助识别噪声大的地表裸地、草本实际操作中训练样本要为每个期望类别准备足够多的像素建议每类至少几百个样本点。Random Forest在GEE中是ee.Classifier.smileRandomForest输入特征直接取CCDC系数图像的波段标签来自样本点属性。4.2 使用预训练系数快速出图如果不想自己重新拟合模型可以直接使用全局CCDC系数数据集这些数据是Gorelick等人为解决CCDC初始系数计算瓶颈而预先生产的。在GEE中加载方式如下// 官方发布的Global CCDC系数集合ID以最新文档为准 var globalResults ee.ImageCollection(projects/GLANCE/global_ccdc) .filterBounds(studyRegion);加载后结合少量本地训练样本训练分类器就能得到土地覆盖分类和变化图。这种方式省去了最耗时的变化检测步骤适合那些只需要分类结果、不关心自定义断点参数的业务场景。但要注意全局系数是用固定参数生成的无法反映你研究区特有的季节节律研究区跨多个气候带时建议仍按第3章流程自己运行CCDC。4.3 用Quick-TSTools目视检查断点质量在进入正式批处理前先用Quick-TSTools做一次目视检查。这是Arevalo等人发布的GEE App在地图上点击任意位置即可看到该像素的SWIR1时间序列和CCDC断点。你可以选择不同波段和时间范围检查断点是否与已知变化事件如伐木年份、火灾年份吻合。典型做法是在哥伦比亚亚马逊区域随机选几个森林转牧场的像素看断点是否出现在砍伐年份再选几个常绿林像素确认没有虚假断点。如果断点普遍滞后于真实事件降低chiSquareProbability到0.98重新试算。若频繁出现短时段的碎片段则增大minObservations到6。4.4 GUI流程与脚本流程的取舍莫桑比克案例展示的是纯图形界面操作在GEE App中设置研究区、时间范围点击运行后自动完成变化检测、分类和制图。GUI方式上手快适合需要在项目会议现场演示、或只做一次快速试算的场景。但GUI模式暴露给用户的参数有限LASSO相关参数和分类器结构都不能自定义。我一般先用GUI试算确认参数区间然后在脚本里固定参数批量运行。两者调用的是同一套ee.Algorithms.TemporalSegmentation.Ccdc实现结果可以互相对照差异只在参数暴露自由度上。4.5 生成变化类型图的完整脚本框架在很多项目中最终交付物是“森林到耕地”“森林到城市”这类变化类型图。脚本需要两步先对每个segment做分类再比较相邻segment的类别。一个简洁的做法是用CCDC输出的系数构建分类特征用训练好的Random Forest对每个segment打分并将最大概率类别写回图像。var training samplePoints; // FeatureCollection with class labels var classifier ee.Classifier.smileRandomForest(100).train({ features: training, classProperty: class, inputProperties: [INTP, amplitude, phase, RMSE] }); var classifiedSegments results .select([INTP, amplitude, phase, RMSE]) .classify(classifier, class);这里inputProperties必须与结果数组中的波段名对齐。如果波段名包含波段编号前缀先用rename重命名。分类后每个segment获得类别标签相邻segment标签不一致的位置即变化像素像元值从源类别编码到目标类别编码就是最终的变化类型图。5. 调参经验与常见坑让结果更稳、跑得更快5.1 先验证云掩膜质量CCDC对输入影像质量要求很高。pixel_qa是fMask算法生成的波段在多云地区仍可能残留薄云或云阴影。运行前先用Map.addLayer叠加一景影像的QA波段和RGB合成目视检查掩膜是否覆盖了典型云块。如果云掩膜明显不足可以对集合额外施加光谱阈值过滤再传入Ccdc。5.2 分块导出后不要直接mosaic分块导出的Array图像合并时要特别小心。直接用mosaic()合并会触发金字塔重采样sample策略会在边缘产生混叠。正确做法是先对每一块执行分类或统计产生标量波段后再mosaic或者保留分块结果在后续处理中按地理位置分别读取。5.3 参数调整的优先级我调试时通常按固定顺序调整参数先保持minObservations4、chiSquareProbability0.99跑通样例区若误检过多优先增加chiSquareProbability到0.995若漏检则降低chiSquareProbability然后再看minObservations是否需要减小。lambda一般不改除非研究发现断点出现频率异常。minNumOfYearsScaler不要低于1否则模型无法积累足够的季节观测。5.4 一个实用的断点验证脚本用实地调查数据验证CCDC断点的正确性一直是项目验收的常见卡点。下面这个脚本把每次变化的时间提取为波段叠加到地图上检查局部一致性var breakYear results.select(tBreak) .reduce(ee.Reducer.min()) .clip(studyRegion); var vis { min: 2000, max: 2019, palette: [darkgreen, yellow, red] }; Map.addLayer(breakYear, vis, First break year);tBreak存的是每个断点对应的小数年reduce(ee.Reducer.min())取每个像素第一次断点年份得到首年变化时间图。这张图可以直接对照已知的伐木或火灾年份做分层随机抽样验证评估CCDC在森林损失和造林再生场景下的时间精度。断点时间精度一般在半年到一年之间受影像观测密度影响不需要追求月级精度。本文还有配套的精品资源点击获取