CN103994737A - 一种超声辐射力弹性成像的时延估计方法 - Google Patents

Abstract

一种超声辐射力弹性成像的时延估计方法，该方法包括步骤1：分别采集材料受辐射力挤压变形前后的超声射频回波信号；步骤2：对采集的材料受辐射力挤压变形前后的超声射频回波信号进行小波去噪处理；步骤3：利用最大似然互相关时延估计算法，求取去噪后的超声射频回波信号的时间延迟；本发明在低信噪比环境下使得超声辐射力弹性成像的时延估计的抗噪能力得以提高，获得了高抗噪能力，在低信噪比环境下能持续保持良好的估算精度、超声辐射力弹性成像的较高分辨率和鲁棒性。

Description

一种超声辐射力弹性成像的时延估计方法

技术领域

本发明属于超声弹性成像技术领域，具体涉及一种超声辐射力弹性成像的时延估计方法。

背景技术

基于机械挤压变形的超声弹性成像技术已日趋成熟，但是这些技术为接触式检测法，应用范围受到极大限制，检测过程容易对待检材料产生损伤。而基于超声辐射力的弹性成像技术有别于上述接触式检测法，利用超声辐射力作为驱动力源，非接触作用于待检材料上使其产生变形，通过检测变形前后的超声回波信号，并应用相应的变形估计方法对这些信号进行分析处理，得到超声辐射力作用后待检材料变形，从而根据胡克定理确定出待检对象弹性模量相对值的空间分布。超声辐射力弹性成像中信号间的时延估计对最终所获得弹性图谱起到决定性作用，精准的时延估计往往可以对相位畸变校正并提高材料变形的检测精度。

时延估计技术在超声辐射力弹性成像中主要用来处理待检测材料变形前后的超声回波信号，以此获得待检测材料的变形量。超声辐射力弹性成像中的时延估计算法，应用较多的有归一化互相关时延估计法(NC)、绝对差总和时延估计法(SAD)、归一化协方差时延估计法(NCov)、非归一化互相关时延估计法(NNC)、平方差总和时延估计法(SSD)、混合符号互相关时延估计法(HSC)、极性重合时延估计法(PCC)以及Meyr-SpieS时延估计法等。

归一化互相关时延估计法(NC)和归一化协方差法(NCov)由于充分考虑到了材料变形前后超声回波信号的能量。因此，具有较高的时延估计精度，受到了学者们广泛认可和青睐。但是这两种算法计算成本太高、效率较低使其应用具有一定的局限性。

综观各种时延估计算法，评价一种算法的好坏主要从算法效率和算法所能达到的精度进行考虑。目前，各种估计算法中，归一化互相关时延估计法(NC)和归一化协方差法(NCov)精度最高，但效率最低；反之，极性重合时延估计法(PCC)和混合符号互相关时延估计法(HSC)则是精度最低，但效率最高。

发明内容

针对现有技术存在的不足，本发明提供一种超声辐射力弹性成像的时延估计方法。

本发明的技术方案：

一种超声辐射力弹性成像的时延估计方法，包括如下步骤：

步骤1：分别采集材料受辐射力挤压变形前后的超声射频回波信号；

设r1(t)、r2(t)分别为材料受辐射力挤压变形前、后的超声射频回波信号；S(t)为源原信号；n₁(t)、n₂(t)为加性噪声；τ为时间延迟；并假设原发射信号及噪声信号为方差为1、均值为0的正态平稳随机过程，且三者之间没有相关性，并设信道衰减因子α＝1；则r1(t)、r2(t)表示为：

r1(t)＝S(t)+n₁(t)         (1)

r2(t)＝α*S(t-τ)+n₂(t)          (2)

r1(t)、r2(t)的互相关函数表示为：

R_r1r2(τ)＝E[r1(t)*r2(t+τ)]＝R_SS(τ-D)          (3)

式中，R_r1r2(τ)为材料挤压变形前后超声射频回波信号的相关函数；R_SS(τ)为源信号S(t)的自相关函数；D为互相关函数R_r1r2(τ)取得最大值时的时间延迟τ，即 $D=\underset{\mathrm{\τ}}{\arg \max }\left[R_{r1r2}\left(\mathrm{\τ}\right)\right];$

步骤2：对采集的材料受辐射力挤压变形前后的超声射频回波信号进行小波去噪处理；

步骤3：利用最大似然互相关时延估计算法，求取去噪后的超声射频回波信号的时间延迟；

最大似然窗如式(4)所示，

$R_{r_1{r}_2}\left(\mathrm{\τ}\right)=\underset{-\∞}{\overset{\∞}{\∫}}\mathrm{\ψ}\left(f\right)G_{r1r2}\left(f\right)e^{j2\mathrm{\πf\τ}}\mathrm{df}---\left(4\right)$

式中，G_r1r2(f)为材料挤压变形前后超声射频回波信号的互功率谱；ψ(f)为最大似然互相关时延估计算法的加权函数，计算公式为式(5)；

$\mathrm{\ψ}\left(f\right)=\frac{1}{\left|G_{r1r2}\left(f\right)\right|}\frac{{\left|r_{r1r2}\left(f\right)\right|}^2}{{|1-|r_{r1r2}\left(f\right)\left|\right|}^2}---\left(5\right)$

式中，|γ_r1r2(f)|为相关系数，G_r1r1(f)和G_r2r2(f)分别为材料变形前后超声射频回波信号的自功率谱。

步骤2包括如下步骤：

步骤2.1：根据采集的材料受辐射力挤压变形前后的超声射频回波信号，从常用的小波基函数中选择最优小波基函数；利用多种常用小波基函数对超声射频回波信号分别进行分解及重构，并将重构的超声射频回波信号与原始的超声射频回波信号进行互相关，得到一系列相关系数，选取各个常用小波基函数对应的最大的相关系数进行比较，选取其中相关系数最大的小波基函数作为最优小波基函数；

步骤2.2：利用选择的最优小波基函数，对采集的材料受辐射力挤压变形前后的超声射频回波信号进行小波去噪处理。

有益效果：本发明的超声辐射力弹性成像的时延估计方法与现有技术相比较有以下优势：

本发明应用的基于小波变换的最大似然互相关时延估计方法，在低信噪比环境下使得超声辐射力弹性成像的时延估计的抗噪能力得以提高，获得了高抗噪能力，在低信噪比环境下能持续保持良好的估算精度、超声辐射力弹性成像的较高分辨率和鲁棒性；

附图说明

图1为本发明一种实施方式的超声辐射力弹性成像的时延估计方法的流程图；

图2为本发明一种实施方式的材料变形前后的超声射频回波信号仿真波形图；

图3为本发明一种实施方式的加噪后的材料变形前后的超声射频回波信号仿真波形图；

图4为本发明一种实施方式的信噪比为-2至-18dB下四种时延估计方法的超声辐射力弹性成像的时延估计比较结果图；

具体实施方式

下面结合附图对本发明的一种实施方式作详细说明。

本实施方式所使用的待检测弹性材料选用自制的与琼脂组织弹性性质相仿的凝胶仿体，采用反复冷冻-解冻的方法来制备聚乙烯醇PVA水凝胶生物仿体。制备生物仿体时，采用的冷冻温度为-20℃，冷冻时间为12h，经过冷冻之后，将生物仿体置于温度为20℃下解冻12h，即，一个冷冻-解冻周期为24h，通过控制冷冻-解冻的循环次数获得不同硬度的聚乙烯醇PVA水凝胶生物仿体。单纯的由聚乙烯醇PVA聚合物和水溶液混合制成的PVA水凝胶，呈现半透明状，为了能够在超声下成像，必须加入散射子来散射超声信号。本实施方式采用的聚乙烯醇PVA聚合物是Polyvinylalkohol聚乙烯醇，散射子使用Sigmacell纤维素。本实施方式中需要使用两种不同硬度的聚乙烯醇PVA水凝胶，其中，外围仿体为5*5*103cm的长方体，由8％浓度的Polyvinylalkohol聚乙烯醇和1％浓度的Sigmacell纤维素的水凝胶制备而成，杨氏模量大约为30kPa。内容物为直径为10mm的圆柱体，由10％浓度的Polyvinylalkohol聚乙烯醇和1％浓度的Sigmacell纤维素的水凝胶制备而成，其杨氏模量的值大约为150kPa，该内容物插在外围仿体的中央即25mm处。

本实施方式的超声辐射力弹性成像的时延估计方法，如图1所示，该方法开始于步骤101，

步骤102：建立材料受辐射力挤压变形前后的超声射频回波信号模型函数；

步骤1：分别采集材料受辐射力挤压变形前后的超声射频回波信号；

设r1(t)、r2(t)分别为材料受辐射力挤压变形前、后的超声射频回波信号；S(t)为源原信号；n₁(t)、n₂(t)为加性噪声；τ为时间延迟；并假设原发射信号及噪声信号为方差为1、均值为0的正态平稳随机过程，且三者之间没有相关性，并设信道衰减因子α＝1；则r1(t)、r2(t)表示为：

r1(t)＝S(t)+n₁(t)          (1)

r2(t)＝α*S(t-τ)+n₂(t)           (2)

r1(t)、r2(t)的互相关函数表示为：

R_r1r2(τ)＝E[r1(t)*r2(t+τ)]＝R_SS(τ-D)          (3)

式中，R_r1r2(τ)为材料挤压变形前后超声射频回波信号的相关函数；R_SS(τ)为源信号S(t)的自相关函数；D为互相关函数R_r1r2(τ)取得最大值时的时间延迟τ，即 $D=\underset{\mathrm{\τ}}{\arg \max }\left[R_{r1r2}\left(\mathrm{\τ}\right)\right];$

时延估计技术在超声声辐射力弹性成像中主要用来处理待检测材料变形前后的超声回波信号，以此获得待检测材料的变形量。还可以将待检测材料挤压变形前后的回波信号r1(t)、r2(t)分别表示为s_r(t)、s_d(t)，而s_r(t)、s_d(t)分别用式(6)和式(7)进行表示，以便于建立材料受辐射力挤压变形前后的超声射频回波信号模型函数：

s_r(t)＝s₁(t)*h₁(t)+n′₁(t)           (6)

s_d(t)＝s₂(t)*h₂(t)+n′₂(t)            (7)

式中，s_r(t)和s_d(t)是两个随机过程，s_r(t)为材料受辐射力挤压变形前的回波信号；t为实时时间，s；s_d(t)为材料受辐射力挤压变形后的回波信号；s₁(t)和s₂(t)分别为超声换能器的发射信号，s₂(t)较s₁(t)有一个时间延迟；h₁(t)和h₂(t)分别为待检测材料挤压变形前后的冲击响应函数；n′₁(t)和n′₂(t)分别为待检测材料挤压变形前后的噪声信号，包括：材料噪声、周边环境噪声及电子噪声。

为了分析噪声的信噪比大小、频率、相关度以及带宽等因素对时延估计算法的影响，建立材料受辐射力挤压变形前后的超声射频回波信号模型函数分别为式(8)和(9)，即超声回波信号可以利用合适的组合基信号来表示。在s₁(t)和s₂(t)之间存在一个时间延迟，为了去除无偏估计的损失，所以往往在评估时延估计算法的性能时通常假定时间延迟为0。

s_r(t)＝s_b1(t)+(1/C)·s_b3(t)        (8)

s_d(t)＝W_1,2s_b1(t)+W_2,2s_b2(t)(1/C)·s_b4(t)         (9)

式中，s_b1(t)，s_b2(t)，s_b3(t)，s_b4(t)分别为超声射频回波基信号，利用sinc函数包络的正弦信号与高斯噪声相卷积产生，sinc函数包络的正弦信号(即点扩散函数)如式(10)所示；C为超声射频回波信号中噪声的信噪比大小，dB；W_i,j是对称矩阵W通过乔里斯基(Cholesky)算法分解后得到的矩阵中的元素；矩阵W是通过式(11)计算得到；

$\mathrm{psf}\left(t\right)=\frac{\sin \left(\mathrm{\πB}{f}_{0}t\right)}{\mathrm{\πB}{f}_{0}t}\sin \left(\mathrm{\πB}{f}_{0}t\right)---\left(10\right)$

式中，B为超声射频回波信号的分数带宽，比特率；f₀为超声射频回波信号的中心频率，赫兹；目前，已经有一些点扩展函数可以用来模拟真实的超声回波信号，但是因为采用sinc函数包络的正弦信号所获得点扩展函数psf(t)的频谱为平带有限谱，与超声回波的实际频谱较为一致。所以，将sinc函数包络的正弦信号作为超声回波信号的点扩展函数得到普遍认可。

$W\·W^T=\left[\begin{array}{cc}1& R\\ R& 1\end{array}\right]---\left(11\right)$

式中，W^T为矩阵W的转置；R为超声射频回波信号期望的相关系数；s_r(t)和s_d(t)通常利用插值的方法使超声回波信号过采样至中心频率的200倍，从而减小由于采样过低带来的误差影响，最终达到提高时延估计的目的。

步骤103：根据材料受辐射力挤压变形前后的超声射频回波信号模型函数，进行材料受辐射力挤压变形前后的超声射频回波信号的仿真,得到材料变形前后的超声射频回波信号；

本实施方式利用Field_II仿真软件在Matlab中进行材料变形前后的超声射频回波信号的仿真。

为了使仿真出的超声射频回波信号能更加贴近实际采集的超声射频回波信号，仿真时设置换能器中心频率为5MHZ,采样率为100MHZ，超声在待检测材料的传播速度为1540m/s。仿真中我们设置材料变形前后的超声射频回波信号的时间延迟为10T，T代表超声射频回波信号的周期，这里为换能器中心频率的倒数，仿真后得到的材料变形前后的超声射频回波信号，如图2所示，其中1为材料变形前的超声射频回波信号的仿真波形，2为材料变形后的超声射频回波信号的仿真波形。本发明是基于一种低信噪比环境下仍能保持良好特性的时延估计方法，所以本实施方式在仿真时，要对不同信噪比下的仿真信号进行处理分析，给超声射频回波信号加的高斯白噪声信号比为-18到-2分贝，加噪后仿真得到的材料变形前后的超声射频回波信号，如图3所示。

步骤104：对仿真得到的材料变形前后的超声射频回波信号进行小波去噪处理；

首先，根据仿真得到的材料变形前后的超声射频回波信号,从常用的小波基函数中选择最优小波基函数；

本实施方式选取了五种常用的小波基函数作为待选择对象，分别为db小波、dmey小波、coif小波、bior小波和sym小波；利用该五种常用的小波基函数对仿真得到的信噪比-18分贝的超声射频回波信号分别进行了分解及重构。根据采样定理，要从抽样信号中从失真的恢复原信号，采样频率应大于2倍信号的最高频率。所以仿真信号的采样频率为100e6HZ，探头固有频率为5e6HZ，所以重构层数应选择第三层0～12e6HZ。

然后将重构信号与原始的超声回波信号进行互相关，得到一系列相关系数，选取该五个小波基函数对应的最大的相关系数进行比较，选取其中相关系数最大的小波基函数为本次试验的最优小波基函数。从表1可以看出，Coif小波在消失矩为3的时候相关系数最大，所以本实施方式选取Coif小波作为本实施方式的最优小波基函数。

表1五种小波基函数的最大相关系数

小波基函数名称	Db	dmey	Coif	Bior	Sym
						消失矩	7	空	3	6	8
最大相关系数	0.9271	0.8032	0.9308	0.9174	0.9199

最后，利用选择的最优小波基函数，对仿真得到的材料变形前后的超声射频回波信号进行小波去噪处理；

即，利用本实施方式选择的最优小波基函数Coif小波对仿真得到的材料变形前后的超声射频回波信号进行小波去噪处理。

步骤105：利用最大似然互相关时延估计算法，求取去噪后的超声射频回波信号的时间延迟；

最大似然窗如式(4)所示，

$R_{r_1{r}_2}\left(\mathrm{\τ}\right)=\underset{-\∞}{\overset{\∞}{\∫}}\mathrm{\ψ}\left(f\right)G_{r1r2}\left(f\right)e^{j2\mathrm{\πf\τ}}\mathrm{df}---\left(4\right)$

式中，G_r1r2(f)为材料挤压变形前后超声射频回波信号的互功率谱；ψ(f)为最大似然互相关时延估计算法的加权函数，计算公式为式(5)；

$\mathrm{\ψ}\left(f\right)=\frac{1}{\left|G_{r1r2}\left(f\right)\right|}\frac{{\left|r_{r1r2}\left(f\right)\right|}^2}{{|1-|r_{r1r2}\left(f\right)\left|\right|}^2}---\left(5\right)$

式中，|γ_r1r2(f)|为相关系数，G_r1r1(f)和G_r2r2(f)分别为材料变形前后超声射频回波信号的自功率谱。

如图4所示，可以看出不同的信噪比下，从普通互相关时延估计法的时延估计结果3、广义互相关时延估计法的时延估计结果4、最大似然四种时延估计相关时延估计法的时延估计结果5和基于小波变换的最大似然互相关时延估计方法的时延估计结果6，可以很清晰的看出，在信噪比大于-11分贝的情况下，四种时延估计方法都保持着良好的估计能力，能准确计算出时间延迟，而且信噪比渐渐降低时，普通互相关时延估计方法和广义互相关时延估计算法的估计精度迅速降低，而在信噪比大于-13分贝的情况下，最大似然互相关时延估计方法和基于小波变换的最大似然互相关时延估计方法仍保持着良好的精度，信噪比再降低时，最大似然互相关时延估计方法失去的估算精度，只有基于小波变换的最大似然互相关时延估计方法仍继续保持的良好的精度。

步骤106：结束。

虽然以上描述了本发明的具体实施方式，但是本领域内的熟练的技术人员应当理解，这些仅是举例说明，可以对这些实施方式做出多种变更或修改，而不背离本发明的原理和实质。本发明的范围仅由所附权利要求书限定。

Claims (3)

1.一种超声辐射力弹性成像的时延估计方法，其特征在于，包括如下步骤：

步骤1：分别采集材料受辐射力挤压变形前后的超声射频回波信号；

设r1(t)、r2(t)分别为材料受辐射力挤压变形前、后的超声射频回波信号；S(t)为源原信号；n₁(t)、n₂(t)为加性噪声；τ为时间延迟；并假设原发射信号及噪声信号为方差为1、均值为0的正态平稳随机过程，且三者之间没有相关性，并设信道衰减因子α＝1；则r1(t)、r2(t)表示为：

r1(t)＝S(t)+n₁(t)          (1)

r2(t)＝α*S(t-τ)+n₂(t)          (2)

r1(t)、r2(t)的互相关函数表示为：

R_r1r2(τ)＝E[r1(t)*r2(t+τ)]＝R_SS(τ-D)        (3)

式中，R_r1r2(τ)为材料挤压变形前后超声射频回波信号的相关函数；R_SS(τ)为源信号S(t)的自相关函数；D为互相关函数R_r1r2(τ)取得最大值时的时间延迟τ，即 $D=\underset{\mathrm{\τ}}{\arg \max }\left[R_{r1r2}\left(\mathrm{\τ}\right)\right];$

步骤2：对采集的材料受辐射力挤压变形前后的超声射频回波信号进行小波去噪处理；

步骤3：利用最大似然互相关时延估计算法，求取去噪后的超声射频回波信号的时间延迟；

最大似然窗如式(4)所示，

$R_{r_1{r}_2}\left(\mathrm{\τ}\right)=\underset{-\∞}{\overset{\∞}{\∫}}\mathrm{\ψ}\left(f\right)G_{r1r2}\left(f\right)e^{j2\mathrm{\πf\τ}}\mathrm{df}---\left(4\right)$

式中，G_r1r2(f)为材料挤压变形前后超声射频回波信号的互功率谱；ψ(f)为最大似然互相关时延估计算法的加权函数，计算公式为式(5)；

$\mathrm{\ψ}\left(f\right)=\frac{1}{\left|G_{r1r2}\left(f\right)\right|}\frac{{\left|r_{r1r2}\left(f\right)\right|}^2}{{|1-|r_{r1r2}\left(f\right)\left|\right|}^2}---\left(5\right)$

式中，|γ_r1r2(f)|为相关系数，G_r1r1(f)和G_r2r2(f)分别为材料变形前后超声射频回波信号的自功率谱。

2.根据权利要求1所述的超声辐射力弹性成像的时延估计方法，其特征在于，步骤2包括如下步骤：

步骤2.1：根据采集的材料受辐射力挤压变形前后的超声射频回波信号，从常用的小波基函数中选择最优小波基函数；

步骤2.2：利用选择的最优小波基函数，对采集的材料受辐射力挤压变形前后的超声射频回波信号进行小波去噪处理。

3.根据权利要求2所述的超声辐射力弹性成像的时延估计方法，其特征在于，步骤2.1选择最优小波基函数的方法为：

利用多种常用小波基函数对超声射频回波信号分别进行分解及重构，并将重构的超声射频回波信号与原始的超声射频回波信号进行互相关，得到一系列相关系数，选取各个常用小波基函数对应的最大的相关系数进行比较，选取其中相关系数最大的小波基函数作为最优小波基函数。