基于改进固定点迭代方法的深层活体量化光声成像

柳芳艳 孟静 司广涛

摘 要：针对限制视图下光声图像的重建伪影问题，提出一种改进的固定点迭代量化光声成像方法。首先，通过传统的反投影重建算法重建由探测器探测到的原始光声压数据，得到原始的光声压图像;接着，利用自适应维纳滤波算法对原始的光声压图像进行滤波去除重建图像伪影;然后，通过光传输模型求解目标成像区域的光通量;最后，进行迭代计算，获得目标组织的光学吸收系数。此外，在求解光通量过程中引入Toast++软件来实现光传输模型的前向求解，

提高量化成像的效率和精确性。仿体和活体实验结果表明，与传统固定点迭代方法相比，所提方法能够获取更高质量的光声图像，重建得到的深层量化光声图像中存在较少伪影;量化重建的深层目标组织的光学吸收系数与浅层目标组织的光学吸收系数的数值非常接近，前者约为后者的70%，能实现深层生物组织光学吸收系数的定量重建。

关键词：量化光声成像;光学吸收系数;深层组织;维纳滤波;Toast++;光子传输模型

中图分类号：TP391.41

文献标志码：A

Abstract： Focusing on the reconstruction artifact of photoacoustic images in restricted view， an improved fixed-point iterative quantitative photoacoustic imaging method was proposed. Firstly， the original photoacoustic pressure data detected by the detector were reconstructed by the traditional back projection reconstruction algorithm to obtain the original photoacoustic pressure image. Secondly， the original photoacoustic pressure image was filtered to remove the reconstruction artifact by adaptive Wiener filtering algorithm. Thirdly， the optical transmission model was used to solve the optical flux of the target imaging region. And finally， iterative calculation was performed to obtain the optical absorption coefficient of the target tissue. In addition， Toast++ software was introduced in the process of solving the optical flux to realize the forward solution of the optical transmission model， which improved the efficiency and accuracy of quantitative imaging. The phantom and in vivo experiments show that compared with the traditional fixed-point iterative method， the proposed method can obtain photoacoustic images with higher quality and there are fewer artifacts in the deep quantitative photoacoustic images reconstructed by the method. The optical absorption coefficient of the quantitatively reconstructed deep target tissue is very close to the optical absorption coefficient of the shallow target tissue， the former is about 70% of the latter. As a result， the quantitative reconstruction of the optical absorption coefficient of the deep biological tissue can be implemented by the proposed method.

Key words： quantitative photoacoustic imaging; optical absorption coefficient; deep tissue; Wiener filtering; Toast++; photon transmission model

0 引言

光聲成像（PhotoAcoustic Imaging， PAI）是近几年发展起来的一种非入侵式和非电离式的新型生物医学成像方式[1]，它融合了光学和超声两种成像模式的特点，兼具纯光学成像的高对比度和纯超声成像的高分辨率优点，已成为目前发展最快的生物医学成像技术之一。光声成像作为一种重要的生物医学成像方式，为研究组织的结构形态，生理功能和新陈代谢指标提供了重要手段，在癌症和心脑血管疾病的早期诊断中具有重要的应用价值[2-3]。光声成像的基本原理是：生物组织经过合适的短脉冲激光照射后，吸收光能，组织内部温度瞬时升高，膨胀产生压力，这种压力以超声波的方式在组织内部传播，最后被放置在生物组织表面的超声换能器探测。根据采集到的超声信号，经过特定的重建方法可以得到生物组织光声图像[4]。

光声图像本质是一种光吸收累积图像，是组织吸收系数和光通量共同作用的结果。在实际应用中，获取组织的光吸收系数分布（称为量化光声成像）具有重要的价值。首先，吸收系数分布与生物组织中吸收发色团的空间分布直接相关，并且不依赖于光通量，这提供了更真实的图像对比度分布，便于深层信号的高对比成像;其次，直接从多波长光声图像获得的光谱成像是有问题的，因为在任何位置的波长依赖性实际上是用整个被照射的组织体积中的所有发色团的光谱特征进行编码的[5]。这种空间光谱串扰的效果会严重破坏感兴趣发色团的光谱特征，损害识别其存在并量化其浓度的能力。然而，若获取组织在多个波长下的吸收系数分布图像可避免上述情况。可见，量化光声成像对组织光谱成像具有重要的意义。结合光声光谱技术，可以定量分析组织中成分的变化，精细地反映极其微小的组织病变及血红蛋白浓度、血氧浓度、氧代谢率等一系列重要的生理参数[6-7]，这在癌症和心血管疾病诊疗方面具有广阔的应用前景。针对量化光声成像问题，国内外很多研究小组开展了大量研究。早期，Rinat O. Esenaliev曾经提出了Temporal Profile模型[8]，主要用于在光声成像中获取生物组织的吸收系数，并用于监测大脑中血氧浓度的变化。随后，Lihong. V. Wang研究小组提出了体外相似组织测量方法[9]，在体外测量经过小鼠脑皮和头后骨的光的分布情況，用来矫正光吸收分布，实现了在小鼠大脑中对氧合血红蛋白浓度的非侵入性成像。2016年，北京大学李长辉教授基于比尔定律矫正光通量，实现了活体人手关节量化成像[10]。以上方法虽然无需求解光传输方程便可得到高精度的活体量化光声成像，但是只适用于浅表组织或者特定深度的断层图像的量化，在临床应用中具有一定的局限性。从2006年开始，许多小组开始开展基于光传输模型的深层量化光声成像方法研究。该类研究方法可分为三类：1）基于混合光声和扩散光学成像系统（Diffuse Optical Tomography， DOT）的量化光声成像方式。2011年Adam Q. Baure与Lihong. V. Wang教授小组合作，应用时域DOT系统和基于超声阵列的光声层析成像系统（Photoacoustic Tomography， PAT），实现了仿体组织的深层量化重建[11]。2012年，Huabei Jiang研究小组设计了在同一个平台上的PAT/DOT双模成像系统，基于这个系统，该小组已经开展了仿体和体外乳腺癌组织成像的研究[12]。这种集成的DOT/PAT的双模成像量化方式，不需要知道组织背景先验信息，可以同时重建组织吸收和散射系数，但是需要两套成像系统的数据采集和信号重建，系统的复杂度高，计算量大。2）单步重建方式。单光声系统下，通过联合求解光声方程和光传输方程实现组织吸收系数重建。Huabei Jiang研究小组在假定背景散射系数和Gruneisien参数已知的情况下，应用有限元方法和传输方程，量化求解人手关节组织的吸收系数[13]。2015年，Sarah Vallelian小组提出了一种用于定量重建光声图像的一步算法，该算法可用于同时从多个声学数据集中恢复吸收系数和介质的超声速度，数值模拟实验说明了所提出重建方法的可行性[14]。同年，Simon Rabanser小组提出了使用一步重建算法进行定量光声层析成像[15]，通过将辐射传输方程与声波方程耦合，将定量光声层析成像的图像重建问题制定为单个非线性逆问题，数值仿真实验验证了算法的可行性。单步量化方法计算量大，实际应用复杂性高，目前大多只用于数值仿真和仿体实验。3）两步重建方式。先通过光声成像系统求解出光吸收累积图像;然后，根据光声系统的光源分布，基于光传输模型求解光通量;最后，通过一定的方法在原始光声图像中去除光通量，得到组织光吸收系数分布。2008年，Banerjee等[16]描述了一种非迭代方法，从边界压力测量值重建的吸收能量图中恢复光学吸收系数分布。该小组通过蒙特卡罗模拟获得边界声压测量值。2017年，剑桥大学的Brochu等[8]基于辐射传输方程，实施了校正光通量的迭代优化，并用于活体小鼠的吸收系数量化成像。伦敦大学的Cox和Bread小组[17]提出了一种简单的固定点迭代模型，在背景组织散射系数均匀且已知的情况下，该模型能够重构组织的光学吸收系数。该方法的局限是需要知道组织背景的光学信息，优点是简单、快速，在软组织量化光声方面具有良好的应用前景。2009年， Cox小组[18]基于固定点迭代模型开展多波长光谱成像方法，通过数值仿真实验，量化求解具有不同吸收系数的载色体浓度。2017年，北京大学李长辉教授小组基于手持式光声成像系统，结合固定点迭代方法，实现了人体乳腺组织的深层活体量化光声成像[19]。

可见，固定点迭代方法受到了研究者的关注。然而，基于该方法的深层活体量化光声成像研究刚刚起步，还有很多问题需要解决。其中，限制视图数据采集模式下重建伪影对深层量化成像结果有较大影响。因此，针对限制视图光声图像的重建伪影问题，本文提出了一种改进的固定点迭代量化光声成像方法。该方法是在迭代工作中对重建光声图像进行伪影滤波，然后进行深层组织量化成像。此外，本文在光传输模型前项求解方面融入了Toast++软件，进一步提高了量化光声成像的效率和准确度。实验中将量化重建算法用于一个新的128阵元商用超声阵列的仿体和活体成像，实验结果验证了本文方法的有效性和准确性。

1 相关研究

1.1 光声成像

根据光声生成理论研究，当激光脉冲的时间宽度足够短时，光声波在无损介质中产生和传播满足以下波动方程[20]：

其中：p0（r）表示由脉冲激光激发的初始声压;c是声速;r是空间位置坐标;t表示时间;表示光声压在位置r上的梯度;δ表示一个与时间成正相关的热激发函数。假设成像区域的组织在声学上是均匀的，可以通过求解以上波动方程来获取生物组织的光声压，因此，p（r，t）可以被表示为：

1.2 量化光声成像

根据光声成像的原理，由于生物组织被短脉冲激光照射，会瞬时使得局部组织的温度升高，在这种情况下，如果局部区域的生物组织通过吸收激光脉冲加热所产生的能量，就会产生声波。如果不考虑其他因素的影响，假设声波在生物组织中先行传播，则声压p的波动方程[17]是：其中：β是体积热膨胀系数;Cp是恒压比热容;Vs代表声速;H是单位时间内在每体积组织中沉积的热量;p和H通常取决于在组织中的位置r和时间t。

当激光脉冲足够短时，热函数可以被建模为以下与时间成正相关的δ函数[17]：

其中：A（r）是在位置r处转化为具体热量的吸收能量密度（即在r处具体的光吸收）。如果在生物组织的局部区域中，在位置r处的吸收系数分布是μa（r），并且在该位置的光通量是φ（r， μa），则在该位置的吸收能量密度A（r）[19]可以表示为：

根据光声效应机理，在短脉冲激光照射条件下，可以将式（3）进行化简，得到在组织激光脉冲照射后产生的声压，即初始的光声压力分布p0，其计算形式如式（6）所示：

其中：Γ是Grüeneisen系数，代表热量转换为压力的效率，可以取0～1的任意常数[21]。在本文中不考虑Γ的影响，假定所有位置的Γ相同，且取值为1（即认为所有的热量全部转换为探测器探测到的声压）。

根据式（5）和（6），可以得出在位置r处的初始压力p0（r）与吸收系数μa（r）和光通量φ（r）的关系式：通常，由于强烈的光散射和吸收，光通量φ（r）在组织中分布不均匀，重建的初始压力p0（r）不能准确地表示吸收系数μa（r）的分布。特别地，当随着组织深度的增大，φ（r）大幅降低时，这种不均匀的光通量导致深层组织中光声信号的对比度减小，影响深层信号成像。

1.3 光通量求解

量化光声成像中，最主要的困难是求解未知的光通量。虽然光在组织中的传播满足麦克斯韦方程，但是由于涉及空间尺度，在浑浊（高散射）介质中传播的光很快就达到了实际计算的极限。扩散方程是光传输模型的一种近似形式，它在高散射介质中可以高精度的表示光子在组织中的传输，扩散方程的数学公式[22]如下：

其中：k（r）=[3（μa（r）+μs（r））]-1是生物组织在某一位置的扩散系数;φ（r，ω）是求解得到的光通量分布; μa（r）是生物组织在某一位置的光学吸收系数， μs（r）是生物组织在某一位置的散射系数;c是光在介质中的传播速;ω是角频率;ζ是在边界处折射率不匹配项;v是外部边界标准;m是在边界处的某一位置;q是在Ω范围内边界上的光源分布; i是一个虚数单位。

扩散方程的求解可以通过有限元方法实现[23]。本文引入Toast++软件，实现了扩散方程的前向求解。首先，根据求解问题得到该问题的变分形式;其次，选择适合成像区域的子空间区域，并根据子空间区域将成像区域进行离散化，在本文中，采用了三角剖分形式;再次，基于Toast++软件中的dotSysmat（）函数，计算单元刚度并组合成总刚度矩阵K;最后，根据公式Kφ=q的变形公式φ=K-1q来求解光通量。

Toast++软件是一个开源软件，是由伦敦大学学院的Martin Schweiger研究小组[22]开发，主要用于解决DOT成像中的正向和逆向问题。该软件包括一组函数库，用于模拟具有复杂边界和不均匀内部参数分布的高散射介质中的近红外光传播，可以对稳态、时域和频域数据采集系统进行建模。在使用过程中，可以直接链接到Toast ++库以编译DOT的应用程序，并且可以快速地调用Toast++函数库中的函数。

2 重建方法与实验设置

2.1 固定点迭代方法

传统的固定点迭代重建算法是一个简单的迭代机制，适用于在重建过程中从生物组织的吸收能量中去除光通量对恢复生物组织的光学系数（吸收系数）的影响，从而根据假设的生物组织的光学系数来重建出生物组织真实的光学吸收系数的分布，得到生物组织真实的光学图像[17]。

传统的固定点迭代算法主要分为以下几个步骤：

1）首先设置初始的光学吸收系数μa=0;

2）根据DOT模型求解光通量分布φ;

3）根據光通量和光学吸收系数计算光声压强;

4）计算原始光声压强与重建光声压强之间的误差;

5）迭代计算μa和φ，直到光声压强误差小于设定阈值。

2.2 改进的固定点迭代方法

传统的固定点迭代方法重建光学吸收系数分布，计算过程简单并且能够得到较为精确的重建结果;但在限制视图光声成像系统中，深层组织的重建光声图像存在较多的重建伪影，严重影响量化光声图像的精度和质量。因此，针对这个问题，本文提出了一种改进的固定点迭代量化光声成像方法，该方法融合了Toast++软件，具体流程如图1所示。

整个流程主要包括三个阶段：

1）重建原始光声图像。在这个过程中，首先，利用传统的反投影重建算法重建经探测器探测到的原始光声压数据得到原始的光声压图像。然后，采用2-D自适应维纳滤波去除原始光声压图像中的噪声和伪影。自适应维纳滤波是一种使噪声图像和其滤波图像之间的均方误差最小的复原方法，并根据局部方差的大小调整滤波效果[24]。在该方法中，使用大小为m×n的邻域（模板）估计局部图像均值和标准偏差。如果邻域选择太小，则只能对噪声滤除产生很小的作用;如果邻域选择太大，则会导致滤波后的图像失真。因此，在本文中通过选择几组不同邻域大小的模板进行实验，最终选取12×12的模板作为维纳滤波在光声图像去噪应用中的邻域大小。

2）光通量求解。在这个过程中，在组织光学散射系数已知的情况下，首先利用Toast++软件中的toastMesh（）函数结合有限元算法对成像区域进行网格剖分;其次，设置初始的光学吸收系数，在本文中，将目标组织的光学系数假设为背景组织的光学吸收系数;再次，根据toast++软件中的dotSysmat（）函数，输入剖分好的成像区域网格，设置的初始光学参数，系统频率等参数，计算每个单元刚度并组合成总刚度矩阵K;最后，根据公式φ=K-1q得出第k次迭代中的光通量φ（k）。

3）更新光学参数分布。首先，根据1.2节中的式（7）求解第k次迭代中的光声压强p（k）c;其次，迭代更新，得到光学吸收系数μ（k+1）a，直到滤波后的光声压图像与计算光声压图像之间的误差小于设定的阈值ε。

单阵元大小 300μm×900μm;3）128-通道数据采集卡。在实验中采用单侧打光的照明方式，用一根1500μm的多模光纤打光，光经过两个透镜在成像组织表面形成一个大小为38mm×5mm的矩形光斑[25]。为实现超声信号的耦合探测，系统准备了装满去离子水的水槽，水槽下部有一个由低密度聚乙烯薄膜封闭的窗口，用于激光照射和信号采集。基于该系统，每一帧二维图像的数据采集需1次激光脉冲照射，三维数据可以通过机械扫描获得。

2.3.2 仿体实验与活体实验

仿体模型为在琼脂块中插入针灸针。实验中，琼脂块中Agar的浓度为1.5g/100ml。在制作过程中，在长方体容器中加入了800ml水，并加入12g Agar，加热到100℃，然后进行自然冷却凝固，凝固后琼脂块的高度为72mm、宽度为80mm，将四根直径大约为0.3mm的针灸针在不同深度斜着嵌入到长方体琼脂块模型中，如图2所示（x是换能器阵列的（横向）方向，y是机械扫描方向，z是深度方向）。图2中四个实心黑点表示要成像的四根针灸针。采用波长为780nm的激光照射扦插钢针的仿体表面，配合机械扫描采集光声信号。

对于活体成像，将激光输出调节至1064nm（一个等吸光点），在该点处氧合血红蛋白和脱氧血红蛋白均匀地吸收光。对于小动物成像，使用Sprague Dawley大鼠（雌性）。

在麻醉的过程中，实施腹腔内注射氯胺酮（85mg/kg）和甲苯噻嗪（15mg/kg）混合物来深度麻醉大鼠，等待一段时间确认大鼠已处于麻醉状态（呼吸、心跳平稳，按压脚趾无疼痛反射）。在实验过程中，首先利用温和脱毛膏除去大鼠腹部皮肤表面上的毛发（为了防止表皮毛发对实验结果的干扰），然后在大鼠腹部表皮上面涂抹一层薄薄的超声耦合胶（避免探测声压直接在空气中传播，提高耦合效率），用动物支架将大鼠固定到成像系统中，并给大鼠通上氧气，在大鼠身下配备加热垫，使其体温始终保持在37±0.5℃。将大鼠腹部皮肤表面的光通量设定为每脉冲约～19mJ/cm2。实验持续1h，在整个实验过程中，本研究确保应用此大小的脉冲能量对组织没有热损伤。实验结束后，动物自然恢复体能，继续放回笼中饲养。本文描述的动物均在中国科学院深圳先进技术研究院和曲阜师范大学动物实验伦理动物委员会批准的顺应性方案中进行。

3 实验结果与分析

3.1 仿体实验结果

图3展示了滤波后的仿体光声图像以及重建后的光学吸收系数分布图像，在本实验中只展示了一帧B-scan仿体实验数据。由于780nm波长下，水是主要的光学吸收物质，琼脂本身既没有散射也没有吸收。因此，在本文实验中，假设琼脂块的吸收系数与散射系数在此波长下与水的吸收和散射系数相同，分别设为μa=0.02346cm-1， μs=1.0cm-1。通过将滤波后的光声图像与光通量图分开来实现光通量补偿从而得到目标物体的光学吸收系数分布图像。图3（a）是经过滤波后的光声压图像，从中可以看出：由于光通量分布不均匀，导致在相同光学系数下获得的位于不同深度的四根针灸针的光学对比度不同;处于较深位置的两根针灸针的对比度较差（图中白色圆圈所示）。通过传统的固定点迭代算法与本文算法，量化重建得到了具有光通量补偿的光学吸收系数图像如图3（b）、（c）所示。图3（b）、（c）是利用传统的固定点迭代算法和本文算法量化得到的光学吸收系数分布图像。从图3（c）标记的两个位置可以看出，在较深位置的两个针灸针与较浅位置的对比度基本相同，也就是说，经过本文算法也可以定量地恢复深层仿体目标的光学吸收系数，获得目标物体相对真实的光学吸收分布图像，实现对深层组织的量化;并且在重建光学吸收系数分布的过程中，通过对原始的光声压信号进行滤波，重建得到的光学吸收系数分布中存在较少的噪声和伪影（图3（b）和（c）中白色方框所示）。

为了定量评估本文算法的效果，用图3（d）、（e）分别展示了位于四个不同位置的滤波后的光声压数值分布与重建得到的吸收系数数值分布。这四个位置分别在图3（a）中用数字标记。从图3（d）中可以清楚地看出，在位置2、3、4这三个位置点上，随着组织深度加深，光声压的数值逐渐减小，并且数值之间差别较大。但是经过本文算法所重建得到的吸收系数之间的数值差别较小（图3（e）所示）。可见，本文算法能有效补偿光通量，能够定量恢复组织的光学吸收系数分布。理论上，相同的物质具有相同的光学吸收系数，针灸针属于钢材料，钢是铁碳合金，主要由铁元素和碳元素构成，且铁元素含量较高，因此在本文中把在780nm波长下铁的吸收系数0.6cm-1当成量化标准。通过图3（e）可以看出虽然本文算法能够在这四个位置上得到较为接近的光学吸收系数，但是每个位置上的吸收系数值仍然存在误差。这是因为在通量补偿的过程中，目前还没有一种有效的算法能够进行精确的通量补偿，本文研究的算法在进行通量补偿时也存在误差，因此本文算法也只是一种相对量化，重建得到的是目标物体光学吸收系数的相对值。

3.2 活体实验结果

由于皮肤下的主要光声信号来源是血管，所以在本文研究中假设除了血管以外的其他背景组织的光学系数都为软组织的光学系数。由于1064nm波长下，软组织中水是主要的吸收物质，因此，组织背景的光学吸收系数设置为1064nm波长下水的光学吸收系数，即μa=0.12cm-1[19]。該波长下软组织的散射系数为μs=7.1563cm-1[26]。

利用传统的固定点迭代算法与本文算法量化重建的实验结果比较如图4所示。图4中展示了四帧不同的B-scan活体实验数据。其中图4（a）是滤波后的光声压图像。在实验结果展示图中，本文标记了两个不同深度的血管（图4中白色虚线圆圈所示）。如前所述，随着深度的增加，由于光通量的影响，重建的光声信号对比度变得越来越弱，当达到一定深度，血管基本上被淹没在背景中。尽管通过传统的固定点迭代算法能够恢复目标组织的光学吸收系数（如图4（b）所示），但是量化重建的结果中存在较多伪影。通过图4（c）可以看出，本文算法能够有效抑制光通量的影响，得到生物组织相对真实的光学吸收分布图像。进一步说明，本文方法可以实现对深层组织的光学系数的量化。通过比较图4（b）和（c）白色虚线方框中的内容，可以看出，当使用改进的固定点迭代方法时，背景中的大部分伪影被移除，提高了图像的保真度，获得较高质量的光学吸收系数图像，实现更加精确的量化重建。

同樣为了评估本文方法的优势，在图5中分别展示了经过滤波后的这四帧不同B-Scan的光声压数值分布与重建得到的吸收系数值分布。其中，主要展示了位于四个不同深度的信号点的数值分布情况，这四个位置分别在图4（a）中用数字标记。图5（a）是经过滤波后的光声压数值分布，图5（b）是利用本文算法量化重建后获得的光学吸收系数数值分布。从数值分布情况可以看出，随着组织深度加深，位于不同位置的光声压的数值之间差别较大（如图5（a）所示），尤其是当目标血管处于较深位置时，光声压值几乎变得非常小。但是通过改进的固定算法重建得到的目标血管的光学吸收系数之间的数值差别较小（如图5（b）所示）。当目标血管的位置较深时，重建获得的光学吸收系数基本上与较浅位置血管的吸收系数相同。可见，通过活体实验，说明本文算法能有效补偿光通量并能够有效地恢复深层生物组织的光学吸收系数分布，在量化深层组织光学吸收系数具有一定的优势。

4 收敛性分析

在本文实验中还计算了当前重建的光声压值与滤波后的光声压值之间的误差。将从当前重建得到的光声压与滤波后的光声压图之间的差异与迭代次数用作迭代停止标准，如2.2节中图1所示。在本文实验中，当重建得到的光声压与滤波后的光声压图之间的差异变得非常小（即小于本文实验所设计的阈值），并且当迭代次数达到最大迭代次数时，迭代过程将停止。

为了验证本文方法的收敛性，针对迭代过程绘制了两个不同B-scan光声图像的滤波后的光声压与重建光声压随迭代次数增加的误差曲线。如图6所示：图（a）是仿体目标滤波后的光声压与重建光声压随迭代次数变化的误差收敛曲线;

图（b）是大鼠腹部血管误差收敛曲线。误差值在前两次迭代中迅速减小，大约在第5次迭代后收敛。

原始声压与重建光声压之间的差异接近于零，表明本文方法能够快速收敛。

5 结语

针对限制视图下光声图像的重建伪影问题，在本文中主要提出了一种改进的固定点迭代量化光声成像方法。对于复杂的深层生物组织，本文方法能够有效抑制量化重建组织光学吸收系数图像中的噪声和伪影，获取高质量的组织光学吸收系数图像。相对于传统的固定点迭代方法，本文方法进一步提高了量化光声成像的效率和准确度。

但是，本文的研究工作只对量化光声成像方法进行了初步的研究，仍然还存在未解决的问题。在今后的工作中，我们将对3-D光声图像进行量化重建，并且也会在本文算法中引入更高级的重建伪影去除方法，进一步提高深层光声图像的量化重建质量。

参考文献（References）

[1] LIN L， HU P， SHI J， et al. Single-breath-hold photoacoustic computed tomography of the breast [J]. Nature Communications， 2018， 9（1）：2352.

[2] FERNANDES D A， KOLIOS M C. Intrinsically absorbing photoacoustic and ultrasound contrast Agents for cancer therapy and imaging [J]. Nanotechnology， 2018， 29（50）： 1361-6528.

[3] MATSUMOTO Y， ASAO Y， SEKIGUCHI H， et al. Visualising peripheral arterioles and venules through high-resolution and large-area photoacoustic imaging， Scientific Reports 8[R]. Nature， 2018： Article number： 14930.

[4] BARDSLEY P， REN K， ZHANG R. Quantitative photoacoustic imaging of two-photon absorption [J]. Journal of Biomedical Optics， 2018， 23（1）：1-11.

[5] LEE C H， FOLZ J， ZHANG W， et al. Ion-selective nanosensor for photoacoustic and fluorescence imaging of potassium [J]. Analytical Chemistry， 2017，89（24）：7943-7949.

[6] DAHAL S， CULLUM B M. Characterization of multiphoton photoacoustic spectroscopy for subsurface brain tissue diagnosis and imaging [J]. Journal of Biomedical Optics， 2016， 21（4）：No.47001.

[7] LIM H T， MURUKESHAN V M. Hyperspectral photoacoustic spectroscopy of highly-absorbing samples for diagnostic ocular imaging applications [J]. International Journal of Optomechatronics， 2017， 11（1）：36-46.

[8] BROCHU F M， BRUNKER J， JOSEPH J， et al. Towards quantitative evaluation of tissue absorption coefficients using light fluence correction in optoacoustic tomography[J]. IEEE Transactions on Medical Imaging， 2017， 36（1）：322-331.

[21] ZHANG J， ZHANG J， LI W. Degradable hollow mesoporous silicon/carbon nanoparticles for photoacoustic imaging-guided highly effective chemo-thermal tumor therapy in vitro and in vivo [J]. Theranostics， 2017， 7（12）： 3007-3020.

[22] SCHWEIGER M， ARRIDGE S. The Toast++ software suite for forward and inverse modeling in optical tomography [J]. Journal of Biomedical Optics， 2014， 19（4）： 040801.

[23] SHAN T， QI J， JIANG M， et al. GPU-based acceleration and mesh optimization of finite-element-method-based quantitative photoacoustic tomography： a step towards clinical applications [J]. Applied Optics， 2017， 56（15）： 4426-4432.

[24] 易三莉，賀建峰.基于BEMD与自适应维纳滤波的图像降噪[J].计算机工程与应用，2013，49（10）：156-158.

（YI S L， HE J F. Image denoising method based on BEMD and adaptive Wiener filter [J]. Computer Engineering and Applications， 2013， 49（10）： 156-158.）

[25] LI M， LIU C， GONG X， et al. Linear array-based real-time photoacoustic imaging system with a compact coaxial excitation handheld probe for noninvasive sentinel lymph node mapping [J]. Biomedical Optics Express， 2018， 9（4）： 1408-1422.

[26] JACQUES S L. Optical properties of biological tissues： a review [J]. Physics in Medicine and Biology， 2013， 58（11）： R37-61.