一种用于超快超声微血流成像的互相关降噪方法
阅读说明：本技术 一种用于超快超声微血流成像的互相关降噪方法 (Cross-correlation noise reduction method for ultrafast ultrasonic micro-blood flow imaging ) 是由 许凯亮 于 2021-09-17 设计创作，主要内容包括：本发明涉及超声波成像技术领域,具体涉及一种用于超快超声微血流成像的互相关降噪方法,旨在解决传统多普勒图像中信噪比低和杂波伪影影响图像质量的问题,该方法将多角度平面波射频回波信号分为两个或多个子组,各子组以相同的方法和参数进行正交解调、相干复合和杂波滤除,得到两组或多组动态血流信号,将每两组信号进行零时延互相关,并将所得数个互相关矩阵进行平均,可得到血流功率多普勒图像,将每两组信号进行固定时延互相关,并将所得数个互相关矩阵进行平均,可得到此时延下的多普勒频偏,继而得到血流彩色多普勒图像。本发明中互相关处理可以有效抑制图像中的随机噪声和非相关杂波伪影,同时保留高度相关的动态血流信号。(The invention relates to the technical field of ultrasonic imaging, in particular to a cross-correlation noise reduction method for ultrafast ultrasonic micro-blood flow imaging, which aims to solve the problems that the signal-to-noise ratio in a traditional Doppler image is low and clutter artifacts affect the image quality. The cross-correlation processing in the invention can effectively inhibit random noise and non-correlated clutter artifacts in the image, and simultaneously retain highly correlated dynamic blood flow signals.)
技术领域
本发明涉及超声波成像技术领域，具体涉及一种用于超快超声微血流成像的互相关降噪方法。
背景技术
微血管检测和成像技术对于评估生理状态和病理疾病有着至关重要的作用，如癌症、慢性肾脏疾病和炎症疾病等。而超声成像是一种安全、低成本且广泛应用的方法，可用于检测和评估体内血流，已在临床实践中应用了数十年。近年来，基于超快超声成像和先进的组织杂波滤波器(如基于特征滤波器、奇异值分解)的超声多普勒微血管成像技术在动物和人体上得到了广泛的应用。
现有公告号为CN111772676A的中国专利公开了一种超快超声多普勒脊髓微血流成像系统，包含硬件和软件两部分，其中硬件部分包含超声波发射和接收模块和实验器材模块，软件部分包含超快超声成像、波束合成、运动校准、杂波滤除、多普勒成像、求差和相关性分析等模块。首先，基于超快超声成像技术和多角度平面波复合成像的理论编写超声平面波发射与接收控制模块，并由计算机软件控制硬件设备发射和接收超声波；对接收到的回波数据进行处理，最终得到多普勒血流图像；并对血流的速度、方向等参数进行分析。该系统还提供了在脊髓加压、受损伤，刺激条件下对脊髓微血流进行成像的模式，可以进行脊髓功能分析和生理病理性分析。
然而，由于非聚焦波在组织深部的穿透能力较弱，基于平面波成像的微血管超快多普勒成像的一个重要缺点是信噪比不足，背景噪声对微血管的识别和成像影响较大，尤其是在深层组织。同时，超快超声微血流成像技术的另一个挑战是由体内强反射、缓慢运动、不均匀的组织散射产生的非相干杂波伪影，传统的杂波滤除方法无法对其进行有效地消除。
发明内容
因此，本发明的目的是提供一种用于超快超声微血流成像的互相关降噪方法，通过将不同偏转角度产生的血流信号进行相关，可以同时抑制超快超声微血流图像中的噪声和非相干杂波伪影，克服传统超快多普勒图像中信噪比低和杂波伪影影响图像质量的缺点。
本发明的上述技术目的是通过以下技术方案得以实现的：
一种用于超快超声微血流成像的互相关降噪方法，包括以下步骤：
S1、确定一个超声平面波的发射和接受序列，一个完整的序列包含多组子序列，在每一组子序列中，通过控制端激励超声换能器向目标成像区域内发射N个偏转角度的平面波；
S2、采样和存储成像区域反射或背散射的RF信号，并设定采样时长为T秒；
S3、对时间T秒内接收的每一组N个角度的RF信号按发射角度分为两个或多个子组；
S4、将每个子组中的RF信号进行正交解调后得到IQ信号；
S5、将每个子组中的IQ信号进行波束合成，重建出B-mode图像；
S6、将每个子组中通过不同发射角度得到的B-mode图像进行相干复合，每组可得到N张相干复合后的B-mode图像；
S7、将每个子组中相干复合后的B-mode图像进行杂波滤除，滤除运动缓慢的组织信号，得到动态血流信号；
S8、将每两个子组中的滤波后信号进行零时延互相关，并将所得数个互相关矩阵进行平均，得到血流功率多普勒图像；
S9、将每两个子组中的滤波后信号进行固定时延互相关，并将所得数个互相关矩阵进行平均，可得到此时延下的多普勒频偏，继而得到血流彩色多普勒图像。
可选地，所述步骤S3中，分组时不要求平均分组，各子组包含若干不同角度数的RF信号，且相邻子组之间的角度差不小于α。
可选地，所述步骤S4中，正交解调包括以下步骤：
S41、将RF信号与一个幅值为1，频率等于超声平面波发射频率的信号相乘进行下混叠；
S42、将下混叠后的信号经过一个低通滤波器，滤除负频谱和期望带宽以外的噪声；
S43、将低通滤波后的信号放大倍，保证信号的能量与原RF信号相同；
S44、数据抽取以得到IQ信号。
可选地，所述步骤S5中，波束合成的方法包括但不限于延时叠加算法、自适应波束合成法、空间复合法、频域-波数域迁移算法。
可选地，所述步骤S7中，杂波滤波的方法包括但不限于高通滤波算法、自适应滤波算法、特征值分解算法、鲁棒主成分分析算法、独立成分分析算法。
可选地，所述步骤S7中，特征值分解算法包括：
S71、将连续多帧的相干复合后图像构建为一个b*k大小的二维矩阵A，并利用特征值分解算法对二维矩阵A进行特征值分解，从而得到特征值矩阵λ和特征向量矩阵U：
E(A*AT)＝λ*U*UT；
式中，特征向量矩阵U为b*b的矩阵，特征值矩阵λ为b*b的对角阵，该对角阵中的对角元素即为矩阵特征值，将所有矩阵特征值由大到小排列为特征值序列；
S72、计算特征向量矩阵U中每个特征向量对应的平均多普勒频移fa：
式中，为第a个特征向量的自相关值，NF为连续多帧的相干复合后图像帧数，ea为第a个特征向量，PRF为阵列超声换能器对应的脉冲发射频率，arg{·}表示求解复数的辐角运算，fa为第a个特征向量对应的平均多普勒频移；
S73、依次判断平均多普勒频移fa是否处于区间[f1，f2]中，判断为否时，将平均多普勒频移fa对应的特征向量置零，从而得到新特征向量矩阵U1；
S74、基于新特征向量U1以及二维矩阵A，对相干复合后图像进行图像重构，从而得到b*k大小的矩阵Y：
S75、将矩阵Y重构为三维图像矩阵，作为连续多帧杂波滤除后的动态血流信号。
可选地，所述步骤S8中，将每两个子组中的滤波后信号进行零时延互相关，并将所得数个互相关矩阵进行平均，得到血流功率多普勒图像，包括如下步骤：
S(x,z,t)表示一组滤波后信号，共有Ng个子组，从Ng个子组中任意选出两个子组i和j，共有个组对，计算滤波后信号的零时延互相关的均值，得到血流功率多普勒图像，公式如下；
式中,PD(x,z)为所得的血流功率多普勒图像，Si(x,z,t)和Sj(x,z,t)分别为选中的两个子组的滤波后信号，NF为相干复合后图像帧数，*表示复数共轭。
可选地，所述步骤S9中，将每两个子组中的滤波后信号进行固定时延互相关，并将所得数个互相关矩阵进行平均，可得到此时延下的多普勒频偏，继而得到血流彩色多普勒图像，包括如下步骤：
S91、S(x,z,t)表示一组滤波后信号，共有Ng个子组，从Ng个子组中任意选出两个子组i和j，共有个组对，计算滤波后信号的固定时延互相关的均值，公式如下；
式中，R(x,z)为求均值后的互相关矩阵，NF为相干复合后图像帧数，lag为固定非零时延，幅值为α，辐角为
S92、根据互相关矩阵计算多普勒速度，公式如下：
其中，fc为发射超声信号的中心频率。
频移与相移有如下关系：
综上，多普勒速度可以由下式计算得到：
式中，CD(x,z)表示血流彩色多普勒图像，多普勒速度大于0表示向着探头方向运动，小于0表示远离探头方向运动。
本发明的有益效果：
本发明通过将不同偏转角度产生的血流信号进行相关，可以同时抑制血流图像中的噪声和非相干杂波伪影，克服了传统超快多普勒图像中信噪比低和杂波伪影影响图像质量的缺点，具有高图像信噪比和高血流/伪影强度比，提升了微血流成像质量。
附图说明
为了更清楚地说明本发明
具体实施方式
或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例的超快超声微血流成像的流程图；
图2为本发明实施例的一帧相干复合后的大鼠脑前囟处冠状面B-mode图像；
图3为本发明实施例采用互相关去噪方法，得到的大鼠脑微血流功率多普勒图像；
图4为本发明实施例不使用互相关去噪方法，得到的大鼠脑微血流功率多普勒图像；
图5为本发明实施例采用互相关去噪方法，得到的大鼠脑微血流彩色多普勒图像；
图6为本发明实施例不使用互相关去噪方法，得到的大鼠脑微血流彩色多普勒图像。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分的实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。
需要说明的是，本申请的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本申请的实施例。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
在本申请中，术语“上”、“下”、“左”、“右”、“前”、“后”、“顶”、“底”、“内”、“外”、“中”、“竖直”、“水平”、“横向”、“纵向”等指示的方位或位置关系为基于附图所示的方位或位置关系。这些术语主要是为了更好地描述本申请及其实施例，并非用于限定所指示的装置、元件或组成部分必须具有特定方位，或以特定方位进行构造和操作。
并且，上述部分术语除了可以用于表示方位或位置关系以外，还可能用于表示其他含义，例如术语“上”在某些情况下也可能用于表示某种依附关系或连接关系。对于本领域普通技术人员而言，可以根据具体情况理解这些术语在本申请中的具体含义。
另外，术语“多个”的含义应为两个以及两个以上。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
一种用于超快超声微血流成像的互相关降噪方法，如图1所示，包括以下步骤：
S1、确定一个超声平面波的发射和接受序列，一个完整的序列包含多组子序列，在每一组子序列中，通过控制端激励超声换能器向目标成像区域内发射N个偏转角度的平面波；
本实施例中优选为200组子序列，即连续发射200组平面波，每组包含-10°～10°之间均匀分布的九个倾斜平面波，并设定超声发射频率为15MHz，其他实施例亦可以选用发散波作为输出方式；
S2、采样和存储成像区域反射或背散射的RF信号，并设定采样时长为T秒；
本实施例中，将目标成像区域选定为大鼠脑前囟处冠状面血流，从而实验前先手术把大鼠颅骨去除，探头选用中心频率15MHz，128阵元的相控阵，另外采样频率为两倍奈奎斯特采样频率，即60MHz，其他实施例亦可以选用不低于2倍的信号发射频率；
S3、对时间T秒内接收的每一组N个角度的RF信号按发射角度分为两个或多个子组；
分组时不要求平均分组，各子组包含若干不同角度数的RF信号，且相邻子组之间的角度差不小于α，是为了降低不同子组之间的噪声及伪影的相关性，所以分组时应尽可能使组别之间角度相差较大；
本实施例中为了方便对比和查看，可分为两个子组，其中一组包含不大于0°的五个角度的RF信号，另一组包含大于0°的四个角度的RF信号，α为45°；
S4、将每个子组中的RF信号进行正交解调后得到IQ信号，以降低回波信号的带宽，同时不丢失基本信息；
具体地，正交解调包括以下步骤：
S41、将RF信号与一个幅值为1，频率等于超声平面波发射频率的信号相乘进行下混叠；
S42、将下混叠后的信号经过一个低通滤波器，滤除负频谱和期望带宽以外的噪声；
S43、将低通滤波后的信号放大倍，保证信号的能量与原RF信号相同；
S44、数据抽取以得到IQ信号，此过程中低通滤波器的截止频率为10MHz，数据抽取原数据的1/6；
S5、将每个子组中的IQ信号进行波束合成，重建出B-mode图像，其中波束合成算法为延时叠加算法、自适应波束合成法、空间复合法以及频域-波数域迁移算法中的一种或多种的组合，在实施例里优选为频域-波数域迁移算法进行波束合成，图2为本实施例的一帧相干复合后的大鼠脑前囟处冠状面B-mode图像；
S6、将每个子组中通过不同发射角度得到的B-mode图像进行相干复合，每组可得到N张相干复合后的B-mode图像；
S7、将每个子组中相干复合后的B-mode图像进行杂波滤除，滤除运动缓慢的组织信号，得到动态血流信号，其中，杂波滤除算法为高通滤波算法、自适应滤波算法、特征值分解算法、鲁棒主成分分析算法以及独立成分分析算法中一种或多种的组合；
本实施例中采用特征值分解算法，具体包括以下步骤：
S71、将连续多帧的相干复合后图像构建为一个b*k大小的二维矩阵A，并利用特征值分解算法对二维矩阵A进行特征值分解，从而得到特征值矩阵λ和特征向量矩阵U：
E(A*AT)＝λ*U*UT；
式中，特征向量矩阵U为b*b的矩阵，特征值矩阵λ为b*b的对角阵，该对角阵中的对角元素即为矩阵特征值，将所有矩阵特征值由大到小排列为特征值序列；
S72、计算特征向量矩阵U中每个特征向量对应的平均多普勒频移fa：
式中，为第a个特征向量的自相关值，NF为连续多帧的相干复合后图像帧数，ea为第a个特征向量，PRF为阵列超声换能器对应的脉冲发射频率，arg{·}表示求解复数的辐角运算，fa为第a个特征向量对应的平均多普勒频移；
S73、依次判断平均多普勒频移fa是否处于区间[f1，f2]中，判断为否时，将平均多普勒频移fa对应的特征向量置零，从而得到新特征向量矩阵U1；
S74、基于新特征向量U1以及二维矩阵A，对相干复合后图像进行图像重构，从而得到b*k大小的矩阵Y：
S75、将矩阵Y重构为三维图像矩阵，作为连续多帧杂波滤除后的动态血流信号；
S8、将每两个子组中的滤波后信号进行零时延互相关，并将所得数个互相关矩阵进行平均，得到血流功率多普勒图像，从而反映的图3为本实施例采用互相关去噪方法得到的大鼠脑微血流功率多普勒图像，图4为不使用互相关去噪方法，得到的大鼠脑微血流功率多普勒图像；
具体包括如下步骤：
S(x，z，t)表示一组滤波后信号，共有Ng个子组，从Ng个子组中任意选出两个子组i和j，共有个组对，计算滤波后信号的零时延互相关的均值，得到血流功率多普勒图像，公式如下；
式中,PD(x,z)为所得的血流功率多普勒图像，Si(x,z,t)和Sj(x,z,t)分别为选中的两个子组的滤波后信号，NF为相干复合后图像帧数，*表示复数共轭，经分析后，采用分组互相关去噪步骤后，图像信噪比提升约10-15dB，信号/伪影强度比提升约8-10dB；
S9、将每两个子组中的滤波后信号进行固定时延互相关，并将所得数个互相关矩阵进行平均，可得到此时延下的多普勒频偏，继而得到血流彩色多普勒图像，从而反映的图5为本实施例采用互相关去噪方法得到的大鼠脑微血流彩色多普勒图像，图6为不使用互相关去噪方法，得到的大鼠脑微血流彩色多普勒图像；
具体包括如下步骤：
S91、S(x,z,t)表示一组滤波后信号，共有Ng个子组，从Ng个子组中任意选出两个子组i和j，共有个组对，计算滤波后信号的固定时延互相关的均值，公式如下；
式中，R(x,z)为求均值后的互相关矩阵，NF为相干复合后图像帧数，lag为固定非零时延，幅值为α，辐角为
S92、根据互相关矩阵计算多普勒速度，公式如下：
其中，fc为发射超声信号的中心频率。
频移与相移有如下关系：
综上，多普勒速度可以由下式计算得到：
式中，CD(x,z)表示血流彩色多普勒图像，多普勒速度大于0表示向着探头方向运动，小于0表示远离探头方向运动，经分析后，采用分组互相关去噪步骤后，大幅度提高了微血流分辨率。
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。