做植被遥感的人多半都有过类似的纠结：想找一个“看得清地块、又覆盖多年、还不用自己从头处理云噪声”的植被指数数据，比想象中难得多。MODIS的250米EVI适合看大趋势，可到了农田地块尺度，一个像元里树、田、路混在一起；Landsat的30米精度能看个大概，但年度单景受云影响严重，一条轨道上能凑出几帧晴空都不容易。这也是我下决心整理这套2019-2024年中国逐年10米分辨率最大值合成EVI数据集的原因——用Sentinel-2的10米表面反射率，按年做最大值合成，把全国连续6年的植被生长状态尽量“压缩”成每年一景、可以直接拿去分析的底图。
这篇内容不是泛泛介绍某个公开产品的宣传稿，而是把这套数据从设计思路、预处理决策、GEE生产管线到质量验证的完整复盘写出来。如果你准备自己做类似的高分辨率植被指数产品，或者想判断这份数据能不能用在自己的研究里，这篇应该能给你省下不少试错时间。
1. 为什么是“10米 + EVI + 最大值合成”这个组合
1.1 EVI比NDVI在10米尺度上到底强在哪
很长一段时间里，NDVI是植被遥感的主流指数，公式简单、计算成本低、对稀疏植被敏感。但它有一个绕不开的问题：当植被覆盖密度上来以后，NDVI会迅速饱和——一片郁闭度0.8的常绿阔叶林和一片郁闭度0.95的针叶林，NDVI可能都挤在0.90附近，几乎拉不开差距。10米像元虽然已经比MODIS那种250米“干净”很多，可遇到南方连片山地森林、高密度农田时，饱和问题依然存在。
EVI的设计初衷就是缓解这个饱和效应。它引入了蓝光波段，配合一个土壤调节项，公式是这样的：
EVI = 2.5 × (ρNIR − ρRed) / (ρNIR + 6 × ρRed − 7.5 × ρBlue + 1)
这个公式里，蓝光的作用是修正大气残留——特别是气溶胶对红光和近红外通道的影响。所以在有大气校正的表面反射率数据下，EVI在高植被覆盖区的梯度保持能力明显好于NDVI。
这一点在做年度最大值合成时尤其关键。最大值合成算法本身就会把指数推向一年中的“偏上限”位置，NDVI在这种场景下很容易整体进入平台期，而EVI还能保留出不同植被类型、不同长势之间的区分度。比如2019年云南某块常绿阔叶林，NDVI已经到0.92，2020年轻度退化后还是0.90，肉眼很难识别；换成EVI，可能是0.68对0.57，变化量级清楚得多。
1.2 为什么是Sentinel-2，为什么从2019年开始起算
10米分辨率这个指标基本锁死了数据源。Landsat系列是30米，MODIS最高250米，Planet类商业星座虽然能到3米，但获取历史长序列的成本不是一般课题组能承受的。真正免费且能回溯到几年前、还有连续表面反射率产品的，只有Sentinel-2。
Sentinel-2的10米波段是B2蓝光、B3绿光、B4红光、B8近红外，恰好覆盖EVI公式需要的三个波段。更重要的是，从2019年开始，A/B双星组网运行稳定，中国大部分地区的重访周期达到5天左右，生长季内能获得足够多的晴空观测。2018年前虽然也有数据，但SR集合的连续性、轨道覆盖的稳定性都差一些，硬要把起点推到2017年，质量检查会多出一堆麻烦。
所以我最终把时间范围定在2019—2024年，这6年正好是Sentinel-2数据质量相对稳定的完整时段，每一年都有足够多的有效观测支撑最大值合成。
1.3 年度最大值合成：一个容易被误用的方法
最大值合成（MVC）是植被指数处理里的老面孔了。它的思路很简单：在一段时间内，对每个像元取所有有效观测的最大值，因为云、云影、气溶胶通常让指数偏低，取最大值能近似恢复“晴空、最绿”的状态。
这套逻辑在月合成、季节合成里很成熟，但直接搬到年度尺度上就必须多留几个心眼。第一，云边界和薄云的残留会造成局部虚假高值，最大值合成反而会把这种噪声当成“生长最旺盛”的像元保留下来；第二，如果时间窗口覆盖了冬季，雪盖会让红光和近红外的反射特征完全走样，EVI数值没有生态学意义；第三，如果把不同气候区的所有月份混在一起取大值，“最大”就不再对应“最绿”，而可能对应“最噪”。
所以我在这个项目里坚持一个原则：预处理和季节窗口的选择，比“取max”这个动作本身重要得多。这也是下一章想详细展开的部分。
2. 生产前必须想清楚的三个预处理决策
2.1 只用表面反射率反演EVI，而不是TOA
做植被指数的人容易在图省事的时候直接用TOA反射率，但EVI这个指数对大气路径比较敏感，尤其是蓝光波段，受瑞利散射和气溶胶影响很大。TOA下算出来的EVI，会随着卫星过境时的太阳天顶角、大气水汽含量产生系统性变化，这种变化不是噪声，是偏差。
我一开始测试过直接用GEE里的COPERNICUS/S2（也就是L1C级TOA数据）跑EVI，结果2019年与2020年之间的数值差异里混入了一部分大气状态差异，没法放心用。后来全部改用COPERNICUS/S2_SR，也就是L2A级表面反射率产品，EVI公式里的蓝光修正项才算真正有物理意义。
需要提醒的是，S2_SR产品自带的波段里已经有SCL分类图、QA60波段，这些辅助信息在GEE里可以直接取用，省掉了自己写云检测算法的大量工作。
2.2 SCL波段的掩膜逻辑：只保留4、5、6
Sentinel-2 L2A自带的SCL波段是一张逐像元的场景分类图，数值含义如下：
| SCL值 | 含义 | 我的处理 |
|---|---|---|
| 0 | 无数据 | 挖掉 |
| 1 | 饱和/缺陷 | 挖掉 |
| 2 | 暗像元（地形阴影等） | 挖掉 |
| 3 | 云阴影 | 挖掉 |
| 4 | 植被 | 保留 |
| 5 | 非植被土地 | 保留 |
| 6 | 水体 | 保留 |
| 7 | 云低概率 | 挖掉 |
| 8 | 云中概率 | 挖掉 |
| 9 | 云高概率 | 挖掉 |
| 10 | 卷云 | 挖掉 |
| 11 | 雪 | 挖掉 |
很多教程会建议保留“云低概率”像元，觉得丢掉太可惜。但我的实测经验是：薄云边缘和卷云附近，SCL给出的低概率类别经常把EVI抬到不合理的程度。最大值合成本身会放大这种残留误差，所以我在掩膜策略上选择宁缺毋滥——只保留4、5、6三个类别，其余全排除。
还有一个隐蔽细节：SCL在GEE的S2_SR集合里，原始分辨率是20米或60米，而EVI计算用的B2、B4、B8都是10米。如果你不重采样就直接拿SCL做掩膜，掩膜边界和EVI像元会错位，结果在云边界处会出现一条条“彩带”伪影。我第一版跑出来的数据就栽过这个跟头。正确做法是把SCL用nearest邻域法重采样到B8的投影网格上再参与掩膜，分类图不能用双线性插值，否则类别数值会被插出毫无意义的中间值。
2.3 时间窗口：全年取max还是生长季取max
理论上“年度最大值”可以直接把1—12月所有合格观测丢进去取max。但我在几个典型区域做了试验，结论很明确：对北方地区来说，冬季太阳高度角低、地表可能有残雪，即使SCL已经把雪掩掉，低角度入射条件下的反射率噪声依然很大，而且冬季植被根本不活跃，EVI的年度最大值理论上不应出现在这个时段。可一旦观测条件差，噪声就会顶上来，形成假的最大值。
对南方常绿区来说，全年取max和生长季取max的结果差异不大，但如果要和北方统一语义，还是得用同一套规则。
我最后选择的处理窗口是每年4月1日到10月31日。这个窗口基本覆盖了我国绝大多数区域的植被生长季，同时避开了北方冬季的积雪和低太阳角噪声。说句实话，这个选择不是为了追求完美，而是为了在全中国范围内用一个统一、可解释的参数设定——窗口内取最大值，近似表示“这一年生长季里植被长到最旺盛时有多绿”。
3. 核心管线拆解：从Sentinel-2集合到逐年合成
3.1 为什么在GEE里做，而不是本地跑
6年、全国、10米分辨率的Sentinel-2原始数据，如果走本地下载流程，存储量是PB量级，个人工作站基本不用想。GEE的COPERNICUS/S2_SR集合把L2A产品在云端管理好，按需取用、边计算边丢弃中间数据，这才让“全国逐年10米EVI”从纸面方案变成能实际跑完的任务。
另外，GEE里做云掩膜和指数计算都是服务端并行，写代码的方式和本地栅格处理很不一样。本地你可能需要遍历几百个瓦片执行GDAL命令，GEE里则是对整个ImageCollection做map和reduce，思路更接近“声明式”——你说明要什么，服务器集群帮你去算。
这个方案当然有门槛，比如需要处理投影一致性、导出分块、配额限制，但这些都有成熟的对策，下面会逐一说明。
3.2 数据准备与逐年合成函数
GEE的JavaScript API里，整个逐年合成可以封装成一个函数。核心逻辑如下：
// 建议使用经过审核的全国行政边界，或者自己上传的高精度国界矢量 var china = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017') .filter(ee.Filter.eq('country_na', 'China')); var roi = china.geometry(); function buildAnnualComposite(year) { var start = ee.Date.fromYMD(year, 4, 1); var end = ee.Date.fromYMD(year, 10, 31); var s2 = ee.ImageCollection('COPERNICUS/S2_SR') .filterBounds(roi) .filterDate(start, end) .map(function(img) { // 将SCL重采样到10米网格 var scl = img.select('SCL') .reproject({ crs: img.select('B8').projection(), scale: 10 }); // 只保留植被/非植被/水体 var valid = scl.eq(4).or(scl.eq(5)).or(scl.eq(6)); var bad = scl.gte(7).or(scl.eq(3)).or(scl.eq(2)).or(scl.eq(0)).or(scl.eq(1)); var mask = valid.and(bad.not()); // 表面反射率计算EVI var evi = img.expression( '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', { 'NIR': img.select('B8'), 'RED': img.select('B4'), 'BLUE': img.select('B2') }).rename('EVI'); return evi.updateMask(mask); }); // 年度最大值合成 var composite = s2.reduce(ee.Reducer.max()).rename('EVI'); // 同时保留有效观测次数，用于质量控制 var count = s2.reduce(ee.Reducer.count()).rename('valid_count'); return ee.Image.cat([composite, count]) .set('year', year) .set('system:time_start', ee.Date.fromYMD(year, 7, 1)); } // 生成2019-2024年集合 var annual = ee.ImageCollection( ee.List.sequence(2019, 2024).map(function(y) { return buildAnnualComposite(ee.Number(y).toInt()); }) );这段代码里有几个要解释的决策点。
第一，我没有在map函数里把图像裁剪到中国边界，只做了filterBounds。原因是在GEE里，对每一景影像做clip到roi会额外增加几何计算开销，导出的阶段再统一处理边界更高效。
第二，SCL的reproject到img.select('B8').projection()，确保掩膜和EVI计算像元完全对齐。这个步骤在跨UTM带拼接时非常关键，少了它，云边界伪影就会回来找你。
第三，reduce之前没有额外设置统一的投影。这意味着集合内各景影像在参与reduce时，GEE会默认采用某种统一的网格进行重采样对齐。我在实践中的做法是，在导出阶段用固定的等积投影统一输出，这样既能保证跨年份可比，又避免输出影像在不同纬度上的面积失真。
3.3 年度合成与投影格网控制
全国范围的10米栅格，如果直接用EPSG:4326经纬度坐标输出，纬度越高每个像元代表的实际面积越小，10米这个“分辨率”会失真。我更推荐使用适合中国区域的双标准纬线等积投影，也就是Albers等积投影。
GEE的Export.image.toDrive的crs参数支持传入PROJ.4字符串，可以这样定义：
var albersChina = '+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'; Export.image.toDrive({ image: annual.filter(ee.Filter.eq('year', 2019)).first().select('EVI'), description: 'CHN_10m_EVI_2019', folder: 'evi_annual', scale: 10, crs: albersChina, region: roi, maxPixels: 1e13, fileFormat: 'GeoTIFF' });但这里必须提醒一个现实问题：全国一整幅10米EVI，单波段float32的GeoTIFF在几百GB量级，直接导出要么触发GEE的完成时间/内存限制，要么导出的文件大到没法正常上传下载。我的习惯是分块导出，而不是追求一个全国整幅的大TIFF。
3.4 分块导出与目录组织
我在本地准备了一个1°×1°的网格矢量（也可以用固定间隔生成Fishnet），每个网格作为一次导出的region。文件名里带上网格的参考经纬度，例如：
data/ 2019/ EVI_2019_N35E105.tif EVI_2019_N35E106.tif ... 2020/ EVI_2020_N35E105.tif ... qc/ valid_count_2019_N35E105.tif这样做的原因是，1°×1°的格子在中国范围大约覆盖100公里见方，单个GeoTIFF在10米分辨率下像素数量适中，无论后续是拼接、裁剪还是按行政区汇总，都能快速定位到对应瓦片。
如果你需要做公开发布，强烈建议把最终产品转成Cloud Optimized GeoTIFF（COG）格式，这样下游用户不需要下载整个大文件，就能在GIS软件里按需读取局部范围。
4. 质量验证与实战中踩过的坑
4.1 云残留导致的“假最大值”
第一版结果出来后，我在青海湖周边和云贵高原做抽查，发现不少像元EVI数值高得离谱，甚至超过1.2，而健康植被的EVI通常在0.9以下。逐景回查后确认，这些异常值来自云边缘的薄云残留区域——SCL虽然标出了高概率云、中概率云，但薄云边缘的过渡带依然漏了一部分进来。
针对这个问题，我做了两道加固。第一，在掩膜基础上再做一次形态学膨胀，把云层掩膜向周边扩大1到2个像元，宁可误杀也不放过过渡区；第二，在最终合成产品上设置EVI大于1.0的像元直接置为空值。这两个措施叠加后，假高值的数量明显下降。
后来我又遇到一个更隐蔽的情况：水体边缘的EVI值被算成了正数。原因是水体在蓝光波段反射率高，EVI公式里的蓝光修正项会让计算值在某些水色条件下变成正值。这也是我在掩膜规则里坚持保留SCL=6（水体）之后还要结合NDWI或必要时再查一遍的原因——水体掩膜不能只依赖SCL分类，浑浊水体、浅水区容易被分到非植被类。
4.2 早期年份观测次数不足造成的空洞
2019年虽然A/B双星都在轨，但当时的重访覆盖密度、数据回传效率都不如现在。加上东南山区春季云量大，个别像元一年里真正“晴空且通过SCL筛选”的有效观测可能只有个位数。如果这个像元唯一一次有效观测恰好发生在残云边缘，最大值合成后的结果自然不可信。
这就是我一直坚持导出valid_count图层的原因。使用这份数据时，我建议对有效观测次数少于3的像元保持高度警惕，在趋势分析里直接把低观测次数区域标记出来，而不是静默插值填补。数据生产可以追求覆盖完整，但算法上不应该掩盖“我们没有足够证据”这件事。
4.3 时间一致性检查的三板斧
每年合成产品单独看都能自圆其说，可放在一起就必须做时间一致性检查。我的常规操作是这三板斧。
第一，用已知物候区域做剖面曲线抽查。比如华北平原冬小麦区，正常年份的EVI时间序列应该是先升、平台、再降的单峰或双峰结构，如果哪一年出现突然跳变，就能反查是哪一步出了问题。
第二，与MODIS植被指数对比。把10米EVI聚合到250米网格，和MOD13Q1的EVI做散点图，相关系数拉到0.8以上才算通过。这里要注意的是，两个数据源的传感器波段设置、大气校正算法并不完全相同，绝对值会有系统偏移，所以对比时看的是空间分布的一致性，不是数值完全对齐。
第三，人工目视抽检。我选了30多个30km×30km的样区，每个样区逐年检查EVI灰度图与真彩色影像的叠加效果，重点看城镇边缘、大型水体边界、云残留区域有没有出现反常的“椒盐”噪声。这一步不能自动化替代，但也是最容易发现隐蔽问题的手段。
5. 这份数据能做什么、不能做什么
5.1 适合的场景
如果研究尺度是“地块级+年际变化”，这份数据基本能直接上手。比如农田长势的连续监测，10米分辨率能区分出不同田块之间的差异，这在250米MODIS数据上做不到；城市绿地与生态修复工程的时序评估，街道尺度的绿化变化也能看出来；还有森林扰动、退耕还林这类地表覆被变化的检测，年度最大值合成可以大幅减少云噪声干扰，作为分类和变化检测的特征输入也很合适。
做土地覆盖分类时，逐年EVI和原始多光谱波段、纹理特征叠加，能明显提升分类器对植被类型的区分能力，尤其是哪些在夏季光谱特征非常相似的类别。
5.2 需要小心的语义边界
必须说清楚的是，年度最大值EVI不等于“全年最高瞬时绿度”的准确测量，它更像一个“生长季内晴空条件下最绿状态的代理指标”。如果研究的问题依赖极端值或特定物候日期，比如想分析春季干旱导致的生长峰值下移，你得回到原始影像去做季节曲线，年度最大值会把这类信息平滑掉。
另外，EVI绝对值在不同数据源、不同产品版本之间存在系统性差异。跨数据源直接比较前，应该先做一致性校准，否则容易把传感器的差异误判为植被变化。
5.3 常用搭配
趋势分析方面，可以对6年逐年EVI做Sen's slope加上Mann-Kendall检验，识别显著变绿和显著变褐的区域。这里建议用年度EVI而不是生长季均值，因为最大值对物候差异不敏感，更稳定。
作物分析方面，可以结合物候期窗口提取EVI峰值时间、积分值和生长速率，再做种植面积或产量估算。比如东北春玉米的播期EVI低值窗口和抽雄期高值窗口，在10米产品上能比较清楚地呈现出来。
如果需要填补合成空洞，我建议保留“有效观测次数”图层，公开发布时也不做静默插值，而是让使用者在自己的分析流程中决定是否需要填补、用什么方法填补。这样既透明，也避免把“算法假设”误当成“数据事实”。
最后，还是想给准备复现这套流程的朋友一个最朴素的建议：不要一上来就追求把全国6年一次跑完。先拿一个省，跑一年，把云掩膜参数、SCL重采样网格、导出投影全部定下来，再做全量。我当时就是在云南先试了2019年，才把SCL掩膜规则和投影统一想清楚的。数据生产这种事，慢就是快。
