袁 霄 陈勇平,2 唐启恒,2 郭文静,2
(1.中国林业科学研究院木材工业研究所 北京 100091; 2.中国林业科学研究院林业新技术研究所 北京 100091)
木材树种识别是古建筑木构件无损检测和安全评估的重要内容(姜笑梅等, 2010; 李鑫, 2017),因为树种不同,材料的物理力学性能不同,其健康评价的基准也不同,若不明确树种信息可能产生一定误判; 同时,木材树种识别也是木结构古建筑维修和保护的基础性工作,通过树种识别获取木材名称和材质参数,可为制定木结构古建筑保护方案及更换木构件提供选材依据,且树种识别结果也可作为档案信息或为科学研究服务(晋宏逵等, 2007; 曹静等, 2019; 李世民等, 2019)。
选取落叶松(Larixsp.)、杉木(Cunninghamiasp.)、硬木松(Pinussp.)、软木松(Pinussp.)、杨木(Populussp.)5种木结构古建筑常用木材树种以及栎木(Quercussp.)、核桃木(Juglanssp.)、水曲柳(Fraxinussp.)3种补充木材树种为研究对象。将以上树种木材加工成100 cm×50 cm×50 cm试样,每树种试样8块,以落叶松为例编号依次为落叶松-No1~落叶松-No8; 采集以上树种且不同批的验证木材,同样加工成100 cm×50 cm×50 cm试样,每树种试样2块,以落叶松为例编号为落叶松-Y1、落叶松-Y2,试样制备示意如图1。选用木材均为气干材,为与实际应用接轨,试样含水率使用便携式含水率测定仪进行检测,落叶松、杉木、硬木松、软木松、杨木、栎木、核桃木和水曲柳试样含水率均在8%~12%之间。
图1 试样制备Fig. 1 Sample preparation
利用Resistograph 4452-P阻抗仪检测试样,检测时微型探针沿木材年轮法线从边材往心材方向探测,获取相应的微损检测曲线。根据现场实际检测情况,为排除一些干扰因素,检测曲线前后各去除5 cm,即检测曲线有效分析范围为探测深度的5~45 cm区间,如图2。
图2 数据采集Fig. 2 Data collection
描述数据曲线主要有2点: 一是曲线数值; 二是曲线走势。
1.3.1 曲线数值量化 曲线数值采用整体均值和曲线占比2种方式量化。整体均值指对应横坐标区间每隔1 cm读取数值100个,5~45 cm区间共读取数值4 000个,计算各对应数值求取平均值,并完善数据和曲线编号,以落叶松为例,编号落叶松-No1修改为落叶松-阻抗仪检测平均值-No1。曲线占比指以颜色直方图方式评估曲线在整幅图像中的占比。颜色直方图在许多图像检索系统中被广泛采用,其描述的是不同色彩在整幅图像中所占比例,并不关心每种色彩所处空间位置。曲线占比一定程度上反映峰谷变化幅度,对所有检测曲线来说画布大小不变,即整幅图像像素点是一定的,而曲线可分成多个像素点,同等条件下峰谷变化越大像素点的接长越长,该颜色占比越高。考虑到颜色直方图基于不同的颜色空间,而本研究检索色彩较少,故采用Gray颜色直方图。
1.3.2 曲线走势量化 曲线走势一定程度上反映图像的纹理特征,故本研究曲线走势采用纹理特征方式量化。纹理特征量化包括特征点选取、特征点主方向确定和特征点描述3个步骤。
1) 特征点选取: 在灰度图像下建立一个Hessian二阶偏导矩阵,求一定范围内像素点的极值,选用二阶标准高斯函数作为滤波器。高斯过滤是一种低通滤波,可除去图像细节而保持整体不变。
建立的Hessian矩阵如下:
H(x,y,σ)=
对该矩阵进行运算:
若det(H)>0,则为极值点,对Hessian矩阵检测到的极值设立一个阈值,若det(H)小于该阈值,则不计算; 若det(H)大于该阈值且比邻近26个点响应值均大,则设为特征点。
2) 特征点主方向确定: 在以特征点为中心、以6s(s为特征点的尺度)为半径的圆形区域内,对图像进行Haar小波(边长4s)响应运算。具体步骤为: 在以 6s为半径的特征点邻域范围内,计算每个像素点x和y方向上的Haar小波响应值,并将其记为hx和hy; 以特征点为中心,对以 6s为半径范围内的像素点进行σ= 2s高斯加权,加权后的值分别表示在水平和垂直方向上的Haar小波响应值,记为Whx和Why; 采用大小为π/3的滑动窗口,计算窗口内水平和垂直方向的响应值总和,记为∑Whx和∑Why; 以步长5°左右旋转滑动窗口,依次计算每个π/3窗口中的∑Whx和∑Why,找到最长矢量作为该特征点的方向。
3) 特征点描述: 以特征点为中心,沿主方向将20s×20s的图像划分成4×4个子区域,每个子区域利用尺寸2s的Haar小波模板进行响应计算,并统计响应值形成特征矢量。具体步骤为: 分别在 16 个子区域中计算每个像素点x和y方向上的Haar小波响应值并进行σ=3.3s的高斯加权平均,计为dx和dy; 将每个子区域中Haar小波响应值dx和dy分别相加,得到每个区域的响应值之和,记为∑dx和∑dy; 将每个子区域中Haar小波响应值dx和dy的绝对值相加,得到每个区域所有像素点的响应值绝对值之和,记为 ∑|dx|和∑|dy|; 将每个区域中的∑dx、∑dy、∑|dx|、∑|dy|组合成一个四维向量,16个子区域就得到一个 16×4=64 的特征向量,即为描述算子。
至此,整体均值、曲线占比和曲线走势量化完成,形成无损检测曲线基础数据库如图3。
图3 基础数据库Fig. 3 Basic database
因已对曲线进行量化,故曲线之间的比较等同于2组数字之间的比较。
1) 曲线占比相似度算法: 假设Q代表待检图像,D代表数据库图像,HQ(k)和HD(k)分别为2幅图像颜色直方图,则2幅图像之间的颜色匹配值P(Q,D)可借助相关系数算法来实现,即:
2) 曲线走势相似度算法: 曲线走势识别即找到曲线中具有典型代表性的局部特征,对该区域进行描述,并建立描述算子,通过描述算子进行相似度搜索。纹理特征比对算法经比较确定为Surf算法。Surf特征识别算法主要先找出特征点建立描述算子,后通过描述算子进行相似度比较。2个局部纹理描述算子的相似性度量,本研究采用欧式距离进行计算:
式中:Xik表示待检图像中第i个特征描述算子的第k个元素;Xjk表示标准图像中第j个特征描述算子的第k个元素;n表示特征向量维数。
对待检图像上的特征点,计算其到标准图像上所有特征点的欧氏距离,得到一个距离集合,对距离集合进行比较运算得到小欧氏距离和次小欧式距离。设定一个阈值0.7,当最小欧氏距离和次小欧式距离的比值小于该阈值时,认为特征点与对应最小欧氏距离的特征点是匹配的,否则没有点与该特征点匹配。逐个度量待检图像中的描述算子,最后用匹配点数除以描述算子总数即曲线走势相似度。为将繁琐的计算比对过程转化为直观的界面,本研究采用自主研发的ImageCV图像识别软件,该软件主要是以1.2中量化后的标准曲线信息为基础数据,植入1.3中曲线数值和曲线走势相似度算法进行待识别曲线的检索和识别。
利用阻抗仪微损检测曲线进行木材树种识别不同于以往木材无损鉴定(刘红清等, 2014; 李敏华等, 2012),该方法主要基于数理统计和相似度分析理论,其树种识别手段为固定区域范围内的选择性比对鉴定方法。识别原理可以简述为: 通过建立木材微损检测曲线中整体均值、曲线占比、区域走势等关键参数与古建筑木构件常用树种信息之间的映射关系,构建古建筑木构件常用树种微损检测信息数据库,而后采集待识别木构件的微损检测曲线,并将该微损检测曲线与微损检测信息数据库中的各微损检测曲线进行比对获得树种识别结果。
不同树种具有不同的阻抗仪检测曲线形式,尤其当木材密度不同时表现得更为明显。一般来说,当树种不同且密度相差较大时,曲线数值和走势差异明显,如图4。
图4 密度差异显著树种的阻抗仪检测曲线Fig. 4 Resistograph testing curves of wood species with significant density differences
这是因为密度与阻抗仪检测值呈正相关,木材密度越大,其实体物质含量越多,微型探针刺入木材内部消耗的能量越大,阻抗仪检测值越高。图4中,栎木木材实测密度约为0.86 g·cm-3,落叶松木材实测密度约为0.61 g·cm-3。
树种不同但密度接近的木材,其阻抗仪检测曲线也有一定差异,如图5。
图5 密度差异不显著树种的阻抗仪检测曲线Fig. 5 Resistograph testing curves of wood species with similar density
由图5可知,虽然二者密度接近(硬木松木材密度约0.46 g·cm-3,软木松木材密度约0.42 g·cm-3),但其阻抗仪检测曲线形式有一定差异,可能是因为除密度外木材软硬等其他材质特性也会影响阻抗仪检测值,这一点从栎木和落叶松的阻抗仪检测值差值也可以看出。
为提高识别精度,本研究采用2种识别判定方法: 第1种识别模式,与待识别曲线比对,相似度最高曲线对应的树种为树种识别结果; 第2种识别模式,与待识别曲线比对,相似度较高8条曲线对应的树种为树种识别备选结果,备选结果中树种占比最高的为树种识别结果,若备选结果中树种占比最高的出现2组或2组以上,则对8条曲线按相似度高低排序,排序号相加最低的为树种识别结果。
从表1可以看出,无论何种识别模式,用于验证的8组16个样本曲线均可被准确识别; 但这并不能说明该识别方法的准确率为100%,首先数据库样本量目前并不足以涵盖所有常用树种以及各树种所有的曲线组合,其次验证样本还需要根据实际情况进行不断补充和完善。表1可以证明,利用阻抗仪微损检测曲线进行树种辅助鉴定的思路是可行的: 从小处看,一个建筑单体用材树种是有限的,取样产地也是有限的,通过前期建立的微损检测曲线数据库对待识别曲线进行比对分析得出鉴定结果是科学可行的; 从大处看,随着数据库逐渐完善,一个地区乃至全国木结构建筑就如同一个建筑单体一样,其比对鉴定原理是相同的。当然,每次应用时根据前期调研结果和部分样品分析,限定比对的适宜树种,可缩小选择范围,提高识别精度; 而且,持续对部分未知树种构件的比对结果进行传统方法验证,也可以完善该辅助技术。
表1 树种识别结果Tab.1 Identification results of wood species
由表1同时可以看出,采用第2种模式进行识别时,相似度较高8条曲线对应的树种并不是单一树种(图6),这一定程度上说明可能会存在误判风险。为了更好地分析,以落叶松-35-Y1为例,分别就整体均值、曲线占比和曲线走势进行详细比对说明。
从图6可以看出,落叶松-35-Y1虽然被准确识别为落叶松,但在颜色相似度比较中,与其他树种并无太大区别,如与落叶松-35-No3颜色相似度为99.99%,与杨木-29-No7颜色相似度为99.35%。这主要是因为检测曲线所在图像中背景色占比均超过90%,曲线占比差异在整体图像中比率较小(图7),故颜色相似度较为接近仅表现为略有差异。
图7 图像背景分析Fig. 7 Background color analysis of image
从图7可以看出,落叶松-35-No3检测曲线所在图像中背景色占比为92.99%,即曲线和辅助坐标占比为7.01%; 杨木-29-No7背景色占比为94.13%,即曲线和辅助坐标占比为5.87%。进一步去除辅助坐标的影响,落叶松-35-No3曲线占比为2.41%,杨木-29-No7曲线占比为1.18%,如图8,2树种曲线占比存在明显差异,该过程相当于曲线拉伸,后续应对颜色直方图算法进行修正。
图8 曲线占比分析Fig. 8 Curves ratio analysis of image
从图6也可以看出,落叶松-35-Y1与水曲柳-128-No5检测曲线颜色相似度接近,这是因为水曲柳-128-No5背景色占比为92.48%,与落叶松-35-No3背景色占比极其接近,但水曲柳在初始的整体均值比较中与落叶松差异较大,如图9所示(图中3条曲线为示意的曲线拉伸效果,并非实际拉伸长度),3树种曲线在实际分析中存在较明显差异。
图9 曲线拉伸Fig. 9 Curves stretch
值得一提的是,硬木松-13-Y1和硬木松-14-Y2虽然均识别为硬木松,但在8条相似度较高曲线中仅有3条识别为硬木松,如硬木松-14-Y2识别中出现杉木、核桃木等树种(图10),不可排除多样品验证时可能出现极大概率的误判。为此,本研究就硬木松木材增加4条验证曲线进行补充分析,其识别比例也多为3~4之间。究其原因主要是上述几种木材现有的阻抗仪检测值略接近且均较低,导致从整体均值、曲线占比和纹理特征难以明显区分。
图10 硬木松识别结果Fig. 10 Identification results of pine
对于此类情况,应对检测曲线的获取或在后期处理上进行适当调整,扩大差异,如图11所示为后期调整的曲线。
图11 曲线调整Fig. 11 Curves adjustment
从图11可以看出,后期调整的曲线走势略为明显,据此本研究重建硬木松、杉木、核桃木和软木松等密度相近木材的树种数据库,并以硬木松-14-Y2为例进行再次对比,但重建后数据库比对效果并无明显改观。单纯依靠后期处理无法解决问题,曲线调整应结合实际检测过程进行,实际检测过程中若部分曲线数值过低,应对仪器检测速度进行实时调整,以增大探针能耗,相对提高阻抗仪检测值,从而得出特征更为清晰的曲线。当然,曲线调整并非是调整个别树种问题,需要综合考虑整体数据库的基准或建立多基准数据库; 此外,对于某些树种是否适用于此类方法也值得进一步探究。
综上所述,利用微损检测曲线进行古建筑木构件用材树种识别是可行的,且不论对于一个特定的木结构古建筑来说,其用材树种是相当有限的。若木材密度相差较大,可分析曲线高低直接区分; 若木材密度接近但曲线峰谷变化差异大,可根据曲线占比和曲线走势区分; 若木材密度和曲线占比均接近,需综合古建筑历史记录、木材取材信息以及其他特征加以确定。实际应用中,结合传统取样鉴定结果并将取样构件的微损检测曲线增列至数据库,可提高识别精度。
1) 以多基准的木材微损检测曲线数据库为基础,通过待识别曲线的比对分析,可实现木结构古建筑常用落叶松、杉木、松木、杨木等木材树种的初步鉴定。
2) 古建筑木构件用材树种现场识别是一个长期性工作,需不断应用-反馈-改进,未必所有用材树种均适合此类分析方法,但随着基础数据库、比对算法、适宜树种、微损曲线特征的不断完善,其鉴定可靠性将得以提高。
3) 古建筑木构件用材树种现场识别是一种区域范围内的选择性鉴定方法,是针对文物建筑特点的一种适宜性保护识别技术,该技术无法取代传统树种微观鉴定方法,但可作为传统鉴定手段的有效补充。