光声成像中的伪影及噪声去除研究硕士学位论文.pdf

西安电子科技大学 硕士学位论文 光声成像中的伪影及噪声去除研究 姓名：张雄斐 申请学位级别：硕士 专业：生物医学工程 指导教师：田捷 201201 塑墨 i l li lrl 1 lll l l ll l li i i i i ii f 一 _ _ - _ _ - _ - _ _ _ _ _ - _ - _ _ _ - - o _ _ _ _ - - _ - - _ - _ _ _ _ _ _ _ _ _ - - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ - - _ - _ _ _ - _ - _ _ - - _ _ _ - _ _ _ _ _ - - _ _ - _ _ _ o - 一 1 _ - _ 一 y 2 0 6 8 3 15 摘要 光声成像技术( p h o t o a c o u s t i ci m a g i n g ) 是一种基于光声效应的新颖的成像方 式。由于在生物组织中，病变结构与周围结构吸收差异较大，而光声成像技术的 任务就是对吸收分布进行重建，所以光声成像技术可以为疾病的判断提供一定的 依据。光声成像技术同时具备了光学成像分辨率高及超声成像成像深度深的优点， 且为无损成像方式，因而已经成为生物医学成像领域的一个研究热点。 光声成像的研究包含前向问题的研究、逆向问题的研究以及成像系统的研究， 前向过程就是组织吸收激光发出光声信号的过程，逆向问题也即重建问题，就是 在计算机中对探测器接收到的光声信号进行重建的过程，成像系统的研究就是对 激励源、探测器及系统其它组成部分优化的研究。 本文主要集中于逆向过程算法的优化以及对光声信号去噪的研究。考虑到滤 波反投影的重建结果存在伪影，因而将各项异性扩散理论引入重建中，提出基于 各项异性扩散理论的光声成像重建方法，使得伪影的影响大大减小。同时本文模 拟了实际光声信号中可能存在的噪声并对利用现有的滤波器对噪声进行处理，之 后提出级联滤波器并对混合噪声进行处理，得到了较好的去噪效果。 关键词：光声成像技术伪影各项异性扩散级联滤波 光声成像中的伪影及噪声去除研究 a b s t r a c t a b s t r a c t p h o t o a c o u s t i ci m a g i n g ( p a i ) i san o v e li m a g i n gm o d eb a s e do nt h ep h o t o a c o u s t i c e f f e c t s i n c et h el e s i o n ss t r u c t u r e sa n dt h eh e a l t h yp a r t sa r eg r e a t l yd i f f e r e n tf r o mt h e a b s o r bd e g r e eo fl a s e ri nt h eb i o l o g i c a lt i s s u e ，p a it e c h n o l o g yc a np r o v i d eab a s i c j u d g m e n to ft h ed i s e a s ef o rt h et a s ko fi t i s t or e b u i l dt h ed i s t r i b u t i o no ft h ea b s o r b p a r t s p a lt e c h n o l o g yh a sb o t ht h ea d v a n t a g eo fh i 曲r e s o l u t i o nf r o mo p t i c a li m a g i n g a n dd e 印i m a g i n gd e p t hf r o mu l t r a s o n i ci m a g i n g ，m o r e o v e r ，i ti san o n i n v a s i v e i m a g i n gm o d e ，s oi ta t t r a c t sm o r ea n dm o r ea t t e n t i o na n dh a sb e c o m ean e wh o t s p o ti n t h ef i e l do fb i o m e d i c a li m a g i n g t h er e s e a r c ho fp a if o c u so nt h r e ea s p e c t s ，t h ef o r w a r dp r o b l e m ，t h ei n v e r s e p r o b l e ma n dt h ei m a g i n gs y s t e m t h ef o r w a r dp r o b l e mi st h ep r o c e s st h a tt h et i s s u e a b s o r b st h el a s e ra n de m i t st h ep h o t o a c o u s t i cs i g n a l ，t h ei n v e r s ep r o b l e ma l s oc a l l e d r e c o n s t r u c t i o np r o c e s si st h ep r o c e s st h a tp h o t o a c o u s t i cs i g n a lb er e c e i v e db yt h ep r o b e a n db eu s e dt or e b u i l dt h ea b s o r p t i o nd i s t r i b u t i o nb yt h ec o m p u t e r , t h er e s e a r c ho f i m a g i n gs y s t e mc o n c e n t r a t e so nt h eo p t i m i z i n go ft h es y s t e mp a r t ss u c ha st h el a s e r s o u r c e ，t h ep r o b ea n ds oo n t h i sp a p e rm a i n l yf o c u s e so nt h eo p t i m i z a t i o no ft h ea l g o r i t h mi nt h ei n v e r s e p r o b l e ma n dt h ed e n o i s i n go ft h ep h o t o a c o u s t i cs i g n a l w h e ns t u d y i n gt h ef i l t e r b a c k - p r o j e c t i o nr e c o n s t r u c t i o na l g o r i t h me x i s t e d ，w ef o u n dt h er e c o n s t r u c t i o nr e s u l t s c o n t a i n sa r t i f a c to fs o m ed e g r e e ，t h e nw ei n t r o d u c et h ea n i s o t r o p i cd i f f u s i o nt h e o r y i n t or e c o n s t r u c t i o n a l g o r i t h m ，a n dp r o p o s et h ea l g o r i t h mf i l t e rb a c k - p r o j e c t i o n r e c o n s t r u c t i o na l g o r i t h mb a s e do nt h ea n i s o t r o p i cd i f f u s i o nt h e o r y , t h i sn e w a l g o r i t h m h a sb e e np r o v e dt og r e a t l yr e d u c et h ee f f e c to fa r t i f a c t f u r t h e r m o r e ，w ec o n s i d e rt h e n o i s et h a tm a yb ec o n t a i n e di nt h ep h o t o a c o u s t i cs i g n a l ，w eu s es o m ef i l t e r sa l r e a d y e x i s t e dt od e a lw i t ht h en o i s ea n df i n dt h er e s u l tn o tg o o de n o u g h ，s ow ep r o p o s et h e m e t h o do fc o m b i n i n gs o m ef i l t e r st o g e t h e ra l s oc a l l e dt h ec a s c a d ef i l t e r , h e r ew e c o m b i n et h r e ek i n d so ff i l t e rt o g e t h e ra n du s et h e mt of i l t e ro u tt h en o i s e ，f i n a l l yw e g e tam u c hb e t t e rr e s u l to fd e n o i s i n g k e y w o r d s ：p h o t o a c o u t i ci m a g i n g a r t i f a c t a n i s o t r o p i cd i f f u s i o ni n t e g r a t e d 矗l t e 鹅 光声成像中的伪影及噪声去除研究 第一章绪论 第一章绪论 1 1 前言 自1 8 9 5 年德国物理学家伦琴发现x 射线之后，第二年就提出把x 射线用于 治疗的设想，在此之后的百年中，x 射线在医学领域发挥了巨大作用，之后陆续 有相关的成像技术以及设备出现，现在比较常用的技术有x 射线成像，c t 成像， m r i 成像，p e t 成像以及超声成像等。这些影像技术现已成为医疗诊断中不可缺 少的手段。这些技术，既能降低诊断成本，又能为病患提供更准确且廉价的诊断 服务，而且对于疾病的及早发现以及及早治疗提供了保证。下面简单介绍一下上 述影像技术。 x 射线成像是作用时间长，对人体危害比较大的一种成像方式。x 射线利用 其穿透性以及荧光效应能在荧屏或者胶片上成像，应用到人体组织时，由于不同 组织结构不同，导致其密度以及厚度有一定的差别，这样在x 射线通过不同组织 时，被吸收的程度也就不同，到达在荧屏或者胶片上的x 射线剂量也就不同，这 样就形成了有差别可以分辨的灰阶图像。但是由于成像的x 射线剂量比较大，对 人体的危害也比较大，同时，由于成像技术本身的限制，x 射线成像的对比度比 较差，噪声以及伪影都比较严重。 c t 成像技术也被称为电子计算机断层扫描成像技术，它利用灵敏度高的探测 器在精确准直的x 射线束照射下，围绕着人体的某一个部位进行一个断面接着一 个断面地扫描，然后通过计算机对探测器接受到的x 射线的衰减信息进行处理， 组成图像矩阵，得到灰阶图像。与x 射线成像相比，c t 成像技术排除了重叠干扰， 同时分辨率也有极大地提高，而且对软组织和实质性器官能进行比较有效地成像。 但是c t 成像技术也与x 射线成像技术类似，对病患危害比较大【1 1 。 m r i 成像技术是利用设备产生的磁场，同时在射频脉冲的照射下，使得置于 磁场中的物体释放出电磁波，通过接受到的电磁波确定物体原子核的位置和种类， 从而准确绘制物体内部的立体图像。m r i 成像技术对软组织的分辨力优于c t 成像 技术，同时对人体没有电离辐射，是一种比较安全的成像技术。但是，m r i 成像 技术对肺部的成像能力比较差，扫描时间长，空间分辨力也不够理想而且价格昂 贵。 p e t 成像技术也即正电子发射型计算机断层成像技术，是当下比较先进的临 床成像技术。它通过在人体代谢所必需的诸如葡萄糖、核酸等物质上标记短寿命 的放射性核素制作成显像剂注入人体，在不同组织不同代谢情况下，标记的核素 2 光声成像中的伪影及噪声去除研究 分布也不同，在核素衰变的过程中会释放出正电子，通过高灵敏度的照相机捕捉， 然后在计算机中进行散射和信息校正，通过对正电子分布的分析，就可以得到生 物体内的三维影像。p e t 成像技术灵敏度高，特异性高，能够进行全身显像同时 安全性也很好，但是p e t 成像技术对于肿瘤和炎症较难分辨，同时需要操作者有 比较丰富的经验而且价格昂贵i l j 。， 超声成像技术是基于人体不同组织声阻抗的差异，利用超声波照射人体，通 过接收超声波在不同声阻抗界面的放射波来进行成像的。超声成像技术能够实时 成像且可以便携化同时没有副作用，但是超声成像技术穿透性比较差，导致成像 深度有限，同时分辨率也不太理想。 从上面的描述我们不难发现，现在的医学成像技术得到了很大程度的改进， 对人体的危害也降低了很多。但是这些成像技术仍然在不同方面有不同程度上的 缺陷，仍然有很多问题亟待解决i 比如分辨率低下以及成像深度不足等等，这就 需要人们开发出更好的成像技术。光声成像技术就是在这种需求下应运而生。 1 2 光声成像技术概述 亚历山大贝尔于1 8 8 0 年在固体中首先发现了光声效应，同时将这种效应报告 给了美国科学进展协会【2 】。但是由于当时技术条件的限制，这种效应并没有得到有 效利用，此后的一百年间，随着各种强光源以及弱信号检测技术的出现，光声效 应又开始重新获得人们的重视，并得到了广泛应用。2 0 世纪8 0 年代，t h e o d o r e b o w e n 等人提出了光声成像技术的概念【3 】【4 】。此后光声成像技术得到越来越多人的 重视，并取得了巨大的进展。 1 2 1 光声成像原理 光声成像是基于光声效应的一种成像技术。所谓光声效应，指的是介质在短 脉冲激光或者调制光源照射的时候，由于内部吸收光导致温度改变从而引起某些 结构和体积的变化，比如体积的涨，而导致形成超声波并往外辐射的一种效应。 光声效应本质上属于能量转换。 在光声成像中，使用持续时间很短的激光脉冲照射组织。此时产生的超声波 的幅度值和被组织吸收的能量成正比，而组织吸收的能量是由光吸收分布决定的。 在组织中，生理特性的变化会导致光吸收物质分布的变化，所以接收到的超声波 就能反映组织中生理特性的变化，进而对病变区域进行成像。光声成原理如图1 1 所示： 第一章绪论 1 2 2 光声成像的数学描述 t r a n s d u c e r s 图1 1 光声成像原理示意图四 仑 8 弓 罂 墨 从光声成像的原理描述中可知，光声成像是基于生物组织对光的吸收，转换 为热能引发膨胀而产生超声波的。而热能的产生是由入射激光脉冲的功率及组织 对其的吸收系数共同决定的。如果用s ( r , o 表示热源，表示吸收系数，l ( r ，t ) 表示 脉冲激光的强度，那么吸收过程可以用公式描述如下： s ( r ，0 = g a i ( r ，t ) ( 1 1 ) 光声成像的基本理论源自于热传导理论，如果激光脉冲的脉宽比热扩散时间 小的多，那么热扩散及可以被忽略，在这种情况下，热传导方程如下： p c p 杀t ( r ，t ) = h ( r ，t ) ( 1 - 2 ) 其中，p 表示密度，c p 表示比热容，t ( r ，t ) 表示吸收脉冲而升高的温度，h ( r , o 表示单位时间单位体积在热源处的加热函数。 本文在研究中主要关注的是在声性质相对均匀的组织中的传播，在声性质均 匀的介质中，声传播由表达式( 1 3 ) 及表达式( 1 4 ) 共同描述【6 】： p 雨u - - u ( r ，0 = - v p ( r ，0 ( 1 - 3 ) v “r ，t ) 一掣+ p t ( r ，t ) 0 - 4 ) 其中，d 表示各向异性体积扩散系数，c 表示声速，u ( r ，t ) 表示声分布，p ( r ，t ) 表 示声压。 将式( 1 2 ) 、( 1 3 ) 及( 1 4 ) 表达式整合在一起，可得到由热源h ( r ，t ) 产生的声压 p c r , t ) f ( ，表达式如下m ： v z p ( r j t ) 一刍薯p ( r ，t ) ：一导袅h ( r ，t ) ( 1 5 ) 4 光声成像中的伪影及噪声去除研究 上述表达式通过格林函数求解【8 l 【9 】可以得到常用的表达式如下： 毗d = 去高掣 m 6 ， 其中t 由如下表达式得到： t ，= t 一必 ( 1 7 ) 在实际应用中，热函数h ( r ，t ) 可以用空间吸收函数a ( r ) 和光强度函数i ( t ) 来表 示，表达式如下： h c r ，t ) = a ( r ) i ( t )( 1 - 8 ) 将式( 1 - 8 ) 式带入式( 1 6 ) 式中，可得声压p ( r ，t ) 的表达式如下： p ( r ，t ) = 去高a ( r ，) i ，( t 7 ) ( 1 9 ) 其中 i ，( t ，) = 掣 ( 1 1 0 ) 1 2 3 光声成像的研究进展 光声成像结合了光学成像和声学成像的优势【1 0 1 。光学成像的优势在于选择性 吸收，就是组织中不同的结构由于成分不同，对照射光的吸收也有很大差异。以 癌变组织为例，由于癌变过程中，组织在进行剧烈地代谢，需要的各种吸收物质 也比较多，所以对光的吸收也更强，因而能很容易地和正常组织区分开来。声学 成像的优势在于高穿透性，就是在组织中传播距离更远，成像深度更深。具体来 说，光声成像是一种无损的成像技术；同时其图像分辨率高，可达到微毫米量级； 而且由于光声成像是基于组织对光的选择性吸收，所以光声成像能在一定程度上 反映组织的结构信息同时对比度较高；此外，光声系统能够实现小型化而且成本 比现行的成像技术更低。因而自2 0 世纪8 0 年代t h e o d o r eb o w e n 等人提出了光声 成像技术的概念至今，光声成像技术得到越来越多人的重视，并取得了巨大的进 展。 目前国外比较有代表性的光声科研小组有：美国华盛顿大学圣路易斯分校光 学成像实验室的l i h o n gv w a n g 【l l 】、美国普渡大学健康科学学院r o b e r t a k r u g e r 1 2 1 、奥地利维也纳大学计算机科学中心光声成像实验室m a r k u sh a l t m e i e r1 1 3 】、英 国伦敦大学物理与生物工程系光声成像实验室p a u lc b e a r d 及b e nc o x 1 4 】、德国慕 尼黑大学v a r i e sn t z i a c h r i s t o r s 【1 5 1 、荷兰屯特大学e e md em u l 【1 6 】以及瑞士伯尔尼 大学m a r t i nf r e n z 【l ”等小组。 第一章绪论 国内比较有代表性的光声科研小组有：华南师范大学激光生命科学研究所邢 达【嵋j 、台湾大学电机工程研究所超声波影像实验室李百祺【1 9 】等小组。 经过科研人员的努力，光声成像技术从基本的成像算法及成像系统研究已经 开始向i 临床方向发展，美国华盛顿大学圣路易斯分校l i h o n gv w a n g 团队研发的与 超声成像系统结合的光声成像系统成像深度已达2 5 c m ，可以开展临床研究【2 0 】。 o p t o s o n t i 公司制造的光声成像系统成像时间已经缩短到1 2 秒，并且能够进行预临 床研究【2 1 1 。相关的设备公司f a i r w a ym e d i c a lt e c h n o l o g i e s 公司【2 2 1 、z o n a r em e d i c a l s y s t e m s 公司【2 3 1 、e n d r a 公司【2 4 】、s e n om e d i c a li n s t r u m e n t 公司【2 5 1 等，也都研制出 了具有各自特色的光声成像设备。 1 3 本文研究内容 尽管经过诸多国内外科研人员的努力，光学成像技术已经接近临床实践，但 是光声成像在成像算法上始终有些不足，而且重建的结果还有伪影及噪声等因素 的影响。本实验室准备在重建的方法以及对结果的处理上做一些工作，然后搭建 一套相关系统并结合实验室已有荧光发光系统构建一个多种功能的大系统。 本文第一章首先介绍了已经应用到临床的一些成像方法，然后对光声成像的 概念、原理以及研究进展进行了简单介绍并说明了光声成像的优势。 第二章介绍了光声成像的前向以及逆向问题所研究的内容，并对逆向问题的 重建算法进行了描述。 第三章提出了一种基于各向异性理论的滤波反投影法改进算法。由于滤波反 投影法本身的局限，伪影在重建的结果中不可避免，为了减小伪影的影响，提高 图像重建的精度，本文将各向异性理论应用到滤波反投影法中，有效地改善了重 建结果。 第四章考虑到实际的成像系统中会有各种噪声对光声成像造成影响，所以在 计算机仿真中对各种可能影响结果的噪声进行了模拟，并使用不同的滤波器进行 仿真处理进而提出组合滤波的方法，使得计算机上的仿真模型能更好地贴近实际 成像系统。 第五章为总结与展望，对前面几章进行总结，并指出其中的问题，提出可以 进一步研究的内容及方向。 6 光声成像中的伪影及噪声去除研究 第二章光声成像前向过程及重建算法的介绍 第二章光声成像前向过程及重建算法的介绍 2 1 光声成像前向过程 描述照射到组织的光被组织吸收，发出声波，进而被探测器接收的整个过程 称为光声成像的前向过程。顾名思义，前向过程的仿真就是对光被吸发出声波并 被探测器接收得到光声信号的过程的仿真，通过建立前向模型可使对其的研究更 加方便。下面将从光的吸收到发出声波以及被探测器接收得到光声信号的逐一进 行描述，进而得到前向模型。 2 1 1 光声信号的产生 光照射到生物组织上，与生物组织的相互作用既有吸收，还有反射和散射， 光声成像中利用的是光与生物组织的热相互作用，即生物组织对光的吸收所导致 的生物组织能量增加进而产生热量的作用，光声效应就是热相互作用的产物，所 以本文只考虑生物组织对光的吸收作用。 描述生物组织光吸收的一个重要参数就是光的吸收系数，所谓光的吸收系数 指的是单位长度长一个光子被吸收的几率，其描述符号为k ，常用单位为c m 。 用公式描述如下： d l1 = 忑t 1 2 - 0 上式中，i 描述的是被吸收前光子的能量，d i 描述的是光传输过程中单位路径 内光子由于被吸收而减少的能量，d ) 【描述的是单位路径。上式描述的是生物组织 对光吸收能力的强弱程度。 光被生物组织吸收后，生物组织能量增加。这部分被生物组织吸收的能量可 以用以下公式描述： a ( r ，t ) = i ( r ，t )( 2 2 ) 上式中，a ( r ，t ) 描述的是被生物组织吸收的能量，k 描述的是光的吸收系数， i ( r ，t ) 描述的是光的入射强度。上式说明在光声成像中，生物组织中热量的产生是 由组织的光吸收系数以及光的入射强度共同决定的。 光被生物组织吸收，生物组织能量增加，就会导致局部温度升高并产生绝热 膨胀，组织间的膨胀收缩导致相互挤压，这种挤压从光吸收部位向其他组织部位 传递，这个过程实际上就是声波的传递过程。 田 光声成像中的伪影及噪声去除研究 在光声成像中，声波主要的描述物理量为声压和声速。声压用p ( r t ) 表示，其 中r 为空间变量，t 为时间变量，也即声压是空间与时间的函数。声速用c 表示， 单位为m s ，在生物组织的光声成像中，不同的组织的声速有差异，声速随组织密 度的增大而增加，为了研究方便，一般把声速c 当作常量，即1 5 0 0 m s 。 经过第一章中数学描述的推导可知，具体光声信号的描述公式如下： p ( r ，t ) = 去f f f d s r ，a h ( r , t ) i t 刮r - r | l 。 ( 2 - 3 ) 其中 h ( r ，t ) = a ( r ) l 。( t ) ( 2 - 4 ) 2 1 2 光声信号的接收 实际应用中，需要用探测器来接受光声信号，本文中主要研究的是单个探测 器的情况。所用的探测器，一般使用超声探测器，单个探测器可以使用水听器。 所谓水听器，是把水下声信号转化为电信号的换能器。之所以要使用水听器，是 因为在实际探测中，为了与重建算法中探测器与生物组织有一定距离且探测器与 生物组织之间的声速保持不变的要求保持一致，需要在生物组织与探测器之间添 加耦合剂，而通常使用的耦合剂为水，因为在声速在水中与生物组织中近似相等， 所以可以使用水听器。 常用的水听器如图2 1 所示： 图2 1 常用水听器示例图 图2 1 所示的水听器为丹麦i 也s o n 公司所生产型号为t c 4 0 1 4 的宽带球形水 听器。水听器的基本技术参数包括频带带宽、接收面积以及探测灵敏度。 此外，由于实际应用中，由探测器接受到的电信号并非是直接可用的光声信 号，而是光声信号与探测器频率响应函数的卷积，所以探测器带宽的选择就尤为 重要，当探测器的带宽范围与光声压谱范围基本吻合的时候，重建出来的图像效 第二章光声成像前向过程及重建算法的介绍 9 果比较好，因为这时重建损失的频率成分会更少闭。 公式( 2 3 ) 描述的是光声信号的通用公式，在声学特性均匀的介质中，如果使 用理想的激光脉冲进行作用，对于球形及圆柱状的吸收体，有如下波g ； 2 7 j ： 图2 2 圆柱状吸收体理论光声信号示意图 上图表明产生的光声信号的形状和光能量的吸收分布相关。 、 g j d i e b o l d 等人给出了对光声信号更为详尽的描述，详情可阅读参考文献 2 7 1 。 2 1 3 光声信号的模拟 毫 为了方便进行研究，需要对光声信号进行模拟仿真，就是对公式( 2 3 ) 进行模 拟仿真。公式( 2 3 ) 实质上是一个偏微分方程，而常用的解偏微分方程的数值方法 有有限微分法1 2 8 】、有限元法1 2 9 】和边界元法【3 0 1 。尽管在很多情况下，这些方法非常 有效，但是对于高频波的时间域模型仿真，这些方法变得慢而且难以处理。因为 这些方法要求每个波形都有很多网格点，然后还要求使用尽可能小的时间步长来 减小不必要的数值离散。有限微分方法中，对模拟过程的时间步骤有要求，时间 步骤太短了容易造成模拟的不稳定性，而且对三维情况。其计算量也非常大，有 限元法以及边界元方法中网格剖分等步骤复杂而且模拟时间较长。 b t c o x 等人将p s m ( p s e u d o s p e c t r a lm e t h o d ) 及k - s p a c e 方法应用到光声前向 模拟中，p s m 可以有效地解决需要很多网格点的问题，k - s p a c e 方法可以有效地解 决时间步长问题1 3 。 ， 在传统的有限微分策略中，区域的微分是通过对网格点值进行线性插值获得 的。使用高阶多项式能获得更多的网格点和更准确的微分值，但是，越多的点， 计算就越复杂。p s m 则将傅立叶序列添加到其中，由于傅立叶成分的幅度可以通 过快速傅立叶变换来计算，而且基函数是正弦函数，每个波长只需要两个点，所 以p s m 方法可以有效地减少计算时间同时减少网格点的使用。 所谓k - s p a c e 就是k 空间，也被称为傅立叶空间或者空间频率空间或者原始数 1 0光声成像中的伪影及噪声去除研究 据空间。k 空间是一个以空间频率为坐标轴的空间坐标系对应的抽象的频率空间。 所谓的空间频率指的是在空间某一方向上单位长度内波动的周期数，单位为h z c m ， 而在一般应用中所谓的频率却是时间频率，即单位时间内波动的周期数，单位为 h z 。空间频率是一个矢量，一般使用类似三维坐标中x ，y ，z 三个相互垂直的分 量来描述，w 胱t q z x 、k y 、k z 来表示，则kw 用从k ( 2 - 5 ) 描述如下： ( 2 5 ) 上式表明，k 空间每一点都a m x 、k y 、k ：三个量来决定，都代表一个固定点 空间频率k ，从上式可以看出，在空间的原点处，频率较低，与三维矢量情况类 似，与原点处空间距离相等的点，具有相同大小的空间频率，但是其方向不一定 相同。 对分布于k 空间的信号而言，k 空间的某一特定点的信号空间频率相同，其 数据的大小则表示信号的强弱，进而言之，对于分布于k 空间的波形，k 空间的 每个数据就代表了一列固定空间频率的波。若对空间频率连续的波进行二维傅立 叶变换，就将得到一个二维空间频率矩阵，矩阵中每一个点对应一定的空间频率 信号，即所谓的二维k 空间。 声波方程可以使用k - s p a c e 方法进行处理，对于声速及密度均匀的介质c 和p 均 为常数，声波方程可用公式( 2 6 ) 描述g a t t 3 2 】： ( 嘉_ c 5 v 2 ) 叫装 一 其中，c 代表声速，p 代表密度，h ( r ，t ) 表示单位时间单位体积在热源处的加热 函数。其中r 为c 2 毒：p 表示各向异性体积扩散系数，c p 表示比热容。 对式( 2 6 ) 进行空间傅立叶变换，并利用格林函数进行求解，可得表达式如下： p ( k ，t ) = f h x j rg ( t _ t ，) 面d h t , “( 2 - 7 ) k 为k 空间中波向量l 【= ( k x ，k y ) 的模，s ( t ) 为r n x 等，其中，格林函数g ( k t ) 的 表达式如下： r 0，t0 g ( k ，t ) = _ - s i n ( c o k t ) ，t 0 ( 2 - 8 ) f c o k 当h t ( t ) = 8 ( 0 ，且瞬时脉冲在t = o 时刻，式( 2 - 7 ) p - - j 以修正为如下表达式： p ( k t ) = f i t xc o s ( c 0 k t )( 2 9 ) 第二章光声成像前向过程及重建算法的介绍 i i 那么根据在t = 0 时刻的式子，可以得到声场中t 时刻的声压表达式如下： p ( x ，0 = f y “ 1 t h x ( k ) c o s ( c o k t ) ) ( 2 1 0 ) 假设时间步长为t ，那么时刻t a t 与时刻t + a t 及二者组合的表达式表示如下： p ( t + a o = f h x ( c o s ( c o 姆) c o s ( c o k a o s i n ( c o k t ) s i nc o s ( c o k a 0 ( 2 - 11 ) p ( t 一t ) = r h x ( c o s ( c o k t ) c o s ( c o ka o + s i n ( c o k t ) s i nc o s ( c o ka t )( 2 - 12 ) p ( t + a t ) + p ( t a t ) = 2 p ( t ) ( 1 2 ( s i n ( c o k a t 2 ) ) z ) ( 2 1 3 ) 再进一步处理可得如下表达式： p ( t + a t ) + p ( t a t ) 一2 p ( t ) = - 4 ( s i n ( c o k a t 2 ) ) z ( 2 1 4 ) 上式就是k - s p a c e 方法得到的时间步长策略，可以应用上式进行传输模拟。 假设有仿体如图2 3 所示，且其中声速c 的值设定为15 0 0 m s ，密度p 的值设定 为1 0 0 0 k g m 3 ，光分布的强度设定为2 。 图2 3 仿体截面吸收分布图 按照k - s p a c e 的要求划分网格，并设定时间步长，在某处探测器接收到的波形 如下图： 图2 4k - s p a c e 方法仿真所得光声信号图 1 2 光声成像中的伪影及噪声去除研究 将图2 4 所示仿真所得的光声信号图与图2 2 所示的理论光声信号示意图以及 文献 2 7 】中研究结果对比，可知上述k - s p a c e 能将传输结果准确的模拟出来。 2 2 光声成像重建算法简介 光声成像的逆向过程，就是描述探测器接收光声信号，而后放大器将所接收 到的光声信号放大，示波器对其进行采样等预处理后信号输入计算机，并使用一 定的重建算法，对光声信号进行反演重建，进而得到初始光吸收分布的整个过程。 光声重建算法在不同的探测器扫描方式下有不同的类别。现在常用的探测器 采集方式有使用聚焦探测器扫描、非聚焦探测器扫描以及超声阵列探测器扫描。 使用聚焦探测器扫描，与b 超扫描类似，可以直接进行成像，不需要成像算法。 l i h o n gv w a n g 研究小组使用聚焦探测器扫描，研制出了类似显微镜的小动物成像 系统p a m ( p h o t o a c o u s t i cm i c r o s c o p ys y s t e m ) 3 3 】。尽管这类系统有非常好的分辨率， 文献 3 3 】中提到的p a m 系统轴向分辨率为1 5 i _ t m ，横向分辨率为4 5 i t m ，但是成像 深度有限，为3 m m 。系统中轴向分辨率由入射激光的脉冲与探测器的频率响应特 性共同决定，而横向分辨率则由探测器的聚焦点直径决定。由于超声衰减系数和 超声频率成正比，系统的轴向分辨率又由探测器的频率相关，为了获得较好的轴 向分辨率，就需要提高探测器的频率，但是这样超声的衰减就会增加，成像深度 就受限制，也即在这种系统中，分辨率的要求会限制成像深度。为了能够接近临 床研究，需要增加成像深度，而使用非聚焦探测器及超声阵列扫描能实现这一点， 所以当前研究人员的研究所使用的主要扫描方式为非聚焦探测器扫描方式及超声 阵列扫描方式。同时这种非聚焦的扫描能获得更大角度范围内的光声信号，此外 还可被视为点探测器，这些优势可以为研究重建算法提供有利条件。 在使用非聚焦探测器的光声成像领域里，研究人员已经研究出非常多的重建 算法，已有的重建算法中比较有代表性的有：基于逆r a d o n 变换的滤波反投影法【蚓； 基于解析解的时域反投影法1 6 】以及基于傅立叶变换的后向投影法1 3 s 。下面将按照时 域及频域对上述方法进行划分并介绍。 2 2 1 时域重建方法 在光声成像中，基于非聚焦探测器的扫描方式一般为圆形扫描，具体如图2 5 所示： 第二章光声成像前向过程及重建算法的介绍 1 3 图2 5 光声成像扫描方式示意图阁 在图2 5 中，虚线的圆弧表示的是光声成像的投影线，这些虚线圆弧都是以p 所在位置为中心，p 所在的位置即探测器所在位置。这与c t 成像中的扫描方式略 有不同，c t 成像中，探测器接收的均为直线投影线。 基于上述扫描方式接收光声信号，光声方程可以做如下处理： v 2 p ( 咖刍嘉p ( r ，t ) = 一芑晏h ( r ，t ) 其中，p 表示各向异性体积扩散系数，c p 表示比热容，h t ) 表示单位时间单 位体积在热源处的加热函数，p ( r ，t ) 表示在位置r ，时间t 的压力。基于格林函数的 解表达式如下： p ( r 加忐高掣k 十l r 一，v 。 倍峋 其中h c r , t ) 为a c r ) l 。( t ) ，而a c t ) 描述的是生物组织吸收的能量密度，如果声波 是由短脉冲8 ( 0 激发的，那么式( 2 1 6 ) 可进一步修正为： ，佃 p e ( r ，t ) = l i e ( t 一曲p 6 ( r ，t ) d t( 2 1 7 ) 其中 p 6 ( r ，t ) = 晏陡肚( r ) d ni i h 协】 ( 2 _ 1 8 ) 式( 2 - 1 8 ) 中，d q 为r 向量对应到r 所在点的角度，p o ( r ) 为初始压力分布。进一 步化简，如图2 5 所示，对于在位置r 0 的点探测器，有表达式如下： p d ( 删= 晏陡肛( r ) d q i 础】( 2 - 1 9 ) 其中，d q 为r 向量对应到r o 所在点的角度。利用式( 2 1 9 ) 就可以根据探测器所 1 4光声成像中的伪影及噪声去除研究 得到的数据p d ( r o ，0 重建出初始压力p o ( r ) 的分布。 为了方便使用逆r a d o n 变换，可以将式( 2 - 1 9 ) 等价变形为： f ( r ，0 = t 4 x ( r o , o d t = p 0 ( r ) d nl 酬= c t ( 2 - 2 0 ) 如果待测目标靠近球形探测面的中心区域，并且其大小远小于探测区域，那 么球面积分可用平面积分近似代替，可以近似使用r a d o n 变换。k r u g e r 等人在文 献【3 4 】中利用r a d o n 变换对上式进行近似处理，首先f c r , t ) 取空间二阶导数，然后 把每个测量位置上的二阶导数进行反投影，之后再将反投影所得到的值进行积分， 就可获取空间吸收分布，进而得到初始压力分布表达式如下： p 。( r ) 一芝1 兀j j f fd r s ；o 【r t o p f d ( r o , t ) + 2 p d ( r 。，t ) 】i t ：l ，一r o l c ( 2 2 1 ) 一般情况下，如果待测目标在探测区域中心或者接近探测区域中心，上述方 法都能得到较好的重建结果，但是如果待测目标远离探测区域中心，或者待测目 标太大，就会出现非常明显的伪影，因为在这些情况下，球面积分不能用平面积 分近似代替。 l i h o n gv w a n g 等通过对光声方程的格林函数解进行解析求解，给出了圆柱及 球扫描方式下光声方程的解析解，对单个探测器及平面超声探测器阵列成像系统 的算法研究有很大的指导价值。 现在给出圆柱扫描示意图如图2 6 所示： 图2 6 圆柱扫描示意图阁 此时假设探测器与样品的距离远大于光声信号的波长，这一点在实际应用中 很容易满足，l i h o n gv w a n g 等人在文献【7 】中给出了详细推导，并得出在圆柱扫描 情况下初始压力分布的表达式如下： m ，z ) = 一去j 5 fd s 0 畦掣i t 川。 ( 2 2 2 ) 第二章光声成像前向过程及重建算法的介绍1 5 上式， 1 = 毒，【n n 。】为权重因子在球面扫描中有示意图如图2 7 所示： lz 厂 、 a f 9 。神 ；广： 笋心一一一 一 ：韵y 7 图2 7 球面扫描示意图阴 也是在探测器与样品的距离远大于光声信号的波长的条件下， 等人经过推导，得到球面扫描情况下初始压力分布表达式如下： a ( r j 咖) 一去f f d s oi 1 t 0 p ( r o , t ) l t = i r - r o i c 上式中，q = 导。 l i h o n gv w a n g ( 2 - 2 3 ) 上述两种方法都是比较经典的方法，但是上述方法要求探测器能够全方位接 收信号，在这种情况下，重建时才能以探测器为圆心进行圆弧反投影得到重建图 像。同时，上述方法要求的数据比较多，需要旋转探测器或者样品来获得多组数 据，因为数据不足将导致反投影的图像伪影很重，这样又会导致扫描时间长，不 利于实时成像，限制了在预临床上的应用。 2 2 2 频域重建方法 k o s t l i 等人在文献【3 5 】等中提出了基于空间二维傅立叶变换的反投影算法。这 种算法是通过对激光触发后每隔一段时间所记录的声强进行傅立叶变换，然后根 据均匀激光激励下，声强分布与空间吸收分布成正比的原理，对空间吸收分布进 行求解，得到初始压力分布。理论上主要经过了三个步骤，第一步是对得到的光 声信号进行空间二维傅立叶变换，第二步是对空间二维傅立叶变换所得结果进行 时间积分，第三步是对积分所得值进行空间三维逆傅立叶变换就可以得到初始压 力分布。 k o s t l i 等人经过推导得到经过空间二维傅立叶变换反投影算法处理的初始压 力分布表达式如式( 2 2 4 ) 所示： 1 6光声成像中的伪影及噪声去除研究 p ( k x ，k y ) = 2 c 4 0 ；o - c 2 k 2f f p ( x ，t ) e x p ( 一i k x x ) c 。s t ) d ) 【d t ( 2 - 2 4 ) 其中， k y = 厣 上式中的p ( k x ，k y ) 可以由信号矩阵经过一次二维傅立叶变换得到，然后初始压 力分布也即初始声强分布p o ( x ，y ) 就可以通过对p ( k x ，k y ) 进行逆二维傅立叶变换得 到。 2 3 小结 本章主要介绍了光声成像的前向形成过程以及其描述方程并对前向模拟的方 法进行了陈述，同时介绍了光声成像的逆向过程构成并对逆向过程中的重建算法 进行了陈述。在本章第一部分内容中，描述了光声成像过程，激光器发出的激光 被组织或仿体吸收，组织或仿体内热膨胀进而发出光声信号，之后这些光声信号 被探测器接收的整个前向过程，同时还给出了对应上述各前向步骤的数学描述， 最后对上述数学描述进行综合整理得出来前向过程传输模型。在本章的第二部分 内容中，首先描述了有限微分及有限元方法对偏微分方程求解的稳定性不足以及 复杂度较高的缺点，然后详细介绍了能较好进行前向模拟的k - s p a c e 方法，并列出 了k - s p a c e 方法的模拟策略，此外还将其应用于实际光声前向模拟中进行仿真，得 到了与理论结果类似的前向模拟结果。在本章的第三部分内容中，首先根据不同 的扫描方式列举出了不同的光声重建方法，然后根据时域及频域的划分，对几种 比较成熟的常用的重建方法进行了描述。本章介绍了基于逆r a d o n 变换的重建