基于 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 等于或稍大于 2nλ(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 软件建立不同形状、尺寸的柱体障碍物的有限元模型，并利用超弱变分公式模拟仿真成像区域中声场的各点声压级，提取各个换能器位置的声压级数值大小可获得超声信号矩阵。

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