乳腺超声断层成像系统的有限元仿真分析∗

张 玉 张国军 裴 毓 张文栋

(中北大学动态测试省部共建实验室 太原 030051)

0 引言

乳腺癌是全世界女性最常见的恶性肿瘤之一。近年来，乳腺癌的发病率呈上升、年轻化趋势，已成为当前社会的重大公共卫生问题，因此，乳腺癌的精确诊断和及时治疗尤为重要。目前，进行乳腺癌筛查的方法主要有X 线钼靶、X 射线计算机断层扫描成像(X-ray computed tomography, X-CT)、磁共振成像(Magnetic resonance imaging, MRI)、多普勒超声[1−2]。然而，这些方法大多具有辐射性或者给患者带来疼痛。近年来，超声计算机断层扫描(Ultrasound computed tomography, USCT)作为一种新型的成像技术[3]，在乳腺癌早期精确诊断方面具有良好的应用前景。超声CT 以超声波在人体不同组织内的声波传播速度差异和人体内的线性衰减系数为基础[4]，利用超声波对物体进行透射得到投影数据，对组织内部的参数进行重建，从而得到乳腺图像。

COMSOL Multiphysics 软件以有限元法为基础，通过求解偏微分方程来实现真实物理现象的仿真。本文利用COMSOL 软件进行乳腺超声断层成像仿真，采用一发多收的扫描方式，得到256个环形阵列换能器的声压分布，提取了接收信号，并利用等角扇束滤波反投影重建算法进行验证。

1 基于COMSOL的仿真分析

1.1 几何建模

通 过COMSOL Multiphysics 中“AC/DC 模块”、“声学模块”和“结构力学模块”接口创建，建立了如图1所示的二维几何模型。仿真测试在半径75 mm 的圆形区域内进行，内置为水，256个压电超声换能器以环形排列的方式均匀分布，每两个换能器间隔1.4◦，置于水中。在(−20 mm, 10 mm)处设置一个半径为12 mm 的圆形区域来代表软组织，在(23 mm,−30 mm)处设置一个半径10 mm 的圆形铁块障碍物来代表肿瘤。由于超声波在软组织和肿瘤中的传播速度接近，本次实验为了直观、清晰地显示超声波在乳腺组织中的声压分布，因此，用铁块来代表肿瘤。复合压电超声换能器采用铝层-粘接层-压电陶瓷(PZT-5H)-粘接层-铝层的结构，粘接层设置的目的是阻抗匹配，如图2所示。

图1 几何模型Fig.1 Geometric model

图2 复合式超声换能器结构Fig.2 Structure of composite ultrasonic transducer

仿真时，采用一个换能器发射、其余255个换能器接收，即一发多收的扫描方式进行，使得接收换能器覆盖面积达到最大，获得足够多的投影数据，以进行高质量的图像重建。发射换能器按顺时针方向变换，依次进行256 次仿真来达到360◦环形扫描。环形分布的换能器相较于线阵分布，数据采集速度更快，避免了因机械旋转而产生的干扰，实验数据更精确。发射换能器发射一个频率为3 MHz的正弦脉冲波，周期为5，幅值为40 V。材料属性见表1。

表1 主要材料属性Table 1 Main material property

1.2 边界条件

波动方程是用来描述声学、光学、电磁学的波动情况的一类重要的偏微分方程，用COMSOL 中“压力声学，瞬态”接口求解瞬态波动方程[5]：

其中，t表示时间，ρ是流体密度，c是声速，pt是总声压，qd是偶极源，Qm是单极源，pb是背景压力波。辐射边界选择平面波辐射，完美匹配层中的典型波速设为水中的声速即1500 m/s，声压级参考压力使用水的参考压力，即Pref,SPL=1 µPa。

固体域和水域的双向自动耦合在水域和换能器边界处实现，边界条件为

其中，n代表表面法线方向，utt是结构加速度，FA是作用在结构上的载荷(每单位面积的力)。

换能器的压电效应在PZT-5域中实现，通过求解线性方程，使得固体力学方程和静电方程相耦合，方程通过将应力和应变与电场和电位移相耦合来模拟压电效应。

1.3 网格划分

在有限元仿真中，需要对所建模型进行网格划分，COMSOL Multiphysics中提供交互式网格划分环境，每一个网格划分操作都将添加到网格划分序列中，当网格划分序列中所有操作完成以后，得到最终的网格。时间步长以及剖分网格的大小和形状很大程度上限制了数值计算的精度。在流体域中，每个波长使用5∼6 个网格单元，通常设置最大单元尺寸为波长/5，即网格尺寸=波长/5。因为此模型中水、软组织和铁块的波长不同，最大单元网格大小分别为0.1 mm、0.101 mm、0.345 mm。完成网格剖分之后网格顶点数为26738，三角形单元数为48102，边单元数为10242，顶点单元数为1796，整个模型网格划分和局部划分如图3和图4所示。

图3 模型网格剖分图Fig.3 Meshing of model

图4 局部网格剖分图Fig.4 Local meshing

2 计算结果与分析

在网格剖分完成之后，选择并行稀疏直接求解器(Multifrontal massive parallel sparse direct solver, MUMPS)进行求解计算，该求解器具有完整的线性系统矩阵，能够实现实数和复数的运算，将求解器中解的最大频率设置为fmax,sol=3 MHz。

确定时间步进的方法有向后差分公式、广义α。在瞬态计算中，通常使用广义α的方法。因为向后差分公式算法会产生散射，波形畸变随计算时间成正比；广义α可以有效规避这种干扰，它用的是前5个时间步长的解，并且可以对下个时间步长的解进行预测[6]。本次计算采用广义α的方法，手动控制时间步长，采用默认的时间步长即周期/60，这样的计算结果精确高，能达到最佳性能。

仿真完成后，对仿真结果进行后处理，得到在乳腺超声断层成像中，超声波在乳腺中的声压分布情况。图5是第一个换能器发射后，时间为27.8 µs、41.7 µs、69.4 µs、97.2 µs时的分布云图。软组织、铁块的声速和密度不同，其声阻抗系数也不同，声速、密度越大，声阻抗越大[7]，超声换能器接收到的信号越小。从图5中可以看出，超声波在水中逐渐衰减，当超声波遇到铁块时，能量有明显的衰减。

由声衍射理论[8]：

式(5)中，k、λ分别代表声波在介质中的波数和波长，a代表障碍物半径。声衍射的强弱与声波波长和障碍物尺寸相关，ka ≫1 时，声衍射现象弱；ka=1时，声衍射现象强。超声波频率为3 MHz 时，为了避免声衍射现象造成的干扰，障碍物尺寸应远大于0.08 mm，所以本模型的设定符合理论，所接收的信号也更加精确。

图5 声压变化云图Fig.5 Cloud chart of sound pressure change

在超声断层成像中，主要是利用渡越时间和幅值衰减进行图像重建，发射信号和接收信号最大幅值处的时间差为渡越时间，最大幅值之比为衰减的程度[9]。仿真结束后，得到每个接收换能器的终端电压和透射原始数据。本次实验利用幅值衰减进行图像重建，一次发射可以得到255组接收信号，记录下每个接收信号前5 个周期中的最大幅值，将幅值衰减取对数得到投影数据，经过256次发射之后，得到一组完整的接收信号矩阵。图6分别为第1 个、第64 个、第128 个和第192 个换能器发射时，正对面5个换能器接收的信号。超声波在不同介质中传播速度不同，在密度高的介质中传播速度快，分析看出，同一个换能器发射，不同位置接收的信号幅值不同，而由于障碍物对超声波的遮挡，使得换能器接收到的信号差别较大。

图6 接收信号Fig.6 Received signal

3 图像重建

等角扇束系统示意图如图7所示，射线源位置S0，探测器分布在圆弧上，各探测器单元性能一致，发射源中心线左右扇面张开的角度都相同，扇形束呈等角度分布，射线S0E的位置由(β,γ)决定，等角扇束重建公式[10]为

其中，重建图像点M(r,ϕ)，S0M长度为L，S0E的投影为pf(γ,β)，D表示S0到旋转中心的距离。

图7 等角扇束系统示意图Fig.7 Schematic diagram of isometric fan beam system

本次实验中256个超声换能器等角度环形分布于乳腺四周，计算得出接收换能器旋转角度增量为1.4◦，发射换能器与任一接收换能器射线间隔λ为0.7◦，符合等角扇束滤波反投影算法的重建规律，可以利用此算法进行图像重建。将256 次发射后得到的接收信号矩阵代入算法中进行重建，结果见图8。可以看出，重建图像精度较高，分析发现与待重建图像即实验设计的模型一致，能分辨出铁块的位置。

图8 等角扇束滤波反投影重建图像Fig.8 Reconstructed image by equal angle fan beam

4 结论

本文利用COMSOL Multiphysics 多物理场仿真软件，模拟了一个256 个换能器呈环形分布的乳腺癌检测环境。采用一个换能器发射、其余255 个换能器接收的方式进行仿真模拟，提取了256 组超声换能器接收信号组成信号矩阵。通过等角扇束滤波反投影算法进行图像重建，根据重建图像，可以较清晰地分辨出软组织和铁块障碍物的位置、形状和大小，经过比对，重建结果与仿真模型一致，且重建图像精度较高。验证了换能器呈360◦排列、进行环形扫描的方式可行。