基于声学散射透镜的低成本光声断层成像方法.pdf

NATURASCIENCEMar.,20232023年3 月JOURNALCOHNANINNIVERSITYVol.59,No.2第5 9 卷第2 期南京大学学报（自然科学）DOI:10.13232/ki.jnju.2023.02.005基于声学散射透镜的低成本光声断层成像方法何君君1，殷杰1.2*,刘广东1,马倩倩1（1.南京科技职业学院，南京，2 10 0 48；2.南京大学声学研究所，南京，2 10 0 93）摘要：光声断层成像需要使用换能器阵列及多通道系统，硬件成本较高，制约了其广泛应用。为降低成本，设计了一种基于声学散射透镜的光声断层成像方法，利用声学多重散射过程中的多路径效应达到增大探测阵列等效数值孔径的作用，从而减少信号探测过程所需换能器数量.通过数值模拟实验，首先分析了散射透镜格林函数的精确测量方法，随后利用四个散射透镜共同作用实现了高质量光声断层成像，并对影响成像质量的相关因素进行了具体研究，提出的散射透镜无须复杂制作工艺或特殊材料，系统只需四个超声换能器即可实现快速成像，可降低光声断层成像系统的硬件成本。关键词：光声断层成像，低成本，散射透镜，格林函数中图分类号：0 42 9文献标志码：ALow-cost photoacoustic tomography based on acoustic scattering lensesHe Junjun,Yin Jiel.2*,Liu Guangdong,Ma Qianqian(1.Nanjing Polytechnic Institute,Nanjing,210048,China;2.Key Laboratory of Modern Acoustics,MOE,Institute of Acoustics,Nanjing University,Nanjing,210093,China)Abstract:Photoacoustic tomography(PAT)is a hybrid biomedical imaging modality based on photoacoustic(PA)effect,itcombines optical excitation with acoustic detection.As a no-ionizing and non-invasive technology,it can provide bothstructural and chemical informations in pre-clinical and clinical applications,including tumor diagnosis,drug deliveringmonitoring,little animal whole-body imaging,brain funtional imaging and so on.As a result,PAT has become one of thefastest growing imaging method in the past decade.One chanllenge that has continuously attracted the attention of researchersand engineerers in biomedical imaging field is to reduce the cost without degradation of performance.PAT system usuallyinvolves transducer array and multi-channels system,which lead to a high cost and bring limitations on widely application.Aswe know,with the known Greens function,multi-paths effect of multi-scattering can enlarge equivalent numerical aperture ofthe detecting array,and PA images can be reconstructed with much less transducers by time reversal algorithm.In this paper,we propose a PAT modality based on acoustic scattering lenses,which utilizing the multi-paths effect of multi-scattering.Firstly,a method for accurate of Greens function measuring is verified,quantitative studies indicate that,the accurracy ofGreens function measuring can be improved by choosing the appropriate coherent gathering zones and increasing the numberof transducers participated in Greens function measuring.Then,we demostrate that high-quality PA images can be achievedwith four scattering lenses,while PA image quality can be improved by controlling the signal length.Scattering lens proposedin this paper needs not precision manufacture or special material,and just four transducers can accomplish quick imaging,withwhich PAT system cost can be reduced.Key words:photoacoustic tomography,low-cost,scattering lens,Greens function基金项目：国家自然科学基金（118 7 42 17），江苏省自然科学基金（BK20181077），江苏高校“青蓝工程”收稿日期：2 0 2 2-12-0 6*通讯联系人,E-mail:232南京大学学报（自然科学）第5 9 卷光声断层成像（PhotoacousticTomography，PAT)是一种基于光声效应的新型非侵人式生物医学成像方法，其具有光学成像对比度高以及声学成像可在深层组织获得高分辨率的特点 1-3 .PAT成像过程中，首先利用非聚焦脉冲激光束照射待测生物组织，生物组织中的光吸收体吸收了光波能量后由于热胀冷缩效应激发出超声波并向周边传播；随后，布置在组织周边的超声换能器接收超声波并将其转化为电信号传输至上位机，上位机通过特定的重构算法形成生物组织内部光吸收系数分布图像.PAT本质上是光学成像，但光学参量的载体是超声波，它打破了光学散射对成像深度的约束，可以在5 6 cm的成像深度获得良好的声学空间分辨率。近年来,PAT在诸多医学领域展现出了广阔的应用前景，如：肿瘤诊断 4-5 、血管疾病诊断 6 、小动物全身成像 7-8 、关节炎诊断 9、药物输运监测 10 以及其他各类医学应用 11-13 早期的PAT系统大多采用单换能器单通道的硬件配置 14，利用步进电机带动换能器沿预先设定的轨道移动，根据合成孔径原理达到阵列成像的效果，这种方式存在鲁棒性较差，扫描时间过长等缺点，难以适应当今临床医学领域对快速成像的需求.近年来，国内外研究机构针对多通道PAT系统开展了相关研究并取得了系列进展 7-8 .然而，多通道系统成本高昂的缺点始终是制约PAT大范围应用的重要因素.为减少换能器及通道数量，国内外学者提出了基于声学透镜的PAT15-16、基于声学腔穴的PAT17等方法，但这些方法亦存在器件需要特殊材质制作、加工难度大、鲁棒性较差等问题.为解决上述问题，本文提出了一种基于声学散射透镜的低成本光声断层成像方法，散射透镜由若干尺寸相同的圆形钢柱构成，单个信号探测换能器固定在散射透镜一侧，成像区位于散射透镜的另一侧，在光声成像的过程中，由于随机散射过程的多路径效应，散射透镜起到了增大换能器等效数值孔径的作用，成像区至换能器的格林函数由超声干涉测量法测得，基于四个散射透镜全包围光吸收体共同作用，利用时间反转算法，可实现高质量的光声断层成像.所提出的散射透镜无须复杂制作工艺或特殊材质，只需四个超声换能器即可实现信号快速探测及成像，进而降低系统硬件成本1理论分析首先考虑如图1所示场景，为避免繁重的数值计算过程，本文采用了2 D模型.多路径效应的产生需要声波的多重散射，因而散射透镜要满足下列条件：（1)散射透镜的厚度要大于声波平均自由程；（2)散射体的尺寸要与声波波长接近.散射透镜核心元件为等距分布的钢制圆柱体，单个超声换能器位于散射透镜左侧ra处，该换能器与散射透镜的相对位置关系保持固定不变.在进行光声成像前，需要在散射透镜右侧放置一个超声换能器阵列，用于测量散射透镜的格林函数.光吸收体位于散射透镜右侧的成像区内，在脉冲激光束的照射下，光吸收体吸收光波能量，并向外辐射光声信号，光声信号的初始声压值与光吸收系数成比例关系，即：q(r,t)A(r)aH(t)/at(1)A（r）为成像区内任一点r处的光吸收系数，H(t)为光脉冲的时域波形，q(r,t)为源，c为声速，光声信号的波动方程为：1 0p(r,t)Vp(r,t)q(r,t)(2)at2部分光声信号经过散射透镜并被位于r处的换能器接收，接收到的光声信号可表示为：pa(rd,t)dtdrq(r,t)G(ra,r,t)(3)R其中，T为信号长度，R为被探测器包围的成像区域,G(rd,r,t)为位置ra与r之间的格林函数.将探测器接收到的信号在时域上反转并发射，反转声场将自适应地汇聚在源点处 18 ,即：p.(r,t)=dtdrpa(ra,t)G(ra,r,t)(4)0R式（4)在均匀介质中可演变为延迟叠加算法(Delay-and-Sum,DAS）,通常情况下,利用式（4)重构光声信号初始声压需要换能器阵列分布的曲面包围R，然而，散射透镜的存在增大了探测器的等效数值孔径，部分原本无法被接收到的信号成分由于多路径效应，经多次散射后被换能器接收，使得利用较少换能器进行光声断层成像成为可233何君君等：基于声学散射透镜的低成本光声断层成像方法第2 期能 19.由式（4)不难看出，成像实现的重点在于格林函数G(rd,r,t)的确定.格林函数的确定可以使用超声干涉测量法，其利用超声信号相关性分析，消除两组信号中共同路径的部分 2 0 .如图1所示，在散射透镜右侧放置一组超声换能器阵列，阵列中的各单元所在位置为rs1,rs2，r s n，散射透镜左侧ra处的换能器向散射透镜发射一个波形为s(t)宽带超声信号，散射透镜右侧换能器阵列各单元接收到的超声信号依次记为u(t),uz(t)，,u（t)，此时，成像区内某点r至超声换能器r.的格林函数可利用下式近似求得：G(ra,r,w)J_F(a)U(w)G(r,rs,)ds(5)其中 G(ra,r,),G(r,rs,w),U()分别是G(rd,r,t),G(r,rst,t),u,(t)的傅里叶变换形式,i=1,2,n;F(w)=S*（w)为波形滤波器,用于在格林函数重构过程中消除s(t)波形对计算的影响.在实际测量过程中，式（5)可离散化为：G(ra,r,a)ZF(a)U()G(r,rs,)(6)1由于模型中成像区和换能器阵列之间的介质为无损均匀介质，故G(r,rs,t)=(7)4元Si在地震干涉测量学等场合，格林函数的重构需要r1,rs，,rs所在的曲面包围ra及r并对其进行各向同性照射 2 0 .由于应用场景不同，本研究采用单侧照射方式，由相关理论可知，在公式（6)中对于格林函数重构起主要作用的叠加信号来自菲涅尔区内的换能器单元 2 0 1,菲涅尔区如图1绿色虚线所示.然而，由于散射透镜的散射效应，必然会造成信号有效相干叠加区发生变化，需要通过定量研究提高格林函数的计算精度.此外，参与叠加的信号数量也与格林函数计算精度相关，故完成一次测量后可将换能器阵列移位，重复一次测量过程,得到一组新的信号ui(t),u2(t),，u(t),并与第一次测量结果共同代入公式(6)中计算，以提高计算精度。u(t)Mm日u(t)有效叠加区Yd测试点换能器成像区移位u(t)Hui(t)Time(us)散射透镜换能器阵列图1散射透镜结构及格林函数测量过程示意图Fig.1Structure of the scattering lens and the illustra-tion for Greens function measuring2数值模拟实验2.1散射透镜的格林函数计算本文利用COMSOL软件建立了相应的数值模型，如图1所示,散射透镜为矩形,中心位置为（一7 mm，0 m m），尺寸设置为3 mm12mm.散射体为直径0.6mm的圆柱,边界条件为硬边界,相邻散射体中心间距为1mm，散射体分布密度为0.9 2 mm-2.理论上散射体排布方式可采用随机分布或规则分布，只要产生足够的多重散射即可，但在随机分布情况下，有可能出现某一区域散射体分布过于密集或过于稀疏的情况，导致多路径效应只集中在散射透镜的某一部分发生，为避免这一情况，根据经验分析，本文中散射体采用上述布置方式。一个超声换能器位于r（一11mm,0mm)处，一个含有2 3 个等间距分布单元的超声换能器线型阵列位于散射透镜右侧,其轴坐标为4mm，y 轴坐标分布范围为一5.5 5.5 mm.ra处的超声换能器向散射透镜发射中心频率为3 MHz的高斯脉冲信号：s(t)=2元 f8(1-1/f.)exp(一 元 fe(t 1/f)(8)s(t)波形如图1插图所示.换能器阵列中各单元接收到的信号即为u(t)，u（t)，u(t)，波形如图2 a所示，换能器阵列接收到的波阵面如图2b所示，从图中可以看出光声信号在时域上存在多次振荡，这是由于散射透镜对光声信号造成了多重散射，多重散射的出现满足了多路径效应的形成条件.234南京大学学报（自然科学）第5 9 卷1.0(a)u(t)0.50.0-0.5-1.051015202530Time(us)(b)Max2015105Min51015202530Time(s)图2(a)经过散射后的超声脉冲波形;(b)换能器阵列探测到的波阵面Fig.2(a)The US pulse scattered by scattering lens,(b)wavefront picked upby the transducerarray为了研究格林函数测量的精确性，我们在成像区中取一点测试点r（1m m，0 m m）,如图1中红点所示.令其发射中心频率为3 MHz波形为s(t)的超声信号，该信号被ra处的超声换能器接收，表示为u(t),可知u(t)=G(rd,rt,t)?s(t),其中G（r d,r t，t)表示r与ra之间实际的格林函数.随后，做r与ra的连线延长线，并与换能器阵列相交，以交点为中心，依次增加交点两侧参与格林函数计算的换能器单元数量，在不同换能器数目的情况下利用公式(6)计算格林函数G(ra,r，t),并令us(t)=G(rd,It,t)s(t),可通过比较us(t)与u（t)波形相似程度确定格林函数计算的精确性，本文选用相关系数比较二者相似程度，相关系数定义为：Couu,(r),u.(r)Ru,(r),u(r)=(9)Varu.(r)Varu(r)使用交点处及左右邻近共5 个换能器单元（=36.87大致为菲涅尔区）数据计算得到的格林函数所对应的u（t)及ut（t)波形如图3 a所示，此时相关系数R=0.839；使用11个换能器数据（=7 9.6 1)计算得到的格林函数所对应的u（t）及u（t)波形如图3 b所示，此时相关系数R=0.935，使用的换能器数目与R值之间的定量关系如图3 c所示.可见，由于散射透镜的出现，对格林函数重构起作用的换能器范围超出了传统的菲涅尔区，信号有效相干叠加区域与重构点的夹角=79.61，当换能器阵列超过该角度时，信号成分不能进行有效的相干叠加，格林函数计算精度随增大而减小：为进一步提高计算精度，将换能器阵列向右侧移动4mm后，重复上述测量过程，阵列中各单元测得一组新的波形数据ui(t),u(t),u(t),与第一次测量结果共同代人公式(6)中计算，在保持=79.61不变的情况下，可将R值提高至0.948，如图3 c中红点所示.可见增加有效叠加区内参与叠加计算的信号数量，亦可在一定程度上提升计算精度.按照上述方法，可根据成像实际需求设定成像区，并计算成像区内各点与ra处换能器之间的格林函数.2.2利用声透镜实现光声断层成像象在确定格林函数的计算方法后，可基于多个散射透镜同时作用实现高质量光声断层成像设置数值模型如图4所示，4个结构相同的散射透镜呈矩形布置，对应的4个换能器坐标分别为（一11mm，0 m m），(0 mm,11 mm),(11 mm,0 mm),(0 mm,一11mm）；四个圆形光吸收体位于矩阵中心位置，用黄色圆球表示，光吸收体的直径为0.6 mm，圆心分别位于（一1mm,1mm），（1m m,1m m),(1mm,一1mm),(-1 mm,一1mm)处,光吸收体光吸收系数设为1.0，周边其余介质设为0，光吸收体及周边介质声速c为15 0 0 ms1，密度为1000kgm-3.设定的成像区用灰色正方形表示，其4个顶点坐标分别为（一2.5 mm，2.5 m m），(2.5 mm,2.5 mm），,(2.5 m m,-2.5 m m),（一2.5 mm，一2.5 mm），随后，利用上文所述计235第2 期何君君等：基于声学散射透镜的低成本光声断层成像方法1.0(a)-U.()0.5.u.(t)0.0-0.5-1.00510152025Time(uS)1.0(b)一u.(t)0.5.u.(t)0.0-0.5-1.0-1.50510152025Time(uS)0.95(c)R0.900.853045607590105120(degree)图3 (a)使用5 个阵列单元测量数据时的u,(t)与u.(t)波形对比(=36.87);(b)使用11个阵列单元测量时的u,(t)与u(t)波形对比(=79.61);(c)格林函数计算精度与角的关系Fig.3(a)Calculated Greens function with 5 elements andtrue value(=36.87),(b)calculated Greens function with11 elements and true value(=79.61),(c)relationshinbetween Greens function computational accuracy and 算方法确定成像区内各点对应的格林函数.为方便展示散射透镜所起到的作用，我们首先将图4所示数值模型中4个散射透镜移除，只留下4个换能器，令激光照射光吸收体并激发出中心频率为2 MHz的光声信号，无散射透镜情况下换能器接收到的光声信号如图5 a所示；随后得到4个换能器在无散射透镜作用的情况下利用DAS算法重构光声信号初始分布场，并采用半峰全宽口Path22PathPath2图4基于散射透镜的光声断层成像数值模型Fig.4 Numerical model for photoacoustic tomographybased on four scattering lenses1.0(a)(b)Max20.51(ww)人00.0-1-0.5-2Min0510152025-2-112Time(us)X(mm)图5(a)无散射透镜时换能器接收到的光声信号；(b)无散射透镜情形下4个换能器利用DAS算法重构的光声断层图像Fig.5(a)Photoacoustic signal received by transducerwithout scattering lens,(b)reconstructed PAT imagesutilizing DAS algorithm by 4 transducers without scat-teringlens法(Full WidthHalfMaximum,FWHM)绘制光声图像，如图5 b所示，4个光吸收体的真实形状和位置用白色圆圈标注，从图中可以看出，由于换能器数目过少，重构的光声图像产生严重的畸变，无法从中得到有价值的信息。随后，将散射透镜复位，重复光声信号激发过程，4个散射透镜换能器接收到的光声信号分别为pal,Pa2，Pa,Pa 4,其波形如图6 a所示，从图中可见信号尾波幅值明显大于图2 a中的信号，这是由于周边散射透镜反射的光声信号在时域的后半段被换能器接收到，为模拟实际测量情况，将光声信号添加权重系数为0.1的高斯白噪声并代人公式（4）236第5 9 卷南京大学学报（自然科学）中，可分别得到时间反转声场pa,Pa,Pra,Pa，最终重构的光声信号初始分布场为：p.=p.十pr2+pra+P-4.采用FWHM法绘制p.，当代人公式（4)进行计算的pa信号长度为2 5 s时，重构图像如图6 b所示，图中白色圆圈代表光吸收体实际形状和位置，图中可见较为明显的伪像和畸变；当信号长度为14s时，重构图像如图6 c所示，4个光吸收体清晰可见，位置、形状均符合实际情况为定量研究信号长度与图像质量之间的关系，引人了峰值信噪比PSNR来衡量图像质量：MAXIPSNR=201g(10)VMSE1MSE:221(1.i)-K(1.)(11)mn式中,MAX,=1.0为光声图像像素灰度相对最大值，m,n分别为光声二维图像在和y方向上的维数，本文m=n=201,即图像像素为2 0 12 0 1;重构图像的像素灰度值标记为I(i，j)；同时，在数值模型中移除散射透镜，使用含12 0 个等间距分布单元、直径为10 mm的圆形换能器阵列对同一成像区域利用DAS算法得到标准光声图像，并将其作为参考图像，参考图像的像素灰度值标记为K(i,j).PSNR值越大代表两幅图像相似度越高，PSNR与信号长度的关系如图6 d所示.从图中可以看出，信号长度为14s左右时，图像重构效果最佳，当信号长度小于14s时，图像质量随信号长度增加而增强，当信号长度大于14s后，随着信号长度的增加，图像质量呈现下降趋势，可见本方法中，光声信号的长度影响着光声图像重构的质量，这是由于利用超声干涉测量法重构的格林函数对应的路径是从成像区直达散射透镜换能器，如图4中红色箭头所示路径1，但当参与重构的探测信号过长时，信号中不但有前行直达成分，还有从周围散射透镜反射的成分，如图4中蓝色箭头所示路径2，这部分信号对应的传播路径未包含在格林函数的计算中，导致图像的畸变.因此，使用本方法重构光声图像时，需选择合适的光声信号长度，以实现最佳成像效果，为了更详细地观察光声初始声场重构效果，对图6 c绘制沿轴及y轴的声压分布曲线，如图7所示，1.0(a)P.(t(b)Max20.54uS(ww)人0.0-0.5-2-1.0Min0510152025-2-1012Time(s)X(mm)2(c)Max40(d)(uw)人10-115210Min-2-1012101214161820X(mm)Signal length(us)图6(a)换能器接收到的光声信号波形；(b)信号长度为25s时得到的光声图像；(c)信号长度为14s时得到的光声图像；(d)PSNR与信号长度之间的关系曲线Fig.6(a)PA signal received by transducers,(b)recon-structed PA images with 25 s signals,(c)reconstructedPA images with 14 s signals,(d)relationship betweenPSNRand signal length1.0(a)1.0(b)0.80.80.60.60.40.40.20.2-2-1012-2-1012Y(mm)Y(mm)1.0(c)1.0(d)0.80.80.60.60.40.40.20.2-2-112-2-112X(mm)X(mm)图7(a)沿x=-1 mm方向的声压分布;(b)沿x=1mm方向的声压分布;(c)沿y=-1 mm方向的声压分布;(d)沿y=1mm方向的声压分布Fig.7(a)Pressure distribution along x=-1 mm,(b)pressure distribution along x=1 mm,(c)pressure distri-bution along y=-1 mm,(d)pressure distribution along y=1 mm从图6 及图7 可以看出，当选择合适的信号长度后，图像重构质量明显优化，光吸收体的形状及位置与实际情况大致吻合，四个光吸收体中心位置声压大致相等.依然存在的少量偏差可能由以下原因导致：（1)格林函数计算还存在一定的误237第2 期何君君等：基于声学散射透镜的低成本光声断层成像方法差；（2)由于对光声信号的长度进行截取，导致尾波中所含的部分信号成分丢失.3结论本文提出了一种基于声学散射透镜的PAT成像方法，开展了理论分析及数值模拟实验，首先，对影响散射透镜格林函数计算精度的相关因素进行了具体研究，结果表明超声干涉测量过程中信号相干叠加区范围的选择以及换能器单元数目对计算精度有直接影响，选择合适的信号叠加区范围，并采用多次测量，可有效提升格林函数计算精度；随后，使用4个散射透镜包围光吸收体同时作用，利用时间反转算法获得了高质量PAT图像；最后，就如何优化图像质量开展了具体研究，结果表明，选择合适的信号长度可避免反射信号的干扰并优化图像质量.由于散射透镜在光声信号探测过程中起到了增大换能器阵列等效数值孔径的作用，可以显著减少PAT所需的换能器数量；同时，本文提出的散射透镜制造简单，无须特殊材质或精密加工，综上所述，本文提出的方法可起到降低PAT系统硬件成本的作用.参考文献1Wang L V,Yao J J.A practical guide tophotoacoustic tomography in the life sciences.NatureMethods,2016,13(8):627638.2Upputuri P K,Pramanik M.Recent advances towardpreclinical and clinical translation of photoacoustictomography:A review.Journal of Biomedical Optics,2016,22(4):041006.3殷杰，陶超，刘晓峻.多参量光声成像及其在生物医学领域的应用.物理学报，2 0 15,6(9)：0 9 8 10 2.（YinJ,Tao C,Liu X J.Multi-parameter photoacousticimaging and its application in biomedicine.ActaPhysica Sinica,2015,6(9):098102.)4Xi L,Grobmyer S R,Wu L,et al.Evaluation ofbreast tumor margins in vivo with intraoperativephotoacoustic imaging.Optics Express,2012,20(8):87268731.5Kang J,Kim E K,Young Kwak J,et al.Optimal laserwavelength for photoacoustic imaging of breastmicrocalcifications.Applied Physics Letters,2011,99(15):153702.6Shang S S,Chen Z J,Zhao Y,et al.Simultaneousimaging of atherosclerotic plaque composition andstructure with dual-mode photoacoustic and opticalcoherence tomography.Optics Express,2017,25(2):530539.7李辉，曹宇，刘红波，等.高速12 8 通道小动物多光谱光声断层成像系统.激光技术，2 0 17,41(5)：6 6 9一674.(Li H,Cao Y,Liu H B,et al.High-speed and128-channel multi-spectral photoacoustic tomographysystem for small animal.Laser Technology,2017,41(5):669674.)8Li L,Zhu L R,Ma C,et al.Single-impulsepanoramic photoacoustic computed tomography ofsmall-animal whole-body dynamics at highspatiotemporalresolution.NatureBiomedicalEngineering,2017,1(5):0071.9Xi L,Jiang H B.High resolution three-dimensionalphotoacoustic imaging of human finger joints in vivo.Applied Physics Letters,2015,107(6):063701.10 Rajian J R,Fabilli M L,Fowlkes J B,et al.Drugdelivery monitoring by photoacoustic tomographywith an ICG encapsulated double emulsion.OpticsExpress,2011,19(15):14335-14347.11 Na S,Russin J J,Lin L,et al.Massively parallelfunctional photoacoustic computed tomography of thehumen brain.Nature Biomedical Engineering,2022(6):584592.12 Guan T P,Shang W T,Li H,et al.From detectionto resection:Photoacoustic tomography and surgeryGuidance with indocyanine green loaded goldNanorodliposome core-shell nanoparticles in liverCancer.Bioconjugate Chemistry,2017,28(4):1221-1228.13 Wu D,Yang J G,Zhang G,et al.Noninvasive in vivomonitoring of collagenase induced intracerebralhemorrhagebyphotoacoustictomography.Biomedical Optics Express,2017,8(4):2276-2286.14 Wang X D,Pang Y J,Ku G,et al.Noninvasive laser-induced photoacoustic tomography for structural andfunctional in vivo imaging of the brain.NatureBiotechnology,2003,21(7):803-806.15 Yang S R,Qin W,Guo H,et al.Design andevaluation of a compound acoustic lens forphotoacoustic computed tomography.BiomedicalOptics Express,2017,8(5):2756-2765.16 Song C L,Xi L,Jiang H B.Liquid acoustic lens forphotoacoustic tomography.Optics Letters,2013,38(15):2930-2933.17 Brown M D,Zhang E Z,Treeby B E,et al.杨贞）（责任编辑238第5 9 卷南京大学学报（自然科学）Reverberant cavity photoacoustic imaging.Optica,2019,6(6):821822.18 Fink M.Time-reversal acoustics in complexenvironments.Geophysics,2006,71(4):SI151-SI164.19 Yin J,Tao C,Liu X J.Dynamic focusing of acousticwave utilizing a randomly scattering lens and a singlefixed transducer.Journal of Applied Physics,2017,121(17):174901.20 Wapenaar K,Draganov D,Snieder R,et al.Tutorial onseismic interferometry:Part 1-Basic principles andapplications.Geophysics,2010,75(5):75A195-75A209.