电动生长锥结合Python树木年轮信息提取研究

常晨, 王秀兰*, 冯仲科, 孙素芬, 张浪, 马文苑

(1.北京林业大学, 精准林业北京市重点实验室, 北京 100083; 2.北京市农林科学院农业信息与经济研究所, 北京 100097; 3.上海市园林科学规划研究院, 上海 200232; 4.上海城市困难立地绿化工程技术研究中心, 上海 200232)

树木年轮信息是重要的林木特征参数，对树木年龄、生长量的确定和树木生长模型的建立有决定性作用。如何快速准确获得树木年轮信息一直是精准林业的研究重点。在传统测树学研究中，通常使用机械生长锥提取树木树芯，一般需要目测识别树芯年轮界限，在其早晚材中间插入大头针来测定年轮数量和年轮宽度，该方法在一定程度上解决了树木年轮信息提取的难题，但效率低，且操作步骤繁琐，无法快速准确获取树木年轮信息[1-2]。

传统的生长锥是非破坏式测定树木年龄和直径生长量的专用工具，由锥管、探针和锥柄三部分组成[3]，其构造简单，体积重量小,但钻入钻出都需要花费极大体力和时间，且钻入时若用力不均匀，极易导致生长锥折断，钻毁树伤[4]。20世纪60年代,Prestemon[5]使用电池驱动的冲击扳手作为动力单元设计了树木动力钻孔机，可钻孔深2～3 in·min-1，效率较低。Johansen[6]使用轻量发电机来操作电动发动机驱动钻头提取树芯，但发电机体积较大，不便于山地携带。 Krottenthaler等[7]发明了新型榫头切割器进行树木树芯提取，使用汽油机作为动力源，在取样高密度热带木材时取得了一定成效。 Kagawa[8]使用Makita TW450D重型电锤作为生长锥动力输出，配合行星减速箱降低转速提高扭矩，扭矩满足绝大部分树种及钻探深度需求，但其成本过高，单电钻就售价超过40 000元，一般林业站或林业工作者无法负担，因此，无法较好推广使用。改进的电动生长锥可全自动进行树木年轮取样，其装置能全自动对任意材质树木任意部位进行树木树芯提取，并可进行速度调控，包括钻头组顺时针钻入和逆时针退出，设备便携可拆装，其动力输出稳定，不存在受力不均匀导致生长锥损坏的情况，且其速度快、效率高、省时省力，日后可在林业行业大范围推广，减少了外业人员的体力劳作，大幅度提高了外业效率。

1 原理与方法

1.1 硬件设计

电动生长锥设备由精准林业北京市重点实验室独立设计，使用电动机代替人力为生长锥钻入树体提供动力，由Haglof树木生长锥(Haglof公司制造生产)、卡头、直流减速电机、电机控制器(信达公司制造生产)、动力锂电池(格氏公司制造生产)、水平助力握把、简易水准仪组成。生长锥需使用卡头与电机相连接，其材料为GCr15专用轴承钢[24]，在进行外壳设计时，为满足其力学特性及便携性需求，选用LY12cz作为壳体主要材料，委托北京三鼎光电仪器有限公司进行加工，精加工完成后进行表面喷砂氧化，提高壳体硬度[25]。经选型后使用信达公司制造的直流减速电机5D209GN-24型，其工作电压为DC24 V，满载功率为200 W，负载转速为2 500 r·min-1，其负载力矩为1.9 N·m,配合减速比150齿轮箱可输出86.8 N·m的堵转扭矩，转速为22 r·min-1。

考虑电机满负载运行时扭矩较大，为稳定固定电机与外壳，在电机外壳主体与电机主体相连接处做加厚处理，在电机顶部设有支撑板，可保证电机在满负载运转时不会发生剧烈震动，导致生长锥或树芯断裂。为保证壳体结构强度减轻壳体重量，壳体底部进行空框处理，在外壳四周及中部均采用加强筋设计，且各个部件均可拆卸，方便野外携带和后期维护维修。电动锥结构如图1所示，其中电动生长锥快拆板可快速进行内部电机及控制器拆装，方便进行调速器及电机的维修更换；在主壳体外壳后侧进行开槽，设有调速器紧固板，用于固定调速器，在壳体两侧分别进行圆形、方形开孔，用于安装调速旋钮和转向开关，壳体中部两侧加强筋处设有螺纹口，用于安装握把仪器外壳上盖为减轻重量采用梯形设计，上盖中部设有简易水准仪，钻取树芯时可根据小气泡状态调整电动生长锥姿态，保证生长锥钻入时垂直于树干；内部组件1起支撑板作用，用于固定电机与仪器主体；内部组件2用于固定控制器与仪器主体；助力握把用于作业时仪器固定；定制卡头用于连接电机主轴与生长锥方形尾槽。电动生长锥具体尺寸如图2所示。

1.2 电动生长锥性能测定

为了检验电动生长锥进行树木树芯提取时的可行性，使用e-Dynamic电子数显扭力扳手代替传统生长锥锥柄对北京市常见树种进行白杨、法桐、旱柳、国槐、梨树进行树芯钻取，分别记录转入和转出时每转动一圈其扭矩峰值，若其峰值大于电动生长锥所能输出最大扭矩，则电动生长锥不适用于该树种的树芯提取，若小于则适用。

1.3 Python年轮自动识别软件设计

1.3.1Python运行环境及函数库调用 获取树芯后，经简单打磨标记处理，使用扫描仪获得树芯年轮影像，在影像树芯区域进行裁剪，获得树芯年轮影像信息。在Windows系统支持下，启动IDLE Python，版本号3.63 Python具有强大的函数库，Python年轮自动识别软件在处理过程中利用其函数库进行图像识别及计算，详情如下：PIL作为Python强大的图像处理库[26]，拥有多种图像处理模块，其中以 Image模块最为重要和基础，图像的调用、打开、显示、波形调整及图像点坐标获取等功能都集成。其中，本文使用from PIL import Image模块进行树芯影像的读取、二值化、像素点坐标读取等处理。numpy库和math库是Python调用函数进行科学计算的基础，提供了Python对多维数组对象运算的支持，使用numpy库和math库进行树芯影像像素点坐标的选择、计算、测距、求值、输出等 为了方便年轮信息的提取后整理使用，本文使用openpyxl库作为数据输出工具，A openpyxl是一个读写Excel文档的函数库， 可将numpy库计算得出的数据按需求汇总至Execl表格中，方便进行读取使用。

1.3.2图像提取及处理 本研究的Python脚本识别算法已命名保存为holesrecog.py，使用者启动Python脚本后，可更改或选择图像识别目录和结果输出目录，将 holesrecog.py所在文件夹内新建文档holesrecog设为工作目录，在该目录下放入需要识别的树芯图像，首先通过程序自动判读，判断相片为jpg格式或png格式定义后续处理格式，使用os.path.join语句与Image.open获取图片路径读取图片信息，为了有效识别每个年轮位置上所标记的黑色圆点，将树芯图像通过image.convert语句将图像转为灰度图，使用self.binarizing(img, 100)将灰度图中灰度值大于100像素点的灰度值清零，小于100像素点的灰度值设置为255，即设定阈值对树芯图像进行二值化处理。则此时经处理后的图像为数个灰度值0的黑色像素点聚集，计算提取每个聚集的坐标便可计算出年轮数量及距离。

在获得二值化树芯图像的基础上，使用np.array语句获得二值化像素点的矩阵，建立自定义语句cordinate进行像素点坐标读取，其函数如①所示。

①def cordinate(self, arr):

w, h=arr.shape #

res=[]

tmp=[]

for i in range(w-1):

iswhileline=True

for j in range(h-1):

if arr[i][j]==0:

tmp.append([i, j])

iswhileline=False

if iswhileline:

res.append(tmp)

tmp=[]

return res

由于图像经过二值化处理，处理树芯时所做的黑色标记形成了黑色像素点聚集，每个黑色像素点聚集之间都被白色像素所分隔，通过使用两个for循环语句循环遍历二值化影像矩阵的值，在循环过程中每当遇到一个灰度值为0像素点,再遇到下一个灰度值为255像素点之前，将这些黑色像素点作为黑色像素点，聚集加入一个临时列表中，在循环结束后得到数个黑色像素点聚集列表，即为年轮处像素坐标聚集列表。为排除因虫洞或斑痕产生的噪声，使用if 50

②point1=list(np.mean(final_arr[0],axis=0)[::-1])

for arr in final_arr:

arr=np.array(arr)

1.3.3年轮距离解算及结果输出 在已知年轮中心点坐标的基础上，使用求距公式可得到相邻两个中心像素点之间的像素距离，其公式[27]如下。

(1)

式中，xn和yn为第n个聚集中心像素点横坐标和纵坐标，xn-1和yn-1为第n-1个聚集中心像素点横坐标和纵坐标，根据对树芯扫描时所设定的扫描分辨率，可使用像素距离换算公式[28]求得实际距离。

(2)

式中，D为n和n-1两点间的像素距离，单位为px；m为扫描时所选用的分辨率；dn为n和n-1两点间的实际距离，单位为mm。

计算之后得到两年轮间距，为了方便统计整理，将数据输出至Excel表格，保存至工作目录下，其计算及输出函数如③所示。

③def out_put(self):

i=0

workbook=openpyxl.Workbook()

wst=workbook.active

wst.cell(row=1, column=1).value=‘图片名称’

wst.cell(row=1, column=2).value=‘x轴坐标’

wst.cell(row=1, column=3).value=‘y轴坐标’

wst.cell(row=1, column=4).value=‘间距’

for im in os.listdir(self.pic_dir):

ext_name=im.split('.')[-1]

pic_name=im.split('.')[0]

lis=[pic_name]

lis.extend(list(np.mean(arr,axis=0)[::-1]))

point2=list(np.mean(arr,axis=0)[::-1])

d=math.sqrt((point2[0]-point1[0])**2+(point2[1]-point1[1])**2)

d0=d*25.4/m

lis.append(d0)

wst.append(lis)

point1=point2

i+=1

try:

workbook.save(self.excel_path)

1.4 操作方法

使用电动生长锥提取树木树芯，由于取样得到的树芯表面粗糙，划痕较多，且潮湿未干燥，需要对树芯进行预处理。首先对取样树芯进行风干处理，风干后使用乳白胶将树芯固定在带有半圆形凹槽木条中，并使用尼龙绳将其绑紧，避免因乳白胶凝固使得树芯变形翘起，阴干2～3 d后，分别依次使用50、200、800目砂带机打磨至表面平整且界线清晰，由于部分树种年轮早晚材在白炽灯或扫描仪光线下很难区分，需使用笔和直尺在树芯上年轮间附加黑色标记点，其标记均在一条直线上且该直线指向树心。

本文使用EPSON Perfection V19扫描仪进行扫描，扫描时根据树芯年轮的早晚材清晰程度和需求精度选择dpi(dot per inch)值，一般选择1 200 dpi进行扫描，清晰度较高速度较快。扫描时将树芯摆放在扫描仪正中央，树皮在左侧，树芯在右侧，扫描时可做暗室处理，确保图像不被外界光线干扰，扫描后保存为jgp格式，如若树芯过长，可分段重叠扫描，扫描后可使用ICE软件进行拼接处理。由于树芯在扫描得到的图像中仅占很小比例，为排除无用像素点干扰，对扫描得到的jgp文件，需要使用Photoshop进行裁剪处理，保存时需注意选择其dpi与原图一致且保存为jgp格式。在软件脚本中可设定工作目录，运行软件脚本，待数秒直至脚本主页面显示“处理完成，共处理了{}张图片，坐标见{}下out_put.xlsx文件”，数据处理完成后可在工作目录中查看结果文件，其格式为xlsx。

1.5 测试方法

测试实验于2018年1月8日至7月20日在海淀区鹫峰林场、大兴区各乡镇进行，以北京市常见树种白杨、悬铃木、旱柳、栾树、梨树为研究目标，选择标准木作为树芯取样样木，标准木的选择采用随机抽样的原则。选择样木40棵，使用450 mm采样长度、4.15 mm采样内径的Haglof树木生长锥作为钻头，分别使用电动生长锥和传统生长锥协同e-Dynamic电子数显扭力扳手钻取树芯，最终取得完整样芯37棵。对比分析不同树种钻取时力矩大小及电动生长锥适用性范围，并分别对普通机械生长锥钻取树芯与电动生长锥式设备进行工作效率对比试验，在试验过程中需要记录树木编号、树种、胸径(diameter at breast high，DBH)等信息，并以Excel数据结果为最终结果，对比和分析不同方法的测量时间，对取得的年轮树芯进行处理，并使用LinTab年轮分析仪与Python年轮自动提取脚本进行对比试验，测量树芯年轮宽度及年轮数量并以LinTab年轮分析仪为标准，评估快速测定方法的测量精度。

1.6 测试精度评价

将LinTab年轮分析仪所提取的年轮信息作为参考值，将并Python脚本所提取的年轮信息作为测量值，对各树种进行精度分析，评价指标包括值偏差(Bias)、相对偏差(Bias%)、均方根误差(RMSE)、相对均方根误差(RMSE%)。

(1)

(2)

(3)

(4)

2 结果与分析

2.1 电动生长锥适用性范围分析

将5种树种的测试扭矩作为参考值，结合生长锥钻入(出)深度值绘制折线图来直观表现扭矩随钻入深度的变化情况如图3所示。

由图3可以看出，生长锥钻入梨树树体时，得到获取树芯所需的最大扭矩为38.8 N·m，电动生长锥可提供最大的堵转扭矩为86.8 N·m，满足生长锥钻入退出需求，且对白杨、栾树、梨树、旱柳、悬铃木适用性良好，即类似材质树种均可使用电动生长锥进行树芯提取。

2.2 作业效率分析

为了分析电动生长锥及Phthon脚本作业效率，在外业实验中同时使用传统生长锥和电动生长锥选择相同树种同等胸径的样木进行树芯获取，在内业中同时使用LinTab年轮分析仪和Python脚本对所有树芯作业所需要的时间，测量结果如表1所示。

表1 作业效率对比Table 1 Comparison of job efficiency

由表1可以看出，使用电动生长锥和传统生长锥进行相同树木对向树芯提取时，传统生长锥组采用2人交替作业的方式进行树芯钻取以保证均衡体力，电动生长锥组仅需1人即可完成20棵树芯提取，作业时间优于传统生长锥组。Python脚本组主要耗时在对处理好的树芯进行标记和扫描处理，其工作流程自动化，无需过多人为干预，其效率优于LinTab年轮分析仪组，且实验后获得的树芯扫描影像可以永久保留，不受储存环境的影响和限制，方便后续继续研究。

2.3 年轮宽度精度分析

为了检验树木年轮自动识别软件进行树芯年轮信息自动提取技术的准确性，对获取得到的37根树芯分别使用树木年轮自动识别软件和LinTab年轮分析仪提取年轮信息，树种及其基本信息如表2所示，其中白杨13棵、悬铃木9棵、旱柳9棵、栾树2棵、梨树4棵，其胸径分布在22.2～43.2 cm 之间，由于不同树种生长速度、木质类型具有差异性，固将不同树种年轮测量精度分开讨论。

表2 取芯样木胸径统计Table 2 DBH statistics core sample

2.3.1年轮测量值拟合分析 分别对树木年轮自动识别软件和LinTab年轮分析仪得到的年轮宽度测量值与参考值进行线性拟合(图4)，其拟合结果均与直线y=x较为相近，其中悬铃木拟合方程为y=0.995 92x-0.005 59，拟合直线的R2为0.998 17；栾树拟合方程为y=1.000 43x-0.007 21，拟合直线的R2为 0.998 64；旱柳拟合方程为y=0.991 93x+0.125 93，拟合直线的R2为0.986 3；梨树拟合方程为y=1.000 43x-0.007 21，拟合直线的R2为 0.998 64；杨树拟合方程为y=0.989 13x+0.018 67，拟合直线的R2为0.993 84；总体拟合方程为y=0.997 90x-0.001 15，拟合直线的R2为0.999 36。测量值与真实值基本重合且分布均匀，具有较高的准确度和可靠性。

2.3.2年轮测量值误差分析 结合表3针对树木年轮自动识别软件测量值对比LinTab年轮分析仪参考值进行误差分析，其中， 对于梨树树种，共取样2棵，由于作为经济林，受人为营林措施干预较多，果树种植目的主要以产果为主，树木生长较为缓慢，经实验共获得227个宽度值。针对5个树种共取样37棵树，获取了782个年轮宽度值，分布范围为0.36～27.77 mm。其年轮宽度值偏差为-0.017 mm，相对偏差为-0.23%，均方根误差为0.139 mm，相对均方根误差为1.87%，满足年轮宽度信息提取基本需求。

表3 年轮测量宽度误差分析Table 3 Analysis of the annual ring width error

3 讨论

本文针对传统测树学中树木树芯获取、年轮信息获取费时费力的问题，并针对树木生长量测量研发了电动生长锥测树装置，通过研究传统生长锥钻取时受力状态，计算其所需扭矩峰值，设计研发电动生长锥，并对其树种适用性进行分析，完善树芯处理流程和电子档案存储，利用Python作为基础语言开发了一种可以快速进行年轮信息提取的软件脚本，该软件基于摄影测量原理、图像处理技术原理和Python函数库的有效利用，有效实现了对年轮信息的快速准确提取。其装置能全自动对树木任意部位进行树木树芯提取，对常见绿化树种如杨树、柳树、栾树等适用性良好，类似材质树种均可使用电动生长锥进行树芯提取，且其成本较低、效率高、省时省力。配合使用的树木年轮信息提取软件软件使用Python语言进行编写，操作简单、无需人为过多干预且无需额外成本，利用扫描仪、电脑即可完成测量，不受专业设备和环境的限制，适用于林业基层进行树木年轮信息获取，且精度较高，有一定的应用价值和推广潜力。

测试实验结果表明本装置能满足生长锥钻入退出需求，且对杨树、悬铃木、旱柳、栾树、梨树适用性良好，类似材质树种均可使用电动生长锥进行树芯提取；使用电动生长锥作业时间优于传统生长锥组，内业年轮提取工作流程自动化，无需过多人为干预，其效率优于LinTab年轮分析仪组，且实验后获得的树芯扫描影像可以永久保留，不受储存环境的影响和限制，方便后续继续研究。

目前，电动生长锥锥头仍依赖由瑞士公司生产的Haglof树木生长锥，随着生长锥国产化的实现，其整体成本会进一步下降；电动生长锥的电机、电池及控制器仍有进步空间，随着电动机、电池和控制器的发展，电动生长锥的重量和体积将会大幅度下降，便携性将进一步提高。本研究中树木年轮信息提取软件的实现仍是建立在获得树木树芯的基础上，而树木树芯获取具有一定的破坏性，对树木生长或多或少都会造成影响，随着近些年医疗设备的进步，便携式CT仪的出现将使无损树木年轮信息提取成为可能，通过对立木胸径层层切片，即可获得其年轮信息，但由于目前成本过高，在林业行业中应用可能存在困难。