木材径切面内部缺陷的应力波成像算法

郑泽宇， 冯海林， 杜晓晨， 方益明

（1.浙江农林大学 信息工程学院，浙江 杭州311300；2.浙江农林大学 浙江省林业智能监测与信息技术研究重点实验室，浙江 杭州311300；3.浙江农林大学林业感知技术与智能装备国家林业和草原局重点实验室，浙江 杭州311300）

木材的内部缺陷会严重影响木材的质量。由于应力波的木材无损检测成本低、携带方便及对木材无损的特性，使得其越来越受到人们的青睐。国内外已开展了应力波速度与木材内部缺陷之间的相关研究。如WANG［1］通过应力波在缺陷木材中的传播规律，表明应力波在腐烂或劣化的木材中比在完好的木材中传播速度慢。徐华东等［2］以校园内的活立木作为样本，使用应力波检测仪器对校园树木进行二维成像，表明在活立木的断层横截面上，应力波横向的波速要小于径向的波速，且应力波在木材内部的波速受制于木材内部的缺陷情况。国内学者也对应力波无损检测仪器的可行性进行了研究。杨学春等［3］利用应力波检测仪器对原木内部腐朽进行了研究，结果表明：应力波测试仪能准确判断不同树种原木内部的腐朽程度，并且能够得到原木内部腐朽基本形状的二维图像。LIN等［4］运用应力波断层成像技术，对不同大小的樟树Cinnamomum camphora人造空洞进行了研究，表明应力波断层成像技术可以反映木材的空洞大小及位置，验证了应力波无损成像技术的可行性。以上研究主要针对应力波传播规律与成像技术的应用领域，但对应力波成像算法的研究不多。FENG等［5］提出了一种基于插值的图像重建算法，利用周围点的值估计未知点的速度情况，进而重建树木内部缺陷图像。DU等［6］提出了基于椭圆的空间插值以及速度补偿的方法对木材的横截面进行成像，能够准确地判定缺陷位置，并得到了良好的成像效果。国内外对于应力波在树木内部的传播规律做了大量的研究［8-10］，但是对于木材径切面传播规律和成像方法研究还比较少。翁翔等［7］利用无损检测仪器对木材的径切面进行测量，对应力波横向速度和径向速度的比值进行了回归，得到了应力波在木材径向上的传播速度模型，其结果的拟合度较高。但在木材径切面的缺陷成像方法研究较少。鉴于此，本研究基于应力波在木材径切面中沿直线上传播的假设，提出一种木材径切面内部缺陷的速度修正成像方法。

1 应力波木材径切面成像原理

1.1 树木径切面模型的建立

木材可以用圆柱体来近似的表示，那么木材的径切面可以被抽象成一个近似的矩形模型，将n个（n一般为12）传感器放置于被测树木两侧（传感器1至n/2位于同一侧）。每个传感器距离地面的高度（或者某一水平高度的距离）分别为（H1，H2，H3， …，Hn）；S1～Sn表示应力波传感器，传感器按照从右向左，呈U字形进行放置；D表示树木直径；θjk为每个传感器两两之间的相对角度。如图1所示：定义一个坐标系表示树木径切面。

图1 木材径切面模型Figure 1 Radial and longitudinal plane model of wood

令Sn/2+1所在的位置为原点坐标，建立传感器的坐标系。

式（1）中：有1≤i≤n，Xi表示传感器位置横坐标，Yi表示传感器位置纵坐标。根据每个传感器的位置数据，可得到每个传感器两两之间的相对角度θjk为：

式（2）中： 1≤j≤n/2，n/2+1≤k≤n。

假定应力波在木材中是沿着直线传播的，通过传感器间的距离和应力波传播的时间计算出传感器间应力波的传播速度［11］。

1.2 速度修正插值方法

为了预测木材离轴的单轴压缩强度，Hankinson公式［12］给出一种数学模型。当应力波在木材径切面上传播时，如果传播方向与垂直于木材纹理角方向所成的角度为θ，那么根据Hankinson公式可以推出如下公式：

式（3）中：V1表示平行于纹理方向上的应力波传播速度；Vr表示垂直于纹理方向上的应力波传播速度；θ表示传播方向与木材纹理方向所形成的角度；Vθ表示在角度为θ时应力波在木材径切面上的传播速度。

翁翔等［7］通过计算得出：

根据式（5），可以得出:

就可以得出速度修正的公式：

式（7）中：Vc（θjk）表示修正后的速度；Vd（θjk）表示为测量得到的应力波传播速度；k是常数；θjk表示2个传感器间与垂直于木材方向纹理所成的角度。

因此经过缺陷区域的Vα经过速度修正后，就可以反映出应力波经过缺陷区域速度呈现出降低的现象，有效地减少了由木材的各向异性［13］导致的应力波在木材不同方向上传播速度不一致现象的影响，为木材内部缺陷成像提供更精确的数据。

将通过传感器之间的线速度VL转换成木材径切面上交点的速度VP：

式（8）中：VL1至VLn分别表示每条经过该交点的应力波传播路径上修正后的线速度；（VL1,VL2,…,VLn）∈Vc（θjk）。

图2 应力波在不同缺陷程度的木材中传播Figure 2 Stress waves propagate in wood of different defect degree

利用速度修正插值法对木材的径切面进行成像。每个已知点都对插值点具有一定的影响，即权重［14］。权重随着已知点和插值点之间距离的增加而减弱，距离插值点越近的已知点的权重越大。而且当已知点在距离插值点一定范围以外时，权重可以忽略不计。在任一待插值点的值是邻域内已知点权重之和。可以表示为：

式（9）中：Vqi为待插值点的预估值；V′Ln为每条经过该交点的应力波传播路径上未修正的线速度，（V′L1,V′L2,…,V′Ln）∈Vd（θjk）；θLn表示每条经过该交点的应力波传播路径与水平面的所成的夹角， （θL1,θL2,…,θLn）∈θjk；di为第i个已知点与待插值点的距离；n表示邻域内参与插值计算的已知点的个数；m为常数。

成像算法步骤：

步骤1：初始化（建立木材径切面的模型，输入传感器的坐标（x,y），计算传感器间的线速度Vd（θjk）；

步骤2： 利用式（7）求出修正后的线速度矩阵Vc（θjk）；

步骤3：传感器两两连线的交点坐标（xc,yc），将邻域内已知属性值的交点的坐标以及每个交点对应的线速度V′Ln代入式（8）计算Vp的属性值；

步骤4：通过将已知点的属性值代入式（9），计算未知点的属性值；

步骤5：重复步骤4，直到计算出所有径切面上的每个点的属性值。

步骤6：将预估点根据属性值的不同进行不同颜色的赋值，并生成二维图像。

2 实验材料与方法

实验用于采集木材应力波数据的仪器是Fakopp。如图3所示：将传感器分别放置在原木径切面两边，每边6个传感器，总共12个。实验过程中使用检测仪器自带的重锤以相同的力度，垂直于传感器方向，对各个传感器分别敲击3～5次，采集总共36～60组数据。将数据用于下文的成像中。

图3 应力波测量装置与传感器的布置Figure 3 Stress wave measuring device and sensor arrangement

本实验样本编号、样本种类、样本周长、样本测量高度等信息如表1所示。

对样本的空洞面积测量方法如下：利用卷尺测量树木空洞的长和宽，假设样本的长为L，宽为W，则对于长方形的空洞缺陷面积SR计算公式为：

假设梯形的空洞，上底长度为A，下底长度为B，缺陷的高度为H，梯形的空洞的面积ST则用公式可以表示为：

表1 实验样本基本信息Table 1 Experimental sample information

假设圆形的缺陷孔洞的直径为R，那么缺陷的面积SY可以表示为：

3 结果与分析

3.1 速度修正插值法（VCI）与反距离加权插值法（IDW）成像结果对比

从图4A可以看出：缺陷为人工挖掘并且形状为长方形。如图4B所示：未修正的线图将经过缺陷区域的线段表示为黄色，未正确区分缺陷区域。图4C表示的是修正后的线图，将经过缺陷区域的线段（如1号和7号、2号和8号等传感器之间的连线）错误地表示为缺陷颜色红色。从图4D可以看出：使用IDW的成像结果，只能粗略反映缺陷位置，对于缺陷大小与形状都不能进行准确成像。对比VCI的成像结果图4E，VCI可以较好还原缺陷位置，可以看出缺陷的位置是处于第3号、5号与8号、10号传感器之间，与图4A所示的真实缺陷位置基本吻合。但是对于缺陷形状的成像不够准确。

图4 雪松样本实验结果Figure 4 Sample experiment results of Cedrus deodara sample

从图5A可以看出：样本缺陷位置处于右边2号至5号传感器的下面3 cm，以及左边11号传感器至7号传感器的下面3 cm的位置，形状为矩形。从未修正线图（图5B）可以看到：仅有少部分经过缺陷区域的线段被表示为缺陷的红色，大部分经过缺陷区域的线段被表示为非缺陷的绿色。修正后的线图（图5C）较好地区分缺陷区域，并以不同的颜色进行区分。可以看到IDW的成像结果虽然能够还原缺陷的位置（图5D），但对缺陷还原的准确度较低，对缺陷位置和大小的还原参考价值较低。如图5E所示：VCI的成像结果可以较好地反映缺陷位置，但存在对缺陷区域成像结果比真实结果偏小的情况，以及对缺陷边缘的重建不够平整，对于样本的缺陷轮廓可以较好还原，可以反映样本的缺陷特征。

图5 泡桐样本实验结果Figure 5 Sample experiment results of Paulownia fortunei sample

从图6A可以看出：树木的缺陷位置位于右边的3号和4号传感器之间，以及左边的9号和10号传感器之间。对比图6B与图6C可以看到：未修正的线图将经过缺陷区域的线段（如7号、1号传感器线段以及6号、12号传感器线段）错误地表示为经过正常区域的线段，呈现绿色，与真实情况不符。而修正后的线图，能够较好区分缺陷部位。从IDW成像结果（图6D）反映缺陷所处的低速度区域，但对于缺陷的形状无法正确成像。缺陷成像结果呈现不规则的块状分布，无法反映样本的缺陷特征。VCI的成像结果（图6E）显示：缺陷形状呈现一个类似圆形的不规则形状，与真实的缺陷形状较为接近。成像结果表明缺陷位于3号、4号传感器与9号、10号传感器之间，与图6A的真实缺陷位置也基本相同。

图6 樟树样本实验结果Figure 6 Sample experiment results of Cinnamomum camphora sample

表2 缺陷与预测的4种组合Table 2 Four combination of real defects and predict defects

3.2 成像结果误差分析

使用混淆矩阵法对成像效果进行定量的误差分析，混淆矩阵包含由分类系统完成的实际和预测的分类信息［15］。 FAWCETT［16］描述了真、 假分类系统中 4 种可能组合（表 2）。

从表2的矩阵可以推出的常用度量值［17］，如准确率或者正确率（A），精确度（P）和查全率（R）。

式（13）～（15）中：准确率A表示对木材状态预测正确的百分比；精确度P表示对木材缺陷预测的可靠性水平高低；查全率R表示对木材缺陷的预测能力高低；C，W，E和F所代表含义见表2。

从表3可见：3个样本的准确率A为72.32%～95.85%。使用VCI方法的3个样本的平均准确率为93.34%，使用IDW方法的3个样本平均准确率为82.63%。IDW方法在1号、2号样本中的准确率均大于85.00%，其原因应该是IDW方法在成像结果中虽然不能很好地反映木材缺陷的特征，但是对木材完好的部位成像结果较好，所以得到较高的准确率。而3号样本只有72.32%的准确率，对比真实缺陷可以看出：IDW方法将部分正常木材错误地成像为缺陷木材，所以导致准确率下降。IDW和VCI的准确度标准差分别为7.38%和2.04%。可见，VCI在准确率上的标准差较小，说明VCI方法在不同样本中的准确率结果更稳定，准确率相较于IDW方法均有所提高。

表3 不同样本的准确率Table 3 Accuracy of different samples

表4 不同样本的精确度Table 4 Precision of different samples

从表4可见：3个样本的精确度P为65.79%～88.35%。使用VCI方法的3个样本的平均精确度为82.26%，使用IDW方法的3个样本平均精确度为70.19%。VCI方法中3号样本的结果要好于其他样本，其原因是3号样本中的预测缺陷大部分与真实缺陷重合，且预测缺陷大部分落在真实缺陷上，所以得到较高的精确度。IDW和VCI方法的精确度标准差分别为5.00%和4.70%。VCI方法在准确度上的标准差较小，VCI方法在精确度的稳定性略高于IDW方法。同时VCI方法在不同样本上的精确度均好于IDW方法。

从表5可见：3个样本的查全率R为17.98%～94.03%。使用VCI方法的3个样本的平均查全率为92.65%，使用IDW方法的3个样本平均查全率为36.64%。VCI方法中的查全率高于IDW方法。可见，VCI方法可以较好地反映木材的真实缺陷结果，且对木材的缺陷位置成像结果更好。在真实缺陷大小不变的情况下，由于VCI方法可以更好地反映木材的缺陷情况，因此得到了更高的查全率，与真实的成像结果相吻合。IDW和VCI方法的查全率标准差分别为：15.59%和1.93%。说明VCI方法的查全率标准差要远小于IDW方法，在查全率的稳定性上要好于IDW方法。

表5 不同样本的查全率Table 5 Recall rate of different samples

4 结论

本研究提出了一种应用于木材无损检测领域的，基于应力波的木材径切面成像方法。通过与IDW的成像效果比较发现：不论是在定量还是定性方面，VCI方法对木材径切面的缺陷成像结果都好于IDW方法。与真实的缺陷进行对比得到较小的误差，可以准确反映木材内部的真实缺陷情况，具有较高的成像精度。利用传感器测量得到的不同应力波数据，通过对不同缺陷的样本进行实验，验证了该算法对于木材的径切面成像效果的可行性。利用混淆矩阵对成像结果进行了定量分析发现：VCI方法的平均准确度为93.34%，平均精确度为82.26%，平均查全率为92.65%，均高于IDW方法，但对于缺陷具体形状的成像结果，其准确度有待进一步提高。