技术领域
本发明涉及医学超声成像领域，特别涉及一种基于自适应加权的双聚焦波束 合成方法。
背景技术
在医学超声成像中，非信源方向上的回波信号引入了广泛的图像杂波，导致 了图像质量的下降。目前，一般采用传统延时叠加方法减弱非信源方向上信号的 干扰，提高信源方向上信号的强度。该方法对不同通道接收到的数据施以不同的 延时，再将其相加，得到期望聚焦点的幅值。单次延时叠加算法简单，易实现， 但其抗噪声能力弱，空间分辨率较低，成像对比度较差。
近年来，虚拟阵元的概念以及自适应波束合成技术已广泛应用于军事、民用 通信以及雷达天线探测系统等阵列信号处理领域。虚拟阵元技术的基本思想是： 在保持真实阵元数目不变的情况下，应用虚拟阵元技术使得阵元数目在虚拟上得 到增加，从而减小波束宽度，提高分辨率，增加探测深度。然而，目前在国内超 声成像领域，将虚拟阵元技术直接应用于超声成像，将使得超声成像在近场的伪 像特别严重，这是虚拟阵元应用于超声成像的一个固有缺点。自适应波束合成的 基本思想是在无失真约束条件下在指定方位或频率上获得最小方差，进而推导出 最优加权向量。由于自适应波束合成方法只适用于远场、窄带非相关信号，而超 声数据具有近场、宽带和强相关性等特点都限制了自适应波束合成算法在医学超 声成像中的应用。
随着临床诊断需求的不断增加，对超声成像系统的成像质量提出了更高的要 求。因此急需一种能够同时提高超声成像分辨率和增加探测深度的波束合成方 法。
发明内容
有鉴于此，本发明所要解决的技术问题是提供一种能够同时提高超声成像分 辨率和增加探测深度的聚焦波束合成方法。
本发明的目的是研究出一种高质量的基于自适应加权的双聚焦波束合成方 法，通过抑制旁瓣伪像，提高波束合成质量，从而全面提高超声图像的质量。
本发明的目的是通过以下技术方案来实现的：
本发明提供的一种基于自适应加权的双聚焦波束合成方法，包括以下步骤：
S1：利用单一固定焦点进行发射和接收超声波；
S2：利用Hanning窗加权函数对回波数据进行加权，第一次聚焦采用了定点 的延时叠加波束合成聚焦；
S3：利用虚拟阵元方法根据第一次聚焦形成的数据来计算第二次动态聚焦的 相关延时参数；
S4：利用最小方差自适应算法根据第一次聚焦形成的聚焦数据计算自适应加 权系数；
S5：根据聚焦延时参数以及相应的动态自适应加权系数进行第二次延时叠加 波束合成聚焦，得到最终成像的扫描线回波数据。
进一步，所述第一次波束合成采用单一固定发射接收焦点，采用Hanning窗 对回波信号进行加权后聚焦，形成的数据存储于缓存器FIFO中，第二次波束合 成采用第一次的聚焦结果，对FIFO中的数据利用自适应加权的方法进行动态波 束合成，形成超声图像的扫描线。
进一步，所述第一次聚焦包括以下步骤：
S21：通过以下公式来计算滑动子孔径中第i个虚拟阵元的位置xi为：
x i = ( i - L + 1 2 ) × d , i = 1,2 , . . . L - - - ( 1 ) ]]>
其中，d为虚拟阵元间距，Zv为虚拟阵元深度，L=Zv/F#/d为滑动子孔径线 性传感器个数，D=L×d为滑动子孔径大小，F#为聚焦系数；
S22：通过以下公式来计算虚拟阵元i的延时τi：
τ i = x i 2 + Z v 2 - Z v c - - - ( 2 ) ]]>
其中c为声速，以滑动子孔径中心作为延迟时间参考点，在虚拟阵元处进行 定点聚焦；
S23：通过以下公式来计算各个通道的回波数据形成扫描线数据：
S ( t ) = Σ i = 1 L w BF 1 ( i ) s i ( t - r c - τ i ) w BF 1 ( i ) = 0.5 [ 1 - cos ( 2 π i L ) ] , 0 ≤ i ≤ L - - - ( 3 ) ]]>
其中，S(t)为N-L+1条扫描线数据，wBF1(i)为Hanning窗加权函数，si(t) 为子孔径中阵元i接收到的回波信号，r/c为声波从像点到子孔径中心的传播时 间，τi为第i个虚拟阵元施加的延时。
进一步，所述第二次聚焦包括以下步骤：
S41：读取缓存器FIFO中的数据并通过以下公式判断数据样本点是否为有 效样本点：
其中，d_x为样本点到相应虚拟阵元的侧向距离，d_z为轴向距离，θa为虚 拟阵元声场的半孔径角，ii为数据样本点编号，jj为虚拟阵元位置点；若Kjj，ii=1， 则样本点位于该虚拟阵元的声场内，该样本点即为有效样本点：
S42：通过以下公式计算半孔径角θa：
θ a = arctan D 2 Z v = arctan 1 2 F # - - - ( 5 ) ]]>
其中，D=L×d为滑动子孔径的大小，Zv表示虚拟阵元的深度，F#=Zv/D为 聚焦系数；
S43：通过以下公式计算样本点ii到虚拟阵元jj的延时参数：
τ jj , ii = 2 Z v ± d _ x 2 + d _ z 2 c - - - ( 6 ) ]]>
其中，“±”表示样本点在相应虚拟阵元的下方或者上方，c为声速；
S44：利用最小方差自适应算法，计算加权系数：
对于M个等间距阵元形成的传感器阵列，在阵列的近场区域存在一些点散 射目标，波束合成的输出可表达为：
y ( t ) = w H ( t ) x d ( t ) = Σ i = 1 M w i ( t ) x i ( t - Δ i ) - - - ( 7 ) ]]>
其中，t为时间系数，w(t)=[w1(t),…,wM(t)]T为加权向量，wH(t)是w(t)的共轭 转置，xd(t)为延时后聚焦的信号，表示为xd(t)=[x1(t-Δ1),…,xM(t-ΔM)]T,Δi为各个通 道的延时量，xi(t-Δi)表示第i阵元通道延时后的信号。
根据第一次聚焦延时后得到的信号S(t)，所述信号S(t)表示为S(t)=[s1(t), s2(t),…,sN-L+1(t)]，其中，L表示每个子孔径传感器阵元的个数，N表示传感器阵 元的总个数，N-L+1表示滑动子孔径的个数。
寻找一个最优的加权向量w，在保持期望方向增益不变的条件下，使阵列的 输出能量最小：
min w w H R i + n w , ]]>约束条件wHa＝1 (8)
其中，公式（8）中的a为取全1方向向量，w为待求的最优加权向量，wH是w的共轭转置，Ri+n是第一次聚焦后输出扫描线数据S(t)的协方差矩阵，可由 公式（9）计算：
S ( t ) = [ s 1 ( t ) , s 2 ( t ) , . . . , s N - L + 1 ( t ) R i + n ( i , j ) = E { [ s i ( t ) - E ( s i ( t ) ) ] [ s j ( t ) - E ( s j ( t ) ) ] } - - - ( 9 ) ]]>
其中，公式(9)中的E(si(t))、E(sj(t))分别表示si(t)、sj(t)的期望值，Ri+n(i,j) 表示协方差矩阵Ri+n中第i行、第j列的元素；结合式（8）、（9），可以以下公式 （10）计算得到最优加权向量为：
w BF 2 = w opt = R i + n - 1 a a H R i + n - 1 a - - - ( 10 ) ]]>
其中，表示Ri+n的逆矩阵，aH表示a的共轭转置，wopt表示最优加权向 量，wBF2表示第2次聚焦时的动态加权系数；
S45：采用自适应加权系数wBF2=[wBF2(1),…,wBF2(N-L+1)]T的逐点接收聚焦 波束合成BF2的第n条扫描线数据为：
H n , ii ( t ) = Σ jj = 1 N - L + 1 w BF 2 ( jj ) K jj , ii S jj , ii ( t - τ jj , ii ) , ii = 1,2 , . . . , P - - - ( 11 ) ]]>
其中，wBF2(jj)为自适应加权系数，P为样本点总数，N为BF1所形成的扫描 线总数，Sjj,ii(t-τjj,ii)为波束合成器BF1所形成的第jj条扫描线上的第ii个样本点， τjj，ii为样本点ii相对于编号是jj的虚拟阵元的延时参数。
本发明的优点在于：本发明将幅度变迹技术、虚拟阵元技术、自适应波束合 成算法引入到超声成像系统，有机融合了传统的延时叠加(DAS)波束合成方法， 利用两次延时叠加实现双聚焦自适应加权波束合成的(dual focusing adaptive beamforming，DFAB)超声成像，能够同时提高超声成像的分辨率和增加探测深 度，有效降低旁瓣等级。该方法兼具虚拟阵元技术和自适应波束合成算法的共同 优点，克服了基于虚拟阵元技术的双聚焦在前后两阶段都采用固定加权系数进行 聚焦，从而导致生成图像在近场存在严重伪像的固有缺点，本发明在第一阶段采 用了固定的Hanning窗加权系数进行定点聚焦，而在第二阶段融合了最小方差自 适应算法，实时动态产生加权系数进行逐点聚焦，有效克服了基于虚拟阵元的双 聚焦波束合成方法在近场产生伪像的固有缺点，最大限度的抑制了旁瓣，并有效 解决了自适应波束合成算法应用于超声成像的关键技术难题，从而大幅度提高了 波束合成质量，使得成像效果获得显著提高。
本发明的其它优点、目标和特征在某种程度上将在随后的说明书中进行阐 述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而 易见的，或者可以从本发明的实践中得到教导。本发明的目标和其它优点可以通 过权利要求书，说明书，以及附图中所特别指出的结构来实现和获得。
附图说明
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明 作进一步的详细描述，其中：
图1为基于自适应加权的双聚焦波束合成DFAB原理框图；
图2为传统聚焦发射声场与虚拟阵元声场对比图，(a)传统发射聚焦，(b)一 个VE的声场，(c)2个VE的声场；
图3为BF1形成扫描线及延时参数的计算示意图，3(a)扫描线示意图，图3(b) 给出了形成第n条扫描线时延时参数τ的计算示意图；
图4为虚拟阵元延时参数计算及声场叠加示意图；
图5为动态孔径接收聚焦(DRF)方法的仿真图像；
图6为合成孔径聚焦(SF)方法的仿真图像；
图7为最小方差自适应波束合成(MV)方法的仿真图像；
图8为虚拟阵元双聚焦波束合成(DFB)方法的仿真图像；
图9为本发明提出的方法(DFAB)的仿真图像；
图10为本发明方法生成的超声图像在85mm处的横向截面图与现有方法的 对比。
具体实施方式
以下将结合附图，对本发明的优选实施例进行详细的描述；应当理解，优选 实施例仅为了说明本发明，而不是为了限制本发明的保护范围。
图1为基于自适应加权的双聚焦波束合成方法DFAB的原理框图；图2为传 统聚焦发射声场与虚拟阵元声场对比图，(a)传统发射聚焦，(b)一个虚拟阵元VE 的声场，(c)2个VE的声场；图3为BF1形成扫描线及延时参数计算示意图，3(a) 扫描线示意图，图3(b)给出了形成第n条扫描线时延时参数τ的计算示意图，如 图所示：本发明提供的一种基于自适应加权的双聚焦波束合成方法，包括以下步 骤：
S1：利用单一固定焦点进行发射和接收超声波；
S2：利用Hanning窗加权函数对回波数据进行加权，通过第一次定点延时叠 加波束合成聚焦，并将第一次波束合成形成的数据存储于缓存器FIFO中；所述 第一次聚焦包括以下步骤：
S21：通过以下公式（1）来计算滑动子孔径中第i个虚拟阵元的位置xi为：
x i = ( i - L + 1 2 ) × d , i = 1,2 , . . . L - - - ( 1 ) ]]>
其中，d为虚拟阵元间距，Zv为虚拟阵元深度，L=Zv/F#/d为滑动子孔径线 性传感器个数，D=L×d为滑动子孔径大小，F#为聚焦系数；
S22：通过以下公式（2）来计算虚拟阵元i的延时τi：
τ i = x i 2 + Z v 2 - Z v c - - - ( 2 ) ]]>
其中c为声速，以滑动子孔径中心作为延迟时间参考点，在虚拟阵元处进行 定点聚焦；
S23：通过以下公式（3）来计算各个通道的回波数据形成扫描线数据：
S ( t ) = Σ i = 1 L w BF 1 ( i ) s i ( t - r c - τ i ) w BF 1 ( i ) = 0.5 [ 1 - cos ( 2 π i L ) ] , 0 ≤ i ≤ L - - - ( 3 ) ]]>
其中，S(t)为N-L+1条扫描线数据，wBF1(i)为Hanning窗加权函数，si(t) 为子孔径中阵元i接收到的回波信号，r/c为声波从像点到子孔径中心的传播时 间，τi为第i个虚拟阵元施加的延时。
S3：利用虚拟阵元方法根据第一次聚焦形成的数据来计算第二次动态聚焦的 相关延时参数；
S4：根据聚焦延时参数进行第二次延时叠加波束合成聚焦；所述第二次聚焦 包括以下步骤：
S41：读取缓存器FIFO中的数据并通过以下公式判断数据样本点是否为有 效样本点：
其中，d_x为样本点到相应虚拟阵元的侧向距离，d_z为轴向距离，θa为虚 拟阵元声场的半孔径角，ii为数据样本点编号，jj为虚拟阵元位置点；若Kjj，ii=1， 则样本点位于该虚拟阵元的声场内，该样本点即为有效样本点：
S42：通过以下公式计算半孔径角θa：
θ a = arctan D 2 Z v = arctan 1 2 F # - - - ( 5 ) ]]>
其中，D=L×d为滑动子孔径的大小，Zv表示虚拟阵元的深度，F#=Zv/D为 聚焦系数；
S43：通过以下公式计算样本点ii到虚拟阵元jj的延时参数：
τ jj , ii = 2 Z v ± d _ x 2 + d _ z 2 c - - - ( 6 ) ]]>
其中，“±”表示样本点在相应虚拟阵元的下方或者上方，c为声速；
S44：利用最小方差自适应算法，计算加权系数：
对于M个等间距阵元形成的传感器阵列，在阵列的近场区域存在一些点散 射目标，波束合成的输出可表达为：
y ( t ) = w H ( t ) x d ( t ) = Σ i = 1 M w i ( t ) x i ( t - Δ i ) - - - ( 7 ) ]]>
其中，t为时间系数，w(t)=[w1(t),…,wM(t)]T为加权向量，wH(t)是w(t)的共轭转置， xd(t)为延时后聚焦的信号，表示为xd(t)=[x1(t-Δ1),…,xM(t-ΔM)]T,Δi为各个通道的延 时量，xi(t-Δi)表示第i阵元通道延时后的信号。
根据第一次聚焦延时后得到的信号S(t)，所述信号S(t)表示为S(t)=[s1(t), s2(t),…,sN-L+1(t)]，其中，L表示每个子孔径传感器阵元的个数，N表示传感器阵 元的总个数，N-L+1表示滑动子孔径的个数。
寻找一个最优的加权向量w，在保持期望方向增益不变的条件下，使阵列的 输出能量最小：
min w w H R i + n w , ]]>约束条件wHa＝1 (8)
其中，公式（8）中的a为取全1方向向量，w为待求的最优加权向量，wH是w的共轭转置，Ri+n是第一次聚焦后输出扫描线数据S(t)的协方差矩阵，可由 公式（9）计算：
S ( t ) = [ s 1 ( t ) , s 2 ( t ) , . . . , s N - L + 1 ( t ) R i + n ( i , j ) = E { [ s i ( t ) - E ( s i ( t ) ) ] [ s j ( t ) - E ( s j ( t ) ) ] } - - - ( 9 ) ]]>
其中，公式（9）中的E(si(t))、E(sj(t))分别表示si(t)、sj(t)的期望值，Ri+n(i,j) 表示协方差矩阵Ri+n中第i行、第j列的元素；结合式（8）、（9），可以以下公式
（10）计算得到最优加权向量为：
w BF 2 = w opt = R i + n - 1 a a H R i + n - 1 a - - - ( 10 ) ]]>
其中，表示Ri+n的逆矩阵，aH表示a的共轭转置，wopt表示最优加权向 量，wBF2表示第2次聚焦时的动态加权系数；
S45：采用自适应加权系数wBF2=[wBF2(1),…,wBF2(N-L+1)]T的逐点接收聚焦 波束合成BF2的第n条扫描线数据为：
H n , ii ( t ) = Σ jj = 1 N - L + 1 w BF 2 ( jj ) K jj , ii S jj , ii ( t - τ jj , ii ) , ii = 1,2 , . . . , P - - - ( 11 ) ]]>
其中，wBF2(jj)为自适应加权系数，P为样本点总数，N为BF1所形成的扫描 线总数，Sjj,ii(t-τjj,ii)为波束合成器BF1所形成的第jj条扫描线上的第ii个样本点， τjj，ii为样本点ii相对于编号是jj的虚拟阵元的延时参数。
S5：得到最终成像的扫描线回波数据，进行坐标变换、后续图像数据处理与 成像显示。
为了验证该算法的有效性，利用Field II对医学成像中常用的点散射目标和 斑散射目标进行成像。Field II是基于线性系统空间响应的原理，它的仿真效果 与实际成像很接近，已被国际上广泛认同为仿真超声系统的标准。点散射目标的 目标散射点共16个，分布在深度为60~95mm，宽度为16mm的区域内；仿真模 型的参数设置为：采用定点发射和分段动态聚焦接收模式，成像采用线性阵列， 阵元总数128，发射信号中心频率3.5MHz，采样频率50MHz，阵元中心间距为 一个波长，声速1540m/s。成像的动态范围为40dB。设置虚拟源深度Zv=20mm， 聚焦系数F#=2，其中，DFB算法采用固定加权(hanning窗)；DFAB算法采用自 适应加权。为提高成像效果，采用了滑动子孔径技术，子孔径阵元数目为48。
图5至图9给出了本发明算法生成的点目标成像图与现有算法生成的点目标 成像图的对比。图5为传统的动态孔径聚焦(DRF)方法；图6为合成孔径聚焦(SF) 方法；图7为最小方差自适应波束合成(MV)方法；图8为双聚焦波束合成(DFB) 方法；图9为本发明所提的方法(DFAB)。图10给出了本发明算法生成的点目标 成像图中在深度85mm处的横向截面图与现有的算法相比。从图5至图9可见： DRF算法横向分辨率和降低旁瓣等级效果均为最差；与DRF算法相比，SF算法 在近场区域虽然有较好的横向分辨率，但是随着探测深度增加，在远场区域横向 分辨率严重下降，且系统复杂度远高于DRF；而MV算法虽然在近场区域横向 分辨率比较均匀，但随着探测深度增加，成像宽度逐渐增加，分辨率逐渐下降， 且远场点目标成像不完整，出现缺失；为了在整个探测深度内得到均匀的横向分 辨率，采用了DFB算法，虽然DFB算法横向分辨率显著提高，但是旁瓣等级太 高，出现了严重的伪像；本发明算法DFAB兼具虚拟阵元技术和自适应波束合成 算法的共同优点，无论是横向分辨率还是旁瓣等级，均得到了很好地改善。因此， 在综合权衡图像的横向分辨率，幅度分辨率以及降低旁瓣等级的能力，DFAB算 法效果最好。
表1给出了各种成像方法在85mm处成像的横向宽度和旁瓣幅度的具体数 据，可以明确表明本发明算法DFAB在提高横向分辨率以及对旁瓣幅度的压制方 面的优越性。因此，在保证相同分辨率的情况下，DFAB算法间接地增加了探测 深度。
表185mm处散射点-40dB时的平均横向宽度及归一化旁瓣幅度
以上所述仅为本发明的优选实施例，并不用于限制本发明，显然，本领域的 技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这 样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之 内，则本发明也意图包含这些改动和变型在内。