在 GEE 里处理 Landsat 8 影像去云几乎是一道绕不过去的坎。尤其当你既看到过“Landsat 8SR”的代码又看到过“T1_L2”的代码很容易被两套名字搞混表面看起来都是“地表反射率产品”但实际用起来波段名不一样、QA 波段不一样、缩放公式也不一样。这篇文章就围绕辽宁省的实测场景把新旧两套数据集的去云区别一条条拆开顺带给出经过验证的 NDWI 水体提取脚本。无论你是第一次在 GEE 里做去云还是正在把旧代码迁移到新数据应该都能直接抄到一份可用的作业。1. 两个数据集到底是什么别被缩写绕晕1.1 从 Collection 1 到 Collection 2命名规则变化很多新手第一次在 GEE 搜索 Landsat 8 数据时会看到一堆非常相似的名字LANDSAT/LC08/C01/T1_SR、LANDSAT/LC08/C02/T1_L2、LANDSAT/LC08/C02/T2_L2再加上C01/T1_RT之类的老古董。这其实不是 GEE 故意为难人而是 USGS 的 Landsat 产品经历了一次代际升级从 Collection 1 切换到了 Collection 2。题主说的“Landsat 8SR”我理解就是指旧版LANDSAT/LC08/C01/T1_SR也就是 Collection 1、Tier 1、Surface Reflectance 产品。而“T1_L2”指的是新版LANDSAT/LC08/C02/T1_L2也就是 Collection 2、Tier 1、Level-2 产品。注意这里 L2 是 Level-2 的缩写表示这一级产品已经完成了大气校正提供的是地表反射率数据不是原始 DN 值。这里有个特别容易踩的坑Collection 1 老产品直接叫T1_SRSR 就是 Surface ReflectanceCollection 2 新产品叫T1_L2L2 其实也是 Surface Reflectance只是官方编号变了。所以很多人看到带 SR 的老代码很亲切看到 T1_L2 反而不认识其实两者都是“表面反射率”但一个是上一代一个是新一代。不搞清楚这个背景后面去云、算 NDWI 时代码很容易张冠李戴。GEE 同时支持这些数据集读取但官方早就停止生成 Collection 1 新产品了存量历史影像还能查新项目不建议继续当成默认数据源。T1 和 T2 的说法也顺便说清楚。Tier 1 指的是几何精度高、辐射质量稳定、适合做时间序列分析的数据Tier 2 是几何控制点不足、精度稍差的数据。标题里的 T1_L2 就是 Tier 1 的 Level-2 产品这也是做时序分析最推荐的数据集。T2 的 QA 波段和 L2 波段结构一样去云代码可以复用但涉及跨日期比较时尽量用 T1。1.2 新旧波段差异速查表两代产品最气人的地方就在这明明都是地表反射率波段名却完全不是一回事。旧版用B1、B2、B3这种简洁命名新版非要在前面加一个SR_前缀变成SR_B1、SR_B2、SR_B3。QA 波段也从pixel_qa改成了QA_PIXEL。我第一次迁移代码时直接用旧代码跑新数据结果 console 里一连串Image.select(...): Band B1 does not exist当场就明白了什么叫版本代沟。为了看起来方便我整理了一张速查表数据项旧版 C1 T1_SR新版 C2 T1_L2GEE 集合 IDLANDSAT/LC08/C01/T1_SRLANDSAT/LC08/C02/T1_L2多光谱波段命名B1-B7SR_B1-SR_B7质量波段pixel_qaQA_PIXEL辐射饱和波段radsat_qaQA_RADSAT反射率缩放DN x 0.0001DN x 0.0000275 - 0.2产品状态已停止更新当前正式版典型真彩色B4、B3、B2SR_B4、SR_B3、SR_B2典型 NDWI(B3-B5)/(B3B5)(SR_B3-SR_B5)/(SR_B3SR_B5)这里必须重点说一下缩放公式。很多人做 NDWI 时直接拿原始 DN 做归一化比值旧版 C1 SR 的反射率缩放系数是 0.0001缩放后数值在 0 到 1 之间DN 直接做比值误差不算大。但新版 C2 L2 的公式是DN * 0.0000275 - 0.2带了一个 -0.2 的偏移量。如果你忽略了后面这个 -0.2算出来的反射率会整体偏亮 0.2NDWI 的分母分子都会受影响暗像元尤其明显。所以正式算指数前一定要先做缩放。1.3 为什么推荐新项目用 C2 T1_L2我个人的建议是新项目无条件选LANDSAT/LC08/C02/T1_L2理由有三条。第一USGS 已经全面切换到 Collection 2C1 数据源停止更新新场景影像不会再进入 C1 集合。如果你做的是持续性的监测项目用 C1 很难跟上数据更新节奏。第二C2 的几何定位精度有提升在辽宁省这种地形起伏不太大但沿海区域较多的地区不同日期影像之间的配准更稳去云后的中值合成不容易出现重影。第三C2 L2 的 QA 波段设计更贴近后续产品比如 Landsat 9 的 L2 数据也是同一套结构代码一次写好L8、L9 都能复用。当然C1 也不是一无是处毕竟大量老教程、老论文、开源代码都基于 C1很多亚区级的长时序实验是在 C1 上做的。如果你只是复现别人十年前的分析C1 历史数据本身没毛病。但我建议做任何“新起点”项目时直接把 C2 T1_L2 作为默认数据源遇到老代码再逐步迁移成本其实很低主要就是改波段名、改缩放、改 QA 波段名这三件事。2. QA 波段与去云原理bitmask 怎么看怎么用2.1 像素级的“合格证”一组二进制的开关去云之前先理解 QA 波段到底存了什么。Landsat 的 QA 波段不是一张普普通通的“云量百分比图”而是把每个像素的多个质量信息压缩进一个 16 位整数里。你可以把它想象成每个像素随身带了一块面板面板上有 16 个指示灯每一位开关表示一种状态。bit 为 1 代表“是/有问题”bit 为 0 代表“否/正常”。问题在于GEE 里我们看到的 QA 值是一个十进制数比如 2500它本身没法直观看出哪几位亮了。所以要去云第一步就是把十进制数拆成二进制然后检查我们关心的那几个位。GEE 里拆位用的是bitwiseAnd和位运算不是空泛的原理而是每行去云代码都要用到的实际工具。比如我们关心 bit 3 是不是 1就可以用qa.bitwiseAnd(8)因为8 2^3 1 3这个值刚好只在 bit 3 上占位。如果结果不等于 0说明 bit 3 已经是 1也就是该像素被标记为云。反过来如果我们希望这个像素不是云就要求qa.bitwiseAnd(8).eq(0)。这个逻辑是整个去云 mask 的核心。很多老手写完这个函数后往往不再解释太多但新手卡住的地方也恰恰在这里到底用 4、8、16 还是 32这些数字对应哪一位下面我分别按新旧两套 QA 波段给出具体对照。2.2 旧版 pixel_qa 的位布局与去云代码旧版LANDSAT/LC08/C01/T1_SR的pixel_qa波段关键位定义是这样的bit 2数值 4卷云Cirrusbit 3数值 8云Cloudbit 4数值 16云影Cloud Shadowbit 5数值 32雪Snow网上很多旧代码喜欢写成var cloudMask qa.bitwiseAnd(4).eq(0);这是不严谨的。因为4对应的是 bit 2也就是卷云位不是普通云位。单独用 4 去掩膜只能过滤掉一部分薄卷云真正的厚云没有处理后期影像上会残留明显白块。我实测跑辽宁夏季影像时这种情况非常明显元数据显示某个 scene 云量只有 10%但局部仍然盖着厚厚的云层只用bitwiseAnd(4)根本去不干净。所以旧版代码我建议至少把云、卷云、云影三个位一起处理function maskL8C1SR(image) { var qa image.select(pixel_qa); var cloudBit 1 2; // 4卷云 var cloudBit2 1 3; // 8云 var shadowBit 1 4; // 16云影 var mask qa.bitwiseAnd(cloudBit).eq(0) .and(qa.bitwiseAnd(cloudBit2).eq(0)) .and(qa.bitwiseAnd(shadowBit).eq(0)); return image.updateMask(mask); }这段代码的思路很简单三个条件都满足也就是既不是卷云、也不是云、也不是云影的像素才保留下来。updateMask是 GEE 里把掩膜“贴”到影像上的方法掩膜为 1 的像素显示掩膜为 0 的变成透明。这算是去云代码里最标准的套路了。2.3 新版 QA_PIXEL 的位布局与去云代码新版LANDSAT/LC08/C02/T1_L2的QA_PIXEL在云、卷云、云影、雪这几个位的定义上和旧版pixel_qa是兼容的。也就是说 bit 2 是卷云、bit 3 是云、bit 4 是云影、bit 5 是雪。所以在 GEE 里写去云函数时代码骨架几乎一样function maskL8C2L2(image) { var qa image.select(QA_PIXEL); var cirrusBit 1 2; // 卷云 var cloudBit 1 3; // 云 var shadowBit 1 4; // 云影 var mask qa.bitwiseAnd(cirrusBit).eq(0) .and(qa.bitwiseAnd(cloudBit).eq(0)) .and(qa.bitwiseAnd(shadowBit).eq(0)); return image.updateMask(mask); }注意这里我加了卷云位。卷云对可见光波段影响不大但 C2 L2 的 QA 里把卷云单独列了出来在辽宁夏季高云多发的时候薄薄的卷云经常让 NDWI 出现很多假水斑。加上卷云位以后虽然可能稍微多去掉一些边缘像素但换来的指数噪声控制很值得。如果你想要更严格的掩膜还可以把 dilated cloudbit 1数值 2加进来。bit 1 表示云周边向外扩展一圈的区域它代表“云附近的像素可能也被云或云影影响”。加上这一位辽宁沿海区域受海雾影响时效果会好不少但也容易导致水体边界附近被过度剔除。我的经验是做水体指数用从严版本做土地利用分类用从宽版本实际根据需求调整。3. 辽宁省案例通吃脚本与 NDWI 应用3.1 研究区与影像筛选辽宁省夏季降水多、云量高而且境内有大伙房水库、辽河干流、盘锦湿地等水体是做 NDWI 和去云测试的天然实验场。更关键的是辽宁纬度高夏季太阳高度角也比较合适水体提取时不容易被地形阴影干扰。为了不依赖本地 ESRI Shapefile我这里用了一个大致覆盖辽宁省的矩形范围来演示。如果你要严格按省界处理可以把自己上传的辽宁省边界作为 geometry 替换掉。var liaoNing ee.Geometry.Polygon([ [[118.8, 38.6], [125.8, 38.6], [125.8, 43.5], [118.8, 43.5]] ]); Map.centerObject(liaoNing, 6);先加载 C2 T1_L2 数据。这里我选择 2021 年夏秋季节云量初筛设在 30% 以下。为什么不用 0%因为辽宁夏季整景云量 0% 的影像少得可怜筛完可能一个场景都不剩。CLOUD_COVER 这个属性只代表整个 scene 的云量不代表你研究区里面就一定晴空所以初始条件放到 20%-30% 都常见后续再用逐像元 QA 掩膜去兜底。var start 2021-06-01; var end 2021-09-30; var collectionC2 ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterBounds(liaoNing) .filterDate(start, end) .filter(ee.Filter.lt(CLOUD_COVER, 30)) .map(maskL8C2L2); print(C2 T1_L2 参与合成的影像数量, collectionC2.size());如果你打印出来的数量是 0 或者只有 1先检查日期和三月份的影像下载情况也可以把云量阈值放宽到 40%或者把时间跨度拉长。辽宁 6-9 月影像数量其实不算少正常打印出来应该有几十到上百条但去云之后能真正合成干净画面的往往就是另外一回事了。3.2 两套数据的去云对比与“翻车”现场做完 C2 集合我们再用同样的条件加载旧版 C1 SR 集合对比一遍。这里的对比目的不是证明谁更高级而是让你直观看到“去云代码必须匹配数据集”这个原则。如果旧代码直接拿到新数据上跑最常见的报错就是波段名不存在。var collectionC1 ee.ImageCollection(LANDSAT/LC08/C01/T1_SR) .filterBounds(liaoNing) .filterDate(start, end) .filter(ee.Filter.lt(CLOUD_COVER, 30)) .map(maskL8C1SR); print(C1 T1_SR 参与合成的影像数量, collectionC1.size());这里要提醒一句C1 毕竟是老数据源2021 年它仍然有影像覆盖用来做对比没问题。但如果你把日期改成最近一两年C1 集合大概率会出现“查不到新数据”的情况。这不是代码 bug而是产品已经停更。遇到这种情况直接换 C2 T1_L2 就好。为了看去云效果我习惯把每套集合分别做一个中值合成再加载真彩色。中值合成的逻辑是同一位置在多张影像中取中间值。它的好处在于即使个别影像残余了小块薄云或云影只要多张影像里该位置在大多数时候是晴空中值就会被“拉回”正常地表反射率水平。这也是为什么单景去云之后还要做合成的核心原因。var compositeC2 collectionC2.median(); Map.addLayer(compositeC2, {bands: [SR_B4, SR_B3, SR_B2], min: 0.01, max: 0.2}, C2 T1_L2 真彩色合成);从肉眼上看如果 C1 和 C2 结果都做得比较干净两者差异不大。真正的“翻车现场”往往出现在三个地方第一旧代码用了pixel_qa新数据里没有这个波段直接报错第二新代码用了QA_PIXEL老数据里没有也报错第三缩放公式没改导致真彩色整体偏亮或偏暗。这些我都踩过属于“表面没报错、结果已歪了”的隐形坑。3.3 用去除云的影像计算 NDWI去云的最终目的是为了得到可靠的地表参数。这里拿 GEE 里最常见的 NDWI 举例。NDWI 有两种常见公式一种是 McFeeters 提出的(Green - NIR) / (Green NIR)对水体敏感城市水体识别常用它另一种是 Gao 提出的(NIR - SWIR1) / (NIR SWIR1)更偏植被含水量和湿地监测。很多人一搜索“gee ndwi”很容易把两个公式混用。这没有绝对对错但要明确自己的目标。在 Landsat 8 C2 T1_L2 里Green 对应SR_B3NIR 对应SR_B5。我习惯先把 SR 波段做缩放再计算指数避免 offset 干扰function maskAndScaleC2(image) { var qa image.select(QA_PIXEL); var mask qa.bitwiseAnd(4).eq(0) .and(qa.bitwiseAnd(8).eq(0)) .and(qa.bitwiseAnd(16).eq(0)); var srBands [SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, SR_B7]; var sr image.select(srBands).multiply(0.0000275).add(-0.2); return image.updateMask(mask).addBands(sr, null, true); }于是在合成后直接计算 NDWIvar ndwiC2 compositeC2.normalizedDifference([SR_B3, SR_B5]).rename(NDWI); Map.addLayer(ndwiC2, {palette: [#a50026, #ffffbf, #313695]}, C2 L2 NDWI (Green-NIR));如果你用旧版 C1 SR只需要把波段名换成B3和B5缩放改成multiply(0.0001)其他逻辑一摸一样。这就是两代数据迁移的精髓表面看是两个数据源实际就是“换名字、换缩放、保持 QA 位思路”三件事。去云不干净的时候NDWI 上的云区往往会出现一片异常高值的“假水”云影区则会出现大片负值凹坑。做完 QA 掩膜再算 NDWI水体边界才会干净许多。4. 实操中的常见问题排查表4.1 为什么去云后仍有薄云或孔洞QA 波段并不是完美的。CFMask 算法对厚云、亮云识别率高但对薄云、半透明云、小碎云经常漏检。辽宁夏季空气湿度大低云和雾常常连成片QA 波段标记出来的云边界也不总是贴合实际。如果你发现去云后仍有白色残影优先做三件事第一把卷云位也加入 mask第二把CLOUD_COVER初筛阈值调低到 20% 甚至 15%第三改用中值合成避免单张影像的残余云直接参与最终结果。还有一个常见误区是看到云影残留就继续增大 shadow bit 的范围结果把大片正常水体也删掉了。C2 的 cloud shadow 检测本身已经带了一定的外扩膨胀再加大很容易误伤。遇到云影比较顽固的情况我更推荐直接剔除这个日期而不是在 mask 上“下猛药”。4.2 为什么新旧代码互相套用报错“报错”还算好的最怕的是不报错但结果错误。旧代码image.select(B4)在新数据集上运行GEE 可能会提示选择失败就算你改了波段名忘了改缩放公式代码也能跑完图片看起来颜色不对劲但你不会第一时间意识到问题是缩放偏移量引起的。所以排查顺序建议是先看波段名再看 QA 波段名再看缩放公式最后再看数据集合 ID。四样对上了代码基本就通。需要特别指出的是新版LANDSAT/LC08/C02/T1_L2里的QA_PIXEL和旧版 C1 里的pixel_qa虽然云、卷云、云影、雪的位定义一样但它们的 High-Level QA 细节存在差异。单纯迁移代码时不要想当然地认为其他位也完全一致尤其是涉及 Snow 和 Water 标记的业务场景一定要针对具体版本重新验证。4.3 冬季辽宁的雪、冰和云如何区分除了云冬季辽宁还有一个大麻烦雪和冰。CFMask 会把雪/冰单独标记在 bit 5数值 32。如果你做夏季水体提取雪不是问题但如果做冬季水体或长时间序列 NDWI积雪区域在绿波段和近红外波段的反射特征会非常接近水体很容易被 NDWI 误判成蓝色水体。我的建议是做水体提取时把 bit 5 也加入掩膜剔除条件写成qa.bitwiseAnd(32).eq(0)。这样至少能防止大批积雪被当成水面。不过要注意这会牺牲冬季雪盖上方的地表信息。如果研究目的里有雪盖范围评估就不要盲目屏蔽雪位而是单独生成一个 snow mask 来记录。冰面也是一样的道理在辽宁水库冬季结冰时NDWI 响应会变得不稳定最好结合时序和温度数据辅助判断。写在最后的一点实操体会我自己的习惯是不管用哪套数据去云函数永远不直接写死在主流程里而是单独拆成一个maskXxx()函数方便切换数据源时复用。在辽宁省这种云量复杂、冬夏季差异极大的地区没有哪个 mask 是“一次写死走天下”的一定要结合研究目标和季节特点去调整 bit。另外每次算指数前我都会随手print一下影像的直方图确认波段数值范围符合反射率语义再决定显示参数和阈值。这套流程看起来笨但正是这些细节让去云结果在一次次项目中保持稳定。