CN116934891A - 一种基于多重幅度加权互相关的自适应滤波超声成像方法 - Google Patents
一种基于多重幅度加权互相关的自适应滤波超声成像方法 Download PDFInfo
- Publication number
- CN116934891A CN116934891A CN202310891567.1A CN202310891567A CN116934891A CN 116934891 A CN116934891 A CN 116934891A CN 202310891567 A CN202310891567 A CN 202310891567A CN 116934891 A CN116934891 A CN 116934891A
- Authority
- CN
- China
- Prior art keywords
- adaptive
- signal
- signals
- filtering
- weighting
- 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.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 68
- 238000001914 filtration Methods 0.000 title claims abstract description 41
- 230000003044 adaptive effect Effects 0.000 claims abstract description 56
- 239000011159 matrix material Substances 0.000 claims abstract description 39
- 238000005070 sampling Methods 0.000 claims abstract description 20
- 238000000034 method Methods 0.000 claims abstract description 19
- 238000012545 processing Methods 0.000 claims abstract description 12
- 230000000295 complement effect Effects 0.000 claims abstract description 4
- 238000004422 calculation algorithm Methods 0.000 claims description 55
- 238000012285 ultrasound imaging Methods 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 230000000903 blocking effect Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000001629 suppression Effects 0.000 description 7
- 238000010521 absorption reaction Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 230000006872 improvement Effects 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000002474 experimental method Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
- A61B8/5269—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving detection or reduction of artifacts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
- Y02D30/00—Reducing energy consumption in communication networks
- Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Animal Behavior & Ethology (AREA)
- Veterinary Medicine (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Pathology (AREA)
- General Physics & Mathematics (AREA)
- Biophysics (AREA)
- Data Mining & Analysis (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Ultra Sonic Daignosis Equipment (AREA)
Abstract
本发明涉及一种基于多重幅度加权互相关的自适应滤波超声成像方法，属于超声成像技术领域；该方法包括：对采样信号进行延时处理以获得超声回波数据，并将希尔伯特变换应用于回波数据；构造多对互补的接收幅度加权函数，对各阵元延时信号进行加权得到两组波束形成信号并计算其两两之间的信号相似程度；对得到的系数阈值化后取平均值，然后对其组成的加权矩阵进行中值滤波；利用中值滤波后的系数自适应确定空间平滑子阵长度进一步估计协方差矩阵，然后利用高分辨率的自适应滤波结构计算权值矢量，并对回波信号进行加权；最后将中值滤波后的系数对改进的自适应滤波结构输出进行加权，得到自适应波束形成的扫描线信号，从而提升超声成像质量。
Description
技术领域
本发明属于超声成像技术领域，涉及一种基于多重幅度加权互相关的自适应滤波超声成像方法。
背景技术
超声成像以其安全性、实时性和低成本等优势在医学临床诊断中被广泛使用。数字波束形成在整个超声成像系统中处于核心位置。超声成像中应用最为广泛的，也是最简单的波束形成技术即延时叠加算法(Delay And Sum，DAS)，它是根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算，然后对延时后的数据对齐叠加。传统DAS算法复杂度低，成像速度快，但由于其采用固定窗函数加权导致主瓣宽度增加，分辨率较低。
最小方差无失真波束形成器(Minimum Variance,MV)能够有效提高成像图像分辨率，通过干扰信号和期望信号的方向性差异对干扰及噪声信号进行滤除。然而其鲁棒性差，不利于实际应用。广义旁瓣相消器(Generalized Sidelobe Canceller,GSC)作为最小方差无失真波束形成器的一种等效结构其对成像的改善效果与其相差不大，能够显著改善横向分辨率，但旁瓣抑制能力仍然较弱，对比度方面未取得相应的改进。除此之外，由于自适应波束形成方法包含大量矩阵运算，自适应波束形成的鲁棒性不及传统DAS算法。
综上所述，目前亟需一种既能保持高分辨率性能，又能提高成像对比度且具有强鲁棒性的波束形成算法。
发明内容
有鉴于此，本发明的目的在于提供一种基于多重幅度加权互相关的高分辨率自适应滤波超声成像算法，能够在保持传统广义旁瓣相消算法良好分辨率性能的情况下，有效提高算法鲁棒性和成像对比度，从而提高算法综合成像质量。
为达到上述目的，本发明提供如下技术方案：
一种基于多重幅度加权互相关的自适应滤波超声成像方法，具体包括以下步骤：
S1：首先对超声阵元接收的回波信号进行放大，AD转换和延时聚焦处理以获得超声回波数据，然后将希尔伯特变换应用于回波数据，得到延时聚焦处理和希尔伯特变换之后的信号x(k)；
S2：构造多对互补的接收幅度加权函数，对各阵元延时信号进行加权得到两组波束形成信号，然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数；
S3：对得到的归一化互相关系数阈值化后取平均值，对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波；
S4：利用中值滤波后的系数确定自适应的空间平滑子阵长度，利用自适应子阵长度估计协方差矩阵，并进行对角加载处理；
S5：利用高分辨率自适应滤波器的结构计算出最优权值矢量，利用权值对超声阵列回波信号进行加权处理；
S6：将中值滤波后的系数作为加权因子乘上改进的高分辨率自适应滤波结构的输出，得到自适应波束形成的扫描线信号，采用优化后的扫描线信号进行最终成像。
进一步，步骤S2具体包括以下步骤：
S21：对于有N个等间距阵元的传感器阵列，通过交替提取M个阵元来构建两个接收幅度加权函数Apod1n和Apod2n，其中n＝1,2...N，j＝1,2...M；
S22：将这两个接收幅度加权函数分别进一步划分为M个子幅度加权函数和其中j＝1,2...M；
S23：使用所得子幅度加权函数对各阵元延时信号进行处理，分别获得M个具有相位多样性的波束形成信号BF1j(k,h)和BF2j(k,h)，
其中，xn(k)为第n个阵元对应的延时回波信号，k为采样时刻。BF1j(k)和BF2j(k)表示一列扫描线信号的第k个采样时刻下的自适应波束信号，BF1j(k,h)和BF2j(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号。
然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数ρj(k,h)，表达式为：
其中，k表示成像区域中第h列的第k个采样时刻，A表示第k个采样时刻分别向上或向下取的样本长度，计算每个像素点时使用的互相关段总长度为2A+1个样本。
进一步，步骤S3具体包括以下步骤：
S31：对得到的M个归一化互相关系数阈值化，并取平均值ρ(k,h)：
其中，δ为阈值，简记为/>
S32：对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波，对矩阵中的任意一点(k,h)将以该点为中心设计中值滤波器窗口ωkh，然后对窗口内包含数值按从小到大的顺序进行排列，取排列中值作为(k,h)点的滤波结果，窗口大小为3λ×10，其中λ为波长。
进一步，步骤S4具体包括以下步骤：
S41：采用中值滤波后的系数自适应判别子阵长度L，定义自适应子阵长度/>为：
其中，N是阵元个数，round(·)是取整函数；整数的取值范围为1到N2；当时，/>为1；
S42：用自适应子阵长度代替固定子阵长度L，处理后的样本协方差矩阵R'(k)的计算公式为：
其中，(·)H表示共轭转置，表示第l个平滑子阵接收到的施加延时后的回波信号，/>
S43：通过以下计算公式对样本协方差矩阵R'(k)进行对角加载，得到对角加载后的样本协方差矩阵R'DL(k)：
R'DL(k)＝R'(k)+ε'·I'
其中，ε'＝Δ·trace(R'(k))，Δ为噪声和信号的比值，trace(·)表示矩阵的迹，I'表示单位阵。
进一步，步骤S5具体包括以下步骤：
S51：通过以下公式计算高分辨率自适应滤波器结构中自适应权矢量w'a：
w'a＝(B'HR'DLB')-1B'HR'DLw'q
其中，非自适应权矢量w'q是的固定权向量，阻塞矩阵B'是/>维矩阵，且满足B'Hw'q＝0，(·)-1表示矩阵求逆运算；
S52：得到改进的高分辨率自适应滤波结构的最优权值矢量w'为：
w'＝w'q-B'w'a
S53：利用改进的最优权值矢量对超声阵列回波信号进行加权求和得到自适应波束形成信号：
其中，y(k)表示计算得到一列扫描线信号的第k个采样时刻下的自适应波束信号，w'H表示w'的共轭转置，w'表示改进的高分辨率自适应滤波算法的自适应加权矢量，表示第l个子阵在第k个采样时刻的回波向量；w'(k)表示第k个采样时刻下的w'。
进一步，步骤S6中，得到自适应波束形成的扫描线信号y'(k,h)，表达式为：
其中，y(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号。
本发明的有益效果在于：本发明相比于现有的波束形成算法，能够在保证其高分辨率性能的同时，提高算法鲁棒性和成像对比度性能。本发明可以显著抑制吸声斑中存在的伪影，能够克服图像高分辨率、高对比度和强鲁棒性之间的权衡困难。
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：
图1为本发明所述超声成像方法的流程图；
图2为基于多重幅度加权互相关高分辨率自适应滤波器整体结构框架图；
图3为4种算法点目标成像结果；
图4为点目标成像60mm处4种算法横向分辨率曲线图；
图5为4种算法吸声斑目标成像结果；
图6为4种算法加入10dB/0dB/-10dB高斯白噪声的吸声斑成像结果，其中，(a)对应10dB；(b)对应0dB；(c)对应-10dB；
图7为4种算法geabr_0数据成像结果。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
请参阅图1～图7，图1为本发明方法流程图，如图1所示，本发明提供一种基于多重幅度加权互相关的高分辨率自适应滤波超声成像算法，具体包括以下步骤：
步骤S1：首先对超声阵元接收的回波信号进行放大，AD转换和延时聚焦处理以获得超声回波数据，然后将希尔伯特变换应用于回波数据，得到延时聚焦处理和希尔伯特变换之后的信号x(k)。
步骤S2：构造多对互补的接收幅度加权函数，对各阵元延时信号进行加权得到两组波束形成信号，然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数，具体包括以下步骤：
S21：对于有N个等间距阵元的传感器阵列，通过交替提取M个阵元来构建两个接收幅度加权函数Apod1n和Apod2n，其中n＝1,2...N，j＝1,2...M；
S22：将这两个接收幅度加权函数分别进一步划分为M个子幅度加权函数和其中j＝1,2...M；
S23：使用所得子幅度加权函数对各阵元延时信号进行处理，分别获得M个具有相位多样性的波束形成信号BF1j(k,h)和BF2j(k,h)，
其中，xn(k)为第n个阵元对应的延时回波信号，k为采样时刻。BF1j(k)和BF2j(k)表示一列扫描线信号的第k个采样时刻下的自适应波束信号，BF1j(k,h)和BF2j(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号。
然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数ρj(k,h)，表达式为：
其中，k表示成像区域中第h列的第k个采样时刻，A表示第k个采样时刻分别向上或向下取的样本长度，计算每个像素点时使用的互相关段总长度为2A+1个样本。
步骤S3：对得到的归一化互相关系数阈值化后取平均值，对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波，具体包括以下步骤：
S31：对得到的M个归一化互相关系数阈值化，并取平均值
其中，δ为阈值，简记为/>
S32：对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波，对矩阵中的任意一点(k,h)将以该点为中心设计中值滤波器窗口ωkh，然后对窗口内包含数值按从小到大的顺序进行排列，取排列中值作为(k,h)点的滤波结果，窗口大小为3λ×10，其中λ为波长。
步骤S4：利用中值滤波后的系数确定自适应的空间平滑子阵长度，利用自适应子阵长度估计协方差矩阵，并进行对角加载处理，具体包括以下步骤：
S41：采用中值滤波后的系数自适应判别子阵长度L，定义自适应子阵长度/>为：
其中，N是阵元个数，round(·)是取整函数；整数的取值范围为1到N2；当时，/>为1；
S42：用自适应子阵长度代替固定子阵长度L，处理后的样本协方差矩阵R'(k)的计算公式为：
其中，(·)H表示共轭转置，表示第l个平滑子阵接收到的施加延时后的回波信号，/>
S43：通过以下计算公式对样本协方差矩阵R'(k)进行对角加载，得到对角加载后的样本协方差矩阵R'DL(k)：
R'DL(k)＝R'(k)+ε'·I'
其中，ε'＝Δ·trace(R'(k))，Δ为噪声和信号的比值，trace(·)表示矩阵的迹，I'表示单位阵。
步骤S5：利用高分辨率自适应滤波器的结构计算出最优权值矢量，利用权值对超声阵列回波信号进行加权处理，具体包括以下步骤：
S51：通过以下公式计算高分辨率自适应滤波器结构中自适应权矢量w'a：
w'a＝(B'HR'DLB')-1B'HR'DLw'q
其中，非自适应权矢量w'q是的固定权向量，阻塞矩阵B'是/>维矩阵，且满足B'Hw'q＝0，(·)-1表示矩阵求逆运算；
S52：得到改进的高分辨率自适应滤波结构的最优权值矢量w'为：
w'＝w'q-B'w'a
S53：利用改进的最优权值矢量对超声阵列回波信号进行加权求和得到自适应波束形成信号：
其中，y(k)表示计算得到一列扫描线信号的第k个采样时刻下的自适应波束信号，w'H表示w'的共轭转置，w'表示改进的高分辨率自适应滤波算法的自适应加权矢量，表示第l个子阵在第k个采样时刻的回波向量；w'(k)表示第k个采样时刻下的w'。
步骤S6：将中值滤波后的系数作为加权因子乘上改进的高分辨率自适应滤波结构的输出，得到自适应波束形成的扫描线信号，采用优化后的扫描线信号进行最终成像：
其中，y(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号。
实验验证：
Field II是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证本发明算法的有效性，利用Field II对超声成像中常用的点散射目标和吸声斑目标进行成像，并利用实际实验数据geabr_0进行成像对比实验。在点目标仿真实验中，设置了纵向间距为5mm，横向间距为6mm的18个点目标，均匀分布在探测深度30mm～65mm、横向距离-9～9mm的成像范围内，并设置图像的成像动态范围为60dB，以此观察各个算法的横向分辨率强弱。在吸声斑仿真实验中，设一半径为4mm的圆形吸声斑，其中心位置在40mm深度处，纵向成像范围在30～50mm之间，横向成像距离为20mm，成像区域内随机分布200000个散射点，背景区域散射点的幅度是吸声斑内部区域的1000倍，设定成像动态范围为65dB。为了验证算法的鲁棒性，对于吸声斑仿真实验，在进行波束形成前对接收数据分别加入了SNR为10dB/0dB/-10dB的高斯白噪声以模拟实际的噪声情况，设定成像动态范围为65dB。仿真都采用发射定点聚焦和接收动态聚焦方式。实验所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为65dB。
对上述四个实验目标采用延时叠加算法(DAS)，多重幅度加权互相关算法(MAX)，广义旁瓣相消算法(GSC)以及基于多重幅度加权互相关的高分辨率自适应滤波算法(MAX-AHRF)进行对比成像实验。
图3给出了4种算法的点目标成像结果，从图3可以看出DAS算法成像质量最差，旁瓣伪影严重，相比于其他算法具有更高的横向伪影和更宽的主瓣宽度，难以区分目标点。MAX与DAS有相似的主瓣宽度，但旁瓣伪影得到一定的抑制。GSC算法在主瓣宽度上明显缩窄，横向分辨率明显提高，但是在旁瓣伪影的抑制上仍然表现不佳。MAX-AHRF与GSC有相似的主瓣宽度，但相对GSC进一步抑制了旁瓣。综上，MAX-AHRF的主瓣宽度和旁瓣等级为4种算法中成像质量最佳。
图4为4种算法点目标60mm处横向分辨率对比图，可以直观的看到MAX-AHRF的旁瓣抑制能力最强，对对比度的提高最高。为了更加直观的对比4种算法成像分辨率，表1给出了不同深度下4种算法-6dB的半峰值宽度(FWHM)数据对比。从表1可以看出，DAS分辨率最低，MAX的分辨率与DAS相似。GSC算法分辨率最高，MAX-AHRF的FWHM在不同深度处几乎与GSC相同，但由图4可知MAX-AHRF的旁瓣抑制能力与GSC相比有明显改进。综上，MAX-AHRF具有更好的综合性能。
表1不同深度下4种算法-6dB的FWHM对比
图5为4种算法吸声斑目标成像结果。从图5可以看出，4种算法对于暗斑轮廓都能较为准确清晰地成像，然而DAS的斑内部存在大量旁瓣伪影，这是由于其分辨率与对比度不足导致的。GSC算法也存在较明显伪影，但相对于DAS算法有了很大的改善，这是因为GSC算法虽然能够有效提高分辨率，但是在提升对比度方面仍然不足。MAX和MAX-AHRF中的吸声斑内部杂波和旁瓣被有效除去，其边缘轮廓更清晰，对于暗斑内的伪影抑制能力和提高对比度方面有显著改善。
表2不同算法成像性能指标比较
表2提供了不同算法的成像性能指标，可以看出DAS内部平均功率过高，即斑内伪影抑制能力弱，因此成像效果不佳。GSC算法的优势在于内部平均功率得到了改善，但是外部平均功率不佳，导致CR降低。MAX-AHRF的内部平均功率比GSC中的低大约50dB，这说明MAX-AHRF可提供有效的伪影抑制，使得吸声斑能被更清晰地识别。MAX-AHRF相对于GSC对比度明显提高，旁瓣和杂波抑制能力明显优于GSC，同时能够保持较好的背景散斑质量。综上，MAX-AHRF相对于其他算法可以得到更好的暗斑成像效果。
为了验证MAX-AHRF算法的鲁棒性，图6给出了4种算法分别加入信噪比为10dB/0dB/-10dB高斯白噪声的成像效果对比。从图6中可以看出，由于加了高斯白噪声，各算法成像对比度均有不同程度的下降。MAX-AHRF算法受噪声信号干扰最小，吸声斑轮廓最为清晰。
表3不同算法成像性能指标比较(10dB/0dB/-10dB)
表3给出了不同算法成像性能指标的对比，以更直观的反映MAX-AHRF算法的鲁棒性改善效果。从表中可以看出其中鲁棒性最差的GSC对比度分别下降了49.9％/69.0％/83.3％，DAS对比度分别下降了26.1％/50.8％/76.0％，MAX对比度分别下降了20.9％/48.3％/69.5％，而MAX-AHRF对比度仅分别下降了15.5％/29.5％/55.4％，且MAX-AHRF在增加噪声信号后成像对比度仍优于另几种算法。综上，MAX-AHRF相对于其他算法有更优越的鲁棒性。
图7给出了4种算法geabr_0数据成像结果。从图7中可以看出，DAS存在明显的暗斑旁瓣伪影，图像受到噪声干扰严重，GSC虽然提高了图像分辨率，但其对比度不佳。MAX对斑内的旁瓣伪影进行了有效的抑制，图像干净清晰，成像质量有显著的提高，但保持了DAS的分辨率。MAX-AHRF在保持GSC窄的点目标宽度的同时，大大提高对比度和鲁棒性，能够克服图像分辨率高、对比度高和鲁棒性强之间的权衡困难。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。
Claims (6)
1.一种基于多重幅度加权互相关的自适应滤波超声成像方法，其特征在于，该方法具体包括以下步骤：
S1：首先对超声阵元接收的回波信号进行放大，AD转换和延时聚焦处理以获得超声回波数据，然后将希尔伯特变换应用于回波数据，得到延时聚焦处理和希尔伯特变换之后的信号x(k)；
S2：构造多对互补的接收幅度加权函数，对各阵元延时信号进行加权得到两组波束形成信号，然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数；
S3：对得到的归一化互相关系数阈值化后取平均值，对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波；
S4：利用中值滤波后的系数确定自适应的空间平滑子阵长度，利用自适应子阵长度估计协方差矩阵，并进行对角加载处理；
S5：利用高分辨率自适应滤波器的结构计算出最优权值矢量，利用权值对超声阵列回波信号进行加权处理；
S6：将中值滤波后的系数作为加权因子乘上改进的高分辨率自适应滤波结构的输出，得到自适应波束形成的扫描线信号，采用优化后的扫描线信号进行最终成像。
2.根据权利要求1所述的自适应滤波超声成像方法，其特征在于，步骤S2具体包括以下步骤：
S21：对于有N个等间距阵元的传感器阵列，通过交替提取M个阵元来构建两个接收幅度加权函数Apod1n和Apod2n，其中n＝1,2...N，j＝1,2...M；
S22：将这两个接收幅度加权函数分别进一步划分为M个子幅度加权函数和其中j＝1,2...M；
S23：使用所得子幅度加权函数对各阵元延时信号进行处理，分别获得M个具有相位多样性的波束形成信号BF1j(k,h)和BF2j(k,h)，
其中，xn(k)为第n个阵元对应的延时回波信号，k为采样时刻；BF1j(k)和BF2j(k)表示一列扫描线信号的第k个采样时刻下的自适应波束信号，BF1j(k,h)和BF2j(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号；
然后计算两组波束形成信号两两之间的信号相似程度，并生成M个归一化互相关系数ρj(k,h)，表达式为：
其中，k表示成像区域中第h列的第k个采样时刻，A表示第k个采样时刻分别向上或向下取的样本长度，计算每个像素点时使用的互相关段总长度为2A+1个样本。
3.根据权利要求2所述的自适应滤波超声成像方法，其特征在于，步骤S3具体包括以下步骤：
S31：对得到的M个归一化互相关系数阈值化，并取平均值ρ(k,h)：
其中，δ为阈值，简记为/>
S32：对成像区域所有像素点对应的系数平均值所组成的加权矩阵进行中值滤波，对矩阵中的任意一点(k,h)将以该点为中心设计中值滤波器窗口ωkh，然后对窗口内包含数值按从小到大的顺序进行排列，取排列中值作为(k,h)点的滤波结果，窗口大小为3λ×10，其中λ为波长。
4.根据权利要求3所述的自适应滤波超声成像方法，其特征在于，步骤S4具体包括以下步骤：
S41：采用中值滤波后的系数自适应判别子阵长度L，定义自适应子阵长度/>为：
其中，N是阵元个数，round(·)是取整函数；整数的取值范围为1到N2；当/>时，/>为1；
S42：用自适应子阵长度代替固定子阵长度L，处理后的样本协方差矩阵R'(k)的计算公式为：
其中，(·)H表示共轭转置，表示第l个平滑子阵接收到的施加延时后的回波信号，/>
S43：通过以下计算公式对样本协方差矩阵R'(k)进行对角加载，得到对角加载后的样本协方差矩阵R'DL(k)：
R'DL(k)＝R'(k)+ε'·I'
其中，ε'＝Δ·trace(R'(k))，Δ为噪声和信号的比值，trace(·)表示矩阵的迹，I'表示单位阵。
5.根据权利要求4所述的自适应滤波超声成像方法，其特征在于，步骤S5具体包括以下步骤：
S51：通过以下公式计算高分辨率自适应滤波器结构中自适应权矢量w'a：
w'a＝(B'HR'DLB')-1B'HR'DLw'q
其中，非自适应权矢量w'q是的固定权向量，阻塞矩阵B'是/>维矩阵，且满足B'Hw'q＝0，(·)-1表示矩阵求逆运算；
S52：得到改进的高分辨率自适应滤波结构的最优权值矢量w'为：
w'＝w'q-B'w'a
S53：利用改进的最优权值矢量对超声阵列回波信号进行加权求和得到自适应波束形成信号：
其中，y(k)表示计算得到一列扫描线信号的第k个采样时刻下的自适应波束信号，w'H表示w'的共轭转置，w'表示改进的高分辨率自适应滤波算法的自适应加权矢量，表示第l个子阵在第k个采样时刻的回波向量；w'(k)表示第k个采样时刻下的w'。
6.根据权利要求5所述的自适应滤波超声成像方法，其特征在于，步骤S6中，得到自适应波束形成的扫描线信号y'(k,h)，表达式为：
其中，y(k,h)表示成像区域中第h列扫描线信号的第k个采样时刻下的自适应波束信号。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310891567.1A CN116934891A (zh) | 2023-07-19 | 2023-07-19 | 一种基于多重幅度加权互相关的自适应滤波超声成像方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310891567.1A CN116934891A (zh) | 2023-07-19 | 2023-07-19 | 一种基于多重幅度加权互相关的自适应滤波超声成像方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116934891A true CN116934891A (zh) | 2023-10-24 |
Family
ID=88390491
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310891567.1A Pending CN116934891A (zh) | 2023-07-19 | 2023-07-19 | 一种基于多重幅度加权互相关的自适应滤波超声成像方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116934891A (zh) |
-
2023
- 2023-07-19 CN CN202310891567.1A patent/CN116934891A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106510761B (zh) | 一种信噪比后滤波与特征空间融合的最小方差超声成像方法 | |
| Synnevag et al. | Benefits of minimum-variance beamforming in medical ultrasound imaging | |
| Nilsen et al. | Beamspace adaptive beamforming for ultrasound imaging | |
| CN102895000B (zh) | 一种基于自适应加权的双聚焦波束合成方法 | |
| CN110501423B (zh) | 一种基于频域分段的高分辨率最小方差超声成像方法 | |
| US20140378834A1 (en) | Continuous Transmit Focusing Method and Apparatus for Ultrasound Imaging System | |
| US20130109968A1 (en) | Ultrasonic imaging device | |
| CN108354627B (zh) | 一种提高帧频的超声波束形成方法 | |
| CN108836389B (zh) | 平面波相关点相干自适应波束合成成像方法 | |
| CN104970831A (zh) | 一种基于特征结构的广义旁瓣相消超声成像波束合成方法 | |
| CN102764139B (zh) | 基于特征空间分析和区域判别的医学超声波束形成方法 | |
| CN111856474B (zh) | 一种基于子阵的空时域条件相干系数超声成像方法 | |
| Li et al. | Blocked element compensation in phased array imaging | |
| US20180038947A1 (en) | Clutter suppression in ultrasonic imaging systems | |
| CN108309352A (zh) | 一种余弦变换域超声成像方法 | |
| CN110501711A (zh) | 一种基于乘幂法的低复杂度最小方差超声成像方法 | |
| CN107137111A (zh) | 一种超声波束形成方法 | |
| CN111278363B (zh) | 超声成像设备、系统及其超声造影成像的图像增强方法 | |
| CN113625286B (zh) | 基于相干特征的强稳健性截断相干系数超声波束形成方法 | |
| CN112120730A (zh) | 一种基于混合子空间投影的广义旁瓣相消超声成像方法 | |
| CN116934891A (zh) | 一种基于多重幅度加权互相关的自适应滤波超声成像方法 | |
| Ylitalo | A fast ultrasonic synthetic aperture imaging method: application to NDT | |
| CN114415187A (zh) | 基于子空间斜投影的强稳健性最小方差超声波束形成方法 | |
| CN109187771B (zh) | 一种融合特征值分解的低复杂度最小方差超声成像方法 | |
| CN113647978B (zh) | 一种带有截断因子的高鲁棒性符号相干系数超声成像方法 |
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 |