基于林内图像的单位面积碳储量估计方法

王雪峰 陈珠琳 管青军 刘嘉政 王 甜 袁 莹

(1.中国林业科学研究院资源信息研究所 北京 100091；2.内蒙古莫尔道嘎林业局 额尔古纳 022357)

森林是陆地生态系统中最大的碳库，森林植被碳储量约占全球植被碳储量的77%，森林土壤碳储量约占全球土壤碳储量的39%，森林每年固定的碳约占整个陆地生态系统的2/3，在全球碳循环和缓冲气候变化过程中具有不可替代的作用(Brownetal.，1993；Johnstonetal.，2004)。据第八次全国森林资源清查结果，我国森林面积现已达2.08亿hm2，森林覆盖率为21.63%，其中，天然林面积1.22亿hm2，占全国森林总面积的58.65%，可见，明确天然林生态系统碳储量和碳库分配特征对更好了解我国森林生态系统碳库和未来森林生态系统碳汇均具有重要意义(魏亚伟等，2015)。

随着温室效应加剧，人们逐渐意识到森林在全球碳循环中的重要作用，但森林究竟是充当碳源还是碳汇角色仍存在很大争议(Caoetal.，1998)。为了解全球碳计量，学者们对世界范围内不同区域、不同类型的森林碳储量进行了估算(Bracketal.，2002；Masellietal.，2009；Rutishauseretal.，2013)，采用的方法有样地实测法、生理生态法和遥感反演法等(Sileshi，2014)。森林碳储量的主要载体为森林生物量，大部分研究直接或间接以生物量的现存量乘以植被的含碳率推算得到森林碳储量。样地实测法(也称收获法)是普遍采用的估算方法，可分为皆伐法、平均生物量法和相对生长法3类(木村允，1976；佐藤大七郎，1973)，其中，皆伐法精度较高，但要耗费大量人财物力且对环境破坏大；平均生物量法样地选择受人为控制，常出现实测结果偏高的现象；相对生长法通过建立林分各部分或整体的生物量回归方程估算碳储量，灵活多变，但需要针对不同地区、不同树种分别考虑模型形式、参数估计和相容性等问题(Kondwanietal.，2018；Chojnackyetal.，2014)。生理生态法基于植物生理因子或生态环境因子建立数学模型，从而估算森林碳储量，其优点是不需要大量工作即可完成碳储量估算，但由于到目前为止人类还不完全了解植物生理生态机制，因此该方法基本停留在研究阶段。遥感反演法(Djomoetal.，2017；Main-Knornetal.，2013；Hoomanetal.，2015)在大范围森林碳储量估算方面具有一定优势，但其总体估计精度偏低，且作为数据源的卫星影像或航片信息也非用户随手可得，同时不同季节、大气效应和天气情况等都会对估计结果产生影响；虽然机载激光雷达(Rossetal.，2004；Patenaudeetal.，2004)、机载高光谱技术以及融合后的多源数据能够解决一些上述问题，但高成本等严重限制了其使用和推广力度。

森林碳储量估计是一个非常棘手的问题，人们一直在探索消耗少、精度高的碳储量估计方法，但是研究进展缓慢。随着图像或视频设备的普及，无论是通过林业现场的图像传感器还是调查员携带的数字相机、智能手机，人们可以非常容易获取反映林木疏密和高矮等信息的林内图像，那么，如果能基于这些图像估测森林碳储量，则必将极大改观现有碳储量的调查方法。鉴于此，本研究从普通图像入手，探寻适合林内图像的目标信息分类提取算法及与碳储量的关系，并尝试提出一种基于林内图像的单位面积森林碳储量估计方法。

1 研究区概况与数据

1.1 研究区概况

研究区位于内蒙古大兴安岭西北麓的寒温带针叶原始林内，该区森林覆盖率94.7%，地理坐标为120°00′26″—121°19′31″E，51°03′16″—52°06′00″N，属大陆性季风气候区的低山丘陵地带，坡度一般小于15°，海拔450～1 405 m，年均气温-5～-7 ℃，冬寒夏温，年降水量414～528 mm，乔木以兴安落叶松(Larixgmelinii)为主，另外有白桦(Betulaplatyphylla)、山杨(Populusdavidiana)、樟子松(Pinussylvestrisvar.mongolica)等；林下植物主要有兴安杜鹃(Rhododendrondauricum)、杜香(Ledumpalustre)、笃斯越桔(Vacciniumuliginosum)、五味子(Schisandrachinensis)等。

1.2 研究数据

沿着防火林带随机设置64块圆形样地，样地中心选择在垂直距离防火带100 m处，样地半径为10.3 m，中心海拔482～1 018 m。记录样地海拔、经纬度、灌草、土壤等基本信息，实测样地内全部树木的胸径、树高。使用分辨率为3 456×5 184像素的佳能单反700D相机，摄取林分纵断面和树冠2类图像。纵断面图像采集从磁北方向开始，每隔45°纵向采集1张，每点8张，拍摄时保持过光心和成像中心的连线与摄影方向地面平行，图1所示为实际样地图像。由于纵断面图像包含林地等信息，为避免林木外因素对预测结果的影响，采用交互式方式剔除地平线以下部分，即仅使用地平线上部分作为待用图像。最终计算的图像密闭度最小值为13.3%、最大值为86.8%，蓄积量分布区间为16～367 m3·hm-2，碳储量基于张俊(2008)模型根据蓄积量计算得到。

图1 样地图像

2 研究方法

林木越密集，图像前景点数越稠密，林木越高大，成像后其所占用的纵向像素数越多，即林分纵断面图像隐式包含林分密度和高度2类复合信息。这种感官上的“密集”“高大”与林分碳储量存在正相关关系。假设某空间点周围存在虚拟环幕，外围林木投影到环幕上，则环幕上林木点与整个环幕的比值就代表该空间的密集大小程度，区间为[0，1]，本研究将其称为“林分密闭度”。如果在该点用摄影机环周摄影得到n帧图像，则此n帧图像的平均密闭度可反映林木对该点周围空间的平均“填充”程度，如此，即将感官上林木对空间的占有程度由定性描述转为定量刻画，其相对客观地反映了林木密度和大小信息，这一参数必定与单位面积林木碳储量存在关系。

2.1 林分密闭度计算

欲得到单位面积林木碳储量，首先要获取林分密闭度，即需要计算某点周围全部林木图像的前景比例均值。

2.1.1 确定全局阈值 找到图像中一个灰度值ts，通过如下方法将图像变为只有前景和背景2个灰度级的新图像Inew：

(1)

如果hg是灰度为g的像素数(即灰度直方图)、a是g≤ts的灰度平均值、b是g>ts的灰度平均值，则容易满足此约束条件。其中，

Brink等(1996)采用对称交叉熵全局阈值，确定使下式dt达到最小的ts就能将图像分割成只有前景和背景2个灰度级的新图像：

(2)

对称交叉熵法给出了一个优化阈值，该阈值相比大津算法和最大熵算法得到的阈值更接近于最优阈值(Li，1993)。但由于整幅图像的等权重考虑，对于光线不均的林内图像常使树体内部出现空洞、非树体部分划归为树体的情况，因此必须兼顾当前像素的周边环境信息，最终决定该像素的归属。

2.1.2 阈值调整与前景抽出 如果g为n×m图像中第x(x=0，…，n×m)个像素的灰度，C++通过hg++遍历全部像素，即可得到该图像的灰度直方图，此时，hg是灰度为g的像素数。需要注意，在遍历前要为从黑(0)到白(255)的每一个灰度级各分配8位内存空间并清零。按照上述对称交叉熵法对图像进行全局性二值化分类，由于各像素是独立考虑的，所以常出现本来是树体像素被判断成背景、而非树体像素被判断成前景的情况。客观需求是，树体内部像素点无论灰度值多少均应归为前景、而树体外部像素点无论灰度大小均应归于背景，一种解决策略是考虑相邻像素属性并采取相应措施。

在判断焦点像素归属时，考虑其上下左右各r个像素，规定以当前焦点像素为中心的方形内，小于阈值的像素数大于γ=r(3r+2)时判为前景，否则判为背景。如图2a所示，当r=2时，a、b、c、d、o、e、f、g、h将全部被判为前景；如果r=1，则a、c被判为前景，其他被判为背景。而对于图2b，不论r=1还是r=2，u、w、v均被判为背景。可以看出，对于物体内部，如果考虑邻近周边像素属性，适当的半径和判断准则对树体有填充效果；而物体外部正好相反，具有抑制作用，这恰好与实际问题相符合。

图2 考虑邻域像素的判断方法在树体内外产生的不同效果

在对称交叉熵法处理基础上，结合邻域像素属性进行分类通常能取得较为满意的结果。图3a为一幅图像的纵剖面图；图3b为基于对称交叉熵法的计算结果，此时得到的前景面积占47.09%；图3c在获取阈值时采用同样方法，但在分类时不单纯针对单个像素，还考虑焦点像素与邻域像素的关系。邻域半径取1，当邻域内的前景像素数大于5归为同类，此时前景面积占全部面积的48.74%，增加1.65%，观察发现，树体内部应填充的部分被填充了，而树体外部无关像素被忽略。

图3 对称交叉熵法处理及在此基础上考虑邻域像素属性的算法比较

为了更加清晰起见，图4给出了兴安落叶松树体部分的一个实际图像计算结果。由于树冠、树叶等相互遮挡，树体内部光线不均，加之噪声存在，致使树体内部不同区域出现过暗或过亮斑点(图4a)；采用对称交叉熵法处理，这些斑点可能区分成背景(图4b)；加入邻域其他像素信息后，误判有所改善(图4c)。这种利用与邻域像素关系以决定当前像素归属的方法具有腐蚀和膨胀双重特性，即当焦点像素处于树体内部时容易将该点归为树体，当焦点像素处于树体外部时容易将该点归为背景，树体内部被膨胀、树体外部被腐蚀，符合客观实际。

图4 对称交叉熵法及与进一步考虑邻域像素的结果比较

需要注意的是，兼顾邻域像素属性的方法对邻域半径大小和区域内多少像素为前景时将焦点像素判为前景的入选准则很敏感，如果这2个参数选择不当，将出现较大误判。试验表明，邻域半径以不大于2为宜，入选准则选择γ∈[3r2+2r,4r2+2r]之间的一个数较为合适。本研究中所用图像取r=1、γ=6，即在以焦点像素为中心的3×3的邻域内，如果有大于6个相似像素出现，则将焦点像素归为该类。

图像被执行前景和背景二值化归类后，统计树体像素比例，即得到该幅图像所代表的林分密闭度，进一步计算旋转一周各图像林分密闭度均值，即得到摄影点的林分密闭度。很明显，该均值算法与图像重叠没有关系，且采集摄影点周围图像越全面，越能代表该摄影点的林分密闭度。

2.2 碳储量估计

由林分密闭度求算单位面积碳储量，最直接有效的方法就是使用模型。林业统计模型众多，从中筛选出效果较好且使用较多的8个模型列于表1。模型中，x为林分密闭度，y为林分单位面积碳储量。式(3)～(6)为2个参数，式(4)～(6)可首先由线性形式给出初始参数，式(7)～(10)为3个参数。需要说明的是，尽管增加参数可能带来更高的模拟精度，但未必是最佳模型；不过，由于生物统计模型中的广泛使用，本研究仍选择使用频率较高的4个3参数模型，以便于比较。

表1 碳储量预测模型选取

首先,根据建模数据计算决定系数、平均残差和均方根误差；然后，利用其他独立数据对模型进行验证，计算的参数包括平均误差、平均绝对误差、平均百分比误差和平均绝对百分比误差等，计算公式参见文献(李凤日，2004)；最后，综合考虑这些数据，评价模型拟合效果。

3 结果与分析

从64块样地中随机抽取50块样地数据用于建模，剩余14块样地数据作为独立数据用于模型检验。模型拟合使用IBM SPSS Statistics 21软件，见表2、3。

表2 模型拟合参数

从表3可以看出，2参数模型中线性模型(3)的决定系数最大，3参数模型中(8)～(10)结果接近，以逻辑斯蒂模型(9)拟合效果最好。从平均误差结果看，线性模型的绝对值最小，表明林分密闭度与林分单位面积碳储量之间可能存在线性关系，但综合考虑其他指标，逻辑斯蒂模型拟合效果相对较优。为了更清晰了解拟合效果，从2参数和3参数模型中分别选出拟合精度较高的前2个模型，绘制实测值和预测值对应的散点图(图5)，粗略地看，几个模型的预测结果接近，但仔细观察发现，逻辑斯蒂模型拟合效果在较大值和较小值的两端表现更加优良。

表3 模型拟合精度

虽然图5的总体趋势较为理想，但各点间较分散，说明还存在没有考虑到的其他因素影响拟合效果。由于大多情况下海拔越高树高越低，从而影响碳储量，为此引入一个虚拟变量k代表海拔因素。

图5 预测值与实测值的对应效果

考虑到大兴安岭北部地区海拔差异不是特别大，同时不至于使模型难于应用，以本次样地中心海拔平均值750 m为界，如果海拔高小于该数值定义k=0，否则k=1。由于不分类时线性模型平均误差最小，而逻辑斯蒂模型在均方根误差和决定系数等指标上表现最优，因此本研究仅测试这2个模型。

尝试将虚拟变量加到模型不同参数位置上，最终显示，表4中模型设置方法得到的模拟精度最高。观察表4结果，考虑海拔因素后，2个模型的模拟精度均得到很大改观，如逻辑斯蒂模型，决定系数从0.859升至0.949，平均绝对误差从9.738降至4.229。

表4 加入海拔后线性和逻辑斯蒂模型的形式、参数和精度

图6 加入虚拟变量前后直线与逻辑斯蒂模型的残差比较

4 讨论

4.1 碳储量模型形式

从本次数据模拟结果看，基于纵向图像预估单位面积碳储量，逻辑斯蒂曲线精度最高，表现出较好的适应性；但是仅有2个参数的线性方程表现出与3参数型曲线很接近的结果，由此推测，图像密闭度与单位面积碳储量之间可能存在线性关系。如果二者确实是线性关系，那么不仅给模型参数估计带来方便，同时也更容易理解，因为固定空间内随着树木不断生长即碳储量不断增加，林分密闭度也不断加大，二者之间呈线性关系最容易理解；当达到土地的最大承载力时，单位面积碳储量达到极限值，此时林分密闭度为1。以本研究模型比如海拔高度大于750 m的样地(k=1)为例，当密闭度达到1时，林地上单位面积碳储量达到133.436 2 t，与逻辑斯蒂曲线的最大值131.952 2 t接近。所以有必要继续测试图像密闭度与单位面积碳储量间的关系，如果确是线性，没必要引用更复杂的多参数方程。

4.2 郁闭度与碳储量的关系

由密闭度一个参数计算单位面积碳储量，应用时仅摄取纵向图像即可，不需要额外参数就可以估测林分碳储量，从而给应用带来方便，这也是本研究的出发点。郁闭度是森林调查时一个常用的调查因子，在遥感应用中，基于郁闭度估测林分蓄积量或林内因子也是常用做法。研究表明，碳储量随郁闭度增加呈上升趋势(刘茂秀等，2011)，因此，如果将郁闭度作为碳储量估计的一个自变量或分类变量能够进一步提高碳储量估计精度，则可以考虑继续建立基于密闭度和郁闭度的“双度”碳储量预测模型，类似于一元和二元材积表，这将给用户提供更加丰富的选择权力。由于在采集林木纵向图像时，将镜头转向冠层就可以获得冠层图像，由冠层图像计算郁闭度简单且精度高，而这不需要增加太多的内外业工作量，因此，继续研究模型中是否包含郁闭度情况下，对林分碳储量估计带来多大影响有重要应用价值。

5 结论

森林碳储量估计工作量巨大。本研究基于实际感知经验，提出了一种基于林内图像的单位面积碳储量估计方法，并以兴安落叶松为例，给出了具体估计方法、步骤和精度分析。

1)基于林内图像能够实现对单位面积森林碳储量的较高精度估计。表面上看，本研究属于图像理解探索内容，但实际上是基于林内图像求算林分密度和高度进而估计碳储量的一种方法。很明显，单位面积森林碳储量与其林木密度和高度直接相关，而林分纵向图像正好隐式包含这2类信息，同时图像方法本身又能相对准确地获取这些信息，因此，估计碳储量也仅仅是一个度量转换问题。

2)采用对称交叉熵法获得初始值，进一步考虑邻域像素属性以决定当前像素所属，能够减弱林内光线不均对图像灰度造成的影响，提高林分分类的准确性。尽管对称交叉熵法是一种较好的图像二值化算法，但是由于算法的整体统一处理，致使计算局部光照不均具有普遍性的林内图像时出现较大偏差。本研究提出利用与邻近像素关系以决定当前像素归属的方法，既保留了全局算法的优点，又增加了对具体像素的处理信息，从而大大增强了适应性，提高了分类精度。

3)将图像作为输入变量，计算出林分密闭度后根据逻辑斯蒂模型就可以估计林地上单位面积碳储量，为碳储量估计提供了一种解决之路。随着林业物联网技术的不断应用，林内图像等各类传感数据不断传回，如何抽取此大量数据中所隐含的实用价值，特别是如何理解其中的海量图像，已成为林业大数据应用的难点之一，本研究方法可以从实时图像中获取监测地碳储量信息，对于林业大数据应用是一个借鉴。

一直以来，林业上使用材积表估计林地蓄积量，对应起见，本研究基于林内图像求算碳储量(或蓄积量)的方法称为图像表法。能够看到，通过建立的图像表也可以实现对单位面积森林碳储量(或蓄积量)的估计，因为该方法通过手机或照相机扫描一周图像即可得到该点为中心的单位面积碳储量(或蓄积量)，效率要远高于目前的测径方法。因此，通过建立某一树种广泛样地的图像表，以此为基础估计该类型单位面积林地碳储量(或蓄积量)可能成为未来一种全新的快速测算方法。