R软件包MtreeRing在生长轮宽度测量中的应用1)

李爽 向玮

(北京林业大学，北京,100083)

生长轮分析作为研究树木生长与环境因子相互作用的定量分析方法，在不同空间尺度下的气候重建中得到了广泛应用[1-2]。此外，生长轮本身指示着林分生长动态，可用于模拟林分生物量[3-4]、林木枯损[5]、碳蓄积[6]，一直是生态学研究的关注热点。随着研究手段的进步，生长轮分析在考古学[7]、地质学[8]、灾害重建[9-10]等领域也取得了丰富的研究成果。生长轮宽度作为树木年代学中最常用的指标，其准确测量是树木年代学研究的前提，直接影响了年表构建和气候重建的可靠性[11]，也对测量方法提出了更高的要求。

目前，生长轮宽度测量多使用WinDENDRO或LinTAB等商业化系统，这类软硬件高度集成的专用工具为研究者提供了成熟、稳定、高精度的分析平台，但售价昂贵，使用门槛较高[12-13]。随着计算机技术的快速发展，数字图像处理技术为生长轮宽度测量带来了全新的思路和启发[14-15]。生长轮数字图像的使用在保证测量精度的前提下，不仅提高了测量效率[16-17]，也显著降低了软硬件成本[12，18-19]。对于早晚材过渡不明显的树种，图像增强技术可以提高人眼的辨识度[20]，而图像拼接技术也为LinTAB难以测量的大体积圆盘样品提供了解决方案[21-22]。

R作为一款开源软件，具有强大的统计和绘图功能，可在Windows和Mac等主流操作系统运行。相比其他软件，R的开源特性使得研究者可以根据自身需求进行二次开发，具有极大的灵活性[23]。近年来，国内外研究者开发了一系列用于生长轮分析的R软件包[17，24-26]。这些R包的出现为树木年代学提供了一个全新的数据分析平台，实现了交叉定年、去趋势、年表构建等功能[27]，是对COFECHA和ARSTAN等传统软件的良好补充，推动了树轮年代学研究的发展。

尽管R软件有许多优点，但命令行的操作方式对使用者编程水平提出了较高要求，这也成为R软件推广普及的主要障碍。目前国内少有介绍R软件在树木年代学中应用的文献，本研究目的是以实例呈现在R软件中使用MtreeRing测量生长轮宽度的流程，对其图像处理算法和误差校正原理进行简单介绍，提出需要注意的问题和使用细节，以期为研究者提供参考。

1 MtreeRing安装

MtreeRing的主要用途是从扫描图像中获取生长轮宽度序列，其R函数提供了以下功能(见表1)：①读取生长轮扫描图像；②使用边缘检测算法自动识别生长轮边界；③手动删除或添加生长轮边界；④根据边界坐标计算和校正生长轮宽度。该软件包在安装目录提供了2张分辨率为1 200 dpi的图像，其中，001.png用于演示生长轮宽度测量，而missing_pith.png用于演示髓心丢失样品的宽度校正，本研究也将使用这两张图片来演示MtreeRing的功能。

用户首先需要安装R软件，其安装包可以从官方网站(https://cran.r-project.org/)免费获取，本研究使用的软件版本为3.6.1。由于MtreeRing的图形界面是基于Shiny框架开发的，因此推荐用户使用RStudio运行该R包，以获得更好的体验。关于R软件及RStudio的安装细节，可参阅《R语言实战》一书[23]。

启动R软件后，在命令提示符“>”后输入install.packages(“MtreeRing”)命令，在弹出窗口选择CRAN镜像站点下载安装包，安装完成后输入library(MtreeRing)命令将MtreeRing加载到当前环境。本研究使用1.4.2版本(https://CRAN.R-project.org/package=MtreeRing)。

表1 MtreeRing软件包中的R函数简介

2 生长轮宽度测量

外业工作采集的样品首先应进行胶装、打磨、扫描等预处理，具体步骤可参考罗春旺[28]与孟盛旺[19]的研究。MtreeRing中一次完整的宽度测量由图像读取、边界检测、宽度计算与修正、数据导出4部分组成。

2.1 图像读取

使用library命令加载MtreeRing后，用户可以使用ring_read函数将图像文件读入R软件，目前支持的文件格式包括png、jpg、tiff、bmp，具体命令如下：

img.path<-system.file(“001.png”,package=“MtreeRing”)；

t1<-ring_read(img=img.path,dpi=1200,plot=TRUE)。

第一行命令system.file将从MtreeRing的安装目录获取图像001.png的完整路径，并以字符串格式赋值到对象img.path。第二行命令ring_read完成图像读取和绘图，并将图像数据保存到了名为t1的对象。其中，参数img指定了图像文件路径；dpi为表示图像分辨率的数值型变量，应根据研究需要选择，通常不低于600 dpi；plot为逻辑值，表示打开新的R图形窗口并绘出生长轮图像。为加快读取和绘图速度，图像不应含有无关内容，如粘贴样芯的木质凹槽；读图时，树皮和髓心应分别放置在R图形窗口的左右两侧，用户可以使用参数rotate旋转原始图像。

2.2 路径创建与边界检测

读图后，用户可以使用ring_detect函数创建测量路径并沿路径识别生长轮边界。该函数是MtreeRing的核心函数，用户可以用鼠标交互方式在R图形窗口中创建路径，并以3种可选方法自动识别边界，ring_detect的重要参数见表2。本研究将使用分水岭算法对图像中的生长轮边界进行识别，具体代码如下：

t2<-ring_detect(ring.data=t1,seg=3,method=‘watershed’,border.color=‘green’)。

该命令将开启一个或若干个R图形窗口，并显示检测到的生长轮边界(图1)。

由于R图形窗口不支持滚动缩放，因此当样品的宽度序列较长时，使用参数seg将原图像分割为等距离的若干段(图1中为3段)，以获得更好的显示效果。ring_detect自动创建了一条从图像中心穿过的虚线路径(单路径模式)，分水岭算法检测到的生长轮边界以绿色突出显示，并标有生长轮编号。该函数的图像处理参数(未在表2列出)适用于常见针叶树种，熟悉图像处理技术的用户可以修改这些参数以适应特定解剖学特征的样品。MtreeRing所使用的形态学运算及分水岭分割、Canny边缘检测等算法的细节可以参考Soille et al.[29]和史景宁等[17]文献。

实际采样中，样芯常出现髓心丢失的现象，导致靠近髓心的生长轮以一定角度倾斜(图2)，此时若直接使用边界之间的水平距离会导致测量结果偏大。这种情况下，用户可以将ring_detect的参数incline设置为TRUE，函数会将两条相互平行的路径同时添加到图像(双路径模式，见图2)，随后的测量将在两条路径上同时进行。

表2 ring_detect函数的重要参数

2.3 生长轮边界修正

当生长轮图像存在噪声(如树脂管道、物理损伤)或边界不清晰时，自动识别会出现过度检测或错过边界的现象。此时用户可以使用ring_modify函数删除或添加新的边界点。以图3为例，编号为21和23的生长轮边界位置错误，用户需要删除这两条边界并使用鼠标单击添加2个新边界，具体代码如下：

t2<-ring_modify(ring.data=t2,del=c(21,23),add=TRUE)。

该命令中，参数del表示要删除的生长轮边界编号，当有多个边界要删除时，使用c函数连接各编号；参数add为逻辑值，表示是否添加新边界。当add=TRUE时，运行该命令后ring_detect创建的R图形窗口会被依次激活(窗口左上角名称显示为ACTIVE)，添加边界的具体步骤如下：①将鼠标移动至被激活的R图形窗口内，此时光标变为十字形；②在待添加的边界与路径交点处单击一次，图像上会显示该点；③如果一个图形窗口有多个边界待添加，重复步骤2；④添加完边界后，单击右键，在菜单栏选择“停止”；⑤如果有多个图形窗口(参数seg大于1)，重复前述步骤即可。当最后一个图形窗口被关闭后，ring_modify会重新打开这些窗口，并为所有边界点重新编号，供用户检查。

2.4 宽度计算与结果导出

当所有生长轮边界都被正确识别后，可以使用ring_calculate函数计算生长轮宽度，代码如下：

rw<-ring_calculate(ring.data=t2,seriesID=“940220”)。

宽度序列rw为数据框格式，参数seriesID为数据框的列名。单路径模式下，根据生长轮边界坐标计算相邻两点之间的距离，并使用图像分辨率(R)将原始单位转换为mm，公式如下：

(1)

(2)

式中：xi和yi分别为第i个生长轮边界的横、纵坐标；Di为第i和i+1个生长轮之间的原始距离；di为实际生长轮宽度(mm)。

在双路径模式下，ring_calculate函数使用生长轮边界倾斜角度αi对相邻边界之间的水平距离Di进行校正。倾斜角度定义为生长轮边界的切线与路径法线之间的夹角，当倾斜超过5°时，应考虑对宽度序列进行误差校正[28]，所用公式如下：

(3)

αi=arctan(mi)；

(4)

(5)

式中：(Xi,Yi)和(xi,yi)分别为第i个生长轮的边界与上、下路径的交点；mi是第i个生长轮边界切线的斜率；Di是根据三角函数关系校正后的生长轮宽度。

需要指出的是，髓心丢失的样品最内侧2～3个生长轮常为圆弧状(图2)，且只能与一条路径相交，因而无法使用双路径模式校正误差。这种情况下，由于最内侧弧线通常为同心圆[30]，用户可以使用pith_measure函数定位髓心，并完成宽度测量与校正。该函数带有插图的示例可以在GitHub下载(https://ropensci.github.io/MtreeRing/articles/pith-MtreeRing.html)。

ring_calculate所得的宽度序列可以直接在R软件中进行交叉定年和去趋势，也可以使用dplR包导出，具体代码如下：

library(dplR)；

write.rwl(rwl.df=rw,fname=“C:Users1.rwl”,format=“tucson”)。

第一行命令将dplR包加载到当前环境，第二行命令将数据框rw以树木年代学中常用的tucson格式保存到用户电脑，参数fname指定了保存路径和文件名。

3 图形用户界面

随着面向对象的发展，图形用户界面逐渐成为软件设计的一种趋势，它不仅减轻了使用者的学习负担，还提供了更友好的交互体验。MtreeRing的图形界面是基于Shiny框架开发的。该框架将各种R函数以下拉菜单、按钮、工具栏等图形对象的形式呈现给用户，用户不需要了解和掌握R软件的数据结构和各类函数，从而可以直观、方便地执行生长轮宽度的测量任务。

加载MtreeRing后，用户可以输入ring_app_launch()命令来启动小程序。该命令会打开系统默认浏览器并加载小程序。小程序中的宽度测量由图像上传、路径创建、边界检测与修正、数据导出4部分构成。与命令行版本不同的是，小程序支持更多路径模式，鼠标与生长轮图像的交互也更加灵活和简便。

3.1 图像上传

小程序开启后停留在图像加载(Image Loading)页面，页面左侧为导航栏，右侧为主窗口。图像上传的具体步骤如下：①主窗口左下角的图像上传(Image Upload)模块中点击浏览(Browse)按钮；②从弹出窗口选择图像，等待上传完成；③点击加载(Load)按钮。

图像上传时要将髓心一侧放在在窗口右端，可以使用图像旋转(Image Rotation)模块调整图像的方向。图像中的无关内容(如粘贴样芯的凹槽)，可使用图像剪裁(mage Cropping)模块删除以加快后续处理速度。

3.2 路径创建和边界

上传图像后，点击左侧导航栏内测量(Measurement)按钮切换页面内容。新的主窗口由工具栏区域和图像区域组成。本研究以默认的单线段路径(single segment path)做演示，其他路径模式的介绍可参考ring_app_launch函数帮助文档内链接。创建单线段路径具体步骤如下：①在路径选项(Path Options)模块内依次输入序列ID、图像dpi、采样年份；②点击主窗口(Main Window)中的路径创建(Path Creation)按钮；③鼠标移至生长轮图像上，在期望的路径起始点双击一次；④鼠标移至路径的结尾点，双击一次完成路径创建。

对于生长轮倾斜的样品，用户可取消勾选水平路径(Horizontal path)，该设定允许沿生长轮生长方向创建一条非水平路径。此外，用户也可以勾选生长轮倾斜(Inclined tree rings)，使用与ring_detect中相似的双路径校正倾斜生长轮带来的误差。

3.3 边界检测与修正

路径创建后，点击Ring Detection按钮切换到边界检测模式，单击绿色的运行检测(Run Detection)按钮完成边界识别过程(图4)。如果自动检测的生长轮边界存在错误，用户可以点击Ring Editing按钮切换至生长轮编辑模式，此时使用鼠标在路径线上双击就可以添加新边界。编辑模式下提供了一个被称为Brush的鼠标操作，用户在图像上按住鼠标左键并拖动鼠标可以创建一个淡蓝色矩形(图4中覆盖2、3号生长轮的淡蓝色区域)，此时点击橙色Delete Border按钮就可以快速删除被矩形覆盖的生长轮边界。

3.4 数据导出

测量完成后，宽度序列可以使用CSV或RWL格式导出，步骤如下：①点击Output模块中的CSV或RWL选项卡；②输入文件名(留空则使用序列ID作为文件名)；③点击蓝色的Download按钮下载文件。导出的生长轮宽度单位为mm，小程序使用的宽度计算和误差校正公式与2.4章节相同。

4 测量结果比较

本研究选取长白落叶松(LarixolgensisA. Henry)和冷杉(Abiesnephrolepis(Trautv. ex Maxim.) Maxim.)对MtreeRing的测量精度进行了对比。两个树种的样芯分别于2016年和2019年在吉林省汪清县金沟岭林场采集。样芯经过干燥、胶装、打磨后使用EPSON 1640XL扫描仪在2 540 dpi分辨率下扫描，结果保存为tiff格式。每个树种随机选取了8个样芯分别使用MtreeRing和WinDENDRO进行自动检测，并对边界识别结果进行目视检查与修改，其中MtreeRing使用分水岭方法和Canny边缘检测分别进行了1次测量。两款程序的图像处理参数均使用默认值，并进行了倾斜生长轮的宽度校正。

图5为3种测量方法的结果对比。由于生长轮宽度的测量结果呈现偏态分布，不满足配对t检验数据总体服从正态分布的前提假设。故本研究使用Wilcoxon秩和检验对2种软件测量的生长轮宽度差异进行了显著性检验，秩和检验使用R软件完成。结果表明，MtreeRing的两种宽度检测方法与WinDENDRO的测量结果差异不显著。对于长白落叶松，分水岭方法：n=302，V=20 836，p=0.498 3；Canny方法：n=302，V=23 988，p=0.347 7；对于冷杉，分水岭方法：n=371，V=32 220，p=0.751 7；Canny方法：n=371，V=34 260，p=0.737 6；其中：n和V分别为样本数量和Wilcoxon统计量。这表明使用MtreeRing测量的生长轮宽度可以满足精度要求。

5 讨论

本研究介绍了使用R软件包MtreeRing测量生长轮宽度的新方法，MtreeRing提供的小程序满足了不同知识背景的用户对生长轮分析的需求。在保证精度的前提下，通过数字图像处理技术加快了测量速度，为不具备LinTAB和WinDENDRO等专用工具的研究者提供了有效的测量手段。

随着计算机技术的发展，数字图像处理技术在生长轮分析中得到了国内外研究者的广泛关注。传统的显微镜分析法不仅需要大量时间和人力，也易受操作者经验影响，由于判读标准不一致而影响测量结果。对圆盘或样芯的数字化有助于样品长期保存与数据分享，也方便了交叉定年后的检查与复测工作。对于生长轮边界不明显的阔叶树种，可以增强边界对比度来改善视觉效果，突出图像所蕴含的解剖学特征。在处理超出扫描仪工作面的圆盘时，可以采取分幅扫描方法，再使用Photoshop或Photoscan进行图像拼接[22]。扫描时应根据研究目的选择适当的分辨率，对于树干解析，800～1 200 dpi较为适宜，而树木年代学分析推荐使用2 540 dpi以达到0.01 mm的测量精度。

相比其他数字图像法[12，19]，MtreeRing借助R软件平台中功能丰富的扩展包，提供了高度集成的测量方案，在单一程序内就可以完成所有步骤。形态学处理和边缘检测技术的使用实现了生长轮的自动识别，对于边界不清晰或解剖学特征复杂的样品，用户也可以采取手动测量作为替代方法。针对实际工作中常遇到的应力木和髓心偏倚现象，MtreeRing提供了多种方法校正由倾斜生长轮带来的测量误差，这也是目前一些已有方法未能解决的。与WinDENDRO的对比分析表明，二者的测量结果具有高度一致性，证明MtreeRing是生长轮宽度测量可靠的替代工具。

尽管MtreeRing对针叶树种取得了良好的识别效果，但在具有复杂解剖学特征的阔叶树图像上，分水岭方法与Canny边缘检测的效果并不理想，生长轮宽度仍依赖人工测量，如何在这类图像上实现高准确率的自动化测量仍需进一步探索。随着图像处理技术的不断发展，使用先进的机器学习技术对生长轮图像进行分析将成为未来研究的发展趋势。Fabijańska et al.[31]使用卷积神经网络(CNN)对环孔材以较高准确率提取了不同形态结构的生长轮边界，在生长轮宽度测量方面表现出良好的应用前景。由于R软件的开源特性，包括CNN在内的机器学习方法已经可以在R中使用，后续研究可将此类新的边界识别方法引入MtreeRing以提高其对复杂特征样品的适应性。综上所述，本研究使用的R软件包MtreeRing在保证测量精度的同时，降低了研究成本，简化操作过程，是对树木年代学分析软件的有益补充。