树木生长轮的测算经常被用于树木年轮分析中。由Douglass A E在20世纪初创立的树轮年代学，在树木年轮分析中常通过研究树木生长层进行定年[1]。当前，树轮年代学的研究成果已经被运用到气候环境、森林生态、环境污染、考古、地质灾害等方面[2]。为了科学测算树木的生长轮，学者们提出了不同的方法。现有的生长轮测算方法和手段有：X射线技术；同位素测定法；CT(computed tomography)技术；生长锥测定法；图像技术；WinDENDRO软件分析和利用GIS(geographic information system)工具[3-9]等。上述树轮年代学的研究对象主要为针叶材。

近年来，众多学者采用不同的图像处理技术来检测树木生长轮。根据研究情况可分成两类。第1类为自动测算类。Vaz等[10]提出生长轮自动识别与测量系统，在环中心识别后，利用最小距离-最小交叉熵原理进行环边界计算。该系统不需要用户进行交互，其缺点是使用者无法针对系统测算结果不满意的问题进行修正。第2类为辅助测算类。Wang等[11]提出自动识别木材微观图中生长轮的方法，采用图像直方图均衡、数学形态学去噪、线性膨胀、优先度优先搜索查找边界，最终取得满意的结果。该方法适用于识别木材的横切面微观图，并对识别木材宏观图中的生长轮具有参考价值。王燕凤等[12]提出的树龄测量方法是对年轮图像双边滤波增强的研究，该方法用双边滤波对图像进行去噪，采用Canny算法进行边缘检测，最后根据多数投票原则和统计方法计算出树龄。宁霄等[13]提出的基于随机森林算法的分类模型可对年轮图像进行分割。杨新华等[14]提出利用小波变换处理图像，通过多阈值分析法识别年轮图像。上述方法都是有关生长轮图像增强的处理方法，但均未分析如何自动计算生长轮数量。

根据现代计算机技术飞速发展的现状，笔者提出一种基于交互模式和图像处理的针叶材生长轮测算方法，以期能够较精确地测算横截面宏观图中的针叶材生长轮数。

1 材料与方法

1.1 试验材料

本研究所用针叶材样本取自西南林业大学标本馆，均可由人工用传统方法数出生长轮数。试验对象为五针白皮松(Pinus squamaia)、杉木(Cunninghamia lanceolata)、翠柏(Calocedrus formosana)、北美红杉(Sequoia sempervirens)。4种样本的生长轮数为19～235，针叶材端面打磨后具有不同的心边材色差。

1.2 试验方法

本研究提出基于交互模式和图像处理的针叶材生长轮测算方法，流程如图1所示。

按照图1所示的流程进行年轮测算前，需做好木材生长轮测算图像的准备工作。准备工作包括锯一段木材、打磨端面、用扫描仪扫描端面、用常用图像处理工具从髓心沿径向截取宽度为N*64(N为正整数)个像素的图像，并保存成JPEG(joint photographic experts group)格式的图片。如果髓心宽度不够，在截图时应通过多截取树皮部分代替髓心。应尽量把端面磨平，否则易导致扫描图不清楚。扫描时，DPI(dots per inch)参数的选择应根据生长轮宽度来定，生长轮越窄，DPI值应越大。

如果对边缘检测结果不满意，就应重新导入图像后再给定不同阈值。对边缘检测结果是否满意，需由人在交互模式下做出判断，并非算法自动实现。例如，边缘检测结果中的轮界线数量少于人的直观判断时，就需再次导入图像，并把边缘检测阈值调小后再进行检测。

2 结果与分析

2.1 图像平滑处理

图像平滑处理是一种可以减少和抑制图像噪声的实用数字图像处理技术。为了减少和抑制图像噪声，通常使用线性平滑滤波器(低通滤波器)，其用滤波模板确定的邻域内像素的平均灰度值代替图像中的每一个像素点值，这种处理可以减少图像的“锐化”变化。模糊程度与邻域半径成正比的邻域平均法虽然会导致图像模糊，但可以有效地抑制噪声[15]。相比灰度图像平滑处理，在处理彩色图像时，采用彩色空间分量代替灰度值[16]。在RGB(红、绿、蓝)空间内，平滑处理的计算公式为：

(1)

式中：(x, y)为平滑后的彩色图像；R(x，y)、G(x，y)、B(x，y)分别表示R、G、B的彩色分量；k为邻域(Sxy)中的像素点数。

木材横截面扫描局部原图和平滑图见图2，经过平滑处理的木材横截面图像显得更加“光滑”、噪声减少。图像平滑处理前后的边缘检测图见图3，经平滑处理的图像，其噪声点明显减少。平滑处理前后的图像直方图见图4，对比可知，平滑处理后的图像直方图灰度值聚集量差异明显。

要准确获取生长轮数，对木材端面扫描图进行平滑操作是有必要的。效果的好坏取决于实际情况，建议用户在交互模式下，通过对有无采用图像平滑处理的边缘检测图效果进行对比，然后决定是否采取图像平滑处理。

2.2 图像直方图均衡化

直方图均衡化又称灰度均衡化，可以使被处理后的图像像素占有尽可能多的灰度级并且均匀分布[15]。在彩色图像处理中，直方图均衡化包括RGB空间分别对每个彩色图像分量进行直方图均衡化，以及HSI(波长、深浅度、强度)空间内仅对亮度进行直方图均衡化2种方法[16]。

RGB平滑均衡化图及平滑均衡化后的直方图分别见图5，6。对比图2b、5可知，RGB均衡化后的图像具有较高的对比度；对比图4b、6可知，RGB均衡化后的图像占有更多的灰度级，并且分布显得更加均匀。

李文书等[15]提出的概率密度函数p(x)范围为[0，1]。本研究所计算的RGB分量值色值范围为[0，255]，该范围内转换前各值概率密度函数(PDF)为：

p(x)， 0≤x≤255

(2)

由概率密度函数的性质可得：

p(x)dx=1

(3)

假设DN(x)为转换后的像素色值，0≤x≤255，则每个像素色值计算公式为：

(4)

本研究对所加载的木材横截面图像进行平滑、均衡化处理后的局部截图见图5。该图已经具有明显的边缘，为下一步边缘检测打下了基础。在交互模式下，由用户选择RGB均衡或HSI均衡。

2.3 图像边缘检测

用阈值方法从图像中自动提取不需要任何形状(空间关系的信息)的低层次特征称为边缘检测，其目的是制作一个线图。最受欢迎的算子有Sobel、Canny和Marr-Hildreth算子。在噪声多时，Roberts交叉算子具有更好的边缘检测效果[17]。彩色图像的边缘检测比灰度图像更复杂，因为不仅要考虑亮度，还要考虑各种不同颜色[16]。

在交互模式下，用户可选择3种不同的彩色图像边缘检测方法。本研究需用边缘检测算法找出生长轮的轮界线，供用户选择使用的检测方法包括：向量法Ⅰ、向量法Ⅱ、RGB分量直接梯度法[16]。本研究推荐使用向量法Ⅱ。

向量法Ⅱ用距离(一般指欧氏距离)来定义彩色梯度，设向量v1和v2分别代表彩色图像中的2个像素点，则这两点之间的欧氏距离DE(v1，v2)定义为：

(5)

式中：下标R、G、B表示RGB彩色空间中的3个分量[16]。本研究以Sobel方法为例，彩色图像的彩色梯度(ER)定义为：

ER=max(DE(va,vb),DE(vc,vd))

(6)

式中，向量va、vb、vc、vd表示图像中3个像素点的组合[16]：

va=[v(x-1,y-1),2v(x-1,y),v(x-1,y+1)]

vb=[v(x+1,y-1),2v(x+1,y),v(x+1,y+1)]

vc=[v(x-1,y-1),2v(x,y-1),v(x+1,y-1)]

vd=[v(x-1,y+1),2v(x,y+1),v(x+1,y+1)]

(7)

在此基础上，根据用户给定的阈值，任何像素点的RGB色值若大于阈值，则赋值255，否则赋值0。因此，经过边缘检测的图像将转换成一幅黑白图像(即二值图像)，如图3所示。本研究的图像为8位灰度图像，黑色为0，白色为255，255对应二值图像中的1。二值图像的获得，将成为后续进行二值图像基本形态学运算的前提。

2.4 图像腐蚀

二值图像的基础形态学运算包括腐蚀运算、膨胀运算、开运算、闭运算和击中击不中变换等[15-16]。基于分析测算针叶材生长轮的需要，本研究采用腐蚀运算和膨胀运算，这两个运算将分别运用于不同的需求。

二值图像腐蚀运算的数学表达式为：

g(x,y)=erode[f(x,y),B]=AND[Bf(x,y)]

(8)

式中：g(x, y)为腐蚀后的二值图像；f(x, y)为原二值图像；B为所处理图像的结构元素[16]，结构元素小于原二值图像f(x, y)。Bf(x, y)定义为：

Bf(x,y)={f(x-bx,y-by),(bx,by)∈B}

(9)

B在二维结构中是一个二维矩阵，b为该结构中的具体元素，即矩阵中的元素[16]，其值为0或255。与运算(AND(x1，…，xn))定义为：当且仅当x1=…=xn=1时，运算结果为1，其他结果为0[16]。

边缘检测后，对所得到的二值图像进行水平方向的腐蚀运算，该运算使被处理图像沿接近垂直轮界线方向形成条状白色区域，该白色区域与边缘检测后所突显出来的生长轮轮界线两端连接，这是下一步尽量消除噪声(图7中的白色小块)的准备工作。在交互模式下，用户需选择与生长轮轮界线一致的方向进行腐蚀。

与图3b相比，图7显示与生长轮相垂直的两侧分别出现了一条白色区域，该白色区域与白色的生长轮轮界相连接。这一运算使得两侧区域与生长轮轮界变成相连通的区域，而不是独立的小区域，但生长轮界之间的小白点并未与其他白色区域连通成大白色区域。

2.5 图像像素间的连通性

在具有0和1值的二值图像中，当2个像素相邻且具有同一灰度值时，则被称为连通[16]。本研究将对边缘检测之后的二值图像进行白色像素点(色值为255)连通性检测，同时计算连通区域的大小(像素点个数)，然后再根据用户给定的区域大小阈值，把小于阈值的连通区域值的像素点值赋为0。

为了进行连通性检测，设p为位于图像坐标(x, y)中的一个像素，则p的4个水平和垂直相邻像素的坐标分别为(x+1, y)、(x-1, y)、(x, y+1)、(x, y-1)，组成了p的4个邻域，用N4(p)表示[16]。像素p的4个对角邻像素的坐标为(x+1, y+1)、(x+1, y-1)、(x-1, y+1)、(x-1, y-1)，用ND(p)表示。ND(p)与N4(p)合称为p的8邻域，用N8(p)表示[16]。

对二值图像进行白色像素点连通性检测时，令V为具有相同灰度值的集合，把具有255值的像素归入连接，则V={255}。如果有2个像素点p和r都在V中取值，且r在N8(p)中，则判定p和r为8连接[16]。假设S1和S2分别为一幅图的2个图像子集，如果存在S1中某些像素与S2中某些像素连接的情况，即可得到这2个图像子集是相连接的结论[16]。如果S1和S2相连接，则把S2合并到S1中，成为新的S1子集。笔者在实现具体算法时，将遍历判断二值图像所有像素点的连接情况，检测出所有白色像素点的连通子集。

通过检测图7的白色连通子集之后，把连通区域像素个数小于给定阈值的连通区域转换为黑色区域的效果见图8。相比图7，图8基本上不存在小的白色区域。本研究采用的方法是把边缘检测结果矩形图像左右两侧的颜色处理成与轮界线一致，并使两侧区域与轮界区域构成连通的区域，再通过连通性判断，把其他噪声基本消除。

2.6 图像膨胀

二值图像膨胀运算的数学表达式为：

h(x,y)=dilate[f(x,y),B]=OR[Bf(x,y)]

(10)

式中:h(x, y)为膨胀后的二值图像[16]。或运算(OR(x1，…，xn))定义为：当且仅当x1=…=xn=0时，运算结果为0，其他结果为1[16]。

本研究需用图像膨胀运算，把生长轮轮界线出现的“双轮界线”清除，把断开的轮界线连接。在交互模式下，用户需选择与生长轮轮界线一致的方向进行膨胀。需通过膨胀运算处理的情况见图9，图9是在对边缘检测图进行去噪、黑白交换之后的效果。

2.7 基于稳定状态多数投票计数法

在使用多数投票计数法计算生长轮数量时，为了尽量膨胀轮界线，需对二值图像的黑白值进行交换。经过交换的二值图，不仅可以达到膨胀时延伸轮界线的目的，而且有利于在交互模式下直接在软件界面标注生长轮序号。

膨胀处理一次，测算一次生长轮数。经过多轮测算，易得到一个稳定的生长轮数。生长轮测算结果散点图见图10，从第9次至第18次测算，10次的测算结果都得到相同的104生长轮数，占总测算次数的50%。测算生长轮数为102的占20%。其余测算生长轮数远小于50%。在交互模式下，用户使用多数投票计数法得到生长轮数为104。

在算法实现时，把散点图转为折线图，然后在折线图上计算斜率为0的最长线段(与X轴平行)。设该折线图最大测算序号为xm，该线段左端点对应的测算序号为xi、对应的生长轮数为yi，右端点对应的测算序号为xj。如果(xj-xi+1)*100%/xm≥50%，则yi为最终生长轮数。

在基于稳定状态多数投票计数法无法取得结果的情况下，需由使用测算系统的人来做出判断。在交互模式下，由人对原图和操作图进行对比辨别。如果操作图基本不存在轮界线断开和“双轮界线”的情况，则停止图像膨胀，然后给生长轮标序并取得最终生长轮数。

3 试验与讨论

作为试验的重要组成部分，测算生长轮的交互系统由作者自主开发。该系统基于.NET框架，用C#语言实现，被命名为“针叶材宏观生长轮测算系统”。其交互性体现在用户不仅可以即时看到每一个处理过程的结果，而且可以在图像平滑、边缘检测、图像腐蚀、图像膨胀等过程选择具体处理算法或设置阈值，如图11所示。用户如果对边缘检测不满意，可以选择重新导入图像。用户如果认为生长轮测算值不稳定，可以选择继续对图像进行膨胀后测算生长轮。“针叶材宏观生长轮测算系统”为用户提供的反馈，除了及时在界面上显示处理效果，还允许用户选择查看直方图(图6)、给生长轮标序号(图12)和保存处理结果。该测算系统具有良好的交互性。

选取4个具有一定代表性的样本端面扫描图(图13)进行试验。其中，图13a、b分别为五针白皮松和杉木经打磨之后的圆盘扫描图，这2个样本圆盘的图像色差小，五针白皮松样本生长轮基本呈圆形，杉木生长轮多数不呈圆形。图13c、d分别为翠柏和北美红杉经打磨之后的端面条状扫描图，这2个样本的生长轮数都超过100，翠柏色差较小，而北美红杉色差较大，且生长轮密集。

本研究提出的测算方法有相对误差小的效果，除了得益于交互模式下用户可给测算方法指定不同阈值，解决针叶材生长轮测算方法无法适应不同针叶材树种、不同样品、不同图像质量的难题，还得益于把边缘检测结果矩形图像左右两侧的颜色腐蚀成与轮界线一致，并使两侧区域与轮界区域构成连通的区域，再通过连通性判断，把其他噪声基本消除。针叶材生长轮经边缘检测之后采用不同算法去噪的对比见图14。由图14可知，本研究提出的去噪方法可以对边缘检测结果图进行进一步去噪，并且去噪效果理想，而其他算法很难在此过程进行理想的去噪。

4个样本的端面扫描截图及其测算结果见图15，由于翠柏和北美红杉生长轮数较多，因此，对其测算结果采用截图方式展示。

4种样本的测算结果见表1。其中，人工计数是人工对样本不同方向上的不含髓心生长轮计数的算术平均值，系统计数为测算系统所测算出的生长轮数。由表1可知，本研究提出的“针叶材宏观生长轮测算系统”对于生长轮测算的相对误差不大于5%。

4 结 论

本研究提出了基于交互模式和图像处理的针叶材生长轮测算方法，该方法是若干个现有图像处理算法的综合应用。通过对针叶材端面的截图，在交互模式下，经过图像平滑、直方图均衡化、边缘检测、腐蚀，以及判断连通性后去噪、膨胀等处理过后，基于稳定状态多数投票计数法结合人工交互判断得到的生长轮数，与人工计数得到的生长轮数相对误差不超过5%。

本研究提出的交互模式允许用户在具体测算过程中给定不同阈值，解决了针叶材生长轮测算方法无法适应不同针叶材树种、不同样品、不同图像质量的难题。本研究的关键创新点是对图像进行腐蚀后再进行连通性判断，最后去除面积小的区域，从而达到消除边缘检测结果图像噪声的效果；计算轮界线之间的区域个数，而不是计算轮界线数量是本研究的另一个创新点。

本研究所提出的测算方法对心边材颜色差异较大的针叶材生长轮图像，需要使用者先对图像手动分段或对生长轮边缘检测结果图像调整大小后再进行后续处理，这是本研究的不足，需要在后续的研究中加以改进。

参考文献：

[1] FRITTS H C. Tree rings and climate[M]. London: Academic Press,1976: 376-412, 534.

[2] 方克艳, 陈秋艳, 刘昶智, 等. 树木年代学的研究进展[J]. 应用生态学报, 2014, 25(7): 1879-1888. DOI: 10.13287/j.1001-9332.2014.0125.

FANG K Y, CHEN Q Y, LIU C Z, et al. Research advances in dendrochronology[J]. Chinese Journal of Applied Ecology, 2014, 25(7): 1879-1888.

[3] SNCHEZ-SALGUERO R, CAMARERO J J, HEVIA A, et al. Testing annual tree-ring chemistry by X-ray fluorescence for dendroclimatic studies in high-elevation forests from the Spanish Pyrenees[J]. Quaternary International, 2018. DOI: 10.1016/j.quaint.2018.09.007.

[4] LOADER N J, McCARROLL D, BARKER S, et al. Inter-annual carbon isotope analysis of tree-rings by laser ablation[J]. Chemical Geology, 2017, 466: 323-326. DOI: 10.1016/j.chemgeo.2017.06.021.

[5] DE MIL T, VANNOPPEN A, BEECKMAN H, et al. A field-to-desktop toolchain for X-ray CT densitometry enables tree ring analysis[J]. Annals of Botany, 2016, 117(7): 1187-1196. DOI: 10.1093/aob/mcw063.

[6] 白雪, 范泽鑫. 哀牢山中山湿性常绿阔叶林水青树年轮宽度对气候变化的响应[J]. 林业科学, 2018, 54(3): 161-167. DOI: 10.11707/j.1001-7488.20180317.

BAI X, FAN Z X. Response of tree ring width to climate change of tetracentron sinensis in humid evergreen broad-leaved forest in the middle Ailao mountains[J]. Scientia Silvae Sinicae, 2018, 54(3): 161-167.

[7] SUNDARI P M, KUMAR S B R, SAHAYARAJ A J. An approach for analyzing the factors recorded in the tree rings using image processing techniques[C]//2017 World Congress on Computing and Communication Technologies (WCCCT), 2-4 Feb. 2017, Tiruchirappalli, India, 2017: 236-239. DOI: 10.1109/WCCCT.2016.65.

[8] WOOD L J, SMITH D J. Intra-annual dendroclimatic reconstruction for northern British Columbia, Canada, using wood properties[J]. Trees, 2015, 29(2): 461-474. DOI: 10.1007/s00468-014-1124-9.

[9] 孟盛旺, 周光, 于健, 等. 利用GIS工具测定树木圆盘年轮的方法[J]. 生态学杂志, 2017, 36(6): 1746-1751. DOI: 10.13292/j.1000-4890.201706.018.

MENG S W, ZHOU G, YU J, et al. Measuring stem disk rings using GIS tools[J]. Chinese Journal of Ecology, 2017, 36(6): 1746-1751.

[10] VAZ C, CARVALHO P, DUARTE F, et al. A vision-based system for automatic growing ring detection and measurement[J]. Computers & Industrial Engineering, 2004, 46(2): 347-354. DOI: 10.1016/j.cie.2003.12.014.

[11] WANG H J, QI H N, ZHANG G Q, et al. An automatic method of tree-rings boundary detection on wood micro-images[C]//2010 International Conference on Computer Application and System Modeling (ICCASM 2010), 22-24 Oct. 2010, Taiyuan, China, 2010: V2-477-V2-480. DOI: 10.1109/ICCASM.2010.5620573.

[12] 王燕凤, 冯海林, 杜晓晨, 等. 一种对年轮图像双边滤波增强的树龄测量方法[J]. 林业工程学报, 2017, 2(5): 109-114. DOI:10.13360/j.issn.2096-1359.2017.05.019.

WANG Y F, FENG H L, DU X C, et al. Improving the tree-age measurement by enhancing the tree-ring image with bilateral filter[J]. Journal of Forestry Engineering, 2017, 2(5): 109-114.

[13] 宁霄,赵鹏.随机森林算法在树木年轮图像分割中的应用[J].林业工程学报, 2018, 3(4): 125-130. DOI: 10.13360/j.issn.2096-1359.2018.04.021.

NING X, ZHAO P . Image segmentation of tree ring based on the random forest algorithm[J]. Journal of Forestry Engineering, 2018, 3(4): 125-130.

[14] 杨新华, 赵娟娟. 基于多阈值分析法的树木年轮特征识别方法研究[J]. 江苏农业科学, 2012, 40(8): 178-180. DOI:10.3969/j.issn.1002-1302.2012.08.071.

YANG X H, ZHAO J J , Study on growth ring feature recognition method based on multi-threshold analysis[J]. Jiangsu Agricultural Sciences, 2012, 40(8): 178-180.

[15] 李文书，赵悦. 数字图像处理算法及应用[M]. 北京：北京大学出版社, 2012.

LI W S, ZHAO Y. Digital image processing algorithms and applications[M]. Beijing: Peking University Press, 2012.

[16] 赵春江. C#数字图像处理算法典型实例[M]. 北京: 人民邮电出版社, 2009.

ZHAO C J. C# examples of typical algorithm of digital image processing[M]. Beijing: Posts & Telecom Press, 2009.

[17] BONNY T, HENNO S. Image edge detectors under different noise levels with FPGA implementations[J]. Journal of Circuits, Systems and Computers, 2018, 27(13): 1850209. DOI: 10.1142/s0218126618502092.