梁 臣， 冯海林， 方益明， 李 剑(1.浙江农林大学 信息工程学院，杭州 311300；2.浙江省林业智能监测与信息技术研究重点实验室，杭州 311300)

无损检测作为一种直接、准确、客观的方法,满足木材生产中非破坏性快速检测和持续检测的需要而越来越受重视。国内外已在此方面做了大量的研究工作,在将超声波、射线、微波等无损检测技术应用到木材物理性质、生长特性及木材缺陷等方面取得了一定的成绩,研究出多种无损检测手段,并通过实验取得了较好的成果[1-3]。但射线、微波检测方法所需设备价格昂贵；相同条件下，超声波在原木中的衰减比应力波更快[4]，所以应力波在原木中的适应性更强。

目前，国内外在使用应力波对树木内部的缺陷已经进行了较广泛的研究。梁善庆等[5]通过应力波传感器波速矩阵重构所得断层二维图像能够直观显示缺陷部位、大小及形状等信息；王立海等[6]发现使用应力波断层成像技术时，若需对原木缺陷进行精确测量，至少需12个传感器才能满足要求；仅需要判断原木是否存在缺陷时，选用6个传感器就能满足要求。余观夏等[7]运用应力波技术分析了原木横截面中的频谱曲线，通过频谱曲线下的面积和一定频率范围的面积可计算木质材料的好坏。

冲击回波法是20世纪80年代中期发展兴起的一种无损检测方法，通过分析应力波在物体内部的反射情况，能根据应力波信号特征检测物体内部缺陷，具有适应性强、反应快、操作简便、造价低等特点。该方法在混凝土、桩基缺陷的检测当中应用比较广泛[8-9]，而在原木检测中的应用较少。

原木中的应力波信号有表面波、纵波、横波等成分，但是各成分之间特征差异较大[10]。本文在现有的应力波检测基础上提出一种基于原木应力波信号分解的冲击回波法。首先将应力波信号分解成多个平稳的、幅度和频率受调制的固有模态函数(IMF)分量，然后找出对空洞敏感的IMF，最后根据该IMF的频域特征实现原木内部空洞径向缺陷定位。

1 基于原木应力波信号分解的冲击回波法

1.1 冲击回波法

应力波在整个原木轴向的传播过程可用一维杆进行分析[11]，原理如图1，小锤敲击物体截面时产生纵波(P 波)、横波(S波)、表面波(R波)。P波又称Primary wave，传播方向与质点震动方向一致，即与小锤的敲击方向一致，P波是物体内传播最快的一种应力波。S波作为一种横波，特点是质点的振动方向与它的传播方向垂直，因此S波的传播时候会形成波峰和波谷。S波的传播速度次于P波。R波是一种在物体表面传播的应力波，R波能量较大，但衰减较快。P波在敲击面和空洞之间来回反射形成瞬态共振，通过传感器获得敲击点附近的振动信号，对振动信号进行FFT运算后可得共振的频率，若知应力波在物体中的传播速度，则空洞与敲击端的距离为

(1)

式中：CP为P波传播的速度；f为P波的峰值频率；b为修正系数，取值为0.96[12]。通过此式可对缺陷位置做出预测。

图1 冲击回波法Fig.1 Impact-echo method

1.2 经验模态分解

冲击回波信号可能包含强烈的表面波,这是一个宽频率范围的间歇信号。然而，经典的傅里叶变换方法只能处理线性和平稳信号，包括应力波在内的许多实际采样信号(模型)是非平稳的或非线性的。短时傅里叶变换、小波变换等分析方法均采用固定的核函数，对信号缺乏自适应性，导致分析结果往往失去实际物理意义而不能正确提取信号本质特征[13-14]。

经验模态分解(Empirical Mode Decomposition，EMD)适合于非线性、非平稳信号分析，其本质是对信号进行平稳化和线性化处理，把复杂的信号分解成有限个瞬时频率有意义、幅度或频率受调制的本征模态分量(Intrinsic Mode Function，IMF)和一个残余函数，即

(2)

式中：r(t)为残余函数，它是一单调函数，Ci(t)为各个IMF分量，并满足以下两个条件：① 在整个数据段，极值点的个数和零交叉点的个数必须相等或相差不能超过一个;② 在任何一点，由局部极大值点形成的包络线和局部极小值形成的包络线的平均值为零。

1.3 基于原木应力波信号分解的冲击回波法

当应力波信号的时间尺度存在跳跃性变化时,对信号进行EMD分解,会出现一个IMF分量包含不同时间尺度特征成分的情况。在分解时为了解决这种情况，法借鉴集合经验模态分解(EEMD)，本文提出基于原木应力波信号分解的冲击回波(流程如图2)步骤如下:

图2 基于原木应力波信号分解的冲击回波法Fig.2 Impact-echo method based on stress wave signal

步骤1 向采集的原木冲击回波信号S(t)中加入正态分布白噪声得到带有白噪声的信号x(t)。

步骤2 确定信号x(t)的所有局部极值点，通过三次样条函数求取其上包络u1(t)和下包络u2(t)的局部均值

(3)

步骤3 令h(t)=x(t)-m(t)，如h(t)不满足IMF条件，则视其为新的x(t)。重复k次得

h1k(t)=h1(k-1)(t)-m1k(t)

(4)

经过n次迭代满足1.2节中两个条件得到的hn(t)即为有效IMF,剩余信号则进入下一轮筛选过程。

步骤4 重复步骤1～步骤3每次加入新的白噪声序列，将每次得到的IMF集成均值作为最终分解结果。

步骤5 剔除IMF集合中的白噪声和表面波。

步骤6 计算各个IMF与S之间的互相关系数R

(5)

RSy(τ)为x和y信号的互相关系数，y为步骤4得到的IMF，T为信号长度，τ为时间间隔。剔除相关系数低于0.3的IMF分量。

步骤7 对剩下的IMF进行傅里叶变换，根据IMF的频域特征结合式(1)计算缺陷位置。

1.4 基于时域特征和分量相关系数的IMF选择方法

由于算法中应力波信号分解加入了高斯噪声，且应力波信号中包含表面波，所以原信号分解后会有一个或多个IMF与噪声或表面波信号特征相似，可从IMF的时域信号特征识别噪声或表面波。

EMD分解不可避免的会出现过多分解现象，即得到的分量个数比原信号实际组成分量多，也即出现所谓的虚假分量。一般来说，信号越复杂，虚假分量的个数可能会越多，对后续信号的时频谱分析和进一步的处理越会产生不利影响。原则上，经EMD得到的各个模式分量ci(t)应分别对应于原信号各实际组分，然而由于前述原因，两者之间总存在误差，其差值即形成虚假分量。分析可知，各分量虽不等同实际组分，但必是后者的近似，即两者存在一定的相关性;而虚假分量由两者差值形成，相关性必然很小。因此，从各分量与原信号相关分析中应能分辨出其真伪。

2 噪声参数选择

加入白噪声的幅度和数量是分解的两个关键参数。通过实验观察，加入足够多的白噪声集合数量能使分解以后信号的噪声变得缓和，并建议在大多数情况下添加噪声的振幅与原信号的标准差比值为0.2[15-17]，但研究内容专注于从振幅相对较小的间歇信号中提取高频与低频连续信号。在冲击回波信号中包含强连续性的纵波信号和能量较强的间歇性表面波，这和上述信号有着不同的特征，也因为表面波拥有的较宽的频率带宽导致纵波的提取比较困难。

为了研究噪声参数设置在冲击回波应用中的效果, 利用有限元建立木材模型，模拟一个信号来表示冲击回波中原木表面质点的平面外位移测试。模拟应力波在木材中的传播已有一些研究，如利用有限元法模拟木材材料性质和缺陷对应力波传播的影响[18]，应力波在木材横切面的传播过程[19]等。本文则假设原木纹理与原木轴向平行，将原木看作一种正交各向异性材料。在有限元仿真软件ANSYS中选取正交各向异性单元Solid185作为建模单元，以松木为建模原型建立三维原木模型。冲击源为一个半径0.02 m的铁球，铁球以10 m/s的速度纵向冲击模型A端，冲击点离模型A端圆心0.07 m，在径向离冲击点0.14 m处采集质点加速度模拟纵波信号，采样频率0.1 MHz，冲击点和信号采集位置如图3,采集的信号如图4(a)。[20]

图3 信号采集和冲击位置Fig.3 Signal acquisition position and impact position

基于有限元分析中冲击回波模型[21]选取表面波函数

y(t)=Aλ(sin(πx/η))3(u(t)-u(t-η))

(6)

A为表面波与纵波的振幅比例，λ为纵波的幅值，x为纵波信号，u(t)为赫维赛德阶跃函数，η为冲击持续时间，取1 ms(如图4(b))。纵波与表面波混合信号如图4(c)。

(a) 纵波信号

(b) 表面波信号

(c) 混合信号图4 有限元仿真模拟原木冲击回波信号Fig.4 The log impact echo signal by finite element simulate

使用分解以后的IMF与纵波的相关系数作为衡量标准。人工设定噪声幅度与冲击回波标准差比例0，0.1，0.2，0.5，0.7，1，1.2，1.5，2，2.5和3；噪声数量从10～200，每次递增10。

添加不同的噪声幅度对相关系数的影响如图5(a)所示，当A=1时，添加较小的噪声比例(0.1)能让分解后的IMF与原冲击回波信号相关系数达到较高的程度(0.95)。当表面波的幅度增加，A=3，5，8时，添加较高的白噪声幅度(与冲击回波信号标准差比例为2)也能使相关系数达到较好的效果，分别为0.937，0.910，0.895。

尽管添加比较多的噪声数量能够影响EEMD分解的效果，但是添加太多数量的噪声也耗费较多的时间。不同的噪声数量对EEMD方法分解效果影响如图5(b)所示，当A=1时，添加噪声数量从20以后，相关系数开始趋向收敛。当A=3，5，8时，添加噪声数量从80以后，相关系数值也趋向稳定。

结果表明，即使在有较强表面波和噪声的情况下依旧能分成很多模态，设置合适的EEMD参数能够使得到的IMF与纵波信号接近。

(a) 噪声幅度与冲击回波标准差比例

(b) 噪声数量

图5 添加不同的噪声幅度和数量对EEMD分解后的IMF与P波之间相关系数的影响

Fig.5 Adding up the amplitude and number of different noises sees the effects on the correlation coefficient between IMF and longitudinal wave

3 试验分析

3.1 试验设备与材料

选取两棵雪松和一棵白杨健康原木，尺寸分别为1#雪松1.9 m，含水率36%，大小端直径分别为19 cm和22 cm；2#白杨2.32 m，含水率42%，大小端直径分别为20 cm和24 cm；3#雪松2 m，含水率40.5%，大小端直径分别为23 cm和26 cm。使用FAKOPP Microsecond Timer微秒计测出两端之间的传播时间，根据原木长度计算出波速。将原木分AB两端置于水平地面，距离1号原木A端1 m、2号A端0.7 m、3号A端0.3 m处凿开空洞，空洞尺寸从0.08×0.08 m、0.1×0.1 m、0.12×0.12 m逐渐扩大。用铁锤敲击原木A端，敲击过程中把握力度保证每次敲击力度差别不大，敲击点离原木髓心0.07 m，在敲击点径向0.14 m处垂直钉上加速度传感器，传感器型号为兰德BZ1105，示波器采样频率为0.25 MHz，B端按同样的方式测试。实验平台、原木模型尺寸、空洞大小与位置分别如图6、图7、图8所示。

图6 试验平台示意图Fig.6 Test platform

图7 原木及各空洞尺寸(cm)Fig.7 Logs and different size holes(cm)

图8 原木模型及空洞位置示意图Fig.8 Models of logs and location of holes

3.2 试验结果分析

分别对1～3号原木两端施加冲击，采集不同尺寸空洞的测试信号，1#原木空洞尺寸0.08 m×0.08 m时，冲击回波信号如图9(a)。对信号进行分解，添加噪声的幅度标准差比例为2，噪声数量为100，分解结果如图9(b)。从图中观察可看出IMF1～2明显为频率较高的白噪声。IMF3幅值较高的波峰均在0.2 ms以内，且衰减较快，与表面波信号特征相似；IMF10为残余分量，在分析时应以剔除。

(a) 冲击回波采集信号

(b) 冲击回波信号分解结果图9 1#原木0.08 m×0.08 m空洞尺寸下A端应力波信号

Fig.9 The 0.08 m×0.08 m hole size log Stress wave signals of No.1 log model in endA

通过各IMF与原信号互相关系数表中(表1)可看出IMF6～9与原信号相关度都在0.3以下，属于虚假分量。IMF4和IMF5与原信号相关性较高，并拥有连续性纵波特征，信息可用于计算。对IMF4和IMF5进行傅里叶变换，结果如图10，然后通过IMF4频域谱峰值2 040 Hz和IMF5频域谱峰值1 120 Hz代入式(1)，计算结果分别为0.957 m对应缺陷位置误差4.3%和1.855 m对应原木末端位置误差2.4%。

图10 1#原木基于时域特征和分量相关系数相结合选择出IMF频域信号

Fig.10 IMF frequency domain signal in No.1 log selects by the combination of the time domain characteristics and component correlation coefficient

2#原木冲击回波信号在1 ms内没有较大的衰减(图11(a))，添加噪声幅度比例0.1、数量100，分解结果如图11(b)。与1#原木情况相似，IMF1～3为白噪声、IMF4为表面波、IMF10为残余分量、IMF7～9的相关系数较低，通过IMF5的频域特征计算缺陷位置为0.71 m，误差1.4%。

由表2可知，利用本文方法计算1#和2#两端、3#B端缺陷位置误差均在10%以内，可准确的判断出空洞位置。当空洞距离检测端0.3 m情况下，3#原木A端检测误差较大，其原因为空洞距离原木末端较近，冲击波在空洞和原木末端发生多次反射形成共振，反射波形成的共振频率大于冲击信号的频率造成，当铁球撞击原木端面时，撞击点振动形成入射波，当入射波的频率小于共振频率后，测试就会不准确[22]，但可从B端检测出空洞位置。并且从表1中可看出对缺陷敏感的IMF与原信号的互相关系数随着空洞的扩大也逐渐增加，这是由于空洞反射面扩大，反射的应力波在采集的信号中占有比例也随之增加。

(a) 2#白杨原木0.08 m×0.08 m尺寸空洞A端应力波时域信号

(b) 冲击回波信号分解结果图11 2#白杨原木0.08 m×0.08 m空洞尺寸下A端应力波信号Fig.11 The 0.08 m×0.08 m hole size log Stress wave signals of No.2 log model in end A

图12 2#原木基于时域特征和分量相关系数相结合选择出IMF频域信号Fig.12 IMF frequency domain signal in No.2 log selects by the combination of the time domain characteristics and component correlation coefficient

表1 各IMF与冲击回波信号相关系数Tab.1 Correlation coefficient between IMF and Impact-echo signal

表2 原木样本检测结果Tab.2 Log sample test results

4 结 论

本文提出了一种基于原木应力波信号分解的冲击回波法，用于原木径向空洞检测。通过有限元仿真模拟原木中应力波信号得出：设置合适的EEMD参数能够使得到的IMF与纵波信号接近，在有较强表面波和噪声的情况下需要添加的噪声数量较多、标准差比例较大。选择空洞位置和大小不同、长度不同的原木样本试验，验证了该方法的可行性。空洞距离检测端较近的情况下(0.3 m)准确率较低，需要进一步研究。对空洞敏感的IMF分量与冲击回波信号的相关系数和空洞大小的变化成正相关。

参 考 文 献

[1] ROSS R J, ZERBE J I, WANG X P, et al. Stress wave nondestructive evaluation of Douglas-fir peeler cores[J]. Forest Products Journal, 2005, 55(3): 38-43.

[2] WACKER J P, WANG X, ROSS R J, et al. Condition assessment of historic wood vessels[C]∥Proceedings of the 15th International Symposium on Nondestructive Testing of Wood. Duluth, MN, 2007.

[3] WANG, X. Effects of size and moisture on stress wave E-rating of structural lumber[C]∥Proceedings of the 10th World Conference on Timber Engineering. Miyazaki, Japan, 2008.

[4] 徐华东，王立海,游祥飞,等. 应力波和超声波在立木无缺陷断面的传播速度[J].林业科学, 2011, 47(4):129-134.

XU Huadong, WANG Lihai, YOU Xiangfei, et al. Propagation velocity of stress wave and ultrasonic wave transmitting on indefectible cross section of standing trees[J]. Scientia Silvae Sinicae, 2011, 47(4):129-134.

[5] 梁善庆，赵广杰,傅峰. 应力波断层成像诊断木材内部缺陷[J]. 木材工业, 2010, 24(5):11-13.

LIANG Shanqing, ZHAO Guangjie, FU Feng. Diagnosis of internal defects of wood with stress wave tomography[J]. China Wood Industry, 2010, 24(5):11-13.

[6] 王立海，徐华东,闫在兴,等. 传感器的数量与分布对应力波检测原木缺陷效果的影响[J]. 林业科学, 2008, 44(5):115-121.

WANG Lihai, XU Huadong, YAN Zaixing, et al. Effects of sensor quantity and planar distribution on testing results of log defects based on stress wave[J]. Scientia Silvae Sinicae, 2008, 44(5):115-121.

[7] 余观夏，张爱珍,史伯章,等. 用应力波频谱分析技术检测原木中的腐朽[J].东北林业大学学报,2007,35(10):20-25.

YU Guanxia, ZHANG Aizhen, SHI Bozhang, et al. Detection of timber decay by stress wave frequency spectrum[J]. Journal of Northeast Forestry University, 2007, 35(10):20-25.

[8] 王靖涛.桩基础设计与检测[M].武汉:华中科技大学出版社, 2005.

[9] 杨学春，王立海. 木材应力波无损检测研究[M].北京:科学出版社, 2011.

[10] WANG X，ROSS R J，BRASHAW B K，et al. Diameter effect on stress-wave evaluation of modulus of elasticity of small-diameter logs[J]. Wood and Fiber Science，2004，36(3)：368-377.

[11] WANG X. Acoustic measurements on trees and logs：a review and analysis[J]. Wood Science and Technology，2013，47(5)：965-975.

[12] SANSALONE M, STREETT W B. Impact echo: nondestructive evaluation of concrete and masonry[M]. Ithaca, NY: Bullbrier Press, 1997.

[13] ZHAO X, PATEL T H, ZUO M J. Multivariate EMD and full spectrum based condition monitoring for rotating machinery[J]. Mechanical Systems and Signal Processing, 2012, 27: 712-728.

[14] SHEN Chenxing, LI Zhixiong, QIN Li, et al. Recent progress on mechanical condition monitoring and fault diagnosis[J]. Procedia Engineering, 2011, 15: 142-146.

[15] ZHANG J, YAN R, GAO R X, et al. Performance enhancement of ensemble empirical mode decomposition[J]. Mech Syst Signal Process, 2010,24:2104-2123.

[16] WU Z, HUANG N E. Ensemble empirical mode decomposition: a noise assisted data analysis method[J]. Adv Adaptive Data Anal, 2008, 1(1):1-41.

[17] GUO W, TSE P W. Enhancing the ability of ensemble empirical mode decomposition in machine fault diagnosis[C]∥ Prognostics and Health Management Conference. IEEE, 2010.

[18] SCHEFFLER M，NIEMZ P，HARDTKE H J. Numerical simulation of sound propagation in wood in presence of defects[J]. Holz als Roh und Werkstoff,2002,60(5)：397-404.

[20] FENG Hailin, FANG Yiming, LI Jian, et al. Detecting Longitudinal of Internal Log Hole Using an Impact-Echo method[J]. Bio Resources,2015, 10(4): 7569-7579.

[21] GIBSON A, POPOVICS J S. Lamb wave basis for impact-echo method analysis[J]. Eng Mech, 2005,131(4):438-443.

[22] CARINO N J. The impact-echo method：an overview[C]∥Proceedings of the 2001 Structures Congress and Exposition. Washington DC：American Society Civil Engineers，2001.