简介:面向GIS与遥感学习者,这份资源聚焦“基于多时相哨兵影像和随机森林算法提取冬小麦”这一实践任务。文档以GeoScene Pro为操作平台,系统记录了哨兵多光谱影像的利用方式、训练重分类器的完整流程,以及训练精度达到0.9995、变量重要性分布等关键输出;同时呈现了生成分类数据集时出现的ERROR 999999异常、并行实例分配及5分20秒运行失败的排错过程,便于读者复现实验并避开同类问题。压缩包内仅1个docx文档,大小1.61MB,内容紧凑、可直接阅读,适合需要快速掌握遥感机器学习分类方法、或正在摸索GeoScene Pro中随机森林工具链的初学者与进阶用户。目前已有203人学习下载,可作为农业监测、作物估产和耕地信息提取类项目的入门参考。对于希望结合多时相影像与集成学习算法解决实际作物制图需求的人来说,这份操作记录能帮助节省环境调试时间。
1. 多时相哨兵影像和随机森林:冬小麦提取为什么绕不开这两样
农业保险定损、种植补贴核验、粮食产量估算,第一件事都是做出准确的冬小麦分布图。用多时相哨兵影像和随机森林算法做冬小麦提取,是当前农业遥感里验证最充分、落到业务也最稳的一条路线。早年用单时相影像分类,冬小麦和油菜、大蒜这些冬季作物光谱长得太像,误分率能到三成以上;后来大家转向多时相思路,抓住冬小麦“秋季播种、冬季休眠、春季返青、初夏收割”的物候节律,配合随机森林的时序特征学习能力,精度才有质的提升。这篇笔记准备把一个在 Gis online 平台上可复现的完整流程讲透:数据窗口怎么选、特征怎么构、随机森林参数怎么设、最容易在哪些地方翻车,以及第二年怎么低成本复用模型。新手能照着跑通第一版,熟手可以拿我这套流程对照自己线上项目的边界。
2. 在 Gis online 上准备哨兵影像:数据选取、时间窗口与云掩膜
2.1 多时相数据怎么选:冬小麦物候窗口的四个关键期
冬小麦从播种到收割跨两个日历年,这是它和春玉米、大豆最大的区别,也是多时相方法能奏效的根本原因。我一般从当年 10 月取到次年 6 月,按物候段划分时间窗口,每个窗口至少保留一景无云影像。下面这个表是我在北方冬小麦主产区常用的窗口划分,直接抄作业问题不大:
| 物候段 | 大致时间窗口 | NDVI 变化特征 |
|---|---|---|
| 出苗—越冬 | 10 月—次年 2 月 | 低位缓降,越冬期呈平台 |
| 返青—拔节 | 3 月—4 月上旬 | NDVI 快速抬升 |
| 抽穗—灌浆 | 4 月下旬—5 月下旬 | NDVI 达到全年峰值 |
| 成熟—收割 | 6 月上旬—中旬 | NDVI 快速回落 |
这四段里,返青拔节和抽穗灌浆两个窗口最重要,它俩决定了冬小麦的“峰值高度”和“峰值出现时间”,这两个特征在随机森林里往往贡献度最高。选数据时每个物候段至少保留一景无云影像,总影像数控制在 6 到 10 景。太少,物候曲线不完整;太多,特征维度膨胀,训练变慢还容易引入噪声。技巧是这样的:先用平台自带的日期排序看每个月的可用影像数,再手动挑云量低的景,别一股脑把整年影像全塞进集合。
提示:多时相不是越多越好。相邻两景影像如果间隔小于 15 天,信息高度冗余,对随机森林区分冬小麦几乎没有新增贡献,反而拖慢训练和导出。
2.2 云量过滤与像素级云掩膜:让影像集合干净可用
哨兵二号影像的云量标识分影像级和像素级两层。影像级看 CLOUDY_PIXEL_PERCENTAGE,像素级看场景分类波段(SCL)。两层都要做,只过滤影像级云量远远不够——一朵云只要遮了研究区的 30%,整景影像就会被过滤掉,你可能把关键的返青期影像全丢了。反过来,云量 10% 以下的影像,研究区边缘也可能挂着碎云,不做像素级掩膜,这些云像元之后会进入特征提取,直接拉低 NDVI,把冬小麦误判成其他地类。
我在 Gis online 上处理冬小麦影像集合的标准流程是这样,下面代码以主流在线遥感平台脚本编辑器的 JavaScript API 为例,不同平台命名空间略有差异,逻辑一致:
// 第一步:按研究区、时间窗、影像级云量粗筛 var roi = table; // 地图上绘制的冬小麦种植区边界 var s2 = ee.ImageCollection("COPERNICUS/S2_SR") .filterBounds(roi) .filterDate("2023-10-01", "2024-06-30") .filter(ee.Filter.lte("CLOUDY_PIXEL_PERCENTAGE", 30)); // 第二步:像素级云掩膜函数 function maskClouds(img) { var scl = img.select("SCL"); var cloudMask = scl.gte(3).and(scl.lte(7)); // 3~7对应云影、低云、中云 return img.updateMask(cloudMask.not()); } // 第三步:逐景应用掩膜 var s2Masked = s2.map(maskClouds);这段代码里,filter lte 30 是影像级初筛,maskClouds 是像素级精细处理。SCL 波段中 3 到 7 分别覆盖云影、低云、中云这几类,gte(3) 和 lte(7) 把一次性圈出来,取反之后 updateMask,云像元就置为无效值。后续算 NDVI 时这些像元不参与统计,不会拖出假的低值或高值。
参数说明:CLOUDY_PIXEL_PERCENTAGE 的阈值我一般取 20 到 30。取太严(比如 5),正值返青期的影像会被大量过滤,因为春季多云;取太松(比如 40 以上),虽然影像多了,但云掩膜会把大片有效像元抠掉,数据利用率太低。SCL 掩码区间不是死的,有的平台 SCL 编码不同,你需要先用控制台打印一景影像的 SCL 分布直方图,确认 3 到 7 确实对应云和阴影之后再固化这个函数。
2.3 构建干净的多波段堆叠:从集合到影像
随机森林训练需要的是“波段 × 样本”的表格结构,所以影像集合要先合成一个多波段影像。常见做法是分物候段取代表影像,计算指数,再合成。我在项目里通常选 4 到 5 个关键日期窗口,每个窗口取云量最低的一景:
// 关键日期窗口:每个窗口取云量最低的一景 var keyDates = [ ["2023-10-01", "2023-11-10"], // 出苗 ["2023-12-01", "2024-01-15"], // 越冬 ["2024-03-01", "2024-04-05"], // 返青 ["2024-04-10", "2024-05-20"] // 抽穗灌浆 ]; var stackList = keyDates.map(function(window) { var img = s2Masked .filterDate(window[0], window[1]) .sort("CLOUDY_PIXEL_PERCENTAGE") .first(); var ndvi = img.normalizedDifference(["B8", "B4"]).rename("NDVI"); var ndwi = img.normalizedDifference(["B8", "B11"]).rename("NDWI"); return img.addBands([ndvi, ndwi]); }); var featureStack = ee.ImageCollection(stackList).toBands();这段代码的关键是 sort("CLOUDY_PIXEL_PERCENTAGE").first(),直接从窗口里挑云量最低的那一景,省去手动逐景检查。合成后的 featureStack 波段数大约等于 4 窗口 ×(10 个原始波段 + 2 个指数)≈ 48 个,每个样本点拿到的是 48 维特征向量。
注意:不同在线平台对 toBands 的波段数上限不一样,我遇到过窗口一多就报 too many bands 的情况。如果报错,把原始波段精简到 B4、B8、B11、B12 这 4 个核心波段,再加 2 个指数,单窗口 6 个波段,4 个窗口 24 维,特征量足够而且不触发限制。
这里有个实操细节:各波段空间分辨率不一致,B2-B4、B8 是 10 米,B11、B12 和 SCL 是 20 米。做堆叠前建议统一重采样。我在函数里显式做一次 nearestNeighbor 重采样,虽然不如 bilinear 平滑,但分类任务里保持原始像元值不容易引入边缘伪像,尤其田块边界很碎的时候更明显。
3. 训练样本与特征工程:让随机森林知道冬小麦长什么样
3.1 样本从哪里来:实采点、历史影像解译与样本量底线
随机森林是监督分类,样本质量直接决定上限。我见过最多的翻车原因是训练样本和验证样本同源,精度报到 95%,换一个县立刻崩。冬小麦提取场景里,样本来源无非两种:一是带手持 GPS 去田头打点,记录冬小麦和非冬小麦;二是基于历史高分影像做目视解译,在自己能确认的田块内取点。第二种时间紧时更常用,但要注意只选“确定是冬小麦”和“确定不是”的像元,拿不准的别勾,宁缺毋滥。
样本量方面,我的底线是每类不少于 200 个像元点,且空间上尽量分散。冬小麦样本要覆盖长势好的、一般的、差的田块,否则模型只认识长势好的,分类图上会出现大面积漏提。非冬小麦样本要覆盖其他作物、林地、水体、建筑、裸土这些常见地类。如果研究区里有油菜,那非冬小麦样本必须给足,它和冬小麦在冬春季光谱太像,给少了模型根本分不开。
样本准备阶段的另一件事是检查样本点的像元值合法性。用 featureCollection 提取特征时会返回每个点在所有波段上的值,有的点在云掩膜区域是 null,这些点要么剔除,要么用相邻时相替代。我习惯在训练前先做一次过滤,把任意特征为 null 的样本点全部删掉,这一步能避免后面训练时报错或者指标虚高。
3.2 特征集怎么定:光谱、指数、时序三件套
随机森林能自动挑特征,但喂什么它学什么。冬小麦提取的有效特征可以分三层理解。
第一层是单期光谱反射率。B4 红波段、B8 近红外、B11 短波红外这三组组合,能区分绿色植被和裸土、水体,但对冬小麦和油菜的区分度很低。第二层是植被指数。NDVI 最核心,NDWI 对水分敏感,能在冬春季把冬小麦和干枯杂草分开;EVI 在生物量高时不易饱和,更适合抽穗期用。我一般在单景影像上同时算 NDVI 和 NDWI,并保留原始光谱,让随机森林自己决定用谁。第三层是时序物候特征,这是多时相方法相对单时相最大的增量:NDVI 全年最大值、最大值出现日期、上升速率、生长季长度。冬小麦的时序特征是峰值高(0.7 以上)且峰值出现在 4 月下旬到 5 月上旬;油菜峰值出现偏早,露地蔬菜则常见多个峰。
实现上,我不用显式把“峰值日期”这类复杂指标算出来,而是把 4 个窗口的 NDVI 直接作为波段堆叠进特征影像,随机森林能靠窗口间的数值关系隐式学到这些物候指标,只要窗口选得对。对于 24 维的小特征集,再显式加物候指标收益有限;但对 48 维的完整特征集,随机森林自己就能完成筛选。
3.3 特征提取与数据集构建:从影像到表格
在 Gis online 上,特征提取最常用的是 sampleRegions,但细节不少。下面是我常用的完整流程,顺便把第 2 章里的指数计算封装成函数:
// 封装指数计算,供多时相堆叠和特征提取复用 function addIndices(img) { var ndvi = img.normalizedDifference(["B8", "B4"]).rename("NDVI"); var ndwi = img.normalizedDifference(["B8", "B11"]).rename("NDWI"); return img.addBands([ndvi, ndwi]); } // 在训练样本点上提取特征值 var samples = featureStack.sampleRegions({ collection: trainingPoints, // 训练样本点集合 scale: 10, // 与影像重采样后的分辨率一致 geometries: true // 保留坐标,便于检查 }); // 过滤掉含 null 特征的样本 var cleanSamples = samples.filter( ee.Filter.notNull(featureStack.bandNames()) ); // 随机打乱并按 7:2:1 划分训练、验证、测试集 var shuffled = cleanSamples.randomColumn("random"); var train = shuffled.filter(ee.Filter.lte("random", 0.7)); var valid = shuffled.filter(ee.Filter.gt("random", 0.7)) .filter(ee.Filter.lte("random", 0.9)); var test = shuffled.filter(ee.Filter.gt("random", 0.9));sampleRegions 的 scale 参数必须和影像实际分辨率一致,否则平台会自动重采样,特征值和影像显示对不上。这里有个很隐蔽的坑:featureStack 如果重采样到 10 米,scale 必须设 10;如果某平台默认 20 米,你就要在重采样后显式指定 scale,否则后续分类的面积统计自带系统性偏差。
数据集划分我用 7:2:1 而不是常见的 8:2,原因有二:一是冬小麦分类样本本身不多,测试集留 10% 也有三四十个点,够算混淆矩阵;二是验证集和测试集分开,能防止你在验证集上调参调到过拟合而不自知。
注意:如果样本点之间空间自相关严重,比如相邻像元都取做样本,7:2:1 随机划分的测试集精度仍然会虚高。严格做法是按田块整体划分,把一块田的全部像元划入训练或测试,而不是随机到像元级。
4. 随机森林建模:参数设定与训练流程
4.1 三个直接影响冬小麦分类效果的参数
随机森林在遥感分类里几乎不用调参,但有三个参数在冬小麦提取场景里值得专门设置。第一个是树的数量 numberOfTrees,默认 100 能跑,但特征维度 30 到 50 时我习惯用 300,树再多精度基本不涨,只增加训练时间。第二个是 variablesPerSplit,每个节点随机选取的特征数,平台默认是特征数的平方根。如果训练集精度极高、验证集明显偏低,说明过拟合,可以把它从默认的 6 调大到 8,让每棵树有更多机会选到关键特征;反过来 OOB 误差一直下不来就适当调小。第三个是 minLeafPopulation,最小叶节点样本数,默认 1 太容易过拟合,对孤立噪声点敏感,我设 5 到 10,分类图的碎斑会明显减少,比事后滤波更治本。
| 参数 | 推荐范围 | 设置逻辑与坑 |
|---|---|---|
| numberOfTrees | 200~500 | 精度 200 棵后饱和,树多了只费时间 |
| variablesPerSplit | 6~8 | 默认取特征数平方根,过拟合时往大调 |
| minLeafPopulation | 5~10 | 抗噪声,减少碎斑块 |
| bagFraction | 0.6~0.7 | 剩余样本用于 OOB 误差估计 |
4.2 训练与分类:一套可以直接复制的流程
第一次跑通,我用下面的完整代码。训练前先取出特征波段名,确保分类器和影像波段严格对应,这个对不上是我见过最频繁的报错来源。
// 取出特征波段名(不包含几何等非波段字段) var featureBands = featureStack.bandNames(); // 创建随机森林分类器并训练 var classifier = ee.Classifier.smileRandomForest({ numberOfTrees: 300, variablesPerSplit: 7, minLeafPopulation: 8, bagFraction: 0.6, seed: 42 }).train({ features: train, // 训练样本集 classProperty: "class", // 样本属性中的类别字段 inputProperties: featureBands // 参与训练的特征波段 }); // 对全研究区进行预测 var wheatMap = featureStack.classify(classifier, "predicted"); // 用验证集算精度 var validationError = valid .classify(classifier, "predicted") .errorMatrix("class", "predicted"); print("混淆矩阵(验证集):", validationError); print("总体精度:", validationError.accuracy()); print("Kappa:", validationError.kappa());这里的 bagFraction 是每次抽样训练一棵树时使用的样本比例,0.6 到 0.7 常用,剩下的样本留作袋外数据估计 OOB 误差。seed 固定为 42,保证每次训练结果可复现,写报告和做对比实验时特别重要,不设 seed 的话每次跑出来面积可能差出几十亩。classProperty 要和你给样本点设置的属性字段名完全一致,常见错误是样本点里叫 typename 而这里写的 class,一训练就报属性不存在。
分类完成后先别急着看面积,第一步看验证集混淆矩阵。冬小麦提取我重点盯三个数字:生产者精度,即真实冬小麦有多少被分对;用户精度,即分出来的冬小麦有多少真正确实是冬小麦;以及错分去向。如果冬小麦大部分错分到油菜,说明返青期特征区分度不够,回去查那个窗口的 NDVI 时序;如果错分到裸土,多半是某景影像有云残留,在 NDVI 上制造了假低值。
4.3 特征重要性:回头看,比盲目调参更有效
训练完强制性做一件事:打印特征重要性。随机森林内部会统计每个特征在分裂中带来的信息增益,Gis online 上通过 explain 输出,注意不同平台返回的键名可能不同,以实际输出为准。
var importance = classifier.explain().get("importance"); print("特征重要性:", importance);拿到排序后,我一般做一次“裁剪重训”:把重要性排后 15% 的特征从 featureBands 里去掉,重新训练,比较验证精度。冬小麦提取里被裁掉的往往是某个窗口用处不大的原始波段,比如 B1 海岸波段,或者越冬期的近红外波段,这类波段对分类是纯噪声,删掉后精度不降,训练速度反而快。如果裁剪后精度下降,说明特征维度已经比较紧,把特征加回去,不折腾。
注意:特征重要性高不代表物理意义就明确。返青期 NDVI 重要性高符合物候预期,但某个短波红外波段重要性高的时候,我会去查这个波段是不是有传感器噪声,别因为模型用了它就笃定写进报告里。
5. 精度验证与避坑:冬小麦提取最容易翻车的位置
5.1 验证数据要独立、指标要交叉,三个数字一起看
精度验证的第一原则是验证样本必须独立于训练样本,而且最好来自不同田块。如果验证样本和训练样本出自同一批田块,空间自相关会让精度虚高。某一次我用同一批田块划分样本,总体精度报 89%,换成独立田块验证直接掉到 72%,这个差距足够说明问题有多严重。
精度汇报不要只报一个总体精度。冬小麦提取必须同时给出生产者精度和用户精度,两者差异本身就是问题信号。生产者精度高、用户精度低,说明模型把大量非冬小麦划进来了,面积偏大;反过来生产者精度低、用户精度高,说明漏提严重,面积偏小。农业统计业务更看重生产者精度,因为漏提会导致补贴漏发;保险定损更看重用户精度,因为错分会把不该赔的地块赔进去。你要按业务目标选主导指标,把另一个当作约束条件。
还有一个常被忽略的是面积偏差校正。随机森林分类是逐像元硬分类,如果混淆矩阵显示错分率超过 5%,面积偏差可能超过实际面积的 10%。这时候只汇报“提取面积 3.2 万亩”是不负责任的,至少要加一句“按验证集误差校正后处于 3.08 到 3.35 万亩之间”,业务方听了反而更有底。
5.2 三个必踩的坑:云残留、物候偏移、样本不平衡
坑 1:云残留制造假“低 NDVI”碎斑。
现象:分类图上冬小麦田块内部出现随机分布的圆形碎斑,被分成油菜或裸土,面积损失 8% 到 12%。
原因:某个窗口的影像 SCL 波段没有掩掉薄云或云影,这些像元在 NDVI 上表现为低值,物候曲线看起来像“冬天就死了”,随机森林把它判给非冬小麦。
解决:第一,打印该窗口实际用的影像 ID,人工看一眼真彩色合成,确认云掩膜是否真的生效;第二,把 SCL 掩膜区间从 3 到 7 扩成 2 到 7,把暗区也算进去,但要让步:这会多掏掉一些有效像元,如果该窗口有效像元占比低于 70%,直接换时间窗重新挑景。
坑 2:不同年份物候偏移,去年的模型今年翻车。
现象:用去年训练的模型分类今年影像,返青期大片漏提,分类面积比官方统计少 20%。
原因:暖冬导致返青提前,或者倒春寒导致滞后,相同日期在今年和去年对应的物候期完全错位。随机森林学的是“返青期 NDVI 要涨到 0.5 以上”这种绝对值,物候偏移几天,特征分布整体移动,泛化就失效。
解决:每年开春先跑一版只用当年 1 到 4 月影像的快速版,和往年结果对比;物候显著偏离的年份,把当年的样本量增加 50% 再重新训练。冬小麦提取模型的保质期通常是一到两个物候季,不要指望一个模型用五年。
坑 3:样本不平衡导致模型“偏心”。
现象:冬小麦样本 500 个,非冬小麦样本 2000 个,模型偏向多数类,冬小麦被大量漏提,生产者精度只有六成。
原因:随机森林的默认目标函数是整体准确率最大化,少数类错误对整体代价影响小,所以模型天然偏向多数类。
解决:一是样本采集时就控制冬小麦和非冬小麦的样本量比值在 1:1 到 1:2 之间,并按田块分层抽样,保证每类代表性地块数相当;二是实在补不了样本,可以在训练时给冬小麦样本加权,把权重提高到非冬小麦的 1.5 到 2 倍。只改混淆矩阵不重训,治标不治本,不建议。
5.3 面积统计与结果导出:让成果能直接交差
分类完成后别直接拿原始分类结果统计面积,先做后处理。我用的顺序是:先做 3×3 焦点众数滤波,把零星碎斑合并进周围多数类;再做最小斑块过滤,把小于等于 3 个像元(10 米分辨率下约 300 平方米)的孤立斑块删除,冬小麦田块不可能碎到单像元级别,这种斑块基本是噪声。
面积统计用平台自带工具,先统计各类别像元数,再乘以单像元面积。注意 10 米像元的实际地面尺寸在不同纬度有差异,全程用平台自带的面积计算函数,别手写固定系数。
// 多数滤波后计算各类面积 var filtered = wheatMap.focalMode(3, "square"); var areaDict = filtered.multiply(ee.Image.pixelArea()) .reduceRegion({ reducer: ee.Reducer.sum().group({ groupField: 0, groupName: "class" }), geometry: roi, scale: 10, maxPixels: 1e10 }); print("各类面积:", areaDict);这段代码 multiply(ee.Image.pixelArea()) 计算每个像元的真实地面面积,按类别分组求和。maxPixels 设到 1e10 是避免大研究区报内存错,如果还报错,改成按区块导出后本地汇总。导出时选 GeoTIFF,别选 PNG,PNG 会丢失坐标系,后续在本地 GIS 里做矢量化和面积复核都麻烦。
6. 把模型用到下一季:迁移与复算的实操习惯
一个容易让人产生错觉的地方是随机森林模型训练一次就能年年用。我前面说过保质期通常只有一到两个物候季,但实际业务里不可能每年从头做一轮完整标注。我的做法是把模型拆成“冻结的部分”和“每年更新的部分”。
冻结的是特征集、随机森林参数、样本选取规则;每年更新的是当年返青期和抽穗期的少量新样本,一般补 60 到 80 个点,约占老样本量的 15% 到 20%。把这些新样本和老样本合起来重新训练,老样本不删除,但按时间衰减给它们打更低的权重。这里有个 trade-off:老样本太多会让模型对往年物候“太恋旧”,老样本太少又丢失了稳定特征,我一般保持老样本和新样本比例在 4:1 到 5:1 之间。
复算的一个快速验证手段是“关键期面积环比”。每年返青期过后,只用当年 1 到 4 月影像先跑一次临时分类,统计冬小麦面积,和去年同期的临时结果对比。偏差在 ±5% 以内就继续用;偏差超过 10% 就说明物候偏移或极端天气影响,启动当年样本补充。这个习惯帮我提前发现问题,而不是等到 6 月收割季才面对一张精度崩掉的图。
我这几年做冬小麦提取,最值钱的经验就是别把时间耗在调随机森林参数上,把精力放在数据选择、物候窗口和样本质量三件事上。随机森林是个黑匣子,但对遥感分类来说它足够稳定、足够快,你喂给它干净的多时相时序特征,它就能给你一张能交差的冬小麦分布图。在 Gis online 上跑通一遍这个流程之后,换数据、换年份、换区域,核心逻辑都不变,变的只有物候窗口和样本。希望帮到你。
本文还有配套的精品资源,点击获取