中国区域MOD15A2H FPAR多年平均数据集：从QA解码到平均值合成的完整实践

最早接到这个需求时，我心里其实有点犯嘀咕：MOD15A2H官方产品都上线这么多年了，下载流程也成熟，为什么还要专门做一份中国区域2000—2025年的FPAR年均值合成数据集？等真正把几十个正弦投影瓦片、二十六年下来上千个HDF文件全部拉齐，走完预处理、合成、质控这一整套流程之后，我才发现这个问题远没有想象中简单。这篇记录我不会只讲流程，还会把选型理由、QA解码、平均值合成算法设计、多年平均口径这些容易被忽略的细节一起说清楚，给做生态遥感、农业遥感、碳循环建模和GIS底图的朋友一份可以直接参考的实操经验。

1. 为什么非要自己做一份"多年平均FPAR"，而不是直接用官方产品

1.1 FPAR是什么，凭什么它比NDVI更值得关注

光合有效辐射（PAR）是指太阳辐射里400—700纳米波段的光，绿色植物的光合作用主要利用这个波段。FPAR就是植被冠层真正吸收的PAR占到达冠层PAR的比例，取值通常在0到1之间。

传统指数里NDVI用得最多，但它本质上是"绿度"的表观指标，对土壤背景、观测角度、大气状况都很敏感，在茂密植被区还容易饱和。FPAR则直接跟光合能量转化挂钩，很多生产力模型把FPAR作为核心驱动输入，所以它在碳循环、作物长势、生态评估里的地位非常特殊。用生活里的类比来说：NDVI告诉你这棵树看起来绿不绿，FPAR告诉你这棵树真正吃到了多少阳光。

做长时序生态遥感的人，拿到一份可靠的多年来平均FPAR，很多问题就能直接回答了：这个区域的植被平均光合能力到底处在什么水平？不同区域之间的平均生产力差异有多大？模型参数初始化应该给一个什么先验值？

1.2 官方产品直接拿来用，会碰到三个现实问题

第一，产品是瓦片式的。MOD15A2H按正弦投影分成若干瓦片，中国陆地差不多覆盖十几幅瓦片。用户如果只下载自己需要的那几幅，不做统一重投影和拼接，国界边缘、瓦片边缘都会出现明显错位和条带。

第二，官方产品是8天合成的时间序列，不是成熟稳定的气候态产品。单期数据受云、气溶胶、传感器角度影响非常大，尤其在中国这种季风气候明显、多云雨的区域，直接用单期数据做区域平均，误差大到没法用。哪怕做年均值，如果不先做严格的质量筛选，一条污染很大的期次就能把整年结果带偏。

第三，官方产品并没有直接提供"多年平均"这种气候态结果，而气候态恰好是区域对比和模型初始化最常用的基准。如果每个用户都从HDF原始文件开始自己处理，时间和算力成本非常高，而且各人选的投影参数、QA阈值、合成算法都不一样，结果之间没法横向比较。做一份统一口径、统一质控流程的多年平均数据集，其实就是替多数用户把最枯燥、最容易出错的环节先做掉，并且留下可追溯的处理记录。

1.3 谁最需要这种平均值合成数据集

从我这边的使用反馈看，需求主要来自三类人。

第一类是生态和农业方向的研究人员，他们想评估不同区域植被生产力的平均水平，比如比较森林、农田、草原在多年尺度上的FPAR差异，或者在碳收支模型里给参数一个合理的先验取值。第二类是GIS和遥感工程师，经常需要往系统里加一张稳定的植被生产力底图，单年数据波动太大，多年平均图更适合做背景图层和分级渲染。第三类是教学场景，老师需要一份能让学生直观理解"遥感时间序列如何变成气候态产品"的示例数据。

针对这些需求，我在产品结构里特意加了两层信息：第一层是全年平均和生长季平均，第二层是多年平均结果及其有效年数，方便不同需求各取所需。

2. 为什么选MOD15A2H：从传感器特点到版本控制

2.1 MOD15A2H的家底

MODIS是搭载在地球观测卫星平台上的中分辨率成像光谱仪，从2000年前后开始提供全球连续观测数据。MOD15A2H是它的L4植被产品之一，空间分辨率500米，时间分辨率8天，每个像元输出FPAR和LAI两个核心变量，同时附有质量评价层（FparLai_QC）。

它比早期1公里版本的同类产品分辨率更高，适合国家尺度甚至省尺度的分析。产品原始数据采用正弦投影，按瓦片分幅，存储为HDF格式。FPAR数据层里的数值实际是0—100的整数，使用时需要先除以100得到0—1的浮点数；质量层则把质量控制信息打包在整数位里，需要解码才能使用。这些看起来是很小的细节，但我在对接用户时发现，忽略它们的人相当多，后面我会专门展开。

提示：MOD15A2H的FPAR数据层单位是0—100整数，使用前务必除以100，否则所有数值会整体偏大一百倍。

2.2 为什么不是AVHRR系列，也不是其他数据源

选择数据源时我也对比过其他候选。

AVHRR的波段条件不如MODIS丰富，反演FPAR更多依赖统计关系，多代卫星之间定标差异大，虽然时间序列能够延伸到80年代，但空间分辨率和一致性不够。国产中高分辨率卫星在光谱设计、连续观测历史、业务化产品方面还有明显短板，尤其不适合做2000—2025这种25年以上统一口径的长时序反演。还有一些全球陆表参数产品虽然也提供几十年的FPAR/LAI，但它们往往把多种传感器的数据做了融合，处理链路更长，对单传感器QA细节的覆盖不足，想针对中国区域做二次质控反而更困难。

MODIS从2000年持续观测至今，产品版本稳定、算法文档清晰、QA层完整，正好覆盖标题里的2000—2025时间段。所以最终选MOD15A2H作为唯一数据源，不是因为它完美，而是它在"长序列、可追溯、分辨率够用"这三件事上最均衡。

2.3 数据获取和版本管理里最容易忽略的事

正式加工前，我把所有期次的MOD15A2H按年份下载归档，文件命名形如MOD15A2H.A2001001.hXXvXX.061.期号.hdf，其中A2001001表示2001年第001天起算的8天合成期。

这里有两个版本陷阱非常值得提醒。第一，MODIS产品有Collection版本之分，我全程只使用061版本，没有混用不同版本。不同Collection对地表反射率校正和反演参数都有更新，会让FPAR产生系统性的数值偏移，混用会让多年平均出现人为突变。第二，8天合成期次的起始日期并不总是落在自然年内，某些期次会从上一年的年末开始，按自然年归档时要按"合成期中心所在的年份"归类，而不是简单看文件夹。

我坚持把原始HDF完整保留，不删除，所有重投影、裁剪、合成的结果都归到独立目录。这样一旦发现处理参数需要调整，就能从原始数据重新开始，不用重新下载。

3. 预处理流水线：把一千多个HDF瓦片变成统一的一张图

3.1 重投影和拼接：为什么先做整景拼接再裁剪

MOD15A2H原始正弦投影是面向全球分幅设计的，直接拼接会产生坐标错位和重叠缝隙。我的做法是先用批处理脚本把每个HDF里的FPAR层和质量层分别提取出来，统一重投影到Albers等积圆锥投影，参数采用中央经线105°E、双标准纬线25°N和47°N、基准面WGS84，像元大小保持500米。

重采样方法上，FPAR层用双线性插值，质量层用最邻近法，避免插值弄乱质量标识。重投影后所有瓦片落在同一网格坐标系下，再用拼接工具合成全中国范围的单张影像。

这里有一个顺序问题：我建议先整景拼接、后统一裁剪，而不是先按瓦片裁剪再拼接。因为瓦片裁剪会引入边缘NoData条纹，后拼接时这些条纹会以不规则形状渗透进结果；整景拼接后再用统一矢量边界裁剪，能得到干净整齐的边界轮廓。

3.2 裁剪和边界像元的处理

裁剪矢量我固定使用一版稳定的中国陆地边界，全流程不换。裁剪结果里，边界外的像元统一赋NoData，边界内的水体和冰川像元先保留，留给最终用户按需求自行掩膜。

边界像元本质上是混合像元，FPAR值会有一定失真，尤其在岸线周边。所以我在产品说明里明确提醒：做区域统计时要使用内部缓冲区，或者结合土地覆盖数据剔除水体。很多人在这一步图省事，直接用原始边界裁剪后做统计，结果往往让沿海地区的平均FPAR被莫名其妙地拉低。

3.3 时间序列组织：如何让26年数据不塌方

每年大约有46个8天合成期次，26年下来单是FPAR层就有上千张栅格。我把它们整理成标准数组结构，维度是"年份×期次×行×列"，按合成期起始日期排序。这样做的好处是后续做月合成、年合成、多年平均时可以直接切片操作。

组织过程里必须完成三件事。一是检查所有影像的行列数和地理范围是否完全一致，重投影参数一旦配错，后面所有统计都会错位。二是记录每个期次的实际日期，而不是只依赖文件名。三是统计每个期次的有效像元比例，筛掉大面积无效的数据期次。这一步如果做不干净，后面合成出来的结果很可能带着一整年的空洞。

3.4 QA质量层的解码和筛选阈值

FPAR质量层（FparLai_QC）把质量控制信息编码在整数里，低三位用于表示数据质量：0表示最佳质量，1表示好质量，2表示可用质量，3表示质量降级，更高值通常表示没有生产或填充值。

解码方法是用位运算取低3位，在Python里就是qc_array & 7，然后把结果跟阈值比较。我在流水线里采用两级策略：默认只使用质量码0和1的像元；如果某个像元在时间序列里有效观测太少，再放宽到质量码2来填补。

这个策略至关重要，因为平均值合成对低质量像元特别敏感。一个被云污染的低值像元如果混进平均，会把整个像元多年平均FPAR明显拉低，而且这种偏差是系统性的，后期很难靠空间滤波彻底消除。QA筛选必须在平均之前完成，而不是先平均完再从结果里找问题。

import numpy as np def extract_fpar_by_qc(fpar_layer, qc_layer, max_qc=1): qc = qc_layer & 7 mask = qc <= max_qc fpar = fpar_layer.astype(np.float32) / 100.0 fpar[~mask] = np.nan return fpar

4. "平均值合成"的算法设计：核心逻辑和两种合成路径

4.1 合成到底在合成什么

遥感产品里的"合成"是指：在一个时间窗口内，把多期观测按某种规则合并成一张代表该窗口状态的影像。常见规则有三种：最大值合成（MVC）、平均值合成、中位数合成。

最大值合成通常用于NDVI这类噪声方向偏低的指数，取窗口内最大值能滤掉云影和部分大气影响，但结果会系统偏高。平均值合成则是把窗口内所有通过QA的像元做算术平均，得到的是窗口内"平均的植被吸收能力"，更贴合气候态语义，缺点是低质量的极低值会污染结果，因此对QA要求更高。

标题里特意强调"平均值合成"，其实就是选择了这种气候态语义。我最后采用了"先月平均、再年平均、最后多年平均"的三级平均结构，每一级都用有效像元算术平均，不用最大值。

4.2 月合成和年合成的时间口径

8天合成期次在一个月内大约有3到4景。月合成时，我先把当月的所有有效期次逐像元求平均，得到月平均FPAR；如果一个像元当月没有任何有效期次，月平均值就设为NoData。

年合成有两种口径。一种是全年平均，把12个月平均值再平均，代表全年的植被平均状态；另一种是生长季平均，取4月到10月的月平均值再平均，代表光合作用旺盛季节的平均状态。

两种口径各有用途：全年平均适合做气候态基线和模型年均参数，生长季平均更适合评估农业和自然植被在生长高峰期的实际生产能力。我把两种结果都输出，避免用户只能硬选一种。

4.3 多年平均的两条路径，我为什么选了"先做年，再做年际平均"

从2000到2025年做多年平均，可以走两条路。

路径A是先算每个月的26年平均，比如1月26年的平均值，再对12个月平均；路径B是先算每个年份的年平均，再对所有有效年份做等权平均。两者在数据完整时数学上接近等价，但数据不完整时差别很大：路径A在某个月份缺失时会把这个月权重变成0，导致年平均不均衡；路径B尊重"有数据的年份"，只要某个像元某年的年内覆盖足够，就把这一年当作一个完整样本。

考虑到2000年数据存在部分期次缺失，以及云雨区经常出现月尺度缺口，我选择路径B。表达式为：

FPAR_clim = (FPAR_year_1 + FPAR_year_2 + ... + FPAR_year_n) / N_valid_years

同时设置两个门槛：单个年份内有效月数少于6个月（生长季口径少于3个月），该年不参与平均；多年平均结果仅当N_valid_years不小于10时才输出，并额外输出一张有效年数图层。这个"有效年数图层"非常实用，用户能一眼看出哪些区域的结果更可靠。

4.4 合成前后的自检和常见异常

合成完成后，一定要看直方图和空间剖面。正常的多年平均FPAR直方图应该是连续单峰分布，峰值集中在0.3—0.6附近；如果看到接近0或接近1的地方出现异常尖峰，就要回头查QA筛选和云污染。空间上要检查不同瓦片拼接处有没有条带，区域边界处有没有不合理的突变。

另一个常见问题是水体像元天然有低值，但FPAR对水体没有明确的生理含义，直接把水体纳入统计会把区域均值往低值拉。做区域统计前一定要结合掩膜剔除水体。这一步不复杂，但能避免给后续使用者造成很大困惑。

5. 质量控制与验证：拿什么证明这份多年平均数据可以信任

5.1 三级质控：从像元到产品

我在整个流程里做了三级质控。

第一级在预处理阶段，就是前面说的FparLai_QC筛选。第二级在时间序列阶段，对每个像元的逐月FPAR做异常值检测：一旦发现数值超出合理范围（小于0或大于1），或者与前后月相比出现突兀变化超过阈值，就标记为异常，逐年统计异常比例。第三级在多年平均完成后，结合有效年数图层、直方图和剖面图检查整体格局。

三级质控的目的不是让数据变得完美，而是把不确定性的来源记录清楚。尤其是在多年平均里，一个像元如果只有两三年有效数据，跟拥有二十年有效数据的像元在置信度上完全不同。如果没有有效年数图层，使用者很容易把二者等量齐观。

5.2 与地面观测和已有产品对比

验证环节我做了两方面比较。

一方面，挑选了有长期冠层PAR观测记录的站点，把站点实测的FPAR换算结果与产品中对应位置3×3窗口内的多年平均FPAR做相关分析。我给自己定的验收线是相关系数不低于0.7、均方根误差不超过0.15，实测结果基本达标。个别山区和农田碎片区误差略大，主要来自站点代表面积与500米像元尺度的不匹配，这个可以接受但不能忽视。

另一方面，把多年平均结果与已有的全球长时序FPAR产品做空间交叉对比，统计全国尺度的均值、标准差和典型生态类型的数值区间。整体一致性不错，差异主要出现在多云雨区域（平均值偏低）和茂密森林区域（FPAR接近饱和，不同反演方法差异被放大）。这些系统偏差属于正常现象，但我会明确写进数据文档，防止用户把本产品与其它产品简单混用。

5.3 版本一致性和可复现要求

数据可复现性不是加分项，而是产品的基本要求。我保留了完整的处理脚本、参数列表和日志，包括投影参数、QA阈值、合成窗口、异常检测规则。所有脚本按顺序编号，参数集中在配置文件中，不散落在代码里。

这个习惯最初是出于偷懒，后来发现它对质量判断的帮助远大于预期。当官方更新数据集合版本，或者我决定把QA阈值从1调整为2时，可以全流程重跑一遍，快速对比参数变更前后的结果差异。没有一套完整的重跑流程，质量控制很难做得彻底。

6. 数据文件长什么样，以及使用时该注意什么

6.1 文件组织、命名和产品规格

交付时我把产品组织成这样的目录：根目录下按年份放逐年平均栅格，根目录再放多年平均栅格和有效年数栅格，另外附一份处理说明和参数文件。

文件命名尽量自解释，例如FPAR_AnnualMean_2005.tif、FPAR_ClimMean_2000_2025.tif、FPAR_ValidYears_2000_2025.tif。格式选了GeoTIFF，方便GIS直接打开；同时另出一份NetCDF版本，便于批量读取和按时间维度切片。核心规格整理成下表：

再次强调源数据除以100这个细节。很多人直接从官方整数当作0—1浮点，得到的结果全部偏大一百倍。做相对比较可能看不出来，但一旦跟站点观测、模型输出做绝对对比，问题就会立刻暴露。

6.2 适用场景、以及必须提醒的不当用法

适用场景方面，这份多年平均FPAR数据集很适合做全国尺度的植被生产力气候态制图、生态区的平均生产力对比、碳循环模型参数初始化、农田和森林长势背景评估，也能用来教学演示长时序遥感产品的处理逻辑。

不当使用的情况也有几条，我全部写进文档。不能把多年平均结果当作某一年的真实值，因为它刻意抹平了年际波动；不要在像元尺度解读局部离散变化，500米像元在山区和破碎地形里混合效应很强；不要把水体、冰川像元纳入统计而不做掩膜；不要在高FPAR区域（大于0.7）硬性比较生产力绝对差异，因为FPAR在茂密植被段会饱和，微小差异更多来自反演误差而不是真实生产力差别。

6.3 如何滚动更新和升级版本

这个数据集不需要每次都全量重算。我采用滚动更新方案：新一年的数据发布后，先按同样流程生成当年的年度平均结果，再用递推公式把多年平均更新一次：

new_mean = (old_mean * old_n + new_value) / (old_n + 1)

同时把有效年数计数加一。这样每年只需要处理一年数据，计算量可控。

但要注意：如果官方发布了新的Collection版本，或者我调整了QA阈值、合成窗口，就不适合滚动更新了，必须从原始HDF开始全序列重处理，否则新旧数据语义不一致，整个序列的一致性会被破坏。我把"什么情况用滚动更新、什么情况必须全量重算"写成决策规则，挂在项目文档首页，避免半年之后自己都拿不准。

7. 做完之后，最想提醒后来人的三件事

7.1 平均合成最怕混入低质量像元

平均值合成没有最大值合成那种天然容错能力。最大值合成因为取大值，会把污染向低值方向压掉；平均值合成只要混入一个被云污染的低值像元，多年平均就会偏低，而且这种偏差经年累月、很难排查。所以做平均合成的第一原则是：QA不达标宁可缺，也不硬凑；缺得太多宁可放宽阈值并在文档里记录，也不要无差别平均。

7.2 时间口径一定要写清楚

同一个FPAR均值，全年平均、生长季平均、月气候态平均，数值差别很大。使用时如果不澄清口径，结果对比就会出问题。我在文档和文件名里都明确标了口径，但还是建议用户在使用前先想清楚自己的科学问题到底需要哪个口径，再决定用哪一层数据。

7.3 脚本和参数要当成产品的一部分来管理

这类数据集的价值不在于算法多新，而在于可复现和口径一致。如果没有版本管理和参数记录，三年后很可能连自己都说不清楚哪个文件用的是哪个QA阈值、哪一次合成用了哪种投影，所谓的质量控制也就无从谈起。把脚本、参数、日志当作产品的组成部分一起归档，是我做完这个项目后最想分享的一条经验。