CN102279394A - 低频窄带高分辨超声探测成像方法 - Google Patents

Abstract

本发明公开了一种低频窄带超声探测成像方法，主要解决低频窄带超声探测纵向分辨率低的问题，其实现步骤是：首先利用超声探测器发射低频窄带信号对远场目标进行探测，将远场目标的回波信号进行采样离散化；然后由离散化的回波信号构造出稀疏基矩阵，并将回波信号写成矩阵形式；最后构造基于CS的高分辨超声探测成像模型，并采用迭代加权最小二乘法求解模型重构出场景目标散射系数，将系数取模值即得到纵向高分辨的低频窄带超声探测成像。本发明相对传统的匹配滤波方法，成像的分辨率显著提高；相对现有的高分辨成像方法，克服了它们局限于高频宽带超声探测成像的不足，可用于人体深部组织、大型工业材料内部伤裂和杂质等探测。

Description

低频窄带高分辨超声探测成像方法

技术领域

本发明属于超声探测领域，特别是一种低频窄带高分辨超声探测成像方法，适合于获取常规超声无法探测到的远场信息，得到远场目标的高分辨率成像，从而为决策者提供更有效的信息。

技术背景

众所周知，超声成像因无创伤，经济且安全已被广泛应用于医疗诊断、工业探伤、海洋水下探测等。为了获得更多的观测对象信息，高清晰度/高分辨率超声影像一直是人们追求的目标。在工业探伤中，高分辨率甚至超高分辨率的超声影像是发现微小有害伤损和杂质的重要保证；在海防目标探测中，高分辨率超声影像是正确识别敌方目标的前提。

为了提高超声探测成像的分辨率，人们提出的一些高分辨超声成像方法，主要包括以下两类：

1、频域法；该方法采用频域外推法实现高分辨。1964年，Harr的文章“Diffractionand Resolving Power”奠定了高分辨成像的基础，他首次采用信号外推的方法成功提高了成像的分辨率。

2、空域法；该方法多采用自适应波束形成方法，如F.Viola等在文献“Time-DomainOptimized Near-Field Estimator for Ultrasound Imaging：Initial Development andresult”提出了一种新的自适应波束形成的方法应用于医学超声，获得了超声成像分辨率的提高。

以上现有的高分辨超声探测成像方法受限于成像系统的带宽和频率，只能用于宽带高频探测机制下的超声高分辨成像。发射宽带高频探测波进行探测的好处是可以获得高的纵向分辨率，但是缺点是宽带高频探测波传播距离有限，其能量随传播距离迅速呈指数衰减很难到达中远场，所以中远场目标不易探测。而在很多应用中，远场高分辨率的需求尤为重要，如医生期望看清人体深部组织信息；工程师需要探测大型工业材料内部伤裂和杂质等。

为了满足远场探测的需求，就必须采用低频窄带探测进行超声成像，但面临的重大科学难题是：基于传统匹配滤波机制下的低频窄带探测超声成像其纵向分辨率非常低，而上述现有的高分辨成像方法又局限于宽带高频探测超声成像，在低频窄带探测超声成像中无法发挥作用。其根本原因是发射低频窄带探测波进行探测的回波数据在频域发生了混叠，所以现有的高分辨成像方法无法从混叠的回波数据中重构出高分辨的影像。

发明内容

本发明的目的在于针对上述已有技术的问题，提出一种低频窄带高分辨超声探测成像方法，以从目标的混叠信息中重构出高分辨率的影像。

实现本发明的技术方案是：将近来兴起的压缩感知理论应用于超声探测成像中，该压缩感知CS理论，是信号处理领域的重大成果，其核心思想是可以从远低于奈奎斯特采样率的有限采样数据中重构信号，而低于奈奎斯特采样率的信号其频域会发生混叠，也就是说CS可以从混叠的采样数据中恢复原始信号，这就为攻克低频窄带探测超声成像纵向分辨率低的科学难题带来了重大的契机。根据这一理论，本发明首先发射低频窄带探测波获取离散的回波数据；然后构建超声回波的稀疏基，并构建基于压缩感知的低频窄带高分辨超声成像模型；最后求解模型得到高分辨率的成像，其具体实现步骤包括如下：

(1)构造离散回波信号模型：

(1a)根据发射的低频窄带波信号为s(t)，将回波信号r(t)表示为发射的低频窄带波信号在不同延时的叠加，即：

$r\left(t\right)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}s(t-{\mathrm{\τ}}_i)+n\left(t\right)$

其中t表示连续时间，τ_i表示第i个目标的延时量，N表示目标散射点的个数，s(t-τ_i)表示第i个目标的回波信号，α_i为表示i个目标的散射系数，n(t)为加性噪声；

(1b)对回波信号r(t)进行采样，令采样间隔为T_s＝1/f_s，其中f_s为采样频率，得到离散化后的回波信号r_k，即：

$r_k=r\left({\mathrm{kT}}_s\right)=\underset{i}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}s({\mathrm{kT}}_s-{\mathrm{\τ}}_i)+n_k$

其中n_k＝n(kT_s)表示离散后的噪声信号，s(kT_s-τ_i)表示第i个目标的离散回波信号，k表示第k个采样点，k＝1，2，L M，M表示采样个数；

(2)将离散化后的回波信号r_k写成矩阵形式，得到回波信号向量r：

r＝[r₁，r₂，...r_k，....r_M]^T＝Aα+n

其中，T表示向量的转置，A是根据不同目标的离散回波信号s(kT_s-τ_i)构建的稀疏基矩阵，α＝[α₁，α₂，...α_k，....]^T为超声探测成像系数向量，n＝[n₁，n₂，...n_k，....n_M]^T为加性噪声向量；

(3)由以上回波信号向量r，构建基于压缩感知的高分辨率超声探测成像模型，即：

$\underset{\mathrm{\α}\∈R^N}{\min }{E}^{\left(p\right)}\left(\mathrm{\α}\right),s.t{\left|\right|r-\mathrm{A\α}\left|\right|}_2^2\≤\mathrm{\ζ}.$

其中min表示最小化，R^N表示N维空间，

ζ为噪声方差；

(4)按如下迭代加权最小二乘法对上述基于压缩感知的高分辨率超声探测成像模型进行求解，得到一维高分辨率的超声探测成像系数向量α：

(4a)初始化：

k＝0，e_k，e_min＝10^-8，

其中k表示第k次迭代，表示A的伪逆，α_k表示当前迭代得到的超声成像系数向量，e_k是一个可调参数，初始设置在0～1之间，e_min表示e_k在迭代中的最小值；

(4b)由α_k构造系数的权值矩阵

p是一个0到1之间的自然数，并更新稀疏基矩阵A_k，得到A_k+1：

A_k+1＝AW_k+1

(4c)由A_k+1更新α_k，得到α_k+1：

${\mathrm{\α}}_{k+1}=W_{k+1}{A}_{k+1}^{+}r$

(4d)判断如下条件是否成立：

${\left|\right|{\mathrm{\&alpha;}}_{k+1}-{\mathrm{\&alpha;}}_k\left|\right|}_2<\sqrt{{\mathrm{\&epsiv;}}_k}/100$

其中

i表示系数向量α_k+1和α_k中第i个元素，α_k+1[i]表示系数向量α_k+1中第i个元素，α_k[i]表示系数向量α_k中第i个元素，

如果不成立，返回步骤(4b)继续迭代；

如果成立，再判断是否满足条件ε_k＜ε_min，如果满足，则停止迭代，得到一维高分辨率的超声探测成像系数向量α＝α_k+1；否则更新ε_k＝ε_k/10，返回步骤(4b)继续迭代，直到满足条件ε_k＜ε_min，得到一维高分辨率的超声探测成像系数向量α＝α_k+1；

(5)将得到的超声探测成像系数α取模值，即得到低频窄带高分辨超声探测成像。

本发明具有如下优点：

A.相比传统的基于匹配滤波的超声成像，由于本发明采用构建基于压缩感知的非凸函数最小化模型和迭代加权最小二乘的工程优化方法，超声成像分辨率得到了质的提高；

B.相比现有的高分辨率成像方法如：频域法、空域法，由于本发明采用压缩感知理论，成功解决了低频窄带探测成像纵向分辨率低的问题。

附图说明

图1为本发明的总流程图；

图2为本发明中迭代加权最小二乘法子流程图；

图3为本发明仿真使用的10个目标点的真实幅值和位置信息图；

图4为不含噪声情况下用现有匹配滤波方法和本发明方法进行探测的结果对比图；

图5为含20dB噪声时用现有匹配滤波方法和本发明方法探测结果对比图。

具体实施方式

参见图1，本发明具体实现步骤如下：

步骤一，构造离散回波数据模型。

(1.1)用超声探测器发射的低频窄带波信号s(t)，对目标场景进行探测，假设目标场景有N个散射点；

(1.2)用超声探测器接收目标的回波信号，根据发射的低频窄带波信号s(t)，将回波信号r(t)表示为发射的低频窄带波信号在不同延时的叠加，即：

$r\left(t\right)=\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{i}s(t-{\mathrm{\&tau;}}_i)+n\left(t\right)$

其中t表示连续时间，τ_i表示第i个目标的延时量，s(t-τ_i)表示第i个目标的回波信号，α_i表示i个目标的散射系数，n(t)为加性噪声；

(1.3)对回波信号r(t)进行离散化采样，令采样间隔为T_s＝1/f_s，其中f_s为采样频率，得到离散化后的回波信号r_k，即：

$r_k=r\left({\mathrm{kT}}_s\right)=\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{i}s({\mathrm{kT}}_s-{\mathrm{\&tau;}}_i)+n_k$

其中n_k＝n(kT_s)表示离散后的噪声信号，s(kT_s-τ_i)表示第i个目标的离散回波信号，k表示第k个采样点，k＝1，2，L M，M表示采样个数。

步骤二，将离散化后的回波信号r_k写成矩阵形式，得到回波信号向量r。

(2.1)由不同延时的离散回波信号s(kT_s-τ_i)构建稀疏基A：

A＝[s(kT_s-Δτ)，s(kT_s-2Δτ)K，s(kT_s-nΔτ)，K，s(kT_s-LΔτ)]

其中Δτ＝T_s，n＝1，2，LL，L＝τ_max/Δτ，这里τ_max＝max(τ_i)，即最大延时量；

(2.2)由稀疏基A重写回波信号模型：

r＝Aα+n

其中r＝[r₁，r₂，...r_k，....r_M]^T表示回波信号向量，T表示向量的转置，α＝[α₁，α₂，...α_k，....α_L]^T为超声探测成像系数向量，n＝[n₁，n₂，...n_k，....n_M]^T为加性噪声向量。

步骤三，由以上回波信号向量r，构建基于压缩感知的高分辨率超声探测成像模型，即：

$\underset{\mathrm{\&alpha;}\&Element;R^N}{\min }{E}^{\left(p\right)}\left(\mathrm{\&alpha;}\right),s.t{\left|\right|r-\mathrm{A\&alpha;}\left|\right|}_2^2\&le;\mathrm{\&zeta;}.$

其中min表示最小化，R^N表示N维空间，ζ为噪声方差；

步骤四，参加图(2)按迭代加权最小二乘法对上述基于压缩感知的高分辨率超声探测成像模型进行求解，得到一维高分辨率的超声探测成像系数向量α。

(4.1)初始化：

k＝0，e_k，e_min＝10^-8，

其中k表示第k次迭代，

表示A的伪逆，α_k表示当前迭代得到的超声成像系数向量，e_k是一个可调参数，初始设置在0～1之间，e_min表示e_k在迭代中的最小值；

(4.2)由α_k构造系数的权值矩阵

p是一个0到1之间的自然数，并更新稀疏基矩阵A_k，得到A_k+1：

A_k+1＝AW_k+1

(4.3)由A_k+1更新α_k，得到α_k+1：

${\mathrm{\&alpha;}}_{k+1}=W_{k+1}{A}_{k+1}^{+}r$

(4.4)判断如下条件是否成立：

${\left|\right|{\mathrm{\&alpha;}}_{k+1}-{\mathrm{\&alpha;}}_k\left|\right|}_2<\sqrt{{\mathrm{\&epsiv;}}_k}/100$

其中

i表示系数向量α_k+1和α_k中第i个元素，α_k+1[i]表示系数向量α_k+1中第i个元素，α_k[i]表示系数向量α_k中第i个元素；

如果条件不成立，返回步骤(4b)继续迭代；

如果条件成立，再判断是否满足条件ε_k＜ε_min，如果满足，则停止迭代，得到一维高分辨率的超声探测成像系数向量α＝α_k+1；否则更新ε_k＝ε_k/10，返回步骤(4b)继续迭代，直到满足条件ε_k＜ε_min，得到一维高分辨率的超声探测成像系数向量α＝α_k+1。

步骤五，将得到的超声探测成像系数α取模值，即得到低频窄带高分辨超声探测成像。

本发明的效果通过以下仿真进一步说明：

1.仿真条件：

运行系统为Intel(R)Core(TM)2 Duo CPU E4500 2.20GHz 2.20GHz的32位Windows操作系统，仿真程序采用MATLAB程序设计语言实现。仿真单频包络的发射信号，信号的中心频率为200KHz，包络形式为高斯窗函数，发射信号持续时间为43us，采样频率为560KHz，p取1，超声波在介质中传播速度设定为1500m/s；

2.仿真内容：

参照图3构造有10个目标的探测场景，其位置分别位于[5 6 7 9 12 17 27 29 3031]*1.25us；其幅值分别为：[0.7 0.7 0.4 0.6 0.7 0.3 0.7 0.5 0.75 1.0]；构造的回波信号带宽为58.24KHz；针对构造的场景分别用传统的匹配滤波的方法和本发明的方法进行目标探测；

(1)在不含噪声的情况下，采用现有的匹配滤波方法和本发明的方法进行探测，探测结果如图4。

(2)在含20dB噪声的情况下，采用现有的匹配滤波的方法和本发明的方法进行探测，探测结果如图5。

3.仿真结果：

从图4可以看出，对于分隔较远的目标，利用匹配滤波的方法可以探测，然而对于距离较近的目标，传统的匹配滤波就不能有效的探测出来；而本发明分方法将所有的目标都精确的探测了出来；本发明方法的分辨率比传统的方法有了明显的提升；

从图5可以看出，在含有20dB噪声的情况下，本发明的方法和传统的方法的探测能力都会受到影响，但是传统的方法对于分隔距离较近的目标仍是不能有效的探测，而本发明的高分辨的方法还是能够保持准确的探测目标。

Claims (2)

1.一种低频窄带高分辨超声探测成像方法，包括如下步骤：

(1)构造离散回波信号模型：

(1a)根据发射的低频窄带波信号s(t)，将回波信号r(t)表示为发射的低频窄带波信号在不同延时的叠加，即：

$r\left(t\right)=\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{i}s(t-{\mathrm{\&tau;}}_i)+n\left(t\right)$

其中t表示连续时间，τ_i表示第i个目标的延时量，N表示目标散射点的个数，s(t-τ_i)表示第i个目标的回波信号，α_i表示i个目标的散射系数，n(t)为加性噪声；

(1b)对回波信号r(t)进行离散化采样，令采样间隔为T_s＝1/f_s，其中f_s为采样频率，得到离散化后的回波信号r_k，即：

$r_k=r\left({\mathrm{kT}}_s\right)=\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_{i}s({\mathrm{kT}}_s-{\mathrm{\&tau;}}_i)+n_k$

其中n_k＝n(kT_s)表示离散后的噪声信号，s(kT_s-τ_i)表示第i个目标的离散回波信号，k表示第k个采样点，k＝1，2，L M，M表示采样个数；

(2)将离散化后的回波信号r_k写成矩阵形式，得到回波信号向量r：

r＝[r₁，r₂，...r_k，....r_M]^T＝Aα+n

其中，T表示向量的转置，A是根据不同目标的离散回波信号s(kT_s-τ_i)构建的稀疏基矩阵，α＝[α₁，α₂，...α_k，....]^T为超声探测成像系数向量，n＝[n₁，n₂，...n_k，....n_M]^T为加性噪声向量；

(3)由以上回波信号向量r，构建基于压缩感知的高分辨率超声探测成像模型，即：

$\underset{\mathrm{\&alpha;}\&Element;R^N}{\min }{E}^{\left(p\right)}\left(\mathrm{\&alpha;}\right),s.t{\left|\right|r-\mathrm{A\&alpha;}\left|\right|}_2^2\&le;\mathrm{\&zeta;}.$

其中min表示最小化，R^N表示N维空间，

ζ为噪声方差；

(4)按如下迭代加权最小二乘法对上述基于压缩感知的高分辨率超声探测成像模型进行求解，得到一维高分辨率的超声探测成像系数向量α：

(4a)初始化：

k＝0，e_k，e_min＝10^-8，

其中k表示第k次迭代，

表示A的伪逆，α_k表示当前迭代得到的超声成像系数向量，e_k是一个可调参数，初始设置在0～1之间，e_min表示e_k在迭代中的最小值；

(4b)由α_k构造系数的权值矩阵

p是一个0到1之间的自然数，并更新稀疏基矩阵A_k，得到A_k+1：

A_k+1＝AW_k+1

(4c)由A_k+1更新α_k，得到α_k+1：

${\mathrm{\&alpha;}}_{k+1}=W_{k+1}{A}_{k+1}^{+}r$

(4d)判断如下条件是否成立：

${\left|\right|{\mathrm{\&alpha;}}_{k+1}-{\mathrm{\&alpha;}}_k\left|\right|}_2<\sqrt{{\mathrm{\&epsiv;}}_k}/100$

其中

i表示系数向量α_k+1和α_k中第i个元素，α_k+1[i]表示系数向量α_k+1中第i个元素，α_k[i]表示系数向量α_k中第i个元素，

如果不成立，返回步骤(4b)继续迭代；

如果成立，再判断是否满足条件ε_k＜ε_min，如果满足，则停止迭代，得到一维高分辨率的超声探测成像系数向量α＝α_k+1；否则更新ε_k＝ε_k/10，返回步骤(4b)继续迭代，直到满足条件ε_k＜ε_min，得到一维高分辨率的超声探测成像系数向量α＝α_k+1；

(5)将得到的超声探测成像系数α取模值，即得到低频窄带高分辨超声探测成像。

2.根据权利要求1所述的低频窄带高分辨超声探测成像方法，其中步骤(2)涉及的根据不同目标的离散回波信号s(kT_s-τ_i)构建的稀疏基矩阵A，是通过如下公式构建：

A＝[s(kT_s-Δτ)，s(kT_s-2Δτ)K，s(kT_s-nΔτ)，K，s(kT_s-LΔτ)]

其中Δτ＝T_s，T_s表示采样间隔，k＝1，2，L M，M表示采样个数，n＝1，2，L L，L表示矩阵A的列数，L＝τ_max/Δτ，这里τ_max＝max(τ_i)，表示最大延时量，τ_i表示第i个目标的延时量。

Images