微钻阻力法测量早晚材密度的可行性试验

2022-09-28姚建峰符利勇宋新宇王雪峰赵燕东郑一力业巧林翟胜楠

林业工程学报 2022年5期

姚建峰，符利勇，宋新宇，王雪峰，赵燕东，郑一力，业巧林，翟胜楠

(1. 信阳师范学院计算机与信息技术学院，信阳 464000； 2. 中国林业科学研究院资源信息研究所，北京 100091；3. 信阳师范学院数学与统计学院，信阳 464000； 4. 北京林业大学工学院，北京 100083；5. 南京林业大学信息科学技术学院，南京 210037)

树木生长除受遗传特性影响外，还受气候、环境和森林经营等综合外界因子的影响[1]。生长轮是树干横截面上的同心圆环，通常每年形成一轮，故又称为年轮。通过建立树木年轮宽度、年轮密度等年轮参数与气候因子之间的数学模型，可重建过去的气候环境变化[2-3]。年轮宽度反映了1年内外界因子对树木生长的综合影响，仅通过年轮宽度指标所获取的树木年轮信息有限。而年轮密度包括每个年轮的平均密度、晚材密度、早材密度、晚材最大密度、早材最小密度等指标，可反映出树木生长季前期和后期外界综合因子对树木生长的影响[4]；同时，年轮密度是表征木材材质综合性变异的重要指标，通过研究年轮密度动态变异规律，可为森林经理、木材定向培育、木材加工与利用等提供理论依据[5-6]。

从20世纪40年代末开始，研究人员就开始尝试各种不同的年轮密度测量方法，其中包括体积比重法、光透视法、硬度测定法、抗张强度测定法、β微粒法等，但这些方法对树木年代学研究而言都存在着不同程度的限制条件，以至于不能被广泛地应用于年轮学研究中[7]。随着X射线技术成功应用于年轮密度的测量中，年轮密度的研究才逐渐兴起[8]。目前，X射线法仍是年轮密度测量中使用最广泛的一种测量方法[8]。但X射线法仪器设备昂贵，维护成本高，体积大、质量大，操作过程烦琐，对所测量样本要求较高，不适合野外作业，因此科研人员一直试图寻找一种可携带、操作简单的年轮密度测量方法来代替X射线法[9]。微钻阻力仪是使用电机控制钻针匀速钻入树木并实时记录钻针阻力的一种仪器[10-11]，钻针阻力与木材密度正相关。当钻针钻入晚材部分时，钻针阻力较大；当钻针钻入早材部分时，钻针阻力较小；当钻针沿径向钻入树木时，钻针阻力呈峰谷交替规律变化。因此，根据钻针阻力曲线图可获取年轮宽度和年轮密度信息[12]。

目前，微钻阻力仪主要由德国Rinntech和IML(Instrumenta Mechanik Labor)公司生产。由于钻针的振动，这2个厂家生产的微钻阻力仪钻针阻力信号中都含有部分的噪声信号，钻针阻力曲线图中的波峰不能与树木年轮的晚材一一对应，年轮识别主要通过专家人工完成，年轮识别精度主要依赖专家经验。2017年，在唐守正院士的指导下，中国林业科学研究院资源信息研究所和北京林业大学工学院成立了微钻阻力仪研制团队，试制了3个微钻阻力仪样机[13-16]，研究了钻针阻力表达方法[17]，提出了4种树木年轮识别算法[18-21]。笔者使用该团队最新研制的微钻阻力仪和德国Resistograph 650-S微钻阻力仪测量辐射松(Pinusradiata)木块的早材密度和晚材密度，以期测试微钻阻力法测量树木生长轮早晚材密度的精度。

1 微钻阻力仪工作原理

微钻阻力仪由2个电机控制：一个是直流电机，控制钻针旋转速度；另一个是步进电机，控制钻针进给速度。微钻阻力仪的机械传动结构如图1所示。钻针与直流电机轴相连，钻针旋转速度与直流电机转速相同，直流电机安装在丝杆滑块上。丝杆滑块在传动丝杆的带动下沿导轨直线移动，步进电机轴与传动丝杆轴相连，传动丝杆的转速与步进电机转速相同。当步进电机旋转时，步进电机带动传动丝杆同步旋转，使丝杆滑块在导轨上直线移动，从而带动钻针移动[16]。

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

当钻针钻入晚材部分时，晚材密度较大，钻针阻力较大，直流电机转速下降，直流电机控制器输出电压上升，使直流电机电流增大，输出力矩增加，直流电机转速上升；当钻针钻入早材部分时，早材密度较小，钻针阻力较小，直流电机转速上升，直流电机控制器输出电压下降，使直流电机电流减小，输出力矩下降，直流电机转速降低。由于直流电机的输出力矩不宜实时测量，且直流电机电流信号中存在大量噪声信号，因此，使用直流电机的电压、测定转速和设定转速来间接表示钻针阻力[17]。钻针阻力(F)计算公式为：

F=Un0/n

(1)

式中：U为直流电机电压，V；n0为直流电机设定转速，r/min；n为直流电机实际转速，r/min。

2 材料与方法

2.1 试验仪器

试验仪器主要有自制的微钻阻力仪和德国Rinntech公司生产的Resistograph 650-S微钻阻力仪。自制微钻阻力仪的钻针旋转速度n0设置为3 500 r/min，钻针进给速度设置为15 cm/min，钻针阻力采样间距为0.002 5 mm，最大测量长度为500 mm。Resistograph 650-S微钻阻力仪钻针进给速度是60 cm/min，钻针阻力采样间距为0.01 mm，最大测量长度为550 mm。该仪器使用直流电机的功率间接表示钻针阻力[10]，但是阻力单位不是功率单位W，而是该公司自定义的单位Resi。因此，仍不能确定该仪器钻针阻力具体的表示方法。

2.2 试验方法

2.2.1 样品加工和处理方法

为了能用体积法测量年轮的早材密度和晚材密度，本试验材料是从木材市场上购买的年轮宽度较宽的新西兰辐射松(Pinusradiata)木材(20 cm×20 cm×200 cm)。新西兰辐射松生长快，心材部分年轮宽度可达2 cm，年轮线清晰，木材强度适中，易于分割年轮的早材部分和晚材部分，是用体积法测量早、晚材密度的理想材料。在原木材两端截取厚度为3.5 cm的板材，每个板材加工成5个宽度为2 cm的长方体木块，共10个长方体木块。长方体的长沿木材径向方向，长方体的宽尽量与大部分年轮线平行，长方体的高与树高方向相同。首先根据年轮线弧度在板材横截面上标记出每个木块的位置，标记方法如图2所示；然后使用手锯沿标记线锯成长方体木块，木材长度为15～18 cm，并使用打磨机修整、打磨长方体木块，使木块高3 cm、宽2 cm；最后把长方体木块置于烘箱中烘烤至绝干状态。为了防止长方体木块受热变形，先把烤箱温度设置为25 ℃，然后每间隔6 h左右增加20 ℃，最后把烘箱温度设置为105 ℃，恒温干燥48 h，使木块干燥至绝干状态，关闭烘箱，使木块自然冷却至常温状态。

图2 木块加工方法Fig. 2 Method of making wood samples

2.2.2 钻针阻力取样方法

分别使用2个微钻阻力仪沿木块径向钻入，为了减小因钻针钻入年轮晚材方向的不同对钻针阻力的影响，每个仪器取2次钻针阻力，一次取样方向是从边材钻入心材，另一次取样方向是从心材钻入边材。为了防止钻针路径重合，每2个取样点在水平和垂直方向上相距1 cm左右。记录钻针阻力序号、钻针路径长度和钻针钻入木块方向。

2.2.3 年轮早材密度和晚材密度测量方法

6.dabusu qami arun-a bui ɡebel ，dalai tenɡkis ba na arun-a要问盐是哪里来的，是从大海和湖泊里来的)

使用游标卡尺测量长方体木块长ai、宽b、高h，精确到0.00 1 cm，使用电子天平测量长方体木块的绝干质量mi，精确到0.01 g。然后使用工具刀从边材方向开始依次切削每个年轮的早材或者晚材，每切削一次，使用砂纸打磨长方体木块切削处，使切削面平整光滑，最后再次测量打磨后的长方体木块长度ai+1和质量mi+1，并根据式(2)计算被切除早材或晚材部分的木材绝干密度ρ。

(2)

式中：mi和mi+1分别为第i次切削木材前后的木块质量，g；ai和ai+1分别为第i次切削木材前后的木块长度，cm；b为木块宽度，cm；h为木块高度，cm。

2.2.4 微钻阻力仪早晚材平均阻力计算方法

根据每次切除的早材或晚材的木材宽度，计算该部分对应的钻针阻力平均值。

1)平滑去噪处理。由于自制微钻阻力仪和德国Resistograph 650-S原始阻力信号中存在高频噪声信号，为了减小噪声信号对年轮密度测量精度的影响，首先对钻针原始数据进行平滑处理，平滑宽度均设置为2 mm。由于自制微钻阻力仪钻针阻力采样间距为0.002 5 mm，因此平滑窗口宽度设置为800个数据采样点；Resistograph 650-S钻针阻力采样间距是0.01 mm，因此平滑窗口宽度设置为200个数据采样点。德国Resistograph 650-S平滑处理前和平滑处理后的阻力曲线图与木块年轮的对比分析见图3。从图3中可以看出，原始阻力数据经过平滑处理后可较清晰地反映出树木年轮信息，峰谷差值较大的波峰与年轮的晚材位置相对应。

2)去钻针针杆摩擦力。尽管微钻阻力仪钻针阻力主要集中在钻针针头上，但是钻针针杆与钻孔之间仍存在一定的摩擦力，且该摩擦力随着钻孔深度的增加而增加。钻针钻入木材之前其针头处于空载状态，钻针阻力主要是微钻阻力仪机电系统空载时的摩擦力；钻针钻出木材之后钻针针头也处于空载状态，钻针阻力主要包括针杆阻力和微钻阻力仪机电系统空载时的摩擦力。为了减小针杆摩擦力对木材年轮密度测量精度的影响，采用以下步骤去针杆摩擦力：①计算钻针钻入木材之前5 mm微钻阻力仪的钻针平均阻力f0；②计算钻针钻出木材之后5 mm微钻阻力仪的钻针平均阻力f1；③按式(3)计算钻针钻入木材深度为d时去针杆摩擦力的钻针阻力F1。

F1=F0-d(f1-f0)/l

(3)

式中：F0为原始钻针阻力；d为钻针钻入木材的深度；l为钻针路径长度。

图4 2个方向阻力曲线及其平均阻力曲线对比Fig. 4 Comparison between resistance curves in two directions and its average resistance curve

2.2.5 数据统计分析、建模与测试

首先对测试数据进行异常数据处理和统计分析，删除测试数据中的异常数据，并分析早晚材绝干密度，以及自制微钻阻力仪和Resistograph 650-S所测早晚材阻力数据的最大值、最小值、平均值和标准差，然后随机抽样4/5的试验数据作为建模数据集，以每个早材或晚材的平均阻力为自变量，以每个早材或晚材的绝干密度为因变量，使用中国林业科学研究院资源信息研究所自主开发的ForStat统计软件建立每个仪器早晚材钻针平均阻力与木材绝干密度的线性模型。最后把剩余1/5的试验数据作为测试数据集，使用式(4)分别计算2个微钻阻力仪预测年轮密度的预测精度ξ。

(4)

3 结果与分析

3.1 异常数据处理结果

由于原始测试数据中可能存在异常数据，为了减小异常数据对数据分析和数据建模的影响，首先使用箱线图法查找原始测试数据中的异常数据，然后使用删除法处理这些异常数据。10个长方体木块被分割成133个早材木片和137个晚材木片(有4个木块两端均为晚材)，因此，共测得133组早材数据和137组晚材数据。分别对早材木片和晚材木片的绝干密度、自制微钻阻力仪钻针阻力和Resistograph 650-S钻针阻力进行归一化处理，把箱线图的延伸倍数设置为1.5，使用R语言分别绘制早材木片和晚材木片归一化数据的箱线图，如图5所示。

图5 箱线图对比Fig. 5 Comparison by boxplot

从图5中可以看出：①早材绝干密度的最大值、自制微钻阻力仪早材木片阻力的最大值和次大值为异常数据；②晚材绝干密度的最大值、自制微钻阻力仪晚材木片阻力的最大值和次大值为异常数据。从原始测试数据中查找对应的异常数据，然后删除包含异常数据的测试记录。

3.2 数据统计结果

对删除异常数据后的测试数据进行统计分析，正常早材木片数据有130组，正常晚材数据有134组，分析各种数据的最大值、最小值、平均值和标准差，统计分析结果如表1所示。

表1 试验数据Table 1 Test data

从表1可以看出：①新西兰辐射松所有晚材木片的平均密度是556.839 kg/m3，所有早材木片的平均密度是386.850 kg/m3，所有晚材平均密度是所有早材平均密度的1.439倍；②自制微钻阻力仪所有晚材部分的平均阻力为4.301 V，所有早材部分的平均阻力是3.706 V，所有晚材部分平均阻力是所有早材部分的1.161倍；③Resistograph 650-S所有晚材部分平均阻力为228.351 Resi，所有早材部分平均阻力是175.577 Resi，所有晚材部分平均阻力是所有早材部分的1.301倍；④单片早材绝干密度的最大值是472.116 kg/m3，单片晚材绝干密度的最小值是377.235 kg/m3。因此，新西兰辐射松晚材密度与早材密度差异较大，2个微钻阻力仪所有晚材的平均阻力均高于所有早材的平均阻力，但由于有少量单片早材的密度高于少量单片晚材的密度，所以有少量单片早材的平均阻力高于少量单片晚材的平均阻力(单个年轮早材部分或晚材部分的平均阻力)。因此，仅使用钻针阻力值不能确定木材的早晚材分界线。

3.3 早晚材钻针平均阻力与木材绝干密度线性模型

随机抽取正常试验数据中104组早材木片数据和107组晚材木片数据作为建模数据集，分别建立自制微钻阻力仪和德国Resistograph 650-S微钻阻力仪所测早晚材钻针平均阻力与木材绝干密度线性模型，2个线性模型的拟合曲线如图6所示，模型参数如表2所示。

自制微钻阻力仪所测早晚材平均阻力与木材绝干密度的线性相关系数为0.657，Resistograph 650-S所测早晚材平均阻力与木材绝干密度的线性相关系数为0.655。从拟合结果可以看出：早晚材钻针阻力与木材绝干密度呈显著正相关；自制微钻阻力仪所测早晚材平均阻力与木材绝干密度的线性相关系数高于Resistograph 650-S。

图6 早晚材钻针平均阻力与木材绝干密度线性模型拟合曲线Fig. 6 Linear fitting curves between average drill resistance and absolute dry wood density of earlywood and latewood

表2 早晚材平均阻力与绝干密度线性模型Table 2 Linear model between average resistance and absolute dry wood density of earlywood and latewood

3.4 模型测量精度

把剩余的54个早晚材木片数据作为测试数据集，分别使用自制微钻阻力仪和德国Resistograph 650-S微钻阻力仪钻针阻力与木材绝干密度之间的线性模型预测木片的绝干密度，2个仪器的测量精度如表3所示。

表3 测量精度Table 3 Measurement accuracies 单位：%

从表3可以看出，2个微钻阻力仪测量早晚材木材绝干密度的平均测量精度高达86%。为了检验2个仪器预测早晚材绝干密度与实际测量绝干密度，以及2个仪器之间预测的木材绝干密度是否存在显著性差异，做了2个仪器预测早晚材绝干密度与实际测量绝干密度间的t检验，检验结果如表4所示。

表4 t检验结果Table 4 Results of t test

从表4可以看出，2个仪器预测早晚材绝干密度与实际测量绝干密度、自制微钻阻力仪和Resistograph 650-S预测的绝干密度在显著性水平0.1上均无显著性差异。

3.5 结果分析

从测试结果可以看出，2个微钻阻力仪测量的树木年轮早晚材绝干密度与实际测量绝干密度在显著性水平0.1上无显著性差异，2个仪器测量的早晚材绝干密度平均精度高达86%，说明微钻阻力法测量树木年轮的早材绝干密度和晚材绝干密度是可行的。但有少量年轮密度测量试验数据精度较低，这主要是由以下几个原因造成的。①钻针振动：微钻阻力仪钻针较细长，当钻针高速旋转时，钻针振动幅度较大，钻针阻力信号中存在部分噪声信号。②钻针针杆摩擦力不稳定：钻针针杆变形、弯曲程度受钻针针头进给方向上阻力的影响。当钻针钻入晚材部分时，木材密度增加，钻针在进给方向上的阻力增加，钻针针杆变形、弯曲程度增加，针杆阻力增加；当钻针钻入早材部分时，木材密度减小，钻针在进给方向上的阻力减小，钻针针杆变形、弯曲程度减小，针杆阻力减小。③树木年轮晚材的柱面不规则，当钻针钻入晚材部分时，钻针针头与晚材的接触面积有时增大，有时减小。

对比分析2个仪器的年轮密度测量结果可知，自制微钻阻力仪测得的早晚材平均阻力与木材绝干密度的线性相关系数和早晚材密度的平均测量精度均高于Resistograph 650-S，说明自制微钻阻力仪机械结构设计合理，钻针阻力表示方式可行。

4 讨论

微钻阻力仪钻针阻力与木材密度正相关[10,22-26]，当钻针沿径向钻入树木时，钻针阻力呈峰谷规律变化，通过钻针阻力可以反映出树木年轮宽度和年轮密度等信息。文献[27]中使用Resistograph微钻阻力仪测量了30余株山杨的年轮宽度，使用COFECHA对获取的年轮宽度进行交叉定年和控制检验，并使用ARSTAN分别建立山杨的标准化年表、差值年表、自回归年表，研究发现山杨年轮宽度指标与该地区月平均气温显著相关(P<0.05)。文献[10]中分别使用Resistograph微钻阻力仪和X射线密度仪测量了冷杉、落叶松、云杉、松树、椴树和杨树等6种树木的阻力数据和密度数据，尽管这2种方式在测量路径、木材含水率、取样精度等方面可能有所相同，但研究发现钻针阻力曲线图与X射线微密度曲线图的变化趋势高度一致。由于钻针阻力信号中含有部分噪声信号，通过钻针阻力曲线图不易确定早晚材分界线，因此，微钻阻力法在测量年轮密度方面尚未得到广泛应用。

笔者最初试图寻找早材密度与晚材密度的分界值，并根据早晚材密度分界值来确定早晚材之间的钻针阻力分界值，但研究发现同一树木不同位置的木材密度差异较大，部分早材密度高于部分晚材密度，因此树木早材密度与晚材密度之间可能尚无统一的分界值。通过对比分析钻针阻力曲线图和木块年轮图，发现钻针阻力曲线中峰谷差值较大的波峰与树木年轮的晚材位置相对应，说明通过钻针阻力可测量年轮密度信息。最后通过人工方式切割早晚材并建立早晚材绝干密度与早晚材钻针阻力平均值之间的线性模型，研究发现2个微钻阻力仪的线性模型R2均在0.4以上，平均早晚材密度测量精度高达86%，微钻阻力法测量结果与体积法测量结果在显著性水平0.1上无显著差异，研究结果表明微钻阻力法测量年轮密度可行。但是从图6中可以看出，钻针阻力与木材绝干密度的散点图比较分散，相关性仍不高，说明钻针阻力还可能受其他因素的影响，如树脂含量、外部推力、仪器晃动、钻针进给速度、钻针旋转速度、钻针锋利程度等[28]。当使用微钻阻力仪测量活力木年轮密度时，还需要考虑木材含水率对钻针阻力的影响。

下一步将研究影响钻针阻力的因素、含水率与钻针阻力之间的关系、早晚材分界线的确定方法等，以期直接通过微钻阻力测量活力木年轮密度，为林业、森林经理学、气候学等学科研究提供一种精确、微损、快速、简便的年轮测量仪器。