这个标题乍一看像数据集仓库里的一条归档记录但内行会立刻意识到背后有一套完整的遥感估算链路NPP净初级生产力是草原碳循环研究的核心指标图门索格特是蒙古国东部的典型草原区1982-1990则是一个很有讲究的时间窗口——正好对应 NOAA AVHRR 全球遥感数据序列的起步阶段R1 又代表了这个产品的第一轮正式版本。当初我拿到这个任务时最想弄清楚的是三件事为什么选这片草原、为什么卡在 1982 到 1990 这几年、以及 R1 版本做了什么质量控制。这篇博文就把我在这套 NPP 草原产品上的研究过程和实操经验完整写出来内容包括区域选择逻辑、数据链路搭建、CASA 模型参数调优、结果时空规律解读以及版本迭代反思适合正在做长时序植被生产力遥感反演、或者想复现蒙古高原草原 NPP 估算的研究者和工程师参考。1. 为什么盯着图门索格特草原NPP研究的区域选择与数据窗口1.1 一个蒙古东部苏木的科学代表性图门索格特Tumentsogt是蒙古国苏赫巴托尔省的一个苏木相当于国内的乡镇级行政区位置大概在东经 112 度附近、北纬 45 到 46 度之间属于蒙古高原东部典型草原亚地带的核心范围。这类区域在遥感研究里经常被视为蒙古国草原的“平均脸”植被组成以克氏针茅、羊草、冷蒿为主草层高度不高盖度中等年降水量在 200 到 250 毫米左右并且集中在 7 到 8 月。生态学上这种典型草原正好处在森林草原向荒漠草原过渡的中间梯度上对气候波动特别敏感NPP 的年际变化幅度远大于湿润区是研究干旱半干旱生态系统生产力响应机制的天然实验场。选择图门索格特作为研究区不是一个随意的决定。我在做区域筛选时主要看重三个条件一是草地类型要有代表性不能选过渡带太窄的特殊点位二是地面资料相对可获取蒙古国 1980 年代的畜牧业统计和气象站点数据虽然稀疏但苏木级别的牲畜量、草场利用记录仍然能够找到这对遥感估产结果的验证非常关键三是像元尺寸与研究区面积的匹配度GIMMS NDVI 数据的空间分辨率是 8 公里如果一个苏木的面积太小混合像元问题会很严重而图门索格特及其周边几百平方公里的连续草原刚好能容纳足够多的纯像元统计意义上的信噪比是够的。这里有一个常被初学者忽略的尺度问题NPP 遥感反演结果并不代表像元中心那一小块草地的生产力它代表的是 8 公里×8 公里范围内的植被生产平均水平。换成图门索格特这样的区域结果的可解释性就体现在“区域尺度生产力”上而不是某个样方尺度的数值。这也是为什么我在后面做验证时不会拿单点实测生物量直接对比栅格值而是先做样方数据的空间聚合。1.2 1982-1990年这个窗口期的生态与数据双重意义1982 到 1990 这个时间段从数据角度讲有非常明确的含义GIMMS NDVI 数据集就是从 1981 年 7 月正式开始发布的1982 年是第一个完整日历年的开始。这是全球范围内最早能够支持长时序植被监测的卫星遥感产品来源是 NOAA 系列气象卫星上的 AVHRR 传感器。MODIS 是 1999 年底才发射的Landsat 虽然更早但重访周期 16 天不要说蒙古国东部这种云雨不稳定区域就是放到全球尺度也没法在 8 年尺度上逐旬拼接出干净的草原 NDVI 时序。换句话说想研究 1980 年代草原生产力动态AVHRR GIMMS 几乎是唯一可用的长时序数据源。从生态学角度看1980 年代蒙古高原草原并不平静。这段时期蒙古国东部经历了明显的干湿波动尤其是 1986 年前后出现过一次较大范围的降水低谷对应草原生产力的显著下降。这样的气候事件在 NPP 时间序列里会留下清晰的痕迹非常适合用来检验模型在干旱胁迫下的表现。更实际的原因是这一时期的牲畜存栏量和草场利用方式有连续的官方统计记录我可以用牲畜量的变化趋势与 NPP 反演结果做交叉验证判断两者是否存在合理的耦合关系——丰水年草多、载畜量上升枯水年草少、牲畜量往往受压这种联动关系如果能在数据里看到说明 R1 产品的时序质量是可信的。关于 R1 这个标识我理解它既是一个数据产品版本号也代表了一整套处理流程的第一轮交付。R1 版本在我这里意味着数据定标、影像合成、模型反演、基础质量控制全部跑通并且产出了完整的 1982-1990 逐年 NPP 栅格序列但还留有后续参数优化和验证迭代的空间。后面第五节我会专门讲 R1 版本做了哪些质量检查以及我在使用过程中发现并修正的问题。2. 从AVHRR原始影像到NPP栅格完整数据链路拆解2.1 数据源选型为什么一定是AVHRR GIMMS NDVI先明确一个概念NPP 不是观测出来的而是模型估算出来的。卫星能直接测到的只是地表反射率经过大气校正和植被指数计算之后得到 NDVI再通过光能利用率模型把 NDVI 换算成干物质产量。所以整个链路的数据源头就是 NDVI 时序数据。在 1982 到 1990 这个区间可选的 NDVI 产品非常有限。GIMMS NDVI 数据集有几个别的数据替代不了的优势时间跨度从 1981 年 7 月一直到 2015 年以后中间没有断层8 公里分辨率和 15 天合成周期足以捕捉草原物候的年际变化最关键的它针对 AVHRR 传感器退化做了交叉定标处理了轨道漂移问题这在长时序分析里是生死攸关的一步——如果不做传感器退化校正NDVI 长期趋势里会混入大量虚假下降信号NPP 的年份间波动就完全失去可信度。我选择的是 GIMMS NDVI3g 版本。需要说明的是NDVI3g 有 R0 和 R1 两个版本子集R1 是更新后的再处理版本对北极地区和干旱区的异常值做了额外修正。我拿到数据后第一步不是直接算 NPP而是先做数据体检把研究区内所有像元的 NDVI 最大值、最小值、缺失比例、时序曲线形态都画出来检查有没有异常跳变、天花板效应和持续的零值。这一步虽然枯燥但能省掉后面大量排查问题的时间。2.2 预处理三件套投影统一、最大值合成、时序去噪先解决投影问题。GIMMS NDVI 原始数据采用的是经纬度网格但后续要做面积统计、与矢量边界叠加、按像元输出逐年统计值经纬度网格在不同纬度上代表的面积不同直接算会有偏差。我统一转成了阿尔伯斯等积圆锥投影Albers Conical Equal Area中央经线取东经 105 度标准纬线取北纬 30 度和北纬 60 度这样研究区范围内的像元面积误差控制在很小的范围内栅格统计和面积加权都更方便。第二步是合成。GIMMS 本身提供的是 15 天最大合成 NDVIMVC即每半个月取该时段内所有可用观测的最大值目的是减少云污染影响。理论上可以直接用这个产品但实际操作中我仍然做了二次筛选把生长季定义为每年的第 129 天到第 273 天约 5 月中旬到 9 月底只保留生长季内的合成值参与 NPP 累积。蒙古草原冬季基本是枯草期NDVI 背景值很低如果不做生长季截取冬季虚假的微弱绿度也会被模型折算成少量 NPP导致年总量被系统性抬高。第三步是时间序列平滑和去噪。虽然最大值合成已经削弱了云的影响但残留的云、气溶胶和传感器噪声仍然会造成 NDVI 时序的异常凹陷。我用了 Savitzky-Golay 滤波做时序重建窗口取 7 个合成期多项式阶数取 3。这组参数不是拍脑袋定的而是反复测试的结果窗口太短去不掉噪声太长又会把真实的物候峰拉平。对蒙古草原这种单峰型生长曲线来说窗口 7 基本能在保留峰值的条件下抹平大部分异常点。平滑之后还有一个兜底操作把 NDVI 小于 0.05 的像元直接标记为无效这些像元通常对应水体、云影或极端退化裸地参与模型计算只会产生噪声。2.3 CASA模型让NDVI变成生产力的核心逻辑CASACarnegie-Ames-Stanford Approach模型是光能利用率模型的经典代表核心公式只有一行NPP APAR × εAPAR 是植被吸收的光合有效辐射ε 是光能利用率。APAR 又可以拆成两个部分入射的光合有效辐射 PAR 和植被吸收比例 FPAR其中 FPAR 是从 NDVI 反推出来的。整套链路可以描述为用 NDVI 通过线性关系计算 FPARFPAR 乘以太阳总辐射和光合有效辐射比例得到 APARAPAR 乘以受温度、水分调节的光能利用率 ε得到该时段栅格的 NPP生长季内逐时段累加得到年 NPP。这套模型的优点在于参数少、输入简单、无需叶片级别的生理参数非常适合大尺度草原应用。但正因为模型结构简单每个参数的取值都直接影响最终结果的绝对量级。比如 FPAR 与 NDVI 的线性关系两端需要确定 NDVI 最小值、最大值对应的 FPAR 最小值、最大值不同的取值会让 NPP 差出百分之十几。我在实际处理时FPARmin 取 0.001对应裸地FPARmax 取 0.95对应浓密植被NDVImin 取研究区内所有像元的 5% 分位数NDVImax 取 95% 分位数这样避免极端像元把两端拉伸得过宽。光能利用率 ε 的核心公式是ε εmax × Tε × WεTε 是温度胁迫系数Wε 是水分胁迫系数两者取值都在 0 到 1 之间干旱严重时 Wε 会显著压低 ε从而让纯气候驱动的“潜在生产力”变成受土壤水分限制的“实际生产力”。这是 CASA 模型在干旱半干旱草原表现好于简单 Miami 模型的关键原因。关于 εmax 和其他参数的取值我在第三节详细展开。注意CASA 输出的 NPP 单位通常是 gC/m²/半月也就是单位面积碳固定量而不是生物量干重。如果后续要折算成牧草产量还需要额外乘以一个碳含量系数并考虑地下生物量分配但作为生态生产力指标直接在碳单位下讨论就可以了。这是我在给非遥感背景的合作者解释结果时反复强调的一点很多人把 NPP 直接理解成“草有多少斤”其实是把概念简化过头了。3. CASA模型在干旱草原的参数调优我反复试验出来的组合3.1 最大光能利用率εmax到底该取多少εmax 是 CASA 模型里影响全局的敏感参数它代表理想条件下植被把光能转化为有机碳的最高效率。Potter 等人在提出 CASA 模型时通常把全球尺度的 εmax 定在 0.389 gC/MJField 等人在后续研究中也沿用类似取值。但我在做图门索格特草原时对这个值做过针对性分析这套模型最初是基于全球植被平均状况标定的包含了森林、农田、灌丛等多种高生产力植被类型而干旱原草在自然状态下存在长期水分亏缺很难达到温带森林那样的光能转化效率。实际操作上我试过三组取值0.389全球默认、0.342部分学者在干旱区推荐的取值、0.30纯干旱区保守值。对比结果很有意思不同取值对 NPP 年总量的影响接近线性0.342 组合得到的多年平均 NPP 大约在 210 gC/m²/年0.389 组合大约在 245 gC/m²/年相差大约 15% 到 17%。由于图门索格特区域没有足够多的实测碳通量数据做绝对标定我最终没有强行修改 εmax而是保留 0.389 作为基准通过 Wε 水分胁迫系数去压低干旱期的实际光能利用率。这样处理的理由是εmax 的原始概念是“潜在最大光能利用率”本身就代表一种理想上限实际生产力应通过胁迫系数来折减而不是预先调低上限。保留统一 εmax 也方便和文献中其他区域的 CASA 结果直接做横向比较。如果你做的是某个具有大量实测通量数据的站点完全可以反过来用实测 NPP 去标定 εmax但在区域尺度、历史时段、无通量塔的背景下保持模型原参数并用胁迫项承担区域差异是更稳妥也更可复现的策略。3.2 温度与水分胁迫系数的简化处理思路温度胁迫系数在 CASA 模型里由两部分构成第一部分 Tε1 反映低温对光合作用的抑制第二部分 Tε2 反映温度偏离最适温度时对光能利用率的影响。蒙古草原生长季平均温度在 15 到 20 摄氏度之间接近 CASA 模型设定的最适光合温度因此我最开始计算的 Tε 数值明显偏大长年接近 0.9 以上。考虑到图门索格特草原生产力变化的主导因子是降水而非温度在 R1 版本里我没有对温度胁迫做过度精细的刻画只保留了一个简化逻辑月均温低于 0 摄氏度的月份直接不参与生长季累积生长季内温度默认接近最适区间Tε 取常数值 1。这样做的本质是把温度胁迫从限制因子里剔除出去让模型对水分变化的响应更加敏感。水分胁迫系数 Wε 的处理才是草原 NPP 反演的重头戏。CASA 原始文档中 Wε 通常用实际蒸散与潜在蒸散之比EET/PET来估算但这个公式依赖逐月气象栅格数据而在蒙古国东部 1980 年代的气象站点非常稀疏插值误差很大直接用蒸散数据反而会引入不可控的不确定性。我采用了一个替代方案利用同期降水量数据计算湿润指数用湿润指数作为 Wε 的输入。具体做法是获取研究区周边气象站点的逐月降水观测用反距离加权插值到每个 NDVI 像元中心计算每月的湿润指数即当月降水与潜在蒸散的比值湿润指数大于 1 时 Wε 取 1无水分限制小于 1 时按线性比例折算。这个替代方案牺牲了一部分物理严格性但换来了两个好处一是降水的观测资料比蒸散资料容易获取得多二是在蒙古草原这种降水与生产力高度相关的区域湿润指数本身已经能解释很大一部分 NPP 波动。当然如果你手上正好有高质量的 ERAS 或类似再分析蒸散数据直接用原始公式会更好R1 版本为了保证可复现性和数据可获取性选择了前者。3.3 针对典型草原植被特点的FPAR修正FPAR 的估算精度直接影响 APAR 的计算而 FPAR 从 NDVI 反演时有一个经典的难题当植被盖度较低时土壤背景反射会显著“污染”NDVI 信号让 FPAR 出现偏低或偏高的系统偏差。图门索格特地区的草原在干旱年份整体盖度可能降到 30% 以下这个问题非常突出。我采用的方法是引入土壤调节植被指数 SAVI 做中间变量而不是直接用 NDVI。SAVI 在 NDVI 基础上增加了一个土壤调节因子 LL 的取值取决于植被密度中低盖度时取 0.5 比较合适。计算 SAVI 后用 SAVI 的线性缩放取代 NDVI 的线性缩放来估算 FPAR这样可以有效降低土壤背景的干扰。另外针对蒙古草原冬季地表裸露比例高、生长季初期的 NDVI 容易受到枯草残留影响的特性我还在 FPAR 计算前增加了一层背景值扣除取每年第 1 到第 2 个合成期的 NDVI 平均值作为该像元的土壤加枯草背景值从生长季 NDVI 中减去后归一化再参与 FPAR 计算。这个处理措施使生长季早期 NPP 累积量更加合理枯草残留造成的伪生产力明显减少。R1 版本里这部分是手动分像元批量处理的虽然效率不高但对结果质量有实质性提升。下面是 R1 版本实际使用的 CASA 模型关键参数汇总参数R1取值说明εmax0.389 gC/MJ保持CASA全球默认值便于横向比较FPARmin/FPARmax0.001 / 0.95对应裸地与浓密植被上限NDVImin/NDVImax5% / 95% 分位数按研究区像元统计确定温度胁迫Tε生长季内取1月均温低于0℃不参与累积简化处理突出水分限制水分胁迫Wε基于降水湿润指数折算替代EET/PET适应稀疏气象站条件土壤调节因子L0.5SAVI计算参数4. 1982-1990年图门索格特草原NPP的时空结果解读4.1 年际波动从干旱年与湿润年的对比看模型敏感度R1 版本输出的 1982-1990 逐年 NPP 空间分布图叠加平均后能看出几个非常典型的草原动态特征。首先是年际波动幅度大多年平均 NPP 大约在 180 到 230 gC/m²/年之间但最湿润年份和最干旱年份的差值常常超过一倍。1985 和 1988 年前后相对湿润NPP 明显抬升1986 到 1987 年出现显著低谷空间上几乎全区域同步下降。这种高度同步的年际波动说明模型对水分信号非常敏感降水的年际变化通过 Wε 和 FPAR 两条途径共同放大了 NPP 响应。我把逐年 NPP 与生长季累计降水做了相关分析决定系数大约在 0.72 左右在生态遥感模型里这样的相关水平已经算比较理想。这让我有底气认为 R1 版本的核心驱动逻辑是健康的。另外还有一个有意思的现象NPP 对前一年降水的滞后响应并不明显。蒙古草原的土壤蓄水能力有限春季融雪加上夏季降水基本决定了当年生产力土壤深层水分的跨年缓冲作用远比黄土高原等深厚土层区域要弱。这解释了为什么当年降水与当年 NPP 的相关性如此直接。4.2 空间格局地形梯度与草地利用的叠加影响空间分布上图门索格特区域的 NPP 呈现明显的西北高、东南低的格局。西北方向靠近肯特省边缘海拔略高地形抬升带来更多地形雨草场长势更好东南方向则逐渐过渡到荒漠草原带降水减少盖度下降NPP 低值区连片出现。R1 版本中这种梯度很稳定几乎在各年份都能复现说明它不是偶然噪声而是区域水热条件的真实反映。另一个值得注意的空间特征是河流谷地和洼地附近的高值带。这些区域地下水埋深浅即便在干旱年份也能维持较好的土壤湿度草地生产力抗波动能力明显强于坡地和梁地。这在实际应用中很有价值如果这个 NPP 产品用于草场承载力量化那么就不能简单按全区域平均值配给载畜量而是需要根据空间异质性分区管理。我在结果整理时还做了一件额外工作按 NPP 数值将研究区分为高、中、低三个等级统计每个等级的面积占比。结果显示1986 年干旱年份高等级面积占比比 1985 年大约下降了 60% 以上这类空间统计对草原管理决策的参考价值比单看平均值直观很多。5. R1版本的质量控制与我的验证心得5.1 R1版本做了哪些质量控制从像元到年尺度的多层检查R1 版本虽然只是第一轮交付但质量控制环节我并没有省。整套检查分四个层次第一层是单像元时序检查。逐个像元输出 1982-1990 年的 NDVI 和 NPP 时序曲线自动标记异常跳变点、负值和长时间零值。对于蒙古草原来说冬季 NDVI 掉到 0 甚至负值是正常的但生长季内出现负值基本可以判定为残留云污染需要回溯原始数据。第二层是空间连续性检查。逐年计算 NPP 栅格的空间自相关并用前后年份的栅格做差值图寻找异常斑块。如果某一年突然出现一个和周围环境完全孤立的低值或高值区很可能对应的是该年份影像质量问题而非真实的草原变化。第三层是年际合理性检查。把每个像元的 9 年 NPP 序列做趋势检验看是否存在突变点和系统性漂移。AVHRR 传感器退化在 1980 年代后期逐渐明显如果不处理NDVI 会在 1988 年以后显示出虚假的下降趋势进而压低 NPP。GIMMS 数据集本身做了定标校正但 R1 版本仍需要检查趋势是否符合降水趋势——如果降水轻微上升而 NPP 持续下降那八成是数据问题而不是生态问题。第四层是与外部资料的交叉验证。当年 NPP 与牲畜存栏量变化的对比方向以及不同草地类型之间的 NPP 差异是否符合生态常识都是判定结果可信度的辅助证据。5.2 与实测和间接数据的对标方法验证 NPP 绝对数值是一件困难的事。图门索格特地区 1980 年代几乎没有连续的涡度通量观测所以我采用了几种间接对标方式。第一种是和区域生物量实测文献值对比相关研究对蒙古东部典型草原的地上生物量测定值通常在 80 到 200 g/m²干重之间把 NPP 换算成干物质时需要除以约 0.45 的碳含量系数两者量级基本吻合。第二种是引用同区域、同模型族的公开发表 NPP 结果做横向对照我的 R1 多年均值与文献中蒙古高原典型草原 150-250 gC/m²/年的区间高度一致。第三种是定性时间对齐检验1986-1987 年的低值期与牲畜量拐点年份有较好的对应关系。这几种验证方式虽然没有通量塔那么“硬核”但在历史数据稀缺条件下已经足够支撑 R1 版本的基本可信度。需要强调的是任何区域尺度 NPP 遥感产品的精度都是相对的验证的本质不是证明绝对正确而是证明“偏差在可接受范围且方向合理”。5.3 我在迭代中踩过的坑与修正记录R1 版本交付后我重新复盘时发现过几个明显的问题。第一个问题是冬季背景值扣除不够彻底研究区最西端靠近河岸林的像元在非生长季仍产生了少量正 NPP导致个别年份年总量系统性偏高大约 8%到12%。这个偏差在高盖度像元上尤其明显后来我在 FPAR 计算中增加了一层冬季 NDVI 最小值背景扣除问题才基本解决。第二个问题是 1986 年部分月份的降水插值在东南方向出现了一个虚假高值区原因是一个远端气象站在该年份记录异常反距离加权插值时把这个异常带到了附近区域导致对应的 NPP 高值斑块失真。后来我对气象数据加入了极值检验和站点一致性检查才避免了同类问题。第三个问题相对隐蔽GIMMS NDVI 的最大合成值在蒙古高原春季容易受到融雪影响雪盖存在时 NDVI 被压得很低但雪融后很快恢复正常。这个信号幅度很小但如果研究区里面有持续到 5 月的积雪区域春季 NPP 会被低估。处理方案是把融雪窗口期大约每年第 121 到第 145 天的 NDVI 替换为前后时段的线性插值而不是直接用原始值参与模型计算。这三个坑让我深刻意识到长时序遥感反演的质量问题很少出在模型公式本身更多是出在数据预处理和外部辅助数据的可靠性上。R2 版本如果再更新我会重点优化降尺度方案、引入更多实测样方数据进行绝对验证以及在水分胁迫系数中尝试融合土壤湿度再分析数据。最后再分享一个小技巧处理这类长时序 NPP 产品时养成保留中间产物的习惯。NDVI 平滑序列、FPAR 栅格、APAR 栅格、逐时段的 NPP 栅格都分目录存放一旦最终结果出现问题回溯定位会快得多。我只在 R1 版本里保存了最终年累积结果排查冬季背景值问题时不得不重跑全流程白白浪费了两天时间。这种细节做过一次的人再也不想做第二次。