这篇笔记是GEE学习笔记的第29篇。前面我写过Sentinel-2的单期NDVI、写过水体指数提取这次要解决一个特别高频的需求把一整年的Sentinel-2影像处理成一张高精度NDVI年均值数据并直接导出下载。所谓“高精度”在这里指的是使用Level-2A表面反射率产品、逐像元云掩膜和10米空间分辨率——不是拿DN值随手一算也不是用粗分辨率搞近似。适用场景很明确比如区域植被年际变化分析、退耕还林评估、农业长势监测、生态红线的本底摸查以及作为后续机器学习分类的输入特征。无论你是刚接触GEE的新手还是正在做时间序列分析的老手这份笔记都可以当作模板直接改。源码我会完整贴出来并把它拆开讲清楚每一步在做什么、为什么这样做。1. 项目背景与整体思路1.1 为什么选Sentinel-2做NDVI年均值这些年大家聊植被指数绕不开MODIS、Landsat、Sentinel-2这三个数据源。MODIS的NDVI产品时间序列长、处理方便但250米分辨率在县域尺度上会有明显的混合像元问题Landsat 30米分辨率可用但单星重访16天再加上云一年内干净影像数量并不算多。Sentinel-2A/B双星组网后5天重访10米分辨率还有Level-2A表面反射率成品非常适合做中高分辨率的年度植被监测。我自己的体感是在南方丘陵地区一块几平方公里的破碎耕地MODIS基本是一片混合值Sentinel-2能把田埂和小水塘都大致分开。所以做区域尺度的NDVI年均值我首选S2。下面是几个常见数据源的对比可以直观看到差距。数据源分辨率重访周期NDVI获取方式一年内有效影像MODIS250 m1-2天现成产品非常多Landsat 8/930 m16天自算少云区域通常够用Sentinel-210 m5天自算一般几十期这里说的“高精度”其实有两层含义。第一层是空间精度10米像元能把地物边界看得比较清楚。第二层是光谱精度我们用的是COPERNICUS/S2_SR地表反射率产品而不是顶层大气反射率COPERNICUS/S2这样可以减少大气散射对红光和近红外波段的影响NDVI本身也会更加可信。1.2 从影像集合到年均值的核心链路整体链路可以概括为四步选数据、洗数据、算指数、汇总导出。选数据是确定时间范围、空间范围、云量阈值洗数据是用QA60波段做像元级云/卷云掩膜算指数是归一化差分汇总导出是把一年内所有合格NDVI像元做均值再按区域裁剪导出。这个流程听起来简单真正决定成果质量的其实是第2步和第4步之间的配合。云掩膜不彻底年均值会被云污染像元拉低掩膜太激进有效观测数不足合成结果出现空洞。所以我在源码里特意加了一个有效观测数图层肉眼检查比事后找原因快得多。你看到红色区域的Valid observations值很低就要警觉了。2. 完整源码与逐段解析2.1 可直接运行的GEE源码先整理成一段能直接跑的完整源码。我默认把研究区设在北京中心点10公里缓冲区年份写的是2023年你只需要替换geometry和year两个变量就能用在别处。// GEE学习笔记 29一键下载 Sentinel-2 高精度 NDVI 年均值 // 平台Google Earth Engine Code EditorJavaScript API // 1. 设置研究区默认示例为北京中心10km缓冲区 // 如果有自己的矢量可以直接写var geometry table; var geometry ee.Geometry.Point([116.4074, 39.9042]).buffer(10000); // 2. 设置年份 var year 2023; var startDate ee.Date.fromYMD(year, 1, 1); var endDate ee.Date.fromYMD(year 1, 1, 1); // 3. Sentinel-2 Level-2A 表面反射率集合 var s2 ee.ImageCollection(COPERNICUS/S2_SR) .filterBounds(geometry) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)); // 4. 云/卷云掩膜基于 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.select([B8, B4]) .updateMask(mask) .divide(10000) .copyProperties(image, [system:time_start]); } // 5. 计算 NDVI function addNDVI(image) { var ndvi image.normalizedDifference([B8, B4]).rename(NDVI); return image.addBands(ndvi); } // 6. 一年内全部合格像元求平均 var ndviAnnual s2 .map(maskS2clouds) .map(addNDVI) .select(NDVI) .mean() .clip(geometry) .rename(NDVI_ year); // 有效像元数统计排查云掩膜过度或影像缺失 var validCount s2 .map(maskS2clouds) .map(addNDVI) .select(NDVI) .count(); // 7. 地图可视化 var visParam { min: -0.2, max: 1.0, palette: [#0000ff, #8b4513, #d2b48c, #ffff00, #008000] }; Map.centerObject(geometry, 11); Map.addLayer(ndviAnnual, visParam, NDVI annual mean year); Map.addLayer(validCount, {min: 1, max: 60, palette: [#ff0000, #ffff00, #00ff00]}, Valid observations); // 8. 导出到 Google Drive大范围、正式结果用这个 Export.image.toDrive({ image: ndviAnnual, description: S2_NDVI_Annual_Mean_ year, folder: GEE_Exports, fileNamePrefix: S2_NDVI_Annual_Mean_ year, region: geometry, scale: 10, maxPixels: 1e13 }); // 小区域也可以直接生成下载链接控制台 print 后点击 print(ndviAnnual.getDownloadURL({ name: S2_NDVI_Annual_Mean_ year, scale: 10, region: geometry }));运行方式很简单打开 Google Earth Engine Code Editor粘贴这段代码点Run。如果一切正常地图上会出现两个图层一个是带颜色的NDVI年均值一个是有效观测数。右侧Tasks面板里会出现S2_NDVI_Annual_Mean_2023的导出任务点击Run即可把GeoTIFF存到你的Google Drive里。2.2 核心函数逐段拆解源码里有几个关键点值得单独拿出来讲尤其是maskS2clouds这个函数很多人直接复制但看不懂一旦结果不对就不知道从哪里排查。先说影像集合筛选。filterBounds(geometry)是空间过滤只保留覆盖研究区的影像filterDate(startDate, endDate)是时间过滤取年初到次年年初filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20))是一个场景级云量过滤意思是整景影像的云量百分比不能超过20%。这个阈值可以根据地区调整后面我会展开讲。再说云掩膜。QA60波段是Sentinel-2的质量波段它用一个整数编码了很多位信息。1 10是1左移10位也就是二进制第10位bitwiseAnd(cloudBitMask)就是把这个位提取出来。如果结果是0说明这个像元没有被云覆盖如果结果是1说明有云。卷云同理用第11位。最后把干净像元保留云和卷云像元变成空值。divide(10000)是因为Sentinel-2 Level-2A产品的表面反射率是以10000倍存储的整数除以10000后得到0到1之间的反射率。这个步骤不是必须的因为normalizedDifference在求比值时上下都同一个尺度最终NDVI不变。但保留这个步骤更规范也方便你后续输出反射率做别的分析。addNDVI里的核心就一行normalizedDifference([B8, B4])公式是(NIR - Red) / (NIR Red)对应Sentinel-2就是(B8 - B4) / (B8 B4)。注意B8是近红外B4是红波段顺序一旦写反NDVI会变成负值。3. 关键参数精讲精度从何而来3.1 场景级云量阈值与像元级QA60掩膜很多初学者只做filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20))以为这样就处理干净了其实这是两个层次的事情。场景级云量是一个全局统计量代表整景影像里有百分之多少的云而QA60掩膜是在像元级别把影像里残留的云和卷云像元单独抠掉。比如一景影像云量10%看起来不高但云可能正好集中在你的研究区上空。如果你不做像元级掩膜那这10%的云区域就会参与NDVI平均直接把结果污染。反过来QA60也不是万能的它对薄云和云阴影的识别不完全冬季高海拔地区的雪也可能误判。这里我给一个经验值干燥地区场景级云量阈值可以放到20到30湿润地区建议10到15。阈值太小容易把整年影像筛没了阈值太大残留云多。如果你所在区域地形复杂、阴影和雪比较麻烦可以考虑用S2的SCL场景分类波段做二次过滤。SCL中的低值通常代表阴影或云但这会让代码更复杂暂不展开。3.2 均值合成、中位数合成与最大值合成怎么选“年均值”这三个字听着简单但合成策略并不唯一。我用同一景区域分别跑过mean、median、max三种方法结果差异其实很大。均值反映全年平均绿度状态对农作物的多茬种植和常绿林地比较友好中位数更稳健能消除残留云和极端值的影响最大值反映的是年内绿度峰值常用于估算植被覆盖潜力。这个任务既然标题叫年均值我默认用mean()。如果你只是想初步看一眼区域植被生产力中位数其实更稳。我个人的习惯是正式跑趋势之前至少同时算一份mean和一份median看看结论是否稳健。如果两者趋势方向一致说明你不是被某个季节的极端影像带偏的。3.3 分辨率、坐标系与导出设置scale参数直接决定导出像元大小。源码里写的10米对应Sentinel-2原始分辨率是“高精度”的关键。但需要认清一个现实10米分辨率是以像素数量成平方增长的。你的研究区如果是一个县10米导出大约会有几万乘几万像元如果是整个省文件几个GB甚至更大导出任务很容易失败。我建议根据范围调整研究区大小scale参数说明10km缓冲区10 m完全无压力典型县域10-20 m可导出注意磁盘空间地级市/省20-30 m推荐平衡精度和体积全国范围100-250 m不建议用10m硬跑坐标系方面源码里我没有显式指定crs这样GEE会使用影像默认投影通常是原始S2瓦片所在的UTM投影。好处是输出像元与原始数据对齐缺点是如果你一个大范围跨多个UTM分带默认投影可能不是你想要的。如果你需要统一的WGS84经纬度坐标可以在导出设置里加crs: EPSG:4326。但注意4326是地理坐标系跨带变形会比UTM明显所以常规做法还是默认投影本地再用GIS转。maxPixels是导出任务的兜底限制。10米分辨率大范围跑的时候一定要把这个值调大否则会报Error: Number of pixels too high。但也不要盲目填一个天文数字先算一下研究区像元数量再乘以1.2倍左右比较稳妥。4. 实操过程从零到一张年均NDVI图4.1 三步快速跑通默认例子拿到代码最直接的方式就是先跑通默认参数不要一上来就换成自己的陌生区域。具体操作分三步。第一步打开GEE Code Editor新建脚本粘贴源码点击Run。如果代码没有语法错误右侧Console会输出NDVI影像的信息和下载链接地图也会自动加载图层。第二步检查地图上的NDVI年均值图层。植被区域应该呈绿色裸地呈黄色或棕色水体和建设用地偏蓝。这个直观检查能提前暴露波段组合或云掩膜的问题。再看Valid observations图层如果一片区域全是红色说明那一年有效观测数非常少后面即使导出结果也没有可信度。第三步点开Tasks面板找到导出任务点击右侧的Run按钮。第一次使用时GEE会要求你授权Google Drive跟着提示走即可。等待任务完成去你的Google Drive里找到GEE_Exports文件夹下载GeoTIFF到本地。下载后可以在QGIS或ArcGIS里打开查看波段属性坐标参考是否正确。4.2 把自己的研究区和年份替换进去默认示例始终只是示例实际项目里肯定要换成自己的区域。最快的方式是用GEE左侧的绘图工具画一个多边形画完后在代码里把开头两句替换掉。如果你画完多边形后GEE生成的是一个变量table代码就写成var geometry table.geometry();如果你是通过Assets上传了Shapefile或GeoJSON代码写成var geometry ee.FeatureCollection(projects/你的用户名/assets/你的矢量路径).geometry();还有一种常见做法是从其他公开数据集里筛出行政边界像我源码里演示的其实是一行注释掉的写法// var geometry ee.FeatureCollection(USDOS/LSIB_SIMPLE/2017) // .filter(ee.Filter.eq(country_na, China)).geometry();注意导出时region参数需要的是一个Geometry对象而不是FeatureCollection。如果控制台打印出来是FeatureCollection记得加.geometry()否则导出任务会报 invalid region。年份修改更简单直接把var year 2023;改成想要的年份。年份区间是自动计算的startDate是当年1月1日endDate是下一年1月1日所以不需要手动改日期。4.3 下载链接方式与Drive导出方式取舍源码里我同时加入了getDownloadURL()它会在Console输出一个可点击的下载链接。这个方案看似方便但不是对所有区域都适用。GEE对于直接下载链接有像素数量限制小范围比如十几平方公里以内没问题一旦研究区变大或者你想要10米高精度链接会直接报错。我的建议是正式项目一律用Export.image.toDrive因为任务式导出没有交互式链接那么严格的像素限制而且可以异步排队不容易断掉。直接下载链接更适合快速验证、临时拿一小块数据。如果你一定要大范围直接用链接下载可以把大范围切成格网循环生成小块的URL。但那样会让代码复杂不少不如直接用Export然后等Drive同步。5. 常见问题与排查技巧实录5.1 导出任务成功但影像全黑或全空白这个现象我遇到不止一次最常见的原因是图层可视化参数和实际数值范围不匹配。比如你把visParam的min设成了0max设成了1但如果年均NDVI在某些高海拔区域普遍偏低低于0就全部显示成蓝色看起来像“空白”。这种情况下先用Map.addLayer配合一个更宽的范围比如min: -1, max: 1试试而不是急着怀疑导出数据。另一个原因是geometry类型不对。Export.image.toDrive的region参数只接受Geometry如果你传入了FeatureCollectionGEE会在任务提交阶段报Invalid region但有时候你看到的错误信息不够明确。我在排查时习惯先加一行print(geometry)看看类型是什么再做转换。还有一种是几何边界太复杂。如果矢量边界有几十万个顶点GEE处理这个边界就需要大量计算结果任务虽然能导出但研究区边缘可能会出现奇怪的黑边或空洞。对复杂边界可以先用geometry.simplify()做简化或者对影像先clip再selfMask()。5.2 NDVI数值范围异常或出现大量空值如果你看到NDVI整片都是负值大概率是波段顺序写反了把normalizedDifference([B8, B4])写成了[B4, B8]。这在所有GEE新手翻车原因里排前三我自己也干过。另一个容易踩的坑是把B8A当成了B8Sentinel-2的窄近红外波段是B8A但NDVI公式里用的是宽波段B8两者数值有差异导出后与别人产品对比时也会不一致。出现大量空值常见原因是云掩膜过于激进。比如你在湿润地区把场景级云量阈值设成5一年下来可能只剩十来景影像再抠掉云和卷云很多像元在所有影像里都是被掩膜状态合成结果自然为空。应对方式是放宽场景级云量阈值或者改用中位数合成。如果只是背景区域空值可以在导出前加一个unmask(-9999)但要注意这样会把没有观测值的区域变成-9999而不是真正的NODATA后续分析前一定要处理。5.3 任务排队时间太长或内存溢出的应对GEE的导出任务排队时间取决于当前服务器负载尤其是下午时段全球用户都在跑任务排队一两个小时很正常。如果等了很久还是Queued不要反复提交同一个任务先把代码停掉等一会儿再看。小任务如果也长时间排队可以换Export.image.toCloudStorage但需要配置云存储不是所有用户都有条件。内存溢出通常是大范围加上10米分辨率导致的。之前有个朋友想导出一个省份的NDVI年均值scale直接写10任务提交后很快报User memory limit exceeded。我的建议是先把maxPixels提高然后看是否还报错还报错就降分辨率比如20或30米。大范围高分辨率的需求最好分块处理不要指望一台浏览器就能合成全国10米数据这不是代码能力问题是地球引擎的资源限制。6. 扩展思路与个人经验6.1 改造成逐月或逐季节NDVI合成年均值是基础实际项目里很多需求要的是物候特征比如生长季峰值、峰值出现时间、春季NDVI增速。把源码改造成逐月合成其实不难用一个map遍历月份即可var months ee.List.sequence(1, 12); var monthlyMean ee.ImageCollection(months.map(function(m) { var start ee.Date.fromYMD(year, m, 1); var end start.advance(1, month); var monthly s2.filterDate(start, end) .map(maskS2clouds) .map(addNDVI) .select(NDVI) .mean() .set(month, m); return monthly; }));得到的结果是一个以“月”为单位的ImageCollection你可以直接导出每一张或者用monthlyMean.toBands()合成一个12波段的影像。后面做物候分析时这个12波段的堆叠影像会非常方便。6.2 用NDVI年均值快速估算FVC植被覆盖度很多研究需要植被覆盖度而不仅仅是一个NDVI值。最常用的遥感估算方式是二分法模型公式是FVC (NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil)其中NDVI_soil是纯裸土NDVINDVI_veg是纯植被NDVI。这两个阈值有地区差异干旱区通常取0.05和0.70湿润地区可以取0.10和0.85。代码实现非常简短var fvc ndviAnnual.subtract(0.10).divide(0.85 - 0.10).clamp(0, 1).rename(FVC);注意不要直接把0和1以外的值扔掉有些区域NDVI比土壤值还低比如水体会被clamp成0高密度植被NDVI超过植被阈值会被clamp成1。这样做出来的FVC在0到1之间后续统计覆盖度百分比会方便很多。6.3 构建多年NDVI序列衔接深度学习和地物分类如果你想做2019到2024年连续6年的植被趋势可以写一个循环把每年的ndviAnnual放到一个列表里最后用toBands()叠加成多波段影像。波段名类似NDVI_2019、NDVI_2020这个堆叠结果可以直接作为随机森林或深度学习的输入特征。我自己的经验是单靠年均NDVI做地物分类区分草地和灌丛的效果往往一般如果加上年内最大值、峰值时间和夏季均值分类精度会有明显提升。这也是为什么我建议先把源码跑通后续扩展的空间真的很大。目前GEE也支持把影像导出成TFRecord格式配合TensorFlow做深度学习前提是你对数据尺度、样本量有清醒的规划。刚开始不用追求太复杂的模型先把手头NDVI时序数据整理干净比任何算法都重要。最后讲一点我自己踩出来的经验。年均NDVI并不是年份越长越适合直接比较不同年份的有效观测数不同云覆盖差异会造成合成结果系统偏低。所以每次导出一张年均值图我都会同时导出有效观测数跑趋势时把观测数小于10的像元剔除。这个习惯帮我避开了很多假信号。源码里的Valid observations图层就是干这个的别嫌它占存储建议一并导出。祝大家跑图顺利。