CN105181805A - 一种基于时反music的多滤波超声成像方法 - Google Patents

Abstract

本发明公开了一种基于时反MUSIC的多滤波超声成像方法。该方法的具体步骤如下：(1)设置实验参数，将探头置于目标信号源的近场进行扫描，对扫描结果进行傅里叶分析得到在一定时间、一定频率下的响应矩阵。(2)由坐标变换对响应矩阵进行数据旋转得到包含噪声信号的两个副矩阵A₁和A₂。(3)根据多散射和单散射的特征，对矩阵A₁和A₂分别在矢量上进行映射，映射结果乘以矢量，分别得到近似单散射的信号矩阵和

Description

一种基于时反MUSIC的多滤波超声成像方法

技术领域

本发明涉及缺陷定位技术领域，特别是涉及一种粗晶材料的高分辨成像方法。

背景技术

粗晶粒材料在现代工业中应用越来越广泛，因此对其安全性检测成为无损检测领域的热点之一。但其因晶粒平均直径大，散射幅度曲线向低频方向偏移，与换能器频带发生较大交叠，散射现象明显，信噪比低，使得超声对缺陷的检出能力严重下降。传统超声相控阵对粗晶材料的缺陷定位手段大致可以总结为三类：1)采用低频窄带探头进行检测2)采用聚焦技术3)采用信号处理技术。其中低频窄带探头、聚焦，由于现有技术的限制，使得很难达到较为理想的状态，现多采用信号处理。目前在超声相控阵领域内对粗晶材料的缺陷定位比较流行的算法是成像算法，它是通过不断迭代，将超声波束的有效能量集中到缺陷所在的位置，从而获可以克服介质不均匀带来的声束畸变和失真。这种方法实现的基础是当信号进入直接中时，大部分只经过一次散射回到阵元，阵列响应矩阵的主对角线对目标之间有决定关系。但当介质进入多散射的环境，这种成像方法就失去了作用。为了克服这种局限性，设计了一种基于时反MUSIC算法克服介质多散射的成像方法，其基本思想是在多散射背景中分离目标的单散射回波信号，通过MUSIC算法成像，从而提高分辨率。

发明内容

本发明的目的在于针对现有技术的不足，提出了一种基于时反MUSIC算法克服介质多散射的成像方法。

本发明的目的是通过以下技术方案来实现的：一种基于时反MUSIC算法克服多散射的成像方法，包括以下步骤：

步骤1：在超声相控阵检测装置主机上安装探头和楔块，对检测装置进行校准，包括声速、楔块延时、灵敏度、编码器的校准，之后对超声相控阵实验试块进行扫描，具体为：

超声相控阵中的一个阵元i首先发射一束脉冲波，之后N个阵元接收并记录回波信号h_ij(t)。在完成第一次接收之后超声相控阵的其他阵元重复上述步骤，得到N²个脉冲响应组成的N×N的响应矩阵H(t)。经过傅里叶变换，时间信号被截断成连续的时间窗，表达式如下：

K_ij(T,t)＝h_ij(T-t)W_R(t)(1)

当t在区间[-Δt/2,Δt/2]时，W_R(t)＝1，否则W_R(t)＝0。在每个时间T下，k_ij的值构成了传输矩阵K。传输矩阵K经过快速傅里叶变换，得到在时间T和频率f下的响应矩阵K(T,f)。

步骤2：由坐标变换对响应矩阵K(T,f)进行数据旋转得到包含噪声信号的两个副矩阵A₁和A₂，具体如下：

对响应矩阵K(T,f)的坐标进行如下变化(x_i,x_j)→(y_u,y_v)

$y_u=\frac{x_i-x_j}{\sqrt{2}}$ 和 $y_v=\frac{x_i+x_j}{\sqrt{2}}---\left(2\right)$

变化之后得到两个副矩阵A₁和A₂，其中A₁的维数L₁是(2M-1)×(2M-1)，A₂的维数L₂是(2M-2)×(2M-2)，其中M＝(N+3)/4。A₁、A₂和响应矩阵K(T,f)的关系如式(3)和(4)：

a₁[u,v]＝k[u+v-1,v-u+2M-1](3)

a₂[u,v]＝k[u+v,v-u+2M-1](4)

步骤3：根据多散射和单散射的特征，对矩阵A₁和A₂分别在矢量上进行映射，映射结果乘以矢量，分别得到近似单散射的信号矩阵和具体为：

矩阵A₁和A₂分别在矢量S上投影，其中矢量S的坐标如式(5)所示：

$s_1=\exp \[jk\frac{y_u^2}{2R}\]{L_1}^{-1/2},s_2=\exp \[jk\frac{y_u^2}{2R}\]{L_2}^{-1/2}---\left(5\right)$

其中R＝cT/2，k＝2πf/c，c为超声波传播速度。投影结果P₁和P₂如式(6)所示：

$P_1=S_1^H{A}_1,P_2=S_2^H{A}_2---\left(6\right)$

对矩阵A₁和A₂进行滤波，滤波后的结果和可以用如下式(7)表示：

$A_1^F=S_1{P}_1=S_1{S}_1^H{A}_1,A_2^F=S_2{P}_2=S_2{S}_2^H{A}_2---\left(7\right)$

最终得到的滤波后的矩阵K^F中的元素如式(8)所示：

步骤4：对信号矩阵K^F进行特征值分解，得到特征值对应的噪声子空间，噪声子空间借助格林函数对目标进行成像，具体为：

对矩阵K^F进行奇异值分解，分解结果如式(9)所示：

K^F＝UΛU^H(9)

Λ是有N个特征值的对角矩阵，其中有M个有效特征值，其余有(N-M)个近似为0的特征值。U为特征值对应的特征向量，其中噪声子空间NS＝{u_M+1(ω),u_M+2(ω),…,u_N(ω)}。

最终的成像函数D(x_p,ω)为：

$D(x_p,\mathrm{\&omega;})={\left(\underset{i=M+1}{\overset{N}{\&Sigma;}}\right|<g(x_p,\mathrm{\&omega;}),u_i^{*}\left(\mathrm{\&omega;}\right)>\left|\right)}^{-1}---\left(10\right)$

其中g(x_p,ω)为格林函数矢量，表达式如下：

g(x_p,ω)＝[G(x_p,α₁,ω),G(x_p,α₂,ω),…,G(x_p,α_N,ω)]^T(11)

X_p表示搜索位置的矢量，a_i代表阵元位置，T表示转置。

本发明的有益效果是：本发明提出了一种包含在强散射介质中缺陷的成像方法，该方法可以克服多散射回波和声速畸变造成的缺陷定位不准确、分辨率不高等问题。并首次将多滤波和时反MUSIC算法结合，为粗晶材料缺陷的定位检测提供了一种新的思路。

附图说明

图1是本发明方法实现的结构框图。

具体实施方式

下面结合附图给出本发明的具体实施例，详细说明本发明的技术方案。

如图1所示，本发明一种基于时反MUSIC算法克服多散射的成像方法，包括以下步骤：

步骤1：在超声相控阵检测装置主机上安装探头和楔块，对检测装置进行校准，包括声速、楔块延时、灵敏度、编码器的校准，之后对超声相控阵实验试块进行扫描，具体为：

超声相控阵中的一个阵元i首先发射一束脉冲波，之后N个阵元接收并记录回波信号h_ij(t)。在完成第一次接收之后超声相控阵的其他阵元重复上述步骤，得到N²个脉冲响应组成的N×N的响应矩阵H(t)。经过傅里叶变换，时间信号被截断成连续的时间窗，表达式如下：

K_ij(T,t)＝h_ij(T-t)W_R(t)(1)

当t在区间[-Δt/2,Δt/2]时，W_R(t)＝1，否则W_R(t)＝0。在每个时间T下，k_ij的值构成了传输矩阵K。传输矩阵K经过快速傅里叶变换，得到在时间T和频率f下的响应矩阵K(T,f)。

步骤2：由坐标变换对响应矩阵K(T,f)进行数据旋转得到包含噪声信号的两个副矩阵A₁和A₂，具体如下：

对响应矩阵K(T,f)的坐标进行如下变化(x_i,x_j)→(y_u,y_v)

$y_u=\frac{x_i-x_j}{\sqrt{2}}$ 和 $y_v=\frac{x_i+x_j}{\sqrt{2}}---\left(2\right)$

变化之后得到两个副矩阵A₁和A₂，其中A₁的维数L₁是(2M-1)×(2M-1)，A₂的维数L₂是(2M-2)×(2M-2)，其中M＝(N+3)/4。A₁、A₂和响应矩阵K(T,f)的关系如式(3)和(4)：

a₁[u,v]＝k[u+v-1,v-u+2M-1](3)

a₂[u,v]＝k[u+v,v-u+2M-1](4)

步骤3：滤波的基础是单散射信号的响应矩阵的每一列对某个函数存在某种已知的关系，而都散射信号不具备这一特征，因此可以根据它们的差异性进行滤波。

根据多散射和单散射的特征，对矩阵A₁和A₂分别在矢量上进行映射，映射结果乘以矢量，分别得到近似单散射的信号矩阵和具体为：

矩阵A₁和A₂分别在矢量S上投影，其中矢量S的坐标如式(5)所示：

$s_1=\exp \&lsqb;jk\frac{y_u^2}{2R}\&rsqb;{L_1}^{-1/2},s_2=\exp \&lsqb;jk\frac{y_u^2}{2R}\&rsqb;{L_2}^{-1/2}---\left(5\right)$

其中R＝cT/2，k＝2πf/c，c为超声波传播速度。投影结果P₁和P₂如式(6)所示：

$P_1=S_1^H{A}_1,P_2=S_2^H{A}_2---\left(6\right)$

对矩阵A₁和A₂进行滤波，滤波后的结果和可以用如下式(7)表示：

$A_1^F=S_1{P}_1=S_1{S}_1^H{A}_1,A_2^F=S_2{P}_2=S_2{S}_2^H{A}_2---\left(7\right)$

最终得到的滤波后的矩阵K^F中的元素如式(8)所示：

步骤4：TR-MUSIC从字面意思上看是一种成像和时间反转技术结合的算法，在实现过程上，需要对时间反转算子进行奇异值分解，并借助于格林函数成像。具体过程如下：

对信号矩阵K^F进行特征值分解，得到特征值对应的噪声子空间，噪声子空间借助格林函数对目标进行成像，具体为：

对矩阵K^F进行奇异值分解，分解结果如式(9)所示：

K^F＝UΛU^H(9)

Λ是有N个特征值的对角矩阵，其中有M个有效特征值，其余有(N-M)个近似为0的特征值。U为特征值对应的特征向量，其中噪声子空间NS＝{u_M+1(ω),u_M+2(ω),…,u_N(ω)}。

最终的成像函数D(x_p,ω)为：

$D(x_p,\mathrm{\&omega;})={\left(\underset{i=M+1}{\overset{N}{\&Sigma;}}\right|<g(x_p,\mathrm{\&omega;}),u_i^{*}\left(\mathrm{\&omega;}\right)>\left|\right)}^{-1}---\left(10\right)$

其中g(x_p,ω)为格林函数矢量，表达式如下：

g(x_p,ω)＝[G(x_p,α₁,ω),G(x_p,α₂,ω),…,G(x_p,α_N,ω)]^T(11)

X_p表示搜索位置的矢量，a_i代表阵元位置，T表示转置。

Claims (1)

1.一种基于时反MUSIC算法克服多散射的成像方法，其特征在于，包括以下步骤：

步骤1：在超声相控阵检测装置主机上安装探头和楔块，对检测装置进行校准，包括声速、楔块延时、灵敏度、编码器的校准，之后对超声相控阵实验试块进行扫描，具体为：

超声相控阵中的一个阵元i首先发射一束脉冲波，之后N个阵元接收并记录回波信号h_ij(t)。在完成第一次接收之后超声相控阵的其他阵元重复上述步骤，得到N²个脉冲响应组成的N×N的响应矩阵H(t)。经过傅里叶变换，时间信号被截断成连续的时间窗，表达式如下：

K_ij(T,t)＝h_ij(T-t)W_R(t)(1)

当t在区间[-Δt/2,Δt/2]时，W_R(t)＝1，否则W_R(t)＝0。在每个时间T下，k_ij的值构成了传输矩阵K。传输矩阵K经过快速傅里叶变换，得到在时间T和频率f下的响应矩阵K(T,f)。

步骤2：由坐标变换对响应矩阵K(T,f)进行数据旋转得到包含噪声信号的两个副矩阵A₁和A₂，具体如下：

对响应矩阵K(T,f)的坐标进行如下变化(x_i,x_j)→(y_u,y_v)

$y_u=\frac{x_i-x_j}{\sqrt{2}}$ 和 $y_v=\frac{x_i+x_j}{\sqrt{2}}---\left(2\right)$

变化之后得到两个副矩阵A₁和A₂，其中A₁的维数L₁是(2M-1)×(2M-1)，A₂的维数L₂是(2M-2)×(2M-2)，其中M＝(N+3)/4。A₁、A₂和响应矩阵K(T,f)的关系如式(3)和(4)：

a₁[u,v]＝k[u+v-1,v-u+2M-1](3)

a₂[u,v]＝k[u+v,v-u+2M-1](4)

步骤3：根据多散射和单散射的特征，对矩阵A₁和A₂分别在矢量上进行映射，映射结果乘以矢量，分别得到近似单散射的信号矩阵和具体为：

矩阵A₁和A₂分别在矢量S上投影，其中矢量S的坐标如式(5)所示：

$s_1=\exp \&lsqb;jk\frac{y_u^2}{2R}\&rsqb;{L_1}^{-1/2},s_2=\exp \&lsqb;jk\frac{y_u^2}{2R}\&rsqb;{L_2}^{-1/2}---\left(5\right)$

其中R＝cT/2，k＝2πf/c，c为超声波传播速度。投影结果P₁和P₂如式(6)所示：

$P_1=S_1^H{A}_1,P_2=S_2^H{A}_2---\left(6\right)$

对矩阵A₁和A₂进行滤波，滤波后的结果和可以用如下式(7)表示：

$A_1^F=S_1{P}_1=S_1{S}_1^H{A}_1,A_2^F=S_2{P}_2=S_2{S}_2^H{A}_2---\left(7\right)$

最终得到的滤波后的矩阵K^F中的元素如式(8)所示：

步骤4：对信号矩阵K^F进行特征值分解，得到特征值对应的噪声子空间，噪声子空间借助格林函数对目标进行成像，具体为：

对矩阵K^F进行奇异值分解，分解结果如式(9)所示：

K^F＝UΛU^H(9)

Λ是有N个特征值的对角矩阵，其中有M个有效特征值，其余有(N-M)个近似为0的特征值。U为特征值对应的特征向量，其中噪声子空间NS＝{u_M+1(ω),u_M+2(ω),...,u_N(ω)}。

最终的成像函数D(x_p,ω)为：

$D(x_p,\mathrm{\&omega;})={\left(\underset{i=M+1}{\overset{N}{\&Sigma;}}\right|<g(x_p,\mathrm{\&omega;}),u_i^{*}\left(\mathrm{\&omega;}\right)>\left|\right)}^{-1}---\left(10\right)$

其中g(x_p,ω)为格林函数矢量，表达式如下：

g(x_p,ω)＝[G(x_p,α₁,ω),G(x_p,α₂,ω),...,G(x_p,α_N,ω)]^T(11)

X_p表示搜索位置的矢量，a_i代表阵元位置，T表示转置。