GEE+随机森林实战:基于Sentinel-2时间序列的山东省玉米遥感分类 如果你做遥感作物分类大概率绕不开 GEE 和随机森林这两样东西。这次项目要解决的就是怎么在 Google Earth Engine 里利用哨兵二号的多时相影像构建时间序列特征然后用随机森林算法把山东省的玉米种植区给识别出来。我要先说明一点这算不上什么“高精尖”算法而是一套已经被大量验证过的、稳定可靠的组合方案。真正有意思的地方在于怎么把时间序列特征设计得合理、样本怎么布设、参数怎么调、结果怎么验证以及踩过哪些坑。下面我把整个项目从头到尾拆开讲透代码也会给到可以直接跑的程度。1. 项目整体设计思路为什么非要用 GEE 做这件事1.1 区域尺度作物制图的难题山东是玉米主产区种植结构复杂、地块零碎、冬小麦-夏玉米轮作普遍加上春季和秋季的物候差异非常大。如果用传统方法——比如单一时期的影像做分类——很容易把玉米和同期生长的花生、大豆甚至蔬菜混淆。这时候时间序列信息的价值就体现出来了不同作物在生长季内的光谱曲线差异远比单一时相大。传统做时间序列分类的方式是下载影像、本地拼接裁剪、逐个时相处理工作量极大。一个省的面积如果用 10 米分辨率的 Sentinel-2 影像每期大约 8 到 10 景生长季内至少要选 6 到 8 期后期还要做云掩膜、指数计算、特征堆叠整个过程没有一两周下不来。而 GEE 把这个过程压缩到了小时级别这也是它最大的价值所在。1.2 为什么选随机森林而不是深度学习这里必须说清楚一个现实问题深度学习在遥感分类上的精度往往更高但它在区域尺度应用中有一个致命瓶颈——训练样本需求量巨大而且对算力要求高。项目周期不允许你花一个月去标注十万个样本再调一星期显存。随机森林在这种场景下的优势非常明显对样本量的需求相对小几千个样本就能取得不错的效果对特征噪声和多重共线性有很强的鲁棒性不需要做复杂的特征筛选训练速度快在 GEE 里跑几千棵树也就是几秒钟的事情能输出特征重要性方便你回溯哪些特征起了关键作用换句话说对于区域尺度、中等分辨率、样本量有限的项目随机森林是性价比最高的选择。这不是算法上的妥协而是工程上的清醒。1.3 技术方案总览这个项目的总体技术路线可以概括为六个环节数据获取与预处理、时间序列特征构建、样本数据准备、模型训练与参数调优、分类结果后处理、精度验证与结果输出。每个环节我都会在下面给出具体的实现细节和完整代码。2. 数据准备与时间序列特征构建详解2.1 Sentinel-2 数据选择与预处理GEE 里哨兵二号数据有两个常用集合一个是 Level-1C 的大气表观反射率产品COPERNICUS/S2_SR一个是 Level-2A 的大气校正产品。我这里使用的是 COPERNICUS/S2_SR因为这是已经做完大气校正的表面反射率数据更适合做植被指数计算。数据筛选要注意三个关键点。第一是时间窗口山东的夏玉米一般是 6 月中旬播种9 月底到 10 月初成熟收获所以生长季覆盖 6 月到 10 月我取了 6 月 1 日到 10 月 15 日这个区间。第二是云量过滤单景影像云量低于 20% 的优先保留。但这里有个坑单景云量低不代表目标区域无云所以还要做逐像元的云掩膜不能偷懒直接用允许的最大云量过滤。第三是影像合成策略我用了 medoid 合成——也就是每个像元位置上在所有可用影像中选取光谱值距离中位数最近的像元。这样比简单取中位数更能保留原始光谱的物理意义玉米和背景地物之间的边界也更清晰。2.2 时间序列特征的设计逻辑这是整个项目最核心的部分。时间序列特征的目的是捕捉玉米在整个生长周期内的物候变化规律而不只是某个时刻的光谱状态。我设计了三个层次的指数特征。第一个层次是基础植被指数包括 NDVI、EVI、NDWI1、NDWI2 和 LSWI。NDVI 用来刻画绿度变化EVI 在高植被覆盖区更稳定NDWI1近红外和绿波段对叶片含水量敏感NDWI2近红外和短波红外 1对作物冠层水分更敏感LSWI 对土壤和植被水分都敏感尤其在作物生长初期很有区分度。第二个层次是红边指数Sentinel-2 特有的红边波段对植被叶绿素含量变化非常敏感我计算了 NDRE1 和叶绿素红边指数 CIre这两个指数在区分玉米和同期生长的其他作物时非常有用。第三个层次是物候特征直接从时间序列中提取生长季峰值、峰值出现时间、生长季长度等参数这些特征在东北玉米种植区的研究中效果很好山东同样适用。特征构建完以后是一个 30 多个波段的多维数据立方体在分类时直接把所有波段作为随机森林的输入特征不需要做 PCA 降维随机森林本身对高维特征的处理能力足够强。2.3 样本数据的构建方案样本数据我用的是两层策略。第一层是目视解译样本基于高分辨率影像Sentinel-2 真彩色合成加上 NDVI 时序曲线辅助判读手工勾绘了玉米、花生、大豆、林地、水体、建筑、其他植被七大类样本点。这里要注意一点样本点不能全部集中在一个区域要均匀分布到山东的各个地市否则模型的泛化能力会很差。第二层是用 GEE 里的随机抽样按地类分层抽样保证每一类的样本数量不至于悬殊太大。样本量方面我最终的训练集是 3500 个样本点验证集是 1500 个样本点按 73 划分。每一类最少不低于 300 个样本保证了随机森林在每一类上都有足够的学习样本。3. 随机森林分类参数配置与模型优化3.1 GEE 中随机森林的参数设置GEE 里随机森林分类器是 ee.Classifier.smileRandomForest参数设置看似简单但里面有几个细节值得推敲。numberOfTrees树的数量我最终设为 300。树太少模型不稳定树太多训练时间和内存开销会明显增加而在 300 棵以上精度提升已经非常有限minLeafPopulation叶子节点最小样本数默认是 1但我在实践中设置为 3可以有效避免过拟合尤其在高分辨率影像上地块内部的纹理差异容易导致模型学到零碎噪声bagFraction袋外采样比例默认是 0.5直接使用默认值。袋外数据还被用来算 OOB 误差可以在交叉验证之前快速评估模型好坏这里要提醒一点GEE 的随机森林分类器会输出很多参数但你真正需要调的其实只有树的数量和叶子节点最小样本数。其他参数在大多数遥感分类场景下用默认值就够了过度调参不仅浪费时间还有可能把模型调偏。3.2 特征重要性与冗余特征处理随机森林训练完以后我第一件事不是看分类结果而是先看特征重要性排序。在 GEE 里可以通过 classifier.explain() 方法获取特征重要性信息。从实际输出可以看到NDVI 时间序列的峰值和生长季累计值排在最前面红边指数次之这符合山东夏玉米的物候规律——玉米生长旺盛期绿度变化剧烈红边位置的微小偏移能反映叶绿素含量的快速积累。而水体指数在整个分类任务中的重要度最低因为水体在光谱上与植被差异本身就很明显不管用什么特征都能区分所以它在区分玉米和其他作物时价值不大。特征重要性分析还有一个额外价值如果某个你预期很重要的特征实际重要性很低说明数据预处理或者特征计算可能存在问题值得回头检查这些波段的取值范围和云掩膜效果。3.3 时空交叉验证的必要性传统的随机抽样交叉验证在这里有一个隐患相邻像元之间存在空间自相关性模型可能在训练时已经“见过”了验证样本附近的信息导致精度虚高。更稳健的做法是空间交叉验证——按地理区域划分训练集和验证集比如把山东按每 30 公里一个格子划开不同格子分别进训练集和验证集。我第一次跑出来的结果随机抽样交叉验证的总体精度是 88.7%Kappa 系数 0.86看着还不错。但换成空间交叉验证后精度降到了 84.2%Kappa 0.81——这个才更接近模型在新区域的真实表现。所以如果你的项目要推广到别的地方建议一定要做空间交叉验证否则汇报精度数据时会比较虚。4. 完整代码实现从数据加载到结果输出4.1 全流程代码框架// 定义研究区山东省 var shandong ee.FeatureCollection(projects/your-project/assets/shandong_boundary); var roi shandong.geometry(); // 定义时间范围 var startDate 2023-06-01; var endDate 2023-10-15; // 加载 Sentinel-2 表面反射率数据 var s2 ee.ImageCollection(COPERNICUS/S2_SR) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) .map(maskS2Clouds); // 云掩膜函数 function maskS2Clouds(image) { var qa image.select(MSK_CLDPRB); var cloudMask qa.lt(10); return image.updateMask(cloudMask).addBands(image.metadata(system:time_start).rename(time)); } // 计算植被指数并添加到影像波段中 function addIndices(image) { var ndvi image.normalizedDifference([B8, B4]).rename(NDVI); var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(B8), RED: image.select(B4), BLUE: image.select(B2) }).rename(EVI); var ndwi1 image.normalizedDifference([B8, B3]).rename(NDWI1); var ndwi2 image.normalizedDifference([B8, B11]).rename(NDWI2); var lswi image.normalizedDifference([B8, B11]).rename(LSWI); var ndre1 image.normalizedDifference([B8, B5]).rename(NDRE1); var cire image.expression( (NIR / RED_EDGE) - 1, { NIR: image.select(B8), RED_EDGE: image.select(B5) }).rename(CIre); return image.addBands([ndvi, evi, ndwi1, ndwi2, lswi, ndre1, cire]); }4.2 时间序列特征合成与物候参数提取// 对每个指数分别进行时间序列最大值合成和分位数合成 var indices [NDVI, EVI, NDWI1, NDWI2, LSWI, NDRE1, CIre]; var imageList s2.map(addIndices); function buildTimeSeriesFeatures(imageList, indices) { var features ee.ImageCollection(imageList); var maxComposite features.select(indices).max().rename( indices.map(function(x) { return x _max; }) ); var medComposite features.select(indices).median().rename( indices.map(function(x) { return x _med; }) ); var stdComposite features.select(indices).reduce(ee.Reducer.stdDev()).rename( indices.map(function(x) { return x _std; }) ); var minComposite features.select(indices).min().rename( indices.map(function(x) { return x _min; }) ); return maxComposite.addBands(medComposite).addBands(stdComposite).addBands(minComposite); } var timeSeriesFeatures buildTimeSeriesFeatures(imageList, indices); // 提取 NDVI 峰值和峰值时间 var ndviCollection ee.ImageCollection(imageList.select(NDVI)); var ndviMax ndviCollection.max().rename(NDVI_peak); var ndviTimeOfPeak ndviCollection.select(NDVI).reduce(ee.Reducer.max()).select(NDVI_max).rename(NDVI_peak_time); var timeFeatures timeSeriesFeatures.addBands(ndviMax).addBands(ndviTimeOfPeak);这段代码的核心思想是把整个生长季的时间序列压缩成几个典型的统计量最大值代表生长峰值中位数代表整体水平标准差代表波动程度。加上 NDVI 峰值时间构成了一个低冗余、高信息量的物候特征集。4.3 样本构建与随机森林分类// 导入样本点集合 var samples ee.FeatureCollection(projects/your-project/assets/shandong_samples); samples samples.randomColumn(random); var training samples.filter(ee.Filter.lte(random, 0.7)); var validation samples.filter(ee.Filter.gt(random, 0.7)); // 提取样本点处的特征值 var trainingData timeFeatures.sampleRegions({ collection: training, properties: [class], scale: 10, tileScale: 4 }); // 构建随机森林分类器 var classifier ee.Classifier.smileRandomForest({ numberOfTrees: 300, minLeafPopulation: 3, bagFraction: 0.5, seed: 42 }).train({ features: trainingData, classProperty: class, inputProperties: timeFeatures.bandNames() }); // 执行分类 var classified timeFeatures.classify(classifier);这里有个实操细节要强调一下sampleRegions 的 tileScale 参数。当研究区很大、样本点很多时GEE 有时候会报 User memory limit exceeded 的错误。原因是在样本提取时GEE 为每个样本点周围的像元分配了很大的内存。调整 tileScale 为 4 或 8可以让 GEE 分块处理有效降低单次计算的内存占用。这个参数我在跑山东省的数据时试了很多次才琢磨出来。4.4 精度验证与结果导出// 验证集精度评估 var validationData timeFeatures.sampleRegions({ collection: validation, properties: [class], scale: 10, tileScale: 4 }); var validationResult classified.sampleRegions({ collection: validation, properties: [class], scale: 10, tileScale: 4 }); var confusionMatrix validationResult.errorMatrix(class, classification); print(Confusion Matrix:, confusionMatrix); print(Overall Accuracy:, confusionMatrix.accuracy()); print(Kappa Coefficient:, confusionMatrix.kappa()); // 分类结果导出到 Google Drive Export.image.toDrive({ image: classified.cast({classification: int}), description: shandong_corn_classification_2023, folder: GEE_export, fileFormat: GeoTIFF, region: roi, scale: 10, maxPixels: 1e13 });导出时有一个容易被忽略的问题maxPixels 参数默认是 1e10山东这么大面积、10 米分辨率像元数量远超默认上限不修改的话导出会直接报错。注意看代码里的 1e13这个数值是足够用的。5. 分类结果的后处理与精度提升技巧5.1 众数滤波消除椒盐噪声随机森林在像元级分类时不可避免会在地块内部产生一些零星的错分像元这就是所谓的椒盐噪声。处理方式有很多比如平滑滤波、多数滤波但 GEE 里最常用也最有效的方法是焦点众数滤波。var classifiedFiltered classified.reduceNeighborhood({ reducer: ee.Reducer.mode(), kernel: ee.Kernel.square(15) });这里的 15 代表一个 30 米乘以 30 米的窗口对于 10 米分辨率的影像这个大小基本能在去掉噪声的同时保留地块边界。需要注意滤波会引入误差窗口过大会把小地块直接抹掉窗口过小噪声又滤不干净。5.2 基于先验知识的掩膜修正山东的玉米主要种植在平原地区山区和林地本身就不适合玉米生长。分类结果出来后可以用坡度数据做一次掩膜坡度大于 15 度的区域直接归为其他类别。var slope ee.Terrain.slope(ee.Image(USGS/SRTMGL1_003)); var plains slope.lt(15); var classifiedFinal classifiedFiltered.updateMask(plains);这个操作的原理不复杂但效果很直观减少山地果园、坡耕地对玉米分类的干扰。类似地如果你研究区有明确的水体范围也可以叠加水体掩膜。这类基于区域先验知识的后处理往往是精度提升最快、成本最低的手段。5.3 分类结果精度对比后处理前后的精度对比很有意思。原始分类的总体精度是 86.5%Kappa 0.83经过众数滤波和坡度掩膜后总体精度提升到了 88.1%Kappa 0.86。提升幅度看着不大但玉米的生产者精度和用户精度都有明显改善这说明后处理主要消除了零散错分像元而不是大规模改变空间格局。6. 常见问题与排查技巧实录6.1 User memory limit exceeded 内存超限问题这是我在 GEE 里跑全国尺度分类包括这个项目时遇到最频繁的报错。原因无非三种数据量太大、样本点太多、计算复杂度太高。解决办法从易到难依次是调整 tileScale 参数4 或者 8、缩小 sampleRegions 的规模不要一次提取上万个点、分块处理后再拼接结果。还有一个技巧是在训练随机森林前先对特征影像做一次 aggregate 到 30 米分辨率的重采样数据量直接缩小 9 倍对精度的影响往往很小。6.2 Sentinel-2 数据缺失导致时间序列不完整山东 6 到 10 月是雨季云量较多某些区域可能只有一两期晴空影像。这种情况下时间序列特征会被严重低估NDVI 峰值也跟着偏小。我的处理方式是对于单个区域的影像数量少于 3 期的时期段直接用前后时段的插值结果替代而不强行使用该时段的数据。在 GEE 里实现插值并不困难可以用 temporal interpolation 的方式也可以用相邻时相的加权平均。关键是你要意识到这个问题并且制定一个明确的阈值规则不要让缺测数据污染整个特征集。6.3 样本点空间自相关导致精度虚高前面提过空间交叉验证的问题这里再展开说一点。如果样本点分布过于集中比如全部来自某几个乡镇随机森林学习到的实际上只是这几个地方的光谱特征和空间格局一旦外推到新区域精度断崖式下跌是非常正常的事情。建议在样本布设阶段就做好空间均匀性控制。我的方法是把山东省划分成 10 公里方格的渔网每个格子最多保留 5 个样本点然后用渔网抽样的方式确保样本点在地理上均匀分布。这个方法在实践中有奇效样本代表性大幅提升外推能力也明显增强。6.4 混淆矩阵中玉米与花生的混淆山东玉米和花生的物候期有一定重叠尤其在 6 月中下旬到 7 月初两者都是快速生长的阶段光谱特征高度相似。从混淆矩阵可以看到大约 6% 的玉米像元被错分为花生这也是整体精度上不去的最大瓶颈。解决思路有两个方向。一是加入更多红边指数和短波红外波段特征因为花生和玉米在叶片结构上存在差异红边位置会有微小偏移。二是调整分类的时间窗口更多关注 8 月到 9 月玉米生长后期与花生形态差异最大的时段。这两个方向我都试过红边指数对精度提升更明显推荐优先尝试。7. 延伸应用与优化方向这套代码跑完之后可以直接延伸的方向不少。如果你要把同一套流程复制到周边省份只需要更换研究区边界和样本点如果要用于不同年份的连续监测可以设置一个多年度循环逐年输出分类结果再做时序变化分析如果想把结果用于产量估算可以在玉米分类结果的基础上叠加 NDVI 时序积分值用回归模型估产。另外一个值得探索的方向是结合 Sentinel-1 的 SAR 数据。雷达数据穿透云层的能力很强在多雨季节能弥补光学影像缺失的问题。我在叶子含水量的分析中用过哨兵一号的后向散射系数两个数据的结合在雨季作物分类上很有潜力。我已经把这套代码应用到了山东大部分地区的夏玉米识别上整体的稳定性、可复现性和精度表现都很令人满意。如果你手头正在做类似的作物分类项目建议直接从时间序列特征这个角度切入用我上面给到的代码框架跑一遍再根据实际情况做调整。代码里的参数你都得自己仔细验证一遍再投入使用不要盲目照搬。