0 引言

在陆地生态系统中, 植被初级生产力包括植被总初级生产力(Gross Primary Productivity, GPP)和植被净初级生产力(Net Primary Productivity, NPP). GPP指绿色植物在单位时间单位面积上由光合作用产生的有机物质总量, NPP指GPP中扣除自养呼吸后的剩余部分^[1-2]. GPP/NPP作为陆地生态系统的关键参数^[3-4], 在预测未来大气中CO₂和其他温室气体含量、认识大气圈与生物圈的相互作用等科学问题研究上发挥着重要作用, 在全球碳循环过程中扮演着重要角色^[5-7].

随着全球气候变化和生态环境恶化, 自20世纪70年代起, NPP成为国际生态学和地学领域的研究热点之一, 国内外学者提出了多种NPP估算模型^[8-11].在众多NPP估算模型中, Carnegie Ames Stanford Approach(CASA)遥感过程模型是基于光能利用率和遥感数据来估算NPP的过程模型, 在全球和区域植被NPP的估算研究中有着十分广泛的应用^[12].例如, 董丹和倪健^[13]利用CASA遥感模型和空间分辨率为1 km的SPOT NDVI遥感数据估算了中国西南喀斯特地区1999—2003年的植被NNP, 指出了该地区植被NPP的时空变化特征; 张镱锂等^[14]利用CASA遥感模型和空间分辨率为8 km的GIMMS遥感数据研究了1982—2009年青藏高原高寒草地NPP的时空分异特征; 李秀芬等^[15]利用CASA遥感模型和空间分辨率为1km的MODIS遥感数据研究了2006年黑龙江省森林NPP的季节性变化规律.对于小区域尺度的NPP估算而言, 空间分辨率为1km的遥感数据因空间分辨率太低而不能发现NPP在区域细节上的变化; Landsat TM(Thematic Mapper)和OLI(Operational Land Imager)数据的空间分辨率为30 m, 有利于提高小区域尺度NPP估算的空间分布精度.如管小彬等^[16]基于Landsat影像监测了武汉市冬季NPP及其时空变化分析; 韩艳飞等^[17]分析了近30年关天经济区植被净初级生产力对土地利用变化的动态响应; 袁甲等^[18]分析了皖江城市带土地利用/覆被变化对区域净初级生产力的影响.

然而, 目前的研究主要利用单传感器反演的NPP来分析区域生态环境.由于MODIS数据和Landsat TM、OLI数据各自在空间分辨率和时间分辨率方面的不足, 导致基于MODIS数据生产的NPP产品难以适用于NPP空间分布规律的分析, 而基于Landsat数据生产的NPP则难以适用于NPP时间变化特征的分析.针对这一问题, 本研究提出了一种综合MODIS和Landsat的NPP产品进行小区域NPP时空变化特征分析的方法.门头沟区是北京市的生态涵养区之一, 具有重要的生态功能, 研究该区的NPP时空变化特征, 可为评价该区乃至北京的生态可持续发展系统提供一个重要生态指标.本文以门头沟区为研究区, 采用CASA遥感模型, 基于30 m空间分辨率的Landsat TM和OLI数据生产门头沟区2003—2014年每年中某一天的NPP产品, 用于分析门头沟区NPP的空间分布研究; 基于2003—2014年MODIS17A3产品(年NPP总量产品)和MODIS17A2产品(8天合成的GPP产品)开展门头沟区NPP的年际变化和GPP年内变化特征研究及驱动因素分析.

1 研究区与数据 1.1 研究区概况

门头沟区是北京市下辖的市辖区, 地处北京西部山区, 位于东经115°25′00″至116°10′07″, 北纬39°48′34″至40°10′37″之间.东西长约62 km, 南北宽约34 km, 总面积1455 km².门头沟区地处华北平原向蒙古高原过渡地带, 地势西北高, 东南低, 山地面积占门头沟区总面积的98.5%, 是北京市唯一的纯山区, 也是北京市西部重要的生态环境屏障, 具有十分重要的生态意义.

1.2 数据来源与预处理

本文选取2003—2014年每年一期无云的Landsat数据.数据选取规则是: ①由于8月份植被状况较快, 且有无云Landsat数据的年份较多, 首选每年8月份的数据; ②如果当年没有8月份的数据, 则选择距离8月份时间最短的数据; ③如果两期数据距离8月份时间相同, 则选择植被状况最好的数据.选取具体数据的日期见表 1.

Landsat TM和OLI影像来自https://earthexplorer.usgs.gov/, 预处理主要包括辐射定标, 大气校正、几何校正、投影转换(即将影像原始的UTM投影转成Albers等积投影)、数据拼接、裁剪等.

MODIS17A2/3产品来自https://earthexplorer.usgs.gov/, 数据预处理主要包括投影转换(即将影像原始投影转成Albers等积投影)、数据拼接、裁剪、降噪等.受云雨天气条件(或其他因素)的影响, 数据存在异常值.例如, 受云的阴影的遮挡, NDVI值会严重偏低, 进而影响GPP值严重偏低.对这些异常值删除后, 采用等距插值方法进行填充.

气象数据为斋堂气象观测站2003—2014年的观测数据, 包括年平均气温和年降水量, 以及对应Landsat影像获取日期的日平均气温.斋堂气象观测站位于门头沟山区, 观测数据几乎不受人为因素的干扰, 能够代表门头沟山区的气候特征.

2 方法和技术路线 2.1 GPP/NPP估算方法

基于Landsat数据和光能利用率原理建立NPP的模型如下公式所述^[19]:

$ \begin{align} {\rm NPP}=PAR\times FPAR\times e^\ast \times F(w)\times F(T)-R_{\alpha}. \end{align} $

(1)

其中: $PAR$ 为光合有效辐射, 单位为MJ/m²; $FPAR$ 为植被层对入射光合有效辐射的吸收分量, 无单位; $\varepsilon ^\ast $ 最大光能利用率, 单位为gC/MJ; $F(w)$ 为水分对光能利用率的影响因子, 无单位; $F(T)$ 为温度对光能利用率的影响因子, 无单位; $R_{\alpha}$ 是植被自养呼吸量, 单位为gC/m².

2.2 估算关键参数遥感反演

按照光能利用率模型, 根据门头沟区的实际辐射、气温、降水、植被指数、植被类型等实际数据反演了NPP模型的各输入参数.各参数的具体计算方法如下.

$PAR$ 是波长400 $\sim$ 700 nm范围内被大气吸收和散射等削弱后穿过大气到达地表的太阳辐射, 其值占入射到地表的太阳总辐射( $R_s $ , MJ/(m²·d))的50%左右, 所以 $PAR$ 可表示为

$ \begin{align} PAR=R_s \times 0.5, \end{align} $

(2)

$ \begin{align} R_s =R_a \times ({a_s +b_s \times n/N}), \end{align} $

(3)

$ \begin{align} N=24\times \omega _s /\pi. \end{align} $

(4)

式中, $n$ 表示实际日照时数, 可由气象站数据插值得到. $N$ 表示最大日照时数. $a_s $ 表示阴天( $n=0)$ 时大气外界辐射到达地面的分量, $b_s $ 表示晴天( $n=N)$ 时大气外界辐射到达地面的分量, 其中 $a_s $ 和 $b_s $ 的取值可以利用模型所推荐的经验参数, 即 $a_s =0.25$ , $b_s =0.5$ . $R_a $ 为大气上界太阳辐射(MJ/(m $^{2}\cdot $ d)), 可由以下公式计算得到:

$ \begin{align} R_a =\dfrac{24\times 60\times G_{sc} }{\pi } \times d_r \times ( {\omega _s \sin \phi \sin \delta +\cos \phi \cos \delta \sin \omega _s }), \end{align} $

(5)

$ \begin{align} d_r =1+0.033\times \cos ({2\pi J} /{365}), \end{align} $

(6)

$ \begin{align} \delta =0.409\times \sin ( {{2\pi J}/{365}}-1.39), \end{align} $

(7)

$ \begin{align} \omega _s =\arccos ( {-\tan \phi \tan \delta }). \end{align} $

(8)

式中, $G_{SC} $ 为太阳常数, 取0.082 MJ/(m $^{2}\cdot $ min), $d_r$ 为日地相对距离的倒数, $\delta $ 为太阳赤纬(弧度), $\omega _s $ 为日落时角(弧度), $\phi $ 为地理纬度(弧度), $J$ 为某一天在一年中的排序数.

$FAPAR$ 是基于Enhanced Vegetation Index(EVI)的线性方程来估算^[19]:

$ \begin{align} FAPAR={\rm EVI}=2.5\times \Big( {\dfrac{\rho _{nir} -\rho _{red} }{\rho _{nir} +6\times \rho _{red} -7.5\times \rho _{blue} +1}}\Big). \end{align} $

(9)

式中, $\rho _{nir} $ 是近红外波段的地表反射率, $\rho _{red} $ 是红波段的地表反射率, $\rho _{blue} $ 是蓝波段的地表反射率.

水分限制因子 $F(w)$ 可以表示为^[20]:

$ \begin{align} F(w)={( {1+LSWI})}/{({1+LSWI_{max }})}, \end{align} $

(10)

$\begin{align} LSWI={( {\rho _{nir} -\rho _{swir} })}/{( {\rho _{nir} + \rho _{swir} } )}. \end{align} $

(11)

式中, $LSWI$ 为陆表水分指数, $LSWI_{max}$ 是植被在生长季内的最大 $LSWI$ , $\rho _{nir}$ 、 $\rho _{swir}$ 分别为近红外和短波红外波段反射率.

本模型 $F(T)$ 按Raich等^[21]建立的陆地生态模型中的方法来计算:

$ \begin{align} F(T)=\dfrac{( {T-T_{min } })\times ( {T-T_{max } })} {( {T-T_{min } } )\times ( {T-T_{max } } )-( {T-T_{opt} } )^2}. \end{align} $

(12)

式中, $T $ 为大气温度, $T_{min}$ 、 $T_{opt}$ 和 $T_{max}$ 分别为光合作用的最低、最适和最高温度.门头沟区以C4植物为主, 对C4植物, 将 $T_{min}$ 和 $T_{max}$ 分别定义为0℃和50℃.假设植物已经适应了其所生长的温度环境的情况下, 最适温度可定义为生长季的长期平均温度.

最大光合利用率基于BGC模型进行模拟, 初始化为MOD17 NPP产品采用的取值.对于C4植物来说, 最大光合利用率 $\varepsilon^{\ast}$ 为一个常数2.76 g $\cdot $ (MJ)^-1.

自养呼吸量( $R_{\alpha})$ 可基于植物地上生物量(Above-ground Biomass, Bag, kg/ha)采用半经验公式计算^{[19, 22]}:

$ \begin{align} R_\alpha =\Big(0.53\times\dfrac{{\rm Bag }}{{\rm Bag} +50}\Big)\times {\rm e}^{(\frac{T_c -T_a }{50})}. \end{align} $

(13)

式中, $T_{c}$ 是植物光合作用的平均气温(℃), $T_a$ 是实测气温(℃), Bag为实测地上生物量 $({\rm g/m^2})$ .

3 结果分析 3.1 NPP时空分异 3.1.1 NPP的空间分布

根据上文的反演方法, 得到2003—2014年门头沟区NPP(30 m分辨率)空间分布图, 如图 1所示.受Landsat数据源限制, 每年尽量选择1景夏季(或接近夏季)的质量较好的影像进行NPP生产, 其中2011年和2012年Landsat数据缺失.对图 1(a)-(j)逐像元计算平均值, 得到近10年门头沟平均NPP的空间分布特征, 结果如图 1(k)所示.

分析图 1, 门头沟区NPP在空间分布上的空间分异特征显著, 西北部山区和北部及南部部分山区单位面积上的NPP值最高, 面积占研究区总面积的7.15%;东部平原区上的NPP值最低, 面积占研究区总面积的1.68%, 主要是因为平原地区建设用地多, 植被稀疏; 研究区中部为山涧河谷地区, 植被较为稀疏, 面积占研究区总面积的18.36%;其他区域面积占研究区总面积的72.81%.如图 1(k)中, 西北部山区和北部及南部部分山区NPP值在12 gC $\cdot$ m^-2· d^-1以上, 东部平原区NPP为4 gC $\cdot$ m^-2·d^-1以下; 中部山涧河谷地区NPP值在4 $\sim $ 8 gC $\cdot$ m^-2·d^-1之间; 其他区域NPP值分布在8 $\sim $ 12 gC $\cdot$ m^-2·d^-1.

总体上门头沟区NPP在空间上呈现四周高、中间低的分异特征, 受地形及人为因素影响, 极低值区主要沿河流谷地及东部河流出山口处分布, 呈现平放的"Y"型特征.

3.1.2 NPP的年际变化

MOD17A3产品为一年内NPP总量产品, 统计门头沟区2003—2014年每年的NPP均值, 结果如图 2所示.据图 2分析, 门头沟区2004年的NPP值为2003—2014年的最大值, 为308 gC $\cdot$ m^-2·a^-1.自2005年以来, 总体呈现平稳态势, 没有明显的增加或降低趋势.但年际变化较大, 2008年和2012年为两次峰值, 2010年为低谷.较2005年以来的均值184 gC $\cdot$ m^-2·a^-1而言, 2008年和2012年NPP的增幅分别为21%和25%, 2010年NPP的降幅为23%, 其他年份的波动幅度在10%以内.从周期的角度分析, NPP的峰值出现的周期为4年, 这一规律的内在原因有待结合其他相关数据进一步研究.

3.2 GPP年内变化

由3.1节分析得知, 门头沟2004、2008和2012年的年NPP高于历史平均水平, 2010年的年NPP为历史最低水平, 其余年份为正常年.随机选择2007年和2011年作为正常年份的代表年.受MODIS GPP数据处理难度的影响, 本文只分析这6年的年内GPP变化特征(其中2004、2008和2012年代表了GPP的高值年份, 2010年代表了GPP的低值年份, 2007和2011年代表了GPP的正常年份), 以期分析不同年份GPP的年内变化特征. GPP年内变化特征如图 3所示.

2004、2008和2012年GPP年内变化呈现正态分布, 其余3年的GPP年内变化呈现偏态分布.第1—121天和265—365天的时期内, 图 3中6个年份GPP的年内变化趋势基本一致, 第121—265天这一时期内GPP呈现较大的分异特征, 表现为极高、较高、正常、较低四种特征.

因此, 年际GPP的差异主要来源于年内第121—265天之间的GPP差异.

4 讨论与结论

本研究利用30 m Landsat TM和OLI影像反演NPP产品, 分析了门头沟区NPP的空间变化特征, 利用MODIS GPP时间序列产品分析了门头沟区的时间变化特征(年际变化和年内变化特征), 为评估该区域生态环境变化提供了数据基础.实验结果表明: ①门头沟区2003—2014年的NPP存在四周高、中间低的空间分异特征; ②自2005年以来NPP总体平稳, 没有明显的增加或减少趋势; ③ 5—9月GPP/NPP的差异对不同年份间GPP/NPP差异的贡献最大.

由于缺少高时空分辨率的NPP产品, 本研究利用30 m Landsat NPP日产品进行NPP的空间分布特征分析, 使用MODIS NPP年产品进行NPP年际变化分析, 各年Landsat数据时相的差异会对分析结果产生干扰, 采用多源遥感数据时空融合模型来获得同一天的多年Landsat NPP产品是未来解决该问题的方法之一^[23-24].地表环境参数遥感反演精度会影响NPP和GPP的反演精度, 如植被指数在叶面积指数大于3的地区存在饱和现象, 这将增加FPAR反演的不确定性.

本文方法存在以下不足: ①本文分析了NPP的时空变化特征, 但由于NPP受多重因素影响, 且各因素间相互影响机理复杂^[19], 本文没有进一步分析NPP时空变化的影响机理, 需要进一步加强相关方面的研究; ②由于Landsat数据的缺失，本文采用MODIS GPP/NPP产品进行了门头沟地区的GPP/NPP时间变化特征分析, MODIS的空间分辨率较低, 其大部分像元是混合像元, 这会降低MODIS GPP/NPP产品的精度.