基于遥感反演的中国森林冠层叶氮浓度 空间格局初步研究

于泉洲，刘煜杰，周 蕾，石 浩，孙雷刚

（1.聊城大学 环境与规划学院，山东 聊城 252059；2.中国科学院 地理科学与资源研究所 生态系统网络观测与模拟重点实验室，北京 100101；3.中国环境科学研究院生态文明研究中心，北京 100012；4.浙江师范大学 地理与环境学院，浙江 金华 321004；5.西北农林科技大学 资源环境学院，陕西 杨凌 712100；6.河北省科学院地理科学研究所，河北 石家庄 050011）

我国森林生态系统在固定大气CO2和减缓气候变暖方面具有巨大潜力[1-4]。森林生态系统氮素水平，尤其是冠层叶片氮浓度（Canopy Nitrogen Concentration,CNC,g/100g）制约着森林的固碳能力[5]。因此，获得准确的森林CNC 空间化数据便成为区域森林生产力和碳汇评估的关键。传统的野外采样方法获取CNC 耗时费力，并且无法获得空间连续的数据，而遥感反演方法则能较好地解决这一问题。当前，对于森林生物量、蓄积量和叶面积指数等参数的遥感反演已较为成熟[6-7]，但对于森林冠层生化参数的遥感反演则正在兴起。

近年来，国外学者基于HyMap、AVIRIS、CASI 以及Hyperion 等多种高光谱数据，在世界多地开展了研究[8-13]，内容包括森林冠层结构、生化组分及树种组成等方面。由于高光谱数据的空间覆盖有限，研究多被限制在小区域。直到Ollinger 等[10]发现北美森林CNC 与近红外反射率（NIR，800～850 nm）呈显著的线性正相关，并 结 合MODIS（Moderate-resolution imaging spectroradiometer）反射率数据，实现了北美森林CNC 的反演。虽然，Knyazikhin 等[14]对CNC 与冠层NIR 的相关关系提出了质疑，并且Wang 等[15-16]在对欧洲混交林CNC 的遥感反演研究中指出，冠层结构与叶化学组分之间的“协同效应”使得CNC 反演效果较好。Yu 等[17]对我国长白山混交林的研究发现，森林CNC 的空间变异主要受森林树种组成的影响，在跨植被功能型的森林生态系统中可以谨慎地使用卫星遥感获得的NIR 来反演森林CNC。因此，Ollinger 等[10]的方法仍被视为一种大尺度森林CNC 遥感估计的可行方法而被广泛应用[13,18]。

国内相关研究较少且集中在小区域。如沈艳[19]、袁金国和牛铮[20]在西双版纳地区开展了Hyperion 星地同步实验，研究了植被覆盖区冠层叶绿素含量、叶氮浓度、碳氮比等。杨曦光等[21]利用Hyperion 数据结合4-scale 模型和神经网络等机器学习模型反演了我国小兴安岭地区森林CNC。Yu 等[22]通过采样调查结合Hyperion 数据分析发现江西中亚热带人工针叶林单一的冠层结构导致NIR 与CNC 相关性降低，限制了区域CNC 的反演。

因此，本研究选取中国生态系统研究网络（CERN）中的千烟洲站和长白山站及其周边的林区为对象，通过森林样方调查采样、测定了438 份叶片氮浓度，并结合Hyperion 高光谱数据及MODIS 的BRDF 反射率数据，反演生长季中国区域森林CNC，同时对反演结果进行验证和空间特征分析。研究旨在揭示我国森林生长季CNC 的空间格局特征，并为国家尺度森林碳循环模拟提供更精确的初始化参数数据，具有重要实践意义。

1 数据与方法

1.1 研究区概况

选取吉林省长白山温带针阔混交林和江西省千烟洲中亚热带人工针叶林为采样区域。长白山温带针阔混交林位于吉林省长白山北坡，是我国保存最完整的自然森林生态系统之一[23]，长白山森林树种类型多样，主要建群种有红松Pinus koraiensis、椴树Tilia tuan、蒙古栎Quercus mongolica、水曲柳Fraxinus mandschurica和色木槭Acer mono等。千烟洲中亚热带人工针叶林位于江西省吉安市内，该地区经过30 多年的综合治理植树造林，森林景观已经形成且高度郁闭[24]，形成了以湿地松Pinus elliottii、马尾松Pinus massoniana杉木Cunninghamia lanceolata为主，兼有木荷Schima superba、枫香Liquidambar formosana等原生阔叶树的人工针叶林生态系统。

1.2 数据处理

1.2.1 遥感数据预处理

1）Hyperion 数据处理

EO-1 Hyperion 数据覆盖范围为7.5 km×180 km的东北西南向条带区域，空间分辨率30 m，光谱范围为355～2 577 nm，光谱分辨率约10 nm。千烟洲地区的Hyperion 影像成像时间为2009年5月初（影像编号为：EO1H1210412009121110PY_PF2_01）。长白山地区的Hyperion 数据获取时间为2012年6月底（影像编号为：EO1H1160302012179110P0_PF2_01）。由于森林为成熟林年际变化较小，因此星地数据可达到季相匹配的要求。选用的Hyperion 数据级别为L1Gst，该级别数据经过辐射校正、系统级几何校正和地形校正[25]，故仅对影像进行大气校正和几何精校正。

2）MODIS 数据

使用的MODIS 数据包括MCD12Q1 和MCD43B4。利用MRT（MODIS Reprojection Tool） 软件进行数据的投影转换和拼接，将正弦投影数据转为Albers 等面积投影，中国陆地区域共利用19 幅MODIS 数据拼接得到。

MCD12Q1 是500 m 分辨率的土地覆盖类型产品。本研究采用其中IGBP（International Geosphere Biosphere Programme，国际地圈生物圈计划）定义的分类系统，将其中的常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林和混交林5 种土地覆被类型合并作为本研究的森林CNC 反演的区域。

MCD43B4是BRDF模型校正的8d合成的地表反射率产品，其空间分辨率为1 km。因我国南方部分林区生长季云量较多，为避免数据大量缺失，故选用2011—2013年间生长季（DOY=177～241）的27个时期的数据进行合成。合成方法是将27 个时期数据按照像元位置逐个时期累加然后取平均。经统计，合成提取的森林地表NIR 像元面积占总林地面积的99.23%，说明3年生长季合成的地表NIR 能基本覆盖我国绝大部分森林区域，可以用于全国森林CNC 的反演。

1.2.2 采样与样品处理

千烟洲地区森林采样时间为2012年4月下旬，长白山地区采样时间为2013年7月上旬，共选取采集了59 个森林样方的438 份叶样品（表1）。

表1 野外采样的树种信息Table1 Tree species information in fieldworks

采样和样品处理步骤如下：

1）在Hyperion 覆盖区域随机选择森林郁闭度高，且地形较平坦的林分作为样地。

2）用样方绳在林分中心设置30 m×30 m 的样方，用GPS 在样方中心测量经纬度坐标。

3）进行样方调查。统计样方中高度达到乔木层的树种及棵数，测量样方中每个树种5 棵代表性立木的胸径。收集冠层上方的叶片样品，取鲜叶100 g 左右，装入封口袋并标记带回实验室。

4）样品在烘箱65 ℃下烘至恒重，利用球磨仪研磨成粉。用德国Elementar 公司的Vario EL cube 元素分析仪测定叶片氮浓度。

1.2.3 叶片氮浓度的尺度扩展

为使星地数据空间匹配，需将叶片尺度的氮浓度扩展到样方冠层尺度。对于纯林样地，以样方内同一树种的多棵树叶片混合测得的叶片氮浓度作为整个样方冠层的叶氮浓度；对于混交林样地，两地区采用不同的叶氮浓度尺度扩展方法。长白山地区采用样方内树种优势度加权的方法进行，即采用某树种的横截面积在整个样方内所有树种总横截面积的百分比作为加权系数[17]。千烟洲地区，采用冠层叶生物量加权的方法。通过收集的优势树种冠层的叶生物量方程（具体见文献[22]），计算样方不同树种的叶生物量，然后将每个树种的叶生物量除以样方总叶生物量，得到每个树种在该样方中的叶生物量比例（加权系数），然后进行不同树种的叶片氮浓度进行加权，最后得到样方CNC。

1.2.4 森林冠层叶氮浓度遥感反演与验证

全国生长季森林CNC 反演的具体流程如图1：

其中，构建的样方加权的森林CNC 与Hyperion 提取的NIR 之间的统计模型如下：

图1 基于MODIS BRDF 数据的中国区域森林CNC 反演流程Fig.1 The flow chart of forest CNC inversion based on MODIS BRDF data in China

由式1 和图2看出，在冠层尺度上CNC 与NIR 具有显著正相关关系（R2=0.75，P＜0.000 1）。

图2 基于森林样方数据的CNC 与NIR 的相关关系Fig.2 Relationship between forest CNC and Hyperion NIR from forest plots

收集了2013年生长季对我国6 个森林站点乔木叶氮浓度的调查数据和王晶苑等[26]对中国4 个森林站点的主要树种叶氮浓度调查数据，将每个站点内的各主要树种的叶氮浓度取平均值作为验证数据。验证站点信息见表2。

表2 CNC 反演验证的森林站点信息Table2 Forest site information for validation of CNC inversion

2 结果与分析

2.1 实测叶片氮浓度特征

调查的59 个森林样方共收集到438 份叶片样品。其中长白山样品中，主要优势种红松、紫椴、糠椴、蒙古栎、水曲柳和白桦的叶片氮浓度均值分别为1.43、2.80、2.84、2.45、2.98 和2.34 g/100 g，长白山所有样品叶片氮浓度均值为2.44 g/100 g。江西千烟洲调查了7 个树种共21 份叶样品，其中优势种湿地松、马尾松、杉木和木荷的叶片氮浓度均值分别为0.91、1.33、1.59 和1.42 g/100 g，千烟洲叶片样品氮均值为1.55 g/100 g。

2.2 中国森林CNC 空间特征

2.2.1 中国森林CNC 空间格局

反演的森林CNC 分布结果具有显著的空间差异特征，可以大致依照胡焕庸线进行划分。首先，在胡焕庸线左侧（西北部）森林分布面积明显少于右侧（东南部），同时左侧的森林CNC 也低于右侧（图3）。具体而言，结合中国植被图[27]对比分析，东北大小兴安岭、川西高原以及藏东南等针叶林占优势的区域CNC 相对较低；对于东北长白山脉、秦岭、云南南部及海南岛尖峰岭等阔叶林为主的区域CNC 较高；而在我国长江中下游、南岭、武夷山脉等阔叶原始林、阔叶次生林和人工针叶林混交分布的地区森林CNC 居中。

图3 反演的中国区域森林CNC 空间格局（1km2 分辨率）Fig.3 Spatial patterns of forest CNC in China (1 km2 resolution)

2.2.2 森林CNC 的分省统计特征

对全国32 个省、直辖市和自治区（将香港和澳门特别行政区的森林CNC 计入广东省）的森林CNC 进行区域统计，包括各省区森林面积、CNC最小值（min）、最大值（max）、均值（mean）、标准差（SD）和变异系数（CV）（表3）。分析发现，森林CNC 均值最高的三个省份分别为海南、吉林和辽宁，CNC 均值依次为2.66 g/100g、2.47 g/100g和2.44 g/100g。而森林CNC 均值最低的3 个省份分别是新疆、山东和天津，CNC 均值依次为 1.71 g/100g、1.68 g/100g 和1.66 g/100g。

将森林CNC 变异系数越高说明区域森林CNC变异性越强，CV 越低说明区域森林CNC 分布越均匀。CV 最大的3 个省份是天津、山东和青海，均是森林分布面积较少的省份，CV 分别为0.277、0.265 和0.253，而CV 最小的3 个省域为贵州、福建和北京，分别为0.069、0.064 和0.061。

2.2.3 不同类型森林CNC 统计特征及与全球数据的对比

反演的CNC 最低值和最高值分别为0.49 g/100g和3.63 g/100g，全国均值为2.236 g/100g，标准差为0.277 g/100g。不同森林类型中，落叶阔叶林CNC 最高，为2.47±0.18 g/100g，其次依次为常绿阔叶林、混交林、落叶针叶林和常绿针叶林。其中常绿针叶林CNC 最低，为1.82±0.28 g/100g（表4）。

根据美国橡树岭实验室基于全球文献755个树木叶片的氮浓度数据（Oak Ridge National Laboratory Distributed Active Archive Center），全球绿叶氮浓度最低值和最高值分别为0.08 g/100g 和5.20 g/100g，叶片氮浓度均值为1.81±0.719 g/100g（表4）。该数据集来自1970—2010年的多项研究，数据覆盖全球绝大多数地区（159.7 W～176.9 E,50 S～68.5 N）[28]。

表3 不同省市区森林CNC 统计特征†Table3 Statistical characteristics of forest CNC in different provincial areas

反演的每个栅格CNC 代表1km2区域内森林树种叶片氮浓度加权的结果。通过对比全球数据集，反演的CNC 数值范围落于全球数据集范围之内，说明反演的中国区域森林CNC 上下限值是合理的。同时，反演CNC 均值比全球均值略高，每种森林植被功能型的CNC 也高于对应树种的平均值，标准差则低于全球数据。由于反演得到的CNC 反映的是林冠中上层叶片的氮浓度信息，而全球数据集中的叶片取样位置和时间没有做明确说明。另外MODIS 的森林分类数据具有一定误差[29]，对反演的CNC 统计特征会产生一定影响。因此，我们认为反演结果是可信的。

2.3 验证与讨论

2.3.1 反演结果初步验证

根据10 个森林验证站点坐标（表2），提取反演结果中对应像元的CNC。比较实测值与反演值发现，两者相关系数为0.52，均方根误差为0.43 g/100g，占反演CNC 均值的19.23%。这说明反演结果能够反映CNC 空间差异，且反演的数值范围可接受。另外，大部分验证点分布在1:1线以上（图4），说明反演结果可能存在一定高估。

图4 中国10 个森林站点的CNC 实测值与反演值的对比Fig.4 Comparison between measurements and inversion values of CNC in 10 China forest sites

2.3.2 讨 论

森林CNC 空间分布受内因和外因两个方面的综合控制。内因指冠层尺度上树种组成的变化，即森林类型变化对于CNC 的影响。植物在长期进化过程中，其各个器官理化性质的功能性收敛使得不同物种的叶氮浓度差异巨大[30]，而与叶氮浓度关联的冠层结构的总体变化直接影响了冠层反射特征[31]。外因主要指外界环境对于森林生态系统的影响，包括气温、降水、纬度、海拔、太阳辐射、土壤养分水分、大气氮沉降、森林扰动（林火、砍伐等）等方面的影响[32]。

由于森林CNC 与林冠结构特征具有协同性，CNC 除其自身的反射特征外，还通过其相应的林冠结构特征影响其冠层反射光谱，故反演的森林CNC 主要反映了不同森林功能型之间的差别[14,17,31]。因此，CNC 的空间变异可以用森林植被功能型的空间差异来解释，而对于同一森林功能型内部的CNC 空间变异应主要体现在环境因素的影响方面，其反演的准确性亦有待进一步评估。

3 结 论

基于野外采样调查和Hyperion 高光谱数据，建立森林CNC 反演的统计模型，并借助MODIS数据，实现了2011—2013年间生长季森林CNC的初步反演和分析，得到如下结论：

1）基于温带森林和亚热带森林的调查数据显示，CNC 与NIR 之间存在显著的线性正相关关系，R2=0.75，P＜0.000 1，建立的模型可用于CNC 反演分析。

2）全国森林CNC 空间分布大致以胡焕庸线为界，东南高西北低。全国森林CNC 介于0.49～3.63 g/100g 之间，全国均值为2.236±0.277 g/100g。分省统计显示，海南、吉林和辽宁的森林CNC 均值最高，而新疆、山东和天津的森林CNC均值最低。

3）不同森林功能型的CNC 之间存在差异。落叶阔叶林CNC 均值最高为2.47 g/100g，其次依次为常绿阔叶林、混交林、落叶针叶林，常绿针叶林CNC 均值最低，为1.82 g/100g。

4）初步验证显示，反演结果在一定程度上反映了我国森林CNC 的空间特征和统计特征，但部分地区可能存在一定高估。另外，森林CNC 的空间变异主要受到森林植被功能型空间变异的影响。

致谢：感谢美国地质勘探局（United States Geological Survey,USGS）网站提供的MODIS 和EO-1 Hyperion 遥感数据！感谢中科院生态网络观测与模拟重点实验室相关老师提供的部分站点叶片氮浓度验证数据！同时对协助调查的同事和野外台站戴冠华、杨风亭副站长及工作人员的帮助深表感谢！