一开始接触 GEE 里的 Landsat 8 数据时很多人会把Landsat 8 SR和T1_L2当成两个差不多的数据集。实际上这俩是两代产品一个停在 Collection 1一个是现在主推的 Collection 2 Level-2。刚好最近在辽宁省做土地覆被和时序分析几乎每一条影像都要做去云。如果只靠网上复制一段去云函数很容易翻车——因为两代数据的 QA 波段、波段名、数值缩放比例都不一样。这篇文章就结合辽宁省的实例把LANDSAT/LC08/C01/T1_SR简称 L8 SR和LANDSAT/LC08/C02/T1_L2简称 T1_L2的去云区别掰开揉碎讲一遍包括代码、坑位和验证方法希望你能少走点弯路。1. 你手上的 Landsat 8 到底属于哪一代SR 与 T1_L2 的身份辨析1.1 Collection 1 与 Collection 2 的来龙去脉先说说产品代际。USGS 在 2017 年发布 Landsat Collection 1而 GEE 里常见的LANDSAT/LC08/C01/T1_SR就是这一代中的“Tier 1 表面反射率产品”。它的关键点在于数据是表面反射率SR已经做过大气校正不是原始的 TOA 辐射亮度。很多老教程里的去云代码默认数据源就是它。大约从 2021 年开始USGS 逐步把 Collection 2 作为正式标准GEE 里对应的就是LANDSAT/LC08/C02/T1_L2。这里的T1是 Tier 1表示几何精度高、适合时序分析L2是 Level-2代表经过大气校正后的地表反射率产品。Collection 1 的影像现在已经停止更新GEE 官方也建议把新流程迁移到 Collection 2。所以从本质上讲你不需要把“Landsat 8 SR”和“T1_L2”理解成两种截然不同的卫星数据它们其实是同一个传感器的两代处理产品。一个是旧家中老底子一个是新家的正式版。1.2 数据集 ID 和波段命名的三大差异在写代码之前最重要的就是认清波段结构。下图只用列表和表格的方式把差异列出来代码里不会踩到波段名不存在的坑。对比项老版 L8 SRCollection 1新版 T1_L2Collection 2GEE 数据集 IDLANDSAT/LC08/C01/T1_SRLANDSAT/LC08/C02/T1_L2红绿蓝近红外波段B2、B3、B4、B5SR_B2、SR_B3、SR_B4、SR_B5QA 波段名pixel_qaQA_PIXEL辐射饱和 QAradsat_qaQA_RADSAT气溶胶 QAsr_aerosolSR_QA_AEROSOL热红外波段B10、B11ST_B10地表温度等这里最容易出问题的是波段命名。在 C1 里你用select(B4)取红波段虽然它其实是表面反射率但波段名和 TOA 数据一样。到了 C2USGS 明显学聪明了直接改成SR_B4从名字上就告诉你这是表面反射率。如果你把 C1 的去云代码直接套到 C2 上select(pixel_qa)这一句就会直接报错偏偏这个报错还不会提示“没有这个波段”只会在后面输出结果时突然全空白很迷惑人。另外还有尺度因子。C1 的 SR 反射率值存放在 16 位整型里转换关系是DN * 0.0001。C2 改成了更精细的DN * 0.0000275 (-0.2)。这意味着两端数据明面上都叫“反射率”实际数值范围完全不一样混用的时候必须手动统一。1.3 什么时候该用 C1什么时候该用 T1_L2我的建议很简单新项目一律用LANDSAT/LC08/C02/T1_L2。C1 已经不再更新2018 年之后的影像你想拿新的也拿不到这等于直接断了时序分析的后路。除非你是在复现一篇很老的论文或者项目里已经沉淀了大量基于 C1 的流程代码否则没必要继续守旧。但这里有个例外如果你做长时序分析需要把 2013 年以来的数据都拼起来那么 T1_L2 本身已经覆盖了完整时间段没必要混合 C1 和 C2。反倒是历史对比分析里C1 的某些像元标识和 C2 有细微差别最好别混着用做逐像元统计。总之记住一个原则去云逻辑要跟着产品代际走不能一套代码吃遍全部数据。2. 去云的核心机制QA 波段位编码的五四三2.1 QA_PIXEL 和 pixel_qa 的位编码对照去云的底层逻辑不是用机器学习识别云而是读取影像自带的 QA 波段。这个 QA 波段里面每个像元是一个 16 位的整数每一位bit代表一种属性。比如第 3 位bit 3如果是 1说明这个像元被判定为云第 4 位是 1说明是云影第 5 位是 1说明是雪或冰。C1 的pixel_qa和 C2 的QA_PIXEL在比特位的定义上有很大的重叠但并不完全一样。C2 在 C1 的基础上增加了卷云置信度、雪/冰置信度等额外信息同时把部分旧标记做了细化。下面是常用位含义的对照表。bit 位C1 pixel_qaC2 QA_PIXELBit 0FillFillBit 1Dilated CloudDilated CloudBit 2CirrusCirrusBit 3CloudCloudBit 4Cloud ShadowCloud ShadowBit 5SnowSnowBit 6ClearClearBit 7WaterWaterBit 8-9Cloud ConfidenceCloud ConfidenceBit 10-11Cloud Shadow ConfidenceCloud Shadow ConfidenceBit 12-13无Snow/Ice ConfidenceBit 14-15无Cirrus Confidence理解位编码之后去云函数无非就是做一件事把某个 bit 为 1 的像元筛掉。GEE 的写法通常是qa.bitwiseAnd(1 3).eq(0)意思是“第 3 位是否为 0”。如果为 0说明不是云保留否则掩膜掉。2.2 为什么 C2 的 QA_PIXEL 比 C1 更适合做云掩膜很多人以为两代数据的 QA 波段只是换了个名字其实算法也在更新。C1 的 QA 数据用的是早期 CFMask 算法对薄云和卷云的识别不算稳定。C2 的 QA_PIXEL 是 USGS 重新处理 Collection 2 时生成的整体上使用更新版本的 CFMask对卷云、云影、尤其是山体阴影附近的误判有所改善。在实际项目里最直观的感受是辽宁省沿海地区春季经常有低云和海雾C1 的 QA 波段会漏掉一部分低云导致 NDVI 时序里突然出现一个很低很假的点C2 因为多检查了 Dilated Cloud 和 Cirrus 位漏检概率更低。当然 C2 也不是神它也会把河谷里的小块雾误判成云只是总体更稳。2.3 SR 数值的范围坑同样叫反射率量纲不一致位编码之外还有一个和去云没有直接关系但经常被忽略的坑SR 波段数值的范围。C1 的反射率波段需要乘0.0001C2 需要乘0.0000275再减0.2。如果你写死一个可视化范围比如min: 0, max: 0.4在 C1 上没问题但在 C2 上DN 值 8000 乘以0.0000275再减0.2出来的反射率是0.02左右中间的暗色像元可能直接黑掉。更严重的是做多期影像数值统计时如果不统一缩放NDVI 会出现系统性偏移看起来像是“地表发生了剧烈变化”其实是数据代际差异。我建议所有基于 C2 的统计、分类、深度学习样本制作都先把整幅影像转成浮点反射率并且把system:time_start属性保留好避免后面的时序处理把时间信息丢掉。3. 以辽宁省为例同一块地方的两种去云操作对比3.1 辽宁区域与数据准备辽宁省地处东北南部地形比较复杂中部是辽河平原东部是山地西部是低山丘陵南部是辽东半岛和沿海。春季和初夏经常有云沿海还容易起雾所以 Landsat 8 16 天重访周期碰上可用晴天影像的概率并不高。要做时序分析去云不是可选项而是必要的预处理。下面这一段代码先框选了辽宁省周边范围并筛选出 2021 年 4 到 9 月之间的影像。为了对比我把 C1 和 C2 两个数据集都加载进来。var liaoNing ee.Geometry.Rectangle([119.5, 38.7, 125.8, 43.5]); var l8sr ee.ImageCollection(LANDSAT/LC08/C01/T1_SR) .filterBounds(liaoNing) .filterDate(2021-04-01, 2021-09-30) .filterMetadata(CLOUD_COVER, less_than, 40); var l8c2 ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterBounds(liaoNing) .filterDate(2021-04-01, 2021-09-30) .filterMetadata(CLOUD_COVER, less_than, 40); print(l8sr.size(), l8c2.size());注意CLOUD_COVER是元数据里的大范围云量估计不能替代逐像元去云。筛选它的主要目的是减少无效计算尤其可以减少大批量处理时的时间消耗。3.2 两套去云函数的逐行解释C1 版本的去云函数如下function maskL8sr(image) { var qa image.select(pixel_qa); var cloud qa.bitwiseAnd(1 3).eq(0); var cloudShadow qa.bitwiseAnd(1 4).eq(0); var cirrus qa.bitwiseAnd(1 2).eq(0); var dilated qa.bitwiseAnd(1 1).eq(0); return image.updateMask(cloud.and(cloudShadow).and(cirrus).and(dilated)); }这段代码很简单逐位检查云、云影、卷云、膨胀云。但 C1 的 QA 波段里没有额外的置信度位来辅助判断所以它对“可能不是云但置信度中等”的像元处理比较粗放。C2 版本的去云函数如下function maskL8c2(image) { var qa image.select(QA_PIXEL); var cloud qa.bitwiseAnd(1 3).eq(0); var cloudShadow qa.bitwiseAnd(1 4).eq(0); var cirrus qa.bitwiseAnd(1 2).eq(0); var dilated qa.bitwiseAnd(1 1).eq(0); var snow qa.bitwiseAnd(1 5).eq(0); // 更严格的一个版本加上云置信度不超过 low 的条件 var cloudCon qa.bitwiseAnd(3 8).rightShift(8); var conMask cloudCon.lt(2); return image.updateMask(cloud.and(cloudShadow).and(cirrus).and(dilated).and(snow)); // return image.updateMask(cloud.and(cloudShadow).and(cirrus).and(dilated).and(snow).and(conMask)); }看起来好像只是波段名从pixel_qa变成了QA_PIXEL其实还可以利用 C2 的置信度位做更严格去云。注释掉的那一行就是开启“云置信度小于 low”的过滤如果你要非常严格地清掉所有可疑像元可以打开它。但要注意辽宁有些林区和草地像元会被误标为低置信度的云开太狠会减少有效像元数量反而影响覆盖度。3.3 实测效果辽东湾沿海一景的掩膜结果我在辽宁省南部选了一景 2021 年 6 月的影像分别用 C1 和 C2 的 QA 波段做去云。两套数据的 QA 对同一块区域的判定并没有完全重合。用 GEE 统计纯云像元的比例C2 检出的云和云影面积通常比 C1 更大一些尤其是薄卷云和低云边缘。实际表现大概是这样在辽东湾靠近海面的区域C1 的pixel_qa会把部分海雾当作水体保留下来导致影像上出现一大片灰白色噪声C2 的QA_PIXEL则会把它们判定为 Dilated Cloud 或 Cirrus直接掩膜掉。对于内陆农田区域的云影C2 的识别也更贴近云的投影方向剩下的云影残边比 C1 少。当然这并不意味着“C2 去云后绝对干净”。如果你把掩膜后的影像放大到像元级依然能看到一些碎小的残云或云影点尤其是在茂密植被和低反差地块上。所以后面我还会专门聊验证和修边的问题。4. 去云代码实战8 个常见坑与绕行方案4.1 把 QA_RADSAT 当成 QA_PIXEL这不是段子我见过很多次。C2 里有两个容易搞混的波段QA_RADSAT和QA_PIXEL。QA_RADSAT表示辐射饱和状态与云完全无关。如果你误选它去做 bit 运算结果是整个掩膜乱七八糟或者根本没有掩掉任何云。绕行方案很简单写代码时先print(image.bandNames())看一眼C2 的QA_PIXEL是全部大写C1 的pixel_qa是全小写。复制粘贴别人的代码时尤其要注意这一点。4.2 位运算的优先级问题在 JavaScript 的 GEE API 里最容易犯错的是位运算的拼接。正确的写法是qa.bitwiseAnd(1 3).eq(0)但有时候你会看到别人写成qa.bitwiseAnd(1 3 0)这种是错的因为优先级高于bitwiseAnd结果完全不对。我的习惯是把位检查单独写成一个个变量再最后合到一起阅读和排查都方便。另外如果你要取某个连续 bit 段比如bits 8-9的云置信度需要用(3 8)和右移rightShift(8)配合。这个在 C2 里很常用别漏了右移否则你得到的还是一个 0-255 范围的数没法直接和 0-3 的置信度等级比较。4.3 只去云不去云影或者把所有雪地都清掉在辽宁省做时序分析春季融雪期和 3、4 月的残雪很常见。如果你只去云不去雪积雪像元会被当成反射率突变的点NDVI 时序里会突然出现一个低值反过来如果直接把所有 snow bit 为 1 的像元都掩掉地表积雪信息就完全丢失做土地利用分类时也会少掉一整个类别。我的建议是根据项目目标来决定如果是做植被 NDVI 时序那么把雪去掉是对的如果是做土地覆被分类可以不急着去雪把它单独作为一个类别或者在下游分类模型里区分。C2 的好处是它有 Snow/Ice Confidence你甚至可以只掩掉高置信度的积雪保留低置信度的融雪像元。4.4 C1 和 C2 混合使用时的缩放问题这个坑我踩过最久。假设你辛辛苦苦把 C1 的数据全部去云然后想跟 C2 的数据合并成一个影像集合做时序如果没做缩放统一C1 的反射率乘的是0.0001C2 乘的是0.0000275再减0.2两组数据的数值实际上不对齐。表面看起来都叫 SR但统计结果差一大截。正确的做法是先对两个集合分别做map把波段数值统一成浮点反射率再合并。同时注意波段名也要统一最好都用SR_B4这样的命名避免B4和SR_B4混在一起。这个工作虽然繁琐但在项目启动时做好后面能节省大量调试时间。4.5 时间范围的“隐形坑”C2 的 Landsat 8 数据从 2013 年开始有但这不意味着每年每个月都有有效影像。辽宁省春季经常连雨夏季对流云多用filterDate(2020-04-01, 2020-06-30)筛选时可能只剩两三景。有人会把范围扩到前一年秋季或后一年春季这也不是不行但要明确跨季节的影像在地表状态上差异很大合并做分析时必须考虑物候漂移。另外如果做的是 2013 年之前的历史分析Landsat 8 的 C2 产品根本不存在只能去用 Landsat 5/7 的 C2 数据。别以为所有数据源都能无缝替换去云逻辑要跟着传感器的波段设计走。4.6 ImageCollection 的 map 里不要 return 外面定义的对象在 GEE 的map回调用我有时会看到有人在函数外定义一个mask someValue然后想在函数内修改它。这是典型的变量作用域混淆。map里面的函数应该只返回处理后的 Image并且所有判断都在 GEE 的服务端对象上进行。不要试图用print在循环里逐景输出那会把浏览器内存撑爆也会让代码非常慢。如果你的影像集合非常大建议先用CLOUD_COVER粗筛再在map里做细去云最后用filter再筛一遍处理后剩余像元比例。这样能显著提升运算效率。4.7 去云后立即统计导致的“空影像”去云后有些影像可能整景都被掩膜掉变成全空影像。如果不做过滤这类影像参与后续合成或统计会影响结果出现大片 NoData。我通常在去云之后再加一个filter用reduceRegion统计有效像元数量把有效像元太少的影像直接丢掉。var validCount image.select(SR_B4).reduceRegion({ reducer: ee.Reducer.count(), scale: 1000, maxPixels: 1e9 });这里的scale可以适当调大加快速度只要能判断大致是否有有效区域即可。4.8 布尔掩膜没有.and()连接在 GEE 里cloud.and(cloudShadow)返回一个新的布尔影像这一步是必要的。有人喜欢用*乘结果虽然 0、1 相乘逻辑上等价但到了浮点影像上容易把掩膜值变成非 0/1导致后面可视化时出现各种半透明像元。去云掩膜逻辑越清晰越好别偷懒。5. 去云结果的验证方法教你判断掩膜到底准不准5.1 先把掩膜可视化出来看很多人去完云直接把数据加到图层里也不看掩膜边界最后合成结果出现空洞才回头找原因。我的习惯是先输出两张图原始影像的真彩色叠加一个掩膜边界掩膜后的影像用半透明背景展示保留像元。在 GEE 里可以用下面这样把 QA 波段渲染出来。var qaVis {min: 0, max: 65535, palette: [black, white]}; Map.addLayer(l8c2.first().select(QA_PIXEL), qaVis, QA_PIXEL raw);然后逐级放大到辽宁沿海、山地、城市等不同地表看看掩膜边缘是否贴合云的轮廓。城市里的高亮屋顶和裸土很容易被 CFMask 误判成云这一步能快速发现。5.2 用 NDVI 时间序列检查漏云掩膜再怎么精细都可能有漏网之鱼。一个很实用的方法是对同一位置的时间序列点做 NDVI 曲线。正常植被 NDVI 曲线在生长季应该是平滑上升、稳定、再下降的如果某个节点突然断崖式下掉然后下一期又恢复正常多半是残云或云影没有掩干净。可以在辽宁省任意选一个农田或林地多边形用 GEE 的ui.Chart.image.series输出 NDVI 曲线。看到异常的负值就回头检查那一期的掩膜。这个方法不需要额外数据操作简单而且能定位到具体是哪一景出了问题。5.3 借助 Sentinel-2 辅助抽检如果项目正好同时使用 Sentinel-2可以在相近时间窗口内对比 Sentinel-2 官方场景图或经过云掩膜处理的影像来验证 Landsat 8 去云结果。不是所有项目都有同期好影像但哪怕一个月内的参考也能帮你判断这块区域到底是“低反射率地表”还是“阴影/云影”。在辽宁省山地、水体交错区域这种抽检尤其有效。需要特别提醒的是不要把所有阴影都当作云影。山地地形阴影和云影在光谱上可能很接近C2 的 QA 波段会把一部分地形阴影标为云影如果你全信 QA相当于白白丢掉了大量阴坡像元。6. 去云之后的事时序合成与深度学习样本的衔接6.1 用去云影像做年度合成和中值合成去云的最大价值在于让后续的合成变得可信。做完去云后最直接的应用就是做中值合成或最大值合成。下面这段代码以 C2 的 T1_L2 为例生成辽宁省 2021 年 5-8 月的 NDVI 中值合成。var c2Masked l8c2 .map(maskL8c2) .map(function(img) { var sr img.select(SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7) .multiply(0.0000275).add(-0.2) .rename(blue, green, red, nir, swir1, swir2); var ndvi sr.normalizedDifference([nir, red]); return sr.addBands(ndvi.rename(ndvi)) .copyProperties(img, [system:time_start]); }); var composite c2Masked.select(ndvi).median().clip(liaoNing); Map.addLayer(composite, {min: 0, max: 10000, palette: [gray, green]}, NDVI median);注意我这里把ndvi放大了10000倍来可视化原因是 NDVI 原始值很小在图层渲染里不放大很难看出层次。如果你直接输出统计要记得还原。6.2 为深度学习模型准备干净样本的关键点近几年“GEE 深度学习”这个词越来越热很多人会把去云后的影像导出为训练样本喂给语义分割模型做土地覆盖分类或云检测。这里有一个容易忽视的问题如果你用 QA 波段生成的掩膜作为训练标签等于默认了 QA 波段完全正确。但 QA 波段本身有漏检和误检特别是在云边界和阴影区域。我从实际操作中得到的经验是在制作深度学习训练集之前先对掩膜做一次形态学收缩比如用focal_min或邻域统计把云边缘向内收缩 1-2 个像元。这样能明显减少“云边缘半透明像元”被当成干净地表的情况。更稳妥的做法是只挑选那些被多个时相互相验证、在所有 QA 中都标记为 clear 的像元作为训练样本把不确定区域留给模型自己去学习边界特征。6.3 从去云到更高阶应用物候、温度与覆盖分类如果你已经能把 C2 T1_L2 影像稳定去云剩下的路就好走多了。辽宁省中部有大量农田春季出苗期如果遇到连续云雨天气不用去云函数你会发现自己根本凑不出连续 4 个时期的影像。去云之后即使不是每一景都完全干净至少能通过多景拼接或插值得到作物物候曲线。还有一点想提醒你地表温度产品ST_B10在 C2 里也是 Level-2 的一部分它与SR_B4等波段一样经常被云层遮挡。去云时不要只掩掉光学波段也要把热红外波段一起 mask 掉否则合成出来的地表温度栅格会出现拼接温度异常值。这个细节我一开始没注意后来处理辽宁省沿海海面温度时才意识到非常影响后续分析。我个人在实际操作中的感受是去云这件事本身不难难的是你永远要对自己的掩膜结果保持怀疑。在辽宁省这种云量随机性很大的区域最稳定的做法就是把 C2 数据、CFMask 的 QA_PIXEL、人工抽检这三者结合形成一个不可分割的流程。如果你现在正准备用 Landsat 8 做时序分析建议直接选 T1_L2它既能覆盖历史时间段也有更好的 QA 质量值得你花时间把去云函数写好。