CN106991708B

CN106991708B - 超声多普勒血流成像的处理方法及处理系统

Info

Publication number: CN106991708B
Application number: CN201710287911.0A
Authority: CN
Inventors: 郭建军; 凌涛
Original assignee: Vinno Technology Suzhou Co Ltd
Current assignee: Feiyinuo Technology Co ltd
Priority date: 2017-04-27
Filing date: 2017-04-27
Publication date: 2020-04-14
Anticipated expiration: 2037-04-27
Also published as: CN106991708A

Abstract

本发明提供一种超声多普勒血流成像的处理方法及处理系统，所述方法包括：获取回波RF信号；并对其做希尔伯特变换以形成正交IQ复信号；在慢时间方向上，对正交IQ复信号做滤波处理，以获得对应的滤波变换信号；在快时间方向，对获得的滤波变换信号做短时傅里叶变换，使其结果形成频域矩阵；获取当前频域矩阵中，每个快时间方向上的频率分量；根据一维自相关算法，估计当前频域矩阵中每列数据在其对应频率分量上的平均速度和能量；对当前频域矩阵中对应频率分量的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。本发明考虑组织衰减对超声信号的影响，基于宽带发射，提供更好的分辨率，进而提升了超声图像的质量。

Description

超声多普勒血流成像的处理方法及处理系统

技术领域

本发明属于医疗超声技术领域，主要涉及一种超声多普勒血流成像的处理方法及处理系统。

背景技术

彩色超声诊断仪(B超机)的彩色血流成像，以其独有的实时动态特性，成为现代医学不可或缺的辅助诊断的手段之一，在临床诊断中成为某些病症的判断标准。

当前彩色血流成像普遍采用自相关技术；根据多普勒效应，散射子运动导致的频移大小正比于发射信号频率和散射子运动速度的乘积。

传统的一维自相关技术，通过在回波信号慢时间方向上的自相关来估算散射子在发射信号频率上导致的相位差，从而估算出散射子的运动速度；该估计算法是建立在窄带信号的模型基础上，也就要求发射波形较长来保证符合窄带模型，其速度估计性能，尤其是对快速运动的散射子的速度估计能力，随发射信号的带宽的增加而显著恶化；同时，当发射波形较长时，血流散射子的成像分辨率会下降。

另外，由于超声在组织内传播时，存在的非线性频率衰减，随着声速在组织内的传播距离增加导致的频率成分不同衰减，传统一维自相关血流成像技术因为是基于固定的发射信号中心频率和固定带宽而进行计算，所以其估算速度的性能也随传播距离的增加而下降。

如上所述，对于传统的一维自相关血流成像技术，由于技术理论基础是基于窄带信号模型而得出，这就需要发射信号脉宽较长，从而影响了图像的分辨率。另外，超声在组织内存在非线性频率衰减现象，随着超声传播距离的增加，其信号的频带带宽和中心频率都会因此而偏移，传统的自相关技术，只能在基于固定的发射信号中心频率进行速度估算，所以随着深度增加时，其估算的速度误差也会增大。

相应的，当前公开的互相关血流成像技术和蝶型搜索成像技术基于宽带信号模型，解决了上述一维自相关技术与分辨率的矛盾，但也存在计算复杂度高，不易实时实现，和不能有效解决频率衰减导致的速度偏差问题。

发明内容

本发明的目的在于提供一种超声多普勒血流成像的处理方法及处理系统。

为了实现上述发明目的之一，本发明一实施方式的超声多普勒血流成像的处理方法，所述方法包括以下步骤：

S1、分别获取每个扫查包下、每条扫查线上、每个采样点对应的回波RF信号；分别对获取的回波RF信号做希尔伯特变换以形成正交IQ复信号；

S2、在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，以获得对应的滤波变换信号；

S3、在快时间方向，分别对获得的滤波变换信号做短时傅里叶变换，使其结果形成频域矩阵；

S4、获取当前所述频域矩阵中，每个快时间方向上的频率分量；

S5、根据一维自相关算法，估计当前频域矩阵中每列数据在其对应频率分量上的平均速度和能量；

S6、对当前频域矩阵中对应频率分量的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

作为本发明一实施方式的进一步改进，所述步骤S3具体包括：

P1、在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量；

P2、对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做L_fft阶短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数。

作为本发明一实施方式的进一步改进，所述步骤S4具体包括：

M1、获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；

M2、以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置；

M3、根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和；

M4、获取M*K频域矩阵中，每列数据对应的频率分量。

作为本发明一实施方式的进一步改进，所述步骤S5包括：根据一维自相关算法，估计M*K频域矩阵中每列数据在频率分量f(k)上的平均速度和能量，并将其集合标记为第一平均速度序列，和第一能量序列；

所述步骤S5和步骤S6之间，所述方法还包括：

N1、根据预设的高速度阈值，低速度阈值，噪声能量阈值和组织能量阈值对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；

若相同频率分量对应的第一平均速度序列中的平均速度值大于所述高速度阈值，对应能量序列中的能量值小于噪声能量阈值；

或相同频率分量对应的第一平均速度序列中的平均速度值小于所述低速度阈值，对应能量序列中的能量值大于组织能量阈值；

则剔除相同频率分量对应第一平均速度序列和第一能量序列对应的中的平均速度值和能量值，形成第二平均速度序列和第二能量序列。

作为本发明一实施方式的进一步改进，所述步骤N1后，所述方法还包括：

N2、获取第二平均速度序列对应的方差值及速度平均值，并根据其对第二平均速度序列和第二能量序列进行第二次预处理；

若第二平均速度序列中的任一平均速度值与第二平均速度序列的速度平均值差值的平方大于第二平均速度的方差值与第一预设常数值的乘积，则剔除相同频率分量对应第二平均速度序列和第二能量序列中的平均速度值和能量值，形成第三平均速度序列和第三能量序列。

为了实现上述发明目的之一，本发明一实施方式提供一种超声多普勒血流成像的处理系统，所述系统包括：信号获取模块，用于分别获取每个扫查包下、每条扫查线上、每个采样点对应的回波RF信号；分别对获取的回波RF信号做希尔伯特变换以形成正交IQ复信号；

滤波处理模块，在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，以获得对应的滤波变换信号；

傅里叶处理模块，在快时间方向，分别对获得的滤波变换信号做短时傅里叶变换，使其结果形成频域矩阵；

处理输出模块，用于获取当前所述频域矩阵中，每个快时间方向上的频率分量；

根据一维自相关算法，估计当前频域矩阵中每列数据在其对应频率分量上的平均速度和能量；

对当前频域矩阵中对应频率分量的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

作为本发明一实施方式的进一步改进，所述傅里叶处理模块具体用于：

在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量；

对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做L_fft阶短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数。

作为本发明一实施方式的进一步改进，所述处理输出模块具体用于：

获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；

以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置；

根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和；

获取M*K频域矩阵中，每列数据对应的频率分量。

作为本发明一实施方式的进一步改进，所述处理输出模块还用于：

根据一维自相关算法，估计M*K频域矩阵中每列数据在频率分量f(k)上的平均速度和能量，并将其集合标记为第一平均速度序列，和第一能量序列；

根据预设的高速度阈值，低速度阈值，噪声能量阈值和组织能量阈值对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；

获取第二平均速度序列对应的方差值及速度平均值，并根据其对第二平均速度序列和第二能量序列进行第二次预处理；

与现有技术相比，本发明的超声多普勒血流成像的处理方法及处理系统，考虑了组织衰减对超声信号的影响，基于宽带发射，提供更好的分辨率，同时，速度估计更稳定准确；提高了超声成像设备的方便性和使用效率，提升了超声图像的质量。

附图说明

图1是传统的成像系统的整体模块示意图；

图2是传统的超声多普勒血流成像的处理系统的模块示意图；

图3是本发明一实施方式中超声多普勒血流成像的处理方法的流程示意图；

图4是本发明一实施方式中超声多普勒血流成像的处理系统的模块示意图；

图5A、图5B是本发明一具体示例中分别采用传统的超声多普勒血流成像的处理方法与本发明的超声多普勒血流成像的处理方法获得的血流图像的对比示意图；

图5C、图5D是本发明一具体示例中分别采用传统的超声多普勒血流成像的处理方法与本发明的超声多普勒血流成像的处理方法获得的速度谱图的对比示意图。

具体实施方式

以下将结合附图所示的各实施方式对本发明进行详细描述。但这些实施方式并不限制本发明，本领域的普通技术人员根据这些实施方式所做出的结构、方法、或功能上的变换均包含在本发明的保护范围内。

需要说明的是，本发明主要应用于超声设备，相应的，所述待测物可为待测组织，在此不做详细赘述。

结合图1所示，多普勒成像系统的模块示意图；脉冲多普勒成像过程中；通过探头向组织中发射脉冲信号，所述脉冲信号经组织中反射形成超声信号经由探头换能器的不同基元转变为电模拟信号，通过前放模块放大，再由A/D数模转换模块转换为数字信号；各个不同基元的数字信号经过波束合成模块，合成为射频信号；射频信号经过固定频率的正交解调后，将正交解调结果I/Q信号送入相应的处理模块。

结合图2所示，传统的多普勒血流(CF)模式成像需经过如下过程：在正交解调模块中，回波射频RF信号先与发射中心频率为f0的正交信号相乘，再通过低通的基带滤波器，获得基带IQ信号，送入后面的壁滤波模块，壁滤波一般为高通滤波器，主要用来虑除信号中低速运动的组织信号。壁滤波后的信号，送入后面的速度能量估计模块计算速度及能量；再通过动态范围压缩模块进行对数压缩，最后进行血流显示。

结合图3所示，图3为本发明一实施方式中超声多普勒血流成像的处理方法的流程图，所述方法包括：S1、分别获取每个扫查包下、每条扫查线上、每个采样点对应的回波RF信号；分别对获取的回波RF信号做希尔伯特变换以形成正交IQ复信号；

本发明一具体示例中，为了方便描述，将每个扫查包下、每条扫查线上、每个采样点分别对应的正交IQ复信号分别以IQ(x,y)表示，所述(x,y)表示采样点的坐标。

需要理解的是，为了方便描述，以下出现的示例中仅以一个扫查包为例做具体介绍。相应的，对于一个扫查包中的第m条扫查线，m取值为[1,2,3……M]，M表示一个扫查包内扫查线的总数量，其对应的正交IQ复信号序列可表示为：{IQ(m,1),IQ(m,2),IQ(m,3),…,IQ(m,N)}，其中，N表示对应任一条扫查线，其在快时间方向上的采样点总数量，需要说明的是，在上述示例的实际应用中，所述扫查包的数量以及每个扫查包中扫查线的数量均可以根据需要具体调节，在此不做详细赘述。

进一步的，本发明一实施方式中，所述方法还包括：S2、在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，以获得对应的滤波变换信号。

所述滤波处理通常为采用高通滤波器对数据进行处理，其用于滤除正交IQ复信号中低速运动的组织信号。接续上述示例，本发明一具体实施方式中，一个扫查包下的，M条扫查线对应的正交IQ复信号序列依次为：

{IQ(1,1),IQ(1,2),IQ(1,3),…,IQ(1,N)}

{IQ(2,1),IQ(2,2),IQ(2,3),…,IQ(2,N)}

{IQ(3,1),IQ(3,2),IQ(3,3),…,IQ(3,N)}

……

{IQ(M,1),IQ(M,2),IQ(M,3),…,IQ(M,N)}；

进一步的，在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，获得每条扫查线上的每个扫查点分别对应的滤波变换信号序列依次为：

{W(1,1),W(1,2),W(1,3),…,W(1,N)}

{W(2,1),W(2,2),W(2,3),…,W(2,N)}

{W(3,1),W(3,2),W(3,3),…,W(3,N)}

……

{W(M,1),W(M,2),W(M,3),…,W(M,N)}，即获取当前数据包中的M*N数据矩阵对应的滤波变换信号。

进一步的，所述方法还包括：S3、在快时间方向，分别对获得的滤波变换信号做短时傅里叶变换，使其结果形成频域矩阵。

本发明一优选实施方式中，所述步骤S3具体包括：P1、在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量，L为偶数。

本发明优选实施方式中，所述L＝k*fc/f0，其中，k为常数，其可以根据需要自行设置，一般其取值可为1、2或3；fc表示系统回波信号的采样频率，f0表示发射信号的中心频率。

接续上述示例，本发明具体实施方式中，当前深度以d表示，则对应获得每条扫查线上的L个扫查点分别对应的滤波变换信号序列依次为：

{W(1,d-L/2+1),W(1,d-L/2+2),W(1,d-L/2+3),…,W(1,d-L/2)}

{W(2,d-L/2+1),W(2,d-L/2+2),W(2,d-L/2+3),…,W(2,d-L/2)}

{W(3,d-L/2+1),W(3,d-L/2+2),W(3,d-L/2+3),…,W(3,d-L/2)}

……

{W(M,d-L/2+1),W(M,d-L/2+2),W(M,d-L/2+3),…,W(M,d-L/2)}，即获取当前数据包中的M*L数据矩阵对应的滤波变换信号。

P2、对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数。通常情况下，短时傅里叶变换的阶数为2的幂，其幂的值可根据需要自行设定，本发明优选实施方式中，其幂取值范围为[5,8]。

接续上述示例，本发明具体实施方式中，做L_fft阶短时傅里叶变换，则获得每个扫查包对应的M*L_fft频域矩阵为：

{S_FFT(1,1),S_FFT(1,2),S_FFT(1,3),…,S_FFT(1,L_fft)}

{S_FFT(2,1),S_FFT(2,2),S_FFT(2,3),…,S_FFT(2,L_fft)}

{S_FFT(3,1),S_FFT(3,2),S_FFT(3,3),…,S_FFT(3,L_fft)}

……

{S_FFT(M,1),S_FFT(M,2),S_FFT(M,3),…,S_FFT(M,L_fft)}。

进一步的，本发明一实施方式中，所述方法还包括：S4、获取当前所述频域矩阵中，每个快时间方向上的频率分量。

本发明优选实施方式中，所述步骤S4具体包括：M1、获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；本发明具体实施方式中，获取M*L_fft频域矩阵靠近中间一行的频域数据，并获取其对应的幅值；本发明具体实施方式中，M为偶数，靠近中间一行的数据为第M/2行或M/2+1行。

继续上述示例，以取第M/2行数据为例做具体介绍。相应的，M*L_fft频域矩阵中第M/2行的频域数据表示为：

{S_FFT(M/2,1),S_FFT(M/2,2),S_FFT(M/2,3),…,S_FFT(M/2,L_fft)}，

其对应的幅值序列为：

{P_FFT(M/2,1),P_FFT(M/2,2),P_FFT(M/2,3),…,P_FFT(M/2,L_fft)}；

进一步的，以当前的幅值序列为搜索基础，在预设幅值范围内搜索幅值的最大值，将其确认为M*L_fft频域矩阵对应的中心频率位置，为了方便描述，将中心频率位置对应的列以N_fd表示。所述预设幅值范围为[N_fL，N_f0]

N_fL＝(f0–B/2)*L_fft/fc；

N_f0＝f0*L_fft/fc；其中，B为预设的信号带宽，f0、fc与上述示例表示相同。

M2、以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置。

本发明具体实施方式中，以当前的幅值序列为搜索基础，以N_fd为搜索中心，向其两侧依次进行搜索，至当前幅值与中心频率对应的幅值的比值符合设置的要求或数据为止，分别获得起止幅值，并将起止幅值分别对应的列以N_frxL和N_frxH表示。

M3、根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和。

本发明具体实施方式中，以M*L_fft频域矩阵为数据基础，获取从第N_frxL列起始至N_frxH列为止的M*K频域矩阵；其中，K表示感兴趣带宽的截止频率所对应列的总和；其取值k＝N_frxL,N_frxL+1,…,N_frxH。

本发明具体实施方式中，为了方便描述，将所述M*K频域矩阵中的任意一个数据以{S_FFT-1(m,k)}表示，k＝[N_frxL,N_frxL+1,…,N_frxH]＝[1,2,3……K],相应的，M*K频域矩阵可表示为：

{S_FFT-1(1,1),S_FFT-1(1,2),S_FFT-1(1,3),…,S_FFT-1(1,K)}

{S_FFT-1(2,1),S_FFT-1(2,2),S_FFT-1(2,),…,S_FFT-1(2,K)}

{S_FFT-1(3,1),S_FFT-1(3,2),S_FFT-1(3,3),…,S_FFT-1(3,K)}

……

{S_FFT-1(M,1),S_FFT-1(M,2),S_FFT-1(M,3),…,S_FFT-1(M,K)}，

上述M*K频域矩阵中的数据值与上述M*L_fft频域矩阵中的第N_frxL列至第N_frxH列的数据值相同，例如：M*K频域矩阵中的S_FFT-1(1,1)，与上述M*L_fft频域矩阵中S_FFT(1,N_frxL)数值相同，上述M*K频域矩阵中S_FFT-1(M,K)，与上述M*L_fft频域矩阵中S_FFT(M,N_frxH)数值相同，在此不做一一赘述。

进一步的，所述方法还包括：M4、获取M*K频域矩阵中，每列数据对应的频率分量。

本发明具体实施方式中，通过S_FFT-1(M,k)对应的具体数值为复数，可获得M*K频域矩阵中每个数据对应的实部和虚部值。如下所示，

S_FFT-1(m,k)＝R_m(k)+i*I_m(k)，其中，R_m(k)表示实部，I_m(k)表示虚部，i为虚部符号。

相应的，M*K频域矩阵中，第k列数据对应的频率分量可表示为：f(k)＝k*fc/L_fft。

进一步的，本发明一实施方式中，所述方法还包括：S5、根据一维自相关算法，估计当前频域矩阵中每列数据在其对应频率分量上的平均速度和能量。

本发明优选实施方式中，根据一维自相关算法，估计M*K频域矩阵中每列数据在频率分量f(k)上的平均速度和能量，并将其集合标记为第一平均速度序列，和第一能量序列；相应的，以公式表示如下：

其中，V0为超声在血流中的传播速度，fprf为血流扫查的脉冲重复频率，V_f(k)、P_f(k)分别表示M*K频域矩阵中第k列数据在其对应的频率分量f(k)上的平均速度和能量。

如此，获得M*K频域矩阵中、K列数据分别在其对应的频率分量上的第一平均速度序列，和第一能量序列，其可以表示为：

{V_f(1),V_f(2),V_f(3),……V_f(K)}

{P_f(1),P_f(2),P_f(3),……P_f(K)}，

由于k＝[N_frxL,N_frxL+1,…,N_frxH]＝[1,2,3……K]，对于M*L_fft频域矩阵，上式可表示为：

{V_f(N_frxL),V_f(N_frxL+1),V_f(N_frxL+2),……V_f(N_frxL)}

{P_f(N_frxL),P_f(N_frxL+1),P_f(N_frxL+2),……P_f(N_frxL)}。

进一步的，本发明一实施方式中，所述方法还包括：S6、对当前频域矩阵中对应频率分量上的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

本发明一优选实施方式中，为了使最终显示的图像更加清晰，所述步骤S6之前，所述方法还包括：N1、根据预设的高速度阈值，低速度阈值，噪声能量阈值和组织能量阈值对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；若相同频率分量对应的第一平均速度序列中的平均速度值大于所述高速度阈值，对应能量序列中的能量值小于噪声能量阈值；或相同频率分量对应的第一平均速度序列中的平均速度值小于所述低速度阈值，对应能量序列中的能量值大于组织能量阈值；则剔除相同频率分量对应第一平均速度序列和第一能量序列中的平均速度值和能量值，形成第二平均速度序列和第二能量序列。

本发明一具体实施方式中，接续上述示例，对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；

若V_f(k)>TV_h且P_f(k)<TP_l，剔除V_f(k)和P_f(k)；

若V_f(k)<TV_l且P_f(k)>TP_h，剔除V_f(k)和P_f(k)，

其中，TV_h表示高速度阈值，TV_l表示低速度阈值，TP_l表示噪声能量阈值，TP_h表示组织能量阈值。

相应的，形成的第二平均速度序列和第二能量序列，表示为：

{V_f-1(1),V_f-1(2),V_f-1(3),……V_f-1(P)}

{P_f-1(1),P_f-1(2),P_f-1(3),……P_f-1(P)}，其中，P小于等于K，第二平均速度序列中的任一平均速度以P_f-1(p)表示，p＝1,2,3……P。

优选的，本发明一实施方式中，所述方法还包括：N2，获取第二平均速度序列对应的方差值及速度平均值，并根据其对第二平均速度序列和第二能量序列进行第二次预处理；若第二平均速度序列中的任一平均速度值与第二平均速度序列的速度平均值差值的平方大于第二平均速度的方差值与第一预设常数值的乘积，则剔除相同频率分量对应第二平均速度序列和第二能量序列中的平均速度值和能量值，形成第三平均速度序列和第三能量序列。

若(V_f-1(p)-V_ave)²＞D_ave*K_d，剔除V_f-1(p)和P_f-1(p)；其中，V_f-1(p)表示第二平均速度序列中的任一平均速度值，P_f-1(p)表示第二能量序列中的任一能量值，V_ave表示第二平均速度序列的速度平均值，D_ave表示第二平均速度序列的方差值，K_d为预设常识值。

相应的，形成的第三平均速度序列，和第三能量序列，表示为：

{V_f-2(1),V_f-2(2),V_f-2(3),……V_f-2(Q)}

{P_f-2(1),P_f-2(2),P_f-2(3),……P_f-2(Q)}，其中，Q小于等于P，第三平均速度序列中的任一平均速度以P_f-2(q)表示，q＝1,2,3……Q。

进一步的，对第三平均速度序列和第三能量序列进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

本实施方式中，用于最终的频谱显示的最终能量和最终速度可分别表示为：

其中，P_d表示最终能量，V_d表示最终速度。

进一步的，将获得的最终能量和最终速度进行动态压缩，以进行输出显示，在此不做详细赘述。

结合图5A-图5D所示，图5A、图5C为采用传统的超声多普勒血流成像的处理方法分别获得的血流图像和速度谱图；图5B、图5D为采用本发明的超声多普勒血流成像的处理方法分别获得的血流图像和速度谱图；

通过比对图5A和图5B可知：对组织中两支距离较近且细的血管进行血流成像过程中，采用本发明获得的血流图像的分辨率更高，且能明显区分出组织中两支血管。

通过比对图5C和图5D可知：其中横轴代表深度方向，纵轴代表速度大小，当发射信号为较短时，对于同样的2个散射子，采用传统的超声多普勒血流成像的处理方法获得的速度估计不稳定，在深度方向分布较宽且误差较大，而采用本发明的超声多普勒血流成像的处理方法获得的速度估计更准确且分辨率更高。

结合图4所示，本发明一实施方式中提供的超声多普勒血流成像的处理系统，所述系统包括：信号获取模块100、滤波处理模块200、傅里叶处理模块300、处理输出模块400。

信号获取模块100用于分别获取每个扫查包下、每条扫查线上、每个采样点对应的回波RF信号；分别对获取的回波RF信号做希尔伯特变换以形成正交IQ复信号；

本发明一实施方式中，滤波处理模块200用于在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，以获得对应的滤波变换信号。

{IQ(1,1),IQ(1,2),IQ(1,3),…,IQ(1,N)}

{IQ(2,1),IQ(2,2),IQ(2,3),…,IQ(2,N)}

{IQ(3,1),IQ(3,2),IQ(3,3),…,IQ(3,N)}

……

{IQ(M,1),IQ(M,2),IQ(M,3),…,IQ(M,N)}；

进一步的，滤波处理模块200在慢时间方向上，分别对获得的正交IQ复信号做滤波处理，获得每条扫查线上的每个扫查点分别对应的滤波变换信号序列依次为：

{W(1,1),W(1,2),W(1,3),…,W(1,N)}

{W(2,1),W(2,2),W(2,3),…,W(2,N)}

{W(3,1),W(3,2),W(3,3),…,W(3,N)}

……

本发明一实施方式中，傅里叶处理模块300用于在快时间方向，分别对获得的滤波变换信号做短时傅里叶变换，使其结果形成频域矩阵。

本发明一优选实施方式中，傅里叶处理模块300具体用于：在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量，L为偶数。

{W(1,d-L/2+1),W(1,d-L/2+2),W(1,d-L/2+3),…,W(1,d-L/2)}

{W(2,d-L/2+1),W(2,d-L/2+2),W(2,d-L/2+3),…,W(2,d-L/2)}

{W(3,d-L/2+1),W(3,d-L/2+2),W(3,d-L/2+3),…,W(3,d-L/2)}

……

傅里叶处理模块300还用于：对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数。通常情况下，短时傅里叶变换的阶数为2的幂，其幂的值可根据需要自行设定，本发明优选实施方式中，其幂取值范围为[5,8]。

{S_FFT(1,1),S_FFT(1,2),S_FFT(1,3),…,S_FFT(1,L_fft)}

{S_FFT(2,1),S_FFT(2,2),S_FFT(2,3),…,S_FFT(2,L_fft)}

{S_FFT(3,1),S_FFT(3,2),S_FFT(3,3),…,S_FFT(3,L_fft)}

……

{S_FFT(M,1),S_FFT(M,2),S_FFT(M,3),…,S_FFT(M,L_fft)}。

本发明一实施方式中，处理输出模块400用于获取当前所述频域矩阵中，每个快时间方向上的频率分量。

本发明优选实施方式中，处理输出模块400具体用于：获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；本发明具体实施方式中，处理输出模块400获取M*L_fft频域矩阵靠近中间一行的频域数据，并获取其对应的幅值；本发明具体实施方式中，M为偶数，靠近中间一行的数据为第M/2行或M/2+1行。

{S_FFT(M/2,1),S_FFT(M/2,2),S_FFT(M/2,3),…,S_FFT(M/2,L_fft)}，

其对应的幅值序列为：

{P_FFT(M/2,1),P_FFT(M/2,2),P_FFT(M/2,3),…,P_FFT(M/2,L_fft)}；

进一步的，处理输出模块400以当前的幅值序列为搜索基础，在预设幅值范围内搜索幅值的最大值，将其确认为M*L_fft频域矩阵对应的中心频率位置，为了方便描述，将中心频率位置对应的列以N_fd表示。所述预设幅值范围为[N_fL，N_f0]

N_fL＝(f0–B/2)*L_fft/fc；

进一步的，处理输出模块400还用于：以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置。

本发明具体实施方式中，处理输出模块400以当前的幅值序列为搜索基础，以N_fd为搜索中心，向其两侧依次进行搜索，至当前幅值与中心频率对应的幅值的比值符合设置的要求或数据为止，分别获得起止幅值，并将起止幅值分别对应的列以N_frxL和N_frxH表示。

进一步的，处理输出模块400还用于：根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和。

本发明具体实施方式中，处理输出模块400以M*L_fft频域矩阵为数据基础，获取从第N_frxL列起始至N_frxH列为止的M*K频域矩阵；其中，K表示感兴趣带宽的截止频率所对应列的总和；其取值k＝N_frxL,N_frxL+1,…,N_frxH。

{S_FFT-1(1,1),S_FFT-1(1,2),S_FFT-1(1,3),…,S_FFT-1(1,K)}

{S_FFT-1(2,1),S_FFT-1(2,2),S_FFT-1(2,),…,S_FFT-1(2,K)}

{S_FFT-1(3,1),S_FFT-1(3,2),S_FFT-1(3,3),…,S_FFT-1(3,K)}

……

{S_FFT-1(M,1),S_FFT-1(M,2),S_FFT-1(M,3),…,S_FFT-1(M,K)}，

进一步的，处理输出模块400还用于：获取M*K频域矩阵中，每列数据对应的频率分量。

本发明具体实施方式中，处理输出模块400通过S_FFT-1(M,k)对应的具体数值为复数，可获得M*K频域矩阵中每个数据对应的实部和虚部值。如下所示，

本发明一实施方式中，处理输出模块400还用于：根据一维自相关算法，估计当前频域矩阵中每列数据在其对应频率分量上的平均速度和能量。

本发明优选实施方式中，处理输出模块400根据一维自相关算法，估计M*K频域矩阵中每列数据在频率分量f(k)上的平均速度和能量，并将其集合标记为第一平均速度序列，和第一能量序列；相应的，以公式表示如下：

{V_f(1),V_f(2),V_f(3),……V_f(K)}

{P_f(1),P_f(2),P_f(3),……P_f(K)}，

{V_f(N_frxL),V_f(N_frxL+1),V_f(N_frxL+2),……V_f(N_frxL)}

{P_f(N_frxL),P_f(N_frxL+1),P_f(N_frxL+2),……P_f(N_frxL)}。

进一步的，本发明一实施方式中，处理输出模块400还用于：对当前频域矩阵中对应频率分量上的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

本发明一优选实施方式中，为了使最终显示的图像更加清晰，处理输出模块400还用于：根据预设的高速度阈值，低速度阈值，噪声能量阈值和组织能量阈值对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；若相同频率分量对应的第一平均速度序列中的平均速度值大于所述高速度阈值，对应能量序列中的能量值小于噪声能量阈值；或相同频率分量对应的第一平均速度序列中的平均速度值小于所述低速度阈值，对应能量序列中的能量值大于组织能量阈值；则剔除相同频率分量对应第一平均速度序列和第一能量序列中的平均速度值和能量值，形成第二平均速度序列和第二能量序列。

本发明一具体实施方式中，接续上述示例，处理输出模块400对M*K频域矩阵对应的第一平均速度序列和第一能量序列进行第一次预处理；

若V_f(k)>TV_h且P_f(k)<TP_l，剔除V_f(k)和P_f(k)；

若V_f(k)<TV_l且P_f(k)>TP_h，剔除V_f(k)和P_f(k)，

{V_f-1(1),V_f-1(2),V_f-1(3),……V_f-1(P)}

优选的，本发明一实施方式中，处理输出模块400还用于：获取第二平均速度序列对应的方差值及速度平均值，并根据其对第二平均速度序列和第二能量序列进行第二次预处理；若第二平均速度序列中的任一平均速度值与第二平均速度序列的速度平均值差值的平方大于第二平均速度的方差值与第一预设常数值的乘积，则剔除相同频率分量对应第二平均速度序列和第二能量序列中的平均速度值和能量值，形成第三平均速度序列和第三能量序列。

{V_f-2(1),V_f-2(2),V_f-2(3),……V_f-2(Q)}

进一步的，处理输出模块400还用于：对第三平均速度序列和第三能量序列进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度。

其中，P_d表示最终能量，V_d表示最终速度。

进一步的，本发明的超声多普勒血流成像的处理系统还包括：动态压缩模块(未图示)以及显示模块(未图示)，所述动态压缩模块将获得的最终能量和最终速度进行动态压缩，所述显示模块用于对经过动态压缩模块处理的数据结果输出显示，在此不做详细赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统中模块的具体工作过程，可以参考前述方法实施方式中的对应过程，在此不再赘述。

综上所述，本发明的超声多普勒血流成像的处理方法及处理系统，考虑了组织衰减对超声信号的影响，基于宽带发射，提供更好的分辨率，同时，速度估计更稳定准确；提高了超声成像设备的方便性和使用效率，提升了超声图像的质量。

为了描述的方便，描述以上装置时以功能分为各种模块分别描述。当然，在实施本申请时可以把各模块的功能在同一个或多个软件和/或硬件中实现。

应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施方式中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

上文所列出的一系列的详细说明仅仅是针对本发明的可行性实施方式的具体说明，它们并非用以限制本发明的保护范围，凡未脱离本发明技艺精神所作的等效实施方式或变更均应包含在本发明的保护范围之内。

Claims

1.一种超声多普勒血流成像的处理方法，其特征在于，所述方法包括以下步骤：

S6、对当前频域矩阵中对应频率分量的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度；

其中，所述步骤S3具体包括：P1、在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量；P2、对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做L_fft阶短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数；

所述步骤S4具体包括：M1、获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；M2、以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置；M3、根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和；M4、获取M*K频域矩阵中，每列数据对应的频率分量。

2.根据权利要求1所述的超声多普勒血流成像的处理方法，其特征在于，

所述步骤S5包括：根据一维自相关算法，估计M*K频域矩阵中每列数据在频率分量f(k)上的平均速度和能量，并将其集合标记为第一平均速度序列，和第一能量序列；

所述步骤S5和步骤S6之间，所述方法还包括：

3.根据权利要求2所述的超声多普勒血流成像的处理方法，其特征在于，所述步骤N1后，所述方法还包括：

4.一种超声多普勒血流成像的处理系统，其特征在于，所述系统包括：

信号获取模块，用于分别获取每个扫查包下、每条扫查线上、每个采样点对应的回波RF信号；分别对获取的回波RF信号做希尔伯特变换以形成正交IQ复信号；

对当前频域矩阵中对应频率分量的平均速度和能量进行线性拟合，获取用于最终的频谱显示的最终能量和最终速度；

其中，所述傅里叶处理模块具体用于：

在同一深度下，对应每个扫查包、在每条扫查线的快时间方向上分别获取L个采样点对应的滤波变换信号，L小于等于N，N表示每条扫查线对应的原始采样点数量；对每条扫查线对应的L个滤波变换信号序列进行加窗补零后，做L_fft阶短时傅里叶变换，使其结果形成对应一个扫查包的M*L_fft频域矩阵，M表示每个扫查包设置的扫查线数量；L_fft表示短时傅里叶变换的阶数；

所述处理输出模块具体用于：获取当前深度下，M*L_fft频域矩阵对应的中心频率位置；以获取的中心频率位置为中心，获取当前深度下，感兴趣带宽的截止频率位置；根据获取的感兴趣带宽的截止频率位置获取其对应的M*K频域矩阵，K表示感兴趣带宽的截止频率所对应的列的总和；获取M*K频域矩阵中，每列数据对应的频率分量。

5.根据权利要求4所述的超声多普勒血流成像的处理系统，其特征在于，所述处理输出模块还用于：

6.根据权利要求5所述的超声多普勒血流成像的处理系统，其特征在于，所述处理输出模块还用于：