孟 琪,孙 正*
(1.华北电力大学电子与通信工程系,河北 保定 071003;2.华北电力大学河北省电力物联网技术重点实验室,河北 保定 071003)
光声层 析成像(Photoacoustic Tomography,PAT)是一种新型的生物医学成像方法,它结合了光学成像的高对比度和超声成像的高分辨率[1],可应用于眼科[2]、胃肠疾病[3]、脑科学[4]、心脏病学[5]和皮肤病学等。其物理基础是生物组织的光声效应,用短脉冲激光照射组织,组织吸收部分光能量后受热膨胀,产生宽带超声波(即光声信号),并迅速向组织表面传播。超声换能器采集光声信号后送入计算机内,通过求解声学逆问题,从声压测量数据中重建出组织表面及内部的初始声压分布图或光吸收能量分布图,以显示组织的形态结构[6]。标准重建方法有延迟求和法[7]、反投影法[8]、时间反演法[9]和傅立叶变换法[10]等。在此基础上,通过求解光学逆问题,还可以定量估算组织的光学特性参数(如光吸收系数和散射系数)及热膨胀系数等的空间分布,实现功能成像,即定量PAT(quantitative Photoacoustic Tomography,qPAT)[11-12]。
为了简化问题,光声成像算法通常采用简单的光子传输模式,或假设入射光在组织表面和内部分布均匀。但在实际应用中,大多数生物组织都是具有高散射和不透明特点的混浊介质,当光在组织内传输时既被吸收又被散射。不同成分的组织结构对光的吸收和散射以及光通量的分布都是不均匀的。此外,随着光在组织中入射深度的不断增加,光的穿透能力随之减弱,光能量的分布随之分散,光通量也会显著降低。除上述原因以外,不合理的照明模式也会造成组织内的光照分布不均匀,只有靠近入射光的部分组织可被相对均匀的照明,进而导致光通量分布不均匀,这不仅会降低重建图像的精度和质量,也会降低成像深度[13]。本文对目前解决PAT 不均匀和不稳定照明问题的主要方法进行归纳和总结,主要包括补偿光通量的变化、精准估计光通量分布和优化照明模式3 种解决方案。
早期研究中,模拟光在组织中的传播时,一般利用假设的光通量分布消除不均匀照明对图像重建精度的影响[14],其主要困难在于需要已知有关组织复杂光学特性的先验知识。目前,补偿光通量变化的主要方法有稀疏分解法、基于成像模型的方法、基于荧光蛋白的方法、光声-声光光谱组合法和光声-超声双模态成像融合法,如表1所示。
表1 补偿PAT 光通量变化的主要方法对比Tab.1 Comparison of main methods of compensating for variations in light fluence in PAT
基于稀疏表示的分解方法是将PAT 图像分解为两个分量:由组织中扩散的光通量引起的低空间频率全局分量和代表光吸收系数变化的高空间频率局域分量。再通过少量基函数(通称为库,通常采用二维离散Haar 小波基和二维离散傅立叶基)之和逼近图像[15]。
根据光吸收能量分布与组织光吸收系数和光辐射通量之间的关系,求得组织光吸收能量分布的稀疏表示,利用正交匹配追踪(Orthogonal Matching Pursuit,OMP)算法找到满足稀疏表达式的最少数量的系数,即实现图像分解,进而计算出光吸收系数和光通量对数的向量表示。
实验结果表明,采用该方法可以有效提高根据光吸收能量分布图重建光吸收系数分布图的精度,如图1(彩图见期刊电子版)所示。该方法不必求解光辐射传输方程,不依赖于选择一个准确的光子传播模型来解释光通量的不均匀性,没有将建模误差传递到光吸收系数的估算结果中,具有较高的量化精度和稳定性。但是,它不适用于具有缓慢变化的光吸收系数的组织成像。此外,当入射光只从成像目标一侧照射时,每个投影都有不同的光照模式,因此,非均匀边界光照也会在光吸收系数的重建结果中引入一定误差,即作为算法输入的光吸收能量分布图不能准确代表光辐射通量与光吸收系数的乘积。为了提高算法的准确性,可以增加库中的元素,寻找更全面的库表示光吸收系数和散射系数,或者使用更精准的分解算法,例如基跟踪。
图1 采用稀疏分解法的PAT 图像重建结果[15]。(a)仿体的几何结构图;(b)光吸收能量分布图;(c)光通量分布图;(d)光吸收系数分布图Fig.1 Results of PAT image reconstruction by using the sparse decomposition algorithm[15].(a)Geometry of the phantom to be imaged;(b)optical absorption distribution;(c)light fluence distribution;(d)optical absorption coefficient
对根据光声信号测量值重建的初始声压分布图进行光通量补偿,可得到光吸收能量分布图,但这个过程可能放大噪声和伪影,特别是在有限视角测量条件下,会造成图像质量下降。针对这一问题,Bu 等[16]提出了一种重建光吸收系数与补偿光通量相结合的方法,将光通量补偿集成到基于模型的图像重建中,使每个测量位置的光通量和光吸收系数分离开,进而估算光吸收系数。
在组织具有均匀声学特性和非黏性介质的前提下,利用格林函数法[17]求解光声波动方程,得到计算光声信号理论值的矩阵表达式,通过使光声信号测量值和理论值之间的L2 范数最小,定量估计光吸收系数的空间分布。该方法引入了矩阵压缩算法[18],可以将矩阵压缩到原大小的1/250,不仅减少了所需的内存,而且加快了重建速度。与常规qPAT 方法相比,该方法在每次迭代中补偿光通量,减少了由于光通量变化所致的重建误差,提高重建图像质量。此外,该方法还具有基于模型的PAT 图像重建方法的优点,例如可减少由于有限角度测量和声速分布不均匀所致的图像伪影和分辨率下降等问题。
该方法的局限性包括:没有考虑超声探测器性能对重建图像质量的影响,超声换能器矩阵的中心频率和带宽较低,导致没有很好地抑制噪声。另外,即使采用了矩阵压缩算法,仍然需要较大内存,使得该方法很难应用于高分辨率、大规模、实时的PAT 图像重建。未来可以考虑利用压缩感知方法进一步提高重建速度,以及采用高质量的优化方法解决大规模优化问题。
基于模型的方法和稀疏分解法均适用于静态成像物体。对于运动器官成像中动态光通量的校正,Deán-Ben 等[19]提出了非时间混合的多光谱光声层析成像(temporally unmixed MultiSpectral Optoacoustic Tomography,tuMSOT)技术,利用可逆光开关荧光蛋白(Reversibly Switchable Fluorescent Protein,RSFP),结合实时采集的光声光谱数据校正散射介质内的光通量分布。RSFP 可以在特定波长下实现暗-亮状态的切换,从而使组织的光物理性质发生变化,产生光声信号。利用荧光分子探针的暗-亮状态之间的转换时间校正组织中光通量分布的不均匀性(图2)。即使成像目标存在运动(如动脉血管或胃肠道),采用该方法仍可实现目标区域中光通量的有效校正。当暗-亮状态切换足够快时,荧光分子探针的可逆性还可用于校正纵向光通量。利用具有波长开光性的光致变色分子,可以在不同波长下校正光通量。
图2 充分激活RSFP 前后的三维光声图像[19]。(a)激活前;(b)激活后Fig.2 Three-dimensional photoacoustic images(a)before and(b)after full activation of RSFPs[19].Reprinted with permission from©The Optical Society.
未来可以通过吸收光谱分析荧光分子探针的空间分布,使RSFP 的分布测量具有更高的时间分辨率。此外,可进一步利用基因编码的蛋白质和人工合成的光致变色分子,实现定量tuMSOT。
基于RSFP 的方法验证了在不同入射光波长下校正光通量的可能性,此后Hussain 等[20]提出了一种光声光谱(Photoacoustic Spectroscopy,PAS)成像中校正不同波长光通量的方法,将光声与声光(Acousto-Optics,AO)成像相结合,在无法预知介质光学特性的情况下,补偿与空间和波长相关的光通量变化。该方法涉及两个光声测量和一个声光测量,用光声和声光信号测量数据来表示光吸收系数,实验装置如图3 所示。
图3 PA-AO 联合成像实验装置原理图[20]Fig.3 Schematic diagram of experimental setup of PA-AO joint imaging[20].Reprinted with permission from©The Optical Society.
声光测量是基于散斑对比度的检测方法[21,22],且在已知入射光波长的情况下,散斑对比度的变化与超声标记光功率成正比,据此可以补偿光通量的变化,如图4(彩图见期刊电子版)所示。
图4 采用PA-AO 光谱组合法得到的光声图像[20]。(a)λ=755 nm 时从side1(左图)和side2(右图)照射介质得到的光声图像;(b)λ=780 nm 时从side1(左图)和side2(右图)照射介质得到的光声图像;(c)λ=755 nm(左图)和λ=780 nm(右图)时补偿光通量后的光声图像Fig.4 PA images obtained by using PA-AO spectral combination method[20].(a)PA images by exciting the medium from side 1(Left)and side 2(Right)whenλ=755 nm;(b)PA images by exciting the medium from side 1(Left)and side 2(Right)whenλ=780 nm;(c)PA images ofλ=755 nm(Left)andλ=780 nm(Right)after light fluence compensation.Reprinted with permission from©The Optical Society.
与基于PAS 的方法相比,该方法不依赖于任何有关组织光衰减特性和光传输模型的假设,而是通过声光测量消除成像质量对组织光学特性的依赖。但是它忽略了由局部散射变化引起的超声光调制效率的变化,其是否影响成像精度还有待进一步研究。声光成像中超声标记体积的大小也会对光吸收系数和散射系数的估计精度产生影响,未来可以考虑利用纳秒脉冲激光进行声光测量。
超声换能器可同时工作在脉冲/回声模式和接收模式,在超声成像装置的基础上增加光声成像单元,二者共用同一个超声换能器和信号采集卡,可形成光声-超声(photoacoustic-ultrasonic,PAUS)双模态成像,同时对组织的光学特性和声阻抗特性进行成像。
对于手持式乳腺PAUS 双模态成像,Zhao 等[23]提出了一种解决深层组织中光通量急剧下降问题的方法,即首先对乳腺进行超声成像,然后根据乳腺中不同组织结构的光学参数,利用有限元法求解光扩散方程计算光辐射通量,最后利用光通量分布图对原始光声图像进行光通量补偿,实验装置如图5(彩图见期刊电子版)所示。
图5 手持式乳腺PAUS 成像装置示意图[23]。(a)成像装置示意图;(b)中央切割平面的光通量Fig.5 Schematic diagram of handheld breast PAUS setup[23].(a)Imaging setup;(b)map of light fluence in the central cut plane
该方法可以改善乳腺PAT 图像的质量,提高深层组织的对比度,但也造成了背景信号的增强(包括深层区域的噪声),从而导致深层组织成像的信噪比严重下降,如图6 所示。
图6 补偿光通量前后的乳腺图像[23]。(a)未补偿光通量的PAT 图像;(b)补偿光通量后的PAT 图像;(c)未补偿光通量的PAUS 图像;(d)补偿光通量的PAUS图像Fig.6 Breast images before and after compensation for light fluence[23].PAT image before(a)and after(b)light fluence compensation;PAUS dual-modal image before(c)and after(d)light fluence compensation
此外,Jin 等[24-25]提出一种基于光声-被动超声双模态成像的光通量补偿方法,超声换能器同时接收组织产生的光声波和被动超声波,得到光声和漫反射(Diffuse Reflectance,DR)双模态图像,再采用反射解耦的方法对两幅图像进行融合,补偿光通量变化,如图7 所示。
图7 光声-被动超声融合成像原理示意图[24]Fig.7 Schematic diagram of photoacoustic-passive ultrasonic fusion imaging[24]
相比于单模态PAT,采用光声-被动超声双模态成像方法可有效减少由于反向散射和漫反射引起的估算误差,改善重建图像质量,如图8(彩图见期刊电子版)所示。但是,该方法的不足是:校准函数相对简单,不适用于所有应用场合;散射层边界处的漫反射图像可能发生畸变;系统分辨率不高;组织中的气泡不满足校准方程,因而无法识别和消除。未来可采用强聚焦激光束(光学分辨率光声显微镜)[26]或聚焦传感器(声学分辨率光声显微镜)[27]实现高分辨率融合成像。
图8 光声与光声-被动超声融合图像对比[24]。(a)光声图像;(b)光声-被动超声融合图像Fig.8 Comparison of(a)photoacoustic images and(b)photoacoustic-passive ultrasonic fusion images[24]
在qPAT 技术中,重建组织光吸收系数分布图时,传统方法通常是将图像重建与光传输模型相结合,利用迭代法求得光吸收系数的最小二乘解。此类方法的估算结果中可能包含与光通量相关的误差,降低成像精度。目前,精准估计光通量的主要方法有扩散光层析成像(Diffuse Optical Tomography,DOT)法、光声-声光信号组合法和表面光增强的方法,如表2 所示。
表2 精准估计PAT 光通量的主要方法对比Tab.2 Comparison of main methods for accurately estimating light fluence in PAT
将PAT 与DOT 相结合,利用低分辨率DOT重构组织的光吸收系数和散射系数的空间分布,进而采用有限差分法求解散射光子密度方程,得到组织表面和内部的连续波光通量,并在此基础上准确重建光吸收系数的空间分布,提高光声图像重建精度和光声光谱分析的精度,校正光声图像中与光通量相关的误差[28-29],如图9(彩图见期刊电子版)所示。
DOT 采用光学方法测量光通量,是一种波长相关的方法,虽然得到的光学特性图是平均量,但仍可用于准确估计光通量,而且可以重建低于DOT分辨率的目标。但其不足之处是由于假设被测组织的Gruneisen 系数是常数,因此只适用于具有相同成分或热力学特性差异不明显的成像目标。同时该方法会增加成像系统及其操作的复杂度,限制了其应用范围,未来仍需加快活体PAT-DOT集成系统的开发。
图9 PAT-DOT 的成像结果[28]。(a)DOT 测量的光吸收系数和散射系数分布图;(b)组织表面的光通量分布图;(c)补偿光通量之前的PAT 图像;(d)补偿光通量之后的PAT 图像Fig.9 Results of the PAT-DOT method[28].(a)The distributions of optical absorption coefficient and scattering coefficient measured by DOT;(b)light fluence on the phantom surface;PAT images(c)before and(d)after compensating for light fluence
将光声和声光信号测量相结合,根据光子传输概率重建组织的光吸收系数和发色团浓度,可提高qPAT 的精度[30-34]。该方法基于两个原理:一个光子通过散射介质时可以按照相等概率沿两个方向运动;采用聚焦超声波可以在介质的一定体积内标记穿过的光子。
该方法在不预先确定介质的局部光通量和光学特性的情况下,仅通过实验测量的方式,利用测量值和实验参数估算局部光吸收系数(图10),可在很大程度上减小光吸收系数估计值与真实值之间的误差。但是光子功率是一个无法直接测量的内部量,需要根据从介质外测量的光声信号重建得到,与直接测量相比,此过程引入了一定的误差。此外,该方法假设光子的平均自由程大于标记体的线性维数,然而实际应用中该假设并不适用于所有光子。未来研究中需重点减小简单离散标记与超声标记之间的差距,并通过活体实验验证方法的可行性。
图10 光子在注入点1、标记体积2 和检测器位置3 之间的轨迹示意图[30]Fig.10 Schematic of photon trajectories between injection point 1,labeling volume 2 and detector position 3[30]
表面光增强的方法是利用扩散近似作为光传输模型,求解使模型输出值(即初始声压理论值)和测量值之间误差最小的最小二乘问题,利用初始声压和表面光测量数据同时估计组织的光吸收系数、散射系数和Gruneisen 系数[35]。与常规qPAT方法相比,该方法可以同时估算出光吸收系数、散射系数和Gruneisen 系数的空间分布,并且提高光吸收系数和散射系数的估算精度,如图11(彩图见期刊电子版)所示。但是最小化是否具有唯一性尚无法确定,而且文献[35]中仅进行了仿真实验验证,需要采用实际成像系统采集真实成像数据进行进一步的可行性验证。
图11 采用表面光增强法的PAY 图像重建结果[35]Fig.11 Results of PAT image reconstruction with the augmented PAT[35].Reprinted with permission from©The Optical Society.
对于PAT,采用不同时间、不同角度和不同位置的照明都会获得不同的光吸收系数分布图。通过改进照明模式,可以改善不均匀和不稳定照明对成像质量的影响,保证成像目标内各个视角都有足够的光穿透性。目前的改进方法可分为外部照明和内部照明两类,如表3 所示。
表3 PAT 照明模式的比较Tab.3 Comparison of two schemes of illumination in PAT
4.1.1 旋转照明
如图12 所示,旋转照明是将成像目标固定在成像系统的中心位置,“I”型光纤照明条产生脉冲激光并围绕目标旋转,照明条与传感器固定于旋转支架,并围绕扫描轴一起旋转,获得断层数据集[36]。一次完整的扫描是由围绕中心旋转512 步的照明条和探测器完成的,每旋转一步,探测器记录声压数据,进而获得目标的三维光吸收能量分布图。
与顶部照明和多侧面照明模式相比,旋转照明不是一次性照亮整个目标,而是通过旋转将入射光覆盖于整个目标,因此具有更高的光穿透性。但同时也在系统方程中引入了断层扫描数据的不一致性,即在每个层析视角记录的声压数据由不同的目标函数产生,目标内部的光通量分布随视角的变化而变化。对于较大的成像目标,该方法采集的数据不一致性程度更高,可能造成重建图像质量的下降。未来可以进一步研究旋转照明与内部照明相结合的方法,以解决大型物体成像深度下降的问题。
图12 采用旋转照明的乳房三维PAT 系统示意图[36]。(a)PAT 系统;(b)乳房假体的二维横截面Fig.12 An illustration of 3D breast PAT system employing the rotating partial illumination design[36].(a)Imaging system;(b)2D cross-section of the numerical breast phantom
4.1.2 非平稳照明
非平稳照明是对旋转照明的优化,将照明条分布在一个环绕成像目标的圆环上,如图13 所示,将入射光投向目标,通过改变照明条的数量与排列方式,得出改进数据不一致性的照明条件[37]。实验结果表明,当照明条数量增加即照明充足且检测面封闭时,可以有效降低测量数据的不一致性[37]。
图13 非平稳照明示意图[37]Fig.13 Schematic diagram of non-stationary illumination[37].Reprinted with permission from©The Optical Society.
该方法的缺点是在组织内部仍存在由于光束的不规则重叠所致的不均匀照明。为了弥补这种不足,Park 等[38]提出利用三维PAT 系统扫描整个乳房,根据比尔朗伯定律估计非均匀照明的极角函数和与深度有关的光衰减并进行补偿,直至达到临界深度。与无补偿成像方式相比,成像深度提高了67%,是目前无压缩三维乳房PAT 中成像深度最高的。
4.1.3 光捕捉器增强照明
组织表面对入射光的反射造成大部分光能量的损失,这种光损耗直接导致光声信号强度的减弱和重建图像质量的下降。Yu 等[39-40]设计了一种针对表面光反射的新型照明增强装置,利用3D 打印技术制作高频光声成像探头,探头内的光捕捉器是一种简单的凹面光反射器,与附在超声探头上的光纤束集成在一起。光捕捉器将收集到的反射光通过凹面镜腔多次随机反射重新分配到组织表面,如图14 所示。
图14 内置光捕捉器的高频光声成像探头示意图[39]Fig.14 Schematic diagram of high-frequency photoacoustic imaging probe with built-in light catcher[39]
采用该装置不仅可提高入射光在组织表面分布的均匀性,而且可增大深层组织内的有效光通量,改善深层组织的成像质量,如图15(彩图见期刊电子版)所示。与增强照明之前相比,光声信号强度可提高约30%,信噪比提高约20%。但是由于组织表面仅反射部分入射光,因此该方法对于成像质量的改善有限。此外,光捕捉器中的反射锥镜采用坚硬的透明材料制成,很难用于表面不平坦的成像物体,限制了其应用范围。未来可以进一步优化材料的选取并提高捕捉器对光的反射效率。
图15 采用光捕捉器增强照明(a)前(b)后的光声图像对比[39]Fig.15 Comparison of photoacoustic images before(a)and after(b)enhanced illumination by using light catcher[39]
4.1.4 优化的三维光声层析成像照明
Mc Larney 等[41]提出一种优化的三维PAT照明方式,将3D 打印技术和定制光纤照明相结合。首先,采用基于光线追踪的方法模拟光纤特性以及超声换能器阵列的几何形状和样品位置等,确定能够在球形表面上提供均匀照明的最佳光纤数量、排列方式和位置等。然后,根据仿真计算结果采用3D 打印制作样品室,样品室采用水密封,最大限度减少光能损失并确保激发出的压力波的声耦合(图16)。
图163 D 打印的样品室、支架、超声换能器阵列和待成像目标的(a)装配图和(b)实物照片[41]Fig.16 (a)Assembly diagram and(b)photograph of 3Dprinted sample chamber,entire holder,ultrasound matrix array transducer and imaged target[41]
该方法采用定制多模光纤束实现了均匀的样品照明,光纤的高数值孔径确保在单脉冲激励下获得高信噪比。与传统的单侧照明相比,该方法不需要校正光照方差,扩大了有效视场,提高了穿透深度,极大地改善了图像的整体质量。此外,采用3D 打印技术可以精确定位样品、光纤和超声换能器阵列,并调整照明位置、探测器阵列的大小和方向以及样品位置等。但是,图像重建精度仍然受到其他因素的影响,例如层析成像检测的覆盖范围有限,不完全层析也会导致中心部分图像产生畸变。另外,该方法适用于孤立的目标,对于不同几何形状的目标,还需要设计专用的探测器阵列。
外部照明模式在成像深度上有很大的局限性,穿透深度一般限制在距离组织表面几厘米的范围内。Li 等[42]设计了一种内部照明的PAT 系统,在样品内嵌入一根带有圆柱形扩散器的定制光纤,光纤的尖端被融合到一个针状的二氧化硅扩散器中,通过光纤照射整个样品,并将光传送到样品内的感兴趣区域,如图17 所示。该方法可以解决外部照明中的低组织穿透性问题,对于体积较大的目标可以实现更均匀的照明,提高成像质量,如图18 所示。
图17 (a)光纤扩散器实物照片;(b)光纤扩散器照明[42]Fig.17 (a)Photograph of the fiber diffuser;(b)fiber diffuse illumination[42]
图18 (a)外部照明与(b)内部照明PAT 图像对比[42]Fig.18 PAT images with(a)external and(b)internal illumination[42]
内部照明模式虽然克服了组织表面的光衰减对穿透深度的限制,但最终的穿透深度仍然受到光纤扩散器与组织之间光衰减的限制。除此之外,当不需要照亮光纤扩散器周围的所有组织时,需要对光纤扩散器做进一步优化,以实现特定方向的照明。未来可以进一步提高光纤扩散器的功率,使其能够适应更大体积的成像物体。
不均匀和不稳定照明对PAT 成像质量的影响主要由照明模式、组织几何形状和光学特性引起的光衰减以及组织内部光吸收和散射的非均匀性造成。本文针对上述3 种原因对现有的解决方法进行了归纳和总结。
目前,对于该问题的解决仍然面临着诸多问题,例如:为了实现光吸收系数和散射系数空间分布图的高分辨率重建,需要不断提高计算机的存储空间和计算速度,因此需要不断更新与完善数据压缩算法;在本文总结的解决方案中大多利用人造仿体或单个活体器官验证方案的可行性,其临床可移植性还有待验证;现有的方法大多解决光在组织内部传输的非均匀性所致的图像质量下降问题,对于相邻组织之间边界处的光传输仍需进一步研究;对于更深层的组织特别是组织中心位置处的成像质量仍需进一步提高。
除此之外,理想条件下的PAT 忽略了诸多影响图像重建质量的因素,例如测量数据集的完备性、组织声学特性的不均匀性、成像目标的运动状态以及超声探测器的性能等。在对非理想条件PAT 的研究中,需要综合考虑上述因素,全面提高成像质量和精度。
近年来,深度学习等人工智能(Artificial Intelligence,AI)技术的广泛应用为实现高质量PAT提供了诸多解决方案。例如:利用深度神经网络准确预测光吸收目标[43]、采用多层小波卷积神经网络实现高光通量光声成像[44]以及基于U-Net的qPAT[45]等,解决了低信噪比、背景噪声以及光辐射通量和横向分辨率随光照深度下降的问题,极大地提高了光声图像的质量。未来,进一步提高效率和准确性是该领域的主要研究方向,此外缺乏实际临床应用、缺乏高质量的训练数据集以及AI算法的完善也是未来研究中需要解决的问题[46-47]。