尧图网页设计 
干我们这行的都知道GEE里坑最多的不是复杂的波段计算反而是样本提取这种“看着很简单”的环节。特别是做全球尺度的森林损失研究你不可能人工去逐点标注几十万个验证样点也没法用简单的随机抽样指望稀疏的损失像元能被抽中。去年我做全球森林变化产品精度评估时被stratifiedSample这个函数卡了整整三天反复研究后才发现问题往往不在函数本身而在你对“分层”这件事的理解。这篇文章就把我从数据准备、类别设计到分层样本提取、质量检查的全过程完整拆出来每一步都给可直接复用的代码和参数希望能让你少走几趟弯路。这个内容的核心是基于Hansen全球森林变化数据集用GEE内置的stratifiedSample分层抽样方法在“森林损失”这个类别分布极不均衡的目标上科学地提取分层样本点。它解决的痛点非常明确——如何在全球动辄几十亿像元的影像中让稀有类别比如某一年发生损失的像元也能被充分抽样同时保持样本的空间代表性和时间完整性。无论你是做森林变化监测、土地覆被产品精度验证还是想学GEE分层抽样的原理这套方案都值得完整看一遍。1. 这个需求背后到底在解决什么问题1.1 为什么偏偏选Hansen全球森林变化数据做森林损失量分析第一步就是选底图。GEE里不缺森林相关的数据集有MODIS的年度森林盖度、有ESA的全球土地覆被还有各种局部区域的林龄产品。但Hansen全球森林变化数据集UMD/hansen/global_forest_change几乎成了默认选项原因无非三点。第一它的空间分辨率是30米而MODIS那种250米到500米的产品对森林损失这种“碎片化”事件来说太粗了。一条道路建设或一次小规模皆伐在MODIS影像上可能只体现为一个混合像元根本没法数清损失了多少。Hansen的30米分辨率可以捕捉到大多数非法的、小规模的森林采伐。第二它提供了从2000年基线到逐年的损失年份lossyear信息。这个“逐年”属性特别关键——你不仅能知道哪里发生了损失还能知道是哪一年损失的这让“全球森林损失量”的分析从静态制图变成了时间序列分析。我在做数据产品精度评价时最看重的就是这一点。第三它在GEE中是以固定资产形式托管的访问速度快、不用自己做复杂的预处理。你要做的只是select出需要的波段剩下的几何校正、大气校正全都不必操心。但这里有个极其重要的实操细节很多人忽略了Hansen的lossyear波段里0代表“2000年至数据版本年份之间未检测到损失”1到数据版本年份才代表对应的损失年份。当你把它直接用于分层抽样时如果不做任何处理0类像元可能占据99%以上的面积你抽取的样点会被无损失像元淹没损失年份样本根本抽不到几条。这正是需要引入“分层”概念的直接原因。1.2 分层抽样 vs 简单随机抽样差在哪说句实话如果目标类别分布是均匀的stratifiedSample和普通的sample没什么太大区别。但在森林损失这个场景下二者差异是决定性的。简单随机抽样是在整个目标区域内完全随机地撒点类别分布和像元分布保持一致。问题在于如果我们研究的区域里只有0.5%的像元属于“2015年发生损失”那么随机抽10000个点理论上抽中2015年损失层的只有50个这50个点根本不足以支撑任何统计分析。分层抽样stratified sampling的思路则是先把研究区按类别规则切成若干个互不重叠的“层”strata然后在每一层内部独立地随机抽样。这样无论2015年损失层有多么稀疏只要你指定给它抽200个点它就一定会给你抽出200个。这就保证了每个类别都有“最低样本保障”同时每一层的样本又保持了随机性。用一个生活类比来说如果你想了解整个城市不同收入人群的消费习惯简单随机抽样就好比从电话簿上随机挑人结果高收入人群可能一个都没被抽到分层抽样则是先把市民按收入分成高、中、低三个名单再从每个名单里分别随机挑人。你所付出的额外成本是必须先知道“怎么分层”但换来的回报是每个层级都有足够的样本数量。在GEE里stratifiedSample函数做的就是这件事。它允许你指定一个分类波段作为分层的依据每一个独特的整数值视为一层然后在每一层内部执行随机采样。它内部会高效地做分桶统计和随机定位比你自己写reduceRegions然后拼接多个filter结果要可靠得多。2. stratifiedSample的核心机制与参数门道2.1 函数签名与每个参数的真实含义直接看GEE里stratifiedSample的调用方式API文档写得比较简洁实际含义得靠用才能体会。基本形式如下var samples image.stratifiedSample({ numPoints: 100, classBand: class, region: roi, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3] });每个参数的真实含义numPoints是每一层内计划抽取的样本点数。注意是“每一层”都抽这个数不是总共抽这个数。如果你的分类波段有5个类别值则总样本数大约是5乘以numPoints实际可能因某个类别像元不足而少于该值。classBand是分层依据的波段名称。它必须是整数类型Int或Byte因为分层就是按整数值归类。如果波段是浮点型你需要先重编码成整数否则函数会报错或者行为不符合预期。region定义了采样范围。可以是一个ee.Geometry、ee.FeatureCollection或ee.Feature。不传的话默认用影像全幅范围但全球30米分辨率逐类采样计算压力极大强烈建议传区域。scale是采样时考虑的空间分辨率。这个值的取舍很值得讲。用30米意味着逐像元级别的采样最精确但对计算资源的需求高用1000米意味着采样的“定位”粒度较粗适合全球快速预览。我的经验是如果研究区是一个国家或大洲级别且你的计算配额不算富余先用250米跑通全流程最后对关键区域再用30米做一次精细抽样。geometries是一个布尔值决定输出的FeatureCollection中每个样本是否携带点几何。如果设为false返回的记录只有属性值没有空间坐标设为true则有GeoJSON点。你做精度验证、人工目视判读或后续叠加高分辨率影像都必须设为true。seed是随机种子。这个参数的意义超过多数人的想象。分层抽样使用了随机数生成器如果不固定种子每次运行生成的样本位置都不同论文和报告的复现性就没法保证。在学术研究中审稿人要求“样本结果可复现”时你只需要亮出固定的seed即可。我建议用日期加版本号比如20240601清晰又不容易冲突。classValues是可选参数用于显式指定要采样的类别值列表。默认会对classBand中的全部类别分别抽样。如果你只关心特定年份的损失层这个参数可以避免多余的计算。2.2 分层变量该怎么构造重新编码成类别波段这是整个流程中最容易被低估的一步。很多人直接把Hansen的lossyear波段塞给classBand然后发现样本数量巨大、运算超时或者返回的结果里损失年份样本出奇地少。原因在于lossyear的取值并不适合直接作为分层依据。Hansen的lossyear取值是0到20的整数但在某些版本和特定区域2000年之前已经有森林覆盖但从未经历损失的大片区域全是0值。直接把0到20每个值当作一层会产生两个问题一是层数太多导致numPoints乘以层数的总量非常惊人二是0层占绝对优势你要重点关注的不同时期损失层被稀释在了过细的分层结构里。解决思路不是放弃lossyear而是在它基础上重新编码成一个“有意义的类别波段”。比如按时间窗口合并类别将森林损失划分为“未损失”、“2001-2005年损失”、“2006-2010年损失”、“2011-2015年损失”、“2016-2020年损失”五类。这样既保留了时间维度上的层次又不会造成层数过多。代码实现如下// 加载Hansen全球森林变化数据集 var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8); // 选择关键波段 var treeCover hansen.select(treecover2000); var lossYear hansen.select(lossyear); // 构建分层类别波段 // 0 未损失森林 (treecover2000 30 且 lossyear 0) // 1 2001-2005年间损失 // 2 2006-2010年间损失 // 3 2011-2015年间损失 // 4 2016-2020年间损失 var forestMask treeCover.gte(30); var classBand ee.Image(0) .where(forestMask.and(lossYear.gte(1).and(lossYear.lte(5))), 1) .where(forestMask.and(lossYear.gte(6).and(lossYear.lte(10))), 2) .where(forestMask.and(lossYear.gte(11).and(lossYear.lte(15))), 3) .where(forestMask.and(lossYear.gte(16)), 4) .clip(roi) .byte();你可能注意到了我把森林定义成treecover2000大于等于30%。这不是拍脑袋定的。Hansen团队官方建议森林/非森林的阈值常用30%这也是许多联合国报告和森林资源评估采用的标准。当然你完全可以根据研究需要改写成10%或50%但必须保持一致并在论文或报告里明确声明阈值。2.3 种子、尺度和几何输出的设计细节在GEE里使用stratifiedSample时还有几个容易踩的细节提前说清楚能省很多调试时间。关于scale的选择逻辑我强烈建议你做一个“两级采样”设计。第一级在全球或大区域尺度用scale参数设置为1000米甚至5000米跑一次分层抽样目的是验证分类波段是否正确、各类别样本空间分布是否合理、总体样本量是否符合预期。这一步成本极低一两分钟就能出结果。第二级在确认分类没问题后再对重点区域用scale: 30做精细抽样。这样可以把90%的调试时间花在低成本的第一级上而不会因为一个低级错误浪费掉几十分钟的计算配额。关于geometries参数还有一个小坑。当你设置geometries: true时每个样本点都带空间几何这个集合如果在GEE的Map面板上直接显示可能因为点数太多导致页面卡顿。建议先用limit或者按类别filter后再显示或者干脆只在导出Data后本地查看。另一个容易被忽略的点是stratifiedSample返回的FeatureCollection不仅包含classBand指定的类别值还会包含参与分层影像的所有波段信息。如果你的分层影像是一个包含多波段的Image比如你把treecover、lossyear都addBands进去了那么每个样本都会带上这些波段的像元值作为属性。这其实是一件好事——你在后续精度验证或者建模时不需要再次去影像里提取属性因为所有属性已经挂在样本点上了。但要注意如果影像中包含了不必要的波段比如原始Hansen数据中的数据质量波段输出的属性表会变得臃肿导出CSV后光字段理解就要花不少时间。所以采样前用select限定波段是个好习惯。3. 全球森林损失样本提取的完整实操3.1 准备数据与目标区域先明确我们要做什么。以东南亚地区为例我想提取覆盖五类森林状态未损失、四个时期损失的样本点每个类别各抽取300个用于后续的高分辨率影像目视判读和损失面积估算。东南亚是森林损失重灾区但也是行政边界复杂、影像噪声较多的地区用它来演示更有说服力。数据准备代码// 行政边界导入可以用FAO的GAUL数据集或全球行政区划数据集 var roi ee.FeatureCollection(FAO/GAUL/2015/level0) .filter(ee.Filter.inList(ADM0_NAME, [Malaysia, Indonesia, Thailand])) .geometry(); // 加载Hansen数据 var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8); var treeCover hansen.select(treecover2000); var lossYear hansen.select(lossyear); Map.centerObject(roi, 5); Map.addLayer(treeCover.updateMask(treeCover.gt(0)), {palette: [#006400, #00ff00]}, Tree Cover 2000);森林盖度图层添加到地图上目的不是看热闹而是确认目标区域内森林的空间分布是否符合预期。比如东南亚的婆罗洲、苏门答腊这些热点区域应该有明显的森林覆盖如果发现城市中心也有大量森林盖度说明数据或区域投影有问题要及早排查。3.2 构造分层类别把连续损失变成离散层接着用之前提到的思路重编码分层类别并把需要携带的属性波段一起放在一个影像里// 构建森林区域掩膜 var forestMask treeCover.gte(30); // 构建分层类别波段 var classBand ee.Image(0) .where(forestMask.and(lossYear.eq(0)), 0) .where(forestMask.and(lossYear.gte(1).and(lossYear.lte(5))), 1) .where(forestMask.and(lossYear.gte(6).and(lossYear.lte(10))), 2) .where(forestMask.and(lossYear.gte(11).and(lossYear.lte(15))), 3) .where(forestMask.and(lossYear.gte(16)), 4) .clip(roi); // 创建用于采样的多波段影像 // classCode是分层依据lossYear和treeCover作为样本属性一并携带 var sampleImage ee.Image.cat([ classBand.rename(classCode), lossYear.rename(lossYear), treeCover.rename(treeCover) ]).byte(); // 给classCode波段起名为class因为stratifiedSample默认以该波段分层 var sampleImage sampleImage.select( [classCode, lossYear, treeCover], [class, lossYear, treeCover] );这里有两个细节。第一我用ee.Image.cat把多个波段拼接成一个多波段影像然后重命名。classCode重命名为class这样stratifiedSample里classBand: class一目了然。第二所有波段最后都转成了byte8位无符号整数。Hansen的lossyear和treecover取值都在0到100之间byte足够存储减少一半以上的内存和网络传输消耗。这在全球尺度的计算中是非常可观的优化。3.3 执行stratifiedSample并导出检查核心语句来了// 设置抽样参数 var numPointsPerClass 300; // 执行分层抽样 var samples sampleImage.stratifiedSample({ numPoints: numPointsPerClass, classBand: class, region: roi, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3, 4] }); print(Sample FeatureCollection, samples); print(Size of samples, samples.size());运行后在Console里观察如果一切正常输出集合大小应该在1500附近5类乘以300但实际情况可能因为某些类别在区域内像元不足而少于1500。比如2001-2005年东南亚的损失相对较少抽取300个点时可能只能抽到280-290个这是正常的说明该类别在给定区域内的有效像元数不够。为了快速检查各类别实际抽到了多少点可以用下面的代码分组统计// 按类别统计样本数量 var classCounts samples.reduceColumns({ reducer: ee.Reducer.frequencyHistogram(), selectors: [class] }); print(Class counts, classCounts);如果某个类别抽到0个样本而你的研究又必须包含该类别就需要放宽region范围或者下调numPoints之外的硬约束。比如可以用多期Hansen版本拼接数据。但多数情况下出现某个类别完全没有样本的情况几乎总是因为分层波段构建逻辑有误。3.4 样本可视化与快速质检拿到样本后第一件事不要急着导出先在地图上可视化检查空间分布。我习惯的配色方案是未损失层用深绿色四个损失时期分别用橙黄红紫四色这样一眼就能看出是否均匀分布在校准框内。// 将样本按类别分割开显示 var palette [#006400, #ff9900, #ff3300, #ff00ff, #0000ff]; Map.addLayer(samples, {color: 000000}, All samples (black));更实用的做法是把样本属性表导出为CSV用统计软件快速查看每类的空间坐标范围、lossyear分布直方图。具体导出代码如下// 导出样本点 Export.table.toDrive({ collection: samples, description: SEAsia_forest_loss_samples, fileFormat: CSV, folder: GEE_exports, selectors: [class, lossYear, treeCover, .geo] });导出后我习惯先在Excel里打开确认几个关键字段class的取值是否只在0到4之间lossYear的取值是否和class相互匹配比如class1时lossYear应该在1到5之间检查是否存在坐标为(0,0)或位于研究区之外的异常点。这一步虽然简单但往往能立刻暴露问题。我在一次实验中发现class42016-2020损失的点里混进了lossYear0的记录追根溯源发现是lossyear构建时分类波段和损失年份的掩膜顺序写反了重新调整处理逻辑后一切正常。4. 实战中一定会遇到的坑4.1 “每个类别的像元数不够”怎么解stratifiedSample最让人头疼的运行报错之一是每个类别要求抽取numPoints个点但某些层实际像元数少于numPoints。常见的提示是Number of pixels in strata number of requested samples。这事我遇到过太多次尤其是当研究区小、且某个损失年份特别稀疏时。解决办法有几种按优先级排列第一种调低numPoints让每个类别的期望样本数小于最少类别的有效像元数。你可以先用reduceRegion或classBand.reduceFrequencyHistogram快速统计各层像元数量再据此设定numPoints。第二种扩大region让稀有类别获得更大的空间范围来补充有效像元。比如把研究区从单一州郡扩大到整个省。第三种调整scale。scale从30米改成100米或250米意味着每个“采样像元”覆盖更大的地面范围类别有效像元数也随之增多。代价是采样点的定位精度下降所以它更适合全球尺度的粗略研究不适合需要精细点位的研究。第四种如果实在无法满足就拆分区域。把研究区分成几个子区域分别执行分层抽样最后用merge合并所有样本。这个方法多花一点代码和时间但能保留30米的分辨率也是我处理非洲多国联合分析时最常用的一招。// 拆分区域后分别抽样再合并 var subregions [regionA, regionB, regionC]; // 每个子区域都要有对应几何 var allSamples ee.FeatureCollection([]); subregions.forEach(function(name) { var subRegion ee.FeatureCollection(rois).filter(ee.Filter.eq(name, name)).geometry(); var subSamples sampleImage.stratifiedSample({ numPoints: numPointsPerClass, classBand: class, region: subRegion, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3, 4] }); allSamples allSamples.merge(subSamples); });需要注意的是GEE的forEach不是同步执行的合并时最好用ee.FeatureCollection的flatten方式或者在ui模式下改用map。上面的代码只是示意实际使用推荐写成函数式合并。4.2 采样点太密内存溢出与地图卡顿一次性生成几千个带几何的样本点在GEE的Map面板上直接显示十有八九会卡顿甚至白屏。这不是程序出bug而是客户端渲染压力过大。我的处理方式是把样本点写到自己的资产里或者导出GeoJSON到本地ArcGIS/QGIS查看。如果必须在GEE里预览部分样本可以这样// 显示每类前50个点共250个点足够预览 var previewSamples samples.limit(250); Map.addLayer(previewSamples, {color: red}, Preview samples (250));另一个内存陷阱是想把样本点和高分辨率遥感影像叠加做目视判读时直接addLayer高分辨率影像同样会让浏览器崩溃。解决方案是使用getThumbURL生成缩略图或者只给单个样本点周围创建一个小缓冲区域再显示。实际操作中我更倾向于用GEE Code Editor的Export.map.toDrive把注意力转到导出任务不在线上做深度交互。4.3 不同版本的Hansen数据波段差异GEE上的Hansen数据集更新过多个版本有v1_7、v1_8、v1_10等。版本之间虽然核心波段设计一致但也有细微差异。比如较新版本增加了datamask、gain等波段损失年份的最大值也不同v1_8损失年份到2020年v1_10则可以到2022年甚至更晚。这意味着如果你在前面代码里写死了lossYear.gte(16)在v1_10中会覆盖所有2020年之后的损失得到的结果维度完全变了。我的建议是研究开始前明确使用哪个版本并检查对应版本的最大损失年份。代码里用动态方式获取版本上限比写死数值更健壮var lossYear hansen.select(lossyear); var maxLossYear lossYear.reduceRegion({ reducer: ee.Reducer.max(), geometry: roi, scale: 1000, maxPixels: 1e9 }).getNumber(lossyear); print(Max loss year in ROI:, maxLossYear);这条命令可以快速确认你所选择的数据版本在目标区域内最晚的损失年份。许多同学跑完整套流程后才发现版本内损失年份的上限和预期不符不得不重算一遍费时费力。在第一步就确认版本和最大年份是一个好习惯。4.4 常见问题速查表为了更直观我把使用stratifiedSample时常见的坑整理成一张速查表建议收藏常见问题可能原因解决方案报错No valid pixels in strataclassBand在ROI内全为掩膜值检查classBand的掩膜范围确保每个类别在ROI内有有效像元某个类别样本数为0类别编码冲突或区域太小检查recode条件扩大ROI或调低scale两个类别的样本点重合分层波段取值定义重叠检查.where条件之间是否有边界重叠逐层用mask验证导出CSV字段过多未限定采样波段采样前做select只保留必要波段每次运行结果不同未设置seed固定seed确保结果可复现Map显示卡死样本点数过多用limit预览或导出本地查看损失年份最大年份异常使用了不同版本Hansen动态读取maxLossYear并检查数据集ID5. 样本拿到之后还能干什么5.1 用样本做面积估算与精度验证分层样本的价值绝不只是“抽出来了就完事”。在森林损失研究中最常见的后续工作是面积估算和分类精度验证。以面积估算为例Olofsson等人提出的良好实践指南里推荐使用分层抽样stratified random sampling得到的样本结合每层的像元数权重来修正面积估计。如果你把整个区域按“未损失”和“损失”两层抽样那么损失面积的无偏估计可以写成损失面积估计 各类别损失像元比例乘以各类别总面积再对各层加权合并。具体公式为A_hat Σ (W_h × A_total × p_h)其中W_h是第h层像元占比p_h是第h层中样本里判定为“损失”的比例A_total是整个区域总面积。通过这样的计算你得到的面积估计值比单纯的“损失像元乘以30×30平方米”要更稳健尤其适合样本量不大、类别不平衡的数据。在GEE中你可以先用reduceRegion拿到每个类别的像元数量和面积权重再用样本点和高分辨率影像人工判读的结果去更新p_h。这两个东西组合起来就能算出带置信区间的损失面积。具体代码实现较多不在这里展开但你要知道stratifiedSample抽出的样本是衔接“遥感分类结果”和“地面真实面积”的桥梁。另外如果你后续要训练一个自定义分类器比如区分火灾、砍伐、风暴导致的损失stratifiedSample产出的样本点还可以作为训练样本。做法很简单把样本点叠加到高分辨率影像切片上人工标注出损失原因再把标注后的点集灌回分类器。这时候你在分层抽样时保留的lossYear和treeCover属性就派上了大用场——它们可以作为辅助特征一起参与模型训练。5.2 结合高分辨率影像做人工审核的流程分层样本落到具体使用中尤其是做精度评估时光靠Hansen的30米分辨率是不够的你需要用更高分辨率的影像Planet、Maxar、Sentinel-2的10米来确认每个点真实的地表状态。实际操作中我的流程是这样的先把stratifiedSample导出的样本在GEE里和Sentinel-2无云影像合成叠加然后导出为PNG切片再放到本地影像管理软件中人工判读。更高效的做法是直接在GEE Code Editor中逐一查看样本点的真彩色影像// 以其中一个样本点为中心显示Sentinel-2真彩色影像 var pointSample ee.Feature(samples.first()); Map.centerObject(pointSample, 15); Map.addLayer(s2Composite, {bands: [B4, B3, B2], min: 0, max: 3000}, S2 RGB); // 查看该样本点的属性 print(Sample attributes:, pointSample);当你人工审核一批样本点后如果发现某些点被Hansen标为“2015年损失”但在Sentinel-2影像里一直是森林状态就说明Hansen在那个区域存在漏检或误检。你可以把所有这类误判样本的坐标记录下来反过来重新评估Hansen数据在特定区域的可信度。这正是“分层样本点”作为独立参考数据的关键价值——它不仅是被动的抽样工具更是主动的质量监测工具。不过人工审核高分辨率影像也有耗时陷阱。我在处理全球尺度样本时几千个样本点逐一点开查看往往要耗费数天时间。后来我用Google Earth Engine的Export.image.toDrive批量导出每个样本点的局部影像快照然后用FastStone之类的看图软件快速浏览效率提升明显。推荐你也把样本点分成小批次批量导出再集中判读既快又不容易漏。结尾想说的几句实在话分层抽样是整个GEE森林遥感分析中最容易被当作“配件”却又最影响结论质量的环节。我现在做任何全球尺度的森林面积或变化分析都会把stratifiedSample放在处理流程的核心位置而不是等分类都做完之后才想起来补样本。这背后的原因很简单样本代表不了总体后面的任何统计推断都是空中楼阁。建议你在跑正式实验前先用小范围的测试区域把分类波段、numPoints和scale这三个参数全部调通确认每层样本数和空间分布都符合预期再放大到全研究区跑。宁可开头多花一个小时调参也不要等到几万样本导出后才发现分层设计出了偏差那才是真正的时间杀手。 GEE分层抽样实战：基于Hansen数据提取森林损失样点
干我们这行的都知道GEE里坑最多的不是复杂的波段计算反而是样本提取这种“看着很简单”的环节。特别是做全球尺度的森林损失研究你不可能人工去逐点标注几十万个验证样点也没法用简单的随机抽样指望稀疏的损失像元能被抽中。去年我做全球森林变化产品精度评估时被stratifiedSample这个函数卡了整整三天反复研究后才发现问题往往不在函数本身而在你对“分层”这件事的理解。这篇文章就把我从数据准备、类别设计到分层样本提取、质量检查的全过程完整拆出来每一步都给可直接复用的代码和参数希望能让你少走几趟弯路。这个内容的核心是基于Hansen全球森林变化数据集用GEE内置的stratifiedSample分层抽样方法在“森林损失”这个类别分布极不均衡的目标上科学地提取分层样本点。它解决的痛点非常明确——如何在全球动辄几十亿像元的影像中让稀有类别比如某一年发生损失的像元也能被充分抽样同时保持样本的空间代表性和时间完整性。无论你是做森林变化监测、土地覆被产品精度验证还是想学GEE分层抽样的原理这套方案都值得完整看一遍。1. 这个需求背后到底在解决什么问题1.1 为什么偏偏选Hansen全球森林变化数据做森林损失量分析第一步就是选底图。GEE里不缺森林相关的数据集有MODIS的年度森林盖度、有ESA的全球土地覆被还有各种局部区域的林龄产品。但Hansen全球森林变化数据集UMD/hansen/global_forest_change几乎成了默认选项原因无非三点。第一它的空间分辨率是30米而MODIS那种250米到500米的产品对森林损失这种“碎片化”事件来说太粗了。一条道路建设或一次小规模皆伐在MODIS影像上可能只体现为一个混合像元根本没法数清损失了多少。Hansen的30米分辨率可以捕捉到大多数非法的、小规模的森林采伐。第二它提供了从2000年基线到逐年的损失年份lossyear信息。这个“逐年”属性特别关键——你不仅能知道哪里发生了损失还能知道是哪一年损失的这让“全球森林损失量”的分析从静态制图变成了时间序列分析。我在做数据产品精度评价时最看重的就是这一点。第三它在GEE中是以固定资产形式托管的访问速度快、不用自己做复杂的预处理。你要做的只是select出需要的波段剩下的几何校正、大气校正全都不必操心。但这里有个极其重要的实操细节很多人忽略了Hansen的lossyear波段里0代表“2000年至数据版本年份之间未检测到损失”1到数据版本年份才代表对应的损失年份。当你把它直接用于分层抽样时如果不做任何处理0类像元可能占据99%以上的面积你抽取的样点会被无损失像元淹没损失年份样本根本抽不到几条。这正是需要引入“分层”概念的直接原因。1.2 分层抽样 vs 简单随机抽样差在哪说句实话如果目标类别分布是均匀的stratifiedSample和普通的sample没什么太大区别。但在森林损失这个场景下二者差异是决定性的。简单随机抽样是在整个目标区域内完全随机地撒点类别分布和像元分布保持一致。问题在于如果我们研究的区域里只有0.5%的像元属于“2015年发生损失”那么随机抽10000个点理论上抽中2015年损失层的只有50个这50个点根本不足以支撑任何统计分析。分层抽样stratified sampling的思路则是先把研究区按类别规则切成若干个互不重叠的“层”strata然后在每一层内部独立地随机抽样。这样无论2015年损失层有多么稀疏只要你指定给它抽200个点它就一定会给你抽出200个。这就保证了每个类别都有“最低样本保障”同时每一层的样本又保持了随机性。用一个生活类比来说如果你想了解整个城市不同收入人群的消费习惯简单随机抽样就好比从电话簿上随机挑人结果高收入人群可能一个都没被抽到分层抽样则是先把市民按收入分成高、中、低三个名单再从每个名单里分别随机挑人。你所付出的额外成本是必须先知道“怎么分层”但换来的回报是每个层级都有足够的样本数量。在GEE里stratifiedSample函数做的就是这件事。它允许你指定一个分类波段作为分层的依据每一个独特的整数值视为一层然后在每一层内部执行随机采样。它内部会高效地做分桶统计和随机定位比你自己写reduceRegions然后拼接多个filter结果要可靠得多。2. stratifiedSample的核心机制与参数门道2.1 函数签名与每个参数的真实含义直接看GEE里stratifiedSample的调用方式API文档写得比较简洁实际含义得靠用才能体会。基本形式如下var samples image.stratifiedSample({ numPoints: 100, classBand: class, region: roi, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3] });每个参数的真实含义numPoints是每一层内计划抽取的样本点数。注意是“每一层”都抽这个数不是总共抽这个数。如果你的分类波段有5个类别值则总样本数大约是5乘以numPoints实际可能因某个类别像元不足而少于该值。classBand是分层依据的波段名称。它必须是整数类型Int或Byte因为分层就是按整数值归类。如果波段是浮点型你需要先重编码成整数否则函数会报错或者行为不符合预期。region定义了采样范围。可以是一个ee.Geometry、ee.FeatureCollection或ee.Feature。不传的话默认用影像全幅范围但全球30米分辨率逐类采样计算压力极大强烈建议传区域。scale是采样时考虑的空间分辨率。这个值的取舍很值得讲。用30米意味着逐像元级别的采样最精确但对计算资源的需求高用1000米意味着采样的“定位”粒度较粗适合全球快速预览。我的经验是如果研究区是一个国家或大洲级别且你的计算配额不算富余先用250米跑通全流程最后对关键区域再用30米做一次精细抽样。geometries是一个布尔值决定输出的FeatureCollection中每个样本是否携带点几何。如果设为false返回的记录只有属性值没有空间坐标设为true则有GeoJSON点。你做精度验证、人工目视判读或后续叠加高分辨率影像都必须设为true。seed是随机种子。这个参数的意义超过多数人的想象。分层抽样使用了随机数生成器如果不固定种子每次运行生成的样本位置都不同论文和报告的复现性就没法保证。在学术研究中审稿人要求“样本结果可复现”时你只需要亮出固定的seed即可。我建议用日期加版本号比如20240601清晰又不容易冲突。classValues是可选参数用于显式指定要采样的类别值列表。默认会对classBand中的全部类别分别抽样。如果你只关心特定年份的损失层这个参数可以避免多余的计算。2.2 分层变量该怎么构造重新编码成类别波段这是整个流程中最容易被低估的一步。很多人直接把Hansen的lossyear波段塞给classBand然后发现样本数量巨大、运算超时或者返回的结果里损失年份样本出奇地少。原因在于lossyear的取值并不适合直接作为分层依据。Hansen的lossyear取值是0到20的整数但在某些版本和特定区域2000年之前已经有森林覆盖但从未经历损失的大片区域全是0值。直接把0到20每个值当作一层会产生两个问题一是层数太多导致numPoints乘以层数的总量非常惊人二是0层占绝对优势你要重点关注的不同时期损失层被稀释在了过细的分层结构里。解决思路不是放弃lossyear而是在它基础上重新编码成一个“有意义的类别波段”。比如按时间窗口合并类别将森林损失划分为“未损失”、“2001-2005年损失”、“2006-2010年损失”、“2011-2015年损失”、“2016-2020年损失”五类。这样既保留了时间维度上的层次又不会造成层数过多。代码实现如下// 加载Hansen全球森林变化数据集 var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8); // 选择关键波段 var treeCover hansen.select(treecover2000); var lossYear hansen.select(lossyear); // 构建分层类别波段 // 0 未损失森林 (treecover2000 30 且 lossyear 0) // 1 2001-2005年间损失 // 2 2006-2010年间损失 // 3 2011-2015年间损失 // 4 2016-2020年间损失 var forestMask treeCover.gte(30); var classBand ee.Image(0) .where(forestMask.and(lossYear.gte(1).and(lossYear.lte(5))), 1) .where(forestMask.and(lossYear.gte(6).and(lossYear.lte(10))), 2) .where(forestMask.and(lossYear.gte(11).and(lossYear.lte(15))), 3) .where(forestMask.and(lossYear.gte(16)), 4) .clip(roi) .byte();你可能注意到了我把森林定义成treecover2000大于等于30%。这不是拍脑袋定的。Hansen团队官方建议森林/非森林的阈值常用30%这也是许多联合国报告和森林资源评估采用的标准。当然你完全可以根据研究需要改写成10%或50%但必须保持一致并在论文或报告里明确声明阈值。2.3 种子、尺度和几何输出的设计细节在GEE里使用stratifiedSample时还有几个容易踩的细节提前说清楚能省很多调试时间。关于scale的选择逻辑我强烈建议你做一个“两级采样”设计。第一级在全球或大区域尺度用scale参数设置为1000米甚至5000米跑一次分层抽样目的是验证分类波段是否正确、各类别样本空间分布是否合理、总体样本量是否符合预期。这一步成本极低一两分钟就能出结果。第二级在确认分类没问题后再对重点区域用scale: 30做精细抽样。这样可以把90%的调试时间花在低成本的第一级上而不会因为一个低级错误浪费掉几十分钟的计算配额。关于geometries参数还有一个小坑。当你设置geometries: true时每个样本点都带空间几何这个集合如果在GEE的Map面板上直接显示可能因为点数太多导致页面卡顿。建议先用limit或者按类别filter后再显示或者干脆只在导出Data后本地查看。另一个容易被忽略的点是stratifiedSample返回的FeatureCollection不仅包含classBand指定的类别值还会包含参与分层影像的所有波段信息。如果你的分层影像是一个包含多波段的Image比如你把treecover、lossyear都addBands进去了那么每个样本都会带上这些波段的像元值作为属性。这其实是一件好事——你在后续精度验证或者建模时不需要再次去影像里提取属性因为所有属性已经挂在样本点上了。但要注意如果影像中包含了不必要的波段比如原始Hansen数据中的数据质量波段输出的属性表会变得臃肿导出CSV后光字段理解就要花不少时间。所以采样前用select限定波段是个好习惯。3. 全球森林损失样本提取的完整实操3.1 准备数据与目标区域先明确我们要做什么。以东南亚地区为例我想提取覆盖五类森林状态未损失、四个时期损失的样本点每个类别各抽取300个用于后续的高分辨率影像目视判读和损失面积估算。东南亚是森林损失重灾区但也是行政边界复杂、影像噪声较多的地区用它来演示更有说服力。数据准备代码// 行政边界导入可以用FAO的GAUL数据集或全球行政区划数据集 var roi ee.FeatureCollection(FAO/GAUL/2015/level0) .filter(ee.Filter.inList(ADM0_NAME, [Malaysia, Indonesia, Thailand])) .geometry(); // 加载Hansen数据 var hansen ee.Image(UMD/hansen/global_forest_change_2020_v1_8); var treeCover hansen.select(treecover2000); var lossYear hansen.select(lossyear); Map.centerObject(roi, 5); Map.addLayer(treeCover.updateMask(treeCover.gt(0)), {palette: [#006400, #00ff00]}, Tree Cover 2000);森林盖度图层添加到地图上目的不是看热闹而是确认目标区域内森林的空间分布是否符合预期。比如东南亚的婆罗洲、苏门答腊这些热点区域应该有明显的森林覆盖如果发现城市中心也有大量森林盖度说明数据或区域投影有问题要及早排查。3.2 构造分层类别把连续损失变成离散层接着用之前提到的思路重编码分层类别并把需要携带的属性波段一起放在一个影像里// 构建森林区域掩膜 var forestMask treeCover.gte(30); // 构建分层类别波段 var classBand ee.Image(0) .where(forestMask.and(lossYear.eq(0)), 0) .where(forestMask.and(lossYear.gte(1).and(lossYear.lte(5))), 1) .where(forestMask.and(lossYear.gte(6).and(lossYear.lte(10))), 2) .where(forestMask.and(lossYear.gte(11).and(lossYear.lte(15))), 3) .where(forestMask.and(lossYear.gte(16)), 4) .clip(roi); // 创建用于采样的多波段影像 // classCode是分层依据lossYear和treeCover作为样本属性一并携带 var sampleImage ee.Image.cat([ classBand.rename(classCode), lossYear.rename(lossYear), treeCover.rename(treeCover) ]).byte(); // 给classCode波段起名为class因为stratifiedSample默认以该波段分层 var sampleImage sampleImage.select( [classCode, lossYear, treeCover], [class, lossYear, treeCover] );这里有两个细节。第一我用ee.Image.cat把多个波段拼接成一个多波段影像然后重命名。classCode重命名为class这样stratifiedSample里classBand: class一目了然。第二所有波段最后都转成了byte8位无符号整数。Hansen的lossyear和treecover取值都在0到100之间byte足够存储减少一半以上的内存和网络传输消耗。这在全球尺度的计算中是非常可观的优化。3.3 执行stratifiedSample并导出检查核心语句来了// 设置抽样参数 var numPointsPerClass 300; // 执行分层抽样 var samples sampleImage.stratifiedSample({ numPoints: numPointsPerClass, classBand: class, region: roi, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3, 4] }); print(Sample FeatureCollection, samples); print(Size of samples, samples.size());运行后在Console里观察如果一切正常输出集合大小应该在1500附近5类乘以300但实际情况可能因为某些类别在区域内像元不足而少于1500。比如2001-2005年东南亚的损失相对较少抽取300个点时可能只能抽到280-290个这是正常的说明该类别在给定区域内的有效像元数不够。为了快速检查各类别实际抽到了多少点可以用下面的代码分组统计// 按类别统计样本数量 var classCounts samples.reduceColumns({ reducer: ee.Reducer.frequencyHistogram(), selectors: [class] }); print(Class counts, classCounts);如果某个类别抽到0个样本而你的研究又必须包含该类别就需要放宽region范围或者下调numPoints之外的硬约束。比如可以用多期Hansen版本拼接数据。但多数情况下出现某个类别完全没有样本的情况几乎总是因为分层波段构建逻辑有误。3.4 样本可视化与快速质检拿到样本后第一件事不要急着导出先在地图上可视化检查空间分布。我习惯的配色方案是未损失层用深绿色四个损失时期分别用橙黄红紫四色这样一眼就能看出是否均匀分布在校准框内。// 将样本按类别分割开显示 var palette [#006400, #ff9900, #ff3300, #ff00ff, #0000ff]; Map.addLayer(samples, {color: 000000}, All samples (black));更实用的做法是把样本属性表导出为CSV用统计软件快速查看每类的空间坐标范围、lossyear分布直方图。具体导出代码如下// 导出样本点 Export.table.toDrive({ collection: samples, description: SEAsia_forest_loss_samples, fileFormat: CSV, folder: GEE_exports, selectors: [class, lossYear, treeCover, .geo] });导出后我习惯先在Excel里打开确认几个关键字段class的取值是否只在0到4之间lossYear的取值是否和class相互匹配比如class1时lossYear应该在1到5之间检查是否存在坐标为(0,0)或位于研究区之外的异常点。这一步虽然简单但往往能立刻暴露问题。我在一次实验中发现class42016-2020损失的点里混进了lossYear0的记录追根溯源发现是lossyear构建时分类波段和损失年份的掩膜顺序写反了重新调整处理逻辑后一切正常。4. 实战中一定会遇到的坑4.1 “每个类别的像元数不够”怎么解stratifiedSample最让人头疼的运行报错之一是每个类别要求抽取numPoints个点但某些层实际像元数少于numPoints。常见的提示是Number of pixels in strata number of requested samples。这事我遇到过太多次尤其是当研究区小、且某个损失年份特别稀疏时。解决办法有几种按优先级排列第一种调低numPoints让每个类别的期望样本数小于最少类别的有效像元数。你可以先用reduceRegion或classBand.reduceFrequencyHistogram快速统计各层像元数量再据此设定numPoints。第二种扩大region让稀有类别获得更大的空间范围来补充有效像元。比如把研究区从单一州郡扩大到整个省。第三种调整scale。scale从30米改成100米或250米意味着每个“采样像元”覆盖更大的地面范围类别有效像元数也随之增多。代价是采样点的定位精度下降所以它更适合全球尺度的粗略研究不适合需要精细点位的研究。第四种如果实在无法满足就拆分区域。把研究区分成几个子区域分别执行分层抽样最后用merge合并所有样本。这个方法多花一点代码和时间但能保留30米的分辨率也是我处理非洲多国联合分析时最常用的一招。// 拆分区域后分别抽样再合并 var subregions [regionA, regionB, regionC]; // 每个子区域都要有对应几何 var allSamples ee.FeatureCollection([]); subregions.forEach(function(name) { var subRegion ee.FeatureCollection(rois).filter(ee.Filter.eq(name, name)).geometry(); var subSamples sampleImage.stratifiedSample({ numPoints: numPointsPerClass, classBand: class, region: subRegion, scale: 30, geometries: true, seed: 20240601, classValues: [0, 1, 2, 3, 4] }); allSamples allSamples.merge(subSamples); });需要注意的是GEE的forEach不是同步执行的合并时最好用ee.FeatureCollection的flatten方式或者在ui模式下改用map。上面的代码只是示意实际使用推荐写成函数式合并。4.2 采样点太密内存溢出与地图卡顿一次性生成几千个带几何的样本点在GEE的Map面板上直接显示十有八九会卡顿甚至白屏。这不是程序出bug而是客户端渲染压力过大。我的处理方式是把样本点写到自己的资产里或者导出GeoJSON到本地ArcGIS/QGIS查看。如果必须在GEE里预览部分样本可以这样// 显示每类前50个点共250个点足够预览 var previewSamples samples.limit(250); Map.addLayer(previewSamples, {color: red}, Preview samples (250));另一个内存陷阱是想把样本点和高分辨率遥感影像叠加做目视判读时直接addLayer高分辨率影像同样会让浏览器崩溃。解决方案是使用getThumbURL生成缩略图或者只给单个样本点周围创建一个小缓冲区域再显示。实际操作中我更倾向于用GEE Code Editor的Export.map.toDrive把注意力转到导出任务不在线上做深度交互。4.3 不同版本的Hansen数据波段差异GEE上的Hansen数据集更新过多个版本有v1_7、v1_8、v1_10等。版本之间虽然核心波段设计一致但也有细微差异。比如较新版本增加了datamask、gain等波段损失年份的最大值也不同v1_8损失年份到2020年v1_10则可以到2022年甚至更晚。这意味着如果你在前面代码里写死了lossYear.gte(16)在v1_10中会覆盖所有2020年之后的损失得到的结果维度完全变了。我的建议是研究开始前明确使用哪个版本并检查对应版本的最大损失年份。代码里用动态方式获取版本上限比写死数值更健壮var lossYear hansen.select(lossyear); var maxLossYear lossYear.reduceRegion({ reducer: ee.Reducer.max(), geometry: roi, scale: 1000, maxPixels: 1e9 }).getNumber(lossyear); print(Max loss year in ROI:, maxLossYear);这条命令可以快速确认你所选择的数据版本在目标区域内最晚的损失年份。许多同学跑完整套流程后才发现版本内损失年份的上限和预期不符不得不重算一遍费时费力。在第一步就确认版本和最大年份是一个好习惯。4.4 常见问题速查表为了更直观我把使用stratifiedSample时常见的坑整理成一张速查表建议收藏常见问题可能原因解决方案报错No valid pixels in strataclassBand在ROI内全为掩膜值检查classBand的掩膜范围确保每个类别在ROI内有有效像元某个类别样本数为0类别编码冲突或区域太小检查recode条件扩大ROI或调低scale两个类别的样本点重合分层波段取值定义重叠检查.where条件之间是否有边界重叠逐层用mask验证导出CSV字段过多未限定采样波段采样前做select只保留必要波段每次运行结果不同未设置seed固定seed确保结果可复现Map显示卡死样本点数过多用limit预览或导出本地查看损失年份最大年份异常使用了不同版本Hansen动态读取maxLossYear并检查数据集ID5. 样本拿到之后还能干什么5.1 用样本做面积估算与精度验证分层样本的价值绝不只是“抽出来了就完事”。在森林损失研究中最常见的后续工作是面积估算和分类精度验证。以面积估算为例Olofsson等人提出的良好实践指南里推荐使用分层抽样stratified random sampling得到的样本结合每层的像元数权重来修正面积估计。如果你把整个区域按“未损失”和“损失”两层抽样那么损失面积的无偏估计可以写成损失面积估计 各类别损失像元比例乘以各类别总面积再对各层加权合并。具体公式为A_hat Σ (W_h × A_total × p_h)其中W_h是第h层像元占比p_h是第h层中样本里判定为“损失”的比例A_total是整个区域总面积。通过这样的计算你得到的面积估计值比单纯的“损失像元乘以30×30平方米”要更稳健尤其适合样本量不大、类别不平衡的数据。在GEE中你可以先用reduceRegion拿到每个类别的像元数量和面积权重再用样本点和高分辨率影像人工判读的结果去更新p_h。这两个东西组合起来就能算出带置信区间的损失面积。具体代码实现较多不在这里展开但你要知道stratifiedSample抽出的样本是衔接“遥感分类结果”和“地面真实面积”的桥梁。另外如果你后续要训练一个自定义分类器比如区分火灾、砍伐、风暴导致的损失stratifiedSample产出的样本点还可以作为训练样本。做法很简单把样本点叠加到高分辨率影像切片上人工标注出损失原因再把标注后的点集灌回分类器。这时候你在分层抽样时保留的lossYear和treeCover属性就派上了大用场——它们可以作为辅助特征一起参与模型训练。5.2 结合高分辨率影像做人工审核的流程分层样本落到具体使用中尤其是做精度评估时光靠Hansen的30米分辨率是不够的你需要用更高分辨率的影像Planet、Maxar、Sentinel-2的10米来确认每个点真实的地表状态。实际操作中我的流程是这样的先把stratifiedSample导出的样本在GEE里和Sentinel-2无云影像合成叠加然后导出为PNG切片再放到本地影像管理软件中人工判读。更高效的做法是直接在GEE Code Editor中逐一查看样本点的真彩色影像// 以其中一个样本点为中心显示Sentinel-2真彩色影像 var pointSample ee.Feature(samples.first()); Map.centerObject(pointSample, 15); Map.addLayer(s2Composite, {bands: [B4, B3, B2], min: 0, max: 3000}, S2 RGB); // 查看该样本点的属性 print(Sample attributes:, pointSample);当你人工审核一批样本点后如果发现某些点被Hansen标为“2015年损失”但在Sentinel-2影像里一直是森林状态就说明Hansen在那个区域存在漏检或误检。你可以把所有这类误判样本的坐标记录下来反过来重新评估Hansen数据在特定区域的可信度。这正是“分层样本点”作为独立参考数据的关键价值——它不仅是被动的抽样工具更是主动的质量监测工具。不过人工审核高分辨率影像也有耗时陷阱。我在处理全球尺度样本时几千个样本点逐一点开查看往往要耗费数天时间。后来我用Google Earth Engine的Export.image.toDrive批量导出每个样本点的局部影像快照然后用FastStone之类的看图软件快速浏览效率提升明显。推荐你也把样本点分成小批次批量导出再集中判读既快又不容易漏。结尾想说的几句实在话分层抽样是整个GEE森林遥感分析中最容易被当作“配件”却又最影响结论质量的环节。我现在做任何全球尺度的森林面积或变化分析都会把stratifiedSample放在处理流程的核心位置而不是等分类都做完之后才想起来补样本。这背后的原因很简单样本代表不了总体后面的任何统计推断都是空中楼阁。建议你在跑正式实验前先用小范围的测试区域把分类波段、numPoints和scale这三个参数全部调通确认每层样本数和空间分布都符合预期再放大到全研究区跑。宁可开头多花一个小时调参也不要等到几万样本导出后才发现分层设计出了偏差那才是真正的时间杀手。