一种基于切比雪夫多项式的复合窗变迹超声波束形成方法
阅读说明：本技术 一种基于切比雪夫多项式的复合窗变迹超声波束形成方法 (Composite window apodization ultrasonic beam forming method based on Chebyshev polynomial ) 是由 王平 梁家祺 王慧悦 何理 柳学功 周瑜 李锡涛 何峰宇 武超 阎鑫龙 陈靖翰 于 2021-08-18 设计创作，主要内容包括：本发明涉及一种基于切比雪夫多项式的复合窗变迹超声波束形成方法,属于超声成像技术领域。该方法首先对超声阵列接收到的回波信号进行处理,得到所需要的超声阵列回波信号；然后输入不同切比雪夫参数分别构造两个切比雪夫窗函数；再将两个切比雪夫窗函数按一定权重进行叠加,得到复合变迹窗,并用复合窗对超声阵列回波信号进行加权；根据近场和远场各自的成像效果分别调整近场窗和远场窗的切比雪夫参数,直到符合预期要求,然后进行后续波束形成和最终成像；该方法能够根据探测深度自行调节谱加权的旁瓣抑制能力,能在近场有效抑制伪像,并改善远场分辨率,从而提高超声成像的整体质量。(The invention relates to a composite window apodization ultrasonic beam forming method based on Chebyshev polynomial, belonging to the technical field of ultrasonic imaging. Firstly, processing an echo signal received by an ultrasonic array to obtain a required ultrasonic array echo signal; then inputting different Chebyshev parameters to respectively construct two Chebyshev window functions; superposing the two Chebyshev window functions according to a certain weight to obtain a composite apodization window, and weighting the ultrasonic array echo signals by using the composite window; respectively adjusting Chebyshev parameters of a near field window and a far field window according to respective imaging effects of the near field and the far field until the Chebyshev parameters meet expected requirements, and then performing subsequent beam forming and final imaging; the method can automatically adjust the side lobe suppression capability of the spectral weighting according to the detection depth, can effectively suppress the artifact in the near field, and improves the far field resolution, thereby improving the overall quality of the ultrasonic imaging.)
技术领域
本发明属于超声成像技术领域，涉及一种基于切比雪夫多项式的复合窗变迹超声波束形成方法。
背景技术
超声成像是一种广泛应用在医学和工业领域的无损检测手段，而延时叠加波束形成器因其结构简单、易于实现的优势，被广泛应用在市面上的超声成像仪器当中。但是传统延时叠加方法得到的波束存在旁瓣水平较高的缺陷，旁瓣是产生伪像的主要原因，较高的旁瓣水平降低了超声成像质量。
抑制旁瓣的技术途径是采用幅度变迹，其实现方法是对接收阵元采用谱加权技术，谱加权技术利用傅里叶变换的性质来改变旁瓣。常见的传统幅度变迹采用的窗函数有三角窗、汉宁窗、海明窗等，都能有效抑制旁瓣高度，防止伪像的产生。然而幅度变迹的引入会带来新的问题，即主波束宽度的增大，尤其在远场成像时，主波束宽度增大所造成的横向分辨率的下降极大地降低了最终成像质量。同时由于超声衰减的存在，远场成像对变迹窗伪像抑制能力的要求没有近场成像这么高，横向分辨率才是影响远场成像质量的主要因素。因此，有学者提出将探测区域按深度分段，不同分段采用不同窗函数的变迹法来改善远场分辨率。但是传统的窗函数难以灵活改变自身旁瓣抑制能力，与分段变迹法无法较好适配，只能将现有的窗函数填入分段中测试成像效果，再从中选择一个效果相对较好的作为该段的窗函数。目前还缺少对于不同分段窗函数的构造与选择的系统性方法。
切比雪夫窗是由多尔夫提出的一种原本应用于雷达阵列定位的变迹窗，它的构造过程基于切比雪夫多项式的性质，且被证明该窗口能得到给定旁瓣高度下最小的主瓣宽度。但是和其他窗函数一样，切比雪夫窗仍存在主波束宽度和旁瓣高度之间的矛盾，即低旁瓣高度会伴随着较大的主瓣宽度。将切比雪夫窗引入到超声成像后，由于超声聚焦环节的存在，实际得到的回波响应和理论上的波束方向图存在偏差，需要通过规范的方法进行调整修正。
综上所述，急需发明一种能够在近场有效抑制伪像的同时，远场又能保持较高横向分辨率的高效波束形成算法，以提高整体超声成像质量。
发明内容
有鉴于此，本发明的目的在于提供一种基于切比雪夫多项式的复合窗变迹超声波束形成方法。该方法能够根据探测深度自行调节谱加权的旁瓣抑制能力，能在近场有效抑制伪像，在远场改善分辨率，从而提高超声成像的整体质量。
为达到上述目的，本发明提供如下技术方案：
一种基于切比雪夫多项式的复合窗变迹超声波束形成方法，该方法包括以下步骤：
S1：对超声阵元接收的回波信号进行放大、AD转换和延时聚焦处理，以获得超声阵列回波数据x(k)，x(k)表示为x(k)＝[x1(k),x2(k),…,xN(k)]，其中N表示超声阵列的阵元个数，k表示为对应采样深度d的采样时刻；
S2：根据切比雪夫多项式性质，输入不同切比雪夫参数AS-near和AS-far，分别构造两个切比雪夫窗函数wnear，wfar；
S3：将两个切比雪夫窗按一定权重进行叠加，叠加权值随探测深度变化，得到复合变迹窗wcomp(d)，其中d表示探测深度，并用复合窗对x(k)变迹后，进行波束形成与成像；
S4：根据近场成像效果调整近场窗wnear的切比雪夫参数直到符合预期；
S5：根据远场成像效果调整远场窗wfar的切比雪夫参数直到符合预期；
S6：完成复合窗调整，进行波束形成与最终成像。
可选的，在所述S2中，根据切比雪夫多项式性质，输入不同切比雪夫参数AS-near和AS-far，分别构造两个切比雪夫窗函数wnear，wfar，具体包括以下步骤：
S21：输入切比雪夫参数AS，单位为dB；定义主波束最大值和旁瓣水平的比值为R，则有：
S22：构造(N-1)阶切比雪夫多项式TN-1(α)，其表达式为：
其中，N为超声探头接收阵元个数；
S23：根据多项式性质，主波束最大值和旁瓣水平的比值等于TN-1(α0)的值，其中α0>1，解得α0的值，并对其进行尺度变换，得到ω：
则波束方向图B(ψ)在ψ空间表示为：
其中，因子1/R用于归一化波束方向图，ψ为空间变量，使得B(0)＝1；
S24：找到波束方向图的原零点ψpo，其位置为：
其在α空间表示为：
变换尺度到ω空间，得到：
最终得到在ψ空间的零点ψp：
S25：构造一个N×N的阵列流形矩阵V(ψ)其表达式为：
V(ψ)＝［ν(0) ν(ψ1) ... ν(ψp) ... ν(ψN-1)]
其中，ν(ψ)表示阵列流形向量，表达式为：
则最终切比雪夫参数AS对应切比雪夫窗的权值为：
w＝[VH(ψ)]-1·[1 0 ... 0]T
S26：按照上述步骤输入不同切比雪夫参数AS-near和AS-far，分别构造两个切比雪夫窗函数wnear，wfar，在参数选择上需要满足AS-near<AS-far，具体数值将在后面步骤进行调整。
可选的，在所述S3中，将两个切比雪夫窗按一定权重进行叠加，叠加权值随探测深度变化，得到复合变迹窗wcomp(d)，并用复合窗对x(k)变迹后进行波束形成，具体包括以下步骤：
S31：将两个切比雪夫窗按一定比例进行复合，得到复合变迹窗wcomp(d)，其表达式为：
其中，dmax和dmin分别表示探测区域的最大和最小深度；
S32：用构造出的复合窗进行对x(k)变迹后，进行波束形成，其波束形成信号表达式为：
其中，SDAS表示波束形成信号，xi(k)表示第i个阵元上对应深度的回波数据，wcomp,i表示第i个阵元上的变迹权值。
可选的，在所述S4中，根据近场成像效果调整近场窗wnear的切比雪夫参数直到符合预期要求，具体包括以下步骤：
S41：在初步成像中选取若干近场目标点，得到其在各自深度下的旁瓣峰值VPSL，计算其与动态成像范围左边界值VDRL的差值Vn,error，并设定一个可接受的阈值Vnear；
S42：若Vn,error大于0，则应减小近场切比雪夫参数AS-near，若差值小于阈值Vnear，则应增大AS-near，改变参数值后重复步骤S3和S4，直到满足条件Vnear<Vn,error<0，则当前的wnear,adj即为最终的近场切比雪夫窗函数。
可选的，在所述S5中，根据远场成像效果调整远场窗wfar的切比雪夫参数直到符合预期，具体包括以下步骤：
S51：在上次成像的最远场选取目标点，得到其幅值AM-far，设定一个最远场幅值的期望值Vfar，该期望值Vfar应略大于动态成像范围左边界值VDRL，满足VDRL<Vfar<1.05VDRL；
S52：若目标点幅值AM-far>Vfar，则AS-far增大10％；若AM-far<Vfar，则AS-far减小10％，重复步骤S3和S5，直到目标幅值AM-far和期望值Vfar的满足0.95|Vfar|<|AM-far-Vfar|<1.05|Vfar|，则当前的wfar,adj即为最终的近场切比雪夫窗函数。
可选的，在所述S6中，生成复合窗，对超声阵列回波信号x(k)变迹后，进行波束形成和最终成像，具体包括以下步骤：
S61：将调整之后的远近场切比雪夫窗重新复合得到wcomp,adj(d)，其表达式为：
S62：用调整后的复合窗进行对超声阵列回波信号x(k)变迹后，进行波束形成与最终成像，其波束形成信号表达式为：
本发明的有益效果在于：本发明将两个分别适用于近场成像的切比雪夫窗wnear和适用于远场成像的wfar按比例进行复合，使得整体复合窗随探测深度从wnear逐渐过渡到wfar，来满足不同深度下对窗函数性能需求侧重点的变化。相比于传统的变迹窗函数如汉宁窗、海宁窗等，本发明的复合窗能在近场有效抑制伪像的同时，在远场能更好的提高横向分辨率，从而提升整体的超声成像质量。
本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。
附图说明
为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：
图1为本发明所述的基于切比雪夫多项式的复合窗变迹超声成像方法的流程图；
图2为不同窗函数变迹后点目标仿真成像结果的比较图；
图3为不同窗函数变迹后点目标仿真结果100mm和200mm深度处切片图像；
图4为不同窗函数变迹后geabr_0数据成像结果比较图。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。
其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。
图1为本发明的算法流程图，本发明提供一种基于切比雪夫多项式的复合窗变迹超声波束形成方法，该方法包括以下步骤：
步骤S1：对超声阵元接收的回波信号进行放大、AD转换和延时聚焦处理，以获得超声回波数据x(k)，x(k)表示为x(k)＝[x1(k),x2(k),…,xN(k)]，其中N表示超声阵列的阵元个数，k表示为对应采样深度d的采样时刻。
步骤S2：根据切比雪夫多项式性质，输入不同切比雪夫参数AS-near和AS-far分别构造两个切比雪夫窗函数wnear，wfar，具体包括以下步骤：
S21：输入切比雪夫参数AS，单位为dB。定义主波束最大值和旁瓣水平的比值为R，则有：
S22：构造(N-1)阶切比雪夫多项式TN-1(α)，其表达式为：
其中，N为超声探头接收阵元个数；
S23：根据多项式性质，主波束最大值和旁瓣水平的比值等于TN-1(α0)的值，其中α0>1，解得α0的值，并对其进行尺度变换，得到ω：
则波束方向图B(ψ)在ψ空间可表示为：
其中，因子1/R用于归一化波束方向图，ψ为空间变量，使得B(0)＝1；
S24：找到波束方向图的原零点ψpo，其位置为：
其在α空间表示为：
变换尺度到ω空间，得到：
最终得到在ψ空间的零点ψp：
S25：构造一个N×N的阵列流形矩阵V(ψ)其表达式为：
V(ψ)＝[ν(0) ν(ψ1) ... ν(ψp) ... ν(ψN-1)]
其中，ν(ψ)表示阵列流形向量，表达式为：
则最终切比雪夫参数AS对应切比雪夫窗的权值为：
w＝[VH(ψ)]-1·[1 0 ... 0]T
S26：按照上述步骤输入不同切比雪夫参数AS-near和AS-far，分别构造两个切比雪夫窗函数wnear，wfar，在参数选择上需要满足AS-near<AS-far，具体数值将在后面步骤进行调整。
步骤S3：将两个切比雪夫窗按一定权重进行叠加，叠加权值随探测深度变化，得到复合变迹窗wcomp(d)，其中d表示探测深度，并用复合窗对x(k)变迹后进行波束形成与成像，具体包括以下步骤：
S31：将两个切比雪夫窗按一定比例进行复合，得到复合变迹窗wcomp(d)，其表达式为：
其中，dmax和dmin分别表示探测区域的最大和最小深度；
S32：用构造出的复合窗进行对x(k)变迹后，进行波束形成，其波束形成信号表达式为：
其中，SDAS表示波束形成信号，xi(k)表示第i个阵元上对应深度的回波数据，wcomp,i表示第i个阵元上的变迹权值。
步骤S4：根据近场成像效果调整近场窗wnear的切比雪夫参数直到符合预期要求，具体包括以下步骤：
S41：在初步成像中选取若干近场目标点，得到其在各自深度下的旁瓣峰值VPSL，计算其与动态成像范围左边界值VDRL的差值Vn,error，并设定一个可接受的阈值Vnear；
S42：若Vn,error大于0，则应减小近场切比雪夫参数AS-near，若差值小于阈值Vnear，则应增大AS-near，改变参数值后重复步骤S3和S4，直到满足条件Vnear<Vn,error<0，则当前的wnear,adj即为最终的近场切比雪夫窗函数。
步骤S5：根据远场成像效果调整远场窗wfar的切比雪夫参数直到符合预期要求，具体包括以下步骤：
S51：在上次成像的最远场选取目标点，得到其幅值AM-far，设定一个最远场幅值的期望值Vfar，该期望值Vfar应略大于动态成像范围左边界值VDRL，满足VDRL<Vfar<1.05VDRL；
S52：若目标点幅值AM-far>Vfar，则AS-far增大10％；若AM-far<Vfar，则AS-far减小10％，重复步骤S3和S5，直到目标幅值AM-far和期望值Vfar的满足0.95|Vfar|<|AM-far-Vfar|<1.05|Vfar|，则当前的wfar,adj即为最终的近场切比雪夫窗函数。
步骤S6：完成复合窗调整，进行波束形成与最终成像，具体包括以下步骤：
S61：将调整之后的远近场切比雪夫窗重新复合得到wcomp,adj(d)，其表达式为：
S62：用调整后的复合窗进行对超声阵列回波信号x(k)变迹后，进行波束形成与最终成像，其波束形成信号表达式为：
Field II是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证所提算法的有效性，利用Field II对超声成像中常用的点目标进行成像并利用geabr_0实验数据进行成像对比实验。在点目标仿真实验中，设置两列横向间隔为10mm，纵向间隔为10mm的40个点目标，均匀分布在30mm～200mm的深度之间，采用发射定点聚焦和接收动态聚焦方式，发射焦点固定在60mm处，设置图像的成像动态范围为65dB，同时设置了超声衰减参数为50dB/m。geabr_0实验所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为50dB。对上述两个实验目标采用矩形窗、汉宁窗、海明窗和本发明提出的复合窗进行变迹成像实验，在点目标成像仿真实验中采用近场-30dB远场-20dB的切比雪夫窗叠加而成的复合窗，在geabr_0实验中采用近场-15dB远场-3dB的切比雪夫窗叠加而成的复合窗。用双点目标的中心幅值(Middle Amplitude,MA)和半高宽(Full-Width at Half-Maximum,FWHM)来评价图像的横向分辨率，判断不同变迹窗的优劣。
图2给出了点目标仿真成像结果的比较图。由图2可知，经过矩形窗变迹的图像在近场存在明显伪像，随着探测深度的增加，超声衰减的影响逐渐加大，伪像幅值逐渐衰减出动态范围从图像上消失；汉宁窗变迹和海宁窗变迹的图像全程不存在明显伪像，但是其在远场的横向分辨率较低，导致同一深度的点目标的像相互混叠在一起无法区分；而复合窗变迹的图像和矩形窗结果相比在近场有更好的抑制伪像的能力，和汉宁窗、海宁窗结果相比在远场有更高的横向分辨率。表1为不同深度下四种窗函数变迹结果的中心幅值和半高宽。由表1可知，复合窗的FWHM平均值为2.404mm，相较于汉宁窗的3.280mm要窄26.6％，相较于海明窗的3.126mm要窄22.4％，在同深度下，复合窗的MA值要小于其他三种传统变迹窗，这意味着复合窗更有利于区分同一深度的两个点目标。
表1不同深度下四种变迹窗函数结果的中心幅值MA(dB)和半高宽FWHM(mm)
图3给出了点目标仿真成像结果中100mm深度和200mm深度处切片图。由图3可知，在100mm深度处，矩形窗变迹结果的旁瓣高度在动态成像范围之内，且与主瓣峰值之间存在谷值，故能形成较明显伪像，且两点中间的旁瓣相互叠加，形成了一个更加明显地伪像，而汉宁窗、海明窗和复合窗变迹结果的旁瓣峰值不在动态成像范围之内，故不会产生明显伪像，且和汉宁窗、海明窗相比，复合窗有更窄的主瓣宽度和更低的中心幅值；在200mm深度处，由于超声衰减的存在，波束幅值减小，矩形窗结果的旁瓣高度衰减处成像范围，故伪像消失，而汉宁窗、海宁窗的结果的主瓣宽度过宽，使其在两点中间混叠出了一个新的峰值，严重降低了成像质量，而复合窗的结果仍保持了相对较窄的主瓣宽度，没有因主瓣混叠产生新的伪像。综上所述，在点目标仿真成像实验中，本发明所提出的复合窗能显著改善最终成像质量。
图4给出了geabr_0成像实验结果比较图。取纵向排列的点目标作为观察对象，将其从近到远进行编号，则部分点目标的半高宽如表2所示。由图4和表2可知，相比与矩形窗、汉宁窗、海明窗的变迹结果，本发明所提出的复合窗横向分辨率更高，点目标横向扩散程度更小，其优势在远场尤为明显。
表2不同窗函数下部分纵向点目标变迹结果的半高宽(mm)
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。