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

在冶金、化工、电力等工业过程中经常出现两相流或多相流介质，例如在冶金行业的选矿厂、冶炼厂和矿山充填工艺生产过程中，常常要求对浆状的混合液中的固体(即金属)量进行测量，以利于成本核算、质量控制和金属回收。

作为一种用以测量多相流流动的过程参数的技术，超声过程层析成像技术(Ultrasonic Process Tomography，UPT)在 20 世纪 80 年代中后期逐渐形成和发展起来，并受到各界的广泛关注。层析成像的理论要求只有在投影角无限多，且每个投影角度下的投影线足够多的情况下才能获得对象的“精确重建”。

目前，扇形束扫描方式由于有效投影数据量大，已经得到了广泛的应用。其主要的实现方式是采用一种半柱面式的压电晶片，换能器的结构如图 1 所示，最大特点是其声辐射/接收的角度范围为180°，其声场可以覆盖检测区域。

但是，超声层析成像技术仍有如下方面的制约：超声波发射/接收同步和扫描切换技术不够成熟，实现超声波快速成像较为困难。此外，若超声波波长与障碍物尺寸相当, 会产生声散射、绕射等现象，且随机发生，给测量分析带来难度。

图1 换能器结构示意图

以本文的工作将利用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 的聚四氟乙烯椭圆柱障碍物。

障碍物聚四氟乙烯的声速 c0=5200 m/s，密度 ρ0= 7865 kg/m3。障碍物外部为介质水，声速 cs=1473 m/s，密度 ρs=997 kg/m3。假定某时刻，根据实际测量过程中采用的扇形束扫描模式，位于最左边的换能器为发射探头，剩下其余的换能器为接收探头。

由于换能器晶片的发射/接收面是半径为 5 mm 的半圆柱面，因此设定成像区域外边界 5 mm 位置，即(?0.06，0) m 处有一柱面超声波源，频率 f=2 MHz，则声波在障碍物区域的波长 λ0=c0/f=2.6 mm，在介质水的区域波长 λs=cs/f= 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 问题，该方法被证明是有效的。

1.3 网格划分

网格的划分对于模拟计算的收敛性和准确性有较大的影响。要给出可接受的准确性，声场模块的超弱变分公式需要一个相对均匀的网格，并通过最大单元尺寸参数 hmax 控制，而不是通过明显的几何细节。

COMSOL 软件的超声散射模型指出，当使用默认的 UWVF 设置时，定义的 hmax 等于或稍微大于 2 nλ(n 为正整数)时，会在收敛性和准确性之间给出较好的平衡，可以获得较好的仿真结果。由于本文建立的模型中介质水和障碍物区域中的波长并不一致，需要在 COMSOL 中对这两个不同区域设置并定义不同大小的 hmax。

根据上述原则，分别对介质水和障碍物两个求解域的 hmax设置为 0.0015 m 和 0.0052 m。网格划分结果如图 3 所示。

图3 模型的网格划分

1.4 求解过程

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

其中：R=[(x?x0)2+(y?y0)2]1/2 指成像区域内任意一点坐标(x，y)离声源位置(x0, y0)的距离；ks=2 π /λs，表示介质中的波数；函数 J 和 Y 分别指第一、二类Bessel 函数，i 为虚部单位。

成像区域的外边界施加辐射边界条件，吸收柱面波。此外，要取得可接受的准确性，声场模块的超弱变分公式采用默认20阶超弱 Helmholtz 单元。对求解域进行求解即可获得声场分布。

2、仿真结果分析

本文算例的模拟结果如图 4 所示，结果以声压级云图表示，显示了超声层析成像系统中液态介质和固态障碍物的声压级分布情况。

图4 声场云图

从图中可以看出，超声波在水中的扩散致使声压沿传播路径逐渐减弱。

另外，每一个柱体对超声波存在明显的遮挡作用，柱体后向声压非常弱，形成阴影区。

考虑到扇形束扫描的几何近似理论：只有当入射波的波长λ 远小于目标物的尺寸时，才能忽略由目标引起的衍射，即 ka >> 1(k=2π / λ，表示声波在介质中的波数；a 是障碍物的等效半径)。所以，满足该条件的 2 MHz 超声波所测障碍物的等效半径须远大于0.12 mm。

由于仿真几何模型中设定的障碍物尺寸均远大于该频率下所测障碍物等效半径的最小值，所以从该图中还可以看到，声波存在的散射现象，非常微弱，符合几何近似理论。

作者在实验中也注意到了类似现象，其表现为障碍物后方对应的传感器仍可以接收到极其微弱的超声信号，可以忽略不计。然后，提取和处理所获取的该声压分布结果，得到理论上阵列传感器的声压分布。

通过对代入二值逻辑反投影算法的信号矩阵进行阈值控制，得到一个信号矩阵 S(m，m)，其中 m 为所使用的换能器个数。那么，第 i 个发射位置对应的第 j 个接收位置的信号为：

3、图像重建算法

目前，采用二值逻辑反投影反演算法采用了上述的几何声学特性，例如基于三角形判断的反投影反演算法，其主要思想如图 5 所示。

当 1 号换能器发射时，4 号换能器被其传播路径中物体遮挡，则该探头的信号经过二值化处理后置 0。这样，由换能器间隔的中点 a、b、c 三点组成的区域表示目标物影响区间，图像矩阵对应区域增加权重“1”。

在此基础上，本文提出了一种基于四边形判断的二值逻辑反投影反演算法，如图 6 所示，当 1 号换能器发射时，4 号换能器被其传播路径中物体遮挡，则该探头的信号经过二值化处理后置 0。

这样，由换能器间隔的中点 a、b、c、d 四个点组成的四边形区域表示目标物影响区间，图像矩阵对应区域增加权重“1”，以此类推，待所有扫描信号处理完毕后通过阈值对图像进行目标和背景分割，获得重建图像。

分别按三角形判断准则和四边形判断准则对其进行反演，得到的结果如图 7 和图 8 所示。

将两种方式的重建结果与仿真对象进行对比发现，声波的散射和边缘效应没有对两种算法的图像重建带来较大影响；障碍物的位置信息、形状信息和尺寸信息，都较为准确。说明了二值逻辑反投影算法所用的几何声学理论在入射波的波长 λ 远小于目标物的尺寸时具备可行性。

图 7 和图 8 虽然可以看出两种判断准则的重建效果大致一样，但仔细对比上述两种算法，可以发现细微的差别。主要表现在：两种准则都能较好地得到圆柱障碍物的截面形状，但是四边形准则得到的圆柱障碍物的尺寸更接近真实模型，而三角形准则测得的偏大；此外，四边形准则得到的椭圆柱障碍物的截面信息也稍微接近真实模型。

4、结 论

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

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

来源：midiio