做了几年遥感数据处理来来回回经手的全国尺度植被数据不在少数。这次整理这份全国省市县NDVI数据2000-2024算是把之前散在各处的活儿集中做了一遍。写这篇东西主要是把数据从原始影像到最后省市级产品的完整路径、踩过的坑、以及一些能直接用的处理思路记录下来给同样需要做长时序植被分析的朋友一点参考。1. 这份全国省市县NDVI数据到底是什么1.1 一句话讲清楚NDVI和这个产品NDVI全称归一化植被指数本质上是近红外波段和红光波段反射率的差值归一化。公式很简单(NIR - Red) / (NIR Red)取值区间大致在-1到1之间。绿色植被在近红外波段反射率高、红光波段吸收强所以NDVI值越高通常代表植被覆盖越好、长势越旺盛。裸土、水体、冰雪这类地物NDVI一般偏低甚至为负。这套数据集覆盖的是2000年到2024年这25年时间跨度刚好跨过了几个重要的植被演变阶段。空间范围是全国行政层级细化到省、市、县三级。我做的产品基本都是月度合成和年度合成两种时间粒度因为它面向的绝大多数应用比如植被长势监测、退耕还林评估、荒漠化趋势分析都是以月、季、年为单位的并没有必要去做日尺度。数据的空间分辨率省级和市级我主要用的是250米县级因为区域小、边界细碎如果还用250米小县的像元数会非常少所以部分地区提升到了100米。这个分辨率选择不是拍脑袋定的纯粹是基于数据源可获取性和实际应用需求折中出来的。1.2 谁在用这份数据、能解决什么问题这三年陆陆续续找我要过类似数据的主要有几类人。第一类是高校和科研院所的师生研究方向集中在土地利用变化、生态屏障区植被恢复评估、历史干旱事件重建。他们最常问的一句话是你们有没有某个县过去二十年每年夏天的NDVI均值因为要做回归分析缺一年都不行。第二类是自然资源和农业相关部门的项目团队用来做植被覆盖度分级、退耕还林地块的成效核查、草原生产力的粗略估算。他们不太关心NDVI的算法细节更关心数据能不能按行政边界一键裁好、年份能不能对齐、字段名是否统一。第三类是一些做碳汇、生态价值核算的咨询公司他们要把NDVI作为绿度指标套进模型里做生态资产测算。这类用户特别看重数据的连续性和可比性——如果某一年的数据因为云污染或者传感器差异出现系统性的偏低整个核算结果都会出问题。这三类用户的共同点是他们都希望拿到的是能直接入库分析的成品而不是一堆还没做大气校正、还没做去云处理的原始影像。所以他们需要的不是原始产品而是经过了时间合成、数据清洗、按行政区划裁剪、统一命名的一套准成品。1.3 我踩过的第一个坑坐标系和分辨率不统一这个数据集最早一批交付的时候我差点在坐标系上栽跟头。某平台的产品默认是正弦投影而用户拿到的行政区划边界通常是WGS84经纬度坐标有一个省级区域的客户直接把数据叠加到国土三调的底图上看边界怎么都对不上。后来统一处理为WGS84地理坐标系投影坐标则按各区域所在分带重新转成适合当地的高斯-克吕格投影或者UTM投影。这里有个经验存档阶段一律用WGS84经纬度不转投影交付给具体用户时再按用户的应用场景转对应的投影。因为不同用户需要的投影五花八门存原稿时如果转了某个投影后续换投影多一次重采样就多一次精度损失。分辨率的问题也类似。虽然250米的数据在很多省级尺度分析中够用但县级区域尤其是东部那些面积只有几百平方公里的县250米分辨率一个县也就一万来个像元统计均值的可靠性尚可但做空间格局分析就捉襟见肘。我后来针对县级产品加入了100米重采样版本代价是数据量膨胀但用户的反馈明显好了很多。2. 2000-2024年长时序NDVI的数据源选型2.1 主流数据源横向对比做2000年至今的NDVI数据绕不开几个主要数据源。MODIS系列产品的优势是时间序列长、波段设置适合植被指数计算、并且有现成的植被指数产品可以直接使用。它的时间分辨率是16天一年23期空间分辨率250米。这个产品我用的最多因为它的质量控制波段做得比较完善可以识别云、阴影、雪等污染像元。Landsat系列分辨率高得多 30米但受重访周期和云的影响单景有效观测非常不稳定。做全国范围、25年尺度的月度NDVI如果全部用Landsat数据数据的空洞问题和计算量都是灾难级别的适合做局部区域的精细分析不适合做全国批量生产。SPOT/VGT系列当年也是植被监测的主力但数据连续性在2014年之后接续情况比较麻烦现在已经不是主流选择了。国产卫星数据比如资源三号、高分系列近年发展很快但2010年以前的历史存档不足没法覆盖完整的2000-2024序列只能作为补充或者验证数据。我最终确定的主数据源是MOD13Q1250米16天合成的植被指数产品辅助数据源是MOD13A1500米产品用于大尺度宏观分析时的快速预览和校验。2010年之前和2015年之后的数据产品版本差异有一些但官方提供了相互兼容的再处理产品基本能够保证全序列的算法一致性。2.2 为什么不能直接拿单期影像做分析很多刚接触遥感的朋友会问既然MOD13Q1已经是16天合成了是不是直接下载下来裁剪就行理论上可以但实际做长时序分析时直接拿单期产品会碰到几个问题。一是日期不对齐。MODIS产品的16天合成周期是按儒略日划分的每年的第一期和最后一期的实际日期会有偏移。要生成月度或年度数据不能简单地把每个自然月内的影像直接叠加因为月与月的边界会把16天合成周期切开。二是云污染残留。虽然MOD13Q1已经做了初步的质量筛选但某些月份如果云覆盖太严重合成期的产品仍有异常低值。如果不做进一步的时间维平滑和异常检测统计出的月均NDVI会明显低于实际值。三是传感器退化问题。MODIS传感器服役多年辐射响应会有缓慢衰减虽然官方发布了校准产品但时序数据里仍然会残留一部分非植被因素造成的趋势成分。做长趋势分析时这个因素如果不处理很可能会得出植被逐年变差或者变好的错误结论。我的处理方式是16天产品先做严格的质量筛选只保留质量标记为良好quality band 等于0或1的像元对连续三期都被标记为云污染的像元在时间维上做插补。月度合成采用最大值合成法这本身就有一定的去云效果因为在一个月内至少有一期是晴空的概率远大于单期。年度合成同样用最大值法保证一年内每个像元都有机会取到植被最繁盛时期的信号。2.3 云平台处理与本地处理的取舍处理全国尺度、25年时序的数据计算量非常大。最初我在本地机器上跑单是下载原始数据就费了很多时间后面处理更是每天跑几百万个像元的运算。后来切换到某云平台的批量运算才顺畅起来。云平台内置了MOD13Q1的原始数据不用自己下载和预处理直接筛选、裁剪、合成导出结果到云端存储再批量下载到本地。这里省下来的不仅是存储空间更重要的是避免了大量重复的I/O操作。不过云平台也不是万能的。自己写脚本的好处是灵活出了问题能随时中断调试在云平台上调试批量任务有时候一个语法错误要等队列调度半天。我的建议是先在本地处理一个县、几个月的测试数据把流程跑通再上云跑全国的批量任务。测试时要注意流程的每一步都有日志输出这样批量跑的时候好定位问题。3. 省市县三级行政区划底图的处理逻辑3.1 行政边界数据从哪来、怎么挑全国省市县三级NDVI数据边界精度直接影响结果的可靠性。要是边界错了所有的省均值县均值都没有意义。我使用的是公开发布的基础地理数据省级和市级边界用了精度较高的全国地理数据库县级边界用了更新到2023年版的行政区划数据。为什么选这一版因为要避免一个很典型的坑有些县在2000到2024年之间发生过撤并、改名或者边界调整。比如某县在2015年撤县设区如果直接用最新的边界回算到2000年早期年份的数据就会混入附近其他县的地块。处理这个问题的办法是建立一套边界版本管理流程。对每个年份使用该年份对应的边界数据进行统计。如果边界调整幅度很小可以忽略如果涉及大范围撤并则提供两套统计口径——一套按现行边界统一裁剪一套按历史实际边界独立保留。使用时要根据分析目的去选择做历史趋势图用历史边界更准做现状监测用最新边界更方便对接土地报批数据。3.2 三级裁切的两种方案各有什么优劣行政区划的裁切有两种思路。第一种是先裁剪省级再在省级范围内按市级边界裁剪最后在市级范围内按县级裁剪。这种嵌套式裁切的好处是逻辑清晰每一级数据都严格包含于上一级范围内不会出现县级边界跨到省外的现象。缺点是文件数量特别多保存和命名管理起来麻烦。第二种是直接用全国范围的市级、县级边界一次性批量裁切。每级数据独立裁剪互不干扰。优点是省去了层级关联的复杂度缺点是如果市级和县级边界在交界处有细微不重合这种情况相当常见比如河流边界在国家级、省级、县级数据里画的位置不完全一致同一块地在市级统计和县级统计里可能出现非常小的面积出入。两种方案我最终结合使用了——市级和县级数据分开独立成库不强行追求每一级数据之间的面积完全吻合而是用面积权重修正的方式把市级汇总值尽量调整到与省级汇总值接近一致。这样做的好处是每一级数据都可以单独交付给用户不用携带父级信息。因为很多业务单位只要某几个县不需要先下载整个省再裁。3.3 文件命名规范和目录结构这套数据的文件命名是一个细节活但很容易被忽略。我最终采用的命名规则是NDVI_省份代码_市代码_县代码_时间粒度_年份_月份_分辨率.tif。举例NDVI_CN_110000_110100_110101_月度_2020_07_250m.tif。文件名里必须有标准的行政区划代码因为后续很多自动化流程要按代码匹配属性表用中文名会有各种莫名其妙的编码问题。目录结构上我按省/市/县三级建目录每个目录下按年归档。年度合成文件单独建一个文件夹月度文件按月份顺序排。设计这个目录结构的时候有一个小的后悔点早期我没有把分辨率写进文件名后来同一地区出现250米和100米两个版本文件名完全一样系统里同时跑两个流程直接互相覆盖了。吃了几次亏之后我强行要求所有文件名必须带分辨率和合成方式月度/年度并且把分辨率、合成方式、坐标系、数据源版本四项信息全部写进文件的元数据属性里确保用户拿到任何一个文件哪怕不看文档也能从文件名读全信息。4. 批量化生产流程的具体实现4.1 整体流程导览整个生产链路是原始数据读取 → 质量筛选 → 时间合成 → 全国拼接 → 按边界裁切 → 异常值清洗 → 重投影 → 单位换算 → 命名归档。听起来不复杂但每个环节展开都有讲究。下面的表格是我实际生产中用到的关键参数处理环节关键参数说明数据源MOD13Q116天合成、250米质量筛选QA波段保留0、1两级剔除2、3级月度合成最大值合成法一个月内所有有效像元取最大年度合成最大值合成法12个月合成后再取最大坐标系WGS84存档统一用经纬度分辨率省级250米 / 市级250米 / 县级100米按需重采样数值范围[-1, 1]存储时转为整型放大10000倍空值处理双线性插值补齐只做时间维不做空间维4.2 最大值合成法与质量控制参数详解在众多NDVI合成算法里最大值合成法一直是生产植被指数产品的主流选择。原理其实很简单在一段时间内取每个像元观测到的最大NDVI值作为该时段代表值。为什么取最大因为云、气溶胶、太阳高度角变化这些因素都会导致NDVI被压低很少会让NDVI虚假偏高。所以在一段时间内取最大值相当于尽量保留了植被生长的真实信号。极端情况下如果整个时间段内全是云覆盖最大值法也会得到一个明显偏低的异常值这就是后面还要做质量控制的原因。我采用的严格质量控制流程是分层的先通过MOD13Q1自带的像元可靠性波段筛掉质量标记为云污染冰雪覆盖大气校正失败的像元如果实在没有合格的像元再用经过校正的相邻观测补位还不行就标记为待插补。影响月度合成结果的一个隐藏问题是生长季和非生长季的NDVI差异巨大。1月的最大值可能只有0.2左右而7月的最大值可能达到0.8。如果做全年12个月的等权平均冬季的低值会把植被全年长势水平拉低不能真实反映植被覆盖上限。所以做年度统计时我固定输出两个版本全年最大值版和生长季均值版生长季定义在多数省份是4月到10月。用户在评估植被覆盖能力和生态系统生产力时用生长季均值版更合理。4.3 裁剪和重采样过程中的经验参数裁剪这一步很多人觉得简单用一个矢量边界去切栅格而已。但全国范围的裁剪有一个参数问题很多人直接用默认设置结果输出栅格的像元大小变了、影像行数也变了每个县的数据像元位置对不齐。我在裁剪时候固定了一个操作规范关闭栅格自动匹配选项强制设定输出分辨率和像元原点保证所有县级文件像元网格完全对齐。这样跨年份对比时同一个像元位置对应的是同一块地而不是因为原点漂移出现半个像元的位置错动。这个问题在做长时序差值分析时尤其致命差值图会出现明显的条带噪声。重采样也有讲究。从250米重采样到100米我选择立方卷积法它比双线性更锐利边缘保留效果更好适合县级小范围。从250米重采样到500米这种降采样的场景则用像元聚合取平均避免信息冗余。绝对不要用最邻近法做NDVI重采样因为NDVI在空间上变化是连续的最邻近法会把边界台阶式放大后期做平滑时麻烦很多。4.4 整型存储与浮点转换的小心思原始NDVI是浮点数范围在-1到1之间。直接存成浮点型TIF每一个像元占32位单精度浮点数据量非常大。全国范围逐月数据几十万个文件直接用浮点存储存储空间根本吃不消。我做了一个换算处理将NDVI加1再乘以5000取整后存储为Int16整型这样值域是0到10000。在分析的时候再把整型转回浮点公式是NDVI 整数值 / 5000 - 1。这个放大倍数的选择也是有讲究的——5000倍的放大在绝大多数应用场景下精度足够而且Int16范围-32768到32767刚好能包容这个取值范围。选10000倍虽然精度更高但细节上没有质的提升存储消耗和文件读取速度却明显变差。做这一步的时候一定要把换算公式写进数据的说明文档里。我见过不止一次下游用户拿到整型数据后不知道要除以5000直接拿0到10000的数值去做分析结果所有植被覆盖度计算全部爆表。因为NDVI计算和覆盖度的公式前提条件就是NDVI值域在-1到1之间用放大的数去算整个结果是错的。5. 数据精度验证和质量控制记录5.1 数值范围与空间分布合理性检查拿到一批裁剪好的栅格我第一件事不是看均值而是先做值域检查——统计每个文件的NDVI最小值、最大值、平均值和标准差。如果某个县某个月的NDVI最小值明显低于-0.2或者最大值达到1以上基本可以断定存在异常像元。时空分布合理性的检查也必不可少。比如冬季的东北、西北地区NDVI普遍应该在0.1左右部分雪覆盖区域可能更低。如果某个县1月份的NDVI均值达到0.6那不用想肯定是数据出了问题。这种检查虽然简单粗暴但往往能最快发现问题。我一般用两个指标来判断一个文件是否异常一是全局平均值的Z分数是否大于3意思是这个县这个月的数据和历年来同月份的数据差异太大二是文件内NDVI大于0.8的像元占比是否超过5%。如果某个月份出现大范围异常高值多半是传感器定标或者大气校正的问题需要回溯原始数据。5.2 与其他数据集的交叉验证在重点区域我会拿这套数据和更高分辨率的数据做交叉验证。做法是选一个县用Landsat 30米分辨率数据计算同期的NDVI然后聚合到250米再与这套数据对比。两套数据之间的相关系数一般能达到0.8以上这表明数据在空间格局上是可靠的。时间序列的验证则是找一个长期定位观测站点的实测植被覆盖度数据把对应像元的NDVI序列提取出来去看年际变化趋势是否和实测数据一致。例如某退耕还林重点县的NDVI应该呈现持续上升趋势且斜率在统计上显著。趋势方向一致说明数据能够反映真实的植被变化。5.3 时间序列连续性的检查长时序NDVI数据集最容易出现的问题是某些年份数据系统性偏低或者偏高。因为传感器老化、定标参数更新、大气校正算法的改进都会造成产品之间的细微差异这些差异在单景影像上很难察觉但放在25年的长时间序列里就会被放大。我用两个办法做时序一致性修正。第一双平台交叉校准在同一时间段用MODIS数据和Landsat数据进行对比计算两者的偏差函数再把偏差从MODIS时序中扣除。第二谐波分析法对每个像元的时间序列做傅里叶分解把年周期信号分离出来然后检查残差项是否存在系统性的年际偏移。如果发现某几年的残差整体偏大我会去查该年份的数据版本说明再用邻近年份做线性调整。这个校准过程不需要对所有像元做——那样计算量太大了。我按生态分区选了大约3000个样本像元覆盖森林、草地、农田、荒漠、湿地等主要生态系统类型用样本像元算出一个全国尺度的校准系数再应用到全体像元上。这样做的好处是高效坏处是生态分区交界处可能出现轻微的不连续。实际生产中这个不连续通常不会影响省级和市级的统计但在县域尺度的空间分析中还是要注意必要时可以对校准系数做空间平滑。6. 常见问题速查表与我的避坑经验现象可能原因解决办法某个县某月NDVI整体异常低云污染未完全清除检查QA波段、改用相邻年份同期数据插补边界与行政区划不匹配坐标参考不一致统一用WGS84重投影后再裁剪文件名乱码或冲突命名规则不统一强制使用行政区划代码命名重采样后边界模糊重采样方法不对升级为立方卷积法时间序列出现阶梯跳变数据源版本混用统一使用同版本产品2000年前后数据差异大传感器型号变化做交叉校准调整统一系数存储空间过大浮点型存储转Int16整型存储县级面积太小、统计不稳定分辨率不足对县级单独提升至100米分辨率6.1 一个我到现在还在用的检查命令如果某个文件出了问题我最先做的就是提取该文件的基本统计信息并和相邻年份的同月份文件做差。原则上每年同一月份的数据差值应保持在合理范围内比如生长季差值一般在正负0.1以内非生长季也不应超过正负0.2。如果差值图出现了一个县内部的明显分界线十有八九是数据拼接时的条带问题而不是真实的植被变化。6.2 关于加值、改名、跨界的历史遗留问题中国的行政区划在过去25年里有过不少调整这是做长时序省市县数据时绝对绕不开的坎。最典型的就是县改区。某县如果变成了某市的一个区那县的名字在后续年份就消失了早期年份的县数据和后面年份的区数据没法直接同名衔接。我的处理方式是在属性表里保留一个现状归属字段不管历年怎么改名字都用最新的区划代码作为主键把历史数据归并到最新边界下。还有一种情况是乡镇级飞地。某些县在市域范围内有飞地虽然面积不大但如果在裁剪时没有把飞地正确归入所属县会导致该县的NDVI均值出现明显的偏移。处理方法是使用带有乡镇级属性的三调数据构建县界这样能最大限度保留飞地的归属关系。6.3 批量处理中如何防止任务静默失败全国批量裁剪过程中最怕的不是报错暂停而是任务静默失败——某个县的某个月文件因为没有有效数据输出了一个全零或者全空的文件但流程却认为是成功的。我的做法是每一个输出文件都立即进行有效性校验并记录日志比如检查行数和列数、像元总数、有效像元比例、值域范围、平均值是否在合理区间。任何一个校验项不过就把该文件标记为异常并且在当轮处理结束时汇总输出一份异常清单。有了这份清单我只需要处理异常文件就行而不用把这个县的数据全部重新跑一遍。这里有一个很小的技巧输出文件名带任务批次号。即使同一个县的数据反复处理了多次也能区分哪一次处理的结果是最终版本避免文件被后来的错误处理覆盖掉。6.4 最终交付时除了数据还应该附什么数据文件本身很重要但交付时如果只给一堆栅格文件用户拿到后大概率会回来问你各种问题。我的交付包一般包含这几个部分文件说明文档README写明数据源、合成方式、坐标系统、分辨率、值域换算公式、文件命名规则。行政区划版本说明写明用了哪一年的区划底图哪些区域在历史上发生过调整不同年份的口径差异。质量报告包括每个文件的统计结果、异常清单、插补列表、精度验证汇总。样例脚本包括读取数据、转换NDVI、计算区域均值等常用操作的示例代码。这些附加材料会增加一倍的工作量但实际使用时的效果会好很多倍。很多用户拿到一个栅格文件之后并不知道里面的整型数值是要换算的也不知道最大值合成和均值合成的区别。有一份说明文档在能避免很多低级误解。7. 最后想补一句实在话全国各地被反复索取NDVI数据的场景其实比我们想象的多得多从生态修复评估到农情长势监测从历史旱情重建到城市热环境分析大家都在用同一个底子上的植被指数做文章。这套25年跨度、省市县三级的数据集做的时候花了大量精力在行政边界怎么配、时间序列怎么对齐、异常值怎么清洗这些看似不起眼的地方。但恰恰是这些隐蔽的细节决定了数据最终分析结果的可信度。我个人的体会是遥感数据的生产八成的精力其实不是花在遥感本身而是花在处理好行政区划、坐标系、时间轴、文件命名这些枯燥琐碎的事情上。这些环节做得越扎实后面分析阶段踩雷的概率就越低。这套数据后续还可以继续往两个方向扩展一是把时间粒度细化到旬甚至周满足物候研究的更高需求二是叠加气象数据做植被与气候因子的联合分析。如果哪个方向有需要我可以接着把做法整理出来分享给大家。