基于COMSOL的超声层析成像仿真及图像重建

顾建飞，苏明旭，蔡小舒



基于COMSOL的超声层析成像仿真及图像重建

顾建飞，苏明旭，蔡小舒

(上海理工大学颗粒与两相流测量研究所，上海市动力工程多相流动与传热重点实验室，上海 200093)

为了验证三角形和四边形准则的二值逻辑反投影(Binary Logic Back Projection，BLBP)算法在超声层析成像重建过程中的可行性，以超声层析成像的重建区域中同时含有圆柱形和椭圆柱形障碍物为例，利用COMSOL对成像区域内存在障碍物的声场进行仿真，得到声场声压分布。然后，通过提取和处理所获取的声压分布，获得信号矩阵。最后，将信号矩阵代入二值逻辑反投影算法进行图像重建。仿真结果表明：当入射波的波长λ远小于障碍物的尺寸时，声衍射的现象并不明显。重建结果表明：两种准则在几何近似理论下可行，且四边形准则的效果优于三角形准则。

超声波；有限元；图像重建；超声过程层析成像

0 引言

在冶金、化工、电力等工业过程中经常出现两相流或多相流介质，例如在冶金行业的选矿厂、冶炼厂和矿山充填工艺生产过程中，常常要求对浆状的混合液中的固体(即金属)量进行测量，以利于成本核算、质量控制和金属回收。作为一种用以测量多相流流动的过程参数的技术，超声过程层析成像技术(Ultrasonic Process Tomography，UPT)在20世纪80年代中后期逐渐形成和发展起来，并受到各界的广泛关注[1-3]。层析成像的理论要求只有在投影角无限多，且每个投影角度下的投影线足够多的情况下才能获得对象的“精确重建”。目前，扇形束扫描方式由于有效投影数据量大，已经得到了广泛的应用。其主要的实现方式是采用一种半柱面式的压电晶片，换能器的结构如图1所示，最大特点是其声辐射/接收的角度范围为180°，其声场可以覆盖检测区域。但是，超声层析成像技术仍有如下方面的制约：超声波发射/接收同步和扫描切换技术不够成熟，实现超声波快速成像较为困难。此外，若超声波波长与障碍物尺寸相当, 会产生声散射、绕射等现象，且随机发生，给测量分析带来难度。

COMSOL Multiphysics是一款大型数值仿真软件，通过有限元方法模拟，在科研和工程中能用偏微分方程(Partial Differential Equation，PDE)描述各种问题。它是集前处理、求解器和后处理于一体的有限元软件，提供友好的图形用户界面和建模工具，简化了建模工作。所以本文的工作将利用COMSOL软件模拟柱面压电换能器在超声层析成像技术中的声场。在层析成像的信号提取阶段，利用COMSOL多物理场耦合软件声学模块模拟声场，得到理论上阵列传感器的声压分布，这对UPT实验系统的搭建具有一定的指导意义，为信号矩阵的精确度提供参照作用。同时也可以对改进的二值逻辑反投影算法做进一步说明和验证。

1 模型建立

1.1 几何模型

本文采用如图2所示的二维几何模型，建立与实验系统相同的测量区域。其中，区域C0安装了16个换能器的阵列装置；区域C1是半径为55 mm的成像区域。区域C2是一个半径为10 mm的圆形区域，圆心位置在(-0.01，0.02) m，表示的是聚四氟乙烯圆柱障碍物。区域C3是一个逆时针旋转45°的长半轴为18 mm、短半轴为12 mm、中点在(0，-0.02) m的聚四氟乙烯椭圆柱障碍物。障碍物聚四氟乙烯的声速0=5200 m/s，密度0= 7865 kg/m3。障碍物外部为介质水，声速s=1473 m/s，密度s= 997 kg/m3。假定某时刻，根据实际测量过程中采用的扇形束扫描模式，位于最左边的换能器为发射探头，剩下其余的换能器为接收探头。由于换能器晶片的发射/接收面是半径为5 mm的半圆柱面，因此设定成像区域外边界5 mm位置，即(-0.06，0) m处有一柱面超声波源，频率=2 MHz，则声波在障碍物区域的波长0=0/=2.6 mm，在介质水的区域波长s=s/= 0.736 mm。这里，计算出的波长对后续设定该区域的网格单元尺寸尤为重要。

1.2 超弱变分公式(Ultra Weak Variational Formulation，UWVF)模块

在固体和流体中，线性时谐波场可以由Helmholtz和Navier方程来描述。但是，由于波的振荡特性难以对波动方程作数值逼近。在能接受的精确度范围内，对Helmholtz或Navier方程问题作近似求解，需要有一个相对密集的空间描述。对于标准的数值计算方法，比如低阶有限差分法(Finite Difference，FD)和有限元法(Finite Element，FE)，根据经验，每个波长取10个点是合适的。然而高频超声，由于声波波长非常短，因此要求将网格划分得更细，这将带来计算复杂性问题。针对高频声波，文中利用COMSOL独有的超弱变分公式可以解决上述问题，和低阶有限元分析相比，UWVF能显著减少计算时间，具有明显的优势。这种优势来自UWVF的单元属性，UWVF单元属性包含了求解自由空间Helmholtz波方程的信息。对于二维、三维Helmholtz问题以及二维Navier问题，该方法被证明是有效的[4]。

1.3 网格划分

网格的划分对于模拟计算的收敛性和准确性有较大的影响。要给出可接受的准确性，声场模块的超弱变分公式需要一个相对均匀的网格，并通过最大单元尺寸参数max控制，而不是通过明显的几何细节。COMSOL软件的超声散射模型指出，当使用默认的UWVF设置时，定义的max等于或稍微大于2λ(为正整数)时，会在收敛性和准确性之间给出较好的平衡，可以获得较好的仿真结果。由于本文建立的模型中介质水和障碍物区域中的波长并不一致，需要在COMSOL中对这两个不同区域设置并定义不同大小的max。根据上述原则，分别对介质水和障碍物两个求解域的max设置为0.0015 m和0.0052 m。网格划分结果如图3所示。

1.4 求解过程

模拟采用了连续波(continuous wave)形式，声源发出的柱面波与第二类零阶Hankel函数成正比：

其中：=[(-0)2+(-0)2]1/2指成像区域内任意一点坐标(，)离声源位置(0,0)的距离；s=2/s，表示介质中的波数；函数J和Y分别指第一、二类Bessel函数，i为虚部单位。成像区域的外边界施加辐射边界条件，吸收柱面波。此外，要取得可接受的准确性，声场模块的超弱变分公式采用默认20阶超弱Helmholtz单元。对求解域进行求解即可获得声场分布。

2 仿真结果分析

本文算例的模拟结果如图4所示，结果以声压级云图表示，显示了超声层析成像系统中液态介质和固态障碍物的声压级分布情况。从图中可以看出，超声波在水中的扩散致使声压沿传播路径逐渐减弱。另外，每一个柱体对超声波存在明显的遮挡作用，柱体后向声压非常弱，形成阴影区。考虑到扇形束扫描的几何近似理论：只有当入射波的波长远小于目标物的尺寸时，才能忽略由目标引起的衍射[5]，即>> 1(=2π /，表示声波在介质中的波数；是障碍物的等效半径)。所以，满足该条件的2 MHz超声波所测障碍物的等效半径须远大于0.12 mm。由于仿真几何模型中设定的障碍物尺寸均远大于该频率下所测障碍物等效半径的最小值，所以从该图中还可以看到，声波存在的散射现象，非常微弱，符合几何近似理论。作者在实验中也注意到了类似现象，其表现为障碍物后方对应的传感器仍可以接收到极其微弱的超声信号，可以忽略不计。然后，提取和处理所获取的该声压分布结果，得到理论上阵列传感器的声压分布。通过对代入二值逻辑反投影算法的信号矩阵进行阈值控制，得到一个信号矩阵(，)，其中为所使用的换能器个数。那么，第个发射位置对应的第个接收位置的信号为：

3 图像重建算法

目前，采用二值逻辑反投影反演算法采用了上述的几何声学特性，例如基于三角形判断的反投影反演算法[6]，其主要思想如图5所示。当1号换能器发射时，4号换能器被其传播路径中物体遮挡，则该探头的信号经过二值化处理后置0。这样，由换能器间隔的中点a、b、c三点组成的区域表示目标物影响区间，图像矩阵对应区域增加权重“1”。在此基础上，本文提出了一种基于四边形判断的二值逻辑反投影反演算法，如图6所示，当1号换能器发射时，4号换能器被其传播路径中物体遮挡，则该探头的信号经过二值化处理后置0。这样，由换能器间隔的中点a、b、c、d四个点组成的四边形区域表示目标物影响区间，图像矩阵对应区域增加权重“1”，以此类推，待所有扫描信号处理完毕后通过阈值对图像进行目标和背景分割，获得重建图像。分别按三角形判断准则和四边形判断准则对其进行反演，得到的结果如图7和图8所示。

图6 四边形判断法

将两种方式的重建结果与仿真对象进行对比发现，声波的散射和边缘效应没有对两种算法的图像重建带来较大影响；障碍物的位置信息、形状信息和尺寸信息，都较为准确。说明了二值逻辑反投影算法所用的几何声学理论在入射波的波长远小于目标物的尺寸时具备可行性。图7和图8虽然可以看出两种判断准则的重建效果大致一样，但仔细对比上述两种算法，可以发现细微的差别。主要表现在：两种准则都能较好地得到圆柱障碍物的截面形状，但是四边形准则得到的圆柱障碍物的尺寸更接近真实模型，而三角形准则测得的偏大；此外，四边形准则得到的椭圆柱障碍物的截面信息也稍微接近真实模型。

4 结论

在超声层析成像领域，采用COMSOL软件建立不同形状、尺寸的柱体障碍物的有限元模型，并利用超弱变分公式模拟仿真成像区域中声场的各点声压级，提取各个换能器位置的声压级数值大小可获得超声信号矩阵。

结合三角形判断和四边形判断准则的二值逻辑反投影算法的反演结果表明，两种判断准则的重建效果大致一样，同时也验证了当入射波的波长远小于目标物时该反演算法的可行性。

[1] Ibrahim S, Yunus M A M, Khairi M T M, et al. A review on ultrasonic process tomography system[J]. Jurnal Teknologi, 2014, 70(3): 1-5.

[2] Abbaszadeh J, Rahim H A, Rahim R A, et al. Ultrasonic tomography system: optimizing the frequency in a metal pipe conveyor[C]//Applied Mechanics & Materials, 2013, 284: 572-576.

[3] Hafiz M, Rahiman F, Ruzairi, et al. Ultrasonic tomography - image reconstruction algorithms[J]. International Journal of Innovative Computing Information & Control Ijicic, 2012, 8(1): 527-538.

[4] Huttunen T. The ultra weak variational formulation for ultrasound transmission problems[M]. Kuopion yliopisto, 2004.

[5] Xu L J, Xu L A. Ultrasound tomography system used for monitoring bubbly gas/liquid two-phase flow[J]. Ultrasonics Ferroelectrics & Frequency Control IEEE Transactions on, 1997, 44(1): 67-76.

[6] 陈孝震, 苏明旭, 蔡小舒. 液固两相系中的超声波层析成像实验研究[C]//中国工程热物理学会. 2012.

CHEN Xiaozheng, SU Mingxu, CAI Xiaosu. Liquid-solid two-phase system of ultrasonic tomographic imaging experiments[C]//Chinese Society of Engineering Thermophysics. 2012.

TheCOMSOL basedsound filed simulation and image reconstruction of ultrasonic tomography

GU Jian-fei, SU Ming-xu, CAI Xiao-shu

(Institute of Particle and Two-phase Flow Measurement, University of Shanghai for Science and Technology,Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering, Shanghai 200093, China)

In order to verify the feasibility of the Binary Logic Back Projection (BLBP) algorithm of triangular and quadrilateral criteria in ultrasonic tomography reconstruction, a typical UPT reconstruction area, which simultaneously contains a cylindrical obstacle and an elliptic cylindrical obstacle,is taken as an example. COMSOL is used to make simulation of sound field when the obstacles exist, and an acoustic pressure distribution is obtained. Then, by extracting and processing the obtained sound pressure distribution, a signal matrix is formed. Finally, the signal matrix is substituted into the binary logic back projection algorithm for image reconstruction. It could be clearly observed that the effect of diffraction is insignificant from the computation results whenthe wavelength of the incident wave is smaller than the size of the obstacle. The reconstruction results show that the two criteria are feasible under the geometric acoustic, and the reconstructed quality of the quadrilateral criterion is better than that of the triangular one.

ultrasound; finite element; image reconstruction; Ultrasonic Process Tomography(UPT)

TB551

A

1000-3630(2016)-03-0231-04

10.16300/j.cnki.1000-3630.2016.03.009

2015-08-31;

2015-11-10

国家自然科学基金资助项目(51176128, 51206113)

顾建飞(1990－), 男, 江苏海门人, 硕士研究生, 研究方向为超声过程层析成像。

苏明旭, E-mail: sumx@usst.edu.cn