去年帮一个做农业遥感的朋友处理一个地级市的农作物种植结构提取他本地那台工作站上挂着好几块硬盘Sentinel-2 L2A 从 2019 年攒到 2023 年光解压和整理元数据就花了两天。等他跑完云掩膜、重采样、叠波段、逐月合成、算 NDVI、再做时序平滑这一整套Python 进程已经挂了三整天中间因为内存爆掉重启了两次最后出来的分类图还被审稿人指出边缘田块碎斑太多。后来我把整套流程搬到 Google Earth EngineGEE上重做了一遍同样的区域、同样的年份、同样一套特征代码从两千多行缩到四百来行跑完一遍完整分类加精度评价不到四十分钟。这篇文章就把这套 GEE 农作物种植结构提取的路子掰开揉碎讲清楚影像怎么选、云怎么掩、时序特征怎么构造、样本怎么来、分类器参数怎么定、结果怎么导出以及我在这几个环节里踩过的、常规教程里基本不会提的那些坑。适合已经有遥感基础、想把种植结构提取从本地搬到云端的人也适合刚接触 GEE、想找一个完整落地案例练手的人。1. 把种植结构提取放到 GEE 上先算清本地处理这笔账1.1 种植结构提取的产出到底是什么东西很多人一上来就把种植结构提取等同于作物分类图这是第一个认知偏差。实际业务里一张合格的种植结构图至少要回答四个层面的问题第一层是作物类型也就是这块地种的是水稻、玉米、小麦还是大豆、棉花、果树第二层是种植制度即复种指数是一年一熟、一年两熟还是两年三熟华北平原的冬小麦—夏玉米轮作和东北的春玉米一熟制在同一张图上的表达方式完全不同第三层是空间格局包括田块尺度的连片程度、破碎化指数、作物种植重心迁移第四层是时间稳定性某些地块今年种玉米明年改种花生这种年际波动在单年分类图里根本看不出来。这四个层面的信息前三个靠单年时序特征就能拿到第四个必须做多年度分类再叠加分析。我在实际项目里见过太多只做一年分类然后强行推演复种指数的做法结果就是把今年休耕误判成一年一熟把临时改种误判成种植制度变化。所以从设计阶段就要想清楚你的分类结果最终要支撑哪个层面的分析这直接决定了特征集和样本方案怎么搭。1.2 本地处理一遍的隐性成本账我拿一个一万平方公里左右的地级市做个粗略估算这样你对自己要面对的数据量有个概念。Sentinel-2 L2A 单景覆盖大约 110 公里见方也就是一万两千多平方公里一景压缩后大约一个 G 上下。这个地级市一年的完整覆盖把双星重访和云覆盖冗余算进去大概需要 60 到 80 景做五年时序就是三四百景压缩包合计三四百 G解压之后算上中间产物轻松突破一个 T。真正折磨人的不是存储是计算。本地做时序合成通常要经历解压、重投影到统一投影、按云掩膜逐像元筛选、重采样到统一网格、逐月做中值合成、算植被指数、时间序列平滑、特征堆叠、模型训练、全区推测、众数滤波。这一条链里每一步都会产生一份中间文件每一步都可能因为一个坐标系不一致或者 NoData 处理不对而重来。而且这套流程天生没有可复现性——半年后你想改一个参数重跑硬盘上那些中间文件早就被删了。更隐蔽的成本是并行度。本地做全区推测只能一块一块地切瓦片、多进程跑切瓦片本身要写不少调度代码进程数受限于物理核数和内存。我这边的经验是一台 32 核 128G 的机器做一万平方公里、五年时序、十个类别全流程跑通加上调试没有三四天出不来。1.3 GEE 在这件事上的能力边界搬到 GEE 之后上面那条链里最耗时的部分——数据存储、重投影、逐像元筛选、并行计算——全部由平台承担。你只需要写清楚要什么不用管怎么调度。这就是它最大的价值把工程复杂度换成了一次性的逻辑设计成本。但 GEE 不是万能的有几条边界必须提前知道。第一交互式计算有时间限制Code Editor 里点一次 Run 超过几分钟基本会超时凡是全区级别的计算都必须走Export用批处理任务在后台跑。第二不适合做需要反复随机读写的深度学习训练虽然现在有 TensorFlow 接口但数据搬运和样本预处理的代价很大涉及深度学习的场景我一般还是把样本导出到本地做。第三没有真正的状态每个脚本每次运行都是从原始影像重新算一遍中间结果要显式地存成 Asset 才能复用这一点和本地工作流的心智模型完全不一样。提示如果你是第一次做种植结构项目建议先用一个县的范围把完整流程走通确认特征和样本方案没问题再放大到全市或者全省。直接上大区域的结果通常是跑了两小时报一个内存溢出然后你连错在哪都找不到。2. 影像源与云掩膜Sentinel-2、Landsat 8/9 到底怎么配2.1 三种数据源在种植结构里的分工GEE 上做作物分类绕不开三类数据源。我把它们的核心参数列成表方便你对着自己的项目周期和精度要求做取舍。数据源空间分辨率重访周期时间覆盖在种植结构里的定位Sentinel-2 SR10m可见光近红外/20m约 5 天双星2017-03 至今主力数据田块尺度分类的唯一现实选择Landsat 8/9 C2 L230m约 8 天双星2013 至今补 Sentinel-2 之前的年限或填补特定时段的空缺MODISMOD13Q1/MCD12Q2250m/500m16 天/年2000 至今长时序物候先验、大区域趋势分析不参与田块尺度分类我的一般组合是Sentinel-2 作为主特征源负责 2017 年之后所有年份如果项目要求 2013 到 2017 年的结果就用 Landsat 8 单独建一套特征两段结果在重叠期做一致性检验MODIS 的物候产品则作为辅助特征塞进特征集尤其是 MCD12Q2 提供的 Greenup、Senescence 这类物候期参数在大区域上能明显提升早稻、冬小麦的区分度。这里有个容易被忽略的点Landsat 和 Sentinel-2 的空间分辨率差三倍混用会引入尺度不一致误差。如果一块田只有十几米宽Landsat 像元里会混进大量非作物地物NDVI 曲线直接被拉平。所以在一年两熟、田块细碎的地区我建议只用 Sentinel-2Landsat 只用在田块面积大、地类单一的区域。2.2 云掩膜的三种做法和各自的残留问题云掩膜是种植结构提取里最容易被低估的环节。掩膜做得不干净时序曲线上就会出现莫名其妙的高值或低值毛刺平滑算法再好也救不回来。GEE 上 Sentinel-2 云掩膜有三条常见路径各自的毛病我列一下。第一条路是 SCL 分类层。Sentinel-2 L2A 自带SCL波段用整数编码区分水体、云、云影、植被、裸土等类别直接筛掉 3、8、9、10、11 这几个值就行。代码写起来最省事但它在薄云和云边缘上处理得不干净经常漏掉半透明的卷云也会把部分云影误判成植被。第二条路是 QA60。这个波段用位运算标记云和卷云是老教程里最常见的写法。但要注意处理基线 04.00 之后发布的数据里 QA60 已经不再更新用它对旧数据没问题对新数据会失效。第三条路是 Cloud Score。这是目前我认为效果最好的一条GEE 里可以直接链接到 Sentinel-2 和 Landsat 的集合上用cs_cdf波段做阈值筛选。它的优点是能识别薄云、云影和霾保留更多有效观测代价是数据要额外做一次链接。// Sentinel-2 云掩膜SCL Cloud Score 双重筛选 var s2 ee.ImageCollection(COPERNICUS/S2_SR_HARMONIZED); var csPlus ee.ImageCollection(GOOGLE/CLOUD_SCORE_PLUS/V1/S2_HARMONIZED); var csBand cs_cdf; var s2Masked s2 .linkCollection(csPlus, [csBand]) // 把 cs_cdf 链接到每景影像上 .filterBounds(roi) .filterDate(2023-01-01, 2023-12-31) .map(function(img) { var scl img.select(SCL); // SCL: 3云影, 8中概率云, 9高概率云, 10卷云, 11雪 var sclMask scl.neq(3).and(scl.neq(8)).and(scl.neq(9)) .and(scl.neq(10)).and(scl.neq(11)); // Cloud Score 概率阈值0.6 是比较稳的经验值 var csMask img.select(csBand).gt(0.6); var mask sclMask.and(csMask); return img.updateMask(mask) .divide(10000) // 反射率缩放务必做 .copyProperties(img, img.propertyNames()); });Landsat 的云掩膜走的是另一套逻辑C2 L2 产品用QA_PIXEL波段的位标记需要按位读取。位 1 是膨胀云、位 3 是云、位 4 是云影、位 5 是雪把这几位判断为 0 才是干净像元。// Landsat 8/9 C2 L2 云掩膜 function maskLandsat(img) { var qa img.select(QA_PIXEL); var mask qa.bitwiseAnd(1 1).eq(0) .and(qa.bitwiseAnd(1 3).eq(0)) .and(qa.bitwiseAnd(1 4).eq(0)) .and(qa.bitwiseAnd(1 5).eq(0)); // C2 L2 的 SR 波段需要乘缩放系数再减偏移量 return img.updateMask(mask) .select([SR_B.*], [B2,B3,B4,B5,B6,B7]) .multiply(0.0000275).add(-0.2); }2.3 跨传感器拼接时的波段对照与缩放系数如果你确实需要把 Landsat 和 Sentinel-2 拼到一起用波段名称和物理量必须对齐否则指数算出来全是错的。物理意义Sentinel-2Landsat 8/9 C2 L2中心波长对应关系蓝B2B2约 490nm vs 480nm绿B3B3约 560nm vs 560nm红B4B4约 665nm vs 655nm近红外B8B5约 842nm vs 865nm短波红外1B11B6约 1610nm vs 1610nm短波红外2B12B7约 2190nm vs 2200nm红波段和近红外的中心波长相隔十几到二十纳米NDVI 会有系统偏差一般做法是对其中一个做线性经验校正或者干脆在特征层面各算各的、让分类器自己去学。缩放系数这块Sentinel-2 SR 是除以 10000Landsat C2 L2 是乘 0.0000275 再减 0.2两个都做完之后才是无量纲反射率才能代入 EVI 这种带常数项和蓝色波段的公式。我见过最典型的错误就是 NDVI 算了缩放、EVI 忘了算结果 EVI 全是大负数还不报错只是分类精度悄悄掉下去。3. 把时间维度的信息榨干从 NDVI 曲线到物候特征3.1 单时相分类在一年两熟区为什么必败先讲一个真实的翻车案例。有位同学做华北平原某县的种植结构选了一景 6 月中旬的 Sentinel-2 影像做监督分类用地物光谱差异硬分结果玉米和小麦混得一塌糊涂。原因很简单6 月中旬正是冬小麦收割、夏玉米刚播种出苗的窗口期小麦地块这时候要么是收割后的裸土加秸秆要么已经播下玉米但还没出苗而常年一熟的春玉米地块这时候刚刚进入快速生长期绿度还不高。两个类别在这个时间点上的光谱几乎重合。这就是种植结构和地表覆盖的根本区别种植结构必须靠时间维度的差异来定义而不是靠某一时刻的光谱。冬小麦的特点是秋季播种、越冬、次年 3 到 5 月快速返青并达到 NDVI 峰值、6 月收割夏玉米是 6 月播种、7 到 9 月达到峰值、10 月收获。把这一年的 NDVI 曲线画出来就是典型的双峰形态而且两个峰值的位置和幅度直接编码了种植制度信息。一熟制春玉米则是单峰峰值在 7 到 8 月。这两种曲线放在特征空间里分得干干净净但用单时相影像去看它们可能完全一样。3.2 平滑重建谐波拟合与 SG 滤波在 GEE 里的实现原始 NDVI 时序一定是有噪声的云残留、观测几何变化、大气校正误差都会造成毛刺。所以第一步是平滑重建目标是求出一条干净、连续、能反映真实物候过程的曲线。GEE 上主流有两种做法。谐波拟合Harmonic Regression是最稳也最好实现的一种。它的思路是把一年的 NDVI 曲线用几个正弦和余弦项叠加去逼近拟合出来的系数本身就是压缩后的物候特征。一阶谐波就能抓住一年一个峰二阶能抓住一年两个峰这对一年两熟区特别关键。// 谐波拟合一阶 二阶numX5 表示 常数项 sin/cos(1阶) sin/cos(2阶) var dependent ndviCollection; // 已经做过云掩膜和月合成的 NDVI 集合 var independent ee.ImageCollection( ee.List.sequence(0, dependent.size().subtract(1)).map(function(i) { var t ee.Number(i).multiply(Math.PI / 6); // 假设是逐月角度步长 π/6 return ee.Image.constant([ 1, // 常数项 t.sin(), // 一阶 sin t.cos(), // 一阶 cos t.multiply(2).sin(), // 二阶 sin t.multiply(2).cos() // 二阶 cos ]).rename([const,sin1,cos1,sin2,cos2]) .set(system:time_start, dependent.toList(1).get(0).get(system:time_start)); }) ); // 把自变量和因变量合成同一个集合前 numX 个波段当 X后面的当 Y var combined dependent.select([NDVI]).combine(independent); // 注意波段顺序 var regression combined.reduce(ee.Reducer.linearRegression(5, 1));拟合完之后从coefficients波段里能直接读出振幅和相位振幅是sqrt(sin^2 cos^2)相位是atan2(cos, sin)一阶振幅代表年内绿度变化强度相位代表峰值出现的时点二阶振幅则直接指示双峰的强度。这三个量放到分类器里比原始 NDVI 的判别能力强得多。Savitzky-Golay 滤波效果也好但 GEE 里没有现成的实现。它的原理是在时间窗口内做局部多项式最小二乘拟合本质是一个固定系数的卷积核。移植到 GEE 上的思路是把时序影像转成数组、用arrayReduce沿时间轴做加权求和但实现起来比较绕而且数组维度和空值处理很容易踩坑。如果只是做一年内的平滑我建议优先用谐波拟合如果一定要 SG可以先用谐波拟合做一次去噪再叠加一个短窗口的移动平均实际效果差距不大。3.3 特征集到底放哪些一份可以直接抄的清单特征不是越多越好冗余特征会稀释随机森林的判别力还会拖慢计算。我把这几年用下来比较稳的一套清单整理出来你可以根据自己的作物类型做增删。特征类别具体内容物理意义逐月中值指数12 个月的 NDVI、EVI、NDWI 中值保留完整物候过程是分类的主体特征分位数全年 NDVI 的 p10、p25、p50、p75、p90描述绿度分布形态对休耕、裸土敏感峰值特征NDVI 最大值、最大值出现日期、生长季积分区分单峰/双峰量化长势差异谐波系数一阶、二阶的振幅和相位压缩的物候形态计算量小判别力强红边指数NDRE (B8-B5)/(B8B5)对叶绿素含量和生育期敏感水稻尤其有用水分指数NDWI、LSWI含短波红外水稻淹水期识别是水稻区别于旱作的关键物候产品MCD12Q2 的 Greenup、Senescence、Dormancy提供先验物候期适合大区域约束地形SRTM 高程、坡度山区区分梯田作物与林地的辅助变量先验地表覆盖Dynamic World 各类别概率波段提供类别先验能显著压低离群误分这套特征加起来大约六十到八十个波段对 Sentinel-2 十五个景左右的月合成集合来说计算量在可接受范围内。有一个坑必须提醒逐月中值合成遇到整月全是云怎么办。如果某个月所有观测都被掩掉了那个月的波段就是全 NoData后续分类时这一景会被整片跳过或者报错。我的处理办法是加一个后备对每个月中值合成后再做一次线性插值用前后月份的值填补空缺代码里用一个简单的focal_mean沿时间维度做确实不现实更靠谱的做法是先判断某个像元的有效观测数少于阈值就用谐波拟合的结果来补。这个细节在大多数教程里都不会写但它在多雨地区是决定成败的。4. 样本是分类精度的天花板采集、迁移与分层4.1 野外调查点、高分目视判读、已有产品迁移这三条路分类器的上限由样本决定这句话我做项目越久越认同。样本来源无非三条路成本和适用场景完全不同。野外实地调查是最可靠的GPS 打点记录作物类型优点是不用怀疑标签真值。但它的限制也很明显一趟野外只能覆盖有限面积成本高而且一个点代表一个田块还是代表一个区域本身就是个需要定义的问题。我的经验是野外点主要用来做锚点用来验证其他方式生成的样本而不是直接当训练集主力。高分辨率影像目视判读是实际项目里最常用的。用天地图、Google 影像或者项目自有的高分影像在 GEE 里逐块画多边形判断作物类型。这里有个技术细节目视判读得到的是多边形不是像元直接拿去训练会引入大量边界混合像元。我的做法是把多边形向内做一次腐蚀在 GEE 里可以用reduceNeighborhood或者直接人工内缩只保留田块内部的核心像元这样训练样本纯度能提升一大截。已有产品迁移是最高效的一条路前提是你有一份可信度较高的历史分类结果或者公开产品。国外有 CDL 这种年度作物分类产品国内有部分公开的耕地和作物制图成果加上 Dynamic World、ESA WorldCover 这类地表覆盖产品。它们的类别体系不一定和你需要的完全一致但拿来生成初始样本、再做人工抽检修正效率比从零画多边形高很多。4.2 样本迁移的具体做法和一致性检验样本迁移是这几年的热门做法核心思想是用上一年的分类结果去生成下一年的训练样本避免每年重复做野外调查。但它有个硬前提地块的种植结构在两年之间没有发生大规模变化。如果第二年遇到大范围改种或者撂荒迁移样本会把噪声直接喂给分类器。我的做法是加一层一致性检验具体分四步。第一步拿上一年分类图对每一类做边界腐蚀只保留田块内部的像元这一步能把混合像元基本清干净。第二步用当年的时序特征在这批像元上计算类内一致性比如看这批像元当年 NDVI 曲线的方差方差过大的整块剔除。第三步用当年已有的少量实地点或者高分目视点对迁移样本做交叉验证如果某个类的用户精度低于一个阈值我一般设 80%就对这个类做降采样或者直接弃用。第四步把最终样本与上一年的样本做空间叠加落在变化热点区的样本单独标记训练时降权或者剔除。// 样本迁移中的边界腐蚀用众数滤波去掉田块边缘像元 var classMap ee.Image(projects/xxx/assets/crop_map_2022); // 先做一次 3x3 众数滤波边缘像元会被邻域主导类替换 var smoothed classMap.focalMode({radius: 1, units: pixels}); // 只保留滤波前后类别一致的像元这些就是田块内部的纯像元 var coreMask classMap.eq(smoothed); var samples smoothed.updateMask(coreMask) .stratifiedSample({ numPoints: 300, classBand: classification, region: roi, scale: 10, geometries: true });4.3 类间不平衡与训练验证集的切分种植结构分类的类别分布天然极不平衡水稻可能占 40%花生可能只占 2%。如果直接用随机采样花生这个类可能在整个训练集里只有几十个点分类器会倾向于把它当成噪声忽略掉用户精度极低。解决办法是按类分层采样并且给少数类设最低样本数下限。我在 GEE 里用stratifiedSample按类别采样每个类至少保证 300 到 500 个样本。样本量不是越多越好超过一定数量后精度提升会趋于平缓反而增加训练时间。规模上一个县域做八到十个类别总样本量在 3000 到 5000 之间基本够用。切分训练集和验证集的时候要注意不能按像元随机切要按空间块切。同一块田里的相邻像元高度自相关如果随机切分验证集里会混进和训练集空间上紧邻的像元验证精度会被系统性高估有时候能虚高十几个百分点。我一般用randomColumn生成随机数然后按 0.7 的比例切但在此之前会先把样本按空间聚类分组保证同一田块的像元落在同一组里。这一步比较麻烦如果不想做退而求其次的做法是训练集和验证集用不同年份的样本天然规避了空间自相关。5. 分类器与后处理随机森林不是按个按钮就完事5.1 随机森林参数的定法GEE 里的smileRandomForest是种植结构分类的默认选择它对特征量纲不敏感、能处理非线性、能输出变量重要性几乎没什么明显短板。但有三个参数需要认真对待。var classifier ee.Classifier.smileRandomForest({ numberOfTrees: 200, // 树的数量 variablesPerSplit: null, // 默认取特征数的平方根 minLeafPopulation: 3, // 叶节点最小样本数 bagFraction: 0.5, // 每棵树的样本抽样比例 seed: 42 // 固定随机种子保证可复现 }); var trained classifier.train({ features: trainingSamples, classProperty: class, inputProperties: featureBands });numberOfTrees设多少我的经验是 100 到 300 之间就够再往上精度提升很小但计算时间线性增长。minLeafPopulation建议设成 3 到 5默认的 1 会让每棵树都长得特别深容易过拟合噪声样本。seed一定要固定否则每次运行结果都不一样你根本没法判断是参数改动起了作用还是随机性在捣乱。bagFraction保持默认 0.5 就行样本量特别小的时候可以调到 0.8 让每棵树看到更多样本。5.2 用变量重要性回看特征集砍掉冗余训练完之后explain()会输出每个特征的贡献度这是被严重低估的一个环节。我一般会把重要性排序打出来看看排在前十的是哪些排在最后二十位的又是什么。var importance trained.explain(); print(importance); // 输出每个特征的 importance 值实际跑下来重要性最高的通常是一阶谐波相位、NDVI 峰值日期、7 到 9 月的 NDVI 中值、LSWI 这几个。有意思的是很多新手喜欢堆一大堆指数什么 SAVI、MSAVI、GNDVI、NDMI 全放进去结果重要性榜单上那些高度相关的指数扎堆排在末尾互相稀释。我的做法是如果两个特征的相关系数超过 0.9只保留重要性更高的那个。这一步通常能把特征数砍掉三成模型训练更快精度几乎不降有时候还会略微上升。5.3 众数滤波、连通域去碎斑和对象化分类分类器输出的原始结果一定是椒盐噪声式的尤其在田块边界和混合像元区域。后处理做得好不好直接决定成果能不能用。最基础的是众数滤波用 3x3 或 5x5 的窗口取邻域众数把孤立的错误像元抹掉。缺点是它也会把真实的细小田块一并抹掉所以窗口不能开太大我做 10 米分辨率一般只用 3x3。进阶一点的做法是对象化分类。先用 SNIC 算法对影像做分割得到一个个空间上连续、光谱上均质的分割对象然后把每个对象内部所有像元的特征做统计均值、中值、标准差用对象的统计特征训练分类器最后把分类结果赋回对象内的所有像元。这样出来的结果天然是田块尺度的不会有碎斑。// SNIC 分割 对象内部众数赋值 var seeds ee.Algorithms.Image.Segmentation.seedGrid(30); // 分割尺度30 像元约 300 米 var snic ee.Algorithms.Image.Segmentation.SNIC({ image: featureStack, size: 30, compactness: 5, connectivity: 8, neighborhoodSize: 128, seeds: seeds }); var segments snic.select(segments); // 把像元级分类结果按分割对象取众数 var objectClass classified.reduceConnectedComponents({ reducer: ee.Reducer.mode(), labelBand: segments });需要注意一点reduceConnectedComponents的效果取决于分割质量分割尺度设得太大会把相邻的不同作物田块合并成一个对象众数一取整片都变了设得太小又和像元级分类没区别。我的经验是分割尺度按当地平均田块大小的 0.5 到 1 倍来设平原地区 10 米分辨率一般取 20 到 40 像元。最后连通域计数可以用来清理小到不合理的图斑。比如一个 3 个像元的水稻图斑在 10 米分辨率下就是 300 平方米比最小田块还小基本可以肯定是误分。// 统计每个连通图斑的像元数小于阈值的并入背景 var count classified.connectedPixelCount({maxSize: 100, eightConnected: false}); classified classified.updateMask(count.gte(20)); // 小于 20 像元的图斑剔除6. 精度评价、导出参数与工程化踩坑清单6.1 混淆矩阵里真正要看的指标很多人做精度评价只看一个总体精度OA然后写一句分类精度达到 92%就交差了。这个数字在类别极不平衡的场景里基本没有参考价值——如果水稻占 80%我把所有像元都判成水稻OA 也能有 80%。真正要看的是三类指标。生产者精度PA反映的是这一类有多少被找出来了PA 低说明漏分严重对应的是地面真实面积被低估。用户精度UA反映的是被判为这一类的有多少是真的UA 低说明错分严重。Kappa 系数在类别严重不平衡时同样会失真我一般把它作为参考不作为主要依据。var validation validationSamples.classify({ classifier: trained, outputMode: CLASSIFICATION }); var cm validation.errorMatrix(class, classification); print(总体精度, cm.accuracy()); print(Kappa, cm.kappa()); print(用户精度, cm.consumersAccuracy()); print(生产者精度, cm.producersAccuracy());如果你的结果要用来估算面积那就必须用分层随机抽样的思路来设计验证样本并且按类别面积比例做加权。裸的像元计数去推面积在类别不平衡时误差能到百分之几十。这套方法在国际上已经有比较成熟的操作规范核心就是用混淆矩阵反推每个类的真实面积区间。6.2 导出参数和投影的选择导出环节的坑基本都出在参数上。我的导出模板大致是这样Export.image.toDrive({ image: classified.toByte(), description: crop_structure_2023, folder: GEE_exports, fileNamePrefix: crop_structure_2023, region: roi, // 建议用矩形范围不要用复杂多边形 scale: 10, // 必须显式指定默认会按影像原始分辨率 crs: EPSG:32650, // 用当地 UTM 分带不要用 EPSG:4326 maxPixels: 1e13, // 大区域必须抬高否则会报错 fileFormat: GeoTIFF, formatOptions: {cloudOptimized: true} // 输出 COG方便后续 QGIS 直接读 });几个要点解释一下。scale必须写清楚不写的话 GEE 会按输入影像的原始分辨率算混合数据源的情况下你不知道会输出什么。crs强烈建议用当地 UTM 分带而不是地理坐标系因为 EPSG:4326 的像元在纬度方向上不是正方形面积统计会偏。region用矩形包络框比用精确边界效率高很多精确边界会导致大量的空瓦片计算。maxPixels在超过一万平方公里的时候必须抬高到 1e13 量级否则会直接报错退出。另外如果导出到 Drive 的影像超过大约 30GGEE 会自动按瓦片拆成多个文件。这时候文件名会带后缀合并的时候要按文件名排序拼接别手工拼错行号。如果区域非常大我更推荐先Export.image.toAsset()存成 Asset再从 Asset 分批导出这样中间结果可复用不用每次重跑计算。6.3 内存、超时和 tileScale我踩过的几个坑这一节是我认为最有价值的部分因为这些都是报错之后才学会的东西。报错信息真实原因解决办法User memory limit exceeded单次聚合计算的数据量超过单用户内存上限加tileScale: 4或更高把计算切成更多小块Computation timed out交互式计算超过时间限制改用Export批处理不要点 RunToo many concurrent aggregations同一次reduceRegion并发聚合过多减小区域、提高 tileScale或拆成多次导出Image.reduceRegion: Parameter reducer is requiredreducer 没传或者传了 undefined检查变量名尤其在 map 函数里导出结果是全 0 或全 NoData类别值超了 byte 范围或者掩膜没清干净用toByte()前确认类别值不超过 255检查掩膜tileScale这个参数值得单独说。它的作用是控制每个计算瓦片的大小值越大瓦片越小单次计算占用的内存越少但总计算时间会略微增加。默认值是 1我一般在大区域reduceRegion或者对象化分类的时候直接调到 4出问题的概率能降一大截。还有一个特别隐蔽的坑reduceRegion的scale参数容易被忽略。如果你的影像集合分辨率是 10 米但reduceRegion没指定 scaleGEE 会按最小分辨率去算可能产生几亿个像元参与聚合内存瞬间爆掉。写代码的时候养成习惯只要是reduceRegion就显式加上scale、maxPixels和tileScale三个参数。最后一个经验是关于脚本结构的。我现在的习惯是把每个阶段的中间结果都导出成 Asset比如云掩膜后的集合、月合成集合、特征堆叠、训练好的分类器各存一份。好处是一旦后面某一步出错你不用从原始影像重算从最近的 Asset 接着跑就行。代价是 Asset 有容量限制一张 10 米分辨率的全市特征堆叠可能好几个 G需要定期清理。但相比每次调试都等四十分钟这点存储成本完全值得。提示调试阶段强烈建议用一个乡镇或者一个县的范围把roi缩小所有计算能在十几秒内出结果这样你一天能试几十次参数。确认方案之后再放大区域效率差好几倍。7. 几个容易被忽略但影响很大的细节前面讲的都是主干流程还有几个细节我觉得值得单独拎出来说因为它们不涉及核心技术但会在你不知不觉中把精度拉低。第一个是作物日历的本地化。不同纬度的同一作物物候期能差一两个月。北方冬小麦 3 月返青南方冬小麦可能 2 月就返青了北方一季稻 5 月插秧南方双季早稻 4 月就开始了。如果你直接套用文献里的特征窗口或者物候参数效果会打折扣。我的做法是先用一批已知类型的样本点把它们当年的 NDVI 曲线都画出来求平均看看不同作物的峰值到底落在哪个月再据此确定特征窗口。这个步骤花不了半小时但能让特征集的针对性提升一个档次。第二个是水体的干扰。水稻田在淹水期NDVI 会急剧下降如果只是机械地按 NDVI 阈值做筛选很容易把正在泡田的水稻判成水体或者裸土。这时候 LSWI 或者 NDWI 就派上用场了它能在 NDVI 低值期提供额外的判别信息。这也是为什么我在特征清单里坚持放水分指数尤其是在水稻种植区。第三个是多云多雨地区的观测数问题。西南和华南的一些区域一年有效观测可能只有二三十次有些月份一个干净观测都没有。这种情况下月合成会大面积出现空缺谐波拟合也会因为样本点太少而失真。我的应对方式是把合成窗口从月改成半月或者旬牺牲一点平滑度换取更多有效观测同时把 Sentinel-2 和 Landsat 合并使用虽然尺度不一致但在观测稀缺的情况下多一份数据比尺度一致性更重要。这个取舍没有标准答案得看你的项目对精度的要求。第四个是分类结果的年际一致性。单年结果看起来不错但把连续几年的图叠起来你可能会发现同一块地每年都在变这显然不符合农业生产的实际情况。解决办法是在多年度分类之间加一层时域一致性约束比如用隐马尔可夫模型或者简单的规则投票把三年里两年都是玉米的地块统一成玉米。这一步在国内的种植结构制图项目里越来越常见因为它直接关系到成果能不能支撑复种指数分析。我在实际使用中最深的一个体会是GEE 把工程门槛降下来之后真正决定成果质量的东西全部回到了遥感本身——数据源的物理特性、作物物候规律、样本的代表性。技术栈只是把试错成本从三天跑一遍降到四十分钟跑一遍让你有更多机会去反复验证这些本质问题。所以别把精力都花在调参和炫技上多花花时间去看每块地的 NDVI 曲线长什么样比什么都管用。