波束域广义旁瓣相消超声成像方法与流程

本发明属于超声成像
技术领域：
，涉及一种波束域广义旁瓣相消超声成像方法。
背景技术：
超声成像具有操作简单，安全性高和成像迅速的特点在医学和工业方面有广泛的应用，目前应用最为广泛的是传统的das成像。它具有成像速度快，但分辨率和对比度较差。为了解决这一问题，有学者提出了自适应波束形成算法，即根据回波数据动态的计算加权系数，最常见的是最小方差(mv)算法，其原理是保持期望方向增益不变，使得阵列输出的能量达到最小，这在一定程度上提高了传统das算法的成像效果。在此基础之上，一些学者相继提出了对角加载(diagonalloading,dl)技术来提高算法的鲁棒性，空间平滑法来消除回波信号之间的相关性，以及特征空间法来进一步提高超声成像的分辨率与对比度。由于超声成像检测对象的不同，声速也存在差异，因为声速不均匀性引起的相位畸变也是导致超声成像质量下降的一个重要来源。为了解决这个问题，asl等提出了相干系数(coherencefactor，cf)作为成像聚焦质量的指标，可以有效压缩旁瓣提高成像的分辨率；li等提出广义相干系数(generalizedcoherencefactor,gcf)作为自适应加权系数应用于成像算法，并且应用最为广泛。广义旁瓣相消(gsc)算法作为mv的等价结构，它分为上下两个通道，上通道为非自适应通道；下通道为自适应通道，通过阻塞矩阵阻止期望信号的进入，最后通过叠加来消除非期望信号，提高成像质量。成像质量与成像时效性一直是广大学者的追求目标，因此减小算法的复杂度也是近年来的研究热点。asl和mahloojifar提出了通过将协方差构造成toeplitz矩阵，来减少矩阵求逆的运算量进而降低算法的复杂度；有些学者也提出了用qr分解算法来实现协方差矩阵求逆的快速运算；也有学者通过householder多级维纳滤波来避免输入信号自相关矩阵的估计进而降低算法的复杂度。基于以上原因，本发明提出了利用离散余弦变换(discretecosinetransform,dct)来构造转换矩阵，将超声阵列回波信号从高维的阵元域转换到低维的波束域，进而减少协方差矩阵的维数降低算法复杂度，变换后的回波信号我们取其低频部分，即可获得原始回波信号中的大部分能量，所以本文提出的算法在成像效果基本不变的情况下，能够明显缩短了成像时间，加快成像速度。为了验证其有效性，本发明对点目标以及吸声斑目标进行了仿真成像，并用geabr_0数据进行了实验成像，将波束域广义旁瓣相消(beamdomaingeneralizedsidelobecanceller,bgsc)与das、gsc，以及esbgsc算法在分辨率、对比度以及成像时间方面进行了比较。技术实现要素：有鉴于此，本发明的目的在于提供一种波束域广义旁瓣相消超声成像方法，该方法能够在成像效果基本不变，有些参数甚至优于传统广义旁瓣相消的情况下，较大地缩短算法的运行时间，有效克服了传统广义旁瓣相消算法由于运行时间较长而限制了其应用范围的问题。本发明的目的是研究一种波束域广义旁瓣相消算法，在降低算法的复杂度的同时，减少超声成像的运行时间。为达到上述目的，本发明提供如下技术方案：一种波束域广义旁瓣相消超声成像方法，该方法包括以下步骤：s1：对超声阵元接收的回波信号进行放大滤波处理与ad转换，获得超声成像所需的阵列回波信号；s2：将接收阵元依次划分为具有一个重叠阵元的子阵列，对相应接收子阵的回波信号进行空间平滑处理，得到阵元域样本协方差矩阵；s3：利用离散余弦变化构造转换矩阵，将空间平滑后的阵元域子阵列回波信号转换为对应的波束域回波信号；s4：将传统广义旁瓣相消算法中的非自适应权矢量和阻塞矩阵利用转换矩阵变换到波束域；s5：利用每个子阵列对应的波束域回波信号获得波束域的样本协方差矩阵，并通过对角加载技术提高算法稳定性；s6：通过获得的波束域非自适应权矢量，阻塞矩阵以及波束域的协方差矩阵，得到波束域广义旁瓣相消算法权矢量的表达式；s7：利用波束域广义旁瓣相消算法权矢量对波束域回波信号进行加权求和，得到波束域广义旁瓣相消算法的输出信号。进一步，所述步骤s2具体包括：s21：将n个阵元分为(n-l+1)个子阵列，其中每个子阵有l个阵元，设表示第l个阵元域子阵列接收的回波信号：其中e表示阵元域下标以区别波束域，k表示第k个采样点，n表示共有n个阵元，l表示第l个子阵列，表示第l个阵元在第k个采样点的回波数据，以此类推和分别表示第l+1个阵元和第l+l-1个阵元在第k个采样点的回波数据，[·]mt为转置运算。s22：利用s21中得到的子阵列回波信号，通过以下公式得到阵元域样本协方差矩阵：其中表示s21中得到的第l个阵元域子阵列接收的回波信号，表示阵元域样本协方差矩阵，矩阵维数为l×l，[·]h表示共轭转置运算。进一步，所述步骤s3具体包括以下步骤：s31：通过以下公式构造(1+p)×l维转换矩阵，须满足(1+p)≤l：矩阵t满足tth＝i，其中i为单位阵，th表示矩阵t的共轭转置矩阵，tm,n表示矩阵t第m行，第n列的值，l为每个子阵的阵元数，p表示协方差矩阵的降秩参数；s32：利用s31中得到波束域转换矩阵后，将s2中得到的阵元域子阵列回波信号转换为波束域回波信号，第l个子阵列按下式进行转换：其中表示波束域第l个子阵列回波数据维数为(1+p)×1的向量，b为波束域下标以区别阵元域，k表示第k个采样点，t(1+p)×l表示维数为(1+p)×l的波束域转换矩阵，表示s2中得到的阵元域第l个子阵列回波数据维数为l×1的向量；s33：利用s32得到的波束域子阵列回波信号，对每个子阵列信号进行空间平滑处理，即利用均匀线阵的平移不变性，将得到的每个(1+p)×1维波束域子阵列进行叠加处理，从而得到波束域回波信号：其中xb(k)为空间平滑后的波束域回波信号，表示波束域第1子阵列的回波信号，和与其类似分别表示第2和第n-l+1子阵列回波信号。进一步，其特征在于：所述步骤s4具体包括以下步骤：s41：波束域广义旁瓣相消算法对应的阻塞矩阵按下式计算：bb＝tb其中，b为传统广义旁瓣相消算法的阻塞矩阵，t为波束域转换矩阵；s42：波束域广义旁瓣相消算法中非自适应权矢量为：wbq＝twq其中，wq为传统广义旁瓣相消中的非自适应权矢量。进一步，其特征在于：所述步骤s5具体包括以下步骤：s51：利用每个子阵列对应的波束域回波信号按下式计算，获得波束域的样本协方差矩阵：即由波束域的样本协方差矩阵代替s22中得到的阵元域样本协方差矩阵矩阵维数由l×l变为(1+p)×(1+p)，其中表示波束域第l个子阵列回波数据；s52：利用对角加载技术增加算法稳定性：其中i为单位矩阵，rb(k)为对角加载后的波束域协方差矩阵，λ为对角加载系数，且满足δ为常数，且满足在本发明中进一步，其特征在于：所述步骤s6具体包括以下步骤：s61：传统广义旁瓣相消算法权矢量的表达式：w＝wq-bwa，s.tbhwq＝0式中，wq是n×1的固定权矢量，b为n×(n-1)的阻塞矩阵，bh是b的共轭转置运算，由于wq为固定的向量，且b与向量wq正交，所以广义旁瓣相消将最小方差算法中的约束优化问题转换为非约束问题，按下式计算得到自适应权矢量的最优解wa：wa＝(bhrb)-1bhrwqs62：波束域广义旁瓣相消算法权矢量wb按下式计算：式中，wbq为波束域非自适应权矢量，bb为波束域阻塞矩阵，wba为波束域自适应权矢量；其中自适应权矢量进一步，所述步骤s7中的自适应输出波束信号按下式计算：式中，为权矢量wb的共轭转置，xb(k)为子阵列波束域回波信号空间平滑后的波束域信号，yb(k)即为波束域的输出波束信号。进一步，所述子阵阵元数目l的取值为l≤n/2。本发明的有益效果在于：本发明采用了一种波束域广义旁瓣相消超声成像算法，通过离散余弦变换构造转换矩阵，降低自相关矩阵的维数，进而降低算法的复杂度。通过点目标，吸声斑仿真以及geabr_0数据进行实验，从成像效果图和分辨率对比图说明bgsc算法与gsc算法在成像结果上基本一致，吸声斑成像的对比度甚至略优于gsc算法，并且由于减少了算法的复杂度，鲁棒性也有所增加；在运算时间方面点目标仿真比gsc在运行时间上缩短23.6％，较esbgsc缩短58.8％；吸声斑仿真较gsc缩短28.1％，相比esbgsc缩短55.7％；最后我们采用michigan大学生物医学实验室提供的geabr_0数据进行实验，得出的结论与点目标和吸声斑仿真结果类似，进一步验证了本发明采用算法的有效性和实用性。附图说明为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：图1为本发明所述方法的流程图；图2为广义旁瓣相消器结构图；图3为空间平滑原理示意图；图4为4种算法点目标仿真成像对比图；图5为轴向距离50mm处点目标侧向分辨率对比图；图6为加噪声情况下4种算法点目标成像对比图；图7为轴向距离50mm处加噪声情况下点目标侧向分辨率对比图；图8为4种算法吸声斑成像结果对比图；图9为4种算法geabr_0数据实验结果对比图；图10为4种算法geabr_0数据实验分辨率对比图。具体实施方式下面将结合附图，对本发明的优选实例进行详细的描述。图1为本发明所述方法的流程图，图2为广义旁瓣相消器结构图，如图1、2所示，本发明提供一种在超声成像中的波束域广义旁瓣相消算法，包括以下步骤：步骤s1：对超声阵元接收的回波信号进行放大处理与ad转换，以获得超声成像所需的阵元域回波数据；步骤s2：将接收阵元依次划分为具有一个重叠阵元的子阵，对相应接收子阵的回波信号进行空间平滑处理；图3给出了前后向空间平滑原理示意图，具体包括以下步骤：s21：将n个阵元分为(n-l+1)个子阵列，其中每个子阵有l个阵元，设表示第l个阵元域子阵列接收的回波信号：其中e表示阵元域下标以区别波束域，k表示第k个采样点，n表示共有n个阵元，l表示第l个子阵列，表示第l个阵元在第k个采样点的回波数据，以此类推和分别表示第l+1个阵元和第l+l-1个阵元在第k个采样点的回波数据，[·]mt为转置运算。s22：利用s21中得到的子阵列回波信号，通过以下公式得到阵元域样本协方差矩阵：其中表示s21中得到的第l个阵元域子阵列接收的回波信号，表示阵元域样本协方差矩阵，矩阵维数为l×l，[·]h表示共轭转置运算。步骤s3：利用离散余弦变化构造波束域转换矩阵，将空间平滑后的阵元域回波信号转换为波束域的回波信号，具体步骤如下：s31：通过以下公式构造(1+p)×l维转换矩阵，须满足(1+p)≤l：矩阵t满足tth＝i，其中i为单位阵，th表示矩阵t的共轭转置矩阵，tm,n表示矩阵t第m行，第n列的值，l为每个子阵的阵元数，p表示协方差矩阵的降秩参数；须满足(1+p)≤l才能有效降低自相关矩阵的维数，进而减小算法的复杂度，在本发明的实例中选取p＝8，子阵列包含的阵元数l＝32；s32：利用s31中得到波束域转换矩阵后，将s2中得到的阵元域子阵列回波信号转换为波束域回波信号，第l个子阵列按下式进行转换：其中表示波束域第l个子阵列回波数据维数为(1+p)×1的向量，b为波束域下标以区别阵元域，k表示第k个采样点，t(1+p)×l表示维数为(1+p)×l的波束域转换矩阵，表示s2中得到的阵元域第l个子阵列回波数据维数为l×1的向量；s33：利用s32得到的波束域子阵列回波信号，对每个子阵列信号进行空间平滑处理，即利用均匀线阵的平移不变性，将得到的每个(1+p)×1维波束域子阵列进行叠加处理，从而得到波束域回波信号：其中xb(k)为空间平滑后的波束域回波信号，表示波束域第1子阵列的回波信号，和与其类似分别表示第2和第n-l+1子阵列回波信号。s4：将传统广义旁瓣相消算法中的非自适应权矢量和阻塞矩阵利用转换矩阵变换到波束域；具体包括以下步骤：s41：波束域广义旁瓣相消算法对应的阻塞矩阵按下式计算：bb＝tb其中，b为传统广义旁瓣相消算法的阻塞矩阵，t为波束域转换矩阵；s42：波束域广义旁瓣相消算法中非自适应权矢量为：wbq＝twq其中，wq为传统广义旁瓣相消中的非自适应权矢量，t为波束域转换矩阵。s5：利用每个子阵列对应的波束域回波信号获得波束域的样本协方差矩阵，并通过对角加载技术提高算法稳定性；具体包括以下步骤：s51：利用每个子阵列对应的波束域回波信号按下式计算，获得波束域的样本协方差矩阵：即由波束域的样本协方差矩阵代替s22中得到的阵元域样本协方差矩阵矩阵维数由l×l变为(1+p)×(1+p)，其中表示波束域第l个子阵列回波数据；s52：利用对角加载技术增加算法稳定性：其中i为单位矩阵，rb(k)为对角加载后的波束域协方差矩阵，λ为对角加载系数，且满足δ为常数，且满足在本发明中s6：通过获得的波束域非自适应权矢量，阻塞矩阵以及波束域的协方差矩阵，得到波束域广义旁瓣相消算法权矢量的表达式；具体步骤如下：s61：传统广义旁瓣相消算法权矢量的表达式：w＝wq-bwas.tbhwq＝0式中，wq是n×1的固定权矢量，b为n×(n-1)的阻塞矩阵，由于wq为固定的向量，且b与向量wq正交，所以广义旁瓣相消将最小方差算法中的约束优化问题转换为非约束问题，按下式计算得到自适应权矢量的最优解wa：wa＝(bhrb)-1bhrwqs62：波束域广义旁瓣相消算法权矢量wb按下式计算：式中，wbq为波束域非自适应权矢量，bb为波束域阻塞矩阵，wba为波束域自适应权矢量；其中自适应权矢量s7：利用波束域广义旁瓣相消算法权矢量对波束域回波信号进行加权求和，得到波束域广义旁瓣相消算法的输出信号。自适应输出波束信号按下式计算：式中，为权矢量wb的共轭转置，xb(k)为子阵列波束域回波信号空间平滑后的波束域信号，yb(k)即为波束域的输出波束信号。本发明中子阵阵元数目l的取值为l≤n/2。为了验证本发明的有效性，在本实施例中，利用fieldii对医学成像中常用的点散射目标和斑散射目标进行成像并采用michigan大学生物医学实验室提供的geabr_0数据进行实验成像。fieldii是基于线性系统空间响应原理，它的仿真结果与实际成像很接近，已被国际上广泛认同为仿真超声系统的标准。在点目标仿真实验中一共设了15个点目标，轴向位置分别在40mm～90mm之间，每隔5mm设立目标点，其中50mm和70mm处分别有三个目标点，侧向位置分别在-2mm、0mm和2mm处，仿真采用定点发射，动态接收聚焦的模式，发射点设在50mm处，接收聚焦点分别位于30mm～60mm之间，成像的动态范围为60db。同时设定中心在25mm处，半径为3mm的圆形区域为吸声斑区域，成像的动态范围也为60db。分别采用动态发射定点接收方式对4种算法进行成像，并比较各种成像算法的分辨率、对比度和算法的运行时间。geabr_0数据实验参数如下：阵元中心频率为3.33mhz，阵元数目为64个，间距为0.2413mm，采样频率为17.76mhz，声速为1500m/s，成像的动态范围设为60db。图4给出了动态发射定点接收下4种算法点目标仿真成像结果。图5给出了轴向距离为50mm处4种算法点目标侧向分辨率。从图4直观来看，传统的das算法的成像效果最差，侧向伪影较为严重，无法清晰的观察出每个独立点目标；esbgsc的成像效果最好，近场区域几乎没有侧向伪影，旁瓣等级最低；而gsc和本文提出的bgsc算法成像效果基本相同，介于esbgsc和das算法之间，这可以从直观上表明在点目标仿真中bgsc在成像效果上相比gsc并没有下降太多。从图5可以看出，das算法的成像分辨率最差，主瓣宽度均大于其它3种算法；esbgsc的成像分辨率明显优于其它算法，而gsc和bgsc的成像分辨率优于传统das算法，但比esbgsc要差，且它们两种算法的分辨率差别不大，分辨率曲线几乎重合。表1给出了4种算法点目标仿真的运行时间，本发明仿真在matlabr2015a下运行，使用的计算机配置如下：inteli74ghzcpu，8gbram，由表1我们可以看出传统das算法由于几乎不涉及矩阵和向量的运算，所以运行时间要比其他3种算法少得多，本发明采用的bgsc算法在相同的情况下，相比gsc和esbgsc算法由于大大降低了矩阵的维数，所以运行时间大大减少。表14种算法点目标成像时间对比表成像算法运行时间/sdas2.66gsc122.04esbgsc226.33bgsc93.23图6给出了动态发射定点接收情况下加白噪声4种算法点目标成像结果，图7给出了轴向距离为50mm加噪声情况下点目标侧向分辨率。从图6的结果可以看出加了信噪比为20db高斯白噪声后，背景噪声明显增强，可以看到较为明显的噪声白斑。从图7中可以看出加了噪声后4种算法的侧向分辨率相比之前不加噪声的情况并没有明显的变化，由此可见，4种算法均具有较强的鲁棒性，在有噪声的情况下也可以保持成像效果，侧向分辨率也没有明显的下降。图8给出了动态发射定点接收下4种成像算法吸声斑仿真结果，表2给出了4种算法吸声斑仿真成像的对比度等参数。从图8可以看出，在动态发射定点接收模式下，与点目标仿真类似，传统das成像效果最差，吸声斑与背景噪声区分度不大，无法较为清晰的看出圆形吸声斑的轮廓；esbgsc的成像效果最佳，吸声斑与背景区分度较大；而gsc与bgsc算法成像效果介于两者之间，只在吸声斑边界区域有较为明显的噪声斑点，但仍能清晰的看出吸声斑区域。对于吸声斑成像的参数指标其中对比度(cr)的定义为外部平均功率与中心暗斑平均功率之差，而背景区域方差是反映算法鲁棒性的参数，一般来说其值越小说明算法的鲁棒性越好。由表2可以看出，传统das算法的背景区域方差最小，说明该算法的鲁棒性最好，bgsc算法次之。而对比度方面，esbgsc相比das、gsc和bgsc分别提高了约8.1db、7.4db、7.4db和7.0db，然而esbgsc算法的稳健性却是最差的。本发明采用的bgsc算法对比度相比das和gsc算法均有所提高，稳健性也高于gsc和esbgsc算法。表24种算法成像指标对比表为了说明算法的运算时间，表3中列出了4种算法吸声斑成像的运行时间。与点目标成像结果类似，das运行时间最短，esbgsc运行时间最长，而本发明采用的bgsc算法相比gsc和esbgsc在相同的情况下分别缩短了10.64s和47.61s。表34种算法吸声斑成像运行时间对比表成像算法运行时间/sdas1.10gsc48.50esbgsc85.47bgsc37.86图9给出了4种算法geabr_0数据实验成像结果，图10给出了4种成像算法在轴向距离为12mm处分辨率对比图。由图9可以看出，esbgsc算法成像效果最好，背景噪声亮度较暗，点目标相对于背景噪声较为明显，对比度较大；而传统das算法成像效果最差，背景噪声亮度较亮，分辨率与对比度均不及其它3种算法，目标点也不够清晰；gsc和本发明采用的bgsc算法成像效果介于das与esbgsc算法之间，可以较为明显的观察出点目标和吸声斑区域。从图10中可以看出，传统das算法不仅对比度不及其它3种算法，分辨率也是最差的，无论从主瓣宽度和旁瓣峰值方面均不及其它算法；esbgsc算法由于加了特征值分解提取出了信号子空间，因此分辨率最优，旁瓣等级也明显小于其它算法；而gsc和本发明采用的bgsc算法成像分辨率基本一致，优于das算法，不及esbgsc算法，其中主瓣宽度与esbgsc算法基本一致，但旁瓣峰值高于esbgsc算法。由此我们可以得出，在实验的条件下，bgsc算法依然可以保证成像效果较于gsc基本不变，但运行所需的时间却有明显的降低，具有一定的实用性。最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。当前第1页12