医学光声成像的物理原理与仿真实现

薛泽鹏 徐家琦 刘炳坤 卫子瑶 张本健 林辉

摘要：光声成像利用生物组织的光声效应进行成像，具有分辨率高、对比度强、穿透性深的优势，在肿瘤医学检测方面具有巨大发展潜力。本工作基于热声转换原理与软件，对圆盘源、图片源及具有高含氧量肿瘤区域的脑部CT图片源的图像重建进行了仿真实现，对比了快速傅里叶变换和时间反转重建算法。研究工作表明两种图像重建算法均可重建复杂图像源。

关键词：光声成像;时间反转;快速傅里叶变换;含氧量

中图分类号：TP18        文献标识码：A

文章编号：1009-3044（2019）23-0277-03

开放科学（资源服务）标识码（OSID）：

光声层析成像（Photoacoustic Tomography）是基于光声效应发展而来的一种新型的无损成像技术。通过短波激光作为光源照射物体，在激励激光脉冲的持续时间远小于物体对激光的吸收和传导时间时，物体内部会发生热膨胀，进而发出球面的超声波[1]。医学光声成像根据生物组织对照射激光吸收的声压分布情况对组织结构进行成像。当生物组织受到短脉冲电磁照射导致瞬间热膨胀，从而发射相应频率和幅值的超声波信号，然后用时间反转算法或快速傅里叶变换算法重建图像，将生物组织内的微波能量吸收分布重建出来。

通过对医学光声图像的分析可用于对患者肿瘤等疾病的诊断[2]。由于癌组织内不同的化学成分具有不同的光学吸收特性，光声成像利用光致超声的成像机制，选择特定的激光波长对癌组织内化学成分进行检测，从而对癌组织进行高对比度的结构成像。肿瘤细胞的一个重要标志是过量的氧燃烧（即高代谢），高代谢导致组织内血红蛋白浓度、血氧浓度和氧代谢率等发生改变[3]。光声成像可对癌组织的这些参数进行定量分析，从而实现癌组织功能成像[4]。因此光声成像技术在软组织肿瘤的早期无损诊断、定位及分析上是行之有效的检测方法，在未来肿瘤的临床应用上具有潜在的应用价值。

本文首先介绍了光声成像的基本原理和重建算法，并采用改编的k-Wave 软件对圆盘源、图片源及具有高含氧量肿瘤区域的脑部CT图片源进行了图像重建，对比了快速傅里叶变换（Fast-Fourier-Transform， FFT）和时间反转（Time-Reversal， TR）重建算法各自的优势。

1 医学光声成像原理

生物组织受到短脉冲电磁照射会瞬时吸热膨胀，从而发射相应频率和幅值的超声波信号，通过超声传感器扫描收集样品发出的热声信号，用图像重建算法就可以把生物组织的微波能量吸收分布重建出来。热声波动方程（1）是热声成像的基本依据：

式中：p（r，t） 是时刻t位置r处的声压;为体膨胀系数;T（r， t） 是媒质由于吸收微波能量而引起的温升;和分别为媒质的密度和声速。

假设在一个小范围内微波能量沉积时间很短，以至于可以忽略热扩散，则可以將温度变化描述为

式中：

将式（2）～（4）代入方程（1）中，得

热声或光声成像本质上就是利用探测到的热声信号p（r， t） 来反演出生物组织中的微波能量吸收分布（SAR）。

2 计算机仿真方法

2.1 k-Wave软件

k-Wave是用于Matlab和C++的开源声学工具箱，由Bradley Treeby和Ben Cox（伦敦大学学院University College London）和Jiri Jaros（布尔诺科技大学Brno University of Technology）开发。该软件设计用于复杂和组织真实介质中的时域声学和超声模拟。仿真功能基于k空间伪谱方法，具有快速、易用的特点。该工具箱包括能够模拟声波传播压力和速度源进行建模的高级时域模型，优化c++代码，最大限度地提高大型仿真计算的性能[5]。

2.2 肿瘤组织的光吸收处理

k-Wave仿真函数主要解决的是超声场的一阶方程组，分别为动量守恒、质量守恒、声压与密度的关系函数。通过分析k-Wave基本方程，可以加深对k-Wave工具箱工作原理的理解，便于进行编程实现相关的模拟仿真。

在运行过程中，由于网格划分问题容易使运行繁琐复杂，同时为了解决数值分散问题，需要设置小时间步长。由于伪谱法采用了有限差分计算方法，并吸收了傅里叶级数的思想，经分析K-space算法能够解决这一问题[6]。采用傅里叶级数具有两大优势，即通过傅里叶变换可以精确计算傅里叶分量的幅值，且由于采用正弦曲线作为基础函数而可以使用更少的节点来保证计算精度。

结合上述相关可视化算法和理论，基于k-Wave仿真工具箱，以工具箱提供的2D-Time- Reversal-Reconstruction-For-A-Linear-Sensor例子为基础，对其进行了代码改编和模型植入。研究了快速傅里叶变换和时间反转两种图像重建算法对圆盘源、图片源，以及具有高含氧量肿瘤区域的脑部CT图片源的重建效果。

3 算例与结果

3.1圆盘源

设水中有两个圆盘声源同时传播，在空间相遇并相互穿越，最后声信号被置于水模体边缘的线性传感器探测到。两个圆盘波源半径分别为2.5mm 和 4mm， 在2D平面上相距42mm。

（b） FFT算法重建图

（c） TR算法重建图 ]

图1给出了圆盘源的原图与重建图，可见TR算法重建图（图1（c）） 相对于FFT算法重建图（图1（b）） 出现了更多的伪影，需要后期加工去除。在重建速度方面，FFT算法是使用傅里叶变化的公式进行运算，而TR算法相当于把传播到水模体边缘的声信号作为新的波源，重新进行一次传播仿真，因此使得TR算法的时间几乎是FFT时间的两倍。对于图像尺寸较小，或图像较为对称（例如点波源、圆盘源），FFT算法较TR算法具有优势。

3.2 图片波源

为校验程序对离散的图片波源的重建效果，本工作还在原程序中植入了图片波源的读入代码，可读入jpg/png/bmp 格式的图片，并将图片网格的灰度值作为各个离散波源的强度值， 重建了二维小熊猫图像。

图2 给出了小熊猫原图与重建图，图2（b）为FFT算法重建的图像，图2（c）为TR算法重建的图像。可见，图2（b）相比于图2（c）出现了更多的伪影。因此对于图片源，即在组织内部为非均匀声强的情况下，TR算法在图像重建方面性能较优。因为相对其他重建算法，TR算法的假定初始条件更少，不需要设定声介质的初始速度。TR算法可以对任意形状的封闭表面进行测量，不受测量表面外的声源影响，甚至在一些假设条件不成立如非均匀的声速分布时，该算法仍可以得到较为理想的结果，因此TR算法被认为是具备“最少约束条件”的光声图像重建算法[7]。

3.3 组织氧含量仿真

光声信号强弱是由光能的吸收和分布决定的。血红蛋白存在于生物组织中，氧的输运需要通过血红蛋白。与氧结合的血红蛋白称为氧合血红蛋白，表现为红色;而脱氧的血红蛋白呈蓝紫色。血红蛋白的含氧量与生物组织的代谢密切相关，在不同的组织中具有不同的氧饱和度。肿瘤区域含氧量高，吸收光能也相对较高。

我们通过改变图片区域的颜色来模拟含氧量的差别，即使用Matlab通过更改图片的RGB256位数值将图片部分区域的颜色更改为蓝色和红色来模拟对光吸收能力的高低。

在原始图像中，由于背景（深色）吸收光的能力要比目标区域高得多。而背景产生的光声信号先于目标区域的信号到达传感器，导致目标区域的光声信号被背景区域的信号所覆盖，导致重建图像严重失真。为了避免这个问题，我们在将原始图像输入代码之前，添加了一个灰度反转程序。

灰度变换是图像增强的一种经典而有效的方法。其原理是将图像中的每一个像素的灰度值，通过一个函数对应到另一个灰度值上去，从而实现灰度的变换。常见的灰度变换有线性灰度变换和非线性灰度变换。对灰度图像进行线性灰度变换能将输入图像的灰度值的动态范围按线性关系公式拉伸扩展至指定范围或整个动态范围。

我们对一脑部肿瘤CT图像进行了图像重建。在模拟重建过程中，利用Matlab 选中原图黑色区域（rgb（0，0，0））替换为蓝色（rgb（0，0，255））（图中为深灰色）（图3（a））和红色（rgb（255，0，0））（图中为浅灰色）（图3（b）），利用灰度值的差异模拟蓝色、红色两种颜色对光能的吸收强弱区别，来模拟低氧浓度和高氧浓度两种情况，然后进行TR算法重建得到图3（c）和（d），观察氧浓度对重建的影响。

可见，浅色组的成像内容多在黑区的一部分，这也证明在光声成像中，光吸收能较强的区域可以更明显地显现出来。因此由于肿瘤细胞含氧量较高，肿瘤区域也可以通过模型很好地重建成像。浅色组重建图像中肿瘤区域部分灰度值要高于深灰色组。重建图像的灰度值与组织氧浓度水平之间呈非线性关系。

4 结论

光声成像利用生物组织的光声效应进行成像，在医学癌症检测方面具有巨大發展潜力。本工作基于热声转换原理与Matlab工具箱k-Wave软件，对圆盘源、图片源及具有高含氧量肿瘤区域的脑部CT图片源的图像重建进行了仿真实现。对比了快速傅里叶变换和时间反转重建算法的优缺点。结果表明两种图像重建算法均可重建复杂图像源，其中对于图像尺寸较小，或图像较为对称（例如点波源、圆盘源），FFT算法较TR算法具有优势。而对于复杂的非均匀强度源TR算法的精度更高。

致谢：本工作得到合肥工业大学大学生创新训练计划项目（编号：2018CXCYS209）的支持。

参考文献：

[1] Bell A G. On the production and reproduction on of sound by light[J]. American Journal of Science， 1880（118）： 305-324.

[2] Fakhrejahani E， Torii M， Kitai T， et al. Clinical Report on the First Prototype of a Photoacoustic Tomography System with Dual Illumination for Breast Cancer Imaging[J]. Plos  One， 2015， 10（10）：e0139113.

[3] Gottschalk S， Fehm T F， Deanben X L， et al. Noninvasive real-time visualization of multiple cerebral hemodynamic parameters in whole mouse brains using five-dimensional optoacoustic tomography[J]. Journal of Cerebral Blood Flow & Metabolism Official Journal of the International Sociely of Cerebral Blood Flow & Metabolism， 2015， 35（4）： 531-535.

[4] Wang L， Maslov K， Xing W， et al. Video-rate functional photoacoustic microscopy at depths[J]. Journal of Biomedical Optics， 2012， 17（10）： 1060071-1060075.

[5] Treeby BE， Cox BT. k-Wave： MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. J. Biomed. Opt. 2010; 15（2）：021314.

[6] 胡文娟.基于k空间的并行磁共振成像算法的研究[D].东北大学.2011

[7] 孙明健. 光声成像系统信号采集与图像重建研究 [D]. 哈尔滨：哈尔滨工业大学，2011：3-5.

【通联编辑：唐一东】