基于有限元和k空间伪谱法的光声成像仿真平台

齐梦雨，严壮志, ， 赵丽丽

1 上海大学通信与信息工程学院，上海市，200444

2 上海大学上海生物医学工程研究所，上海市，200444

光声成像技术是一种新型无创在体生物医学成像技术。光声成像结合了超声和光学成像的优点[1]，兼具超声成像的高分辨率和纯光学成像的高对比度特性，因而成为了生物医学成像领域的研究前沿和热点之一[2]。光声成像的基本原理是光声效应[3]：使用短脉冲激光照射生物组织，组织体由于吸收能量导致的瞬时热膨胀进而辐射超声波。通过检测和采集这些携带生物组织特性的超声信号，就可以利用一定的重建算法重建出生物组织内部的光学吸收分布图像。由于生物组织的光吸收特性与组织的生理特征、结构形态特征、病变特性、代谢状态甚至神经活动密切相关，因此光声成像技术已经在肿瘤早期诊断[4-5]、生物分子成像[6]、生物组织功能成像[7]等生物医学领域中取得了重大的进展。

由于光声成像技术是一种多模态的混合成像技术，该技术的研究常常需要昂贵的实验设备，以及随之而来的繁杂的实验操作步骤。如果能将仿真技术应用在光声成像技术中，不仅可以避免真实实验所带来的开销，方便地调节光学、声学以及各项系统参数，同时可以有选择地忽略真实实验的系统噪声和测量误差，为理论研究提供有利条件。因此，光声成像仿真技术对于光声成像的理论研究及其在实际中的应用有重要的指导意义，也是验证实验设计性能和算法性能的首选方案。

目前，计算机仿真技术已经被广泛应用于光声成像领域的研究中。例如，TREEBY等[8]基于Matlab平台使用k空间伪谱法（k-space Pseudospectral Method）对光声信号在有损介质中的传播进行了模拟仿真；YUAN等[9]使用有限元方法在Matlab平台上对光声成像的后向过程进行了模拟。这部分研究中的仿真主要是针对声场仿真或重建过程仿真，并不涉及光源和光场部分的仿真，因此可扩展性和适用性都受到了较大限制。对于包含光场和声场的全过程光声成像仿真，一般使用商用多物理场仿真软件COMSOL Multiphysics或者有限元分析软件ANSYS搭建仿真平台进行仿真。例如，MOON等[10]使用ANSYS对金属元素颗粒进行了全过程的光声成像仿真；SOWMIYA等[11]使用COMSOL对光声断层成像的全过程进行了仿真。而COMSOL或ANSYS软件不仅价格高昂，专业性强，而且由于商业软件的非开源性，内部代码和运行机理无法得知。

针对上述问题，本文根据有限元法和k空间伪谱法[12]，借助开源工具箱k-Wave[13]有限元计算工具箱TOAST++[14]，设计并实现了光声成像仿真平台和图形用户界面(GUI) ，完成了对光声成像全过程的模拟和仿真。光声成像仿真平台是基于Matlab搭建完成的。Matlab是由美国MathWorks公司出品的商业数学软件，也是目前使用最普遍的数学软件之一。Matlab内部集成了极为完善的科学计算工具和图形显示工具，并且支持python、C/C++等其它扩展语言。因此基于Matlab搭建的仿真平台将更具扩展性和普及性。

1 光声成像理论基础

在光声成像技术中，一般使用周期性的脉冲激光照射目标组织，吸收了激光能量的目标组织会由于升温而发生周期性的热胀冷缩，从而辐射出超声波。一般而言，为了有效产生超声波，激励源需要满足热学限制条件和压力限制条件：① 脉冲激光持续时间小于组织的热扩散时间；② 脉冲激光的脉宽小于声波在目标组织中传播所需的时间。在满足这两个条件的情况下，目标组织由于激光照射而产生热量的过程可以描述为：

其中，ρ是组织密度，Cp是比热，T(r, t)表示位置r处、t时刻生物组织吸收光能量产生的温度变化，λ是热导率，H(r, t)是激光的能量沉积。由于热学限制条件，这个过程可以当作绝热过程，所以简化为：

如果假设目标组织内部的密度和声速为常数，那么可以通过以下公式来描述声场：

其中，矢量函数u(r, t)表示声位移。同时，由于温度上升导致的热膨胀可以用热膨胀方程来描述：

其中β为热膨胀系数，c是声速。联立(2)～(4)式即可得到均匀生物组织中的光声波动方程：

特别地，由于实际应用中激光的脉宽很短，所以可以将激光的时间包络看作是阶跃函数δ(t)，此时H(r, t)=A(r)δ(t)，A(r)是光吸收系数分布函数。那么在位置r，初始时刻的t0的光声压强可以写作：

其中，Γ=c2β/Cp称为格日尼森系数（Grüneisen Coefficient），表示的是光能量沉积转化为光声波的效率。Γ是一个常数，室温下大部分生物组织的Γ在0.11左右[13]，这说明光声效应产生的初始声压正比于光吸收系数分布。据此，重建生物组织光吸收分布的过程可以简化为重建初始光声压强的过程。而重建初始光声压强的过程，则是通过采集到光声信号后，通过特定的重建算法完成的。至此完成了光声成像的全过程。

2 光声成像仿真平台的建立

仿真平台主要可以分为四个部分：光学仿真模块、声学仿真模块、图像重建模块以及数据输入输出模块。仿真平台的工作流程如图1所示。

图1 光声仿真平台工作流程图Fig.1 Flow chart of the photoacoustic simulation platform

2.1 光学仿真模块

光学仿真模块不仅可以模拟和计算光在组织中的扩散，同时可以设置光学仿真过程中的各项参数，包括光源的属性、目标组织和介质的光学参数等。由于光在组织中的传输和扩散属于光声成像的前向过程，因此正确计算并得到吸收光的光场分布至关重要。光学仿真模块使用有限元方法（FEM）对光扩散方程进行求解从而得光场分布。所以在使用光学仿真模块进行计算前，需要先对目标组织进行离散化和网格划分，得到标准的网格文件（*.msh），这一过程可以利用AutoCAD、Gmsh等专业制图软件来完成。值得注意的是，在有限元方法中，网格划分的好坏是影响求解精度和速度的关键因素。网格划分应首先考虑网格的数目，一般来说，网格数量越多计算精度越高，但计算速度也就越慢，所以在确定网格的数目时应该有所权衡综合考虑。此外，网格划分还有一些其它原则，例如：在硬件条件能达到要求的情况下，最好使用形状规则的网格单元；在计算时数据变化梯度较大的部位（比如组织的边界）需要采用比较密集的网格；对于研究中需要重点关注的部位（例如待重建的目标组织）使用更为精细的网格；当几何模型的结构和形状对称时，网格也应划分成对称网格等。

2.2 声学仿真模块

光声成像理论中，在满足压力和热力学限制的条件下，利用格日尼森系数可以将组织的光吸收分布转化为初始声压分布，进而对声在组织体内的传播进行建模从而得到声场分布。为了对声场进行较精确的仿真，如果继续使用有限元方法则需要进一步提高网格的数目和精度，除了会提高仿真平台的运行时间外，也会使仿真平台占用更多的硬件资源。因此，声场仿真的数值方法采用k空间伪谱法。k空间伪谱法用于模拟复杂介质中光声信号的传播，和传统的仿真方法，例如有限元和有限差分法相比，不仅速度快而且计算结果也较为精确[15]。利用开源的Matlab声学仿真工具箱k-Wave，可以较方便地在Matlab上使用k空间伪谱法对光声信号的传播进行计算。

由于k空间伪谱法在计算时要求网格必须是均匀网格，而光学仿真模块中使用的有限元网格为非均匀网格，这就导致了两种网格在网格形状和节点数目上并不一致。所以想要在声学仿真模块中使用光学仿真模块的结果，首先要进行网格的统一[16]。相比于k空间伪谱法使用的均匀网格，有限元网格更加精细，因此需要选取均匀网格中每个网格节点所对应的多个有限元网格的节点，取它们的平均值作为均匀网格中该网格节点的实际值。这样就可以将光学仿真模块的计算结果施加在声学仿真模块建立的模型上。

此外，声学仿真模块的另一个功能是定义声场仿真中需要的各种参数，例如组织介质和目标组织的声学相关属性，不同组织区域的声速和密度等。

2.3 图像重建模块

利用声学仿真模块采集到光声信号后，使用图像重建模块及其图形用户界面，可以非常便捷地对采集到的时域光声信号进行解析法重建，例如快速傅里叶重建算法和时间反演重建 (Time Reversal Reconstruction, TR)等。此外，图像模块集成的重建算法还有：用于求解大型线性方程组的经典迭代算法LSQR算法；基于L1范数正则化的稀疏重建算法YALL1算法[17]。特别地，针对实际应用中使用广泛的其它迭代算法，该模块虽然不能直接进行重建，但是可以生成并导出这些迭代算法在重建时所需要的系统测量矩阵和探测器数据以供使用。

2.4 数据输入输出模块和图形用户界面

数据的输入输出也是仿真平台中非常重要的一个部分。这部分主要包括模型和网格的导入，仿真数据的导出等，这些操作都可以通过图形用户界面来完成。与此同时，用户图形界面还负责对各个模块进行控制以及进行交互显示等重要工作。图形用户界面不仅使仿真平台的工作过程更加直观和简洁，同时也极大地提高了仿真平台的易用性和交互性，有利于平台的进一步开发和推广。借助于Matlab的图形界面工具集GUIDE，完成了对仿真参数设置、仿真过程控制、仿真结果展示以及数据导入导出等操作的图形化和窗口化。图2所示为仿真平台的主界面以及光学仿真模块的操作界面。

图2 仿真平台主界面以及光学仿真模块控制界面Fig.2 Interfaces of main program and the optical module

3 仿真平台的可靠性验证和重建结果分析

3.1 验证实验设计

通过设计仿真实验，分别从前向过程、重建过程两方面对仿真平台的可靠性进行验证。

前向过程的可靠性验证。使用仿真平台的光学和声学仿真模块进行光声成像的前向过程仿真，得到光声信号时域波形。通过对比仿真波形和理论计算波形的形状结构来验证仿真平台前向过程的可靠性。

重建过程的可靠性验证。使用光声成像仿真平台的图像重建模块进行图像重建。通过对比重建结果和目标物的实际位置来验证重建结果是否准确。

3.2 前向过程的可靠性验证

通过对比理论计算结果和仿真结果的时域波形结构，可以验证光声仿真平台的前向过程的仿真是否可靠。理论研究表明，使用短波脉冲光源照射球状目标组织产生的时域光声信号具有“N”型结构[18]。因此实验中首先建立球状目标组织的几何模型，进行离散化和网格划分后导入仿真平台进行仿真，并将采集得到的光声信号进行对比。图3分别展示了光声波形的理论计算结果和仿真平台的仿真结果。

图3 时域光声信号的理论计算结果和仿真结果对比Fig.3 Comparison of theoretical temporal waveform and simulation result of photoacoustic signal

3.3 重建过程的可靠性验证

图4 所示为验证实验中使用的仿体几何模型和离散后的网格模型。仿体分为两部分，大圆柱体为组织介质，其底面半径为15 mm，高为60 mm；内部的小圆柱体为目标物，底面半径为4 mm，高为6 mm。目标物和介质的光吸收系数的比值设置为20:1。共16个点光源均匀分布在以大圆柱体中心为球心，半径为30 mm的球体表面。32个超声探测器组成超声探测阵列，均匀设置在大圆柱体的表面。为了直观显示目标物的重建结果和真实位置，选取了了z=30 mm平面的重建结果，如图5所示。

3.4 验证结果分析

图3展示了一个周期内光声信号时域波形的仿真值和理论值。其中图中虚线为理论计算结果，实线为仿真平台得到的仿真值。当目标物体为实心球状物时，仿真得到的光声信号波形和理论计算值基本重合，不但呈现出了“N”型的结构，而且从信号出现到消失的持续时间也均为4 μs左右。这验证了仿真平台前向过程仿真的可靠性。

图4 实验模型Fig.4 Experimental model

图5 重建结果和目标物的真实位置（z=30 mm平面）Fig.5 Reconstructed results and real position of target (plane z=30mm)

图5 展示的是重建结果和目标物的真实位置。图中的亮团是仿真平台的重建结果，目标物的真实位置用实线进行了标注。重建过程将大圆柱体离散成为长方体，所以仿体的真实边界在图中也以实线标出。由图5可知，重建出的亮团虽然存在着一定程度的伪影，但是整体呈现出较为规则的圆形形状，这和目标物本身的形状相符。此外，亮团最亮的部位和目标物的真实位置基本重合，这说明仿真平台的重建结果是准确的，因此证明了仿真平台后向过程的可靠性。

4 总结与展望

本文阐述了基于Matlab的光声成像仿真平台的设计和实现。对图形用户界面以及各模块的工作原理进行了介绍。通过仿真实验，从前向过程和后向过程两个方面证明了仿真平台的可靠性，同时提供了利用该平台进行光声成像仿真和重建的简单演示。与使用COMSOL，ANSYS等商业软件相比，该仿真平台可以大大降低仿真的成本，在非复杂实验的场景下提高仿真效率。同时，该仿真平台的开源特性也使其具有更好的可移植性和可维护性。为光声成像技术的理论研究提供了一个较为可靠和高效的研究工具。但是距离成为一个成熟的仿真工具，该仿真平台仍需进一步地改进和完善。针对该仿真平台的改进工作将会主要集中在以下几个方面：①集成几何建模和网格离散的功能，使仿真平台更具独立性；②提供更多的可设置系统参数，进一步细化仿真过程以应对更加复杂的仿真场景；③为不同类型的真实实验提供更加丰富的重建算法并提高重建精度。