基于聚焦声场模型的光声层析成像时间延迟快速校正反投影方法*

罗晓飞 王波 彭宽 肖嘉莹

(中南大学基础医学院,生物医学工程系,长沙 410083)

光声层析成像是一种发展迅速的成像技术,其可提供生物组织的结构和功能信息,结合了光学成像高光学对比度与声学成像高穿透深度的优点.然而,由于现有的反投影成像算法通常将围绕目标扫描的超声换能器等效为一个点探测器,导致非中心成像区域图像的切向模糊,严重影响了图像质量.本文提出一种新的光声层析成像算法,其采用聚焦声场等效模型,可以快速有效地克服换能器孔径效应所造成的声场畸变,恢复非中心成像区域的切向分辨率.仿真结果表明,该方法对直径5 mm,距离旋转中心6 mm 的目标,切向分辨率提升至少达2 倍.实验结果表明,该方法可以有效地恢复边缘图像的切向模糊,使得复杂目标的微小结构能被清晰探测.这种新方法为传统的反投影方法提供了一种有价值的替代选择,对基于圆/球扫描的光声层析成像系统的设计具有重要的指导作用.

1 引言

光声层析成像(Photoacoustic tomography,PAT)结合了超声高穿透深度和光学成像高光学对比度的优点,能够在光学扩散极限的深度范围内对生物组织进行高分辨率的结构和功能成像[1-4].目前,它已被广泛应用于各种生物医学成像和检测[5],如脑血管可视化[6]、小动物全身成像[7]、乳腺癌检测[8]和手指关节成像[9]等.

在PAT 中,通常对整个样本使用宽场照明,假设每个换能器覆盖整个目标图像域,并通过围绕样品的单换能器扫描或使用换能器阵列来采集超声信号.然后利用图像重建方法一次性恢复整个区域的光吸收分布[10].PAT 扫描的几何形状可以是线性的、圆形的或球形的,也可以是上述几何形状的组合.与线性扫描[11]相比,圆形(二维成像)[9]和球形(三维成像)[10]扫描可以覆盖更完整的目标视角,以提供更好的成像质量.除了优化扫描几何结构外,快速准确的图像重建算法也是PAT 研究的热点[12].

在各种已有的PAT 重建算法中,反投影法(Back-projection,BP)因其稳定性好、计算量小而被广泛采用[13,14].然而,包括BP 方法在内的大多数当前的PAT 重建算法都简单地将超声换能器建模为点探测器,然而,实际上这些换能器的探测面通常是平面的,尺寸上可以达到几毫米大小[13,15].在重建算法中,超声换能器探测面的这种几何模型与实际情况的不匹配导致了“有限孔径效应”[15,16].由于这种影响,在基于圆/球扫描的PAT 中,当目标远离旋转中心时,切向分辨率会迅速下降,从而导致目标物成像结果的切向模糊.理论分析表明,其切向分辨率与距旋转中心的距离成正比,在换能器探测面处与换能器孔径大小相等[17].

目前,已有改善被拉伸的切向分辨率的方法.很多研究人员提出了在基于圆周扫描的PAT 中使用聚焦换能器并结合虚拟点探测器的方法[18].但这些方法并不能可靠地重建目标的高频部分[19].为了克服“有限孔径效应”,一些研究者采用了基于各种最优化模型的图像重建方法,结合换能器的几何形状和频率响应对图像进行重建,但其计算成本太高.反卷积作为一种图像后处理方法也有可能提高切向分辨率,但作为一种典型的逆方法,它会引起较强的图像噪声[20,21].因此,仍然需要发展新的PAT 重建算法来克服基于圆形/球面扫描的PAT中有限换能器尺寸的影响.

为此,本文提出一种新的PAT图像重建算法.采用聚焦声场等效模型,可以快速有效地克服换能器孔径效应所造成的声场畸变,恢复非中心成像区域的切向分辨率.本文将给出该算法的理论模型,并对其进行数值仿真以及实验验证.

2 算法模型

2.1 算法理论模型

图1 是基于环形扫描的PAT 扫描结构示意图.其中,脉冲激光照射在目标物上.换能器以目标物为圆心,沿着环形轨迹对目标进行扫描,最终通过一定的成像算法来重建扫描区域中的光吸收能的分布.

图1 基于二维环形扫描的PAT 系统Fig.1.2D circular-scanning-based PAT system.

其中,反投影法具有快速、简单、稳健性高、计算量小、易于并行的优点,是当前PAT 中最常用的算法.其核心是先测量像素与每个传感器之间的时间延迟 Δt,然后由相应时间延迟后的传感器信号S(i,Δt(i,r) 之和得到像素值:

其中,N为换能器扫描过的总切向角度数,r为要重建的像素坐标,S(i,t) 为换能器在第i个位置接收到的光声信号,t为时间.Δt由所求像素在换能器声场中的位置,即像素与换能器的相对位置所决定.其中,使用当前的反投影算法时,一般将换能器当作点探测器,其中心位置当作模型中点换能器的位置,如图2(a)所示.在这种情况下,投影线(或者更确切地说是等时延迟线)是一组以探测器位置为中心的同心曲线(如图2(a)所示的蓝线),并且对于位于成像域中 (x,y) 处的任意像素,其值为

图2 不同重建模型的反投影算法示意图 (a)点状换能器模型;(b)无限大小换能器模型;(c)聚焦声场换能器模型Fig.2.Schematics of the back-projection algorithm of different reconstruction models:(a) The model for point-like ultrasonic transducer;(b) the model for ultrasonic transducer with infinite element size;(c) the model for ultrasonic transducer with focused sound field.

其中R是换能器扫描轨迹的半径,M是换能器总数,θi是第i个换能器的角坐标v,是介质中的声速.但是这个模型忽略了换能器的探测面形状和大小的影响.这种模型误差,会导致随着目标远离扫描中心,其切向分辨率急速上升.理论表明,在换能器探测面处,目标的切向分辨率变得与换能器的尺寸相等.

图2(b)是另一种常见的PAT 重建模型.该模型中将换能器视为无限大的平面探测器,其反投影线是与探测器平面平行的一组直线.这种情况下的图像重建与Radon 逆变换非常相似,其点 (x,y) 处的像素值如下:

然而,该方法仅适用于换能器直径相对于重建区域较大,以及换能器声场准直性非常高的有限几种情况.而在常规的PAT 中,超声换能器的探测面通常为3—10 mm 的有限尺寸大小,其声场在距离换能器较近时表现为直线向前传播,因此其反投影线近似于有限尺寸大小换能器模型中(图2(b))中的直线,而在远场情况下其声场表现为一定角度的扩散模型,其反投影线又近似于图2(a)中那样的一组同心圆,如图2(c)所示.为了使得反投影中换能器近场和远场条件下都有较小的相位误差,本文这里采用一种基于聚焦声场的模型,来近似有限尺寸超声换能器的声场.其公式为

其中,v为声速,在聚焦探头中,w0和z0分别为声场的束腰半径和半景深,其由换能器的数值孔径与中心波长决定.在平探头中使用聚焦声场的方法来近似平探头的声场,把聚焦声场的束腰半径(通常被称为焦宽)w0定义为换能器孔径l的一半,式(4)中a和b分别表示所求像素与换能器探测面中心点之间的轴向以及横向距离.另外,根据平探头的发散角θ可以得到该聚焦声场的焦深,其中fc为换能器的中心频率,即有:

可以看出,当a无限小时:

因此,表达式接近近场直线传播模型.

当a很大时,Δt的表达式为

这时接近远场同心圆模型.

2.2 数值模拟与仿真

在以下模拟中,在x轴上有7 个点均匀分布在—6—6 mm 之间.平面换能器的带宽为60%,直径为5 mm.对点目标进行半径为20 mm 的二维环形扫描,坐标原点设在旋转中心,因此超声换能器相对旋转中心的坐标为(—20 mm,0 mm).全扫描共有360 个探测器,角度间隔为1°.在这些仿真中,测试了不同的换能器中心频率,分别为1 MHz,3 MHz,5 MHz和10 MHz.采样频率为100 MHz.本文中除特殊说明外,一般首先对原始信号进行希尔伯特变换来得到复光声信号,重建后采用像素值的模作为最终输出图像.另外,这里还提取了模拟的点目标的切向剖面的半高宽(full width half maximum,FWHM)作为切向分辨率.在噪声评估时,采用了1000 个试验来模拟信号的随机涨落.在每一次试验中,使用Box-Muller 方法对数据添加高斯白噪声,标准差为旋转中心点状目标信号幅度的5 %.将每个目标的1000 个重建幅值的平均值作为信号As,将这些重建幅值的标准差作为噪声An,计算信噪比(signal to noise ratio,SNR)为

图3 显示了不同中心频率下的重建结果,以探讨该方法的适应性.不同的行表示不同算法的重建结果,不同的列表示不同的中心频率.从结果看出,中心频率为1 MHz 的时候,点探头模型相较于无限大模型能将7 个点目标很好地重建,这是因为中心频率越小,波长越大,根据单缝隙干涉原理,此时超声换能器的发散角就越大,因而有限尺寸聚焦声场模型就越接近点探头模型.

但是从图3(a)—3(d)可以清楚地看到,随着频率的增高,换能器的发散角增大,有限换能器孔径的影响就越显著.虽然第一行中,在各频率下x=0和2 mm 处的目标都能较好重建,但是对于x=4或者6 mm 处的目标,随着频率的增大,如在3 MHz和5 MHz 下,其切向模糊就会增大.对于较高的10 MHz 频率则该现象更为显著.从图3(e)—3(h)可以看出,这种切向分辨率增大的现象仍然存在.但是,通过聚焦声场的方法进行时间延迟补偿,所有频率下的偏心目标的切向模糊都得到了很好的恢复.

图3 不同重建模型下换能器中心频率变化的仿真结果 (a)—(d)当中心频率分别为1 MHz,3 MHz,5 MHz和10 MHz 时,用传统点模型方法得到的结果;(e)—(h)无限大平面模型的相应结果;(i)—(l)聚焦声场模型的相应结果.所采用的平面换能器直径为6 mm,7 个点的位置分别为—6 mm,—4 mm,—2 mm,0 mm,2 mm,4 mm和6 mm.在模拟数据中加入5%的噪声,所有图像的振幅都归一化为1Fig.3.Simulation results with different reconstruction methods when the transducer central frequency varies:(a)—(d) The results with the conventional point-like method when the central frequency are 1 MHz,3 MHz,5 MHz and 10 MHz,respectively;(e)—(h)the corresponding results with the model of the infinite element size;(i)—(l) the corresponding results with the model of focused sound field.The diameter of the employed planar transducer was 6 mm,and the locations of the seven points are —6 mm,—4 mm,—2 mm,0 mm,2 mm,4 mm,and 6 mm,respectively.The simulated data were added with 5% noise,and the amplitudes of all the images were normalized to 1 for comparison.

为了更好地说明基于聚焦声场模型对该换能器有限孔径效应的修正,图3 为中心频率为3 MHz的4 个目标(0 mm,2 mm,4 mm和6 mm)的实数重建值切向截面被提取,并在图4(a)—4(d)中进行比较.另外,表1 列出了图3 中x=6 mm 处最右侧目标的SNR和FWHM.从表1 可以看出,与传统的点目标模型方法相比,基于聚焦声场的方法将x=6 mm 处目标的切向分辨率至少提高了2 倍,同时SNR 也至少提高了0.5 dB.

表1 Table 1 不同算法和换能器频率下偏离中心6 mm 目标的FWHM和SNRTable 1.The FWHM and SNR results of 6 mm off-center target with different methods and transducer frequencies.

图4 由3 MHz 中心频率超声换能器模拟得到的4 个目标的横向剖面,红线、蓝线和绿线分别表示点状换能器模型、无限大小换能器模型和聚焦声场换能器模型的结果 (a) 0 mm;(b) 2 mm;(c) 4 mm;(d) 6 mmFig.4.The simulated lateral profiles for the four targets obtained from 3 MHz central frequency ultrasonic transducer.The locations of the four points are 0 mm (a),2 mm (b),4 mm (c),and 6 mm (d),respectively.The red,blue and green lines represent the results from the conventional point-like model,the infinite element size model,and the focused sound field model,respectively.

3 点目标实验与动物实验

为更全面考察本文提出的新算法对环形扫描的光声层析成像的改善效果,通过点目标和鼠脑等复杂目标的实验验证了本章方法的有效性.本文采用典型的单个超声换能器的环形扫描系统[22],系统装置示意图如图5(a)所示.从光学参量振荡器激光器(SpitLight 600 OPO-532,Innolas)发出的脉冲740 nm 激光被反射镜引导,通过一块凹透镜给目标提供均匀的扩散照明.激光在样品表面形成的直径大约为3 cm,光强约为5 mJ/cm2,重复频率为20 Hz.二维扫描总步数为360 步,角度步长为1°.旋转中心到换能器检测表面的距离约为20 mm.信号首先由脉冲发生器/接收器(DPR500,Ultrasonics)放大,然后在计算机中用采集卡(NI-5124,12 bit,100 MHz)进行数字化.整个系统通过激光进行同步,并将数据采集到硬盘上进行后续处理.本文中所采用的平面换能器(Blatek Industries,5-3)中心频率为5 MHz,带宽为80 %,压电单元直径为3 mm.脉冲发生/接收器滤波器设置为1—10 MHz 的带通,此时用一根10 μm 钨丝测量得到的换能器光声脉冲和频谱响应分别如图5(b)和5(c)所示.

图5 PAT 系统和超声换能器的传输特性示意图 (a)二维PAT 系统图;(b)超声换能器光声脉冲响应信号;(c)超声换能器光声频谱图Fig.5.Schematic of the PAT system and transmit characteristic of the transducer:(a) Schematic of the 2D PAT system;(b) the photoacoustic impulse waveform of the transducer;(c)spectral diagram of the transducer.

点目标仿体实验中所用琼脂仿体的直径为3 cm,散射系数为1 mm—1,吸收系数为0.007 mm—1.将8 根铅笔头(0.5 mm 厚) 垂直插入琼脂仿体作为点目标,如图6(a)所示.用传统点探头模型的方法、无限大平面换能器的方法和有限尺寸的聚焦声场模型方法重建的复数幅值图像分别在图6(c)—6(e)中示出.与照片图像6(a)对比可以看出,这3 幅图都能可靠地揭示出8 个目标的位置.然而,由于有限换能器孔径的影响,采用点换能器模型时远离旋转中心的7 个目标在切向上严重变形(如图6(c)).采用无限大平面换能器模型时,其切向分辨率也提高有限(如图6(d)).而通过聚焦声场模型对时间延迟进行补偿,这些目标的切向分辨率得到了很好的恢复,如图6(e)所示.

图6 不同重建算法对点目标成像的实验结果 (a)仿体照片;(b)8 个点目标的切向分辨率分布图;(c)点状模型BP 算法;(d)无限大小模型BP 算法;(e)聚焦声场模型BP 算法Fig.6.The experimental results of point target imaging with different reconstruction algorithms:(a) The photograph of phantom;(b) tangential resolution distributions of eight point targets;(c) the BP method of point-like model;(d) the BP method of infinite element size model;(e) the BP method of focused sound field model.

本文对各重建模型下8 个目标的切向分辨率进行了测量,结果如图6(b)所示.在3 条点线图中,红色折线图代表点模型重建算法,蓝色折线图代表无限大平面换能器模型的重建结果,绿色折线图代表基于聚焦声场模型的重建结果.从可以看出,对于距离旋转中心为8 mm 的6 号目标,基于聚焦声场模型的方法的FWHM 比传统点模型和无限大平面模型的方法分别提高了约5.5 倍和2 倍.这些结果与我们在模拟研究中得到的结论是一致的.

鼠脑成像实验是为了验证本文所提的算法对复杂活体目标成像的能力.选用KM 小鼠(35 g,3 周龄,湖南SJA 实验动物有限公司),该动物实验方案已获得中南大学动物伦理委员会的批准.在成像前,用脱毛膏去除小鼠头部的毛发.用戊巴比妥钠(20 mg/kg,腹腔注射)麻醉小鼠,使其在整个实验过程中保持不动,用一个自制的动物支撑架来固定鼠头,在水和鼠脑之间用一层透明的薄膜来进行超声耦合.图7(a)是当时实验结束后拍摄的掀开小鼠头皮后的脑部血管照片图,图7(b)—7(e)是不同重建算法对小鼠脑成像的实数值重建结果,对比可看出,图7(b)中箭头1 所标记的上矢状窦血管远端被明显拉伸,这是因为点目标模型不符合实验中所用的探测面达9 mm 的超声换能器.相比之下,无限尺寸换能器模型图7(b)和有限尺寸换能器模型图7(c)均能给出较好的重建结果.这两种算法重建的小鼠的脑部血管网络差别也都较小,都和血管照片图7(a)有较好的对应,如箭头2到7 所标记的其他较小的血管所示.

图7 不同重建算法对小鼠模型成像的实验结果 (a)小鼠脑图片;(b)点状模型BP 算法;(c)无限大小模型BP 算法;(d) 聚焦声场模型BP 算法Fig.7.The experimental results of mouse brain imaging with different reconstruction algorithms:(a) The photograph of the mouse brain;(b) the BP method of point-like model;(c) the BP method of infinite element size model;(d) the BP method of focused sound field model.

4 结论

本文提出了一种基于聚焦声场模型的反投影重建法来对重建图像进行时间延迟校正,以克服二维圆扫描PAT 中的“有限孔径效应”.与传统的点换能器模型和无限大平面模型相比,本文方法通过施加精确的时间延迟,可以很好地恢复偏心区域退化的切向分辨率,从而恢复这些区域垂直轴向的图像特征.基于聚焦声场的方法继承了传统反投影方法的简单性、稳定性和低计算量的优点,是一种有价值的替代方法.该方法对基于圆/球扫描的PAT系统设计具有重要的指导意义.