基于灰度图像的峰谷年轮识别方法研究

姚建峰，郭旭展，宋新宇，王静贤，祁 柯，代琪琦

(信阳师范学院 a.计算机与信息技术学院； b.数学与统计学院， 河南 信阳 464000)

0 引言

在温带和亚热带地区，树木在生长季早期形成的木材颜色较浅，通常把这部分木材称为“早材”；在生长季晚期形成的木材颜色较深， 通常把这部分木材称为“晚材”[1]。当年的晚材与次年的早材之间界限分明，呈现纹轮，称该纹轮为年轮线。相邻年轮线之间的距离称为年轮宽度。树干基部的年轮个数记录了树木年龄，年轮宽度记录了气候、环境和森林经营等综合外界因子对树木生长的影响[2-3]。通过研究年轮宽度与气候因子之间的模型，可重建过去的气候变化情况，甚至可预期未来的气候变化规律[4]。同时，树木年龄是森林资源调查的一项基本指标，是林业经营管理中重要的时间指标，在确定林分龄级和龄组、调整林分龄级配置、评价森林碳汇潜力、制定森林经营方案等方面具有重要意义[5-8]。

目前测量树木年轮的方法主要有显微镜法、X射线法和图像法[9]。

显微镜法是使用显微镜或者放大镜，人工测量经过打磨的圆盘或者生长锥木芯获取树木年轮信息。该方法是基于实物的测量年轮方法，对样品的存储环境要求较高、检查或者复测样品困难。

X射线法是使用X射线扫描经过打磨的圆盘或者生长锥木芯的图像，然后使用X射线图像获取树木年轮信息[10]。该方法可获取准确的年轮宽度和年轮密度信息，但是该方法设备昂贵，使用成本和维护成本较高，操作复杂[11]。

图像法是使用扫描仪获取经过打磨的圆盘或者生长锥木芯图像，然后使用年轮图像处理软件测量树木年轮信息[12]。年轮图像处理软件有商业化的软件和科研人员自主开发的年轮识别软件。目前最具有代表性的商业化年轮图像处理软件是加拿大Regent公司开发的WinDendro软件[13]，该软件可进行年轮标记、获取年轮个数和年轮宽度等数据，但是，该软件价格较高，国内仅少数科研机构可以承担。

随着计算机图像技术[14-16]的发展，越来越多的科研人员使用图像处理技术获取树木年轮信息[17-19]。图像处理技术涉及图像滤波算法、早晚材边界分割算法和早晚材边界增强算法等，不同圆盘图像需要设置不同的参数才能达到较高的年轮识别精度。因此，基于图像处理技术的年轮识别方法要求使用人员具有较好的数字图像处理技术基础，熟悉相关算法的原理，掌握相关参数调试方法。由于大部分基层林业工作人员不熟悉数字图像处理技术，所以，该方法不适合广大基层林业工作人员使用。潘虹等[20]提出了一种基于微钻阻力的峰谷年轮识别算法，该算法只需设置1个参数(峰谷差值)的阈值,便可识别出钻针阻力信号中的树木年轮信号，使用简单。本文使用峰谷年轮识别算法识别圆盘图像中的年轮信息，并计算年轮识别的精度，为自动识别圆盘图像年轮信息提供一种简单、可行的方法。

1 灰度图像峰谷年轮识别原理及实现过程

年轮图像是由多个明暗相间的同心圆环构成，圆环中心是树木髓心。早材颜色较浅，图像灰度值较大；晚材颜色较深，图像灰度值较小。因此，圆盘灰度图像中经过树木髓心的每行灰度数据的曲线呈峰谷交替规律变化。在基于微钻阻力的峰谷年轮识别算法中，钻针阻力波峰对应年轮晚材，钻针阻力波谷对应年轮早材，与灰度曲线中的对应关系相反。

为了便于使用峰谷年轮识别算法处理年轮灰度曲线，本文取灰度值的补码进行年轮识别处理。由于圆盘图像中存在树脂斑点、打磨擦痕等噪声信号，如果只取1行灰度数据值进行年轮识别处理，这些噪声信号对年轮识别结果影响就较大，年轮识别精度不高。

为了减小噪声信号对年轮识别的影响，本文选取经过髓心、宽为3～5 mm的圆盘图像，然后取每列灰度补码的平均值作为该列在行方向上各像素点的灰度数据，从而把行、列两个方向上的灰度数据转化为行方向上的灰度数据。

为了减小行方向上噪声信号对年轮识别的影响，首先在年轮识别前对行方向上的灰度数据进行平滑处理，然后再对平滑后的灰度数据进行年轮识别处理。

峰谷年轮识别算法是根据波形图中波峰与波谷之间的差值来识别波形图中的波峰是否为年轮信号。该算法首先设置波峰与波谷之间差值的阈值Det。当波峰与相邻波谷之间的差值大于或者等于设定的阈值时，则把该波峰记为树木年轮信号；当波峰与相邻波谷之间的差值小于阈值时，则把该波峰记为噪声信号。将年轮信号的个数作为选取圆盘图像部分的年轮数，将年轮信号波峰极大值点作为年轮之间的分界点，将相邻年轮信号波峰极大值点之间的距离作为年轮宽度。

峰谷年轮识别函数的处理过程如下：

(1)获取图像灰度数据：根据图像名Name值，读取圆盘图像的数据，并去掉背景色的数据，然后将彩色图像数据转换成灰度图像数据。

(2)计算每列像素点的灰度平均值，并对灰度平均值取补码，得到行方向像素点的一维灰度补码数据，使用平滑算法对一维灰度补码数据进行平滑处理，平滑窗口长度为Det，并记录行方向上的像素点的个数m。

(3)根据图像分辨率，计算行方向上各像素点的位置，并存储在二维数组Value的第1列中；每个像素点的灰度补码值存储在数组Value的第2列中。

(4)比较数组Value的第2列中前2个数值的大小，确定像素数据变化趋势trend是上升段还是下降段。若Value [1][2]< Value[2][2]，则数据变化趋势trend处于上升阶段，把trend赋值为1;若Value [1][2] ≤Value [2][2]，则阻力变化趋势trend处于下降阶段或者稳定阶段，把trend赋值为0。把数值较小点的数据值和位置作为局部极小值点的数据值vally_y和位置vally_x，把数值较大点的数据值和位置作为局部极大值点的数据值peak_y和位置peak_x；

(5)当数据变化趋势处于上升阶段时，若新数据点的值大于局部极大值，则把新数据点作为局部极大值数据点；若新数据点的值比局部极大值点小，且它们之间的差值d小于设置的阈值Det，则认为该数据点是一个噪声信号；若新数据点的值比局部极大值点小，且它们之间的差值d大于或者等于设置的阈值Det，则把该数据点作为本年轮结束的标志，数据变化趋势从上升阶段变为下降阶段或者稳定阶段，把局部极大值点作为一个年轮信号进行处理，有效波峰的个数n加1，把新数据点的值和位置分别作为局部极小值点的值vally_y和位置vally_x；

(6) 当数据变化趋势处于下降阶段或者稳定阶段时，若新数据点的值比局部极小值点的阻力值小，则把新数据点作为局部极小值数据点；若新数据点的值比局部极小值大，且它们之间的差值d小于设置的阈值Det，则认为该数据点是一个噪声信号；若新数据点的值比局部极小值大，且它们之间的差值d大于或者等于设置的阈值Det，则把该数据点作为该年轮晚材部分的开始阶段，数据变化趋势从下降阶段或者平稳段变为上升阶段，把局部极小值点作为一个年轮信号的谷值点进行处理，把新数据点的值和位置分别作为局部极大点的值peak_y和位置peak_x。

灰度图像峰谷法识别年轮的流程图如图1所示。

图1 峰谷算法流程图Fig. 1 Flow chart of peak-valley algorithm

使用Matlab语言编写峰谷年轮识别函数实现峰谷年轮识别算法。识别函数的定义如下：

[Value, Value_Key, Wide, n]=Peak_Vally(Name, Resolution, Det)。

在形式参数中，参数Name是一个字符串，用来传递灰度图像的文件路径和文件名；参数Resolution用来传递图像的分辨率，单位为dpi；参数Det用来传递波峰与波谷之间差值的阈值。函数返回值中，参数Value用于存储转化后的灰度图补码数据，是一个二维数组，第1列用于存储像素点的位置，单位为mm，第2列用于存储像素点像素值的补码；参数Value_Key用于存储峰谷识别后的关键点数据，即年轮信号波峰极大值点和波谷极小值点的位置和像素数据，是一个二维数据，第1列存储像素点的位置，第2列像素值的补码，第3列存储波峰或者波谷的标志，若该像素点是波峰，则该数组元素为1；若该像素点是波谷，则该数组元素为0；参数Wide是一维数组，用于存储年轮宽度；参数n是树木年轮信号的个数，即数组Value_Key中波峰的个数。

2 材料与方法

2.1 试验材料及处理方法

试验材料是2020年9月在信阳师范学院校内的天然次生林中采样的12棵近2个月内枯死的马尾松，在树干高度为1.3 m处附近截取树干通直、无明显缺陷、厚约10 cm的圆盘，并记录圆盘北向方向。使用打磨机打磨圆盘，使年轮线清晰可见。在圆盘光面东西、南北方向两个方向上选取宽为3～5 mm木材，两个方向的木材都要经过圆盘髓心，并用记号笔做好标记。使用Lintab6年轮分析仪测量每个圆盘东西、南北2个方向上的年轮线个数和年轮宽度。把扫描仪分辨率设置为400 dpi，使用扫描仪扫描圆盘获取圆盘彩色图像，使用画图软件选取圆盘图像上标记区域的木材部分，并做好记录。

2.2 阈值确定方法

阈值设置值要合理，若阈值设置值过大，则把波峰与波谷差值较小的年轮信号判别为噪声信号；若阈值设置值过小，则把波峰与波谷差值较大的噪声信号判别为年轮信号。本文根据3个圆盘图像在不同阈值下识别结果与圆盘图像的对比分析图，选择能正确识别大部分年轮的阈值作为所有圆盘图像的阈值。

2.3 年轮识别误差计算方法

对照圆盘图像和像素数据年轮识别图，若像素数据年轮识别图中的年轮线与圆盘图中的晚材位置相对应，则把该年轮线标记为正确的年轮线；若像素数据年轮识别图中的年轮线与圆盘图中的早材位置相对应，则把该年轮线标记为错误的年轮线；若圆盘图的晚材位置没有对应的年轮线，则把该年轮标记为未能识别的年轮。

保存圆盘图和像素数据年轮识别图的匹配图，分别记录每个像素数据的预估的年轮个数m1、误判的年轮个数m2和漏判的年轮个数m3。每个取样方向上以Lintab6测量的晚材个数为m为年轮数的真值，每个图像数据的年轮误判率e1、年轮漏判率e2、年轮识别错误率e和年轮个数测量精度ε1,计算公式分别如下:

e1=m2/m,

(1)

e2=m3/m，

(2)

e=e1+e2,

(3)

ε1=1-|m1-m|/m。

(4)

2.4 年轮宽度测量精度计算方法

对照圆盘图和像素数据年轮识别图的匹配图，首先分析圆盘图像测量的年轮宽度数据，把错误年轮线所对应的年轮宽度进行合并，计算相邻的2个正确识别年轮线之间的宽度w1。然后分析Lintab6测量的年轮宽度，若2个正确识别年轮线之间没有漏判的年轮，则把Lintab6测量该年轮的宽度作为该年轮宽度的真值w2；若2个正确识别年轮线之间存在漏测的年轮，则把2个正确识别年轮线之间所有年轮的宽度和作为这2个正确识别年轮线之间的年轮宽度真值w2，最后使用式(5)计算每相邻2个正确识别年轮线测量年轮宽度的精度ε0:

ε0=1-|w1-w2|/w2。

(5)

把每个圆盘每2个正确识别年轮线测量年轮宽度的精度ε0的平均值，作为每个圆盘图像测量年轮宽度的平均精度ε2。

3 结果与分析

根据3个圆盘图像在不同阈值下的识别结果和圆盘图像的对比分析图可知，当阈值设置为6时，能正确识别圆盘中大部分年轮信息。图2是其中一个圆盘图像的年轮识别结果与圆盘彩色图像和灰度图像的对比图。从图2中可以看出，尽管原始图像中存在色斑、霉变、不清晰的年轮线等问题，但是峰谷年轮算法仍可正确识别圆盘中大部分年轮信息。

图2 年轮识别结果与圆盘图像对比Fig. 2 Comparision between tree-ring recognition results and disc image

通过峰谷年轮识别结果与圆盘图像对比，在圆盘图像晚材位置标记漏判的年轮，在年轮识别图中的年轮分界线上标记误判的年轮，然后计算年轮识别的漏判率、误判率、错误率和年轮宽度测量精度，各图像数据的年轮识别结果如表1所示。

表1 年轮识别结果

续表1

从表1中可以看出，峰谷年轮识别算法的平均漏判率为10.084%，平均误判率为4.059%，平均错误率为14.143%，年轮个数的平均测量精度为90.096%，年轮宽度的平均测量精度为89.263%。从总体识别结果来看，峰谷年轮识别算法预测树木年龄精度和年轮宽度精度高达90%左右。

部分图像数据的误判率或者漏判率较高，可能有以下原因：(1)设定的阈值对部分圆盘图像不合适：本文对所有数据都使用相同的阈值进行识别，对有些圆盘图像阈值偏大，导致早材与晚材颜色差异较小的年轮信号被识别为噪声信号，使漏判率偏高；对有些圆盘图像阈值偏小，导致把部分峰谷差值较大的噪声信号识别为年轮信号，使误判率偏高。在实际使用过程中，可通过人工方式标记不清晰的年轮线、对不同圆盘图像设置不同的阈值等方式来进一步提高年轮识别精度。

4 结语

首先使用列方向求均值法，将二维圆盘灰度数据转化为一维灰度数据，然后对一维灰度数据使用峰谷年轮识别算法获取年轮信息。经测试，峰谷年轮识别算法预测树木年龄精度和年轮宽度精度高达90%左右。该算法易于理解，使用简单，可为林业工作人员提供一种简单、可行的年轮识别方法。