1 引言

光声成像作为一种新型的无损医学成像技术,近年来在国内外受到了广泛关注^[1-3]。在脉冲激光的照射下,生物组织吸收光能之后由于热弹性效应会产生压力,进而释放出声波,这种现象即为光声效应。光声成像技术通过超声换能器探测光声效应产生的声波信号,反演光吸收的强弱分布,从而可以获得生物组织的二维或三维图像。这种成像方式使得光声成像技术具备光学成像对比度高和超声成像穿透深^[4]的优点。光声成像技术可以对血管^[5]、脑功能^[6],甚至细胞核^[7]进行成像,在肿瘤检测^[8]等生物医学领域具有广阔的应用前景。

光声断层成像是使用非聚焦激光照射的光声成像方式,这种成像方式的成像深度能够达到50 mm左右。光声断层成像需要在多个角度探测光声信号,通过一定的反演算法来重构图像。目前,主要的光声重建算法有滤波反投影算法^[9-10]、时间反演法^[11-12]和傅里叶变换重建算法^[13]等。实际的光声成像往往存在采样数据不足的问题,从而严重影响了上述算法的图像重建效果。基于压缩感知理论的光声图像重建算法可以提高稀疏角度采样时的光声图像重建质量^[14-15],然而,这类重建算法需要建立的投影矩阵会随着重建图像大小的增加而变得非常大,远远超过了一般计算机内存空间所能存储的范围^[14]。另一种提高稀疏角度采样时光声图像重建质量的方法是对光声信号进行超分辨率重建。将低采样率的光声数据进行超分辨率重建后,再以滤波反投影等算法进行图像重建,最终可以获得较高的成像质量。

近几年迅速发展起来的基于字典学习的稀疏表示方法在图像超分辨重建^[16]、图像分类^[17]等领域应用广泛。在核磁共振成像^[18]和计算机断层扫描(CT)成像^[19-20]中,以基于字典学习的稀疏表示项作为正则化项,对重建图像的解空间进行约束,可以提高低采样率和低信噪比情况下的图像重建质量。利用字典学习的稀疏表示方法在对欠采样的CT正弦图进行超分辨率重建之后,可以改善滤波反投影法的重建质量,大幅减少伪迹^[21]。在光声成像领域,基于字典学习的稀疏表示方法可以用来去除光声显微镜中的回波信号^[22]以及对重建图像去噪^[23]。本文将基于字典学习的稀疏表示方法应用于光声信号的重建,恢复稀疏采样时缺失的信号数据,并通过仿真实验来验证经过重建的光声信号可以显著改善图像的重建质量这一结论。

2 光声信号的恢复模型

光声成像的基本方程^[9]为

∇2p(r,t)-1c2∂2p(r,t)∂t2=-βCpA(r)∂I(t)∂t,(1)

式中:p(r,t)为位置r处的声压;t为时间;A(r)为电磁波吸收分布;I(t)为激光脉冲能量函数;c为生物组织中的声速;β和C_p分别为生物组织的等压膨胀系数和比热容。假设声速恒定,求解上述方程后可以得到在r₀位置处t时刻超声换能器的声压^[9]为

p(r0,t)=β4πCp∂∂t∯r'-r0=ctA(r')td2r',(2)

式中:r'为声源的位置向量。

在实际的光声成像中,通过超声换能器阵列或者使用单个换能器进行扫描,可以获得多个位置处的光声信号。但是受硬件条件和扫描时间的制约,实际成像过程中的采样点数往往远小于理想情况。此外,超声换能器有一定的接收信号面积和角度,这使得光声信号之间会产生混叠和模糊。欠采样数据Y与理想采样数据X间的关系为

Y=SX+ε,(3)

式中:Y∈R^M×T,M为实际中不足的采样点数(M<N),T为信号时间长度;X∈R^N×T,N为理想的采样点数;S∈R^M×N为下采样和模糊矩阵;ε为噪声。从欠采样信号中恢复出完备信号是一个病态问题,一般需要利用信号本身的先验知识对解空间进行约束才能精确求解。本研究通过引入基于字典学习的稀疏表示项,对(3)式中的理想采样信号X进行求解。

如果一个信号的非零元素个数远小于信号的维度,就可以认为该信号是稀疏的。对于自然信号而言,虽然其本身并不稀疏,但是可通过变换域或者过完备字典来对其进行稀疏表示。如果信号y∈R^m可以被一个m×n(m≪n)型矩阵D中的k个列向量线性表达,并且k≪n,那么称矩阵D就是能将信号y进行稀疏表示的过完备字典,D中每个列向量称为一个原子。稀疏表示的目标函数^[24]为

α∈Rnminα‖0 s.t. y=Dα,(4)

式中:α为稀疏系数;‖α‖₀为α的l₀范数,即α中非零元素的个数,用来表示α的稀疏度。信号y在过完备字典D稀疏表示下,稀疏系数α可以用正交匹配追踪算法(OMP)来求解。

为了对光声信号进行最稀疏表示,本研究通过已知光声信号集对初始随机字典基于K-SVD算法进行训练,即进行字典学习。利用包含N组信号的信号集{x₁,x₂,…,x_N}进行字典学习,其目标函数为

minD,αiNi=1∑i=1N‖xi-Dαi‖22 s.t. ‖α‖0≤k0,1≤i≤N,(5)

式中:α_i为第i组信号的稀疏系数向量;‖· ‖22为l₂范数。

在光声信号重建中,以引入基于学习得到的字典的稀疏表示项作为约束,可以有效缩小可行解空间。在光声信号重建中,为了减小计算量,将M×T大小的光声信号Y交叠分成n×n大小,并重排为n²维的列向量y_i。对于每个y_i的重建,其目标函数为

minx,α∑j=1n2([sxi]j-[yi]j)2+λ‖x-Dα‖22+‖α‖0,(6)

式中:x_i∈R^K×1表示y_i对应的待重建的信号;s为下采样矩阵;λ为正则化参数,用以平衡数据保真项与稀疏表示项的影响;x为待重建的信号矩阵。由于(6)式中有2个变量x_i和α,因此需要采取交替求解的策略,即固定其中一个变量对另一个变量进行求解。当稀疏向量α固定时,目标函数变为

minx∑j=1n2([sxi]j-[yi]j)2+λ‖xi-Dα‖22。(7)

利用抛物线分离的优化方法^[25]对方程(7)进行求解,迭代过程为

xikt=xikt-1-∑j=1n2sjksxi]t-1]j-[yi]j+2λ[[xi]t-1-Dα]k∑j=1n2sjk∑l=1Ksjl+2λ,k=1,…,K,(8)

式中:t为迭代次数。

字典学习与光声信号重建的算法流程如图1所示,其中ε为残差阈值。

图 1. (a)字典学习算法流程;(b)基于字典学习的光声信号重建算法流程

Fig. 1. Flow charts of (a) dictionary learning algorithm and (b) photoacoustic signal reconstruction algorithm based on dictionary learning

下载图片 查看所有图片

3 仿真实验

在验证所提算法有效性的模拟实验中,用K-wave工具箱^[26]来模拟光声成像过程中光声信号的产生。计算区域大小为10 mm×10 mm,即划分为256×256的网格,模拟区域边界使用完美匹配层,以使得光声信号向前传播。模拟中设置声速为1500 m/s,模拟光声信号时间间隔为10 ns,每个位置共采集1207个点。为了验证算法的稳健性,在模拟中加入1%的噪声,即使得信噪比为40 dB。

首先,为了通过学习的方式得到对光声信号稀疏表示的字典,模拟生成理想采样情况下的光声信号。在超声换能器所在的圆周上均匀地选取160个角度采集光声信号,如图2所示。得到的理想采样光声信号数据空间大小为160×1207(160个采样点,1207个时间点),可视化结果如图3所示,其中

图 2. 模拟区域结构

Fig. 2. Schematic of simulation

下载图片 查看所有图片

C方向为换能器序列,T方向为时间序列。

为了计算方便,将理想采样点数时模拟得到的160×1207大小的光声信号重叠分块,重排成若干64×1大小的列向量(每个信号小块均保留了空间和时间两个维度的信息)作为样本集。去掉样本集中方差较小的样本,将剩下的约20000个样本归一化后,通过图1(a)所示的算法流程对初始随机字典进行训练。表1所示为在算法实现中选取的关键参数。图4所示为初始随机字典和从光声信号样本中学习后得到的字典,字典中每个小方格即为字典的一个原子。由图4可知,与初始字典中每个原子的结构都杂乱无章不同,学习得到的字典的绝大部分原子都展现出了明显的结构特征,这说明字典已经学习到了光声信号的先验信息,在对光声信号稀疏表示时可以获得更好的稀疏度。

图 3. (a)模拟出的光声数据;(b)图3(a)中白线处的信号

Fig. 3. (a) Simulated photoacoustic data; (b) signal at where marked with white line in Fig. 3 (a)

下载图片 查看所有图片

图 4. (a)初始随机字典;(b)从样本集中学习得到的字典

Fig. 4. (a) Initial random dictionary; (b) learned dictionary from sample set

下载图片 查看所有图片

学习得到字典后,通过图1(b)所示的基于字典学习的光声信号重建算法流程处理欠采样光声信号(在160个采样点的光声信号均匀时采样得到的40个采样点的信号)。使用图5(a)所示的图像作为成像目标,模拟产生160个采样点的光声信号,如图5(b)所示。分别将160个采样点的光声信号、下采样之后的光声信号、经所提算法重构恢复的光声信号用时间反演法重建光声图像,重建结果如图5(c)和 图5(d)所示。从图5(c)和图5(d)所示的图像重建结果中可以看到:下采样之后的光声图像出现了大量伪影,细节处(如图5(c)虚框内部分)难以分辨;经所提算法处理的重建图像不仅伪影显著减少,而且在一些细节处(如图5(d)虚框内部分)也实现了足够的分辨。

图 5. (a)原始图像;(b) 160个采样点下的光声重建图像;(c) 40个采样点下的光声重建图像;(d) 40个采样点下基于所提算法的光声重建图像

Fig. 5. (a) Original image; (b) photoacoustic image reconstructed by measurement with 160 detector locations; (c) photoacoustic image reconstructed by measurement with 40 detector locations; (d) photoacoustic image reconstructed by measurement with 40 detector locations based on proposed algorithm in this paper

下载图片 查看所有图片

为了定量表征图像的重建质量,引入峰值信噪比(PSNR)和结构相似度(SSIM)^[27]。均方误差(MSE)、峰值信噪比、结构相似度的表达式分别为

eMSE=1MN∑m=1M∑n=1N(umn-vmn)2,(9)rPSNR=10lg(IMAX2/eMSE),(10)αSSIM(u,v)=2μuμv2σuv+c2)(μu2+μv2+c1)(σu2+σv2+c2),(11)

式中:e_MSE为均方误差;I_MAX为图像最大强度;u和v分别为重建图像和原图像,大小均为M×N;r_PSNR越大,表明图像重建效果越好;μ_u和μ_v分别为图像u和v的平均强度;σ_u和σ_v分别为图像u和v的标准差;σ_uv为u和v的协方差;c₁和c₂为常数;α_SSIM(u,v)越大,表明图像重建效果越好。

在采样点数分别为40和80的情况下进行2组模拟实验。表2所示为在不同采样点数和信噪比下,所提算法、常用的插值方法和不进行信号恢复时重建图像的峰值信噪比、结构相似度。由表2可知,与原始情况和常用的插值方法相比,所提算法的峰值信噪比分别提高8 dB和3 dB。可见,所提算法具有良好的稳健性。

4 结论

针对光声成像中采样角度稀疏、采集数据不足的问题,本课题组提出了一种光声信号的恢复重建算法。通过从光声信号样本集中学习得到的字典,对欠采样的光声信号稀疏表示进行迭代求解,从而恢复采样过程中缺失的信号。仿真实验结果表明,所提算法可以显著减少重建图像的伪迹,分辨出图像的细节,图像的峰值信噪比约提高了8 dB(与不进行信号恢复时重建图像的峰值信噪比相比)。为了检验所提算法的抗噪能力,在40 dB、30 dB、20 dB等信噪比下进行仿真,结果表明:所提算法具有良好的稳健性。然而,实际的光声采样过程更加复杂,影响光声信号质量的因素更多,结合实验来修正所提算法中光声信号降质模型的工作还有待进一步开展。