2019-2024中国逐年10米EVI合成数据集：生成技术与应用

做植被遥感的小伙伴应该都遇到过这种尴尬拿到的NDVI产品动不动就是250米、500米分辨率在平原地区凑合能用一旦到了破碎地块、山区沟壑一个像元里混着林地、农田、裸岩好几类地物指数值四不像根本没法支撑精细分析。这两年哨兵二号Sentinel-2数据普及开来10米分辨率的植被指数终于有了盼头但很多人卡在了“如何把一整年的影像处理成一张干净的逐年合成图”这一步——单景影像有云、有阴影时间序列上还有缺失处理不好做出的图全是噪点。这篇内容要聊的就是一套2019到2024年中国范围逐年10米分辨率最大值合成EVI数据集。简单说它把每年全国所有Sentinel-2影像按最大值合成MVC的思路压缩成一张无云遮挡、反映当年植被最旺盛状态的EVI影像空间分辨率10米逐年独立。对做农业长势监测、林业资源调查、生态修复评估、城市绿地变化研究的人来说这套数据可以直接拿来当底图或者做时间序列分析省掉自己从原始影像开始处理的巨大工作量。文章会把这套数据集的生成思路、核心处理流程、关键参数设置和踩过的坑完整梳理一遍不管是想直接用数据还是想自己动手做一套类似的产品都能找到参考。1. 为什么偏要攒一套逐年10米的EVI而不是继续用NDVI1.1 所有植被指数里为什么选EVI而不是NDVI先解决一个绕不开的问题EVIEnhanced Vegetation Index增强型植被指数和NDVI归一化植被指数到底差在哪。NDVI算得简单近红外减红光除以近红外加红光公式朴素、理解容易但它有两个老大难毛病。第一个是饱和问题当植被覆盖度很高、LAI叶面积指数超过3或者4的时候NDVI的增长曲线会变得非常平缓哪怕实际生物量还在明显增加NDVI数值却几乎不动了。对森林、茂密农田这种高覆盖区域NDVI跟实际长势的相关性会大打折扣。第二个是大气和土壤背景干扰红光通道对大气气溶胶敏感阴霾天或者有薄云时红光被散射NDVI值会被明显压低地表有裸露土壤时NDVI也容易受土壤亮度干扰。EVI就是冲着这两个问题去的。它的标准公式是EVI 2.5 × (NIR − Red) / (NIR 6 × Red − 7.5 × Blue 1)公式里多了一个蓝光波段用来做大气残留散射的校正因为蓝光对气溶胶的响应跟红光有相关性可以据此把大气影响消掉一部分。而系数6和7.5是经验调出来的权重用来压土壤背景的影响系数2.5是增益项把数值范围拉伸到适合显示的区间。L1的常数项则是土壤调节参数。这套设计让EVI在高植被覆盖区不容易饱和对大气噪声也相对鲁棒整体动态范围更大更适合做逐年对比——不同年份之间气象条件不一样大气状态不一样用EVI做时间序列分析比NDVI稳得多。当然EVI也不是没有缺点。它对输入数据的质量更挑剔蓝光波段本身信噪比不如红光和近红外如果大气校正做得不到位EVI反而会出现一些奇怪的异常值。所以实际生产EVI产品时大气校正这一步绝对不能省这也是很多自己取数算指数的朋友容易忽视的地方。1.2 10米分辨率意味着什么为什么说“每年一张全国图”很值钱10米分辨率放在全国尺度上看是一个相当奢侈的粒度。MODIS的植被指数产品是250米到500米Landsat的植被指数是30米而Sentinel-2的10米分辨率相当于把Landsat的一个像元细分成了9个。这意味着原来一个混合像元里可能包含的农田、道路、林缘、水体现在能分开识别了。在实际应用中这个差别很直观。比如做华北平原的冬小麦面积提取30米分辨率下一条3米宽的田埂和5米宽的机耕路会跟麦田混在一起提取出的面积普遍偏高10米分辨率下田块边界清晰多了田埂和道路能被单独分出来。再比如做南方丘陵地区的果园监测地块普遍小而碎30米影像里一个像元经常同时包含柑橘树、梯田坎和杂草指数被“平均”掉之后果园和草地的区分度很差但10米像元基本可以做到纯果园像元。每年一张、覆盖全国、没有云污染的10米EVI图层价值就在于“连续可比”。单独一年的影像好找但六年连续、逐年的最大值合成产品不好攒因为每一年的有效观测次数、云覆盖情况、物候偏移都不一样处理流程稍微不统一出来的数据就没法做趋势分析。1.3 2019到2024这个时间窗口的特殊意义2019年是Sentinel-2A和2B双星组网稳定运行、全球覆盖频次达到5天一回访的成熟期。再早的年份单星运行时期重访周期10天加上云遮挡部分地区全年有效观测次数可能不够支撑最大值合成。所以2019年作为起点是数据源条件决定的合理选择不是随意定的。而2019到2024这六年正好覆盖了“十四五”期间我国生态文明建设推进、双碳目标落地、耕地保护和粮食安全战略强化这几个关键阶段。农业上这六年里大豆扩种、玉米大豆复合带状种植等政策调整在空间上都有体现生态上天然林保护、退耕还林还草工程的成效评估也需要一个稳定的高分辨率本底。有了逐年可比的10米EVI就能回答很多“哪里变了”“变了多少”“趋势是向好还是恶化”这类问题。2. 生成一套逐年最大值合成EVI的完整技术链路2.1 数据源选型为什么只能选Sentinel-2以及它的局限性做10米分辨率植被指数目前来说Sentinel-2基本是唯一现实选择。Landsat系列最高30米达不到10米要求高分系列数据不开放获取没法做全国范围的长时间序列Planet、Spot这些商业数据价格高昂做六年全国覆盖成本根本不敢想。Sentinel-2凭借10米空间分辨率、5天重访、免费开放的策略成了这类任务的事实标准。但Sentinel-2本身也有一些绕不开的局限生产数据集时必须心里有数。第一它2015年才发射2019年之前的历史数据覆盖稳定性不足。第二国产化处理中L1C产品是大气表观反射率没有做大气校正必须先用Sen2Cor或者专门的大气校正模块处理否则算出来的EVI数值会偏高或偏低而且随季节和纬度变化误差模式还不一样。第三Sentinel-2的重访周期虽然是5天但在中国南方尤其是四川盆地、贵州、湖南这些多云地区一年里能拿到无云有效观测的次数可能只有十几次甚至更少这会直接影响最大值合成结果的可靠性。2.2 最大值合成MVC的原理以及它为什么能“去云”最大值合成英文叫Maximum Value Composite简称MVC是遥感时间序列处理里最经典、也最耐用的方法之一。核心逻辑特别简单在一段时间内对同一地理位置的所有有效观测按指数值取最大用这个最大值代表该时段内的植被状态。这个方法的理论基础是云、阴影、大气浑浊等因素只会降低植被指数值不会抬升它。云在可见光波段反射率高但在近红外波段也亮算出来NDVI或EVI会是很低甚至负值云阴影直接压暗所有波段指数也会掉下来气溶胶散射对NDVI是压制对EVI的影响相对小但也以降为主。所以取最大值等于在时间维度上做了一次“择优录取”——把每一天里那个最接近真实地表植被状态的值挑出来。MVC也不是万能的。它有一个先天弱点会把物候上的异常值误当成“最佳值”。比如某一年某块农田因为干旱出现早衰本该在7月达到峰值但7月的观测恰好被云污染了而6月有一次观测植被正处于旺盛期数值也很高MVC就会把那一次的值当成全年代表掩盖了干旱的影响。所以做逐年最大值合成时不能只丢一张全年影像了事最好结合物候窗口做分月合成再按植被生长季而非自然年去取年最大值。2.3 去云与大气校正的关键细节在实际处理Sentinel-2 L1C数据时第一个绕不开的步骤是大气校正和云检测。推荐做法是直接用官方的Sen2Cor插件把L1C处理成L2A地表反射率产品。Sen2Cor不仅做大气校正还会同时输出场景分类图SCL其中包含云、云影、雪、水体等类别可以拿来做后续的像素级云掩膜。这里有个很多新手容易踩的坑Sen2Cor对薄云的检测能力有限尤其对高海拔地区的卷云经常漏检。处理全国数据如果全靠Sen2Cor的云掩膜南方多云地区还是会出现不少残留云污染。更稳的做法是在Sen2Cor的基础上叠加一个基于蓝光波段阈值的补充检测。因为薄云在蓝光波段的反射率异常高而且云越厚蓝光越亮可以设定一个蓝光反射率上限我实测下来0.15到0.2这个区间比较合适超过阈值的像元直接标记为云再做一次膨胀dilation处理把云的边缘过渡带也抹掉。2.4 逐年合成的分月窗口策略直接拿全年所有影像做最大值合成从理论上讲是可行的但实际做出来会发现一些问题。比如北方的冬天地表被积雪覆盖时EVI会很低不会干扰年最大值但春季融雪季节土壤湿度大、部分地表还有残雪某些像元的蓝光波段反射异常若处理不当会产生刺眼的异常亮点。更稳妥的做法是先把一年分成若干个月份窗口对每个月的影像先做一次MVC合成得到12个月的月度EVI再从12个月度合成值里取年最大。这样做的优势有两个一是月度合成本身相当于又做了一次时域滤波能进一步压制残留的大气噪声二是为后续做物候分析保留了中间产品如果需要按生长季比如4月到10月取最大直接对月度数据做筛选就行不用重跑整个流程。2.5 技术平台选型与计算量预估全国范围、六年、10米分辨率这个数据量级单机处理不现实。Sentinel-2的10米波段全国一景原始数据大约在1GB左右全国每年大概需要4000到5000景影像覆盖考虑云量和重叠直接下载原始数据存储压力巨大。最合适的处理平台是Google Earth EngineGEE它的公共数据目录里已经归档了全球Sentinel-2 L2A数据可以在云端直接取用不需要本地存储原始影像而且有并行计算能力。GEE上做这件事的计算量也很可观。全国范围逐月MVC合成一次大约需要跑几十TB级别的像素计算个人用户免费额度勉强能撑住单年单次处理但六年连续处理建议用付费配额不然跑不完会被限流。如果是国内网络环境不方便用GEE也可以考虑用PIE-Engine或者AI Earth这类国产遥感云平台它们同样归档了Sentinel-2数据对国内用户来说访问更稳定。3. 实操流程从GEE上构建逐年EVI最大值合成数据集3.1 核心代码框架与流程分解整套流程在GEE上大概分为五步加载Sentinel-2 L2A影像集合、逐影像计算EVI、按月份分组做月度最大值合成、从月度合成取年最大、最后做全国分幅导出。下面给出一个可以直接参考的代码框架。// 1. 定义研究区范围此处示意为中国边界 var china ee.FeatureCollection(USDOS/LSIB_SIMPLE/2017) .filter(ee.Filter.eq(country_co, CH)); // 2. 加载Sentinel-2 L2A影像集合并限定时间范围 var s2 ee.ImageCollection(COPERNICUS/S2_SR) .filterBounds(china) .filterDate(2020-01-01, 2020-12-31) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 80)); // 3. 对每一景影像计算EVI并附加云掩膜 function addEVI(image) { var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(B8), RED: image.select(B4), BLUE: image.select(B2) }).rename(EVI); return image.addBands(evi).updateMask(image.select(SCL).gte(4).and(image.select(SCL).lte(7))); } var s2evi s2.map(addEVI); // 4. 按月份分组做月度最大值合成 var months ee.List.sequence(1, 12); var monthlyMax months.map(function(m) { var monthImage s2evi.filter(ee.Filter.calendarRange(m, m, month)) .select(EVI) .max() .set(month, m); return monthImage; }); // 5. 从12个月度合成中取年最大 var annualMax ee.ImageCollection(monthlyMax).reduce(ee.Reducer.max()) .rename(EVI_max) .clip(china);这段代码的思路很清晰通过SCL波段过滤掉云和阴影然后按月份把相同月份的影像叠在一起取EVI最大值最后把12个月的最大值再做一次max归约。注意代码里我用了CLOUDY_PIXEL_PERCENTAGE小于80这个过滤条件目的不是真的要全云景影像参与计算而是把明显废掉的景筛掉减少无效计算量真正的云掩膜靠的是SCL波段。3.2 参数选择的解释和调整建议这里有几个参数值得仔细抠一抠。SCL波段筛选条件SCL是Sen2Cor输出的场景分类取值范围中4代表植被5代表非植被6代表水体7代表未被云遮挡的土地。我把筛选条件设为gte(4).and(lte(7))相当于把云3、云影8、9、雪11全剔掉了。但实际运行时你会发现Sentinel-2 L2A产品里SCL等于0或1的像元也不少代表无数据和饱和这些在计算EVI时会变成空值所以最好再加上SCL.gt(1)的条件。月份窗口的分组方式上面的代码是严格按自然月分组的。但中国地域辽阔东部和西部的植被物候差异极大。同样是5月海南的植被已经进入旺盛期而东北的森林可能刚刚展叶。做全国尺度合成时分组窗口需要更灵活。可以把1到12月改成生长季窗口比如亚热带区域可以取4到10月温带区域取5到9月但这样会增加处理复杂度需要在代码里做区域分块。如果偷懒按自然年取最大得到的结果会更偏向夏秋季因为全国绝大多数区域植被峰值都出现在这个时段。CLOUDY_PIXEL_PERCENTAGE的阈值设到80意味着云量超过80%的整景影像直接丢弃。这个阈值看起来宽松但配合SCL像素级掩膜、月度最大值合成最终质量是有保障的。如果你把阈值设到20能大幅减少计算量但要小心一些多云地区可能连一景合格影像都留不下导致大片数据空洞。3.3 导出与后处理坐标系、缩放级别和瓦片化计算完成之后导出环节也有讲究。GEE上直接导出全国范围的单张影像文件会非常大而且10米分辨率下全国范围逐像元输出数据量轻松超过几十GB导出经常超时。更实用的做法是分省导出或者按经纬度网格分块导出GeoTIFF然后再本地拼接。导出时的坐标系建议使用UTM投影分带。全国范围横跨多个UTM带如果强行用Web Mercator或统一的Albers投影导出容易在投影变换中产生锯齿和错位。分带导出的优点是每个分带内的变形最小缺点是拼接时边缘需要做羽化处理不然会出现明显的接边痕迹。缩放级别上建议导出GSD地面采样距离为10米的原始分辨率而不要为了图方便直接按缩放级别导出瓦片。瓦片化虽然适合Web展示但做定量分析时精度会受重采样影响。3.4 质量检查样本点验证和空间完整性检查数据出来之后不能直接拿去用至少要做两轮检查。第一轮是数值范围检查。EVI在理论上取值范围是-1到1但实际地表上健康植被的EVI通常落在0.2到0.8之间。如果一张合成图里出现了大面积的负值或超过1的异常峰值说明处理过程有问题多半是云掩膜没有盖干净或者蓝光波段受残留气溶胶影响。最简单的检查办法是在GEE里直接看影像直方图确认绝大多数像元落在合理区间。第二轮是空间完整性检查。把每年的合成图叠加在行政边界上统计各省份的有效像元覆盖率。南方多云省份如果有效覆盖率低于80%那这一年这个区域的数据就不能直接用于变化检测需要补充其他数据源或者更换合成窗口。我个人的习惯是每年随机抽取50个地面样本点跟Landsat 30米EVI做一次交叉验证计算两个产品之间的相关系数。理论上Sentinel-2和Landsat的EVI在同一区域应该高度相关如果某个时间点突然出现相关性骤降那大概率是Sentinel-2那边有未剔除的云残留或者合成窗口内物候异常。4. 制作这套数据集时踩过的坑与解决实录4.1 南方多云地区的“数据黑洞”最开始做2020年数据时其他省份都很顺利唯独四川盆地和贵州一带合成的EVI图出现了大片空洞部分地区甚至整县都是空值。排查下来发现两个原因一是这些地区常年云量高全年有效观测次数确实少二是代码里把CLOUDY_PIXEL_PERCENTAGE设在80可这些区域的影像达到这个阈值的比例很高大量景被提前过滤掉了。解决办法有两步。第一步是把场景过滤阈值放宽到95让更多“半云半晴”的影像进入候选池靠SCL像素级掩膜去剔除云区这样至少能保住那些仅局部晴天的影像中的有用部分。第二步是调整合成窗口把按月合成改成按两月合成比如2到3月、4到5月这样相当于给多云地区增加了每个时间窗口内的有效观测样本量。代价是时间分辨率降低了但对于年度最大值这种产品来说影响不大。4.2 SCL分类错误导致的“消失的森林”第二次明显翻车是在2022年的数据上云南省一块区域的EVI值突然比周边低了一大截而且空间形态上正好是山脉走势。反复排查后发现问题出在SCL分类错误上——高海拔地区的针叶林在阴影下被Sen2Cor误分成了云影SCL8结果整片森林被当成云影遮掉了。这类错误的隐蔽性极强因为SCL是官方产品自带的分类结果很多人会默认信任它。实际上Sen2Cor在复杂地形下、尤其是高差大的山区误分类率并不低。解决办法是增加一个规则对被掩膜为云影但近红外反射率仍然较高说明下面有植被的像元做一次“赦免”处理把这类像元从云影掩膜中恢复出来。4.3 传感器差异引起的年度系统性偏移做时间序列分析时细心的朋友可能会发现2019年和2020年之间的EVI均值有一个小的跳变。这不一定是真实地表变化更可能是传感器原因。Sentinel-2A和2B两颗卫星的传感器虽然设计上一致但实际辐射响应有细微差异同一个地物在2A和2B影像上计算出的EVI可能会有0.01到0.03的系统性偏差。做六年连续数据集时这种系统性偏移会直接影响趋势分析的可信度。建议在正式的逐年合成之前先做一个交叉定标找一块大面积、辐射特征稳定的区域最好是沙漠或者大湖对比同期2A和2B的EVI值算出一个校正系数把2B的数据统一校正到2A的辐射水平上。这样能最大程度消除传感器间差异对趋势分析的影响。4.4 坐标系和投影的隐形坑还有一个不太容易发现但影响很大的坑在全国分幅导出时如果每一幅的边界是经纬度规则的但投影是UTM的那么相邻分幅在空间上会有一定的重叠或缝隙。拼接时如果不做严格的边缘匹配会出现明显的数据接缝。正确做法是在导出时给每幅图设置一个固定的缓冲区导出后再用gdal_merge或者ArcGIS的镶嵌工具选“最大值”作为重叠区像素取值规则这样接缝处的数值过渡会自然很多。5. 这套数据集到底能干什么、不能干什么5.1 最适合的应用场景农业领域可以做逐年耕地长势对比。华北平原冬小麦拔节到灌浆期的EVI峰值年际对比能直观反映出墒情差异和政策调控的效果。如果再结合物候窗口单独提取春季峰值还能辅助判断播种面积的变化。生态领域这套数据非常适合做森林恢复和退化评估。10米分辨率下人工林和天然林的边界、林窗、采伐迹地都看得比较清楚连续六年的EVI趋势可以帮助区分是长期恢复还是在退化。对南北两山的绿化工程、退耕还林区域的成效评估来说是很好的底图数据。城市生态领域可以分析城市绿地的年际变化。10米分辨率足够识别出城市公园里单棵大树级别的植被变化对海绵城市、公园城市建设的绩效评估有实际价值。5.2 不适合或不建议的用法不要用年最大EVI去做物候提取。年最大值只能反映当年植被最旺盛的那个瞬间体现不了返青期、枯黄期这些关键物候节点。要做物候得回到月度合成数据而不是年度数据。也不建议直接拿年最大EVI做不同土地覆盖类型间的横向比较。森林的年最大EVI和农田的年最大EVI数值体系不同农田在行播期和封垄期之间变化剧烈单靠年最大值很难准确区分作物类型和森林类型。对小区域的高精度总量估算比如估产或者生物量反演10米分辨率仍然不够而且EVI本身是光谱指数不等于生物量物理量需要做回归建模或辐射传输反演。5.3 和现有公开产品的配合使用建议用这套数据的时候最好跟MODIS MOD13Q1的250米NDVI/EVI产品做一次联合使用。MODIS数据时间长2000年至今虽然分辨率粗但可以用于做长时间趋势标定和异常年份识别。先用MODIS跑一遍全国趋势锁定变化显著的区域再用10米的逐年EVI数据做细节核查和边界勾绘这种“粗筛细查”的双尺度方案非常高效。6. 下一步的扩展方向6.1 从年合成扩展到季节合成年最大值合成解决了“一年一张干净图”的问题但很多研究需求其实是“关键物候期各来一张”。比如对水稻产区6月和9月的EVI分别对应移栽分蘖期和抽穗灌浆期前后的EVI差异能反映长势变化。把年度流程改成季节窗口合成输出的产品会更有分析价值。6.2 补充植被趋势斜率与显著性检验有了六年逐年数据一个很自然的衍生品是Theil-Sen趋势斜率和Mann-Kendall显著性检验。这两个方法对短时间序列特别友好不需要数据服从正态分布对异常值也有很强的鲁棒性。生成一张全国10米分辨率的EVI趋势图等于把六年变化压缩到一张图上对主管部门或者保护区管理者来说直观性远胜于六张单年图。6.3 考虑加入Landsat 8/9做时间回溯2019年之前的数据空白是数据源限制导致的但Landsat-8从2013年就有数据虽然分辨率是30米但通过数据融合方法比如时空自适应反射率融合模型可以把Landsat的长时间序列和Sentinel-2的高分辨率结合往前回溯到2013年左右。代价是精度略降但趋势连续性会明显改善。我把这套流程从数据获取到合成输出完整跑通之后最大的感受是10米分辨率并没有想象中那么遥远难的不是算法而是数据治理——从影像筛选、去云、大气校正到传感器定标每一步的小误差都会像滚雪球一样叠加到最终产品里。做这类长时序数据集真正的功夫都花在看不见的地方。最后分享一个小技巧做完第一版之后别急着马上发布把它跟已有的Landsat长时序数据做一次空间相关性回归如果相关系数低于0.7的区域占比超过10%说明数据里还有系统性噪声需要回头查一下合成窗口和云掩膜的参数。我一直在用这个标准来验收自己做的每一版产品算是一个低成本高回报的质量控制手段。