这个系列写到第十二篇我们来聊生态遥感里绕不开的一个话题NPP反演。NPP的英文全称是Net Primary Productivity净初级生产力指的是绿色植物在单位时间、单位面积内通过光合作用积累的有机物质总量扣掉植物自身呼吸消耗后真正留下来的那部分。说白了就是一片区域在单位时间内被植物“固定”下来的碳。这个参数在碳循环研究、区域生态系统健康评估、作物长势监测里都是底层基础数据所以很多同学做毕业设计或者科研项目时都会碰到“反演NPP”这个需求。NPP的遥感反演方法很多接触最多的还是光能利用率模型而这个家族里最经典的一支就是CASA模型全称Carnegie-Ames-Stanford Approach。经典CASA模型直接用空间尺度和参数取值往往跟实际研究区对不上所以现在的主流做法都是做“改进CASA模型”。这篇文章我把自己在ENVI里基于改进CASA模型反演NPP的完整流程整理出来包含原理拆解、数据预处理、波段运算表达式、参数取值和踩坑记录适合正在做植被生产力、碳收支或者生态遥感参数反演相关研究的同学参考。我用的数据是Landsat 8 OLI空间分辨率30米研究区尺度在县域到市域之间这套流程换影像源、换参数也能直接平移。1. 先把原理啃透经典CASA模型的骨架与短板1.1 经典CASA模型的数学框架CASA模型的出发点是光能利用率理论。它认为NPP等于植被吸收的光合有效辐射APAR乘以实际光能利用率ε用公式表示就是NPP APAR × ε这里APAR的全称是Absorbed Photosynthetically Active Radiation即光合有效辐射中被植被吸收的部分单位通常是MJ/m²·月ε是实际光能利用率单位gC/MJ。两者相乘NPP的单位就是gC/m²·月。APAR这部分可以再拆成三步APAR SOL × FPAR × 0.5其中SOL是太阳总辐射FPAR是植被冠层对光合有效辐射的吸收比例0.5这个系数表示在400到700纳米的可见光波段里真正能被叶绿素吸收的光合有效辐射大约占太阳总辐射的一半。这个0.5不是随便拍的而是光合有效辐射在太阳短波辐射中的经验比例。FPAR在经典CASA里通常用两种植被指数来约束一种是NDVI另一种是比值植被指数SR分别算出两个FPAR后再取平均FPAR_NDVI (NDVI - NDVI_min) / (NDVI_max - NDVI_min) × (FPAR_max - FPAR_min) FPAR_minFPAR_SR (SR - SR_min) / (SR_max - SR_min) × (FPAR_max - FPAR_min) FPAR_minSR (1 NDVI) / (1 - NDVI)FPAR (FPAR_NDVI FPAR_SR) / 2这里FPAR_max和FPAR_min分别取0.95和0.001这两个数来自前人研究中对植被冠层吸收效率上下限的估计几乎是CASA模型里的固定值。NDVI_max和NDVI_min则要根据研究区植被类型和植被指数取值范围来定这一步也是后面“改进”的重点之一。实际光能利用率ε在经典CASA里是一个被温度和水分“打折”后的量ε ε_max × T_ε × W_εε_max是最大光能利用率也就是理想条件下植物把光能转化成有机碳的效率T_ε是温度胁迫系数W_ε是水分胁迫系数。两个系数的取值范围都在0到1之间分别表示低温和干旱对光合作用的削减程度。1.2 经典CASA模型在实际应用里的三个短板经典CASA模型很经典但直接照搬到自己的研究区至少有三个现实问题。第一个问题是数据源太粗。经典CASA当初是基于AVHRR NDVI数据构建的空间分辨率8公里甚至更粗做全国、全球尺度的NPP估算没问题但拿到县域、市域这种中小尺度几乎没法用。现在大家手上都是Landsat、Sentinel-2这种10米到30米分辨率的影像拿8公里分辨率的结果去跟30米土地利用数据做叠加分析颗粒度完全不匹配连图都画不到一块。第二个问题是水分胁迫系数W_ε在经典模型里依赖区域实际蒸散量和潜在蒸散量的比值也就是W_ε 0.5 0.5×E/E_p。这个式子本身没问题但E和E_p需要气象站点的逐日观测数据来插值而且蒸散量计算通常要套Penman-Monteith公式里面又涉及温度、湿度、风速、日照时数等一堆参数。研究区范围一缩小气象站点往往只剩一两个插值出来的空间异质性很差等于变相削弱了遥感影像的空间优势。第三个问题是最大光能利用率ε_max是一个“均质化”参数。经典CASA在全球范围取一个统一的ε_max比如0.389 gC/MJ但不同植被类型的光能利用率差异很大针叶林、阔叶林、草地、农田根本不在一个水平上。如果用全球统一值硬套一个区域结果偏差会非常明显。搞清楚这三个短板“改进”的方向也就自然出来了。2. 改进方案怎么设计三个改动点逐一说明2.1 改进方向一用Landsat 8 OLI影像替代AVHRR数据源这次反演用的是Landsat 8 OLI。Landsat 8有11个波段其中红波段是Band4波段范围0.64到0.67微米近红外是Band50.85到0.88微米短波红外Band6是1.57到1.65微米Band7是2.11到2.29微米。这个波段配置对NPP反演来说非常理想因为算NDVI只需要Band4和Band5算LSWI水分指数只需要Band5和Band6这些在经典AVHRR数据上要么没有要么波段范围太宽怎么说都绕不过去。空间分辨率方面Landsat 8可见光和近红外波段是30米做县域尺度的NPP反演颗粒度很合适既不会像MODIS 250米那样丢失细节又比无人机影像覆盖范围大很多。时序上Landsat 8是16天重访周期单景影像基本能覆盖一个中型研究区做月尺度或者生长季尺度的NPP估算完全够用。高分辨率数据带来的另一个好处是可以把NPP结果直接跟土地利用分类、植被覆盖度这些产品叠加分析做空间归因。我后面验证结果时把NPP异常值点叠到土地覆盖图上能很直观地看出水体、建设用地和农田的NPP差异这种操作放在8公里分辨率的经典数据上想都不敢想。2.2 改进方向二水分胁迫系数W_ε改用遥感指数计算这是“改进CASA模型”里最核心的一个改动。经典模型里的W_ε需要算蒸散比E/E_p麻烦不说空间精度还吃亏。改进方案里用陆地表面水分指数LSWI来替代蒸散比公式变成W_ε (1 LSWI) / (1 LSWI_max)LSWI (ρ_NIR - ρ_SWIR) / (ρ_NIR ρ_SWIR)ρ_NIR是近红外波段反射率ρ_SWIR是短波红外波段反射率。从物理意义上讲水分子在短波红外波段有很强的吸收特征植物叶片含水量越高、冠层越湿润短波红外的反射就越弱同时近红外反射率和叶片结构、生物量密切相关。所以LSWI这个比值能很好地反映土壤-植被系统的水分状况用它来代替气象站点插值得到的蒸散比既合理又方便。这个改进带来的直接收益是完全绕开了蒸散量计算只需要Landsat影像的近红外和短波红外两个波段就能得到逐像元的水分胁迫系数空间精度和影像分辨率完全一致。在ENVI里算LSWI只需要一步Band Math(float(b5)-float(b6))/(float(b5)float(b6))前提是输入影像已经做好大气校正、反射率定标。LSWI_max的取值可以直接统计研究区内全部有效像元LSWI的累计95%或99%分位数不需要查表这也是本地化改进的一部分。2.3 改进方向三关键参数本地化标定本地化标定主要体现在两个参数上最大光能利用率ε_max以及FPAR公式里的NDVI_max和NDVI_min。ε_max按植被类型分开设。经典CASA用0.389这个全球平均值区域应用时偏差很大。根据国内大量文献在区域尺度的研究结论几种常见植被类型的取值为针叶林0.389阔叶林0.485灌丛0.429草地0.542农田0.542其他植被类型按0.389左右处理。这些值不是拍脑袋编的是很多研究者在不同区域用站点观测NPP校准后得到的结果范围。如果你的研究区有通量塔实测NPP数据最优做法是用实测数据反算ε_max并做对比没有实测数据就用文献值并写清楚来源。NDVI_max和NDVI_min我的做法是在计算出的NDVI栅格上做统计取累计概率1%和99%的分位数作为上下限。这样做的好处是不需要额外准备植被类型图也避免了直接套用0.9和0.02这种固定值带来的误差。选1%和99%还有一层考虑能剔除云污染、水体残留等异常像元的干扰比直接用影像最大值和最小值稳健得多。实际上水体像元NDVI经常在0以下云和云影像元有时候高达0.9以上这些极端值如果不处理会把FPAR的线性拉伸整体带偏。3. 动手前的数据准备与预处理流程3.1 数据清单影像、气象、辅助数据一样都不能少反演NPP不是一张影像就能搞定的完整的数据清单大致分三块。第一块是遥感影像。我用的是Landsat 8 OLI Level-1产品可以从USGS EarthExplorer或国内的一些地理空间数据平台下载。尽量选择研究区云量低于10%的时相如果做生长季平均NPP最好挑6月到9月之间2到3期无云影像分别反演后再求平均。还要注意如果研究区横跨Landsat两条轨道相邻轨道之间有重叠区要提前规划好是拼接还是裁剪别等到算NDVI的时候才发现两景影像范围对不上。第二块是气象栅格数据。CASA模型需要太阳总辐射SOL、月平均气温T_mean和最适温度T_opt。SOL可以采用已有的辐射产品或者气象站点插值T_mean可以用气象站点月平均气温做空间插值反距离权重或样条插值都可以T_opt的物理含义是植物光合作用最适宜的温度一般取研究区生长季内NDVI最大月份的月平均气温也可以近似取研究区6到8月月均温的最高值。气象数据最终都要重采样成和Landsat影像相同的投影、相同30米分辨率、相同行列数这一步在ENVI里用Resample Data工具完成。第三块是辅助数据。研究区矢量边界用来裁剪和最终掩膜比如行政区边界或流域边界如果要做分植被类型设定ε_max还需要土地覆盖数据比如GlobeLand30或MODIS的MCD12Q1产品。3.2 辐射定标与大气校正这一步偷懒后面全是白做NPP反演对反射率精度要求很高因为所有植被指数、FPAR、LSWI都基于反射率计算反射率偏一点最终NPP偏差会被链路上的公式逐级放大。所以预处理里辐射定标和大气校正必须做扎实。辐射定标在ENVI里的操作路径是Toolbox → Radiometric Correction → Radiometric Calibration把DN值转成表观反射率或辐亮度。定标时注意两个选项一是Output Interleave选BIL或者BSQ都行但后面处理要保持一致二是Scale Factor要选对Landsat 8 Collection 2产品的定标系数里带有缩放因子设置不对会出现所有像元值被放大10000倍的情况这是新手最容易踩的坑。定标后先看一眼影像统计值如果最大值大得离谱基本就是Scale Factor没设对。大气校正推荐用FLAASH模块路径是Toolbox → Radiometric Correction → Atmospheric Correction Module → FLAASH Atmospheric Correction。FLAASH输入要求是辐亮度数据格式必须是BIL单位是μW/(cm²·nm·sr)如果不是这个单位要先做转换。大气模型根据研究区纬度和季节选择比如中纬度夏季选Mid-Latitude Summer气溶胶模型选Rural能见度初值设40公里左右如果有实测能见度就用实测值。这一块计算量比较大30米影像全幅跑FLAASH可能要十几分钟到半小时属于正常情况耐心等就好。FLAASH做完后还要进行裁剪和掩膜。裁剪用Subset Data from ROIs掩膜用Build Mask把水体、云影、建设用地这些非植被像元排除掉。特别是水体LSWI在水体区域会异常高如果不掩膜W_ε会被带偏。掩膜时务必保证掩膜范围和原影像严格对齐投影、行列数都要一致否则后面Band Math里两个变量维度不匹配会直接报错。3.3 气象栅格数据的插值与重采样气象数据处理的要点是“先插值后重采样再对齐”。先把气象站点数据通过ArcGIS或ENVI里的插值工具生成连续表面然后统一投影到WGS84 UTM或Albers等积投影最后重采样成和Landsat影像一致的30米像元大小和行列数。顺序不能反如果先把站点数据重采样再插值会引入无意义的空洞。T_opt这个参数要特别说一句。它的字面意思是最适温度但CASA模型里的T_opt更偏向“让植被光合能力达到最强的那个温度”。实际操作中我用的是研究区生长季6到9月NDVI最高的那一期影像所对应月份的月平均气温把它做成全区域一个常量栅格。如果研究区地形高差很大比如山地丘陵区不同海拔的温度差异明显可以把T_opt按海拔做分区插值平原地区通常没必要这么精细。4. 核心环节实现在ENVI中一步步算到NPP4.1 第一步NDVI与SR准备好大气校正后的反射率影像先算NDVI。在ENVI的Toolbox里找到Band Math输入表达式(float(b5)-float(b4))/(float(b5)float(b4))变量b5对应近红外波段b4对应红光波段。计算完成后检查结果范围正常情况下NDVI应该在-0.1到0.95之间含云、水、雪像元时会出现负值。如果整体数值异常比如最大值超过1或者很多像元在1.5以上基本可以断定是辐射定标或大气校正出了问题回头查别继续往下算。接下来算SRSR (1NDVI)/(1-NDVI)在Band Math里输入这个表达式时变量可以命名为NDVI选择之前算好的NDVI文件。NDVI等于1或接近1的像元会导致分母为0但正常情况NDVI不会达到1如果有异常像元可以先对影像做取值范围裁剪把超过合理范围的像元去掉。4.2 第二步FPAR_NDVI、FPAR_SR与FPARFPAR_NDVI用线性拉伸公式计算。假设通过统计得到NDVI_min0.05NDVI_max0.85FPAR_max0.95FPAR_min0.001则表达式为(NDVI-0.05)/(0.85-0.05)*(0.95-0.001)0.001FPAR_SR的表达式类似先算SR再代入(SR-SR_min)/(SR_max-SR_min)*(0.95-0.001)0.001SR_min和SR_max可以由NDVI_min和NDVI_max换算SR_min(10.05)/(1-0.05)约等于1.105SR_max(10.85)/(1-0.85)约等于12.333。当然最稳妥的方式是直接在SR影像上统计累计1%和99%分位数不依赖NDVI的换算。最后取平均FPAR (FPAR_NDVI FPAR_SR)/2FPAR的结果范围应该在0到0.95之间。如果出现负值或明显大于1的值多半是NDVI超出了统计范围导致线性拉伸越界把NDVI_min和NDVI_max重新统计一遍就好。4.3 第三步温度胁迫系数T_ε温度胁迫系数由两个子系数相乘得到T_ε T_ε1 × T_ε2。T_ε1衡量的是高温或低温对光合作用的抑制公式是T_ε1 0.8 0.02 × T_opt - 0.0005 × T_opt²同时有两个限制条件当月均温T_mean≤-10℃时T_ε1直接取0因为极端低温下光合作用基本停止当月均温高于最适温度10℃以上时T_ε1也取0这是高温胁迫导致光合效率趋近于零的经验阈值。在ENVI的Band Math里T_ε1可以写成T_e1 (0.80.02T_opt-0.0005T_opt*T_opt) * (T_mean gt -10) * (T_mean le T_opt10)其中(T_mean gt -10)在ENVI里是一种布尔比较运算满足条件返回1不满足返回0直接作为掩膜乘法使用这样就把分段条件嵌进了表达式不需要额外做掩膜。T_ε2描述的是温度偏离最适温度时的递减效应公式稍长T_ε2 1.1814 / ((1exp(0.2*(T_opt-10-T_mean))) * (1exp(0.3*(-T_opt-10T_mean))))在Band Math里T_opt和T_mean都用栅格变量exp函数直接写exp。这个计算本身不复杂但括号嵌套多写表达式时建议反复检查。算完先看统计值正常情况下T_ε2应该在0到1之间。T_ε就是两个结果相乘分开算更方便排查如果你处理的是月尺度数据每个月对应不同的T_mean影像每个月都要算一组T_ε。4.4 第四步水分胁迫系数W_ε改进CASA模型里水分胁迫系数用LSWI来计算。先用Band Math算LSWI(float(b5)-float(b6))/(float(b5)float(b6))这里b6对应Landsat 8的短波红外1波段。算完后统计LSWI影像的累计95%分位数作为LSWI_max。然后计算W_ε(1LSWI)/(1LSWI_max)在ENVI里把统计出来的LSWI_max作为数值常量直接代入表达式即可省事。如果研究区较大、不同区域水分背景差异明显也可以用生长季最大LSWI合成影像作为变量但那样对年际间气候波动比较敏感。我自己用下来的体会是取累计95%分位数当常量最稳。如果个别湿润像元算出来W_ε略大于1可以在统计或最终合成前用Where函数把大于1的像元压回1但如果是大范围超过1说明LSWI_max取得偏低应改用累计99%分位数。4.5 第五步APAR、光能利用率ε与NPP合成到这一步所有中间变量都齐了按公式逐级合成。先算APARAPAR SOL × FPAR × 0.5SOL是月太阳总辐射栅格单位MJ/m²·月要确保和FPAR是同一月份或同一生长季尺度。再算光能利用率ε ε_max × T_ε × W_εε_max按研究区主导植被类型选定比如以农田为主的研究区可取0.542以阔叶林为主取0.485。最后算NPPNPP APAR × ε得到的NPP单位是gC/m²·月。如果想得到生长季总量把各月NPP相加如果需要换成吨碳每公顷tC/ha再做单位换算。合成完NPP后用研究区矢量边界做一次掩膜把边界外像元去掉再检查最大值和最小值。正常情况下县域尺度生长季NPP月值在几十到两百多gC/m²·月之间如果出现几千甚至几万的异常值先检查掩膜把云、阴影、水体混进来的像元剔除干净。5. 实操中的参数取值与问题排查5.1 常用参数取值速查表写这篇总结的一个重要动力就是把自己踩过的参数坑记下来。下面这张表是我这次流程里用到的核心参数可以直接抄参数取值/方法说明NDVI_max影像累计99%分位数约0.78~0.90按实际统计结果NDVI_min影像累计1%分位数约0.02~0.08按实际统计结果SR_max(1NDVI_max)/(1-NDVI_max)由NDVI范围换算SR_min(1NDVI_min)/(1-NDVI_min)由NDVI范围换算FPAR_max0.95文献常用FPAR_min0.001文献常用ε_max农田0.542、草地0.542、阔叶林0.485、针叶林0.389、灌丛0.429按植被类型选T_opt生长季NDVI最高月平均气温区域常量或分区插值LSWI_maxLSWI累计95%分位数按实际统计结果APAR系数0.5光合有效辐射占太阳总辐射比例这里再强调一遍ε_max的取值对最终NPP结果影响非常大。如果有条件最好用研究区实测NPP数据反推校准没有实测数据至少要在论文里写清取值来源和依据不能闷头拍一个数字。5.2 常见问题与排查技巧整理几个实际操作里最高频的问题每一个我都遇到过。第一个问题NDVI结果整体偏低。大概率是大气校正没做好的锅。Landsat原始DN值直接算NDVI大气程辐射会让红光波段被散射增强、近红外被吸收削弱NDVI系统性偏低。解决方案只能回到FLAASH重新检查大气模型、气溶胶模型和能见度参数。另一个常见原因是影像本身云量多或者研究区有雾霾这种情况只能换时相。第二个问题NPP结果出现大量负值。负值往往不是来自NPP公式本身而是APAR里的SOL辐射数据有问题。如果气象站太阳总辐射插值后某些像元出现负值或异常小的值乘到APAR里就会出负值。排查时先单独查看SOL栅格统计再依次查看FPAR、T_ε、W_ε中间变量的取值范围用排除法定位。第三个问题Band Math报错“Variable has no value”或者维度不匹配。这种问题九成是输入影像的行列数或投影不一致。解决方法是先对所有输入栅格统一执行Resample Data或Reproject确保投影、像元大小、行列数完全一致。最好在做任何一个波段运算前先查看所有文件的基本信息右键View Metadata确保数值一致再开工。第四个问题FLAASH大气校正后影像变暗或者出现异常条纹。多半是辐射定标单位没转对。FLAASH要求的辐亮度单位是μW/(cm²·nm·sr)但ENVI辐射定标默认输出可能是W/(m²·μm·sr)两者数值大约差10倍。这种量级差异足以让大气校正结果面目全非。解决办法是定标后先查看辐亮度影像统计值如果最大值在0到20之间通常接近μW单位如果最大值在0到200之间通常是W单位然后按需求做换算。第五个问题NPP空间分布和土地利用类型明显不符。比如农田区域的NPP反而不如建设用地高这种违背常理的情况先检查掩膜是否把建设用地像元剔除干净。其次是研究区如果混合了多种植被类型而全局只用了同一个ε_max结果会偏。这种情况应引入植被类型图分类型用不同ε_max计算最后再按类型合成。关于结果验证我常用的做法是拿MODIS的MOD17A3 NPP产品做对照。MOD17A3空间分辨率500米时间尺度是年它本身也是基于光能利用率模型算出来的虽然分辨率不一样但数量级和空间分布可以作为交叉验证的参考。把MOD17A3重采样到30米统计研究区内NPP总量和平均值跟自己的结果对比如果误差在20%以内说明整个流程基本可信如果差了好几倍优先检查ε_max和气象辐射数据。这种验证方式虽然不是绝对真值检验但能很快暴露出流程里明显的系统性错误。最后分享一个小技巧每次算完一个中间变量我都顺手在ENVI里做一次直方图统计把最小值、最大值、均值、标准差记下来。CASA链路是按顺序相乘的任何一个中间变量取值异常最终NPP都会有明显病征。把所有中间结果的统计值保存下来排查时一路往回查比对着结果图瞎猜高效得多。这个系列写到这里CASA反演NPP的完整流程也算交代清楚了。从原理上讲改进CASA模型并没有多高深本质就是把经典模型里不适合区域尺度应用的参数换成遥感可观测的替代量从操作上讲真正花时间的也从来不是公式本身而是数据预处理和参数标定这些“磨刀”的活。希望这篇总结能让你少走一点我走过的弯路。