一种基于频率特征的木材应力波无损检测方法

摘要

本发明公开了一种基于频率特征的木材应力波无损检测方法，在原木周围均匀地安装压电式加速度传感器，用脉冲锤敲击其中编号为0的传感器，利用数据采集卡完成传感器输出信号的采集，保存；对采集到的信号做K点快速傅立叶变换，然后求出实测频率响应函数；根据实测频率响应函数及健康木的频率响应函数，构建观察矩阵；作基于二阶统计量的盲源分离，估计缺陷点与观测点之间的频率响应函数；最后以聚类结果为依据判别木材内部有无缺陷、缺陷数量以及缺陷大小。该根据木材频率响应函数的应力波无损检测方法，不受反射波、折射波信号的干扰，检测结果更加准确，能自动检测出木材内部是否存在缺陷、缺陷大小的信息，检测过程简便，实用性强。

说明书

技术领域

本发明属于木材质量检测技术领域，尤其涉及一种基于频率特征的木材应 力波无损检测方法。

背景技术

木材应力波无损检测技术的基本原理如图1所示，利用脉冲锤敲击木材， 使其内部产生应力波的传播，通过安装在原木周围的传感器，可采集到各观测 点上的应力波信号，对应力波信号进行时域、频域分析，确定木质材料的性质(如 弹性模量和缺陷等)。

基于应力波的木材无损检测技术具有实现简单、使用方便等优点，得到了 广泛应用，但通常用于反映木材内部缺陷检测的特征是应力波在木材中的传播 时间或传播速度，这种方法不能判断出木材中腐朽的位置，而且测量结果分散 性较大，检测结果很难达到高度准确性。

频率响应函数是系统输入信号和输出信号的傅立叶变换之比，它包含了木 材内部物理参数的所有信息。正是基于这一性质，中外众多研究者都提出了利 用木材频率响应特征的无损检测技术。随后出现了基于频域特征的应力波无损 检测技术，Bozhang等人研究了木材频率响应特性，健康木中应力波信号高频 成分衰减较小，可以应用这种方法有效识别木材中的早期腐朽，余观夏等人提 出应力波信号频谱的谱线面积与木材的成分有对应关系，可利用频谱面积大小 来表示对应木材成分的多少，并通过对已知内部腐朽情况的原木试样进行测试， 证实了该方法是可行的，而且具有测试结果非常直观等优点。

但是这些方法都没有考虑木材中应力波信号的反射、折射、透射等现象， 在实际检测过程中，传感器观测到的应力波信号中混杂了大量的反射波、折射 波等信号，给检测结果带来了误差。

如图2所示为一原木的横截面，A点为脉冲锤敲击点，B点为观测传感器 的安装位置，由于腐朽点F的存在，应力波除了沿直线AB向B点传播外(假 设这部分应力波信号为s₁(t))，B点还接收到反射信号s₂(t)。

根据频率响应函数的定义：

$H_{\mathrm{AB}}\left(\mathrm{j\omega }\right)=\frac{X_B\left(\mathrm{j\omega }\right)}{X_A\left(\mathrm{j\omega }\right)}=\frac{S_1\left(\mathrm{j\omega }\right)+S_2\left(\mathrm{j\omega }\right)}{X_A\left(\mathrm{j\omega }\right)}=\frac{S_1\left(\mathrm{j\omega }\right)}{X_A\left(\mathrm{j\omega }\right)}+\frac{S_2\left(\mathrm{j\omega }\right)}{X_A\left(\mathrm{j\omega }\right)}=H_1\left(\mathrm{j\omega }\right)+H_2\left(\mathrm{j\omega }\right)---\left(1\right)$

式中：X_A(jω)、X_B(jω)分别是A点、B点观测信号的傅立叶变换；

S₁(jω)、S₂(jω)分别是信号s₁(t)、s₂(t)的傅立叶变换。

显然，真正反映AB间木材物理特性的应该是H₁(jω)，若简单地以观测信号 傅立叶变换之比H_AB(jω)表征AB间木材的物理特性，必然会导致较大误差。

发明内容

本发明提供了一种基于频率特征的木材应力波无损检测方法，旨在解决现 有技术提供的利用木材频率响应特征的无损检测方法，都没有考虑木材中应力 波信号的反射、折射、透射现象，传感器观测到的应力波信号中混杂了大量的 反射波、折射波信号，检测结果的误差较大的问题。

本发明的目的在于提供一种基于频率特征的木材应力波无损检测方法，该 方法包括以下步骤：

步骤一，在原木周围均匀地安装I个通过信号线缆连接至数据采集卡的压 电式加速度传感器，用脉冲锤敲击其中编号为0的传感器，然后通过数据采集 卡完成压电式加速度传感器输出信号的采集，保存为x₀(n)，x₁(n)，...，x_i(n)...， x_I-1(n)；

步骤二，对采集到的信号做K点快速傅立叶变换，得到X₀(k)，X₁(k)，...， X_i(k)，...，X_I-1(k)，k＝0，1，ΛK-1，然后求出第0个传感器至其他各观测传感器 之间的实际频率响应函数：

$H_i\left(k\right)=\frac{X_i\left(k\right)}{X_0\left(k\right)},i=1,\mathrm{\Lambda I}-1;---\left(2\right)$

步骤三，定义相同条件下无缺陷木材中第0个传感器至其他各观测传感器 之间的频率响应函数为H_0，i(k)，然后将H_i(k)减去H_0，i(k)，得到：

Y_i(k)＝H_i(k)-H_0，i(k)，i＝1，ΛI-1；             (3)

步骤四：将Y＝[Y′₁(k)，Y′₂(k)，Λ，Y′_I-1(k)]′作为观察矩阵，对其作基于二阶统计量 的盲源分离，估计出混合矩阵W＝[G′_0，1(k)，G′_0，2(k)，ΛG′_0，J(k)]′和源向量G_j，i(k)， i＝1，ΛI-1，j＝1，ΛJ，其中G_0，j(k)为第0个传感器至木材内部第j个缺陷点之间 的频率响应函数，G_j，i(k)为第j个缺陷点至第i个观测点之间的频率响应函数， 详细步骤为：

(1)估计Y的相关矩阵

${\hat{R}}_Y\left(0\right)=\frac{1}{K}\underset{k=0}{\overset{K-1}{\Sigma }}Y\left(k\right)Y^{\prime }\left(k\right)---\left(4\right)$

(2)对作特征值分解(EVD)

${\hat{R}}_Y\left(0\right)=U_Y{\Sigma }_Y{V}_Y^{\prime }=V_Y{\Lambda }_Y{V}_Y^{\prime }=V_S{\Lambda }_S{V}_S^{\prime }+V_N{\Lambda }_N{V}_N^{\prime }---\left(5\right)$

其中(I-1)×J维矩阵V_M＝[v₁，v₂，Λv_I-1]是与J个按下降顺序排列的主要特征值 Λ_S＝diag{λ₁≥λ₂Λ≥λ_J}相对应的特征矢量；(I-1)×(I-1-J)维矩阵V_N包含(I-1-J) 个噪声特征Λ_N＝{λ_J+1≥Λ≥λ_I-1}对应的噪声特征矢量，并且λ_J＞λ_J+1；

(3)白噪声方差估计成(I-1-J)个不重要特征值的均值；

(4)进行稳健的预白化变换：

$\bar{Y}\left(k\right)={\hat{\Lambda }}_S^{-1/2}{V}_S^{\prime }Y\left(k\right)=\mathrm{QY}\left(k\right)---\left(6\right)$

其中${\hat{\Lambda }}_S=\mathrm{diag}\{({\lambda }_1-{\hat{\sigma }}_K^2),({\lambda }_2-{\hat{\sigma }}_K^2),\Lambda ({\lambda }_1-{\hat{\sigma }}_K^2)\};$

(5)对于给定的p≠0，估计矢量的协方差阵，并进行协方差阵的奇 异值分解：

${\hat{R}}_{\bar{Y}}\left(p\right)=\frac{1}{K}\underset{k=0}{\overset{K-1}{\Sigma }}\bar{Y}\left(k\right){\bar{Y}}^{\prime }(k-p)=U_{\bar{Y}}{\Sigma }_{\bar{Y}}{V}_{\bar{Y}}^{\prime }---\left(7\right)$

(6)对于给定的p，检查对角阵所有奇异值是否不同，如果相同，对 于不同的时滞p重复步骤(5)，如果奇异值不同的，而且彼此远离，则估计出 源向量：

${\hat{G}}_j\left(k\right)=U_{\bar{Y}}^{\prime }{\hat{\Lambda }}_S^{-1/2}{V}_S^{\prime }Y\left(k\right),j=\mathrm{1,2},\Lambda ,J;---\left(8\right)$

步骤五，对估计出的J个进行k均值聚类分析，具体做法为：

(1)对于给定大小为J的数据集，令L＝1，选取3个初始聚类中心Z_i(L)， i＝1，2，3；

(2)计算每个数据对象与k个聚类中心的距离j＝1，2，Λ，J，i＝1，2，3，如果满足$D({\hat{G}}_j\left(k\right),Z_i\left(L\right))=\min \{D({\hat{G}}_j\left(k\right),Z_i\left(L\right)),j=\mathrm{1,2},\Lambda ,J\},$则 ${\hat{G}}_j\left(k\right)\in C_k;$

(3)重新计算k个新的聚类中心i＝1，2，3，以及平方误差 准则函数的值：$E(L+1)=\underset{m=1}{\overset{k}{\Sigma }}\underset{p\in C_m}{\Sigma }{\left|\right|p-Z_i\left(m\right)\left|\right|}^2;$

(4)判断：若||E(L+1)-E(L)||＜ε，则算法结束，否则L＝L+1，返回继续执 行(2)。

步骤六，根据聚类的结果判别木材内部有无缺陷，类别数量决定缺陷区域 的多少，每类中数据对象的数量表示该缺陷区域的大小。

进一步，在步骤三中，H_0，i(k)反映的是应力波信号沿球面方向的传播特征， 与缺陷信息无关，可由健康木材内的频率响应函数得出。

进一步，在步骤五中，G_j，i(k)为第j个缺陷点至第i个观测点之间的频率响 应函数，包含了大量的关于缺陷的信息。

进一步，对于健康木材，可视为匀质木材，各G_j，i(k)的数值大小相等；对于 有缺陷的木材，缺陷点上的G_j，i(k)数值与其他区域不同，相近的缺陷点数值大致 相同。

本发明提供的基于频率特征的木材应力波无损检测方法，通过基于二阶统 计量的盲源分离方法，获取缺陷点与观测点之间的频率响应函数，并以缺陷点 与观测点之间的频率响应函数为依据，对木材缺陷进行检测；该根据木材频率 响应函数的应力波无损检测的方法，不受反射波、折射波信号的干扰，检测结 果更加准确，并能自动检测出木材是否存在缺陷、缺陷大小的信息，检测过程 简便，实用性强，具有较强的推广与应用价值。

附图说明

图1是现有技术提供的利用应力波法确定木质材料性质的基本原理示意 图；

图2是现有技术提供的利用木材频率响应特征的无损检测方法的基本原理 示意图；

图3是本发明实施例提供的基于频率特征的木材应力波无损检测方法的原 理示意图；

图4是本发明实施例的接线图。

具体实施方式

本发明的主要思想如下：

参见图3，为一有缺陷的原木截面图，图中阴影部分表示木材腐朽空洞。 当脉冲锤敲击编号为0的传感器时，应力波会在木材内部传播。各观测传感器 接收的应力波信号x_i(n)中，除了从敲击点直接传播过来的应力波信号外，还包 含有从腐朽空洞反射过来的应力波信号。把大的腐朽空洞看作J个小缺陷点 f₁，f₂...，f_J的集合，根据傅立叶变换具有线性性的特点，各传感器接收到的信 号频谱为：

$X_i\left(k\right)=X_0\left(k\right)H_{0,i}\left(k\right)+X_0\left(k\right)\underset{j=1}{\overset{J}{\Sigma }}{G}_{0,j}\left(k\right)G_{j,i}\left(k\right)---\left(9\right)$

式中，X₀(k)为第0个传感器观测到的应力波信号频谱，H_0，i(k)为假定第0 个传感器与第i个传感器之间为健康木材时的频率响应函数，G_0，j(k)为第0个传 感器至第j个缺陷单元之间的频率响应函数，G_j，i(k)为第j个缺陷单元至第i个传 感器之间的频率响应函数。

式(9)中右边第一项X₀(k)H_0，i(k)为从健康木材中沿球面传播至第i个传感 器的应力波信号频谱，等式右边第二项为各缺陷点反射至第i 个传感器的应力波信号频谱。

将式(9)左右两边同时除以X₀(k)，再将H_0，i(k)移至等式另一边，得到：

$Y_i\left(k\right)=\underset{j=1}{\overset{J}{\Sigma }}{G}_{0,j}\left(k\right)G_{j,i}\left(k\right)=\frac{X_i\left(k\right)}{X_0\left(k\right)}-H_{0,i}\left(k\right)=H_i\left(k\right)-H_{0,i}\left(k\right)---\left(10\right)$

如前所述，G_j，i(k)为第j个缺陷单元至第第i个传感器之间的频率响应函 数，包含了大量的关于缺陷的信息。对于健康木材，可以视为匀质木材，各 G_j，i(k)的数值大小相等；对于有缺陷的木材，缺陷点上的G_j，i(k)数值与其他区 域不同，而且相近的缺陷点数值大致相同。因此可以通过G_j，i(k)识别木材中有 无缺陷，是单个区域出现缺陷，还是多个区域出现缺陷，缺陷区域的大小。

因此木材缺陷检测的问题转换为求取G_j，i(k)，在式(10)中，G_0，j(k)和G_j，i(k) 均为未知信息，唯一可以利用的信息就是Y_i(k)。本发明对Y_i(k)作基于二阶统 计量的盲源分离方法，进而估计出G_j，i(k)，获得木材内部缺陷信息。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实 施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施例 仅仅用以解释本发明，并不用于限定发明。

本发明的一种具体实施方式为：一种基于频率特征的木材应力波无损检测 方法，其具体做法是：

步骤一，在原木周围均匀地安装6个压电式加速度传感器，加速度传感器 输出信号经电荷放大器接入数据采集卡，接线图如图4所示。参见图3，用脉 冲锤敲击第0个传感器，6个传感器将接收到的应力波信号转换成电信号输入 数据采集卡。在上位机软件的控制下完成数据的转换等工作，并记录信号x₀(n)， x₁(n)，x₂(n)，x₃(n)，x₄(n)，x₅(n)。选择抽样周期为1μs，记录信号长度为10ms。

步骤二，对x₀(n)，x₁(n)，x₂(n)，x₃(n)，x₄(n)，x₅(n)作快速傅立叶变换，变 换点数选择为1024，得到频谱X₀(k)，X₁(k)，X₂(k)，X₃(k)，X₄(k)，X₅(k)；然 后根据式(2)求第0个传感器至其他5个传感器之间的频率响应函数：

$H_1\left(k\right)=\frac{X_1\left(k\right)}{X_0\left(k\right)},$$H_2\left(k\right)=\frac{X_2\left(k\right)}{X_0\left(k\right)},$$H_3\left(k\right)=\frac{X_3\left(k\right)}{X_0\left(k\right)},$$H_4\left(k\right)=\frac{X_4\left(k\right)}{X_0\left(k\right)},$$H_5\left(k\right)=\frac{X_5\left(k\right)}{X_0\left(k\right)}.$

步骤三，根据同树种健康木材中应力波信号沿球面方向传播特征，得到第 0个传感器至其他5个传感器在无缺陷情况下的传递函数，H_0，1(k)，H_0，2(k)， H_0，3(k)，H_0，4(k)，H_0，5(k)；再将步骤二中的实际传递函数与之相减，得到：

Y₁(k)＝H₁(k)-H_0，1(k)

Y₂(k)＝H₂(k)-H_0，2(k)

Y₃(k)＝H₃(k)-H_0，3(k)

Y₄(k)＝H₄(k)-H_0，4(k)

Y₅(k)＝H₅(k)-H_0，5(k)

步骤四，步骤三中的Y_i(k)也是第0个传感器至缺陷点间频率响应函数G_0,j(k) 与缺陷点至各观测点之间的频率响应函数G_j，i(k)的瞬时混合，即： 在G_0，j(k)与G_j，i(k)都未知的情况下，可以采用基于二阶统计 量的盲源分离方法进行盲源估计。将Y₁(k)，Y₂(k)，Y₃(k)，Y₄(k)，Y₅(k)当作观测 向量，采用基于二阶统计量的盲源分离方法，估计出G_0，j(k)与G_j，i(k)。具体做法 为：

(1)估计矩阵Y＝[Y′₁(k)，Y′₂(k)，Y′₃(k)，Y′₄(k)，Y′₅(k)]′的相关矩阵为：

${\hat{R}}_Y\left(0\right)=\frac{1}{1024}\underset{k=0}{\overset{1023}{\Sigma }}Y\left(k\right)Y^{\prime }\left(k\right)$

(2)对作特征值分解(EVD)：

${\hat{R}}_Y\left(0\right)=U_Y{\Sigma }_Y{V}_Y^{\prime }=V_Y{\Lambda }_Y{V}_Y^{\prime }=V_S{\Lambda }_S{V}_S^{\prime }+V_N{\Lambda }_N{V}_N^{\prime }$

其中5×J维矩阵V_M＝[v₁，v₂，L v₅]是与J个按下降顺序排列的的主要特征值 Λ_S＝diag{λ₁≥λ₂L≥λ_J}相对应的特征矢量；5×(5-J)维矩阵V_N包含(5-J)个噪声 特征Λ_N＝{λ_J+1≥L≥λ₅}对应的噪声特征矢量，并且λ_J＞λ_J+1；

(3)白噪声方差估计成(5-J)个不重要特征值的均值；

(4)进行稳健的预白化变换：

$\bar{Y}\left(k\right)={\hat{\Lambda }}_S^{-1/2}{V}_S^{\prime }Y\left(k\right)=\mathrm{QY}\left(k\right)$

其中${\hat{\Lambda }}_S=\mathrm{diag}\{({\lambda }_1-{\hat{\sigma }}_{1024}^2),({\lambda }_2-{\hat{\sigma }}_{1024}^2),L({\lambda }_5-{\hat{\sigma }}_{1024}^2)\};$

(5)对于给定的p≠0，估计矢量的协方差阵，并进行协方差阵的奇异 值分解，${\hat{R}}_{\bar{Y}}\left(p\right)=\frac{1}{1024}\underset{k=0}{\overset{1023}{\Sigma }}\bar{Y}\left(k\right){\bar{Y}}^{\prime }(k-p)=U_{\bar{Y}}{\Sigma }_{\bar{Y}}{V}_{\bar{Y}}^{\prime };$

(6)对于给定的p，检查对角阵所有奇异值是否不同，如果相同，对于 不同的时滞p重复步骤(5)，如果奇异值不同的，而且彼此远离，则估计出源 向量${\hat{G}}_j\left(k\right)=U_{\bar{Y}}^{\prime }{\hat{\Lambda }}_S^{-1/2}{V}_S^{\prime }Y\left(k\right),$j＝1，ΛJ。

步骤五，对估计出的J个进行k均值聚类分析，具体做法为：

(1)对于给定大小为J的数据集，令L＝1，选取3个初始聚类中心Z_i(L)， i＝1，2，3；

(2)计算每个数据对象与k个聚类中心的距离j＝1，2，Λ，J，i＝1，2，3，如果满足$D({\hat{G}}_j\left(k\right),Z_i\left(L\right))=\min \{D({\hat{G}}_j\left(k\right),Z_i\left(L\right)),j=\mathrm{1,2},\Lambda ,J\},$则 ${\hat{G}}_j\left(k\right)\in C_k;$

(3)重新计算k个新的聚类中心i＝1，2，3，以及平方误差 准则函数的值：$E(L+1)=\underset{m=1}{\overset{k}{\Sigma }}\underset{p\in C_m}{\Sigma }{\left|\right|p-Z_i\left(m\right)\left|\right|}^2;$

(4)判断：若||E(L+1)-E(L)||＜ε，则算法结束，否则L＝L+1，返回继续执 行(2)。

步骤六，根据聚类的结果判别木材内部有无缺陷，类别数量决定缺陷区域 的多少，每类中数据对象的数量表示该缺陷区域的大小。

本发明实施例提供的基于频率特征的木材应力波无损检测方法，通过基于 二阶统计量的盲源分离方法，获取缺陷点与观测点之间的频率响应函数，并以 缺陷点与观测点之间的频率响应函数为依据，对木材缺陷进行检测；该根据木 材频率响应函数的应力波无损检测的方法，不受消除反射波、折射波信号的干 扰，检测结果更加准确，并能自动检测出木材是否存在缺陷、缺陷大小的信息， 检测过程简便，实用性强，具有较强的推广与应用价值。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发 明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明 的保护范围之内。