GEE实战：Sentinel-2 Level1C指定地点NDVI计算与云掩膜全流程

1. 项目缘起与整体设计思路

第一次接触遥感植被指数计算的朋友，十有八九是从NDVI入门的。原因很简单，归一化植被指数（Normalized Difference Vegetation Index）的公式就一行：近红外波段减去红光波段，再除以两者之和。但真正落到工程上，尤其是要针对某个指定地点、用Sentinel-2 Level1C数据、在云端批量跑出来，中间的门道远比公式本身多得多。

我这次要拆解的项目，核心目标很明确：在Google Earth Engine（后面统一简称GEE）平台上，对用户指定的一个地理坐标或区域，调用Sentinel-2 Level1C影像，完成NDVI的计算与可视化输出。听起来像是一个教学Demo，但实际做下来，它涉及影像筛选、云掩膜、大气校正的取舍、时间窗口设定、坐标系与分辨率对齐、结果导出等一系列环节。任何一个环节处理不当，最后出来的NDVI要么被云污染得没法看，要么数值范围异常，要么导出后跟GIS软件里对不上。

这个内容适合谁？如果你是遥感、地理信息、生态监测、农业估产方向的在校学生或刚入行的工程师，手头有GEE账号但还没系统跑通过一个完整流程，那这篇就是写给你的。如果你已经用过GEE但一直搞不清楚Level1C和Level2A到底该选哪个、云掩膜为什么总是不干净，那也能从里面找到答案。我会把每一步的“为什么这么做”讲透，而不是只丢一段代码让你复制。

整体设计上，我采用的思路是“先定区域、再定时间、后定处理链”。很多人习惯一上来就写ee.ImageCollection筛选，结果区域边界没定清楚，后面裁剪、导出全乱套。正确的顺序应该是：先把研究区（ROI）确定下来，可以是一个点加缓冲区，也可以是一个矢量边界；然后根据研究目标确定时间窗口，比如要监测某年生长季，就锁定6到9月；接着才是影像集合的筛选、云处理、指数计算；最后是可视化参数设定和导出。这个顺序不是拍脑袋定的，而是因为GEE是惰性计算（lazy evaluation）的，你写的每一行只是构建计算图，真正触发计算是在getInfo、Export或者地图渲染的时候。顺序清晰，调试起来才不痛苦。

还有一个关键决策：用Level1C还是Level2A。Sentinel-2的Level1C是经过几何精校正的大气顶层反射率（TOA），而Level2A是经过大气校正的地表反射率（BOA）。理论上做NDVI应该用Level2A，因为它去掉了大气散射和吸收的影响，数值更接近真实地表。但Level2A在GEE里的覆盖时间比Level1C晚，早期数据可能没有，而且有些区域Level2A会有条带或缺失。所以标题里指定Level1C，我理解是有意为之——要么是教学演示简化流程，要么是目标区域Level2A不可用。我的处理策略是：以Level1C为主，但在文中会说明如果Level2A可用，如何切换，以及两者NDVI的差异大概在什么量级。

2. 核心概念与关键细节解析

2.1 Sentinel-2波段与NDVI的对应关系

Sentinel-2 MSI传感器有13个波段，做NDVI只需要两个：B4（红光，中心波长约665nm）和B8（近红外，中心波长约842nm）。在GEE的COPERNICUS/S2数据集中，这两个波段的名称就是B4和B8。注意，GEE里还有COPERNICUS/S2_SR，那是Level2A的地表反射率产品，波段命名一样，但数值已经过大气校正。

NDVI的计算公式是：

NDVI = (NIR - RED) / (NIR + RED)

代入Sentinel-2就是：

NDVI = (B8 - B4) / (B8 + B4)

这个公式的输出范围理论上是-1到1。健康茂密植被的NDVI通常在0.6到0.9之间，稀疏植被在0.2到0.5，水体、云、雪通常是负值。这里有个容易踩的坑：Level1C的TOA反射率数值范围是0到10000（被放大了10000倍存储为整数），而Level2A的BOA反射率范围是0到10000甚至更大。如果你直接拿原始DN值算NDVI，公式本身是比值，放大倍数会约掉，所以结果不受影响。但如果你要做阈值判断，比如NDVI大于0.3算植被，那就要确保你算出来的确实是归一化后的值，而不是被缩放过的。

2.2 云掩膜：NDVI质量的生命线

云是光学遥感的天敌。一朵云飘过去，NDVI直接从0.8掉到-0.5，如果你不做云掩膜，时间序列曲线会像心电图一样乱跳。Sentinel-2 Level1C自带一个QA60波段，里面用位掩码标记了不透明云和卷云。GEE官方示例里常用的云掩膜函数就是基于QA60的：

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).divide(10000); }

这段代码的逻辑是：QA60是一个16位整数，第10位表示不透明云，第11位表示卷云。1 << 10就是1024，bitwiseAnd之后如果结果不为0，说明该位被置1，即有云。eq(0)就是取没有云的位置。最后divide(10000)是把DN值转成反射率，方便后续可视化。

但这里有个大坑：QA60的云检测并不完美，尤其是薄云和云边缘，经常漏检。我实测下来，在南方多云多雨地区，光靠QA60掩膜，生长季能用的影像可能只剩两三成。这时候可以考虑结合COPERNICUS/S2_CLOUD_PROBABILITY数据集做更精细的云概率掩膜，或者干脆用Level2A自带的SCL场景分类波段，里面直接标了云、云影、雪、水等类别。不过标题指定Level1C，我们就以QA60为主，在注意事项里会提替代方案。

2.3 时间窗口与合成策略

指定地点做NDVI，通常有两种需求：一种是看某一景影像的瞬时NDVI，另一种是看一段时间内的合成NDVI。前者简单，筛选出单景影像直接算；后者需要做时间合成，比如取一个月内所有可用影像的中值（median）或最大值（max）。中值合成能有效抑制云和异常值，最大值合成则倾向于保留植被最茂盛的状态，常用于生长季峰值监测。

时间窗口的设定要看研究目标。如果是农业监测，一般按旬或按月合成；如果是物候分析，可能需要逐日或逐周。这里没有标准答案，但有一个原则：窗口内至少要有3到5景可用影像，否则合成结果不稳定。我一般会先打印出筛选后的影像数量，如果太少，就放宽时间范围或者扩大区域。

2.4 坐标系与分辨率

Sentinel-2的空间分辨率是10米（B4和B8都是10米），GEE默认使用WGS84地理坐标系（EPSG:4326）进行显示和计算。但如果你要导出到本地做面积统计，最好投影到UTM或者Albers等面积守恒的投影，否则在高纬度地区面积会失真。GEE的Export.image.toDrive可以指定crs和scale参数，我通常会用crs: 'EPSG:32650'（对应UTM 50N）或者直接用image.projection()继承原影像投影。

还有一个细节：GEE里reduceRegion做区域统计时，默认的scale是影像的原始分辨率，但如果你不指定，它可能会用最大分辨率（比如10米），在大区域上计算量巨大甚至超时。所以做区域NDVI均值统计时，我一般会把scale设成30或100米，牺牲一点精度换速度。

3. 完整实操流程与代码实现

3.1 研究区定义：点加缓冲区还是矢量边界

第一步是定义ROI。最简单的方式是定义一个点，然后加一个缓冲区：

var point = ee.Geometry.Point([116.391, 39.907]); // 示例坐标 var roi = point.buffer(5000); // 5公里缓冲区

如果你有矢量边界（比如Shapefile或GeoJSON），可以先上传到GEE的Assets，然后ee.FeatureCollection('users/yourname/yourasset')加载。注意，GEE的buffer默认单位是米，但如果你在地理坐标系下操作，缓冲区会按球面计算，精度足够日常使用。

提示：缓冲区大小要结合研究目标。做城市绿地监测，1到2公里够了；做流域尺度，可能几十公里。缓冲区太大，计算量指数上升，而且会混入不同地类，NDVI均值意义不大。

3.2 影像集合筛选与云掩膜

接下来筛选Sentinel-2 Level1C影像：

var s2 = ee.ImageCollection('COPERNICUS/S2') .filterBounds(roi) .filterDate('2023-06-01', '2023-09-30') .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) .map(maskS2clouds);

这里我加了CLOUDY_PIXEL_PERCENTAGE小于20的过滤条件，这是影像级别的云量元数据，能快速剔除云量过高的整景影像。但注意，这个百分比是整个影像的，不是ROI内的。如果ROI很小，整景云量20%可能意味着ROI内全是云。所以更稳妥的做法是：先按云量粗筛，再在ROI内做像素级云掩膜，最后统计ROI内有效像素比例。

云掩膜函数前面已经给了，这里再补充一个细节：divide(10000)之后，影像的反射率范围是0到1左右。但有些像素可能因为云边缘或地形阴影出现异常高值，我一般会再加一个updateMask把反射率大于1的像素去掉：

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).divide(10000) .updateMask(image.select('B8').lt(1)); }

3.3 NDVI计算与合成

有了清洗后的影像集合，就可以逐景计算NDVI，然后做合成：

var addNDVI = function(image) { var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI'); return image.addBands(ndvi); }; var withNDVI = s2.map(addNDVI); var ndviComposite = withNDVI.select('NDVI').median().clip(roi);

normalizedDifference是GEE内置的归一化差值函数，等价于(B8-B4)/(B8+B4)，但更简洁且自动处理了除零问题。median()是中值合成，你也可以换成max()、mean()或qualityMosaic('NDVI')。qualityMosaic会按NDVI值从高到低排序，取每个像素上NDVI最高的那景影像的值，适合做生长季峰值合成，但会偏向于选择植被最茂盛的日期，可能掩盖物候变化。

注意：clip(roi)一定要加，否则后续导出或统计会包含ROI外的区域，浪费计算资源。

3.4 可视化参数设定

NDVI的可视化调色板很关键。我常用的是从棕色到深绿的渐变：

var ndviVis = { min: -0.2, max: 0.8, palette: ['#a50026', '#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a', '#66bd63', '#1a9850', '#006837'] }; Map.centerObject(roi, 12); Map.addLayer(ndviComposite, ndviVis, 'NDVI Median');

min和max的设定会影响视觉效果。如果你把min设成-1，max设成1，那大部分植被区域的NDVI在0.3到0.8之间，颜色对比就不够明显。我一般会根据研究区实际情况调整，比如城市区域min设-0.1，max设0.6；农田区域min设0，max设0.9。

3.5 区域统计与导出

如果你要统计ROI内的NDVI均值、最大值、标准差：

var stats = ndviComposite.reduceRegion({ reducer: ee.Reducer.mean().combine({ reducer2: ee.Reducer.max(), sharedInputs: true }).combine({ reducer2: ee.Reducer.stdDev(), sharedInputs: true }), geometry: roi, scale: 30, maxPixels: 1e9, bestEffort: true }); print('NDVI Stats:', stats);

scale: 30是为了加速，maxPixels: 1e9防止像素数超限，bestEffort: true让GEE自动调整分辨率。如果你要导出NDVI栅格到Google Drive：

Export.image.toDrive({ image: ndviComposite, description: 'NDVI_Export', folder: 'GEE_NDVI', region: roi, scale: 10, crs: 'EPSG:32650', maxPixels: 1e9 });

导出时scale: 10保持原始分辨率，crs指定投影。导出任务不会立即执行，需要到GEE的Tasks面板手动点击Run。

4. 常见问题与排查技巧实录

4.1 NDVI值异常：全为负或全为1

这是新手最常遇到的问题。原因通常有三个：一是波段选错了，比如把B4和B8写反了，或者用了B3和B4；二是云掩膜把有效像素全掩掉了，导致median()返回空值；三是divide(10000)没加，DN值直接算NDVI，虽然比值不变，但如果后续做了阈值判断就会出错。

排查方法：先打印影像集合的大小print(s2.size())，如果是0，说明筛选条件太严；再打印单景影像的波段列表print(s2.first().bandNames())，确认B4和B8存在；最后用Map.addLayer分别加载B4和B8，看看数值范围是否正常。

4.2 云掩膜后影像数量骤减

在南方地区，生长季云量高，QA60掩膜后可能只剩几景。这时候可以考虑：放宽CLOUDY_PIXEL_PERCENTAGE到30或40，但配合更精细的像素级掩膜；或者改用Level2A的SCL波段，它把云影也标出来了，掩膜更彻底；再或者用Sentinel-1的SAR数据做辅助，但那是另一个话题了。

4.3 导出后与GIS软件对不上

常见原因是投影和分辨率不一致。GEE默认导出可能是WGS84地理坐标系，而你在ArcGIS或QGIS里用的是UTM投影，叠加时会有偏移。解决办法是在导出时明确指定crs和scale，并确保GIS软件里也设置相同的投影。另外，GEE导出的GeoTIFF可能带有noData值，在GIS里要设置透明或忽略。

4.4 计算超时或内存不足

大区域、长时间序列、高分辨率三者叠加，很容易触发GEE的计算限制。我的经验是：先用小区域和短时间窗口测试代码逻辑，跑通后再逐步扩大；区域统计时把scale设大一点；导出时用maxPixels: 1e9和bestEffort: true；如果还是不行，就分块导出，比如按月份或按子区域拆分。

4.5 Level1C与Level2A的NDVI差异

我做过对比，同一地点同一天，Level1C的NDVI通常比Level2A高0.02到0.05，因为TOA反射率包含了大气散射的贡献，尤其是在蓝光波段，但NDVI只用红光和近红外，大气影响相对小一些。在植被茂密区，两者差异不大；在稀疏植被或水体边缘，差异会明显一些。如果你的研究对精度要求高，建议用Level2A；如果只是做趋势分析或教学演示，Level1C够用。

5. 进阶技巧与个人实操心得

5.1 用qualityMosaic做生长季峰值合成

如果你要提取一个区域在生长季的NDVI峰值，qualityMosaic('NDVI')比max()更好用。max()是逐像素取最大值，但不同像素的最大值可能来自不同日期，合成结果在空间上不连续。qualityMosaic是先按NDVI排序，然后取每个像素上NDVI最高的那景影像的完整波段，这样合成结果在时间上是一致的，适合做物候分析。

var ndviPeak = withNDVI.qualityMosaic('NDVI').select('NDVI').clip(roi);

5.2 加入时间序列图表

GEE的ui.Chart可以快速画出ROI内NDVI的时间序列：

var chart = ui.Chart.image.series({ imageCollection: withNDVI.select('NDVI'), region: roi, reducer: ee.Reducer.mean(), scale: 30 }).setOptions({ title: 'NDVI Time Series', vAxis: {title: 'NDVI'}, hAxis: {title: 'Date'} }); print(chart);

这个图表能直观看出植被的物候变化，比如返青期、峰值期、枯黄期。如果曲线有异常跳变，说明云掩膜没做好，需要回头检查。

5.3 批量导出多个月份的NDVI

如果你要按月导出NDVI，可以用循环：

var months = ee.List.sequence(6, 9); months.getInfo().forEach(function(m) { var start = ee.Date.fromYMD(2023, m, 1); var end = start.advance(1, 'month'); var monthly = s2.filterDate(start, end).map(addNDVI) .select('NDVI').median().clip(roi); Export.image.toDrive({ image: monthly, description: 'NDVI_2023_' + m, folder: 'GEE_NDVI', region: roi, scale: 10, crs: 'EPSG:32650', maxPixels: 1e9 }); });

注意，getInfo()是同步操作，会阻塞直到返回结果。如果月份很多，建议用ee.List的map配合Export，但导出任务本身是异步的，需要在Tasks面板手动运行。

5.4 个人踩坑记录

第一个坑：早期我用filterBounds(roi)之后直接median()，忘了clip(roi)，结果导出的影像比ROI大一圈，在GIS里叠加时才发现。第二个坑：QA60的位掩码我一开始写成了1 << 10和1 << 11，但实际测试发现有些影像的QA60波段是空的，后来查文档才知道，Sentinel-2 Level1C从2022年1月25日之后，QA60波段被移除了，改用MSK_CLDPRB和MSK_SNWPRB。所以如果你的时间窗口在2022年之后，云掩膜函数要改成基于MSK_CLDPRB的：

function maskS2clouds(image) { var cldProb = image.select('MSK_CLDPRB'); var mask = cldProb.lt(20); // 云概率小于20%视为无云 return image.updateMask(mask).divide(10000); }

这个变化很多人不知道，网上很多老教程还在用QA60，跑2023年数据就会报错。我建议你拿到数据后先print(s2.first().bandNames())看一眼，确认有哪些波段可用。

第三个坑：reduceRegion的scale参数。我一开始没设，默认用了10米，结果一个5公里缓冲区的统计跑了十几分钟还超时。后来改成30米，几秒钟就出结果。精度损失可以接受，因为NDVI均值对分辨率不敏感。

5.5 关于Level1C的再思考

标题指定Level1C，我猜测可能是教学场景或者目标区域Level2A缺失。但如果你在实际项目中，我强烈建议优先用Level2A（COPERNICUS/S2_SR），因为大气校正后的NDVI更准确，尤其是在气溶胶含量高的地区。切换方法很简单，把COPERNICUS/S2换成COPERNICUS/S2_SR，云掩膜改用SCL波段：

function maskS2sr(image) { var scl = image.select('SCL'); var mask = scl.neq(3).and(scl.neq(8)).and(scl.neq(9)).and(scl.neq(10)); return image.updateMask(mask).divide(10000); }

SCL里3是云影，8是中概率云，9是高概率云，10是卷云。这样掩膜比QA60干净得多。

最后再分享一个小技巧：如果你要对比不同年份的NDVI，一定要用相同的合成方法和时间窗口，否则物候差异会混入年际变化。我一般会固定用生长季（6到9月）的中值合成，这样年际可比性最好。