CN106063711A - 一种骨导波检测与分离、识别及重构方法 - Google Patents

Abstract

本发明提供一种骨导波检测与分离、识别及重构方法：设计制作了开放式“发射‑气隙‑接收”骨导波远端检测装置；利用该骨导波远端检测装置得到多模骨导波；联合平滑伪维纳时频分布、Rayleigh–Lamb频散曲线以及自适应窗宽的高斯函数平滑方法对检测所得多模骨导波进行模式分离识别；利用分离识别得到的各单模骨导波对骨表入射点导波进行重构，有效避免了入射波、反射波对骨导波检测的混叠干扰。

Description

一种骨导波检测与分离、识别及重构方法

技术领域

本发明涉及超声造影成像技术领域，特别涉及基于开放式“发射-气隙-接收”骨导波检测装置与分离、识别及重构方法。

背景技术

超声成像已广泛用于骨表肌肉组织系统成像，特别对肌骨系统肿瘤、骨表皮瓣的诊断监控。当超声作用于骨表时，骨表皮质将产生导波，这种骨导波已经引起研究者广泛关注并用于骨质定量检测。同时，骨导波向周边软组织层泄漏，这将干扰软组织回波信号的信噪比，进而可能影响到骨表图像质量。骨导波因其非线性传播而具有多模式频散特性，并已用于检测或操控纳米液滴、无包膜微泡和亚微粒子。但是这种检测仅利用了导波的线性特征，而其非线性频散特性的作用被忽略了。但在实际研究骨导波的非线性频散特性时，需要提取排除其他杂波干扰的骨导波，然而骨导波中往往包含着入射波和反射波的混叠干扰，这使得提取骨导波十分困难。

同时，骨导波的分离方法主要有各类改进的FFT变换、小波多尺度分析等模式分离方法，如2D-FFT、S-FFT等。2D-FFT、S-FFT等方法可快速提取导波第一接收波(FAS)速度，但此类方法更适用于具有固定带宽比的非平稳信号分析，对于多阶次频变非平稳导波信号的分析，严重受到时频交叉项的干扰、低信噪比且对高阶导波模式因其在骨皮质的快速衰减而难以准确分离。为此，Kailiang X.等将锥形核算子引入Wigner-Viler分布形成Zhao-Atlas-Marks分布。另外，非线性调频小波变换(Chirplet)可准确分离各频率分量、有效抑制混叠干扰与虚假频率分量，但少有肌骨系统导波分析的应用案例。

发明内容

针对在提取骨导波过程中入射波与反射波的混叠干扰难以滤除的困难，本发明提出了一种骨导波检测与分离、识别及重构方法。

为实现上述目的，本发明采用了以下技术方案：

(1)利于基于发射-气隙-接收模型的装置将骨导波检测中的入射波和反射波限制在发射区；

(2)联合平滑伪维纳时频分布、Rayleigh–Lamb频散曲线以及自适应窗宽的高斯函数平滑方法对单模骨导波进行分离识别；

(3)经过步骤(2)后，利用各单模骨导波的声压波形、衰减系数和传播速度对入射点骨导波进行重构。

所述基于发射-气隙-接收模型的装置包括两个之间留有气隙的水箱，两个水箱相对侧的壁面上设置有孔，两孔上覆盖有橡胶片，两橡胶片上设置有用于供骨样本穿入对应侧水箱的孔，该孔内设置有密封固定部件；两个水箱的内壁设置有吸声材料层；其中一个水箱内设置有用于发射入射角为θ、发射频率为f的入射波的线阵探头，另一个水箱内设置有接收角为的水听器。

所述气隙厚度<3mm。

所述步骤(2)具体包括以下步骤：

(2.1)调整入射角θ和接收角以及骨导波传播距离d，在骨导波检测中对混叠的多模骨导波在时域进行分离；

(2.2)对多模骨导波时域波形进行SPWV时频分析，得到多模骨导波SPWV时间-频率能量分布；

(2.3)依据Rayleigh–Lamb频散方程，计算测试骨样本的各单模导波的频厚积-相速度频散曲线；

(2.4)依据测试骨样本皮质材料信息、入射波频率f和各单模骨导波的频厚积-相速度频散曲线，计算各单模骨导波的频厚积-群速度频散曲线；

(2.5)将频厚积-群速度频散曲线进行坐标翻转变换得到群速度-频率频散曲线，将群速度-频率频散曲线与所得的多模骨导波SPWV时间-频率能量分布一一耦合匹配，从而识别各单模骨导波的对称与反对称模式；

(2.6)若多模骨导波在时域能够分离，则选择窗宽与各单模骨导波持续时间相适应的高斯函数进行移动平滑点乘，按到达检测点时间先后顺序分离各单模骨导波声压波形p。

若多模骨导波在时域难以分离，则依据SPWV时间-频率能量分布，确定各单模骨导波的中心频率与带宽，依次对各单模骨导波进行带通频域滤波以分离各单模骨导波声压波形p。

所述步骤(3)具体包括以下步骤：

(3.1)改变骨导波传播距离d，测得以下模式为j的单模骨导波的非线性传播参数，包括：衰减系数α_j、传播速度v_j：

${\mathrm{\&alpha;}}_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{V_{jk}-V_{j1}}{d_{jk}-d_{j1}}/(N-1)$

$v_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{d_{jk}-d_{j1}}{t_{jk}-t_{j1}}/(N-1)$

其中，N为测试次数，d_jk为在第k个检测点测得的模式为j的单模骨导波传播距离；V_jk为在第k个检测点测得的模式为j的单模骨导波的峰值强度；t_jk为在第k个检测点测得的模式为j的单模骨导波的到达时间，其中到达时间t为单模骨导波首次幅值超过5％峰值强度对应的时刻；

(3.2)依据模式为j的单模骨导波的声压波形p_j、衰减系数α_j、传播速度v_j与骨导波传播距离d_j，重构出入射点i处各单模骨导波声压波形P_ji：

P_ji(t-d_j/v_j)＝(1+α_jd_j/Δp_j(t))×p_j(t-d_j/v_j)

其中，Δp_j为单模骨导波声压波形p_j经传播距离d_j后的幅值改变；

(3.3)依据重构得到的入射点各单模骨导波声压波形P_ji，重构出入射点i处多模骨导波的声压波形P_Li。

$\begin{array}{c}{P}_{Li}\left(t\right)=\underset{j=1}{\overset{n}{\&Sigma;}}{P}_{ji}(t-d_j/v_j)\\ =\underset{j=1}{\overset{n}{\&Sigma;}}\&lsqb;(1+{\mathrm{\&alpha;}}_j{d}_j/{\mathrm{\&Delta;p}}_j(t\left)\right)\&times;p_j(t-d_j/v_j)\&rsqb;\end{array}$ 。

本发明的优点如下：

1、设计了一种骨导波的远端检测装置，用一空气隙隔离了入射波、反射波的混叠干扰，装置简单，开放式的多点检测方式操作简便。

2、提出了骨导波分离、识别与重构的成套方法，准确快捷地对骨导波所含有的诸多单模式导波进行分离、识别并重构，为骨导波非线性频散特性的研究提供了提取骨导波的方法。

附图说明

图1为基于发射-气隙-接收模型的骨导波检测实验原理框图；

图2为发射-气隙-接收模型的结构图，其中，(a)为模型正视图，(b)为同径环状凸体侧视图；

图3为基于发射-气隙-接收模型测得的骨导波；

图4为基于SPWV时频分析-频散曲线的单模骨导波单模式分离识别；

图5为基于单模骨导波的水-骨表界面入射点的导波重构波形；

图中：1为左侧水箱，2为气隙，3为右侧水箱，4为同径环状凸体，5为骨样本，6为支撑固定部件，7为吸声材料，8为排水口，9为凹槽，10为橡胶薄片，11为密封固定部件。

具体实施方式

下面结合附图和实施例对本发明做详细叙述。

本发明主要提出在入射波和反射波混叠干扰情况下，提取出骨导波，其关键在于在排除入射波和反射波干扰的同时，提取出骨表纯导波。为此，本发明设计了基于发射-气隙-接收模型的实验装置，将入射波和反射波限制在发射区，并在接收检测区检测多模骨导波，通过联合SPWV时频分布、Rayleigh–Lamb频散曲线以及自适应窗宽的高斯函数平滑方法对检测到的多模骨导波进行分离、识别，最后利用各单模骨导波相关信息对入射点骨导波进行重构。

步骤一：设计基于模块化设备的开放式“发射-气隙-接收”骨导波远端检测装置，其结构部件如下，参见图1：

(1)该装置的主体包括：左右两侧水箱以及两水箱之间的窄小空气隙2(厚度<3mm)，参见图2；

(2)左侧水箱1的右壁，以及右侧水箱3的左壁均按最大可测试骨样本直径开孔；

(3)所述孔向箱壁两侧设有含一凹槽9的同径环状凸体4(与孔直径相同)，松弛的橡胶薄片10覆盖在孔表面，橡胶薄片10四周被固定在凹槽9槽体中，测试骨样本5贯穿两橡胶薄片，并被两侧水箱内的支撑固定部件6固定；

(4)两侧水箱内壁固定了一定厚度(3mm～10mm)的吸声材料7，测试骨样本5被纯水淹没；

(5)左侧水箱1为发射区，Sonix-Touch的线阵探头以一入射角θ发射频率为f的入射波，在该侧水箱激励测试骨样本产生骨导波；

(6)右侧水箱3为骨导波接收检测区，直径为0.5mm的针式水听器(PrecisionAcoustics Ltd.,Dorchester,UK)以一接收角在该侧水箱检测从骨表传来的骨导波，接收检测区的针式水听器与发射区的线阵探头的距离为d(即骨导波传播距离)；

(7)骨导波的时域波形经过辅助放大器放大后被数据采集卡CS14100(GageApplied,Inc.Lachine,QC,Canada)采集记录；

(8)入射波发射与骨导波接收同时被数字化超声成像系统Sonix-Touch同步控制；

(9)入射角θ和接收角由角度仪调控，骨导波传播距离d由三维扫描系统5800PR调控。依据Snell定律，调整入射角θ、接收角以及骨导波传播距离d，检测得到在时域基本分离的多模骨导波。

步骤二：联合平滑伪维纳(smoothed-pseudo Wigner–Ville，SPWV)时频分布、Rayleigh–Lamb频散曲线以及自适应窗宽的高斯函数平滑方法对单模骨导波进行分离识别，其步骤如下：

(1)对检测得到的多模骨导波时域波形进行SPWV时频分析，得到多模骨导波的SPWV时间-频率能量分布；

(2)依据Rayleigh–Lamb频散方程，计算测试骨样本的各单模导波的频厚积-相速度频散曲线；

描述导波相速度的方程即为Rayleigh-Lamb频散方程：

对称(S)模式

$\frac{{\mathrm{tank}}_{s}b}{{\mathrm{tank}}_{l}b}=-\frac{4k_0^2{k}_l{k}_s}{{(k_0^2-k_s^2)}^2}$

反对称(A)模式

$\frac{{\mathrm{tank}}_{s}b}{{\mathrm{tank}}_{l}b}=-\frac{{(k_0^2-k_s^2)}^2}{4k_0^2{k}_l{k}_s}$

其中，

$k_l^2={\left(\frac{\mathrm{\&omega;}}{c_l}\right)}^2-k_0^2,k_s^2={\left(\frac{\mathrm{\&omega;}}{c_s}\right)}^2-k_0^2$

式中：k₀为沿骨表水平方向的波数，b为1/2骨表厚，ω为角频率，ω＝2πf，c_l为纵波速度，c_s为横波速度，相速度c_p＝ω/k₀；

(3)依据测试骨样本皮质材料信息(厚度，泊松比，以及杨氏模量)、入射波频率f和各单模骨导波的频厚积-相速度频散曲线，计算骨表各单模导波的频厚积-群速度频散曲线；

群速度可表示为：

c_g＝dω/dk₀

将k₀＝ω/c_p代入上式，有

$c_g=d\mathrm{\&omega;}{\&lsqb;d\left(\frac{\mathrm{\&omega;}}{c_p}\right)\&rsqb;}^{-1}={\left(c_p\right)}^2{\&lsqb;c_p-\mathrm{\&omega;}f\frac{{\mathrm{dc}}_p}{d\mathrm{\&omega;}}\&rsqb;}^{-1}$

式中：k₀为沿板水平方向的波数，ω为角频率，ω＝2πf，c_p为相速度；

(4)将频厚积-群速度频散曲线进行坐标翻转变换得到群速度-频率频散曲线，与测得各单模骨导波的SPWV时间-频率能量分布一一耦合匹配，准确识别各单模骨导波对称(S)与反对称(A)模式；

(5)如检测所得多模骨导波在时域基本分离，则选择窗宽与各单模骨导波持续时间相适应的高斯函数进行移动平滑点乘，按到达检测点时间先后顺序准确分离各单模骨导波声压波形p；如检测所得多模骨导波在时域难以分离，则依据SPWV时间-频率能量分布，确定各单模骨导波的中心频率与带宽，依次对各单模骨导波进行带通频域滤波以分离各单模骨导波声压波形p。

步骤三：利用各单模骨导波相关信息对入射点骨导波进行重构，其步骤依次如下：

(1)在开放式“发射-气隙-接收”骨导波远端检测装置中，改变骨导波传播距离d，测得的模式为j的骨表单模导波的非线性传播参数，包括：衰减系数α_j、传播速度v_j；计算公式如下：

${\mathrm{\&alpha;}}_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{V_{jk}-V_{j1}}{d_{jk}-d_{j1}}/(N-1)$

$v_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{d_{jk}-d_{j1}}{t_{jk}-t_{j1}}/(N-1)$

其中，N为测试次数，d_jk为在第k个检测点测得的模式为j的骨表单模导波传播距离；V_jk为在第k个检测点测得的模式为j的骨表单模导波的峰值强度；t_jk为在第k个检测点测得的模式为j的骨表单模导波的到达时间。其中到达时间t为骨表单模导波首次幅值超过5％峰值强度对应的时刻。

(2)依据检测点测得模式为j的单模骨导波的声压波形p_j、衰减系数α_j、传播速度v_j与距离d_j，重构出入射点i处各单模骨导波声压波形P_ji；计算公式如下：

P_ji(t-d_j/v_j)＝(1+α_jd_j/Δp_j(t))×p_j(t-d_j/v_j)

其中，Δp_j为单模骨导波声压波形p_j经传播距离d_j后的幅值改变。

(3)依据重构得到的入射点各单模骨导波声压波形P_ji，重构出入射点i处多模骨导波的声压波形P_Li

$\begin{array}{c}{P}_{Li}\left(t\right)=\underset{j=1}{\overset{n}{\&Sigma;}}{P}_{ji}(t-d_j/v_j)\\ =\underset{j=1}{\overset{n}{\&Sigma;}}\&lsqb;(1+{\mathrm{\&alpha;}}_j{d}_j/{\mathrm{\&Delta;p}}_j(t\left)\right)\&times;p_j(t-d_j/v_j)\&rsqb;\end{array}$

n表示模式数量。

应用举例

选择牛胫骨皮质(厚度2mm，泊松比0.34，杨氏模量19.76GPa)作为入射波激励对象，即骨样本，在接收端获得骨导波波形。橡胶薄片10固定于同径环状凸体4上的环状凹槽9，橡胶薄片10的中间有一小孔，使得骨样本可以穿过，小孔内设置密封固定部件11，具有使骨样本与橡胶薄片紧密结合并阻止水箱中的纯水进入气隙2的作用。支撑固定部件6具有支撑固定骨样本并使骨样本两端保持水平的作用，水箱内壁粘贴的吸声材料7为吸声海绵。依据Snell定律，调整入射角θ为34°、接收角为27°以及骨导波传播距离d为30mm，检测得到在时域基本分离的多模骨导波。检测到的接收端的骨导波时域波形如图3所示。本发明涉及算法均在matlab平台上编程实现，各个模式骨导波SPWV能量分布与其频散曲线的匹配识别如图4所示，可以看出，按时间顺序，骨导波模式分别为S0和A1。最终，根据各单模骨导波重构出入射点骨导波，如图5所示。从而实现了入射点骨导波的重构提取。

Claims (6)

1.一种骨导波检测与分离、识别及重构方法，其特征在于：包括以下步骤：

(1)利于基于发射-气隙-接收模型的装置将骨导波检测中的入射波和反射波限制在发射区；

(2)联合平滑伪维纳时频分布、Rayleigh–Lamb频散曲线以及自适应窗宽的高斯函数平滑方法对单模骨导波进行分离识别；

(3)经过步骤(2)后，利用各单模骨导波的声压波形、衰减系数和传播速度对入射点骨导波进行重构。

2.根据权利要求1所述一种骨导波检测与分离、识别及重构方法，其特征在于：所述基于发射-气隙-接收模型的装置包括两个之间留有气隙(2)的水箱，两个水箱相对侧的壁面上设置有孔，两孔上覆盖有橡胶片，两橡胶片上设置有用于供骨样本(5)穿入对应侧水箱的孔，该孔内设置有密封固定部件(11)；两个水箱的内壁设置有吸声材料层；其中一个水箱内设置有用于发射入射角为θ、发射频率为f的入射波的线阵探头，另一个水箱内设置有接收角为的水听器。

3.根据权利要求2所述一种骨导波检测与分离、识别及重构方法，其特征在于：所述气隙(2)厚度<3mm。

4.根据权利要求1所述一种骨导波检测与分离、识别及重构方法，其特征在于：所述步骤(2)具体包括以下步骤：

(2.1)调整入射角θ和接收角以及骨导波传播距离d，在骨导波检测中对混叠的多模骨导波在时域进行分离；

(2.2)对多模骨导波时域波形进行SPWV时频分析，得到多模骨导波SPWV时间-频率能量分布；

(2.3)依据Rayleigh–Lamb频散方程，计算测试骨样本的各单模导波的频厚积-相速度频散曲线；

(2.4)依据测试骨样本皮质材料信息、入射波频率f和各单模骨导波的频厚积-相速度频散曲线，计算各单模骨导波的频厚积-群速度频散曲线；

(2.5)将频厚积-群速度频散曲线进行坐标翻转变换得到群速度-频率频散曲线，将群速度-频率频散曲线与所得的多模骨导波SPWV时间-频率能量分布一一耦合匹配，从而识别各单模骨导波的对称与反对称模式；

(2.6)若多模骨导波在时域能够分离，则选择窗宽与各单模骨导波持续时间相适应的高斯函数进行移动平滑点乘，按到达检测点时间先后顺序分离各单模骨导波声压波形p。

5.根据权利要求4所述一种骨导波检测与分离、识别及重构方法，其特征在于：若多模骨导波在时域难以分离，则依据SPWV时间-频率能量分布，确定各单模骨导波的中心频率与带宽，依次对各单模骨导波进行带通频域滤波以分离各单模骨导波声压波形p。

6.根据权利要求1所述一种骨导波检测与分离、识别及重构方法，其特征在于：所述步骤(3)具体包括以下步骤：

(3.1)改变骨导波传播距离d，测得以下模式为j的单模骨导波的非线性传播参数，包括：衰减系数α_j、传播速度v_j：

${\mathrm{\&alpha;}}_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{V_{jk}-V_{j1}}{d_{jk}-d_{j1}}/(N-1)$

$v_j=\underset{k=2}{\overset{N}{\&Sigma;}}\frac{d_{jk}-d_{j1}}{t_{jk}-t_{j1}}/(N-1)$

其中，N为测试次数，d_jk为在第k个检测点测得的模式为j的单模骨导波传播距离；V_jk为在第k个检测点测得的模式为j的单模骨导波的峰值强度；t_jk为在第k个检测点测得的模式为j的单模骨导波的到达时间，其中到达时间t为单模骨导波首次幅值超过5％峰值强度对应的时刻；

(3.2)依据模式为j的单模骨导波的声压波形p_j、衰减系数α_j、传播速度v_j与骨导波传播距离d_j，重构出入射点i处各单模骨导波声压波形P_ji：

P_ji(t-d_j/v_j)＝(1+α_jd_j/Δp_j(t))×p_j(t-d_j/v_j)

其中，Δp_j为单模骨导波声压波形p_j经传播距离d_j后的幅值改变；

(3.3)依据重构得到的入射点各单模骨导波声压波形P_ji，重构出入射点i处多模骨导波的声压波形P_Li。

$\begin{array}{c}{P}_{Li}\left(t\right)=\underset{j=1}{\overset{n}{\&Sigma;}}{P}_{ji}(t-d_j/v_j)\\ =\underset{j=1}{\overset{n}{\&Sigma;}}\&lsqb;(1+{\mathrm{\&alpha;}}_j{d}_j/{\mathrm{\&Delta;p}}_j(t\left)\right)\&times;p_j(t-d_j/v_j)\&rsqb;\end{array}.$