定量光声层析成像的研究进展

孙 正,郑 兰

定量光声层析成像的研究进展

孙 正*,郑 兰

(华北电力大学 电子与通信工程系，河北 保定 071003)

光声层析(Photoacoustic tomography,PAT)成像结合了超声成像的高分辨率和光学成像的高对比度的优势，是一种新型的生物医学成像模式。PAT成像算法包含两个逆问题，即根据组织产生的光声信号构建初始声压分布图(即图像重建)以及在此基础上估算成像区域的光学特性参数。后者是一个非线性的不适定问题，通常称为定量光声层析(Quantitative photo-acoustic tomography，qPAT)成像。本文在介绍光声成像原理的基础上，对主要的qPAT算法进行综述，讨论各自的优势和不足，并对未来可能的发展方向进行展望。

定量光声成像； 图像重建； 逆问题； 光吸收系数； 散射系数； Gruneisen系数

1 引 言

光声层析(Photoacoustic tomography,PAT)成像是近年来发展起来的一种非电离式的新型生物医学成像方法。它采用脉冲激光对生物组织进行激发，获取组织的光吸收系数分布情况，因此继承了光学成像的一系列优点，如对组织损伤小、对比度高、可进行功能成像和分子成像等[1]。PAT探测的是组织产生的超声波信号，因此它继承了超声成像对深层组织成像的高分辨率的特点[2]。与X射线成像、超声成像和磁共振成像等传统的医学成像技术相比，PAT技术的优势主要体现在：采用非电离波段，成像过程中不改变生物组织的属性；较容易界定组织产生的光声信号和组织的生理状态之间的关系；可与超声或光学成像技术相结合，以获得更多的诊疗信息；可根据实际应用的需要调整成像深度和成像分辨率等[3]。

目前，对生物PAT成像系统的改进主要体现在三个方面：第一，探测组织产生的光声信号的方法，除了采用传统的基于压电效应的超声换能器或者PVDF 膜的水听器之外，还可利用光学方法实现光声信号的探测。例如：探测光声效应引起的组织折射率的变化[4]、探测压力传感器的变形[5]或者采用附着在光纤束末端的FP(Fabry-Perot)聚合物薄膜作为超声传感器[6]；第二，提高成像速度，例如：通过提高脉冲激光光源的重复频率，加快光声信号的产生过程，缩短信号采集的等待时间[7]，或者提高光声信号的采集速率从而节省数据的采集时间[8]等；第三，成像分辨率从超声分辨水平发展到了光学分辨水平[9]。

PAT成像的逆问题包括声学和光学两方面：声学逆问题是指根据探测器采集到的光声信号(本质是超声波)重建组织内部的初始声压分布或空间光吸收能量密度，即一般意义上的光声图像重建[10-12]；光学逆问题是指运用合适的光传输模型与优化算法，根据探测到的光声信号或者光吸收能量密度，估算出组织的光学特性参数(包括光吸收系数和散射系数)的空间分布[13]和热膨胀特性参数(Gruneisen系数)等[14-16]，得到对组织光学特性的定量评价，即定量光声层析成像(Quantitative photoacoustic tomography,qPAT)。由于声波在组织中传播的时间比激光照射组织以及组织吸收光的时间长约3个数量级，因此可以将声学逆问题和光学逆问题分离开来分别求解[16]。声学逆问题的重建结果，即光吸收能量密度是由局部的光吸收系数和光子数分布共同决定的，并不能反映组织的光吸收系数分布，而光吸收系数与组织的化学成分密切相关，正常和病变组织的光学特性参数之间通常有较明显的差异，因此qPAT可为疾病的早期诊断提供更加准确可靠的信息。

目前对qPAT的研究已成为光声成像领域的热点之一，尤其是同时重建光吸收系数和散射系数的空间分布。本文在简介光声成像和qPAT原理的基础上，对目前已提出的主要qPAT算法进行总结和归纳，简介其中典型算法的原理，并讨论各自的优势和不足。

2 光声成像的原理

生物PAT成像是一种以超声波作为媒介、以组织的光吸收系数和散射系数作为成像参数的生物光子成像方法，其物理基础是生物组织的光声效应[3]。成像系统包括3个组成部分：光声信号的产生单元、接收单元和后处理单元。一束短脉冲(～10 ns)激光经过光学元件后照射到生物组织上，在组织内部形成与组织光学特性参数相关的能量沉积分布。由于激光脉宽很窄，组织吸收的光子能量不能在短时间内释放，导致组织温度瞬间的不均匀提升。周期性热流使周围的介质产生热膨胀从而激发出宽频带的超声波，并迅速向组织边界传播。基于这种超声波信号的特殊产生机理，为了区别于其他的超声信号，通常称其为光声信号[17]。超声换能器接收到光声信号之后，经计算机处理即可重建出组织内部光能量沉积的分布图，揭示病变组织的内部信息。

如图1所示，从组织内的吸收体吸收光子到探测器测得声压时间序列的整个过程称为光声成像的正问题，可分为光学正问题和声学正问题两部分：前者的结果是光吸收能量密度，后者的结果是探测器测得的声压信号。假设采用单一波长λ的光源照射组织，组织的光吸收系数和散射系数分别为μa和μs，待测组织区域Ω内一点r处的光吸收能量密度H(r,λ)为

H(r,λ)=μa(r,λ)Φ(r,λ;μa(r,λ),μs(r,λ)),

(1)

其中，H(r,λ)定义为单位体积单位时间内的热能转换；r∈Ω，Ω⊂Rn(Rn为n维实数空间，n=2,3)为有界域，Φ是光能流率，则r处的初始声压为

(2)

图1 光声成像正问题框图Fig.1 Photoacoustic imaging block diagram of the forward problem

3 定量光声成像的原理

3.1 声学逆问题

假设介质的声学特性均匀，在理想激光脉冲的均匀照射下，被照组织产生的三维光声信号的幅值与脉冲激光的幅值成正比，光声信号的特性由光能量的吸收分布决定。因此，可以根据探测器测量到的声压时间序列p(r,λ,t)重建初始声压的空间分布p0(r,λ)，进而得到光吸收分布H(r,λ)，即PAT声学逆问题，也就是通常所说的光声图像重建。

如果超声波在密度均匀的无损介质中的传播速度是均匀的，并且光激发可被看作是瞬时的，那么超声波传播的物理模型可由以下3个方程表示[9]：

(3)

式中，p(r,t)为r∈Ω⊂Rn处、t∈R+时刻的声压；u(r,t)为介质的振动速度；c为声速；ρ(r,t)为位置r处、时刻t时的声学密度；ρ0是介质密度。

从不同的应用角度出发，借鉴其他成像技术，研究者已提出了多种图像重建算法，例如滤波反投影法(Filtered back-projection,FBP)、时间反演法(Time-reversal,TR)、相控聚焦算法、基于傅里叶变换的重建算法、反卷积重建算法和迭代重建算法等[12]。

3.2 光学逆问题

PAT光学逆问题是指由初始声压分布p0(r,λ)估算组织的光吸收系数和散射系数的空间分布。在早期的研究中，通常假设组织的光散射系数是已知的，那么采用递归法[18]或者非递归的方法[19]即可重建出光吸收系数的分布。但是，该假设在多数情况下都是不成立的，更为通用的方法是基于误差最小化的方法。由式(2)可知，若待测组织内部的光吸收系数和散射系数是有界常数且满足Lipschitz连续，且已知Gruneisen系数，光吸收能量密度的测量值为

(4)

那么，光学逆问题可表述为如下的非线性最小二乘问题[20-21]

C·H(r,λ;μa(r,λ),μs′(r,λ))‖2)， (5)

其中，μs′是组织的约化散射系数

μs′=μs(1-g)，

(6)

光能流率Φ是未知的，而且它与组织的光吸收系数和散射系数有关，同时PAT成像是三维高分辨率成像，所涉及的数据量极大，因此由光吸收能量密度的测量值重建组织的光学特性参数是一个大规模的非线性不适定问题，特别是当同时重建光吸收系数和散射系数时[22]。

求解光学逆问题需要解决两个关键问题[20]：第一，建立光在组织中传输的前向模型，模型的输出就是光吸收能量密度的理论值；第二，选择适当的优化策略，使光吸收能量密度的测量值与前向模型的输出值之间的误差最小，进而得到光学参数的估计值。

3.3 光在组织中传输的前向模型

通常采用辐射传输方程(Radiative transfer equation，RTE)描述光子在混浊介质中的迁移过程[22]。RTE是积分-微分方程，求解时常需要在空间域和角度域内对方程进行离散化，步骤较为繁琐，因而通常对其进行扩散近似(Diffusion approximation，DA)。相比于DA，RTE能够更准确地描述光子在组织中的迁移过程，尤其是在非扩散区域，但是其复杂的求解过程和较高的运算成本限制了它的广泛应用。

3.3.1 辐射传输方程

RTE描述了在特定控制体(Control volume)内的能量守恒：当光在待测区域内沿某一方向传输时，组织对光子的吸收、偏离光传输方向上光子的散射和超过待测区域的光子外流出都会导致能量的损失，而沿光传输方向的光子散射或介质中的其他光源又会造成能量的增加。

(7)

(8)

假设除了光源Γs⊂∂Ω处以外，在待测区域边界∂Ω处光子不会向内传输，那么方程(7)的边界条件是

(9)

(10)

其中，Φ是光能流率，也称为光辐射通量。

3.3.2 扩散近似

假设组织的约化散射系数远大于吸收系数，即光在组织中发生散射的几率远高于吸收(对于人体的大多数组织，该假设是成立的)，并且光在组织中的传输和分布近似于各向同性，则可将RTE简化为扩散方程(Diffusion equation,DE)，它是用球谐函数表示的RTE的一阶相位近似。

在DA模型中，光子密度被简化为只与空间变量r有关，设有界域Ω的光滑边界为Γ，则满足Robin边界条件的时间独立的DE可写作[14,20]

-·(D(r)Φ(r))+μa(r)Φ(r)=S(r),r∈Ω，

(11)

其中，S(r)是源项，D(r)是扩散系数：

(12)

为了求式(11)的数值解，需要确定适当的边界条件。例如可采用Dirichlet边界条件[20]，即在外推边界处(距离介质边界2D处)，令Φ=0。计算出Φ(r)后，则

H(r)=μa(r)Φ(r),

(13)

DE的求解方法相对较简单，只需在空间域内对方程离散化即可，因而它是目前qPAT中最常用的前向模型[14]。在较为理想的情况下，例如各向同性的浑浊介质中嵌入柱状的不均匀物质，求解DE可以得到光吸收能量密度的解析解[24]。对于实际的生物组织，则需要采用有限差分法或者有限元法得到DE的数值解。

DA仅在扩散域内有效，即在距离光源几个迁移平均自由程(Transportmeanfreepaths)的区域内μa=μs′，在该区域内光由于发生多次散射，因而失去方向性而扩散开。但是对于PAT成像，接近光源的区域和待测区域的边界构成了图像的重要组成部分，通常包含诊疗疾病所需的重要信息，在这些区域内，光的传输是高度各向异性的，因此DA是不成立的[23]。

3.4 主要的优化方法

为了得到光在组织中传输的前向模型的数值解，一般需要先对其进行有限元离散化，相比一般的单网格方法，采用双网格方法可在保证重建精度的前提下，明显缩短重建时间，提高计算效率[20,25]。在此基础上求解式(5)的非线性最小二乘问题，使光吸收能量密度的测量值与其理论值之间的均方误差最小，即可求得光吸收系数和散射系数的空间分布。最常用的方法是Levenberg-Marquardt(LM)方法[26]，即L2范数Gauss-Newton法[27]，通过计算Hessian矩阵的逆矩阵，可迭代地调整(μa,μs)的值。但是该计算过程非常耗时，因此出现了近似计算逆Hessian矩阵的算法，主要包括基于Jacobian矩阵的线性方法[21,28]、非线性梯度法[23,14,29-32]、Bregman迭代法[15]等，那么光学逆问题解的精度将很大程度上取决于所选数值模型的精度。

3.4.1 基于Jacobian矩阵的线性方法

考虑到Gauss-Newton法存在运算成本过高的问题，可采用吸收系数和散射系数的Jacobian矩阵构建Hessian矩阵的近似矩阵，用于在Gauss-Newton迭代中更新(μa,μs)的值。此类方法的优点是原理简单，易于实现，而且它保留了Gauss-Newton法的二阶收敛性。但是，PAT是三维成像，数据量巨大，因此得到的Jacobian矩阵是大型密集矩阵，对计算速度和存储空间的要求都很高，而且计算该密集矩阵的逆矩阵也非常耗时，可能降低此种方法的实用性[23]。文献[27]提出一种无矩阵(Matrix-free)求解近似Hessian矩阵的方法，即只计算矩阵向量积，可降低对存储空间的要求，可是计算量仍然很大。

3.4.2 非线性梯度法

此类方法采用基于梯度的拟牛顿(Quasi-Newton)最小化法[33-35]求解式(5)的最优化问题，避免计算Jacobian矩阵和Hessian矩阵，而只需计算式(5)中误差函数的梯度，从而迭代更新光吸收系数和散射系数的值。通常需要引入正则项，降低问题的病态性，使之成为良态问题，如L2正则化。采用多位置激光照射的方法采集多个光声信号数据集，在已知Gruneisen系数的前提下，采用此类方法可同时、唯一地确定光吸收系数和散射系数。与基于Jacobian矩阵的方法相比，该类算法具有超线性收敛性，且计算速度快，所需的存储空间小，对三维图像的重建更具实用性。但是其计算较复杂，一般需要较多的迭代次数。

3.4.3 正则化方法

对于病态问题，采用正则化方法(即引入正则项)可降低问题的病态程度，把非线性问题转化为某种条件下的线性问题。因此在解决qPAT这一非线性病态逆问题的过程中，选取合适的正则化方法可使算法更加简单，保证稳定近似解的存在[36]。对正则化方法的选取依赖于光吸收系数和散射系数的先验知识。

(1) L1和L2正则化法

(2) Tikhonov正则化方法

基于变分原理的Tikhonov正则化方法[39]是研究不适定问题的最重要的正则化方法之一，即在求解过程中结合解的某些已知信息对解进行限制。PAT是对不同的组织区域进行成像，因此对组织结构光滑性的约束也不同，这一信息可用于设定Tikhonov方法中解的光滑性假设，反映光学参数的光滑性[36]。对于光滑(连续)变化的光学参数，可采用标准Tikhonov正则化方法[39]，设定解的光滑性，反映光学参数的光滑性。PAT非常适合于对血管这种非光滑性的结构进行成像，这种情况下，将组织的光吸收系数和散射系数看作是分段连续的函数更为适合[40-42]，文献[36]设计了一种基于标准Level-set[43]的正则化方法，求解分段连续的光吸收系数和散射系数。此外，还可以考虑采用迭代的正则化方法[44]。

(3) 基于全变差(Total-variation,TV)的非光滑正则化方法

对于待求光学特性参数是分段常数的情形，带约束的全变差(或全变分)正则化方法比L2正则化法更为适合[15,21]，全变分项可以良好地约束信号的平滑性。如果将TV与Bregman方法相结合，则解的精度和计算效率都会得到大幅提升。

3.4.4 Bregman迭代法

Bregman迭代法采用基于次梯度[45]的最小化法求解式(5)的最优化问题，基本模型是ROF(Rudin-Osher-Fatemi)模型。该方法可以解决非线性的非凸优化问题，把非凸问题转换成凸问题去逼近，通常需结合其他优化方法。

目前，Bregman迭代法在qPAT中的应用还处于初步研究阶段。基本思路是：设定光吸收系数和散射系数的初值为0，在每次迭代过程中通过更新次梯度极小化Bregman距离[45]，达到目标函数最小化，最终得到光吸收系数和散射系数的最优解。该算法的计算效率高、收敛速度快。引入正则化后，正则化参数是不变的，且在迭代过程中系统稳定。

在qPAT中常采用以下两种Bregman迭代法求解最优化问题[15]：

(1)基于Jacobian矩阵的线性Bregman迭代法

采用吸收系数和散射系数的Jacobian矩阵构建Hessian矩阵的近似矩阵，用于在Bregman迭代法中迭代更新吸收系数和散射系数的值。在每次迭代过程中求解的是关于Jacobian矩阵的误差函数最小值，比直接计算Jacobian矩阵简单。也可将TV正则化与Bregman迭代法结合，在目标函数中增加由当前迭代过程产生的误差(其初始值为0)，用于修正下一次迭代，可提高结果的精度，但会增加迭代次数。

(2)基于梯度的Bregman迭代法

若式(5)可微，则可对式(5)采用L-BFGS(limited-memory BFGS)算法[15]得出光吸收系数和散射系数的初值，再引入正则化，进而进行Bregman迭代求解。L-BFGS算法通过对吸收系数和散射系数加权，克服同时重建吸收系数和散射系数对其梯度数据的敏感度问题，所需要的存储空间小。

4 定量光声成像的主要算法

4.1 按照PAT成像光源的不同分类

根据PAT成像时所用激光光源的不同可将qPAT方法分为单光源、多光源和光谱qPAT 3类。

(1)单光源qPAT

采用单一波长的单个激光光源照射组织，在已知待测组织的光散射系数分布的情况下，可以重建出光吸收系数的空间分布。若未知组织的内在散射特性，则不能唯一、准确地同时重建出光吸收参数和散射系数的分布[15]。

(2)多光源qPAT(Multiple-source，MS-qPAT)[14,28,34,46-51]

采用相同波长、不同位置的多个激光光源照射组织，采集多个光声信号数据集，进而得出多个初始声压数据集。在已知待测组织的Gruneisen系数的情况下，可唯一、准确地同时重建出待测组织边界和内部的光吸收系数和散射系数的空间分布。

(3)光谱qPAT[16,19,30,52-54]

采用多波长的激光光源在不同位置照射生物组织，获得多个初始声压数据集，在已知待测组织的光散射特性与入射光波长之间关系的前提下，可唯一地同时重建出光吸收系数、散射系数和Gruneisen系数的空间分布。其中求解光吸收系数和散射系数是非线性问题，求解Gruneisen系数是线性问题。此类方法的缺点是在每次不同波长的激光照明时都需要重新测量初始声压分布图，计算繁琐。

4.2 按照所用数据源的不同分类

根据重建光学特性参数时所用数据源的不同，可将qPAT算法归结为两步算法和一步算法两大类，其中对两步算法的研究和优化是目前的主流。

4.2.1 两步算法

两步算法指在得到声学逆问题结果的基础上，再求解光学逆问题。即假设从探测器采集到的声压时间序列中重建出的光吸收能量密度或初始声压分布是已知的，在此基础上求解光学参数的空间分布。

目前已提出的多数qPAT算法都属于两步算法，其基本思想如§3.2～§3.4中所述，首先建立光在组织中传输的前向模型(如RTE或者DA[21,47,55])，并对其进行有限元离散化，然后采用适当的优化算法使光吸收能量密度的测量值与其理论值之间的均方误差最小，求得光吸收系数和散射系数的空间分布。通过对光吸收能量密度进行尺度变换(如对数化处理)[21]，可以大大提高优化算法的收敛速度。但是qPAT测得的光强度的动态范围可能非常大，需要优化尺度变换的比例，以确保最优化问题数值解的稳定性。为了得到声速不均匀情况下的解，可采用不同波长的激光对组织进行照射，即MS-qPAT，利用采集到的多个光声信号数据集重建光吸收系数，在这个过程中超声波的传播速度是可以求解的。

矢量场法[48,53-54]也属于两步算法，其基本思想是：用相同波长的激光光源从两个不同位置照射组织，收集两个光声信号数据。假定声速均匀，待测区域内组织的Gruneisen系数以及边界处的光吸收和散射系数是已知的，待测组织内部的光吸收系数和散射系数是有界常数且满足Lipschitz连续，组织吸收的光子能量全部转化为初始声压。通过建立基于Dirichlet边界条件的DA模型，将根据两个初始声压数据集构建的向量场代入光扩散方程，通过等量变换，求解出光吸收系数和散射系数。该算法的优点是：把一个非线性问题分解成两个线性问题，非迭代地求解光学参数，计算速度快且高效。其局限之处在于：需在待测区域的边界光散射系数已知的情况下才能使用；只能解决满足Dirichlet边界条件的光散射问题；不能同时确定地重建光吸收系数、散射系数和Gruneisen系数，但假定待测区域的某一系数值已知时，可以唯一地确定另外两个系数。

混合数值重建法[53-55]是在矢量场法的基础上优化而来的，原理是：假定组织的Gruneisen系数是一个不确定的常数，声速均匀，被测组织吸收的光能量全部转化为初始声压。通过扩散光学层析(Diffuse optical tomography,DOT)成像模型获取待测区域的边界光吸收系数和散射系数。建立基于Dirichlet边界条件的DA模型，采用非线性最优化方法[36,49,55]更新边界光吸收系数和散射系数。采用矢量场法求解出区域内部的光吸收系数和散射系数，并迭代更新，直至通过RTE和DA模型生成的数据与实测数据一一匹配。该算法除了具有矢量场法的优点以外，由于用非线性优化问题来处理PAT逆问题，因而减小了未知空间，而利用矢量场法可以减小优化范围；结合DOT成像，求解出待测区域的边界系数值，克服了矢量场法的缺陷。

此外，可以利用DOT成像测得光能流率[56]，或者利用已知光吸收谱的光学对比剂来定量测量待测组织上的局部光能流率[57]。将PAT与声光调制结合起来，亦可去除初始声压测量中光能流率的影响[58]。

两步算法的主要不足在于：根据探测器接收的光声信号重建光吸收能量分布是存在误差的，所以在此基础上重建出的光学参数也是不准确的；超声波在组织内的传播速度是变化的，并不能准确获得声速，而两步算法仅在声速已知的情况下重建光学参数的误差较小；光学噪声和光散射问题也是不可忽略的。

4.2.2 一步算法

一步算法指根据超声探测器接收到的组织产生的光声信号直接重建光学特性参数，包括一步反演算法[59]、单级法[60]和免校正法[61-62]等。

一步反演算法的基本思想是采用相同波长、不同位置的光源照射待测组织，采集多个声压数据集 。假定已知组织的光散射系数和Gruneisen系数，将求解光吸收系数分布的逆问题分为线性和非线性两种情况分别进行分析。对于线性逆问题，假定已知背景组织的光吸收系数，采用Born近似法[63]和Landweber迭代法[64]，求解待测组织和背景组织的光吸收系数之间的差值，进而重建光吸收系数分布。这种方法一般适用于声速变化且未知、光吸收系数的变化相对较小的情况。对于非线性逆问题，采用最优控制法，设定一个目标函数，用Levenberg-Marquardt方法使目标函数最小，求解光吸收系数的最优解。在不同组织的边界处，超声波传播速度和组织光学系数的相对变化一般较小，因而可以看作是线性逆问题。该算法的优点是：可在超声波速未知的情况下求解光吸收系数，同时还可以求解出超声波的传播速度。但是不能同时重建光吸收系数、散射系数和Gruneisen系数[52]。

单级法的基本思想是在单一波长、不同位置的激光照射条件下[53]，采集超声探测器产生的声压数据集。假设已知Gruneisen系数，且声速均匀，使用RTE前向模型，采用有限元离散化；在包含成像目标的封闭区域内，使用广义Tikhonov正则化方法[39]求解出RTE算子参数；使用近端梯度算法使Tikhonov函数最小化重建光吸收系数。在已知光散射系数的情况下，利用该方法可以唯一地重建出光吸收系数。也可在不同波长的激光照明条件下，利用该算法重建光学参数。

免校正法的原理是：假设待测区域内的介质均匀且已知光散射系数，采用有限元离散化的方法，使用DA前向模型，根据待测边界的声压测量值和计算结果得出光吸收系数的测量值；然后，使用牛顿法和Marquardt-Tikhonov正则法[42]迭代更新光吸收系数，寻找目标函数最小化的最优值，重建出光吸收系数分布图。该算法的优点是它仅利用声压的归一化边界测量值和入射光的强度，不依赖于任何校正程序；采用标准光声数据，消除了超声探测器灵敏度的影响。此外，还可以将该算法应用于光谱光声测量中。

4.2.3 小结

两步法需要首先根据光声信号重建组织内部的光吸收能量分布或光声压分布图，再据此重建组织的光学参数。在早期的研究中，多采用一种波长的激光照射组织进行光声成像。在后续的研究中发现，不同波长的激光照射组织形成的光吸收能量分布和光声压分布图不同，这样会增加重建光学参数时的变量。对于光谱qPAT问题，也可以采用两步算法解决，分别使用多种波长的激光照射组织，再重复多次重建过程，但这样会使步骤过于繁琐，延长重建过程。目前对两步算法的优化主要是针对光谱qPAT或者MS-qPAT。由于使用不同光源照射时，生物组织的光学参数是不会改变的，因而对于MS-qPAT，采用一步算法重建光学参数可解决上述两步算法的重建过程中变量增加的问题，相对较准确地重构组织的光学参数。此外，探测器采集的光声信号是通过测量区域的边界部分获得的，据此重建出的初始声压分布图是不准确的，进一步重建出的光学参数也存在较大误差，采用一步算法则可以避免重建初始声压分布图时产生的误差，重建出相对较准确的光学参数。

5 结 论

PAT成像为研究生物组织的形态结构、生理和病理特征以及代谢功能等提供了重要的手段。生物组织的光吸收特性与其结构功能和病理特征等紧密相关，qPAT根据生物组织产生的光声信号恢复其光学特性参数，是一个大规模、非线性、病态的逆问题。

目前，qPAT面临的挑战主要包括：(1)声学方面，进一步提高重建图像的精度。(2)光学方面，在多波长激光照射下实现参数重建，考虑未知波长的依赖性影响，以及全面求解光学和声学逆问题的规模。(3)目前的研究工作多是侧重在相对理想的条件下重建组织的光学参数，需进一步考虑实际情况的复杂性，例如组织Gruneisen系数的变化、超声波在组织中的非匀速传播、光在组织表面的非均匀分布以及噪声等，使重构出的光学参数分布图更接近真实情况。(4)将二维图像重建算法扩展到三维，恢复出光学参数的三维分布图。(5)从数学计算的角度来讲，要实现对光吸收系数和散射系数空间分布的高分辨率重建，所需的计算成本非常高，尤其是重建光学参数的三维空间分布时，对计算机的存储空间和运行速度都有很高的要求。

除了光学参数之外，光声成像还能利用其他成像参量进行多尺度成像，例如：化学参量成像、黏弹参量成像、温度参量成像、速度参量成像、声学功率谱参量成像和物化谱参量成像等[65-66]，为疾病的临床诊治提供更为丰富的组织结构和功能信息。

[1] YAO J,WANG L V.Breakthroughs in photonics photoacoustic tomography [J].IEEEPhoton.J.,2014,6(2):0701006.

[2] NAM S Y,CHUNG E,SUGGSL J,etal..Combined ultrasound and photoacoustic imaging to noninvasively assess burn injury and selectively monitor a regenerative tissue-engineered construct [J].TissueEng.Part C:Methods,2015,21(6):557-566.

[3] WANG L V,GAO L.Photoacoustic microscopy and computed tomography: from bench to bedside [J].Annu.Rev.Biomed.Eng.,2014,11(16):155-185.

[4] NUSTER R,GRATT S,PASSLER K,etal..Photoacoustic section imaging using an elliptical acoustic mirror and optical detection [J].J.Biomed.Opt.,2012,17(3):99-106.

[5] LAUFER J,JOHNSON P,ZHANG E,etal..Invivopreclinical photoacoustic imaging of tumor vasculature development and therapy [J].J.Biomed.Opt.,2012,17(5):449-450.

[6] ANSARI R,ZHANG E,MATHEWS S,etal..Photoacoustic endoscopy probe using a coherent fiber-optic bundle [J]SPIEBios.,2016,142 (10):97080L

[7] BILLEH Y N,LIU M,BUMA T.Spectroscopic photoacousti microscopy using a photonic crystal fiber supereontinuum source [J].Opt.Express,2010,18(18):18519-18524.

[8] MALONE E,POWELL S,COX B T,etal..Reconstruction-classification method for quantitative photoacoustic tomography [J].J.Biomed.Opt.,2015,20(12):126004.

[9] SONG H,WANG L V.Optical-resolution photoacoustic microscopy auscultation of biological systems at the cellular level [J].Biophys.J.,2013,105(4):841-847.

[10] LUTZWEILER C,DEN-BEN X L,RAZANSKY D.Expediting model-based optoacoustic reconstructions with tomographic symmetries [J].Med.Phys.2014,41(1):013302.

[11] KUCHMENT P,KUNYANSKY L.Mathematics of photoacoustic and thermoacousic tomography [J].Eur.J.Appl.Math.,2009,2(19):817-866.

[12] 孙正,韩朵朵,王健健.血管内光声成像图像重建的研究现状 [J].光电工程,2015,42(3):20-27.SUN Z,HAN D D,WANG J J.Review of image reconstruction for intravascular photoacoustic imaging [J].Optoelect.Eng.,2015,42(3):20-27.(in Chinese)

[13] SONG N,DEUMIC C,DAS A.Considering sources and detectors distributions for quantitative photoacoustic tomography [J].Biomed.Opt.Express,2014,5(11):3960-3974.

[14] GAO H,OSHER S,ZHAO H.MathematicalModelinginBiomedicalImagingⅡ [M] Berlin/Heidelberg: Springer,2012 131-158.

[15] GAO H,ZHAO H,OSHER S.Bregman methods in quantitative photoacoustic tomography [J].Cam.Rep.,2010,30(6):3043-3054.

[16] COX B,LAUFER J G,ARRIDGES R,etal..Quantitative spectroscopic photoacoustic imaging: a review [J].J.Biomed.Opt.,2012,17(6):ID 061202.

[17] 孙正,苑园，王健健.血管内光声成像的研究进展 [J].中国生物医学工程学报,2015,34(2):221-228.SUN Z,YUAN Y,WANG J J.Progress of intravascular photoacoustic imaging [J]Chin.J.Biomed.Eng.,2015,34(2):221-228.(in Chinese).

[18] COX B T,ARRIDGE S R,KOSTLI K P,etal..Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method [J].Appl.Opt.,2006,45(8):1866-1875

[19] BANERJEE B,BAGCHI S,VASU R M,etal..Quantitative photoacoustic tomography from boundary pressure measurements: non-iterative recovery of optical absorption coefficient from the reconstructed absorbed energy map [J].J.Opt.Soc.Am.A.2008,25(9):2347-2356.

[20] LI S,MONTCEL B,YUAN Z,etal..Multigrid-based reconstruction algorithm for quantitative photoacoustic tomography [J].Biomed.Opt.Express,2015,6(7):2424-2434.

[21] TARVAINEN T,COX B T,KAIPIO J P,etal..Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography [J].InverseProbl.,2012,28(8):1067-1079.

[22] TARVAINEN T.ComputationalMethodsforLightTransportinOpticalTomography[D].Kuopio: University of Kuopio,2006.

[23] SARATOON T,TARVAINEN T,COX B T,etal..A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation [J].InverseProbl.,2013,29(7):75006-75024.

[24] LI S,MONTCEL B,LIU W,etal..Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging [J].Opt.Express,2014,22(17):20500-20514.

[25] 肖嘉莹.定量光声成像技术及在骨关节炎诊断的研究 [D].长沙:中南大学,2011.XIAO J Y.QuantitativePhotoaeoustieTomographyandItsApplicationtoTheDiagnosisofOsteoarthritis[D] Changsha: Central South University,2011.(in Chinese)

[26] KUCHMENT P.TheMathematicalLegacyofLeonEhrenpreis[M].Milan: Springer,2012:183.

[27] SCHWEIGER M,ARRIDGE S R,NISSILA I.Gauss-Newton method for image reconstruction in diffuse optical tomography [J].Phys.Med.Biol.,2005,50(10):2365-2386.

[28] COX B T,TARVAINEN T,ARRIDGE S.Multiple illumination quantitative photoacoustic tomography using transport and diffusion models [C]ProceedingsofInternationalConferenceonTomographyandInverseTransportTheory,USA,MathematicalSoc.,2011.

[29] COX B T,ARRIDGE S R,BEARD P C.Gradient-based quantitative photoacoustic image reconstruction for molecular imaging [C].ProceedingsofSPIEInternationalConferenceonPhotonsPlusUltrasound:ImagingandSensing.SPIE-IntSoc..OpticalEngineering，SanJose,California,UnitedStates,2007.

[30] BAL G,REN K.On multi-spectral quantitative photoacoustic tomography in diffusive regime [J].InverseProbl.,2012,28(2):025010-25022.

[31] MAMONOV AV,REN K.Quantitative photoacoustic imaging in radiative transport regime [J].Commun.Math.Sci.,2012,12(2):201-234

[32] SOONTHORNSARATOON T.Gradient-basedMethodsforQuantitativePhotoacousticTomography[D].London: University College London,2014.

[33] NOCEDAL J.NumericalOptimization(SpringerSeriesinOperationsResearchandFinancialEngineering)[M].New York: Springer,2006:134.

[34] BAL G,Uhlmann G.Inverse diffusion theory of photoacoustics [J].InverseProbl.,2009,6(8):1407-1426.

[35] KLOSE A,HIELSCHERA H.Optical tomographic image reconstruction with quasi-Newton methods [J].InverseProbl.,2003,19(2):387-409.

[36] CEZARO A D,CEZARO F T D.Regularization approaches for quantitative photoacoustic tomography using the radiative transfer equation [J].J.Math.Anal.Appl.,2015,429(1):415-438.

[37] GAO H,ZHAO H.Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization [J].Opt.Express，2010,18(3):1854-1871.

[38] GAO H,ZHAO H.Multilevel bioluminescence tomography based on radiative transfer equation Part 2: total variation and l1 data fidelity [J].Opt.Express,2010,18(3):2894-2912.

[39] GRASMAIR M,GROSSAUER H,Haltmeier M,etal..VariationalMethodsinImaging[M].New York: Springer,2009：294.

[40] NEATER W,SCHERZER O.Quantitative photoacoustic tomography with piecewise constant material parameters [J].SIAM.J.ImagingSci.,2014,7(3):1755-1774.

[41] BERETTA E,MUSZKIETA M,NAETAR W,etal..A variational method for quantitative photoacoustic tomography with piecewise constant coefficients [J].PreprintArXiv,2015:1509.04982.

[42] CEZARO A D,LEITO A,TAI X C.On piecewise constant level-set (pcls) methods for the identification of discontinuous parameters in ill-posed problems [J].InverseProbl.,2013,29(1):015003.

[43] JIANG H,YUAN Z,YIN L,etal..Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media [C].SPIE8thConferenceonBiomedicalThermoacoustics,Optoacoustics,andAcousto-optics.PhotonsPlusUltrasound:ImagingandSensing,SPIE-IntSoc.OpticalEngineering,SanJose,California,UnitedStates,2007.

[44] KALTENBACHER B,NEUBAUER A,SCHERZER O.Iterative regularization methods for nonlinear ill-posed problems [C].RadonSeriesonComputationalandAppliedMathematics，Berlin,2008.

[45] ELVETUN O L,NIELSEN B F.The split Bregman algorithm applied to PDE-constrained optimization problems with total variation regularization [J].Comput.Optim.Appl.,2016,64(3):699-724.

[46] AMMARI H,BOSSY E,JUGNON V,etal..Reconstruction of the optical absorption coefficient of a small absorber from the absorbed energy density [J].SIAMJ.Appl.Math.,2011,71(3):676-693.

[47] BERGOUNIOUX M,SCHWINDT E L.On the uniqueness and stability of an inverse problem in photo-acoustic tomography [J].J.Math.Anal.Appl.,2015,431(2):1138-1152.

[48] ZEMP R J.Quantitative photoacoustic tomography with multiple optical sources [J].Appl.Opt.,2010,49(18):3566-3572.

[49] BAL G,REN K.Multi-source quantitative photoacoustic tomography in a diffusive regime [J].Inverseprobl.,2011,27(7):75003-75022.

[50] BAL G,UHLMANN G.Inverse diffusion theory for photoacoustics [J].InverseProbl.,2010，26(8):25021-25032.

[51] SHAO P,HARRISON T,ZEMP R J.Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data [J].Biomed.Opt.Express,2012,3(12):3240-3249

[52] COX B T,ARRIDGE S R,BEARD P C.Estimating chromophore distributions from multiwavelength photoacoustic images [J].J.Opt.Soc.Am.A,2009,26(2):443-455.

[53] BAL G,JOLLIVET A,JUGNON V.Inverse transport theory of photoacoustics [J].InverseProbl.,2010,26(2):25011-25045.

[54] BAL G,REN K.On multi-spectral quantitative photoacoustic tomography in diffusive regime [J].InverseProbl.,2012,28(2):25010-25022.

[55] REN K,GAO H,ZHAO H.A hybrid reconstruction method for quantitative PAT [J].SIAMJ.ImagingSci.2013,6(1):32-55.

[56] BAUER A Q,NOTHDURFT R E,ERPELDING TN,etal..Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography [J].J.Biomed.Opt.,2011,16(9):096016.

[57] RAJIAN J R,CARSON P L,WANG X.Quantitative photoacoustic measurement of tissue optical absorption spectrum aided by an optical contrast agent [J].Opt.Express,2009,17(6):4879-4889.

[58] DaOUDI K,HUSSAIN A,HONDEBRINK E,etal..Correcting photoacoustic signals for fluence variations using acousto-optic modulation [J].Opt.Express,2012,20(13):14117-14129.

[59] DING T,REN K,VALLELIAN S.A one-step reconstruction algorithm for quantitative photoacoustic imaging [J].InverseProbl.,2015,31(9):65003-65023.

[60] HALTMEIER M,NEUMANN L,RABANSER S.Single-stage reconstruction algorithm for quantitative photoacoustic tomography [J].InverseProbl.,2015,31(6):65005-65028.

[61] YUAN Z,JIANG H B.A calibration-free,one-step method for quantitative photoacoustic tomography [J].Med.Phys.,2012,39(11):6895-6899.

[62] JIANG H B,ZHANG Q Z,YUAN Z,etal..Simultaneous reconstruction of acoustic and optical properties of heterogeneous media by quantitative photoacoustic tomography [J].Opt.Express,2006,14(15):6749-6754.

[63] RONDI L,SANTOSA F.Enhanced electrical impedance tomographyviathe mumford-shah functional [J].EsaimContr.Optim.,2000,6(6):517538.

[64] KIRSCH A.AnIntroductiontoTheMathematicalTheoryofInverseProblems[M].New York: Springer,1996:234

[65] 殷杰，陶超，刘晓峻.多参量光声成像及其在生物医学领域的应用 [J].物理学报,2015,64(9):098102.YIN J,TAO C,LIU X J.Multi-parameter photoacoustic imaging and its application in biomedicine [J].ActaPhys.Sinica,2015,64(9):098102.(in Chinese)

[66] 苗少峰,杨虹,黄远辉,等.光声成像研究进展 [J].中国光学,2015,8(5):699-713.MIAO S F,YANG H,HUANG Y H,etal..Research progresses of photoacoustic imaging [J].Chin.Opt.,2015,8(5):699-713.(in Chinese)

孙正(1977-)，女，河北保定人，博士，教授，硕士生导师，2004年于天津大学获得博士学位，主要从事生物医学信号处理、医学成像和图像处理等方面的研究。

E-mail: sunzheng_tju@163.com

Review on Progress of Quantitative Photoacoustic Tomography

SUN Zheng*,ZHENG Lan

(DepartmentofElectronicandCommunicationEngineering,NorthChinaElectricPowerUniversity,Baoding071003,China)

Photoacoustic tomography (PAT),an emerging medical imaging modality,combines the high resolution of ultrasonic imaging and high contrast of optical imaging.Current research on PAT includes two inverse problems,i.e.,constructing the distribution of initial acoustic pressures according to the photo-acoustic signals generated by the tissues and estimating the optical absorption and scattering coefficients of the tissues within the imaging region based on the results of the first inversion.The latter,known as quantitative photoacoustic tomography (qPAT),is in general a nonlinear ill-posed problem.This paper summarizes current algorithms for solving the qPAT inversion.Related advantages and limits as well as perspective studies are discussed.

quantitative photoacoustic tomography (qPAT); image reconstruction; inverse problem; optical absorption coefficient; diffusion coefficient; Gruneisen coefficient

1000-7032(2017)09-1222-11

2017-03-24；

2017-04-27

国家自然科学基金(61372042)； 中央高校基本科研业务费专项资金(2014ZD31)资助项目

O439

A

10.3788/fgxb20173809.1222

*CorrespondingAuthor,E-mail:sunzheng_tju@163.com

Supported by National Natural Science Foundation of China(61372042); Special Fund for Basic Scientific Research Business in Central Universities(2014ZD31)