CN110501423B - 一种基于频域分段的高分辨率最小方差超声成像方法 - Google Patents

一种基于频域分段的高分辨率最小方差超声成像方法 Download PDF

Info

Publication number
CN110501423B
CN110501423B CN201910757165.6A CN201910757165A CN110501423B CN 110501423 B CN110501423 B CN 110501423B CN 201910757165 A CN201910757165 A CN 201910757165A CN 110501423 B CN110501423 B CN 110501423B
Authority
CN
China
Prior art keywords
frequency domain
signal
sub
window function
window
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910757165.6A
Other languages
English (en)
Other versions
CN110501423A (zh
Inventor
王平
杜婷婷
王林泓
孔露
李锡涛
柳学功
孔美娅
田训
梁家祺
王慧悦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chongqing Caojie Shipping Power Development Co ltd
Chongqing Mostag Energy Management Co ltd
Chongqing University
Chongqing College of Electronic Engineering
Original Assignee
Chongqing Mostag Energy Management Co ltd
Chongqing University
Chongqing College of Electronic Engineering
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chongqing Mostag Energy Management Co ltd, Chongqing University, Chongqing College of Electronic Engineering filed Critical Chongqing Mostag Energy Management Co ltd
Priority to CN201910757165.6A priority Critical patent/CN110501423B/zh
Publication of CN110501423A publication Critical patent/CN110501423A/zh
Application granted granted Critical
Publication of CN110501423B publication Critical patent/CN110501423B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/02—Analysing fluids
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
    • G01N29/0654—Imaging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4454—Signal recognition, e.g. specific values or portions, signal events, signatures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/04—Wave modes and trajectories
    • G01N2291/044—Internal reflections (echoes), e.g. on walls or defects

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Acoustics & Sound (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Signal Processing (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

本发明涉及一种基于频域分段的高分辨率最小方差超声成像方法，属于超声成像领域。该方法首先对阵元接收的采样信号进行延时处理，获得超声聚焦所需的回波数据；其次根据STFT中自适应窗函数的最大集中度测量准则，选取频域分段最优窗函数；利用STFT将超声回波信号转换为窄带子信号；利用共轭对称性，前一半窄带子信号经过共轭对称处理生成另一半窄带信号；将接收阵列依次划分为具有重叠阵元的子阵，对频域信号进行前后向平滑和对角加载处理，获得样本协方差矩阵；最后利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号。该方法可以显著提升超声成像分辨率，提高对比度，可以在整体上提高超声成像的质量。

Description

一种基于频域分段的高分辨率最小方差超声成像方法
技术领域
本发明属于超声成像技术领域，涉及一种基于频域分段的高分辨率最小方差超声成像方法。
背景技术
超声成像以其安全性、无创性、实时性和低成本等优点被广泛应用于无损检测与诊断领域，而波束形成技术是超声成像的关键技术，直接决定超声成像的图像质量。延时叠加算法(DelayAnd Sum，DAS)，是超声成像中应用最为广泛的，也是最简单的波束形成技术。它根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算，然后对延时后的数据对齐叠加。传统DAS算法复杂度低，稳健性好，成像速度快，但由于其采用固定窗函数加权导致主瓣宽度增加，因此分辨率和对比度较低。
近年来，为了提高波束形成算法的对比度和分辨率，自适应算法得到越来越多的研究。1969年Capon提出的最小方差(MinimumVariance,MV)波束形成算法是目前使用最为广泛的自适应算法。该方法依据保持期望方向增益不变，并使阵列输出能量达到最小的原则，通过动态计算聚焦延时后的回波信号加权矢量，再将该矢量与回波信号相乘完成自适应加权，提高了图像对比度和分辨率。但该算法的缺点是在算法涉及矩阵运算，复杂度高，严重影响了成像实时性，并且算法稳健性不如传统DAS算法；并且MV算法最初是针对窄带、非相关信号设计的，而超声信号具有宽带和强相关特性，并不满足MV算法的适用条件。因此，最小方差算法在成像分辨率、对比度、鲁棒性和成像效率都还有很大的提升空间。
为了提高MV算法的性能，对角加载方法和空间平滑方法分别用于提高算法的鲁棒性和降低超声信号的强相关性。在MV算法的基础上，提出了基于特征空间的最小方差算法(ESBMV，Eigenspace-based Minimum Variance)，该方法对样本估计协方差矩阵进行特征值分解，并根据特征值的大小将信号子空间与噪声子空间分离。ESBMV虽然进一步提高了分辨率和对比度，但由于需要计算大量的特征值和特征向量，因此具有很高的复杂度和很低的运算效率。
综上所述，急需发明一种适用于超声回波信号特点，从本质上提高图像分辨率和对比度，并且保持算法运行效率和稳健性的自适应波束形成算法，以全面整体提高超声成像质量。
发明内容
有鉴于此，本发明的目的在于提供一种基于频域分段的高分辨率最小方差超声成像方法，该方法能显著提高算法的成像分辨率和对比度，并同时提高算法效率和波束形成稳健性，有效克服了宽带、强相关的超声回波信号不满足传统最小方差自适应波束形成算法窄带、非相关的应用条件的矛盾问题，从而全面提高了超声成像的质量。
为达到上述目的，本发明提供如下技术方案：
一种基于频域分段的高分辨率最小方差超声成像方法，该方法包括以下步骤：
S1：对超声阵元接收的回波信号进行放大、AD转换和延时处理，以获得超声回波数据；得到延时处理之后的信号x(τ)＝[x1(τ),x2(τ),...xN(τ)]，x1(τ)...xN(τ)分别表示各阵元接收的回波信号，N表示超声阵元数，τ表示为对应深度的采样时刻；
S2：根据STFT中的自适应窗函数的最大集中度测量准则，选取适合超声回波信号的最优窗函数；
S3：依据S2所选窗函数，对各阵元的超声回波信号进行STFT频域分段处理，获得等间距的窄带子信号；
S4：利用STFT的共轭对称性，前一半窄带子信号经过共轭对称处理生成另一半窄带信号；
S5：利用窗函数无信号重叠特性，对同一阵元的窄带子信号按窗函数滑动顺序进行重构，生成各阵元新的频域信号；
S6：将接收阵列依次划分为一个具有重叠阵元的子阵，然后对相应接收子阵的频域信号进行前后向平滑和对角加载处理，以获得频域的样本协方差矩阵；
S7：根据线性约束最小方差原则，计算得出频域分段最小方差波束形成权值；
S8：利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号。
进一步，在步骤S2中，根据STFT中的自适应窗函数的最大集中度测量准则，选取适合超声回波信号的最优窗函数，具体包括以下步骤：
S21：采用基于超声回波信号x(τ)的自适应窗函数的STFT结果S(k,ω)为：
Figure GDA0002219367240000021
其中，ω＝0,1,...,W-1，W是窄带子信号的长度，zj(k,ω)是需要求取的自适应窗函数，j(k,ω)是用于确定窗函数的时刻k和频率ω的索引函数，i是虚数变量；
S22：根据STFT的最大集中度测量准则选择适用于超声回波信号的最优窗函数，最大集中度测量准则表示为：
Figure GDA0002219367240000031
其中，jMC(k,ω)应用最大集中度测量准则确定窗函数时刻k和频率ω的索引函数，argmax是对集合范围内求最大值函数，Θω是包含矩形窗、三角窗、布尔曼窗、汉明窗和汉宁窗的窗函数集合；Cp(k,ω)是最大集中度测量值；Sp(τ,q)是超声信号x(τ)使用自适应窗函数p的STFT结果；q表示对应子频带的采样频率，D(k,ω)独立于时刻变量k，是频率变量ω的低通加权函数：
Figure GDA0002219367240000032
其中，zp(τ-k)是选择自Θω窗函数集合的自适应窗函数。
进一步，在步骤S3中，依据S2所选窗函数，对各阵元的超声回波信号进行STFT频域分段处理，获得等间距的窄带子信号，具体包括以下步骤：
S31：超声回波信号x(τ)通过STFT实现的窄带频域分段如下式所示：
Figure GDA0002219367240000033
其中，z(τ)是通过步骤S2选取的窗长和短时傅里叶变换点数均为64，无信号重叠的汉宁窗；
S32：通过STFT，将每个传感器阵元的超声回波信号转换为若干个独立的等间隔窄带子信号，第n个阵元上的第m个窄带的子信号Sn(m,ω)表达式为：
Sn(m,ω)＝[Sn(W·(m-1)+1),...,Sn(W·(m-1)+W-1),Sn(W·m)]
其中，m＝1,2,...,M，M是窗函数的长度，等同于分段窄带数；ω是窄带子信号频率变量，ω＝1,2,...,W，W是窄带子信号的长度。
进一步，在步骤S4中，利用STFT的共轭对称性，前一半窄带子信号经过共轭对称处理生成另一半窄带信号，经过共轭对称处理后的窄带子信号Sn(m,ω)的维度简化为：
Figure GDA0002219367240000041
其中，fix(·)是取整函数，length(·)是求取信号长度的函数。
进一步，在步骤S5中，利用窗函数无信号重叠特性，对同一阵元的窄带子信号按窗函数滑动顺序进行重构，生成各阵元新的频域信号Xn(Ω)，其表达式为：
Xn(Ω)＝[Sn(1,ω),...,Sn(M-1,ω),Sn(M,ω)]
其中，Ω＝1,2,...MW是重构频域子带信号的长度，该信号在频域的波束形成y(Ω)为：
Figure GDA0002219367240000042
其中，w(Ω)＝[w1(Ω),w2(Ω),...,wN(Ω)]是需要计算的频域自适应加权向量，(·)H表示共轭转置运算，重构频域信号X(Ω)＝[X1(Ω),X2(Ω),...,XN(Ω)]。
进一步，在步骤S6中，将接收阵列依次划分为一个具有重叠阵元的子阵，然后对相应接收子阵的频域信号进行前后向平滑和对角加载处理，以获得频域前后向样本协方差矩阵，具体包括以下步骤：
S61：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵Rl(Ω)，然后根据以下公式计算频域前向平滑估计协方差矩阵R(Ω)：
Figure GDA0002219367240000043
其中Xl(Ω)＝[Xl(Ω),Xl+1(Ω),...,Xl+L-1(Ω)]表示第l个子阵的频域前向平滑向量，Xl(Ω)...Xl+L-1(Ω)是第l个子阵中的各重叠阵元的平滑向量，l＝1,2,...,N-L+1，Xl(Ω)H为Xl(Ω)的共轭转置；
S62：通过以下计算公式对频域前向估计协方差矩阵R(Ω)进行对角加载处理，得到对角加载后的协方差矩阵
Figure GDA0002219367240000044
Figure GDA0002219367240000045
其中，ε＝trace(R(Ω))·δ，trace(R(Ω))为信号的等效功率，trace(·)是求矩阵迹的函数，δ为空间噪声与信号功率之比，
Figure GDA0002219367240000051
I为单位矩阵；
S63：通过以下计算公式由频域前向估计协方差矩阵求得后向估计协方差矩阵，并进行求和平均得到频域前后向估计协方差矩阵
Figure GDA0002219367240000052
Figure GDA0002219367240000053
其中，J是I的反对角矩阵，
Figure GDA0002219367240000054
是
Figure GDA0002219367240000055
的共轭转置。
进一步，在步骤S7中，根据线性约束最小方差原则，计算得出频域分段最小方差波束形成权值：
Figure GDA0002219367240000056
其中，wSTFTMV为频域自适应波束形成权值，a为全1的方向向量，aH是a的共轭转置，
Figure GDA0002219367240000057
是协方差矩阵
Figure GDA0002219367240000058
的逆矩阵。
进一步，在步骤S8中，利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号，具体包括以下步骤：
S81：在步骤S7获取频域分自适应波束形成权值的基础上，基于频域分段的的最小方差算法在频域的波束形成输出ySTFTMV(Ω)为：
Figure GDA0002219367240000059
其中，
Figure GDA00022193672400000510
是频域自适应波束形成权值wSTFTMV的共轭转置；
S82：通过快速傅里叶逆变换计算得到最终时域自适应波束形成信号ySTFTMV(k)为：
Figure GDA00022193672400000511
其中，IFFT(·)是快速傅里叶逆变换函数。
本发明的有益效果在于：
本发明提供了一种基于频域分段的高分辨率最小方差超声成像方法，该方法首先根据超声回波数据特点和最大集中度原则选取STFT中的频域分段最优窗函数，并将回波信号转换为窄带子信号以满足MV算法的窄带要求，显著提高了成像分辨率；其次利用STFT的共轭对称性，可减少一半数据处理量以提高算法效率；同时应用前后向空间平滑方法降低信号强相关性并提高算法稳健性；因此，本发明的方法显著改善了图像分辨率和对比度，同时提高了成像效率以及算法稳健性，有效克服了宽带、强相关的超声回波信号不满足传统最小方差自适应波束形成算法窄带、非相关的应用条件的矛盾，不能显著提高图像质量和效率等问题。
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：
图1为本发明所述方法的流程图；
图2为前后向空间平滑算法示意图；
图3为4种算法点目标成像结果；
图4为45mm、60mm和75mm处4种算法横向分辨率曲线图；
图5为4种算法不同深度处横向分辨率曲线；
图6为4种算法吸声斑目标成像结果；
图7为4种算法geabr_0数据成像结果；
图8为geabr_0实验75mm处散射点横向截面图。
具体实施方式
下面将结合附图，对本发明的优选实施例进行详细的描述。
图1为本发明所述方法的流程图，图2为前后向空间平滑算法示意图，如图所示，本发明提供一种基于时频分段的高分辨率最小方差超声成像方法，具体包括以下步骤：
步骤S1：对超声阵元接收的回波信号进行放大、AD转换和延时处理，以获得超声回波数据；得到延时处理之后的信号x(τ)＝[x1(τ),x2(τ),...xN(τ)]，x1(τ)...xN(τ)分别表示各阵元接收的回波信号，N表示超声阵元数，τ表示为对应深度的采样时刻；
步骤S2：根据STFT中的自适应窗函数的最大集中度测量准则，选取适合超声回波信号的最优窗函数，具体包括以下步骤：
S21：采用基于超声回波信号x(τ)的自适应窗函数的STFT结果S(k,ω)为：
Figure GDA0002219367240000061
其中，ω＝0,1,...,W-1，W是窄带子信号的长度，zj(k,ω)是需要求取的自适应窗函数，j(k,ω)是用于确定窗函数的时刻k和频率ω的索引函数，i是虚数变量。
S22：根据STFT的最大集中度测量准则选择适用于超声回波信号的最优窗函数，最大集中度测量准则表示为：
Figure GDA0002219367240000071
其中，jMC(k,ω)应用最大集中度测量准则确定窗函数时刻k和频率ω的索引函数，argmax是对集合范围内求最大值函数，Θω是包含矩形窗、三角窗、布尔曼窗、汉明窗和汉宁窗的窗函数集合；Cp(k,ω)是最大集中度测量值；Sp(τ,q)是超声信号x(τ)使用自适应窗函数p的STFT结果；q表示对应子频带的采样频率，D(k,ω)严格而言独立于时刻变量k，是频率变量ω的低通加权函数：
Figure GDA0002219367240000072
其中，zp(τ-k)是选择自Θω窗函数集合的自适应窗函数。
步骤S3：依据S2所选窗函数，对各阵元的超声回波信号进行STFT频域分段处理，获得等间距的窄带子信号，具体包括以下步骤：
S31：超声回波信号x(τ)通过STFT实现的窄带频域分段如下式所示：
Figure GDA0002219367240000073
其中，z(τ)是通过步骤2选取的窗长和短时傅里叶变换点数均为64，无信号重叠的汉宁窗。
S32：通过STFT，将每个传感器阵元的超声回波信号转换为若干个独立的等间隔窄带子信号，第n个阵元上的第个m窄带的子信号Sn(m,ω)表达式为：
Sn(m,ω)＝[Sn(W·(m-1)+1),...,Sn(W·(m-1)+W-1),Sn(W·m)]
其中，m＝1,2,...,M，M是窗函数的长度，等同于分段窄带数；ω是窄带子信号频率变量，ω＝1,2,...,W，W是窄带子信号的长度。
步骤S4：利用STFT的共轭对称性，前一半窄带子信号可经过共轭对称处理生成另一半窄带信号；经过共轭对称处理后的窄带子信号Sn(m,ω)的维度简化为：
Figure GDA0002219367240000081
其中，fix(·)是取整函数，length(·)是求取信号长度的函数。
步骤S5：利用窗函数无信号重叠特性，对同一阵元的窄带子信号按窗函数滑动顺序进行重构，生成各阵元新的频域信号Xn(Ω)，其表达式为：
Xn(Ω)＝[Sn(1,ω),...,Sn(M-1,ω),Sn(M,ω)]
其中，Ω＝1,2,...MW是重构频域子带信号的长度，该信号在频域的波束形成y(Ω)为：
Figure GDA0002219367240000082
其中，w(Ω)＝[w1(Ω),w2(Ω),...,wN(Ω)]是需要计算的频域自适应加权向量，(·)H表示共轭转置运算，重构频域信号X(Ω)＝[X1(Ω),X2(Ω),...,XN(Ω)]。
步骤S6：将接收阵列依次划分为一个具有重叠阵元的子阵，然后对相应接收子阵的频域信号进行前后向平滑和对角加载处理，以获得频域的样本协方差矩阵，具体包括以下步骤：
S61：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵Rl(Ω)，然后根据以下公式计算频域前向平滑估计协方差矩阵R(Ω)：
Figure GDA0002219367240000083
其中Xl(Ω)＝[Xl(Ω),Xl+1(Ω),...,Xl+L-1(Ω)]表示第l个子阵的频域前向平滑向量，Xl(Ω)...Xl+L-1(Ω)是第l个子阵中的各重叠阵元的平滑向量，l＝1,2,...,N-L+1，Xl(Ω)H为Xl(Ω)的共轭转置。
S62：通过以下计算公式对频域前向估计协方差矩阵R(Ω)进行对角加载处理，得到对角加载后的协方差矩阵
Figure GDA0002219367240000084
Figure GDA0002219367240000085
其中，ε＝trace(R(Ω))·δ，trace(R(Ω))为信号的等效功率，trace(·)是求矩阵迹的函数，δ为空间噪声与信号功率之比，
Figure GDA0002219367240000091
I为单位矩阵。
S63：通过以下计算公式由频域前向估计协方差矩阵求得后向估计协方差矩阵，并进行求和平均得到频域前后向估计协方差矩阵
Figure GDA0002219367240000092
Figure GDA0002219367240000093
其中，J是I的反对角矩阵，
Figure GDA0002219367240000094
是
Figure GDA0002219367240000095
的共轭转置。
步骤S7：根据线性约束最小方差原则，计算得出频域分段最小方差波束形成权值：
Figure GDA0002219367240000096
其中，wSTFTMV为频域自适应波束形成权值，a为全1的方向向量，aH是a的共轭转置，
Figure GDA0002219367240000097
是协方差矩阵
Figure GDA0002219367240000098
的逆矩阵。
步骤S8：利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号，具体包括以下步骤：
S81：在步骤7获取频域分自适应波束形成权值的基础上，基于频域分段的的最小方差算法在频域的波束形成输出ySTFTMV(Ω)为：
Figure GDA0002219367240000099
其中，
Figure GDA00022193672400000910
是频域自适应波束形成权值wSTFTMV的共轭转置。
S82：通过快速傅里叶逆变换计算得到最终时域自适应波束形成信号ySTFTMV(k)为：
Figure GDA00022193672400000911
其中，IFFT(·)是快速傅里叶逆变换函数。
在本实施例中，为验证所提算法的有效性，利用Field II对超声成像中常用的点散射目标和吸声斑目标进行成像并利用实际实验数据进行成像对比实验。Field II是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。在点目标仿真实验中，设置两列横向间隔为2mm，纵向间隔为5mm的18个点目标，深度分布在40mm～80mm之间，采用发射定点聚焦和接收动态聚焦方式，发射焦点固定在60mm处，并在接收回波中加入一定噪声，设置图像的成像动态范围为50dB。同时，设一中心在45mm，半径为3mm的圆形区域吸声斑，外部随机分布着100000个散射点，在接收回波中加入一定噪声，并设定成像动态范围为50dB。实验采用密歇根大学生物医学超声实验室提供的完备数据集gearb_0，所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为60dB。对上述三个实验目标采用延时叠加算法(DAS)，最小方差算法(MV)，特征空间最小方差算法(ESBMV)，基于时频分段的最小方差算法(STFTMV)进行对比成像实验。
图3给出了4种算法点目标成像结果，从图3中可以看出DAS算法成像质量最差，分辨率最低，相比于其他4种算法横向伪影最多，难以区分相邻目标点。MV算法较DAS算法旁瓣有所降低，在焦点处散射点已基本能够区分，但在其他深度处横向伪影仍较多，分辨率有待提高。ESBMV算法在整个深度范围内可以明显分辨出相邻目标点，成像质量较MV有一定提高。STFTMV算法成像质量最优，分辨率和对比度最高，对噪声鲁棒性最好，对点目标的分辨能力最好。
图4给出了55mm焦点处4种算法横向分辨率曲线图，图5给出4种算法不同深度处横向分辨率曲线，其中(a)为-6dB点目标处分辨率，(b)为-20dB点目标处分辨率。从图4中可以看出，DAS算法成像分辨率最差，主瓣宽度最宽且旁瓣等级最高。MV算法较DAS算法成像有所提高，其主瓣宽度和旁瓣等级都有所改善。ESBMV算法与DAS和MV算法相比，主瓣宽度以及旁瓣等级已经改善明显，STFTMV算法主瓣最窄，其主瓣宽度较MV算法和ESBMV算法减小了50％左右，旁瓣等级最低，较MV算法和ESBMV算法分别降低10.78dB和6.21dB，图像对比度最高。从图5可以看出4种算法横向分辨率随着深度的增加呈降低的趋势，由于焦点在60mm处，因此，在焦点处分辨率有所改善，分辨率曲线出现拐点；由图5可以得出结论：在不同深度处，STFTMV算法点目标分辨率均优于DAS、MV和ESBMV算法。
图6给出4种算法吸声斑目标成像结果，表1给出4种算法对比度。从图6中可以看出，DAS算法相比于其他算法成像效果最差，噪声抑制能力最弱，吸声斑内部存在噪声干扰严重。MV算法和ESBMV算法对噪声的抑制较DAS有所改善。STFTMV算法噪声含量最小，算法旁瓣抑制能力最强。由表1可见，DAS算法对比度最低，由于其只进行简单的叠加成像，计算复杂度低，因而背景方差小，算法稳健性好。MV算法提高了中心暗斑平均功率，但其外部平均功率也同时提高，对比度较DAS算法上升约2dB。ESBMV以及算法中心暗斑及背景功率分别在MV基础上有所提高，对比度较MV算法上升约2dB。其中，STFTMV算法中心平均功率上升最多，对比度较DAS，MV，ESBMV算法分别提高了5.41dB、3.25dB、1.22dB，且背景区域方差最低，对比噪声比最高，算法稳健性最好。
表1 4种算法对比度
Figure GDA0002219367240000111
表2给出4种算法的点目标和吸声斑仿真结果的运行时间。由于点目标深度扫描范围为40～80mm，吸声斑深度扫描范围为40～55mm，因此点目标的仿真成像时间比吸声斑长。从表2中可以看出，传统DAS算法由于计算简单，复杂度低，所以其点目标和吸声斑的成像时间最短，所需时长是MV算法的10％左右。在自适应算法中，ESBMV算法由于获取大量特征值和特征向量，其点目标和吸声斑成像时间最长，是MV算法的2.5倍左右；STFTMV的成像时间最短，约占ESBMV的35％左右，且优于MV算法；因此所提STFTMV算法在自适应算中成像效率最高。
表2 4种算法成像时间
Figure GDA0002219367240000112
图7给出了4种算法geabr_0数据成像结果；图8给出geabr_0实验75mm处散射点横向截面图。从图7中可以看出，传统DAS算法成像效果最差，近场点目标受背景噪声干扰最为严重，采用自适应算法成像较DAS算法都好，其图像分辨率和对比度都有所改善，其中STFTMV算法分辨率最高，对比度改善明显。从图8可以看出，ESBMV与MV算法分辨率相当且都高于传统DAS算法，而ESBMV算法进一步降低了旁瓣等级，提高了对比度。STFTMV算法分辨率和对比度最高，其主瓣宽度最窄，旁瓣等级最低，因此所提基于频域分段的高分辨率最小方差算法成像效果最佳，与点目标和吸声斑仿真结论基本相符。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (6)

1.一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：该方法包括以下步骤：
S1：对超声阵元接收的回波信号进行放大、AD转换和延时处理，以获得超声回波数据；得到延时处理之后的信号x(τ)＝[x1(τ),x2(τ),...xN(τ)]，x1(τ)～xN(τ)分别表示各阵元接收的回波信号，N表示超声阵元数，τ表示为对应深度的采样时刻；
S2：根据短时傅里叶变换(Short-time Fourier Transform，STFT)中的自适应窗函数的最大集中度测量准则，选取适合超声回波信号的最优窗函数；根据STFT中的自适应窗函数的最大集中度测量准则，选取适合超声回波信号的最优窗函数，具体包括以下步骤：
S21：采用基于超声回波信号x(τ)的自适应窗函数的STFT结果S(k,ω)为：
Figure FDA0003210686070000011
其中，ω＝0,1,...,W-1，W是窄带子信号的长度，zj(k,ω)是需要求取的自适应窗函数，j(k,ω)是用于确定窗函数的时刻k和频率ω的索引函数，i是虚数变量；
S22：根据STFT的最大集中度测量准则选择适用于超声回波信号的最优窗函数，最大集中度测量准则表示为：
Figure FDA0003210686070000012
其中，jMC(k,ω)应用最大集中度测量准则确定窗函数时刻k和频率ω的索引函数，argmax是对集合范围内求最大值函数，Θω是包含矩形窗、三角窗、布尔曼窗、汉明窗和汉宁窗的窗函数集合；Cp(k,ω)是最大集中度测量值；Sp(τ,q)是超声信号x(τ)使用自适应窗函数p的STFT结果；q表示对应子频带的采样频率，D(k,ω)独立于时刻变量k，是频率变量ω的低通加权函数：
Figure FDA0003210686070000013
其中，zp(τ-k)是选择自Θω窗函数集合的自适应窗函数；
S3：依据S2所选窗函数，对各阵元的超声回波信号进行STFT频域分段处理，获得等间距的窄带子信号；具体包括以下步骤：
S31：超声回波信号x(τ)通过STFT实现的窄带频域分段如下式所示：
Figure FDA0003210686070000021
其中，z(τ)是通过步骤S2选取的窗长和短时傅里叶变换点数均为64，无信号重叠的汉宁窗；
S32：通过STFT，将每个传感器阵元的超声回波信号转换为若干个独立的等间隔窄带子信号，第n个阵元上的第m个窄带的子信号Sn(m,ω)表达式为：
Sn(m,ω)＝[Sn(W·(m-1)+1),...,Sn(W·(m-1)+W-1),Sn(W·m)]
其中，m＝1,2,...,M，M是窗函数的长度，等同于分段窄带数；ω是窄带子信号频率变量，ω＝1,2,...,W，W是窄带子信号的长度；
S4：利用STFT的共轭对称性，前一半窄带子信号经过共轭对称处理生成另一半窄带信号；
S5：利用窗函数无信号重叠特性，对同一阵元的窄带子信号按窗函数滑动顺序进行重构，生成各阵元新的频域信号；
S6：将接收阵列依次划分为一个具有重叠阵元的子阵，然后对相应接收子阵的频域信号进行前后向平滑和对角加载处理，以获得频域的样本协方差矩阵；
S7：根据线性约束最小方差原则，计算得出频域分段最小方差波束形成权值；
S8：利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号。
2.根据权利要求1所述的一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：在步骤S4中，利用STFT的共轭对称性，前一半窄带子信号经过共轭对称处理生成另一半窄带信号，经过共轭对称处理后的窄带子信号Sn(m,ω)的维度简化为：
Figure FDA0003210686070000022
其中，fix(·)是取整函数，length(·)是求取信号长度的函数。
3.根据权利要求2所述的一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：在步骤S5中，利用窗函数无信号重叠特性，对同一阵元的窄带子信号按窗函数滑动顺序进行重构，生成各阵元新的频域信号Xn(Ω)，其表达式为：
Xn(Ω)＝[Sn(1,ω),...,Sn(M-1,ω),Sn(M,ω)]
其中，Ω是重构频域子带信号的长度，Ω＝1,2,...MW，该信号在频域的波束形成y(Ω)为：
Figure FDA0003210686070000031
其中，wH(Ω)表示w(Ω)的共轭转置，w(Ω)是需要计算的频域自适应加权向量，w(Ω)＝[w1(Ω),w2(Ω),...,wN(Ω)]，重构频域信号X(Ω)＝[X1(Ω),X2(Ω),...,XN(Ω)]。
4.根据权利要求3所述的一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：在步骤S6中，将接收阵列依次划分为一个具有重叠阵元的子阵，然后对相应接收子阵的频域信号进行前后向平滑和对角加载处理，以获得频域前后向样本协方差矩阵，具体包括以下步骤：
S61：把N个阵元依次划分为阵元数目为L的子阵，并分别计算各个子阵的样本协方差矩阵Rl(Ω)，然后根据以下公式计算频域前向平滑估计协方差矩阵R(Ω)：
Figure FDA0003210686070000032
其中Xl(Ω)＝[Xl(Ω),Xl+1(Ω),...,Xl+L-1(Ω)]表示第l个子阵的频域前向平滑向量，Xl(Ω)～Xl+L-1(Ω)是第l个子阵中的各重叠阵元的平滑向量，l＝1,2,...,N-L+1，Xl(Ω)H为Xl(Ω)的共轭转置；
S62：通过以下计算公式对频域前向估计协方差矩阵R(Ω)进行对角加载处理，得到对角加载后的协方差矩阵
Figure FDA0003210686070000033
Figure FDA0003210686070000034
其中，ε＝trace(R(Ω))·δ，trace(R(Ω))为信号的等效功率，trace(·)是求矩阵迹的函数，δ为空间噪声与信号功率之比，
Figure FDA0003210686070000035
I为单位矩阵；
S63：通过以下计算公式由频域前向估计协方差矩阵求得后向估计协方差矩阵，并进行求和平均得到频域前后向估计协方差矩阵
Figure FDA0003210686070000036
Figure FDA0003210686070000041
其中，J是I的反对角矩阵，
Figure FDA0003210686070000042
是
Figure FDA0003210686070000043
的共轭转置。
5.根据权利要求4所述的一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：在步骤S7中，根据线性约束最小方差原则，计算得出频域分段最小方差波束形成权值：
Figure FDA0003210686070000044
其中，wSTFTMV为频域自适应波束形成权值，a为全1的方向向量，aH是a的共轭转置，
Figure FDA0003210686070000045
是协方差矩阵
Figure FDA0003210686070000046
的逆矩阵。
6.根据权利要求5所述的一种基于频域分段的高分辨率最小方差超声成像方法，其特征在于：在步骤S8中，利用快速傅里叶逆变换对频域分段最小方差波束形成权值进行处理，得出最终时域自适应波束形成信号，具体包括以下步骤：
S81：在步骤S7获取频域分自适应波束形成权值的基础上，基于频域分段的的最小方差算法在频域的波束形成输出ySTFTMV(Ω)为：
Figure FDA0003210686070000047
其中，
Figure FDA0003210686070000048
是频域自适应波束形成权值wSTFTMV的共轭转置；
S82：通过快速傅里叶逆变换计算得到最终时域自适应波束形成信号ySTFTMV(k)为：
Figure FDA0003210686070000049
其中，IFFT(·)是快速傅里叶逆变换函数。
CN201910757165.6A 2019-08-15 2019-08-15 一种基于频域分段的高分辨率最小方差超声成像方法 Active CN110501423B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910757165.6A CN110501423B (zh) 2019-08-15 2019-08-15 一种基于频域分段的高分辨率最小方差超声成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910757165.6A CN110501423B (zh) 2019-08-15 2019-08-15 一种基于频域分段的高分辨率最小方差超声成像方法

Publications (2)

Publication Number Publication Date
CN110501423A CN110501423A (zh) 2019-11-26
CN110501423B true CN110501423B (zh) 2021-10-08

Family

ID=68588191

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910757165.6A Active CN110501423B (zh) 2019-08-15 2019-08-15 一种基于频域分段的高分辨率最小方差超声成像方法

Country Status (1)

Country Link
CN (1) CN110501423B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111208213A (zh) * 2020-02-25 2020-05-29 重庆大学 融合交替乘子迭代的谱寻求子带最小方差超声成像算法
CN112035788A (zh) * 2020-06-12 2020-12-04 中国航天科工集团第二研究院 一种提高超声系统成像质量的方法
CN113296064A (zh) * 2021-04-13 2021-08-24 武汉卓目科技有限公司 一种基于Frank码的SAR雷达接收通道时延校准方法及系统
CN113218560B (zh) * 2021-04-19 2022-05-17 中国长江电力股份有限公司 一种超声螺栓预紧力实时估计方法
CN113422630B (zh) * 2021-06-17 2023-02-07 长安大学 一种自适应聚焦宽带波束形成方法及系统
CN113569192B (zh) * 2021-08-05 2024-03-12 阳光学院 一种多相位分级的嵌套阵列天线波束合成方法
CN113647977B (zh) * 2021-08-18 2023-10-10 重庆大学 一种基于切比雪夫多项式的复合窗变迹超声波束形成方法
CN114268381A (zh) * 2021-12-15 2022-04-01 重庆大学 融合维纳后置滤波器的子频带分割最小方差超声成像方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000201352A (ja) * 1999-01-08 2000-07-18 Matsushita Electric Ind Co Ltd 画像圧縮符号化装置
CA2727415A1 (en) * 2008-06-10 2009-12-17 Uti Limited Partnership Signal processing with fast s-transforms
US9739861B2 (en) * 2012-04-19 2017-08-22 University Of Southern California Method for reduced field of view MRI in an inhomogeneous field with rapid outer volume suppression
US11116474B2 (en) * 2013-07-23 2021-09-14 Regents Of The University Of Minnesota Ultrasound image formation and/or reconstruction using multiple frequency waveforms
CN104361894A (zh) * 2014-11-27 2015-02-18 湖南省计量检测研究院 一种基于输出的客观语音质量评估的方法
CN106510761B (zh) * 2016-12-12 2019-06-21 重庆大学 一种信噪比后滤波与特征空间融合的最小方差超声成像方法
CN108761530B (zh) * 2018-05-22 2019-10-11 闽南师范大学 一种地震信号谱分解方法
CN108880586B (zh) * 2018-06-28 2019-11-08 中国人民解放军战略支援部队信息工程大学 一种宽带弱信号增强方法与装置

Also Published As

Publication number Publication date
CN110501423A (zh) 2019-11-26

Similar Documents

Publication Publication Date Title
CN110501423B (zh) 一种基于频域分段的高分辨率最小方差超声成像方法
Synnevag et al. Benefits of minimum-variance beamforming in medical ultrasound imaging
Nilsen et al. Beamspace adaptive beamforming for ultrasound imaging
Asl et al. Eigenspace-based minimum variance beamforming applied to medical ultrasound imaging
CN109490850B (zh) 主瓣干扰下宽带阵列自适应波束形成方法
CN111208213A (zh) 融合交替乘子迭代的谱寻求子带最小方差超声成像算法
CN108836389B (zh) 平面波相关点相干自适应波束合成成像方法
Diamantis et al. Experimental performance assessment of the sub-band minimum variance beamformer for ultrasound imaging
Holfort et al. P2b-12 minimum variance beamforming for high frame-rate ultrasound imaging
Agarwal et al. Improving spatial resolution using incoherent subtraction of receive beams having different apodizations
WO2005008280A1 (en) Corrections for wavefront aberrations in ultrasound imaging
CN108761466B (zh) 波束域广义旁瓣相消超声成像方法
CN105760892A (zh) 一种改进的最小方差超声成像方法
CN111856474B (zh) 一种基于子阵的空时域条件相干系数超声成像方法
CN102764139B (zh) 基于特征空间分析和区域判别的医学超声波束形成方法
Kim et al. Efficient array beam forming by spatial filtering for ultrasound B-mode imaging
Salari et al. User parameter-free minimum variance beamformer in medical ultrasound imaging
CN110501711A (zh) 一种基于乘幂法的低复杂度最小方差超声成像方法
Kozai et al. Optimization of frequency and plane-wave compounding by minimum variance beamforming
CN109187771B (zh) 一种融合特征值分解的低复杂度最小方差超声成像方法
CN111278363B (zh) 超声成像设备、系统及其超声造影成像的图像增强方法
CN114415187A (zh) 基于子空间斜投影的强稳健性最小方差超声波束形成方法
Wang et al. Generalized sidelobe canceller beamforming method for ultrasound imaging
Izadi et al. Weighted Capon beamformer combined with coded excitation in ultrasound imaging
Okumura et al. Stabilization techniques for high resolution ultrasound imaging using beamspace Capon method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20210705

Address after: 400044 No. 174 Shapingba street, Shapingba District, Chongqing

Applicant after: Chongqing University

Applicant after: CHONGQING College OF ELECTRONIC ENGINEERING

Applicant after: Chongqing mostag Energy Management Co.,Ltd.

Address before: 400044 No. 174 Shapingba street, Shapingba District, Chongqing

Applicant before: Chongqing University

Applicant before: CHONGQING College OF ELECTRONIC ENGINEERING

GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20211210

Address after: 400044 No. 174 Shapingba street, Shapingba District, Chongqing

Patentee after: Chongqing University

Patentee after: CHONGQING College OF ELECTRONIC ENGINEERING

Patentee after: Chongqing mostag Energy Management Co.,Ltd.

Patentee after: Chongqing Caojie shipping power development Co.,Ltd.

Address before: 400044 No. 174 Shapingba street, Shapingba District, Chongqing

Patentee before: Chongqing University

Patentee before: CHONGQING College OF ELECTRONIC ENGINEERING

Patentee before: Chongqing mostag Energy Management Co.,Ltd.