基于CASA模型和SEVI指数的福建省植被NPP遥感估算与分析

江 洪，虞嘉玮，蒋世豪，黄贝莹，李玉洁

（1.福州大学数字中国研究院，福建 福州 350108；2.福州大学空间数据挖掘与信息共享教育部重点实验室，福建 福州 350108；3.福建省基础地理信息中心，福建 福州 350003）

植被净初级生产力（Net Primary Productivity，NPP）是指绿色植被在单位时间、单位面积上通过光合作用同化CO2并扣除自养呼吸消耗后的剩余有机物［1］，是评估陆地生态系统碳汇和调节过程的主要因子，在全球的碳平衡中发挥重要的作用［2］.目前，遥感反演植被NPP的模型主要有统计模型、过程模型和光能利用率模型.统计模型通常使用叶面积指数（Leaf Area Index，LAI）或者归一化植被指数（Normalized Differential Vegetation Index，NDVI）进行估算［3-4］.过程模型常见的有Biome-BGC［5］，BEPS［6］，TEM［7］，InTEC［8］等.光能利用率模型以太阳辐射为基础，结合植被的吸收、转化能力估算植被NPP，例如CASA［9］，SEBAL［10］，Nnpp［11］，C-Fix［12］，Glo-PEM［13］等模型.其中CASA模型综合考虑了太阳辐射、水分、温度、不同植被类型的影响.根据研究区的差异，通过参数的改进，CASA模型可以较好地应用于不同地区.例如朱文泉［14］等利用误差最小原则，根据中国的植被NPP实测数据模拟中国各种植被类型的最大光能利用率参数进而估算中国的植被NPP.在山地丘陵等地形崎岖的区域，遥感影像上光学辐射信号受地形影响会发生畸变，影响植被NPP估算精度.结合数字高程模型（Digital Elevation Model，DEM）数据对CASA模型中的参数进行改进可以达到一定的效果［15］，但是基于DEM的地形校正模型只对地形本影校正效果良好，对地形落影的校正效果欠佳［16］.近年提出的阴影消除植被指数（Shadow Eliminated Vegetation Index，SEVI）可以有效地消除地形本影和落影对山区植被信息的影响，并且无需DEM数据的支持［17］.因此，笔者将SEVI应用到CASA模型，改进其中的植被光合有效辐射吸收比率因子（Fraction of Photosynthetically Active Radiation，FPAR）计算模型，提高山区植被NPP遥感估算精度，并计算分析福建省2005年与2015年的植被NPP.

1 研究区与数据

1.1 研究区概况福建省（23°33′～28°20′N、115°50′～120°40′E）地处我国东南沿海，下辖9个地级市和1个综合实验区（如图1所示），陆地面积12.4×104km2，山地、丘陵占全省总面积的80%以上，地势总体上西北高东南低.福建省是我国南方重点林区之一，气候温和，冬短夏长，降水充沛，自然条件优越，植被种类多达5 000余种，森林覆盖率达66.8%［18］，居全国第一.

图1 研究区与验证样本示意图

1.2 数据源及预处理遥感影像采用2005年和2015年的MODIS数据，包括空间分辨率250 m的MOD13Q1数据和空间分辨率500 m的MCD12Q1数据和MOD17A3H数据（来源于美国国家航天局，https：//ladsweb.modaps.esodis.nasa.gov/）.MOD13Q1数据时间分辨率为16 d，用于CASA模型中光合有效辐射（APAR）的计算.MCD12Q1数据时间分辨率为1 a；用于光能利用率ε计算中植被生态系统的划分.MOD17A3H数据为NPP产品，时间分辨率为1 a，用于本文植被NPP估算结果的对比.气象数据采用福建省所有气象站点的月值观测数据，包括降水、气温、湿度、日照时数等（来源于中国气象数据网，http：//data.cma.cn/）.

MODIS数据的预处理包括数据格式转换、投影转换以及裁剪，并利用最大值合成法将年内多景MOD13Q1数据合成12个月的月值数据.DEM数据的处理包括格式转换、投影转换、镶嵌以及裁剪.气象数据的预处理包括：数据的清洗，用往年的平均值代替异常值或者缺失值；数据插值，采用克里金插值法将站点气温、降水、湿度等点数据转换成栅格数据.

2 研究方法

CASA模型主要通过光能利用率ε和植被吸收的光合有效辐射（Absorbed Photosynthetic Active Radiation，APAR）计算植被NPP，

其中，x表示像元，t表示时间；N PP(x，t)表示像元x在t月份的净初值生产力，单位：g·m-2；ε(x，t)表示像元x在t月份的实际光能利用率，单位：g·MJ-1；APA R(x，t)表示像元x在t月份获得的光合有效辐射，单位：MJ·m-2.

2.1 光能利用率的计算光能利用率指植被将实际吸收的太阳有效辐射转换成自身有机物的效率，通过温度胁迫因子T、水分胁迫因子W以及最大光能利用率εmax计算获得，见式（2）.温度胁迫因子的计算见式（3）和（4），水分胁迫因子的计算见式（5），最大光能利用率参数如表1所示［14］.

表1 最大光能利用率参数表 %

其中，ε(x，t)表示像元x在t月的光能利用率，Tε1(x，t)表示植被在低温与高温情况下受到的限制作用，由式（3）计算获得；Tε2(x，t)表示植被所处的环境温度，由公式（4）计算获得；Wε(x，t)表示水分条件对光能利用率的影响，由公式（5）计算获得；εmax表示在理想条件下植被转换成A P AR的效率.

其中，Tε1(x，t)表示像元x在t月的第一温度胁迫因子；T opt(x)表示像元x年内生长的最适温度，当像元x在t月的温度小于-10℃时，Tε1(x，t)=0［19］.

其中，Tε2(x，t)表示像元x在t月的第二温度胁迫因子；T(x，t)表示像元x在t月的温度，比最适温度高10℃或者低13℃时，Tε1(x，t)等于最适温度T o pt(x)时Tε2(x，t)的一半［19］.

其中，Wε(x，t)表示像元x在t月的水分胁迫因子；E ET(x，t)表示像元x在t月的实际蒸发量，采用区域实际蒸散模型求取［20］；P ET(x，t)表示像元x在t月的潜在蒸散量，根据Thornthwaite植被-气候关系模型求取［21］.

2.2 光合有效辐射的计算光合有效辐射由太阳辐射总量（SOL）以及植被光合有效辐射的吸收比率（FPAR）决定，

其中，SO L(x，t)表示像元x在t月的太阳辐射总量，单位MJ·m-2，由天文辐射总量和日照百分率线性拟合而成［22］，

其中，Q S为月太阳辐射总量；Q0为月天文辐射总量，起伏地形下的天文辐射总量通过DEM数据计算获得，见式（8）和（9）［23］；s为日照百分率；a，b为经验系数，引用陈惠［24］等利用实测数据拟合得到的福建省的经验系数，如表2所示.

表2 各月经验系数a和b

其中，Q是地形起伏下日天文辐射量，单位：MJ·m-2；E sc是太阳常数；r0/r是当天日地距离订正系数；n为可照时角的离散数目，δ为太阳赤纬角；ωsi和ωri分别是可照时间段的起始和终止太阳时角（°）.

其中，u、v、w是与地形、地理有关的特征因子，φ、α、β分别为当地的纬度、坡度和坡向.

2.3 FPAR改进植被光合有效辐射吸收比率（FPAR）是CASA模型估算植被NPP过程中主要的驱动因子.目前基于植被指数反演FPAR的方法众多［25］，常用的植被指数有NDVI和比值植被指数（RVI）等，由NDVI反演FPAR的计算公式

其中，FPA Rmax和F PA Rmin的取值分别为0.95和0.001，ND V I(i，max)和ND V I(i，min)的取值分别为第i种地类N DV I的95%及5%的百分位数［14］.

由于山区植被受地形阴影影响显著，根据文献［26］研究结果表明：未经地形校正的NDVI反演得到的FPAR在本影和落影的数值远低于相邻非阴影阳坡，相对阳坡的误差大于70%；经C模型或者FLAASH+C模型校正的NDVI反演得到的FPAR，本影的相对误差小于10%，但是落影的相对误差仍超过24%；通过SEVI反演的FPAR，本影、落影与阳坡的数值十分接近，相对误差小于3%.因此采用能有效消除地形阴影影响的SEVI替代NDVI，改进FPAR的反演模型

其中，SE VI(i，max)和SE V I(i，min)的取值分别为第i种地类SE VI的95%及5%的百分位数.

SEVI计算公式［17］

其中，ρnir和ρred分别代表近红外波段和红光波段的反射率，f(Δ)为地形调节因子，通过相关系数法进行优化选取.以0.001为步长，令f(Δ)从0到1.000，计算f(Δ)对应下的SEVI与RVI的相关系数r1，SEVI与SVI的相关系数r2，当r1和r2接近或相等时，确定最优的调节因子（f opt），计算公式

其中，x，y1，y2分别为遥感影像SEVI，RVI以及SVI的计算结果，n为SEVI，RVI以及SVI影像的像元数.

2.4 FPAR结果验证分别利用MODIS数据计算的NDVI和SEVI反演FPAR，选取相邻阴阳坡的FPAR样本，统计分析阴影区域和对应阳坡区域的FPAR均值以及相对误差，验证不同植被指数计算的FPAR精度.相对误差的计算公式

其中，FP A Rshadow为阴影处的植被光合有效辐射吸收比率，FPA Rsunny为相邻非阴影阳坡的植被光合有效辐射吸收比率.

250 m空间分辨率的MODIS影像相较于文献［26］中使用的30 m空间分辨率的Landsat8 OLI数据，难以区分本影和落影，但在地形复杂的山区仍有可以目视观察的阴影.根据地理学第一定律，一般情况下，地物间的距离越接近，其相关性也越大；同时，也有相关研究表明地处亚热带的福建省山地的阴坡水热条件充沛，阴坡与阳坡的天然林长势状态极为接近［27-28］.因此，在Google Earth影像参考下，通过目视解译，手动选取了45组相邻阴阳坡植被长势接近的样本，如图1所示.

3 结果与分析

3.1 FPAR估算结果估算结果显示采用NDVI反演的FPAR会造成阴影区域的数值偏低（如图2a所示），而采用SEVI反演的FPAR在阴阳坡的数值接近（如图2b所示）.样本FPAR统计结果显示采用NDVI反演的FPAR在阴影区域的相对误差为7.85%，而采用SEVI反演的FPAR在阴影区域的相对误差降至0.53%（如表3所示）.这说明采用MOD13Q1数据反演复杂地形山区地表参数时，有必要进行地形校正，这与文献［26］中基于Landsat8 OLI数据研究得到的结论一致.

图2 样本F P A R折线图

表3 FPAR在阴影、阳坡样本处的均值M以及相对误差E

3.2 植被NPP估算结果2005年福建省植被NPP估算结果平均值为861.9 g·m-2·a-1，不同植被类型NPP均值大小顺序依次为常绿阔叶林、针阔混交林、常绿针叶林、灌丛、草原和农用地，均值分别为1 082 g·m-2·a-1、908.2 g·m-2·a-1、863.7 g·m-2·a-1、852.7 g·m-2·a-1和711.2 g·m-2·a-1.2015年福建省植被NPP均值为855.7 g·m-2·a-1，不同植被类型年均NPP均值依次为1 107.9 g·m-2·a-1、875.5 g·m-2·a-1、831.7 g·m-2·a-1、796.7 g·m-2·a-1、764.8 g·m-2·a-1和758.0 g·m-2·a-1，大小顺序与2005年保持一致.

由于植被NPP实测结果较难获取，采用交叉验证法分析本文估算结果的合理性.查阅文献表明已有多位学者对2005年福建省植被NPP进行了估算（如表4所示），均值范围在788.3～811.5 g·m-2·a-1，常绿阔叶林、针阔混交林、常绿针叶林、灌丛、草原和农用地的NPP范围分别为710.7～1 600 g·m-2·a-1、639.0～1 100 g·m-2·a-1、1 022.2～1 300 g·m-2·a-1、800～900 g·m-2·a-1、700～800 g·m-2·a-1、600～700 g·m-2·a-1.基 于MOD17A3H数据估算的2001～2010中国东南部植被NPP均值为511.4 g·m-2·a-1，其中福建省的年均植被NPP范围在500～1 000 g·m-2·a-1，不同植被类型NPP大小顺序为常绿阔叶林、针阔混交林、常绿针叶林、草地、农田植被［29］.福建、广东、广西等长江以南各省区，常绿阔叶林实测NPP范围在910～1 340 g·m-2·a-1［30］.利用CASA模型估算的武夷山区一带亚热带森林年均NPP超过500 g·m-2·a-1［31］.BEPS模型估算的2004年福建省森林生态系统的NPP全年均值为578.97 g·m-2·a-1，阔叶林为780.0 g·m-2·a-1［32］.本文估算结果高于2004年BEPS模型估算的结果，主要是由于模型不同、年份不同，以及本文研究涵盖的植被类型更多，因此，结果更高.对比其他学者2005年的研究结果，除针阔混交林和灌木外，本文估算的各植被类型NPP结果与已有研究结果吻合.本文估算的福建省植被NPP结果高于Biome-BGC模型的估算结果［33］，也略高于其他学者的CASA模型估算结果［33-35］，可能的原因是SEVI相对NDVI提高了FPAR的估算值，进而提高了植被NPP的估算结果.

表4 2005年福建省植被NPP估算结果比较 g·m-2·a-1

3.3 植被NPP空间分布特征福建省年均植被NPP的空间分布总体呈现内陆向沿海递减的趋势，如图3所示.按照行政区划统计各地级市年均植被NPP，结果与图3相似（如表5所示）：2005年植被NPP最大的为龙岩市，其次为三明市、南平市、宁德市、漳州市、福州市、莆田市、泉州市和厦门市；2015年全省植被NPP第一、二位和第八、九位保持不变，第三到七位分别为漳州市、南平市、宁德市、莆田市和福州市.

图3 福建省2005年和2015年年均植被NPP空间分布

表5 不同行政区的年均植被NPP及变化

以100 m为步长，对福建省0～1 800 m海拔划分为18个等级，统计不同海拔等级植被NPP.结果显示，植被NPP随海拔变化的趋势在2005年和2015年基本一致.低海拔区域（500 m以下）植被NPP数值较低，随海拔升高NPP逐步增加；中海拔区域（500～1200 m）植被NPP达到峰值并稳定，在700 g·m-2·a-1附近小幅度波动；高海拔区域（1 200 m以上）植被NPP逐步下降，如图4所示.这一结果表明，福建省植被NPP的空间分布与海拔高度相关性强，沿海城市的平均海拔处于第一级，人类活动密集，森林覆盖率低，因此平均植被NPP较低；而内陆城市的平均海拔高度处于第二级到第四级，森林覆盖率高，人类对森林的干扰较少，因此平均植被NPP较高.

图4 不同海拔等级对应的植被NPP

3.4 植被NPP月度分布特征福建省月均植被NPP与温度、降水、日照百分率3个水热因子统计结果表明福建省植被NPP月平均值在时间序列上的分布呈现倒“V”型特征，如图5所示.2005年1～2月，全省植被NPP均值在20 g·m-2以下；从3～7月，植被NPP均值总体上升，7月达到最大值148.95 g·m-2；8～12月，全省植被NPP均值稳步下降，12月再次低至20 g·m-2·以下.2015年福建省植被NPP变化趋势与2005年一致，1～6月份的植被NPP均值高于2005年同月的结果，8月份全省植被NPP达到最大值123.84 g·m-2·，9～12月植被NPP均值逐渐减少.

图5 2005年和2015年月均植被NPP、温度、降水和日照百分率

2005年月平均温度1月最低，为8.01℃，7月出现最大值27.95℃，总体呈现先逐步上升，后缓慢下降的趋势.月均降水在6月出现最大值401.41 mm，12月出现最小值23.81 mm，上半年波动明显，下半年逐步下降.月均日照百分率也出现明显的波动特征，7月出现最大值58.65%，2月出现最小值11.82%.与2005年相比，2015年月平均温度走势保持不变，月均降水和日照百分率的走势区别较大，月均降水在8月出现最大值350.98 mm，1月出现最小值47.20 mm.日照百分率在10月出现最大值41%，12月出现最小值17.29%.综合月均植被NPP和水热因子的走势，可以看出植被NPP月度分布特征与温度因子的走势最为相似，相关性分析结果显示，二者的相关系数在2005年和2015年分别达到0.96和0.95.

3.5 2005～2015年植被NPP变化2005～2015年福建省植被NPP变化结果显示，福建省植被NPP总体呈现略有减少的趋势.结合图6和表5可以看到，植被NPP减少的城市主要为南平市、宁德市和三明市，其中宁德市减少27.45 g·m-2·a-1，减幅达3.54%.植被NPP增加的为福州市、莆田市、泉州市、厦门市、漳州市和龙岩市，其中漳州市增加117.95 g·m-2·a-1，增幅达15.78%.

图6 2005～2015年福建省植被NPP变化结果

4 讨论

采用SEVI指数改进CASA模型中FPAR参数的反演，降低了山区阴影和阳坡的FPAR相对误差.利用改进的CASA模型估算福建省2005年和2015年的植被NPP，估算结果略高于其他学者基于NDVI反演的结果，可能是由于SEVI降低了崎岖山区地形效应的影响从而提高了FPAR的估算结果，进而提高植被NPP的反演结果.

福建省植被NPP的空间分布特征在2005年和2015年保持了相同的趋势，龙岩市、南平市和三明市均高于800 g·m-2·a-1，而沿海地区的厦门市、泉州市和福州市均低于750 g·m-2·a-1，说明植被NPP的空间分布特征较为稳定.这种空间分布特征可能与森林覆盖率等因素有关，龙岩市、南平市和三明市森林覆盖率较高，同时山区人类活动较少，植被生长状态好.厦门市、泉州市和福州市地处福建省沿海，城镇化水平高，人口密度大，植被覆盖率较低.不同海拔对应的植被NPP差异明显，海拔较低时，人类活动明显，对植被的生长影响较大，而在中海拔区域，较少的人类活动以及合适的地理环境均有利于植被的生长，植被NPP达到最大值，而过高的海拔不适宜植被的生长.因此，植被NPP在海拔1 200 m以上出现明显的下降趋势.

福建省植被NPP在2005年和2015年表现出相同的月序特征，夏季（6～8月）出现最大值，冬季（12～2月）出现最小值，主要是由于夏季平均气温高、日照充足、降水少，适合植被的生长，而冬季的水热条件较差，植被生长缓慢，导致植被NPP数值较小.

2005～2015年福建省植被NPP总体呈现略有减少的趋势.根据中科院地理科学与资源研究所公布的中国土地利用现状遥感监测数据，统计了福建省2005年到2015年土地利用类型的变化，发现有林地、灌木林、疏林地的占比均下降，而城乡、工矿、居民用地占比均上升.由于林地的植被NPP高于建设用地的植被NPP，因此契合福建省植被NPP略有减少的趋势.此外，统计空间分辨率500 m的MODIS NPP产品（MOD17A3H）也表明2005年到2015年福建省行政区域内的植被NPP均值由925.91 g·m-2·a-1略降为922.59 g·m-2·a-1.还有其他相关研究也得出福建省近年来年均植被NPP呈现下降趋势的结论［29，36-37］.

5 小结

采用SEVI指数改进CASA模型中FPAR参数的计算，可以提高CASA模型对山区植被NPP的估算精度.福建省2005年和2015的植被NPP反演结果表明：2005～2015年福建省植被NPP出现小幅度的减少，但是各植被类型NPP以及时空分布特点保持不变，不同植被类型NPP均值排序依次为常绿阔叶林＞针阔混交林＞常绿针叶林＞灌丛＞草地＞农用地.

福建省植被NPP空间分布上呈现内陆高沿海低的特征，并且随海拔升高呈现先增加后减少的特点，在500～1 200 m的中海拔区域达到最大值.植被NPP月度分布特征为7月和8月达到最大值，而12月和1月出现最小值，与月平均温度的走势相关性强，2005年和2015年月均植被NPP和月均温度的相关系数分别为0.96和0.95.