基于微钻阻力的树木年轮测量方法研究

姚建峰 赵燕东 符利勇 宋新宇 卢 军 黎沙沙

(1.信阳师范学院计算机与信息技术学院，信阳 464000；2.北京林业大学工学院，北京 100083；3.中国林业科学研究院资源信息研究所，北京 100091；4.国家林业和草原局森林经营与生长模拟重点实验室，北京 100091；5.信阳师范学院数学与统计学院，信阳 464000)

0 引言

在温带和亚热带地区，树木生长受季节影响比较大。在生长季早期，树木生长较快，所形成木材的密度较小，颜色较浅，通常把这部分木材称为“早材”。在生长季晚期，所形成木材的密度较大，颜色较深，通常把这部分木材称为“晚材”。当年的晚材与次年的早材之间界限分明，呈现纹轮，称该纹轮为年轮线。树干任意一处从髓心到形成层的年轮线数目就是该处树干生长的年龄，相邻年轮线之间的距离称为年轮宽度[1]。因此，树木年轮不仅记录了树木的年龄，还记录了气候、环境和森林经营等综合外界因子对树木生长的影响[2-3]。通过建立树木年轮与气候因子之间的数学模型，可重现过去的气候环境变化，甚至可能预期未来的气候变化规律[4]。同时，树木年龄是森林资源调查的一项基本指标，是林业生产中重要的时间指标，在预估蓄积量和生长量、确定林分龄级和龄组、调整林分龄级配置、评价森林碳汇潜力、制定森林经营方案等方面具有重要意义[5-8]。目前测量树木年轮常用生长锥法。生长锥法取样木芯速度较慢，且对树木生长有一定的负面影响，因此，不能使用生长锥大范围地测量树木年轮。科研人员试图寻找一种无损、快速的测量方法来代替生长锥法[9]。微钻阻力仪是使用电机控制钻针匀速钻入树木并实时记录钻针阻力的一种仪器[10]。钻针阻力与木材密度正相关。当钻针钻入晚材部分时，钻针阻力较大；当钻针钻入早材部分时，钻针阻力较小。当钻针沿径向方向钻入树木时，钻针阻力呈峰谷交替规律变化。因此，根据钻针阻力峰谷特征可获取树木年轮信息[11]。针头宽度仅3 mm，对树木损伤很小[9]，是取代生长锥成为测量树木年轮的最佳仪器。

目前世界上使用的微钻阻力仪主要由德国Rinntech和IML(Instrumenta Mechanik Labor)公司生产。由于这两个厂家生产的微钻阻力仪钻针阻力中含有大量的噪声信号，钻针阻力中的波峰不能与树木年轮的晚材一一对应，年轮识别主要通过人工完成，年轮识别工作量大。唐守正院士团队研究了微钻阻力仪工作原理及其在树木年轮测量方面的应用，并取得一系列成果[12-18]。但是，由于钻针振动，钻针有时偏离直线方向钻入树木，钻针阻力波峰位置与年轮晚材位置存在错位情况，所以，根据阻力波形图不能获取精确的年轮宽度信息。本文在前期研究基础上，对微钻阻力仪的钻针形状、钻针支撑挡片和峰谷年轮识别算法做进一步改进，以提高所研制的微钻阻力仪测量树木年轮的精度。

1 微钻阻力仪机械部件改进设计

微钻阻力仪钻针由2个电机控制：直流电机控制钻针旋转速度，步进电机控制钻针进给速度。微钻阻力仪的机械传动结构如图1所示。钻针通过钻针夹与直流电机轴相连。传动丝杆与步进电机轴相连。丝杆滑块中心有一个螺纹孔，丝杆滑块通过螺纹孔嵌套在传动丝杆上。直流电机固定在丝杆滑块上。当步进电机旋转时，步进电机带动传动丝杆同步旋转，使丝杆滑块在直线导轨上移动，从而带动钻针移动[18]。

图1 机械传动结构图Fig.1 Mechanical transmission structure diagram1.步进电机 2.联轴器 3.后丝杆支撑座 4.直流电机 5.钻针夹 6.传动丝杆 7.直线导轨 8.钻针 9.钻针套头 10.前丝杆支撑座 11.丝杆滑块

钻针针头宽度为3 mm，钻针针杆直径为1.5 mm[10]，针头宽度是针杆直径的2倍，因此，钻针阻力主要集中在钻针针头上。德国Rinntech公司生产的钻针形状如图2所示。

树木晚材呈圆弧形状，钻针距离树木髓心越近，钻针所切削晚材部分的弧度越大。晚材部分密度较大，钻针阻力较大；早材部分密度较小，钻针阻力较小。由于钻针针尖长度仅0.3 mm，钻针针尖过短，当钻针针尖两边的阻力差异较大时，钻针针尖不能固定钻针进给方向，导致钻针易向阻力较小的早材部分倾斜，使钻针不能沿直线方向钻入树木。因此，仅通过钻入点和钻出点的位置无法确定钻针路径。如果钻针针尖过长，当年轮宽度较窄时，钻针针头可能同时跨越多个年轮，使钻针阻力不能清晰地反映出树木年轮信息。为了使钻针尽量沿直线方向钻入树木，同时使钻针阻力可测量年轮宽度大于0.5 mm以上的年轮信息，本设计将钻针针尖的长度改为0.5 mm。同时使用日本生产的高强度铬钼钢SCM435制作钻针，以减小钻针的振动幅度和弯曲程度。改进前后钻针针头形状如图3所示。

由于钻针又细又长，当钻针高速旋转时，钻针易振动、弯曲、变形，因此，需要设计钻针支撑挡片减小钻针振动幅度。文献[14]中提出的钻针支撑挡片设计方案所使用的金属片长达1.5 m，不易加工，且加工精度难以得到保证。因此，本文提出一种新的钻针支撑挡片设计方案。该设计方案中使用最长的金属片仅0.8 m，易加工，有效提高了钻针支撑挡片的加工精度。单个钻针支撑挡片的结构如图4所示。

钻针从前横向支撑片中间的钻针孔穿过，直线导杆从前横向支撑片和后横向支撑片两边的直线导杆孔穿过。钻针支撑挡片可在直线导杆上移动。多个长度和高度不同的钻针支撑挡片可叠加在一起。钻针支撑挡片组装示意图如图5所示。

图5 钻针支撑挡片组装图Fig.5 Assembly diagram of baffle plate for supporting drill needle

当直流电机座在起始位置时，钻针支撑挡片尾端的横向支撑片叠加在一起，钻针支撑挡片前端的横向支撑片均匀分布在直线导杆上。当直流电机座向前移动时，直流电机座推动钻针支撑挡片前端的横向支撑片从而推动钻针支撑挡片向前移动，使钻针支撑挡片前端的横向支撑片逐步叠加在一起。当直流电机座移动到直线导轨的最前端时，所有钻针支撑挡片前端的横向支撑片叠加在一起，钻针支撑挡片尾端的横向支撑片均匀分布在直线导杆上。当直流电机座后退时，直流电机座推动钻针支撑挡片尾端的横向支撑片从而使钻针支撑挡片向后移动。

当钻针钻入晚材部分时，晚材密度增大，钻针阻力增大，直流电机转速降低，直流电机控制器增加直流电机电压，从而使直流电机电枢电流和功率增加，使电机转速上升；当钻针钻入早材部分时，早材密度减小，钻针阻力减小，直流电机转速升高，直流电机控制器减小直流电机电压，从而使直流电机电枢电流和功率减小，使电机转速降低。因此，直流电机的电压、电流和功率与钻针阻力正相关，可用直流电机的电压、电流或功率来间接表示钻针阻力[14]。由于直流电机电枢回路中存在大量的噪声信号，因此，直流电机电流信号和功率信号中也存在大量的噪声信号。前期研究结果表明：使用转速修正后的直流电机电压表达钻针阻力比较合适[18]。钻针阻力计算公式为

U1=Un0/n

(1)

式中U1——修正后的直流电机电压，V

U——直流电机电压，V

n0——直流电机设定转速，r/min

n——直流电机实际转速，r/min

2 峰谷年轮识别算法改进设计

自制微钻阻力仪使用峰谷年轮识别算法识别树木年轮[14-15]。峰谷年轮识别算法是根据钻针阻力波形图中波峰与波谷之间的差值来识别钻针阻力波形图中的波峰是否为年轮信号。该算法需要设置波峰与波谷之间差值的阈值。当波峰与相邻波谷之间的差值大于或者等于设定的阈值时，则判别该波峰是树木年轮信号；当波峰与相邻波谷之间的差值小于阈值时，则判别该波峰是噪声信号。把树木年轮信号的个数作为测量路径的年轮数，把树木年轮信号波峰极大值点作为早晚材之间的分界点，把每个树木年轮信号波峰极大值点与下一个相邻的树木年轮信号波峰极大值点之间的距离作为年轮宽度。使用C语言编写峰谷年轮识别函数实现峰谷年轮识别算法。峰谷年轮识别函数的定义如下：

int Peak_Vally(floata[][2],intk,floatb[][3],floatd[],floatdet)

峰谷年轮识别函数中，形式参数a是二维数组，用于存储钻针阻力的测量数据，第1列存储阻力的序号，第2列存储阻力；形式参数k存储钻针阻力的个数；形式参数b是二维数组，用于存储波峰和波谷的阻力数据，第1列存储阻力的序号，第2列存储阻力，第3列存储波峰或者波谷的标志，如果该阻力数据是波峰，该数组元素为1，如果该阻力数据是波谷，该数组元素为0；形式参数d是一维数组，用于存储年轮宽度；形式参数det是波峰与波谷之间差值的阈值Δ；函数返回值是树木年轮信号的个数，即数组b中波峰的个数。

阈值设定值要合理，如果阈值设定值过大，则把部分波峰与波谷差值较小的有效年轮信号判别为噪声信号；如果阈值设定值过小，则把部分波峰与波谷差值较大的噪声信号判别为有效的年轮信号。本文以每个测量数据的识别年轮误差最小的阈值作为该测量数据的最优阈值，把每个测量数据最优阈值的平均值作为所有测试数据的最优阈值，然后以该阈值预估每个测量数据的年轮数。具体方法如下：

(1)设置最优阈值的寻优范围：把测量数据波形图与圆盘图像对比，初步确定哪些波峰是树木年轮信号，哪些波峰是噪声信号，把树木年轮信号中波峰与波谷差值的最大值作为最优阈值的最大值Δmax，树木年轮信号中波峰与波谷差值的最小值作为最优阈值的最小值Δmin。当波峰与波谷之间的差值大于Δmax，把该波峰作为一个有效的年轮信号，当波峰与波谷之间的差值小于Δmin，把该波峰作为一个噪声信号。

(2)设置寻优步长：为了使获得的最优阈值较精确，需要尽量减小寻优步长。但是，当寻优步长过小时，寻优时间过长。因此，本设计把寻优范围进行200等分，即步长设置为Δmax与Δmin差值的0.5%。

(3)寻找每个测量数据的最优阈值：以每个测量数据和寻优范围内的每个阈值作为参数调用峰谷年轮识别函数，以预测年轮精度最高的阈值作为该测量数据的最优阈值。

(4)计算每个测量数据最优阈值的平均值，将该值作为所有测试数据的最优阈值Δ。

(5)以每个测量数据和Δ作为参数调用峰谷年轮识别函数。

使用峰谷年轮识别算法处理阻力数据的流程如图6所示。

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

使用自适应低通滤波器对原始钻针阻力数据进行滤波处理，然后根据钻针路径长度和树皮厚度，选取钻针路径上木质部部分的钻针阻力数据，调用峰谷年轮识别函数获取每个年轮的波峰和波谷点的数据和年轮宽度，最后使用R语言绘制每个阻力数据的年轮识别图。

3 性能试验

3.1 试验仪器

试验仪器主要有自制的微钻阻力仪、德国Rinntech公司生产的Resistograph 650-s型微钻阻力仪和Lintab 6型高精度树木年轮分析仪。自制微钻阻力仪的钻针旋转速度设置为3 500 r/min，钻针进给速度设置为30 cm/min，钻针阻力采样间距为0.005 mm，测量长度为500 mm。Resistograph 650-s型微钻阻力仪钻针进给速度是60 cm/min，钻针阻力采样间距为0.01 mm。Lintab 6型高精度树木年轮分析仪的分辨率是0.01 mm，测量长度为560 mm。

3.2 试验方法

3.2.1试验样品加工方法和处理方法

2020年9月在信阳师范学院校内天然次生林中采样9棵近2个月内枯死的马尾松作为试验材料，在树高1.3 m附近无明显缺陷位置截取10 cm厚的圆盘，并记录圆盘北向方向。使用打磨机打磨圆盘直至年轮线清晰可见为止。在圆盘光面分别使用红色、蓝色的铅笔标记东西、南北方向，两个方向都要经过圆盘髓心。使用扫描仪扫描圆盘图像，使用Lintab 6测量每个圆盘东西、南北2个方向上的年轮数和年轮宽度。

3.2.2钻针阻力取样方法

分别使用2个微钻阻力仪沿圆盘东西、南北2个方向钻入圆盘，尽量使钻针路径经过圆盘髓心。在同一个方向上，2个仪器的钻入点在垂直方向上相距2 cm左右，以防止2个仪器钻针路径重合。记录钻针路径长度和钻针钻入圆盘、钻出圆盘处的树皮厚度。

3.2.3Resistograph 650-s型微钻阻力仪年轮宽度测量方法

使用Resistograph 650-s型微钻阻力仪自带的DECOM软件识别树木年轮。首先根据钻针路径长度和树皮厚度，选取钻针路径上木质部部分的钻针阻力数据，然后使用DECOM软件自动识别树木年轮，最后保存每个测量数据的年轮识别结果图和年轮宽度数据。

3.2.4微钻阻力仪年轮识别误差计算方法

对照钻针路径圆盘图和钻针阻力年轮识别图，如果钻针阻力年轮识别图中的年轮线与圆盘图中的晚材位置相对应，则把该年轮线标记为正确的年轮线；如果钻针阻力年轮识别图中的年轮线与圆盘图中的早材位置相对应，则把该年轮线标记为错误的年轮线；如果圆盘图的晚材位置没有对应的年轮线，则把该年轮标记为未能识别的年轮。保存钻针路径圆盘图和钻针阻力年轮识别图的匹配图，并分别统计每个测量数据识别正确的年轮线个数m1、识别错误的年轮线个数m2和未能识别晚材的个数m3。每个测量路径上使用Lintab 6测量的晚材个数为m，每个阻力测量数据的年轮误判率e1、年轮漏判率e2和年轮识别错误率e计算式为

e1=m2/m×100%

(2)

e2=m3/m×100%

(3)

e=e1+e2

(4)

3.2.5微钻阻力仪年轮宽度测量精度计算方法

对照钻针路径圆盘图和钻针阻力年轮识别图，在圆盘图像晚材位置标记漏判的年轮，在年轮识别图中标记正确识别的年轮分界线和错误识别的年轮分界线。如果钻针阻力图中相邻2个正确识别的年轮分界线之间正好对应圆盘图中的1个年轮，则把这2个正确识别年轮分界线之间的距离作为微钻阻力仪测量的年轮宽度w1，把Lintab 6测量的年轮宽度作为微钻阻力仪所识别年轮的宽度真值w2；如果钻针阻力图中相邻2个正确识别的年轮分界线之间与圆盘图中多个年轮相对应，即正确识别的2个年轮线之间存在漏测的年轮，则把这2个年轮分界线之间的距离作为微钻阻力仪测量的年轮宽度w1，把Lintab 6测量的多个年轮宽度和作为微钻阻力仪所识别年轮的宽度真值w2。相邻的2个正确识别年轮线测量年轮宽度的精度ε0计算公式为

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

(5)

把每个测量数据每2个正确识别年轮线测量年轮宽度精度ε0的平均值作为每个测量数据的年轮宽度的测量精度ε。

4 试验结果与分析

4.1 自制微钻阻力仪年轮测量结果

根据钻针阻力曲线与圆盘年轮图像的对比分析，确定最优阈值的最大值Δmax为2.01 V，最优阈值的最小值Δmin为0.01 V，寻优步长设置为0.01 V，各测试数据的最优阈值如表1所示。所有测试数据最优阈值的平均值为0.68 V。把阈值设置为0.68 V，调用峰谷年轮识别函数，记录各数据的测试结果，绘制年轮识别图。通过与圆盘年轮图像对比，在圆盘图像晚材位置标记漏判的年轮，在年轮分界线上标记误判的年轮，然后计算年轮识别的漏判率、误判率、错误率和年轮宽度测量精度。表1为各测试数据的年轮识别结果，图7为其中一个测试数据的年轮识别匹配图。从图7中可以看出，峰谷年轮算法可正确识别圆盘图中约90%的年轮。

图7 自制微钻阻力仪年轮匹配图Fig.7 Tree-rings matching diagram of self-made micro drill resistance instrument

表1 自制微钻阻力仪年轮测量结果Tab.1 Results of tree-rings measurement measured by self-made micro drill resistance instrument

4.2 Resistograph 650-s型微钻阻力仪年轮测量结果

使用DECOM软件年轮识别模块自动识别各测量数据的年轮，记录各数据的测试结果，保存年轮识别图。通过与圆盘年轮图像对比，在圆盘图像晚材位置标记漏判的年轮，在年轮分界线上标记误判的年轮，然后计算Resistograph 650-s型微钻阻力仪年轮识别的漏判率、误判率、错误率和年轮宽度测量精度。表2为Resistograph 650-s型微钻阻力仪年轮识别结果，图8为其中一个测试数据的年轮识别匹配图。从图8中可以看出，DECOM软件可正确识别圆盘图中约75%的年轮。

表2 Resistograph 650-s型微钻阻力仪年轮测量结果Tab.2 Results of tree-rings measurement measured by Resistograph 650-s

图8 Resistograph 650-s型微钻阻力仪年轮匹配图Fig.8 Tree-rings matching diagram of Resistograph 650-s

4.3 结果分析

从测试结果可以看出，2个微钻阻力仪所识别的树木年轮分界线绝大部分都与树木年轮的晚材位置相对应，2个仪器的平均年轮测量宽度的精度达85%左右，说明微钻阻力仪测量树木年轮个数及其宽度是可行的。但有些测试数据年轮识别错误率较高，主要有以下原因：①部分年轮早晚材密度差异较小，钻针阻力峰谷差值较小，特别当年轮宽度小于1 mm时，钻针阻力不能清晰地反映出树木年轮信息，在识别树木年轮时，易把这些峰谷差值较小的波峰判别为噪声信号，如图9所示，自制微钻阻力仪漏判的3个树木年轮所对应的波峰的峰谷差值均较小。②部分树木中存在伪年轮，部分伪年轮的钻针阻力峰谷差值较大，容易把伪年轮的波峰识别为一个正常年轮信号，图9中左边第7条年轮分界线对应的即为伪年轮。③钻针钻入方向偏离树木径向方向，由于树木年轮一般不是一个规则的圆形，且在实际操作时，钻针很难通过树木髓心，因此，钻针钻入晚材的方向有时与年轮线切线方向不垂直，由于钻针形状是扁平的，导致钻针钻过晚材时，钻针针头与晚材的接触面积有时增大，有时减小，因此，这部分钻针阻力峰谷差值较小，特别当树木年轮宽度小于1 mm时，钻针针头可能同时跨越多个年轮，易把这部年轮信号识别为噪声信号，如图10中漏判的树木年轮信号大部分都是由该原因引起的。

图9 包含早晚材密度差异较小的年轮、窄年轮和伪年轮的年轮匹配图Fig.9 Tree-rings matching diagram included small difference in density between early-wood and late-wood,narrow tree-rings and pseudo tree-rings

图10 钻针方向偏离树木径向方向的年轮匹配图Fig.10 Tree-rings matching diagram of drilling direction deviated from tree radial direction

对比分析2个仪器的测量结果，自制微钻阻力仪年轮识别错误率为14.60%，年轮宽度的平均测量精度是86.40%，Resistograph 650-s型微钻阻力仪年轮识别错误率为27.30%，年轮宽度的平均测量精度为84.93%，自制微钻阻力仪的年轮识别错误率比Resistograph 650-s型微钻阻力仪低12.7个百分点，年轮宽度测量精度比Resistograph 650-s型微钻阻力仪高1.47个百分点，说明自制微钻阻力仪机械结构设计合理，机械部件加工精度、特别是钻针加工精度达到现有微钻阻力仪的加工精度，年轮识别方法可行。

5 讨论

前人研究表明，微钻阻力仪钻针阻力主要由木材绝干密度决定[10,19-23]，当钻针沿径向方向钻入树木时，钻针阻力可以反映出树木年轮信息。文献[24]对比分析了生长锥法和微钻阻力法测量135棵土耳其松的年轮结果，发现两种方法测量的年轮宽度的线性相关系数高达0.97，大部分年轮边界一致。文献[15] 给出了峰谷年轮识别算法的具体过程，并根据生长锥木芯测量的树木年龄，确定每个阻力数据所对应的峰谷差值的阈值，再使用峰谷年轮算法预估每个测量数据的树木年龄，研究发现，峰谷年轮算法预估的树木年龄与生长锥法测量结果相接近，但该文献没有分析峰谷算法年轮宽度的测量精度。文献[25]使用中值法去除钻针阻力信号中的伪波峰，即把所有波峰的峰谷差值的平均值作为阈值，把峰谷差值小于平均值的波峰作为噪声信号，选取部分年轮数据分析中值法的测量精度，研究结果表明中值法测得的年轮宽度与真实年轮宽度差异很小。本文在前人研究的基础上，对微钻阻力仪的机械结构和峰谷年轮识别算法做进一步改进，使用所有测量数据最优阈值的平均值作为峰谷算法的阈值，提高了微钻阻力仪测量树木年轮宽度的精度。由于钻针振动、伪年轮、钻针方向偏离树木径向方向等原因，使得部分树木年轮信号的峰谷差值小于部分噪声信号的峰谷差值，且不同树木年轮信号的峰谷差值差异较大，很难找到一个合适的阈值精确识别所有树木年轮信号。下一步将进一步研究树木年轮信号和噪声信号在哪些特征上存在显著差异，并建立这些特征变量与峰谷差值阈值的模型，使阈值随着年轮特征的变化而变化。

6 结论

(1)阐述了钻针针尖长度对钻针路径弯曲程度和可识别年轮最小宽度的影响，并对钻针形状进行了改进，有效减小了钻针路径的弯曲程度。

(2)提出了一种钻针支撑挡片设计方案，提高了钻针支撑挡片的加工精度，减小了钻针振动幅度。

(3)对峰谷年轮识别算法进行了改进，提出了年轮宽度的计算方法。自制微钻阻力仪的年轮识别错误率比Resistograph 650-s型微钻阻力仪低12.7个百分点，年轮宽度测量精度比Resistograph 650-s型微钻阻力仪高1.47个百分点，说明本文峰谷年轮识别方法可行。