MODIS 这摊子事玩遥感的几乎没有不知道MOD44B的。这个产品全名叫 Vegetation Continuous Fields中文圈一般叫“植被连续场”或者干脆叫“植被覆盖百分比”。我这次做的是把 2000 到 2020 年中国区域的 MOD44B 数据整理成一套干净可用的植被覆盖百分比数据集周期长、瓦片多、坑也不少。整个过程走下来从下载、拼接、投影处理到时间序列分析每一步都有最容易被忽略的细节。这篇文章我把整套操作思路、遇到的问题和现在的处理方案都写出来想碰这套数据的人应该能少走不少弯路。先说清楚这东西适合谁看。如果你是要做长时序植被变化、荒漠化监测、生态工程评估或者只是想给论文补一张像样的植被覆盖趋势图那么 MOD44B 会是一个比 NDVI 阈值法更省心也更靠谱的选择。当然如果你从来没接触过 HDF 文件、MODIS 正弦投影和瓦片拼接这篇文章也足够带你走完一个完整流程不需要额外翻太多文档。1. 为什么拿 MOD44B 当主角给中国植被覆盖长时段监测做底稿1.1 一个植被覆盖研究中绕不开的产品用遥感做植被覆盖有一条经典路线是拿 NDVI 算覆盖度公式现成、工具也多。但这条路有个让人头疼的地方需要自己定阈值、定像元二分模型参数不同区域不同传感器结果还不一样。MOD44B 的逻辑完全不同它用训练好的回归树模型直接输出“每个像元里树木覆盖、非树木植被、裸地分别占百分之多少”把覆盖度这个本来要靠估算的变量变成了一种标准化的观测值。这个产品来自 NASA 的 MODIS 传感器搭载在 Terra 卫星上。卫星从 2000 年起持续回传数据所以它天然就是做 2000 年至今长时序分析的骨干数据源。空间分辨率 250 米一个像元相当于地面 250×250 米的范围对省域、国家级别的生态环境监测来说很合适既不会粗糙得看不了细节也不像 30 米 Landsat 那样需要海量存储和计算资源。MOD44B 最核心的一点是它把“植被覆盖”这个问题分段处理不再笼统地说“这里植被多那里植被少”而是把覆盖结构拆成树覆盖、低矮植被和裸地。这种分层在西南林区、西北荒漠区、华北农田区之间做对比时特别有用比如同样是 60% 的覆盖度纯草地和混交林在结构和生态功能上完全不同但如果只看总数就很容易被误导。1.2 MOD44B 对中国区域意味着什么中国的地形和植被类型跨度大从热带雨林到温带荒漠都有分布这对任何一个全球尺度覆盖产品都是一个考验。MOD44B 在很多地区的表现比较均衡尤其在大尺度森林监测中已经被广泛验证。它对中国区域的覆盖是完整的只要把覆盖中国的 MODIS 瓦片下载齐全再经过拼接裁剪就能得到一张完整的 250 米分辨率的全国图。时间范围方面2000 到 2020 正好跨越了 21 年足够支撑趋势分析和突变检测。这些年里中国经历了退耕还林、天然林保护等一系列重大生态工程也经历了部分区域的城市扩张和荒漠化。用 MOD44B 算出的树覆盖趋势图可以直接看到哪些地方在变绿、哪些地方在退化不用再靠定性描述或者小区域样地推断。我从实际处理中得到的一个直接经验是MOD44B 非常适合当做“底图”或者“主变量”来用。你可以把它和气象数据、土壤数据、夜间灯光数据叠加做驱动因子分析也可以把它简化成“植被覆盖度大于多少”的图来制作宏观生态区划。它的语义清晰、可解释性强基本上任何后续分析都能接得住。2. MOD44B 数据集结构与背后算法深度拆解2.1 四个子数据集就是四张“体检表”打开一个 MOD44B 的 .hdf 文件里面不是一张单波段影像而是一堆子数据集。一般你能看到四个主要的层树木覆盖百分比Percent_Tree_Cover、非树木植被覆盖百分比Percent_NonTree_Vegetation、裸地/非植被百分比Percent_NonVegetated以及一个质量控制层Quality。数值范围是 0 到 100代表百分比。要注意的是这三个覆盖百分比加在一起并不一定严格等于 100因为算法本身存在不确定性和噪声加上云污染、传感器视角等影响像元内会有残差。处理的时候不要强行对三个波段做归一化加和否则会人为引入误差。质量控制层是很多人下载完就忽略的但它恰恰是最有价值的一个波段。质控层用位标志Bit Flag记录了像元质量、水体、云、雪等情况。比如某些像元是水体、某些像元被云覆盖如果不加过滤就直接参与统计最后出来的趋势图会被严重干扰。这一点后面分析处理时我会专门展开讲。算法层面MOD44B 使用的是回归树模型每个树节点会根据光谱反射率、纹理信息等特征把训练样本分成不同类别最终输出的是连续覆盖度而不是离散类别。这种做法的好处是抗噪能力强训练样本中混合像元的覆盖度能被更平滑地估计出来比单纯的最大似然分类要稳定得多。2.2 投影、瓦片与文件命名的几个硬知识点MOD44B 默认采用 MODIS 正弦投影Sinusoidal并且按瓦片tile组织全球数据。整个全球网格横向从 0 到 35纵向从 0 到 17一共 36×18 个瓦片每个瓦片约 10 度经度、约 8 度纬度。中国区域大概覆盖了 h23 到 h28、v03 到 v07 这个范围的瓦片具体到不同研究区还需要再确认。文件名也有一套固定的逻辑比如MOD44B.A2000001.h25v05.006.2015167203617.hdf。拆开看就是产品名 A 年份加儒略日 瓦片编号 数据版本 生产时间。这里有个容易被忽视的点虽然文件名里写的是 A2000001表示 2000 年第 1 天但 MOD44B 实际是年度合成产品不是每天的影像你展开来看会发现这个“1 月 1 日”其实是代表一整年的结果。所以做年度趋势分析的时候直接按年份来组织文件就行不要纠结具体日期。下载前建议先用 MODIS 瓦片分布图或者工具算出自己研究区涉及的瓦片编号。别凭感觉猜我之前就吃过亏一开始只下了 h25v05 一个瓦片结果图一拼出来发现缺了一大块。工具方面可以直接用 MODIS 官方的 Tile Calculator或者打开一个全球瓦片索引矢量文件把你的研究区 Shapefile 叠上去一目了然。2.3 与其他植被覆盖产品对比选它的理由市面上能拿来做植被覆盖的数据不少我试过几类简单说下体验基于 NDVI 的覆盖度提取灵活但需要自己标定像元二分模型阈值换一个区域就得重新调不好做横向比较也不适合大范围自动处理。MODIS 土地覆盖产品 MCD12Q1分成离散类型比如常绿针叶林、落叶阔叶林、草地但没有连续的覆盖比率生态学计算中往往需要转成连续量转换过程容易损失精度。GLASS 植被覆盖度产品是另一个不错的连续覆盖产品时间序列也长但下载路径和多数据源融合流程对初学者不够友好。Landsat 系列空间分辨率高但 30 米数据覆盖全国需要处理太多景影像2000 年以前的还好说做 20 年长序列的存储和算力需求比较大。MOD44B 的 250 米分辨率、连续百分比输出、单一来源卫星持续观测让它在大尺度生态监测里成了一个“性价比”很高的选择。尤其适合做国家级或省域级别的时空变化分析以及配合气候变量做驱动分析。如果你是做小流域精确制图的那 MOD44B 的分辨率不够还是回到 Landsat 更实用但如果你需要快速拿到一套中国的长时序植被底图它就是当前最省事的方案之一。3. 2000-2020 中国区 MOD44B 数据获取全流程3.1 前置准备账号、工具和网络节奏去 NASA 的 Earthdata 官网下载 MODIS 数据第一步是注册一个 Earthdata 账号。这个账号是免费申请的但别想着绕过登录直接抓链接现在官方数据分发基本都要求带登录凭证。注册时填好邮箱和机构信息即可选学科方向时可以选“Ecology”或者“Earth Science”审批一般很快我那次是几分钟之内就收到了确认邮件。工具方面建议备齐这几个一个支持 HDF4 读取的 GDAL 环境一个支持断点续传的下载器或脚本再加一个影像拼接与裁剪工具。GDAL 可以通过 Anaconda 安装也可以用 OSGeo4W 包。国内玩家如果在 Windows 环境下用 pip 安装 GDAL 遇到兼容问题建议直接装gisinternals或osgeo的预编译版本省去自己编译 HDF4 驱动的痛苦。下载节奏也是经验之谈。NASA 服务器在国内的响应速度不稳定高峰期经常几十 KB 每秒一个瓦片传半天。我后来基本选择凌晨或上午下载成功率明显更高也建议用支持断点续传的下载脚本别断了就从头再来。3.2 官方下载路径的选择与细节MOD44B 的官方下载路径主要有两条一条是 NASA Earthdata Searchsearch.earthdata.nasa.gov另一条是 AppEEARSappeears.earthdata.nasa.gov。Earthdata Search 适合下瓦片原文件。搜索框里输入 MOD44B然后在时间面板选择 2000-2020再在地图上画出中国大致的范围系统会把涉及的所有瓦片列出来。你可以添加到购物车后批量下载。缺点是它给的是原始 HDF 瓦片后面还需要自己做投影转换和拼接。AppEEARS 更适合做点采样和区域提取。你上传一个 Shapefile 或者直接在图上画范围选择要的产品和年份系统会自动裁剪、重投影可以指定 WGS84输出 GeoTIFF。这种方式的用户体验好但局限是如果研究区特别大、时间序列特别长提交任务后等待时间会很长而且返回的数据组织方式有时候需要再整理。真正追求批量自动化的人我建议直接用脚本下载。Earthdata 提供了wget或curl配合 token 的方式也有现成的 Python 第三方库比如pyModis和modis-tools。我个人现在的习惯是写一个小脚本自动读取瓦片列表、逐条下载并校验文件大小如果下载失败就自动重试三次这比在网页上手工点上万次鼠标要靠谱得多。3.3 批量下载规划瓦片清单与命名规则中国区域要下载全 21 年如果使用瓦片方式一年大概需要 15 到 25 个瓦片全部下来文件总数可能会接近四五百个。单个 HDF 文件几十到上百 MB整套数据完整落盘可能要几十 GB。下载前最好按照“年份/瓦片”的目录结构存放后面处理时检索会快得多。我整理瓦片清单的办法是先准备一个全国 Shapefile然后转成 MODIS 正弦投影和全球瓦片网格叠加把相交的瓦片编号提取出来导出成tiles.txt。再用一个循环逐一拼 URL 下载。这种办法虽然第一次准备有点麻烦但以后你再下别的 MODIS 产品可以直接复用同一份瓦片清单一步到位。下载完一定要抽查文件是否完整。HDF 文件损坏时GDAL 打开会报错但有些文件表面能打开、内部却缺层。稳妥的做法是写一个检查脚本遍历所有.hdf尝试用 GDAL 打开并统计子数据集数量小于预期数量的文件直接标记为异常重新下载。别小看这一步我在东北地区的瓦片中真的碰到过几个文件下载中途断了但大小看着正常结果一打开全是空数据。4. 数据处理与质量控制的核心操作4.1 数据可以这样处理拼接、裁剪、重投影拿到一堆瓦片以后处理链路一般固定为“拼接瓦片→重投影→裁剪到研究区→按年份合成时间序列”。拼接这一步用 GDAL 的gdalbuildvrt很方便它不会真的拷贝像素而是生成一个虚拟栅格把多个瓦片链接在一起。等做完投影转换后再输出成实际文件可以省掉一次中间文件的转换时间。举个例子把 2001 年的多个瓦片拼接并转换到 WGS84可以用这几条命令gdalbuildvrt mod44b_2001.vrt MOD44B.A2001001.*.hdf gdalwarp -t_srs EPSG:4326 -te 73 18 135 54 -tr 0.0025 0.0025 mod44b_2001.vrt mod44b_2001_wgs84.tif这里-te指定的是范围中国大致经纬范围-tr指定输出像元大小。0.0025 度大概对应 250 米左右和原始分辨率保持一致。裁剪到行政区时可以直接用-cutline参数加一个中国省级边界矢量gdalwarp -cutline china.shp -crop_to_cutline mod44b_2001_wgs84.tif mod44b_2001_china.tif需要注意gdalwarp在处理多波段 HDF 子数据集时要先拆出你需要的那个覆盖度层。HDF4 里的子数据集不能像普通 GeoTIFF 一样直接整个文件扔给 GDAL得用HDF4_EOS:...子数据集路径或者用 Python 先把子数据集转成 GeoTIFF 再继续处理。这也是新手最容易卡住的一步。4.2 质量图层和数值缩放怎么用才不掉入陷阱MOD44B 的覆盖度层虽然是 0-100 的百分比但在原始 HDF 中存储类型通常是整型可能还带缩放因子。有些处理软件会自动应用缩放有些不会这就要你看数据说明文档确认scale_factor。如果确实有缩放使用前要乘以系数别把原值直接当百分比否则时间序列会整体偏低或偏高。质量层的使用更要谨慎。它不是简单的“0好1坏”而是位标志编码不同 bit 位代表不同含义有的是像元质量等级有的是水体掩膜有的是云/雪覆盖有的是气溶胶程度。想过滤像元需要把质控值转换成二进制按位读取对应标志位。我见过不少人在 GEE 里直接写quality 0这样很可能把大量边界区域或山区像元误删。我的做法是把质量层读出来然后根据产品说明把需要屏蔽的像元集合云、水、雪、低质量统一设为 NoData。对 MOD44B 来说一般会去掉水体和云污染的像元保留陆地和非云区。这样处理后年际间像元覆盖更加一致趋势分析时不会因为某一年的云污染而出现局部异常。为了不丢数据处理时建议用一个栅格数组保存数据一个布尔数组保存掩膜。每一年先把未通过质控的像元填充为 NaN再参与后续统计。否则直接把 NoData 记成 0会让裸地区域的比重被系统性抬高整张趋势图直接失真。4.3 时间序列统计和趋势计算要留意的细节21 年序列的统计工作最经典的是算每像元的趋势斜率。常见方法是最小二乘线性回归更稳健的选择是 Theil-Sen 估计因为它对异常值不那么敏感。下面这段 Python 思路可以复用from osgeo import gdal import numpy as np from scipy import stats slope np.ones((rows, cols)) * np.nan for i in range(rows): for j in range(cols): vals series[:, i, j] mask np.isfinite(vals) (vals 0) (vals 100) if mask.sum() 12: res stats.theilslopes(vals[mask], years[mask]) slope[i, j] res[0]在写循环之前记得先把所有年份数据读入一个三维数组维度是年份行列并且统一所有年份的 NoData 表达。还有一个容易忽略的问题是不同年份的边界范围可能因为裁剪和投影精度略有偏差必须把所有年份重采样到完全相同的网格相同的范围、行列数、分辨率否则数组逐像元比较时会错位。做显著性检验时我建议同时输出 p 值。p 值小于 0.05 的像元才标记为显著变化否则很多弱趋势会被图上颜色放大成“大片变绿”误导读者。实际写论文、画图时把“显著变绿”“显著退化”“无显著变化”三类的面积统计出来比只扔一张斜率图更有说服力。5. 被问爆的相关 MODIS 产品气象信息系列与 NPP5.1 MODIS 里能拿到的“气象相关”数据集有段时间很多人在后台问我MODIS 到底有哪些产品能获取气象相关的信息这里统一说清楚。严格意义上MODIS 不是气象卫星它是遥感对地观测平台很多产品更偏“地表状态”而不是“气象观测”。但以下产品和气象研究密切相关的变量直接挂钩MOD07_L2大气剖面产品包含温度廓线、湿度廓线、大气可降水量等是做大气校正和天气过程分析的常用辅助数据。MOD05_L2近红外可降水量产品能反演大气中的水汽总量对研究降水过程和水循环也很有用。MOD06_L2云产品包含云顶温度、云顶气压、云相态、云光学厚度等是天气分析里常用的参数。MOD11A1 / MYD11A1地表温度产品LST虽然叫地表温度但和近地表气温关系密切很多山区气温插值研究都会用 LST 做协变量。MOD13A1植被指数产品NDVI/EVI虽然不是直接的气象变量但做干旱监测时经常和降水、温度数据联合使用。这些数据同样可以通过 Earthdata 和 AppEEARS 获取处理链路和 MOD44B 相似都是 HDF 瓦片、正弦投影需要重投影和裁剪。如果你做植被-气候相互作用分析一般会拿 MOD44B 做植被覆盖变量再配 MOD11A1 的地表温度和 MOD07 的水汽变量作为气候侧的解释变量效果很好。5.2 MOD17A3H NPP 下载实操记录MOD17A3H 是年度的 NPP净初级生产力产品空间分辨率 500 米版本现在以 MOD17A3HGF 为主。很多做碳循环、植被生产力评估的朋友都来问“MOD17 NPP 怎么下载”其实流程和 MOD44B 一模一样只是产品代码不同。走 AppEEARS 的步骤很简单登录 AppEEARS点新建任务在数据列表里搜MOD17A3HGF选年份比如 2000-2020上传你的研究区 Shapefile或者在地图上框选范围指定输出格式为 GeoTIFF、投影为 WGS84提交任务等邮件通知下载。要注意的是MOD17A3HGF 子数据集不止一个有 GPP、NPP、NPP_QC 等需要按需勾选。NPP 的单位是 kg C/m²/year下载后要先确认数值范围遇到特别大的值需要查看是否包含填充值。产品说明里通常给了scale_factor不按缩放直接出图的话数值会有量级偏差。另外MOD17 产品经常更新有不同版本老版本的数据可能在服务里已经下线。如果你用的是 Earthdata 直接搜 HDF 瓦片注意把版本号对清楚.061和.006在文件名里都会有体现。5.3 一份能少走弯路的 MODIS 数据获取清单很多新手一上来就埋头下载结果后期分析时发现缺少气候数据、缺少辅助验证数据又得重新补。这里分享一份我常用的清单核心植被覆盖MOD44B树覆盖、非树植被、裸地植被生产力MOD17A3HGFNPP/GPP植被指数MOD13A1 或 MOD13Q1NDVI/EVI地表温度MOD11A1逐日或 MOD11A28 天合成大气水汽MOD05_L2土壤湿度SMAP 产品不是严格 MODIS但与 MODIS 常配合使用夜间灯光NPP-VIIRS做城市扩张和人类活动影响评估时用这套组合基本覆盖了“自然人为”两类驱动因子的主流变量。做 2000-2020 的中国植被覆盖研究我建议先把这个清单里的变量下载到本地后续不管做回归分析还是机器学习驱动因子识别都不用再为数据来源发愁。6. 实操过程问题排查与避坑记录6.1 下载和处理中途最容易暴雷的地方把实际操作里遇到最多的几个问题整理成一张速查表遇到相似问题可以直接对症下药现象原因解决思路HDF 文件打不开GDAL 没有编译 HDF4 驱动安装完整版 GDAL或改用 gisinternals 包下载到一半就断NASA 服务器连接不稳定使用断点续传脚本设置失败重试机制数据全是 0 或 255没做 NoData 和缩放处理检查产品文档里的 _FillValue 和 scale_factor拼出来的图有黑缝瓦片边界有重叠或填充值没统一在 gdalbuildvrt 里设置-srcnodata拼接后统一处理不同年份边界对不上重投影范围不一致所有年份在同一个-te和-tr下重采样植被覆盖百分比加起来超过 100算法本身存在残差不要强制归一化保留原始输出做分析趋势图上出现大片异常没有用质量层过滤云和水按位解析 Quality 层建立统一的掩膜这些问题的核心共同点是你当它是普通图片但它是科学观测数据。每个产品都带着自己的坐标系、缩放系数、填充值和质量标志拿到新数据的第一件事绝对不是开心地画图而是先读产品文档把元数据搞明白。6.2 一个典型的全年数据提取流程这里给一个可以直接参照的全年数据提取流程它综合了下载、处理和质量控制的关键步骤。第一步确定研究区和瓦片号。我一般直接用全国范围的 Shapefile叠加 MODIS 瓦片网格自动提取出要下载的瓦片列表。第二步写下载脚本从 Earthdata 拉取 HDF每个文件下载后自动校验大小小于阈值的重新下载。第三步将每个 HDF 的子数据集拆出来树覆盖、非树植被、裸地、质量层分别存成 GeoTIFF。第四步逐年用 gdalbuildvrt 拼接再统一重投影到 WGS84 并裁剪到研究区。第五步根据质量层生成掩膜把云、水体、低质量像元设置为 NaN。第六步把所有年份堆叠成一个三维数组做趋势分析、统计面积和画图。这六步听起来挺线性实际执行中会有很多循环和调试。比如拆子数据集那一步HDF4 的路径经常因为文件名里版本号变化而改变需要写正则表达式去匹配又比如掩膜处理时不同年份有效像元数差异较大后续趋势计算要把有效像元数少于一半的年份像元剔除掉否则回归结果不稳定。6.3 我反复踩过、也反复修改的几条心得第一不要迷信“一个产品走天下”。MOD44B 确实好用但它是年度的只反映一年内覆盖状态的综合不能看出植物猛然变绿或者叶片凋落的季节变化。如果要做物候分析需要补充 MOD13Q1 等 16 天合成的植被指数产品。第二下载数据时要带着“验证需求”去下。如果你只是为了构图下载三五年就够用了但做趋势分析建议把 2000-2020 的完整序列都下齐因为后续要补某一年的数据是非常痛苦的NASA 的数据档案虽然保留完整但下载服务器有速率限制临时补一两年数据的等待时间比一次下载完整序列更闹心。第三质控图一定要留档。很多论文要求提供像元有效数量和掩膜情况的说明如果处理时没留质控图后面补绘会非常麻烦。我在处理时会给每一年都出一张“有效像元占比”图做成一个逐年变化小动画既能检查数据质量又能给报告添彩。第四处理中间文件尽量选择 GeoTIFF 而不是继续用 HDF。HDF 文件虽然保留了最完整的信息但读取速度慢而且很多软件对它的支持不完全。拆成 GeoTIFF 后以后用 QGIS、ArcGIS、Python 都能顺畅操作。当然原始 HDF 还是要备份一份防止处理过程出错后无法回头。最后再分享一个小技巧做别把全中国放进一幅图里出趋势结果。这个大区域里气候差异太大同一套趋势参数在东南湿润区和西北干旱区解释力完全不同。我实际操作时是按八大生态区或者省域分区分别统计的效果比全国一刀切要清楚得多。这个处理方式在写结论时也更有说服力。