基于数字影像的北京松山森林物候模拟及其与气象因子的关系*

李润东 田文东 于海群 李鑫豪 靳 川 刘 鹏 查天山 田 赟

(1.北京林业大学水土保持学院 北京 100083；2.上海勘测设计研究院有限公司 上海 200335；3.北京林业大学水土保持国家林业和草原局重点实验室 北京 100083；4.北京市房山区林果科技服务中心 北京 102400；5.北京市园林绿化规划和资源监测中心(北京市林业碳汇与国际合作事务中心) 北京 100013)

植被物候一般指植被生长发育过程中，如叶片展开、开花和叶片衰老等周期性事件的发生时间和持续时间(Zhouetal.，2013)，不同空间尺度上的植被物候可观测提供有关生物圈与大气之间能量交换、碳水通量动态变化的重要信息(张学霞等，2003)。植被物候变化是通过光、温度、土壤湿度和养分有效性等的微妙变化来协调的，对生态系统生产力具有控制作用，高精度模拟气候变化条件下的植被物候变化至关重要(代武君等，2020；王连喜等，2010)。目前，在不同时空尺度上，植被物候与气候变化之间的关系仍然存在不确定性，传统植被物候监测依赖于人类对离散物候事件或物候期的直接观察，虽然比较客观准确，但费时费力，成本较高(高琪等，2019)；遥感观测具有连续观测、覆盖度高等优点，但获得的物候指标值与地面实际观测值存在偏差，同时空间分辨率也较低。随着图像处理技术逐步成熟，数字相机影像作为物候监测的新方法，弥补了人工观测和遥感观测的不足，通过连续拍摄生态系统高时空分辨率影像，已经应用于全球多个生态系统中(Campilloetal.，2008；Browningetal.，2017)。

与此同时，涡度相关技术可以直接测定生态系统-大气间的碳交换，具有准确、连续、非破坏性等优点，监测数据可用于检验各种基于碳过程的模型，大力推动了陆地生态系统碳循环研究的发展(Hinko-Najeraetal.，2016；纪小芳等，2019；徐丽君等，2011)。近年来，植被物候与生态系统碳循环之间的关系成为生态学家们关注的热点，物候与碳循环的季节变化密切相关，被认为是气候变化的一项重要生物学指标，人们逐渐利用数字相机技术进行物候监测和生态系统碳循环研究(Keenanetal.，2014；Richardsonetal.，2013；周磊等，2012a)。以往研究表明，数字重复摄影能够量化落叶阔叶林、草地、农作物的光合作用期持续时间和生态系统总初级生产力(gross primary productivity，GPP)(Zhangetal.，2019)，数字相机衍生的相对绿度指数(relative green chromatic coordinate，Gcc)数据与碳通量数据高度相关(Kangetal.，2016)。对于落叶林，Gcc对生长季开始(start of growthing season，SOS)时间的估计值与野外人工观测和卫星衍生SOS的估计值非常匹配，与森林生态系统的卫星衍生数据相比，数字相机绿度指数衍生的物候数据与野外人工观测更相关(Kurcetal.，2010)。通过数字相机获得的绿度指数，能够准确描述森林物候年际变化，并与GPP有较强相关性，这说明数字相机不仅可提供在生态系统尺度上较为准确的植被物候信息，而且还可用来测算碳通量精度，同时也为研究植被物候的环境控制因子提供了新方法(Ahrendsetal.，2009；Richardsonetal.，2009)。有学者将数字相机数据与气象数据相结合研究植被物候与气象因子之间的关系，结果发现植被物候与环境因子关系密切，温度、辐射等均对植被物候具有重要影响，降雨对植被物候具有“触发”作用(周惠慧等，2016；周磊等，2012b)。本研究以北京松山落叶阔叶林生态系统为例，基于2019年数字相机时间序列数据，提取森林植被生长期中的关键指标，并将其与通量观测数据拟合估算出的物候指标进行对比分析，以检验物候模拟精度；同时结合观测站点环境因子数据，分析林地物候与环境因子间的关系，比较不同环境因子对植被物候的影响程度，以期提高植被固碳模型和区域碳固定模拟的准确性。

1 研究区概况与研究方法

1.1 研究区概况

研究地点位于北京松山国家级自然保护区内(115°47′11″E，40°30′48″N)，占地面积4 671 hm2，海拔1 165 m，距北京市区约104 km。属大陆性季风气候，是暖温带与中温带、半干旱与半湿润的过渡地带，年均气温8 ℃，年降水量424.6 mm，主要集中在7、8月，占年降水量的62%。年日照时数2 726 h，日照率62%，植物生长期160天左右，生长季叶面积指数3.65。土壤类型为棕壤，pH 6.34，土壤有机质含量149.75 g·kg-1。林分平均冠层高度4 m，结构为复层结构，乔木层包括核桃楸(Juglansmandshurica)、大果榆(Ulmusmacrocarpa)、华北五角枫(Acertruncatum)、大叶白蜡(Fraxinusrhynchophylla)、蒙古栎(Quercusmongolica)等；灌木层包括绣线菊(Spiraeasalicifolia)、铁线莲(Clematisflorida)等；草本层包括异穗薹草(Carexheterostachya)、等齿委陵菜(Potentillasimulatrix)、牛扁(Aconitumbarbatumvar.puberulum)等(李润东等，2020)。

1.2 数据获取

本研究主要数据来源于近地表相机记录的Gcc、涡度相关技术观测的通量数据、归一化植被指数(normalized difference vegetation index,NDVI)数据以及微气象观测系统观测的空气温度(air temperature，Ta)、土壤温度(soil temperature，Ts)、光合有效辐射(photosynthetically active radiation，PAR)、饱和水气压差(vapor pressure deficit，VPD)和降水量(precipitation，P)。

数字相机安装在观测塔上高度20 m处，图像采集于2019年，从当地时间早上9点到下午4点每30 min采集1次图像，并自动存储在数据采集器(CR-1000X)中，格式为JPG。松山观测站拍摄图像样例如图1所示。Gcc数据依据所采集图像提取，计算过程详见数据处理部分，最终获取Gcc数据样本365个。NDVI数据由SRS光谱反射传感器(安装在观测塔上高度6 m处)测得的入射光和反射光计算得到。涡度相关通量数据采用2019年GPP数据，由SmartFlux在线计算模块计算每30 min数据的平均值。微气象观测数据包括空气湿度和温度、降雨、土壤含水量(soil volume water content，SWC)、土壤湿度和温度、光合有效辐射、风向、风速等。各种气象因子传感器每分钟采集1次数据，由采集器(CR-1000X)记录，各种数据每周下载。最终计算空气温度、土壤含水量、光合有效辐射等日均值以及降雨量日累计值，得到各环境因子数据样本365个。

图1 北京松山落叶阔叶林样品摄影(a)及森林生长动态(b)Fig.1 Sample photography(a)and broad-leaved forest growing dynamic(b)in Songshan,Beijing红色多边形表示感兴趣区域，用于计算相对绿度指数。Red polygon represents the ROI (region of interesting)used to calculate relative green chromatic coordinate (Gcc).

1.3 数据处理

1.3.1 相对绿度指数(Gcc)计算 首先筛选采集的图像，剔除因降雨等原因造成的质量较差的图像，采用相邻日期图像进行插补。构建数字图像时间序列后，选取一张质量较高的图像作为模板，在模板图像上定义合适的感兴趣区域(region of interesting，ROI)。利用ROI对所有图像进行掩蔽，计算每个像素的RGB颜色通道信息(Zhou，2019)。本研究计算白天9:00 AM—4:00 PM每张图像内ROI的Gcc，公式如下：

(1)

式中：Gdn、Rdn和Bdn分别为ROI中像素点R、G、B光谱波段的平均值。

为了降低高频Gcc数据噪声并检测季节周期，采用3天移动窗口法对Gcc时间序列进行平滑处理，称为第90百分位法(Richardsonetal.，2018)。对每30 min的Gcc结果取平均值，得到每日的值。

1.3.2 时间序列数据处理 日时间序列环境因子数据包括Ta、Ts、SWC和PAR等，均采用9天窗宽的移动平均法进行过滤(Zhou，2019)。为了确定Gcc与环境因子间的时间序列关系，Gcc的每日值也采用相同移动平均法进行平滑处理。通过绘制生长曲线和检验Pearson相关系数分析Gcc与环境因子间的相关关系。

1.3.3 时间序列拟合与物候指标计算 基于平滑插值的时间序列数据，采用Klosterman方法对NDVI时间序列数据、Gcc数据以及GPP数据进行平滑和插值，以拟合森林生长季节轨迹(Klostermanetal.，2014)。地面通量实测数据定义当生长初期第1次出现连续3天总GPP小于夏季最大碳固定的5%时，为生长季的开始日期；而当生长末期第1次出现连续3天GPP小于夏季最大碳固定的5%时，为生长季的结束日期(Zhaetal.，2009；Xieetal.，2016)。从GPP中推导出的样本生长曲线和提取的主要物候指标如图2。

图2 从GPP中推导出的样本生长曲线和提取的主要物候指标Fig.2 A sample growth curve and key phenological indicators derived from GPPSOS：生长季开始Start of growing season；POP：峰值位置Position of peak value (maximum)；MSP：春季平均值Mean spring value；MGS：生长季平均值Mean growing season value；RSP：春季生长速率Rate of spring green up；RAU：秋季衰老速率Rate of autumn senescence；EOS：生长季结束End of growing season；MAU：秋季平均值Mean autumn value；LOS：生长季长度Length of growing season ；PEAK：峰值Peak value (maximum).下同The same below.

1.3.4 数字相机数据与地面通量数据和NDVI数据的差异 本研究认为，实测GPP时间序列数据提取的物候指标为真实物候指标，对Gcc、NDVI和GPP时间序列数据采用Klosterman方法拟合(周玉科，2018)，并计算各自生长季物候指标，进而计算所提取物候指标的相对差异，用于评价基于数字相机和SRS-NDVI数据的物候指标精度，公式如下：

(2)

式中：P_VIi表示Gcc与NDVI时间序列数据提取的第i个物候指标的值；P_GPPi表示GPP时间序列数据提取的第i个物候指标的值。

1.3.5 统计分析 植被指数能够反映植被生长状况，而植被生长发育受多种环境因子影响，这些环境因子对植被生长的影响程度需要准确估计。采用Pearson相关分析方法判定Gcc与各环境因子间的关系，看变量间有无一定程度的相关性。将松山观测站环境因子数据利用主成分分析法(principal component analysis，PCA)进行分析。主成分分析是一种多元统计技术，其采用正交变换法将一组相关变量转换成一组正交的、不相关的轴，具有降低数据空间维数的作用，是识别高维变量线性组合最常用的技术(Yaoetal.，2012)。为了定量分析Gcc与各环境因子间的关系，采用全子集回归分析法筛选Gcc与环境因子间的变量方程，利用贝叶斯信息准则(Bayesian information criterion)和调整决定系数(adjust coefficient of determination)对全子集筛选效果进行评价，最终建立Gcc与环境因子间的回归方程(谭丞轩等，2020)。

2 结果与分析

2.1 相对绿度指数与环境因子的关系

Gcc的变化趋势与Ta、Ts的变化趋势相似，Pearson相关系数分别为0.88和0.86(P<0.01)(图3a、b，表1)。Gcc与SWC相关性较低(图3c，表1)(R=0.29,P<0.01)，可能是降雨突发性影响土壤水分变化所致。对比降雨量(图3f)与土壤水分发现，土壤水分对降雨量有正响应关系，如7—8月时SWC随降雨增加而增加。PAR的季节变化也表现出与Gcc的一致性(图3d，表1)，Pearson相关系数R=0.45(P<0.01)。VPD与Gcc的季节变化趋势一致(图3e，表1)(R=0.65，P<0.01)。松山森林生态系统的Gcc与降雨虽然相关性较低(图3f，表1)(R=0.25，P<0.01)，但是降雨对森林变绿期具有“触发”作用，即降雨能够使土壤中微生物呼吸加强，促使土壤有机物分解，释放出植被生长发育所需养分，进而促进植被生长。Gcc随降雨变化呈现波动状态，这表明基于数字相机的Gcc具有可实时捕捉到森林生长动态变化的能力。

表1 相对绿度指数与环境因子的相关性①Tab.1 Correlation relationships between Gcc and environmental factors

图3 相对绿度指数与环境因子的关系Fig.3 Relationships between Gcc and environmental factors

2.2 相对绿度指数与环境因子的回归分析

从上述得知，Ta、Ts、PAR、P、SWC等均对植被生长发育存在不同程度影响，但任何一个环境因子都不能单独对植被生长起决定性作用，各环境因子信息相互叠加，存在相关关系，作为一个整体影响植被生长发育。利用主成分分析法(PCA)分析松山观测站环境因子之间的相关性(各变量数据样本365个)，根据特征值大于1的原则选取2个主成分，2个主成分对总变异的解释贡献率分别为57.03%和19.81%，累计贡献率为76.83%。各环境因子在2个主成分中的荷载如表2所示，荷载越高表示贡献率越高。由表2可知，在第1主成分中各指标载荷均为正值，说明环境因子与第1主成分具有一定相关性。土壤体积含水量在第2主成分中的载荷最大且为正值，说明第2主成分主要代表土壤体积含水量。

表2 环境因子在2个主成分中的荷载Tab.2 Environmental factors’load in the two principle components

采用全子集回归分析法，选择Ta、Ts、PAR、P、SWC为自变量，Gcc为因变量进行分析(各变量数据样本365个)。筛选不同数目自变量的随机组合方式，根据贝叶斯信息准则和调整决定系数筛选出最优组合，结果如表3所示，当组合变量包括Ta、Ts、VPD、PAR时调整决定系数最大，此时方程为最优模型，最终得到Gcc与环境因子间的回归方程如下：

表3 全子集筛选最佳组合统计结果Tab.3 Best combination result after total subset selection

Gcc=0.34+(1.68Ta-0.39Ts+

0.015PAR-9.89VPD)×0.001。

(3)

2.3 Gcc、GPP和NDVI数据提取的物候指标对比

物候学是全球生态学和陆地生态系统碳循环研究的新线索，在环境监测和管理中具有重要应用价值(Berraetal.，2021)。本研究比较基于数字相机时间序列数据提取的Gcc数据、地面实测的通量数据和NDVI数据的关键物候事件时间序列，关键物候指标不仅包括生长季开始和结束时间，还包括春季生长速率和秋季衰老速率。采用Klosterman方法对Gcc时间序列数据、NDVI时间序列数据和GPP数据进行拟合并计算物候指标，结果如图4所示，可以看出，各时间序列数据提取的物候指标差异明显。通过物候模型拟合提取物候指标的相对差异对比结果见表4。Gcc与GPP提取物候指标对比，除生长季结束时间和秋季衰老速率的相对差异小于0.1外，其余物候指标的相对差异均大于0.1，说明基于Gcc提取的物候指标与GPP差异较大，其中相对差异最大的是春季生长速率，为12.5。这是因为春季植被生长初期光合作用已经进行，但叶子较小，导致数字相机拍摄图像绿色波段值较低。随着植被快速生长发育，绿度值快速增长，基于Gcc提取的生长季开始时间较GPP晚，春季生长速率较大，但是秋季物候指标提取结果十分相近，说明数字相机对植被叶子变色过程捕捉敏感，Gcc提取的生长末期物候指标精度大于生长初期。NDVI与GPP提取物候指标对比，生长季开始时间相差6天，相对差异较小，说明SRS-NDVI 测量仪测定的NDVI数据能够敏感捕捉到植被生长状态变化。生长季结束时间与GPP的相对差异为0.021，说明提取的生长末期物候指标精度大于生长季初期，与Gcc提取物候指标结果一致。

表4 Gcc、NDVI和GPP提取物候指标的相对差异Tab.4 Relative difference between key phenological indicators derived from Gcc,NDVI and GPP

图4 利用Gcc(a)、GPP(b)和NDVI(c)反演森林生长曲线和主要物候指标Fig.4 Gcc (a),GPP (b)and NDVI (c)were used to invert forest growth curve and key phenological indicators

3 讨论

3.1 相对绿度指数与环境因子的关系

本研究中，Gcc与温度、光合有效辐射等环境因子间存在较强相关性，一些研究结果也表明，植被生长对温度、光合有效辐射等环境因子具有正响应关系(刘鑫等，2019；Seyednasrollahetal.，2020)，降雨对Gcc具有“触发”作用(周磊等，2012b；Berraetal.，2019)。在逐步回归分析中，温度是Gcc的主要解释变量，这可能是因为温度升高会使生长季开始时间提前(Huangetal.，2018)。SWC与植被生长间存在复杂的相互作用，随降雨发生出现明显波动，说明降雨可以有效补充土壤水分。土壤温度对提高植被生产力具有积极影响(徐满厚等，2013)，与本研究逐步回归结果一致，说明土壤温度在控制植被绿度方面起着关键作用。饱和水气压差可表征空气干燥程度，VPD与Gcc的Pearson相关系数为0.65，说明VPD对Gcc变化具有重要作用。在过去半个世纪中，温带落叶林物种对气候变化表现出广泛的物候响应(Elmoreetal.，2012；Chenetal.，2020)，但未来温度变化如何影响关键的物候期事件，如植被叶子出现和衰老的时间，还存在相当大的不确定性。

3.2 Gcc、NDVI和GPP数据提取物候指标的差异

Gcc与GPP的季节变异模式相似(图4)，即春季升高，初夏达到峰值，随后略有下降，秋季急剧下降。Gcc时间序列数据存在不确定性，如光照变化、被分析的ROI有其他物体造成背景效应、图像分辨率、气象扰动等，这些不确定性影响植被SOS日期(Ahrendsetal.，2008)。通过Gcc数据确定的其他物候指标日期均与GPP提取的物候指标日期相差不大，说明Gcc能够敏感捕捉到植被生长动态变化。生长初期与生长末期相比，数字相机对生长末期的物候指标提取比生长初期更加准确。之前研究中，通常利用遥感数据提取的植被指数估算生态系统物候指标日期，虽然遥感数据存在时间间隔长、云层遮挡等不确定性，但基于遥感数据提取的植被指数仍是物候研究中使用最广泛的指标(Duetal.，2019；Wangetal.，2016)。本研究中，NDVI同样具有类似的季节变化，与GPP提取的物候指标相比差异较小，生长季开始与结束日期的相对差异均小于0.1，说明NDVI数据能够有效捕捉植被生长动态变化，与以往利用遥感数据提取的植被指数监测结果一致(Melaasetal.，2016)。Gcc与NDVI相比，季节变化趋势相同，均为先升高，达到峰值后略有降低，随后急剧下降。卫星植被指数与数字相机提取的物候指标比较发现，Gcc与卫星植被指数具有显著相关性，Gcc能够敏感捕捉生态系统植被生长变化(Wangetal.，2017)，说明SRS-NDVI测量仪测定的NDVI数据用于物候指标估算的可行性较高。

3.3 局限与展望

本研究仅针对一个生态观测站，且时间序列数据较短，后续应继续构建补充数据集。物候学研究需要跨植物物种和不同时空尺度的植被状态长期观测(Sonnentagetal.，2012)，随着物候观测网络的进一步发展，更多生态观测站将配备物候相机，因此，更多的环境变量数据和更多的全球观测点物候图像会显著改善物候变化和生态系统碳循环研究(Whartonetal.，2016)。未来将收集更多站点的环境变量数据和物候图像，并添加遥感图像(如Phenocams、Modis、Landsat 8和Sentinel 2)的长期观测数据，以评估自然因素对不同植被类型、不同生态系统物候转换日期的影响，提高物候模拟精度。

4 结论

本研究以北京松山落叶阔叶林生态系统为例，利用数字相机拍摄图像提取物候信息，结合当地微气象观测系统的环境因子数据，分析二者间的相关性；将数字图像提取的Gcc、NDVI与GPP数据进行对比，分析其差异性，得到以下结论：1)在松山落叶阔叶林生态系统中，基于数字相机提取的Gcc与Ta、Ts、PAR、VPD等环境因子具有很强相关性，温度是影响植被生长变化的关键因子，降雨对Gcc具有“触发”作用；2)与GPP数据提取的物候指标相比，基于Gcc数据能够敏感捕捉到植被生长动态变化，且生长末期物候指标捕捉精确度高于生长初期；NDVI数据也可以敏感捕捉植被生长动态变化，与GPP提取物候指标的相对差异较小。