干遥感的人,十有八九遇到过这种场景:影像下载了几十个G,本地解压半天,处理到一半发现研究区有一大片被云盖住,旁边还拖着阴影;好不容易拼完了,想算个NDVI,结果发现数值不是负得离谱就是超出正常范围,折腾一通才发现是波段系数没处理。这些问题在Google Earth Engine(GEE)里其实是同一套流程里的三件事——影像加载、去云、波段系数转换。我刚转GEE那阵子也是把这三件事分开学的,后来把流程串起来才发现,很多所谓“翻车现场”的根源就是数据链路没打通。这篇文章就把这三个环节放在一起讲透,从原理到可直接跑的脚本,把每一步背后的逻辑说清楚,适合刚接触GEE的人,也适合被云困扰、搞不清DN值和反射率关系的同行参考。
1. 核心思路拆解:为什么把“加载、去云、转换”当成一个整体
1.1 遥感影像处理的“三关”
拿到任何一期卫星影像,处理链路其实可以拆成三关。
第一关是“请数据进门”。在GEE里也分两层:一是从上百个公开数据集中找到合适的数据源,二是把研究区范围内的影像筛选出来。这一步看上去是“查字典”,但它决定了后面所有工作的数据基础。
第二关是“剔除杂质”。光学遥感里最讨厌的杂质就是云和云影。云挡住的地方底下地物信息完全丢失,云影则会让地表光谱变得诡异。如果不去云,哪怕波段系数转换做得再精确,出来的结果也是错的——因为输入就带了脏数据。
第三关是“统一度量衡”。传感器记录的原始DN值是无量纲的整数,不能拿两景不同时间的影像直接做比较,更不能直接塞进NDVI这类波段运算公式里。必须把DN值转换成反射率或辐亮度这样的物理量,才能真正用于定量分析。
这三关有严格的先后顺序:先加载拿到影像,再去云去掉坏像元,最后做系数转换给“干净”的数据配上统一的物理刻度。很多人习惯只盯着其中一步,比如费劲学了各种去云算法,结果加载的波段没区分好,或者转换系数没做,最后分析照样翻车。把这作为整体来看,问题就简单多了。
1.2 为什么选Google Earth Engine而不是本地软件
我早年做遥感处理基本是ENVI加QGIS,流程是:下载影像、解压、逐景打开、拼接、裁剪、大气校正、波段运算。这套流程不是不行,但一遇到大范围、多时相的需求就非常痛苦。一个30景的Sentinel-2区域,下载恐怕就要一晚上,本地预处理脚本动不动报内存不足。
GEE的思路是完全反转的:影像数据本身就在服务器上,不用下载;处理脚本在云端并行执行,你面对的不是某一个文件,而是整个数据集合。我需要全北京市半年内的影像时,筛选加载往往就几秒钟。这种“数据不动代码动”的模式特别适合批处理和长时间序列分析。
但也要说句公道话,GEE不是万能的。如果你只是单景影像做高精度的辐射定标、或者做复杂的几何精校正,本地软件反而更可控。GEE的强项是快速浏览、批量预处理、指数计算和大范围空间统计。搞清楚边界,选型就不会出错。
1.3 这套方案适合谁、解决什么问题
这套“加载+去云+转换”三板斧能直接解决三类需求:土地覆盖分类前的数据准备、植被长势监测和指数反演、水体提取与动态监测。凡是需要用光学影像算NDVI、EVI、NDWI这些指数的人,都绕不开这三步。
以我自己的经验来说,做区域尺度的植被监测,这套流程处理完的SR反射率影像,直接往分类器或回归模型里扔,基本不会因为数据源的问题翻车。所以我建议任何刚开始用GEE的人,都先把这三步固化成一个模板,后续不管接入什么项目,都从这套模板开始改。先能在小范围内跑通再说加需求,别一上来就写几十行的复杂流程。
2. 影像加载:数据筛选与首屏可视化
2.1 GEE数据集的基本认识:Image、ImageCollection和常用地表数据
GEE里的影像有两种基本对象:Image是单景影像,ImageCollection是一个影像集合,可以理解成按时间排列的序列。日常处理很少直接用单景,一般是从ImageCollection里筛选一批,再合成或选代表。
常用的光学数据集我列一个表,方便大家按场景选择:
| 数据源 | 数据集ID | 空间分辨率 | 主要特点 |
|---|---|---|---|
| Sentinel-2 L2A | COPERNICUS/S2_SR | 10米(可见光/近红外) | 地表反射率产品,波段多,适合精细植被与土地利用分析 |
| Sentinel-2 L1C | COPERNICUS/S2 | 10米 | 大气表观反射率,需自行除10000 |
| Landsat 8/9 L2 | LANDSAT/LC08/C02/T1_L2、LANDSAT/LC09/C02/T1_L2 | 30米 | 长时间序列最佳选择,L2已经是SR产品 |
| MODIS | MODIS/061/MOD09GA | 500米 | 大尺度快速监测、植被指数产品成熟 |
这个表里最需要注意的是数据集版本。GEE这几年把不少数据集的波段命名统一改了,比如COPERNICUS/S2_SR的光学波段在新版里叫SR_B1到SR_B12,而老版叫B1到B12。很多在教程里能跑通的代码,你拿到自己账号上会报波段不存在,多半就是版本差异造成的。
2.2 一个能跑的加载脚本:筛选、裁剪与显示
先写一个最基础的加载脚本,不做任何高级处理,目标是让影像出现在地图上。
// 定义研究区:这里用矩形示意,实际项目中可以替换为矢量边界 var aoi = ee.Geometry.Rectangle([116.3, 39.8, 116.5, 40.0]); // 从Sentinel-2 L2A数据集中筛选 var s2 = ee.ImageCollection('COPERNICUS/S2_SR') .filterBounds(aoi) .filterDate('2022-05-01', '2022-10-01') .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)); // 取影像集合中的第一景 var image = s2.first().clip(aoi); // 加载到地图 Map.centerObject(aoi, 11); Map.addLayer(image.select(['SR_B4', 'SR_B3', 'SR_B2']), {min: 0, max: 3000}, 'S2真彩色(未缩放)');这段脚本做了四件事:定义研究区、筛选时空范围、按云量过滤、取一景裁剪显示。每一行都可以单独替换成自己的需求。CLOUDY_PIXEL_PERCENTAGE这个属性是元数据里的整景云量估计值,在加载阶段先粗筛一遍,能省不少后续工作量。
2.3 加载阶段最容易踩的坑
首先就是波段名。新版S2_SR数据用SR_B4,老版用B4。我习惯的做法是在筛选后先打印一次波段名:
print(s2.first().bandNames());看到真实波段名再往下写,比瞎猜靠谱得多。
其次是可视化参数的设置。S2_SR产品如果不做系数缩放,DN范围大概在0到10000以上,真彩色合成时min设0、max设3000左右一般能看个大概。如果你用的是我已经统一缩放过的影像,可视化的min/max就要修改(后面完整脚本会体现)。
最后是空间范围问题。filterBounds用的是影像元数据的空间覆盖范围,不等于影像裁剪到你的AOI。要真正裁剪,必须调用clip(aoi)。很多人加载出来的图还是整个条带,就是因为漏了这一步。这些坑单个看都不大,但串在一起就会让新手卡很久。
3. 去云处理:从质量波段到多时相合成
3.1 云为什么让定量分析“破功”
云和云影是光学遥感最大的天敌。云不仅直接遮蔽地表,云影还会让原本绿色的植被在影像上呈现出深暗甚至偏蓝褐色的混乱光谱。如果我们只做目视解译,云的干扰还能靠经验脑补;但一旦进入NDVI、分类模型这类定量流程,云区像元就会变成异常值,拉偏统计结果,甚至让整个影像的分类精度崩塌。
我在一次小麦长势监测中就遇到过典型情况:研究区连续两周多云,能用的几景影像里,有一景恰好一大块云压在研究区核心位置。当时去云函数写得粗糙,NDVI结果在那个区域出现了完全不合理的负值,差点影响结论。这件事之后我下决心把去云逻辑研究透彻,不再幻想运气好能躲开云。
3.2 单景去云:QA60位掩码和SCL分类
去云的低层逻辑是:每个影像都带有质量评估信息,记录每个像元是不是云、云影、冰雪。在Sentinel-2上,最传统的方式是用QA60波段做位掩码,取出其中第10位(不透明云)和第11位(卷云)的状态:
function maskS2clouds(image) { var qa = image.select('QA60'); var cloudBitMask = 1 << 10; // 不透明云 var cirrusBitMask = 1 << 11; // 卷云 var mask = qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); return image.updateMask(mask); }这段代码的思路是:如果对应位被置为1,说明这个像元是云,掩膜设为false,后续计算自动排除。这是GEE官方示例的标准写法,优点是简单可靠,缺点是位运算对新手不太友好,而且QA60对云影无能为力。
另一个更好理解的方式是用SCL场景分类波段。SCL波段是Sentinel-2预处理时生成的土地覆盖分类结果,其中类别3是云影,类别8、9、10分别是中概率云、高概率云和极大概率云。可以做如下掩膜:
function maskS2cloudsSCL(image) { var scl = image.select('SCL'); var cloudShadow = scl.eq(3); var cloud = scl.eq(8).or(scl.eq(9)).or(scl.eq(10)); var mask = cloudShadow.or(cloud).not(); return image.updateMask(mask); }SCL的优势很直观:不需要背位掩码表,而且能处理云影。不过SCL也有自己的毛病,比如容易把高山积雪和高亮建筑物误判为云,水体的深阴影有时也会被归为云影。所以选择哪种方式,要结合研究区地物情况。
3.3 多时相合成去云:median合成与“以时间换空间”
单景去云有一个天然局限:这一景里云盖住的地方永远缺数据,只能留空。如果研究区没有完全无云的时相,单景不管怎么去云都没用。
这时就该用多时相合成思路。原理很简单:同一个位置,在N期影像中大概率有一些期是干净的。我们把每一期的像元取中值,云通常是高亮尖峰,中值非常稳健地把尖峰“刺”掉。这就是median()合成:
var composite = s2 .filterBounds(aoi) .filterDate('2022-05-01', '2022-10-01') .map(maskS2clouds) .median();我用过很多次median()合成的效果,只要原始影像期数够多(一般5期以上),合成结果往往相当干净,而且不用写复杂的拼接逻辑。这也是为什么我在完整脚本里最终选择中值合成而不是单景去云:它把去云问题从“单帧清洁”变成了“多帧统计消除”,稳健性高一大截。
3.4 去云的边界:误判与调参
去云最怕的不是去不干净,而是误删。
第一个典型误判是雪和云。北方的冬春季影像里,SCL经常把雪区分成云,结果把好端端的地面变成空洞。这种情况可以在掩膜函数中加入NDSI(归一化雪指数)判断,把NDSI高的像元重新保留。
第二个误判是云影吃掉水体。特别是山区深水、峡谷阴影,SCL会认为这就是云影。你要是过度依赖自动掩膜,可能出现整个水库区域被挖空的现象。我一般在掩膜后目视对比一下原影像和掩膜结果,发现异常了就调整类别条件,或者临时增大buffer范围做局部修正。
第三个问题是薄云和碎云。薄云在QA60里可能识别不全,尤其是轻薄的卷云。如果研究区常受薄云影响,建议把蓝波段的反射率阈值加进去做二次过滤,因为蓝波段对大气散射最敏感,薄云会明显抬高蓝波段值。这类问题没有一劳永逸的答案,需要根据影像实际效果微调参数。去云这件事,本质上是经验活。
4. 波段系数转换:DN值、TOA与地表反射率
4.1 转换到底在转什么
先理清三个概念。DN值是传感器记录的原始“像素亮度”,不同传感器、不同时相里同一地物的DN值可能完全不同,它没有物理单位。TOA反射率是大气层顶的反射率,经过了辐射定标和太阳高度角校正,理论上可以跟其他传感器比较了,但没有去除大气的影响。地表反射率(SR)是经过大气校正后的产品,反映了地表真实反射特性。
生活里的类比:DN值是相机RAW文件里未经任何处理的原始像素值,TOA像是相机自动白平衡导出的JPG,SR则是你在Photoshop里认真校正过白平衡和色温后的成片。我们做定量遥感,目标当然是拿“成片”来分析。
很多人在GEE里算NDVI时犯的一个经典错误,就是直接拿不分幅、不缩放的原始影像套公式。结果两条几乎一样的影像,同一个小麦地块的NDVI一个0.5一个0.2,看着让人抓狂。问题就出在DN值没有可比性,必须先转换。
4.2 常见转换方式和GEE实现
不同数据集的转换方式不太一样,但原理都是基于元数据里的定标参数。Landat 8 Collection 2 Level-1产品中,可见光近红外波段的TOA反射率转换公式是:
ρ = DN × 0.0000275 - 0.2在GEE里写成:
var toa = image.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7']) .multiply(0.0000275).add(-0.2) .rename(['blue', 'green', 'red', 'nir', 'swir1', 'swir2']);Sentinel-2 L1C则简单得多,官方已经做了部分定标,反射率等于DN值除以10000:
var toa = image.divide(10000);而如果我们直接选用L2级别的SR产品,比如COPERNICUS/S2_SR或LANDSAT/LC08/C02/T1_L2,它们本身已经完成了大气校正和系数转换,是地表反射率,理论上不需要再做完整定标转换了。
这里我要特别提醒:SR产品虽然也是整数存储,但它的数值含义跟原始DN完全不同。S2_SR的光学波段实际反射率需要除以10000,也就是说你用代码跑出来一个像元值是2345,真实反射率就是0.2345。如果忘了除这一步,直接用原始值去算NDVI或做可视化,数值同样会乱掉。
4.3 转换之后才能做的指数运算
搞清楚系数转换后,很多指数计算就顺理成章了。NDVI的定义是近红外与红波段之差除以两者之和:
var ndvi = composite.normalizedDifference(['SR_B8', 'SR_B4']).rename('NDVI');在GEE里,normalizedDifference会自动处理波段顺序,第一个参数是被减数(近红外),第二个是减数(红波段)。做完系数转换后,NDVI的值域一般在-0.1到0.9之间,水体为负值,裸土和建设用地接近0,茂密植被在0.6以上。如果你算出来的NDVI大量出现超过1或低于-1的值,别怀疑公式,先回去检查波段系数转换。
除了NDVI,EVI的公式对大气和土壤背景做了更多校正,系数组合更复杂,但在GEE里也是一行:
var evi = composite.expression( '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', { 'NIR': composite.select('SR_B8'), 'RED': composite.select('SR_B4'), 'BLUE': composite.select('SR_B2') }).rename('EVI');这些公式本身不难,难的是确保输入波段是“同一种物理量”。所以我建议指数运算一律基于SR产品做,不要用L1C的TOA去算精确指数,省下来的大气校正步骤,最后都会以误差形式还给你。
5. 完整脚本:加载-去云-转换-指数计算一锅端
5.1 脚本结构设计
前面分别讲了三个环节,现在把它们组合成一个完整的、可直接运行的脚本。结构分成五段:定义AOI与时间范围、定义去云与系数转换函数、筛选加载并处理影像集合、多时相合成、计算并可视化指数。先跑通再优化,是我一直坚持的习惯。
这里要说明一下波段命名:新版GEE的S2_SR光学波段是SR_B1到SR_B12,老版是B1到B12。下面的脚本用SR_B前缀编写。如果你的账号是老版,只需把代码里的SR_B4替换成B4即可,其他逻辑不变。
5.2 完整代码
// ===== 1. 定义研究区与时间范围 ===== var aoi = ee.Geometry.Rectangle([116.3, 39.8, 116.5, 40.0]); var startDate = '2022-05-01'; var endDate = '2022-10-01'; // ===== 2. 定义去云函数(基于SCL) ===== function maskS2clouds(image) { var scl = image.select('SCL'); var cloudShadow = scl.eq(3); // 云影 var cloudHigh = scl.eq(8).or(scl.eq(9)).or(scl.eq(10)); // 中/高/全云 var mask = cloudShadow.or(cloudHigh).not(); return image.updateMask(mask); } // ===== 3. 加载并处理S2 SR影像集合 ===== var s2 = ee.ImageCollection('COPERNICUS/S2_SR') .filterBounds(aoi) .filterDate(startDate, endDate) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 30)) .map(maskS2clouds) .map(function(img) { // SR产品光学波段除以10000得到真实反射率,同时只保留光学波段 var sr = img.select('SR_B.*').divide(10000); return sr.copyProperties(img, ['system:time_start', 'system:index']); }); // ===== 4. 多时相中值合成 ===== var composite = s2.median().clip(aoi); // ===== 5. 计算指数 ===== var ndvi = composite.normalizedDifference(['SR_B8', 'SR_B4']).rename('NDVI'); var evi = composite.expression( '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', { 'NIR': composite.select('SR_B8'), 'RED': composite.select('SR_B4'), 'BLUE': composite.select('SR_B2') }).rename('EVI'); // ===== 6. 可视化 ===== Map.centerObject(aoi, 11); Map.addLayer(composite.select(['SR_B4', 'SR_B3', 'SR_B2']), {min: 0, max: 0.3, gamma: 1.2}, '真彩色合成(SR)'); Map.addLayer(ndvi, {min: -0.1, max: 0.9, palette: ['blue', 'white', 'green']}, 'NDVI'); Map.addLayer(evi, {min: -0.2, max: 1.0, palette: ['brown', 'white', 'darkgreen']}, 'EVI');整个脚本的核心是.map()函数。它对影像集合里的每一景影像依次执行去云和系数转换,然后再把这些处理后的影像集合用median()合成。这段流程最大的好处是:你不需要关心每一景影像各自的质量差异,集合级操作会自动把所有符合条件的像元融合起来。
5.3 结果验证与导出建议
脚本跑通之后,千万别看完颜色说“真好看”就完事了。要验证结果,我的习惯是打印几个关键统计值:
print('NDVI stats', ndvi.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: aoi, scale: 10, maxPixels: 1e9 }));如果NDVI统计最小值低于-1或最大值高于1,先回去检查是否漏了divide(10000),这是最容易被漏掉的步骤。另一个常用的验证方式是把掩膜前的真彩色图层和合成后的真彩色图层叠加上去,目视对比云的去除效果。如果合成图层上还有大片版块状异常颜色,大概率是SCL误判或期数太少。
导出方面,如果你需要把结果存到本地,可以加一句:
Export.image.toDrive({ image: ndvi, description: 'NDVI_result', scale: 10, region: aoi, maxPixels: 1e9 });导出前先想清楚分辨率。S2数据虽然原始分辨率是10米,但Export时设置20米或30米可以大幅减少导出时间和存储占用。研究尺度不需要10米精度时,别为了“全分辨率”死磕,导出会慢到让你怀疑人生。
6. 常见问题与排查速查表
6.1 问题速查表
这里把我在GEE实操中高频遇到的问题整理成表,方便随时对照:
| 问题现象 | 可能原因 | 解决办法 |
|---|---|---|
| 图层渲染全黑 | 可视化min/max设置不合理,或大量像元被掩膜 | 打印影像统计值,重新设置拉伸范围 |
| 真彩色效果灰蒙蒙 | 没做系数转换直接可视化 | SR产品除以10000,L1C除以10000后再显示 |
| NDVI大量超过1或为负 | 漏掉波段系数转换,或用了TOA数据的原始DN值 | 检查缩放;改用SR产品计算指数 |
| 波段名找不到 | GEE账号版本与教程不同,S2_SR的波段名有SR_B与B差异 | print(bandNames())确认,替换波段名 |
| 云没去干净 | QA60对薄云/卷云识别有限 | 改用SCL掩膜;增加云量阈值过滤 |
| 水体或雪被误删 | SCL将深水阴影、雪山误判为云 | 调整SCL类别条件,加入NDSI或NDWI判断 |
| 导出报错超时 | 区域太大或导出分辨率过高 | 提高scale、限制maxPixels,或分块导出 |
| 合成图有空洞 | 所有时相在该区域都被云覆盖 | 放宽云量过滤,拉长时间窗口,或用Radar数据补充 |
6.2 高频问题详解
第一个高频问题是NDVI超值域。我见过不少于五个新手把divide(10000)写在去云之前,结果QA60也被除了,后面掩膜逻辑一塌糊涂。顺序很重要:先updateMask做掩膜,再select波段做转换,最后才计算指数。
第二个高频问题是波段名报错。这也难怪,GEE更新数据集常常“偷偷”改命名规则。我自己从老版B2切到新版SR_B2时也懵了一阵。遇到这种问题别慌,print一下波段名,照着实名来写就行。所有教程里的代码都只是参考,真实环境以你的数据集为准。
第三个高频问题是云的残留形态。很多人去云后看合成图,发现大块云是没了,但图像上还有一道道“拖影”或模糊的亮暗斑块。这些往往是云影边缘、薄卷云残留,SCL和QA60都无法完美识别。我的经验是:对时间序列数据再做一次“每像元可用观测数”统计,如果某像元只有1-2期可用数据,合成结果大概率不稳定,可以通过ee.ImageCollection.count()检查有效观测数,把观测数太少的像元也掩掉。
6.3 几条实操心得
这套流程我重复跑了不下几十遍,踩过的坑比写出来的还要多几次。有两个心得特别想分享。
第一,调试时永远保留一个没有去云的图层作为对照。比如把原始真彩色放在图层管理器里,和去云后的结果来回切换。这样能快速判断是云没去掉,还是去得太狠把地物误删了。一上来就只盯着去云图层看,很容易误判。
第二,所有中间结果都打print。不要觉得打印占内存丢人,print(composite.bandNames())、print(ndvi.reduceRegion(...))这些日志在排查时比任何文档都管用。GEE的IDE就是你的实验室笔记本,保留实验记录能让你在跑更大的项目时少很多返工。
第三,不要在脚本一开始就写大而全的复杂函数。先用一个简单AOI、短时间窗口、一景影像验证逻辑,确认没问题后再扩展到多时相和大区域。这个习惯帮我省掉了无数“莫名其妙的上限超时”。
老实说,加载、去云、波段系数转换每个单拎出来都有人讲过,但把它们串成一条完整链路后,你会发现很多后续分析的“数据底子”问题都迎刃而解。我现在拿到任何区域的光学影像需求,都会先把这套三板斧跑一遍,确认数据干干净净、物理量统一,才敢放心往下一步的分类、回归或时序分析走。你也不妨把这套流程固化成一个自己的GEE模板,下次打开编辑器,直接从模板开始改AOI和时间窗口就行了。