做全球森林损失量相关的研究或业务在 Google Earth EngineGEE里有个函数你大概率绕不开——stratifiedSample。我第一次见到它的时候以为它跟ee.FeatureCollection.randomPoints差不多无非是多一个“按类别抽”的选项。真正用它去抽 Hansen 全球森林变化数据里的样本点才发现门道比想象中多得多类别波段怎么构造、参数怎么配、为什么抽出来的点数不是期望值、全球范围为何老超时……这篇就把我在实际项目里用stratifiedSample提取分层样本点的完整过程、完整代码和踩坑记录一起整理出来方便后来人直接抄作业。1. 为什么森林损失量评估绕不开“分层”这两个字1.1 一个“翻车”案例简单随机抽样抽不中损失区先说一个我早期做森林损失验证时的真实经历。当时目标是在东南亚某国评估 2001-2020 年的森林损失面积我随手用ee.FeatureCollection.randomPoints在边界内撒了 2000 个点兴致勃勃地去看这些点落在 Hansen 数据的lossyear波段上是什么值。结果一统计98% 以上的点都是 0也就是说几乎全部落在“二十年间没发生过林冠损失”的像元上。这个结果并不意外。全球森林年损失率本身就很低通常只占森林总面积的零点几个百分点到几个百分点如果把损失年份再细分每一个单年份的比例就更小了。简单随机抽样在样本量固定的情况下抽到每个类别的点数基本等于“按面积比例碰运气”稀有类别很容易被抽成个位数甚至为 0。样本量一少后面无论是做面积估算的精度验证还是训练一个分类器去识别损失年份都会变得非常不稳定。后来我在项目里换成了stratifiedSample把无损失区和各个损失年份分别作为“层strata”每一层都抽取固定的点数。同样 2000 个样本损失区每个年份段都能拿到几百个点。训练和验证结果完全不一样。1.2 分层抽样的本质先给每个“稀有类别”发一张入场券很多人把分层抽样理解成“按类别随机抽点”这只是表象。它真正的核心是把研究区按照目标变量划分成若干个互不重叠的类别然后在每个类别内独立地进行随机抽样并且可以控制每个类别的样本量。放到森林损失场景里目标变量就是“2001 到 2020 年间森林是否损失、以及在哪一年损失”。我们先把影像变成一张整数值的“分层标签图”例如 0 表示无损失、1 表示 2001-2005 年损失、2 表示 2006-2010 年损失然后在每个标签上分别抽固定数量的点。这样一来哪怕某个损失时段只占研究区总面积的 1%只要它还有 1000 个像元你也能在它上面抽到足够多的样本点。这个过程带来的直接收益有两个训练样本的类别结构不再被面积比例绑架稀有类别也能拥有足够样本量如果后续要做分层面积估算比如用 Olofsson 2014 那套无偏估计器分层样本可以提供更小的方差比同样本量的简单随机抽样更高效。理解这层逻辑之后再看 GEE 的stratifiedSample函数你就会明白它为什么“按类别抽”而不是“按面积比例抽”。2. Hansen 全球森林变化产品底层数据口径必须吃透2.1 核心波段逐个拆解题目里的“全球森林损失量”在 GEE 中最常用的数据源就是 Hansen 团队发布的全球森林变化数据集。在代码里加载是这样var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8);这个影像的波段比较多真正做分层抽样时常用的只有这几个波段名含义取值范围treecover20002000 年的林冠覆盖度百分比0-100loss2001-2020 年间是否发生林冠损失0 或 1lossyear损失发生的年份编码0-20值 1 对应 2001 年值 20 对应 2020 年gain2000-2012 年间是否发生林冠增加0 或 1datamask数据有效掩膜0 表示无数据1 表示有效观测几个容易搞混的点先说清楚第一lossyear里的 0 不代表“2010 年”之类的年份它代表“2001-2020 年间在该像元上没有探测到林冠损失”。所以不能直接拿lossyear当时间轴使用。第二treecover2000是百分比不是树的数量。业界做森林损失相关研究时普遍把treecover2000 10作为森林掩膜的常用门槛但你也可以按自己的研究需求调整成 20、30没有绝对标准。第三datamask和treecover2000的掩膜要分清。datamask是数据本身是否有有效观测treecover2000的阈值筛掉的是“不是林地区域”的像元。分层抽样的类标签最好只建立在“有效观测且是森林”的像元上否则你可能会抽到大量无数据或非林地区域的点。2.2 从lossyear到“分层标签”的合并方案明确了数据口径之后就要开始做分层标签波段了。这里有一个关键决策是直接使用原始lossyear作为分层标签还是把它合并成粗粒度的时间段我两种方式都用过直接说结论局部区域、样本需求精细直接使用lossyear原始值0-20 共 21 类每个年份单独抽点。优点是可以逐年做变化检测缺点是类别数多同样的目标总样本量会被摊薄如果某些年份面积特别小还可能出现某类像元不足导致抽不满的情况。全球或大洲尺度、业务需求不关心单年份把lossyear合并成几个时间段例如无损失、2001-2005、2006-2010、2011-2015、2016-2020共 5 类。缺点是损失时间精度变粗优点是类别少、每个类样本量足、运行稳定。我的经验是做全球森林损失量这种尺度非常大的题目时优先用合并时间段方案。20 个年份类看起来信息更全但全球范围下stratifiedSample的运算量和对每个类别的像素要求都更高很容易超时或内存溢出。先把分层跑通再考虑要不要细化。下面这段代码是把lossyear构建成时间段的经典写法var lossYear hansen.select(lossyear); var treeCover hansen.select(treecover2000); // 森林定义2000 年树冠覆盖度不低于 10% var forestMask treeCover.gte(10); // 分层标签0无损失12001-200522006-201032011-201542016-2020 var strata ee.Image(0) .where(lossYear.gte(1).and(lossYear.lte(5)), 1) .where(lossYear.gte(6).and(lossYear.lte(10)), 2) .where(lossYear.gte(11).and(lossYear.lte(15)), 3) .where(lossYear.gte(16).and(lossYear.lte(20)), 4) .updateMask(forestMask) .toInt();这里的关键点在于先用ee.Image(0)打底再逐层用where赋值。这样做可以保证无损失类标签 0在森林掩膜内一直存在而不会因为原始lossyear在无林地区也有 0 而把标签搞混。最后一定要.toInt()因为stratifiedSample要求类别波段最好是整型。顺便说一句如果你确实想用逐年标签代码反而更简单var strataByYear lossYear .updateMask(forestMask) .toInt();但是你要做好心理准备0-20 共 21 个类别当numPoints设置为 50 时返回的样本总量会是 21 × 50 的级别后面我会专门讲这个数量关系。3.stratifiedSample逐参数拆解与完整可跑代码3.1 API 参数逐项说明ee.Image.stratifiedSample()是 GEE 在 Image 对象上的方法接收一个参数字典。我实际项目里最常用到的参数如下参数作用我的建议numPoints每个分层类别要抽取的目标点数先设为 50-200 跑通再按任务调整classBand存放分层标签的波段名一定用整型波段值域为 0-20 这种离散整数region抽样范围必须显式指定否则可能按影像全球范围计算scale采样尺度米全球用 500-1000国家用 100局部用 30seed随机种子固定数值可复现不设则每次结果不同geometries是否返回点几何要看后续用途决定后面细讲这里特别提醒一个容易误解的地方scale不是“点与点之间的最小距离”而是“抽样的像元尺度”。你可以把它理解为在多大的网格上逐像元判断类别值。把scale从 30 改成 1000不一定能直接控制抽样点的密度但它会显著影响参与计算的像元数量和运行速度。全球范围跑不动时第一步就是把scale调大。还有seed我建议无论是在写文档还是写实验记录时都固定下来。分层抽样用于精度验证时一个确定的随机种子意味着别人可以复现你的样本点而做训练/验证样本分离时你可以用两个不同的seed分别跑两次产生两套互不重叠的样本集。3.2 完整代码以柬埔寨森林损失分层抽样为例以下代码我以一个东南亚国家为例完整演示从加载 Hansen 数据到生成分层样本点的全过程。你可以把roi换成自己的研究区// 1. 研究区柬埔寨也可以直接用全球多边形但更建议分区域处理 var roi ee.FeatureCollection(FAO/GAUL/2015/level0) .filter(ee.Filter.eq(ADM0_NAME, Cambodia)) .geometry(); // 2. 加载 Hansen 数据 var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8); var lossYear hansen.select(lossyear); var treeCover hansen.select(treecover2000); // 3. 构建分层标签波段 var forestMask treeCover.gte(10); var strata ee.Image(0) .where(lossYear.gte(1).and(lossYear.lte(5)), 1) .where(lossYear.gte(6).and(lossYear.lte(10)), 2) .where(lossYear.gte(11).and(lossYear.lte(15)), 3) .where(lossYear.gte(16).and(lossYear.lte(20)), 4) .updateMask(forestMask) .toInt(); // 可视化分层结果 Map.centerObject(roi, 7); Map.addLayer(strata, {min: 0, max: 4, palette: [green, yellow, orange, red, darkred]}, strata); // 4. 核心分层抽样 var samples strata.stratifiedSample({ numPoints: 100, classBand: strata, region: roi, scale: 100, seed: 42, geometries: true }); // 5. 查看每个类实际抽到了多少点 var counts samples.reduceColumns({ reducer: ee.Reducer.count().group(1, strata), selectors: [strata] }); print(counts); // 6. 加载到地图上观察 Map.addLayer(samples, {color: black}, stratified samples);直接把这段代码复制到 GEE 编辑器里只要网络没问题基本可以一次跑通。打印出来的counts会按strata值分组显示每个类别的实际点数正常情况下 5 个类别都应该是 100 上下。某些类别如果像元数极少实际点数会少于 100这是正常的后面会解释。3.3 结果检查抽样数量为什么和预期不一样刚跑完代码很多人会发现一个现象打印出来的实际样本量并不是简单的numPoints也可能不是恰好 5 × 100 500。第一个原因numPoints是“每个类别的目标点数”不是总点数。类别数乘目标点数才是理论总样本量。5 个类、每个类 100 点理论总量 500 点这没错。但如果你直接拿lossyear原始值当类别波段那就是 21 个类别numPoints 100时理论总量就是 2100 点。这个乘数关系很容易被忽略。第二个原因某些类别的像元数不足。比如一个典型的 2001-2005 年损失时段在大面积国家尺度通常够用但如果你把区域缩小到某个小流域该类别的总像元可能只有 30 个那么 GEE 只抽 30 个不会凭空造点。第三个原因掩膜没有处理好。如果forestMask把无损失区也误掩了类别 0 的像元可能很少甚至为 0。所以在构建strata时一定要先想清楚“每一层到底有哪些像元能参与抽样”。为了下一步能放心我习惯把分层结果和样本点叠在地图上逐个类别目视检查一遍看看点是否落在了应该落的位置。机器不会骗人但逻辑没写对时机器也会老老实实地把错的结果跑给你看。4. 实战中的坑规模爆炸、类别不足、坐标问题逐个击破4.1 全球范围跑不动scale和分块策略标题里写了“全球森林损失量”所以我猜不少读者最终是想在全球尺度上抽取样本点。这里必须提前打预防针在全球范围直接对一个 30 米分辨率的影像做stratifiedSample哪怕类别只有 5 个也非常容易遇到Computation timed out或User memory limit exceeded。我的解决办法是分级降尺度加分区域。具体来说把scale从 30 米调到 500 米甚至 1000 米。分层抽样不要求点必须落在原始 30 米像元中心500 米尺度下抽出来的点照样可以作为区域级样本。按国家或大洲分块执行。全球多边形不是不能做但一次性把所有陆地区域都塞进一个任务里效率和稳定性都会很差。我通常是按 FAO 的 GAUL 行政边界拆成几十个国家循环提交导出任务最后在本地合并。利用导出任务的重试机制。GEE 导出任务偶发失败很正常尤其是大尺度任务。设置好seed之后失败重跑也能得到相同结果这恰恰是固定seed的另一个好处。代码层面分块抽样可以这样设计var countries ee.FeatureCollection(FAO/GAUL/2015/level0); // 只演示前 10 个国家的循环 var nationList countries.limit(10).toList(10); for (var i 0; i 10; i) { var nation ee.Feature(nationList.get(i)); var samples strata.stratifiedSample({ numPoints: 100, classBand: strata, region: nation.geometry(), scale: 500, seed: 100 i, geometries: true }); Export.table.toDrive({ collection: samples, description: stratified_samples_ i, folder: gee_export, fileFormat: SHP }); }注意这个循环是把多个导出任务提交到 GEE 任务列表里而不是在代码里同步等待。这也是 GEE 的推荐思路——计算在服务端执行客户端只需要排队。4.2 类别波段的数据类型整型 or 浮点结果差很多stratifiedSample是按classBand的像素值进行分组计数的。如果类别波段是浮点型比如lossyear在某些派生数据里变成了浮点那么 5.0 和 5.9999 可能被当作两个不同类别导致类别数量爆炸、每个类别的像元数被严重稀释。这通常是“抽样结果莫名其妙”的头号原因。所以我的习惯是无论原始波段是什么类型在传给stratifiedSample之前一律显式.toInt()。如果你担心取整导致标签偏移可以先.round().toInt()不要省略这一步。4.3geometries: true与geometries: false的选择这个参数很多人不重视但在后续使用中影响很大。geometries: true表示返回的每个样本是一个带点几何的Feature。适合你在本地直接Map.addLayer目视检查或者导出成 SHP、GeoJSON 继续做外业验证。geometries: false返回的样本点不包含真实点几何longitude、latitude会以属性字段的形式存在记录里。这样做的最大好处是大幅减小数据体积。当你在全球范围抽取几万个点、准备导出成 CSV 时用false会让导出和后续处理轻松很多。一个实用经验先把geometries设为true跑一个小范围确认类别逻辑没问题正式大规模抽点时切成false把输出结果压缩成表格需要坐标时再从属性里取。4.4 与randomPoints的分工不要混用ee.FeatureCollection.randomPoints是另一个高频抽样函数但它是完全的均匀随机撒点没有任何类别配额。如果你只是随机撒点做土地利用覆盖验证randomPoints够用但只要你关心“每个损失类别都要有足够样本”就必须用stratifiedSample。我在项目里还见过一种错误做法先用randomPoints撒了几万个点然后做空间连接去提取lossyear值试图事后筛选每个类别。这样不是不行但样本量效率极低尤其是稀有类别几万个随机点筛下来可能只剩几十个。与其事后补救不如一开始就用分层抽样按类别配好额度。5. 验证与进阶从样本点到森林损失面积估算5.1 分层代表性的检查看一眼类别分布抽样跑完不能直接信任至少要做一个分布检查。前面代码里的counts只能告诉你“每个类别抽了多少点”还不够。我还习惯把分层结果里的类别面积占比也算出来和样本量占比做一个对比确保分层抽样的样本结构符合预期。// 计算每个分层类别的面积占比 var areaImage strata .multiply(ee.Image.pixelArea()) .reduceRegion({ reducer: ee.Reducer.sum().group(1, strata), geometry: roi, scale: 100, maxPixels: 1e10 }); print(areaImage);把areaImage和样本点数量的分组结果放在一起看你就能直观地感受到一个事实类别 0无损失的面积可能占了 95%但它的样本量只有 1/5这正是分层抽样的设计意图——牺牲无损失区的冗余样本换取损失区的有效样本。5.2 从样本到森林损失面积估算的链路题目里的“全球森林损失量”最终往往要落到面积估算上。分层样本在这条链路中的关键作用是提供“类别比例”的无偏估计。你把每个分层类别的样本点落回原始影像统计每个类别的样本比例再乘上该类别的实际面积就能得到损失面积的分层估计值。公式不复杂第 h 类的面积A_h由该类别的像元数乘像元面积得到第 h 类面积占总面积的估计比例p_h用该类的样本比例代替总体损失面积估计为各时段类别面积比例加权求和。这个过程中stratifiedSample负责提供每一类有代表性的样本比例。抽样时需要注意的一点是如果某个类别的样本量太少面积比例估计的置信区间会很宽。所以我在实操中都会把重点关注类别的numPoints调得高一点牺牲一点无损失区的样本量。5.3 这些样本如何继续用于分类器训练除了面积估算分层样本点还有一个高频用途——作为监督分类的训练样本。抽完点之后通常要这样处理把样本点与多波段遥感影像做采样得到每个点在各个光谱波段上的值用样本点的strata值作为标签没有损失为 0损失年份时段为 1-4将特征和标签喂给随机森林、CART 或深度学习模型。GEE 里做第一步最常用sampleRegionsvar trainingImage hansen .select([treecover2000, loss, gain]) .addBands(strata.select(strata).rename(label)); var trainingData trainingImage.sampleRegions({ collection: samples, scale: 30, geometries: true }); // 拆分训练集和验证集 var split trainingData.randomColumn(split, 123); var train split.filter(ee.Filter.lt(split, 0.7)); var test split.filter(ee.Filter.gte(split, 0.7));这里的randomColumn相当于又一次抽样但它是在已有的分层样本内做二次切分不会改变每个类别的样本配额结构。如果你的分层样本量很紧张我更建议用两个不同seed分别生成训练样本集和验证样本集而不是事后切分。5.4 时间维度的延伸如果要做逐年变化检测前面提到过合并时间段的方案更适合全球尺度粗粒度评估。但如果你确实要分辨每一年损失比如分析某年森林火灾的突然集中爆发那你需要回到逐年方案。这时的样本量需求会急剧上升因为 21 个类别中面积小且像元少的年份类很容易抽不足。我的做法是先做一个快速估算数一数每个lossyear类在目标区域的像元数再决定numPoints。如果某个年份像元数特别多numPoints可以适量加大像元数特别少的年份接受它“抽不满”的现实后续在模型里可以对该年份类做类别加权或过采样。分层抽样解决的是“每种类型至少有机会被抽中”它解决不了“某种类型在地理上本身就极少”的问题这一点必须心里有数。最后再提醒一个细节记录好seed和所有构建类别的阈值不只为了论文复现。分层样本点一旦导出后续外业调查或人工目视判读可能要花掉大量时间没有种子和阈值记录重新生成一套“一模一样”的样本会非常麻烦。我在实际项目中吃过这个亏现在每一次跑stratifiedSample都会在注释里写清楚用的是什么 Hansen 版本、什么森林阈值、什么随机种子这比任何花哨的参数调优都更救命。