基于米氏理论的蒙特卡洛光学相干层析成像系统模拟研究

1、第27卷第4期August,2005文章编号:1005-5630(2005)04-0049-06基于米氏理论的蒙特卡洛光学相干层析成像系统模拟研究黄亚达,朱红毅,廖然,吕俊󰀁(浙江大学现代光学仪器国家重点实验室光及电磁波研究中心,浙江杭州310027)摘要:光学相干层析成像(opticalcoherencetomography,OCT)是近年来快速发展的一种生物组织高分辨力实时成像技术。由于生物组织的复杂性,故一直没有完美的OCT理论模型。现利用结合米氏理论(Mietheory)的蒙特卡洛算法研究光在生物组织(高散射介质)内的行为,并以此为基础建立了一种OCT理论模型。然后利用

2、该模型模拟重建了多层介质的图像,分析了空间滤波和偏振光对图像质量的影响,获得了较理想的结果。关键词:光学相干层析成像;蒙特卡洛模拟;米氏理论;高散射介质中图分类号:O431.1文献标识码:AOpticalcoherencetomographysimulationwithMonteCarlomethodbasedonMietheoryHUANGYa-da,ZHUHong-yi,LIAORan,LVJun(CenterforOpticalandElectromagneticResearch,StateKeyLaboratoryofModernOpticalInstrumentation,Zheji

3、angUniversity,Hangzhou310027,China)Abstract:Opticalcoherencetomography(OCT)isafastdevelopingtechnologyformedicaldiagnosisinvivo.Buttillnow,therestillisntexistwell-acceptedphysicalmodelsofOCT.ThispaperusesMonteCarlosimulationandMietheorytobuilduptheOCTmodel.Oursimulationbasedonthismodelfittheexperime

4、ntaldataverywell.Thismethodprovidesusagoodandreliablewaytostudytheopticalcoherenceimaging.Keywords:opticalcoherencetomography;MoneCarlosimulation;Mietheory;highscatteringmaterial1引言光学相干层析成像1(opticalcoherencetomography,OCT)是近年来快速发展的一种生物组织高分辨力实时成像技术。目前OCT系统已经成功地应用到青光眼研究,神经外科与神经科学研究,皮肤癌诊断,神经外科指导,癫痫脊椎手术

5、,提供临床牙科诊断图,牙外科,牙龈,牙黏膜疾病诊断,密度物质研究等很多方面2。但是由于光在生物组织中的行为的复杂性,故OCT至今没有一个令人满意的理论模型。生物组织是一种高散射介质。研究光在高散射介质中的传输的光学特性及光强的分布对于许多医学应用都有重要意义3。最早人们想通过蒙特卡洛方法来模拟光在高散射介质中的随机过程,WangLH在󰀁收稿日期:2004-11-29(:(,50光学仪器第27卷1995年提出的光在多层介质中的蒙特卡洛模型(MonteCarloinmultiplayers,MCML)4成功模拟了光子在多层复杂高散射介质中的行为。但MCML只是把光当作一个个的光子包

6、,其中并不包含光的相位、偏振等信息,无法直接应用到OCT系统。后来,SchmittJM和LuQ等分别将扩展Huygens-Frense3,Mie理论等物理模型引入蒙特卡洛方法,成功地解决了这一问题,更精确地描述了光在多层高散射介质中的行为。通过结合米氏理论的蒙特卡洛方法数值模拟了光在分层高散射介质中的行为,并以此为基础建立了OCT理论模型。然后根据该模型模拟实现了对多层介质的层析成像,并讨论了空间滤波,偏振等现有OCT技术对层析成像结果的影响。2理论与模型2.1OCT系统5图1OCT系统示意图光学相干层析成像(OCT)系统是一种类似迈克尔逊干涉仪的具有纵横向扫描功能的系统(图1)。通过使用宽带

7、光源,使得样品反射的光走过的光程必须与参考光走过的光程相等时才相干加强,而其他样品光与参考光光程不相等的光信号则被抑制。通过调节参考光的光程可以选择出从样品不同深度返回的光信号,从而达到层析成像的目的。OCT理论模型中最重要也最难解决的是光在样品中的行为,现采用结合米氏理论的蒙特卡洛模拟方法计算。数值模拟过程中,每个光波走过的光程都被记录下来,通过分别对返回界面时光子走过的不同光程对应的光强的统计,可以获得样品内部不同层面的情况。在模拟中记录的不同的光程即相当于实验中参考臂运动到不同位置时对应的输出,用于统计光程的分辨力即OCT中的相干长度,也就是OCT中的轴向最小分辨力。2.2蒙特卡洛算法蒙

8、特卡洛方法4是一种通过对大量随机过程进行统计获得稳定的物理量输出的方法。由于生物组织是高散射介质,OCT系统中进入样品的每个光子在生物组织中的行为具有一定的随机性,因此通过蒙特卡洛方法追踪并对数以亿计的光子进行统计会在样品表面获得稳定的后向散射光场分布,把这个光场作为OCT系统的样品臂返回光场与参考臂返回光场相干成像,比传统的不考虑散射直接计算样品分层间的反射的OCT模型更加符合实际。蒙特卡洛模拟介质中光传播主要是模拟单个光子在介质中的随机碰撞的过程,记录光子在每个散射点的散射和吸收,通过跟踪大量光子的这种行为可以得到介质中光子传播和散射的统计结果。蒙特卡洛算法把光在高散射介质中的行为分为两个

9、具有一定随机性的过程:传播与散射。在传播过程中两次散射间的传播距离即自由程具有随机性,可随机采样取为L,所有L的平均值需等于光在这种介质中的平均自由程。散射过程中每次散射的散射角󰀁和方位角󰀂也具有随机性,需随机采样,这对角决定了散射后光继续传播的方向。各种不同的蒙特卡洛模拟采用不同的󰀁,󰀂取样方法,最简单的是分别在(0, )和(0,2 )中均匀采样。MCML中对󰀁利用Henyey-Greenstein相函数采样,对󰀂则是(0,2 )中均匀采样。现采用的是利用米氏理论推导出米氏相函数对󰀁

10、;和󰀂采样,不但Henyey-Greenstein相函数已被证明是米氏相函数的在各向同性下的近似5,6,而且米氏理论是有精确解析解的电磁理论,它具有其他理论没有的第4期黄亚达等:基于米氏理论的蒙特卡洛光学相干层析成像系统模拟研究512.3米氏理论米氏理论是计算经小球散射后的电磁场的解析解的理论7。样品是一个悬浮着很多小球的均匀介质。这种物理模型与大多数悬浊液以及OCT仿真实验中经常采用的聚乙烯小球水溶液等高散射介质完全一致,在生物组织中这个假设也被证明是合理的。在模拟过程中,为每个光子增加一个新的参数Stokes量。其表达式:IQUV*EE*+EE8S=EE-EE*EE*+EE

11、*(1)*i(EE-EE)*用Stokes量来表示该光波的偏振状态,通过Stokes量也很容易得到E矢量的大小方向。光波在直线传播过程中并不会改变其偏振态,所以在蒙特卡洛模拟中,只需在散射过程中更新Stokes量,其它过程与MCML是一致。图2散射前后局部坐标系关系每次散射,坐标系绕k方向转过一个󰀂,然后绕e转过一个󰀁(图2),新的Stokes量和原Stokes量的关系:Sn+1=M(󰀁)R(󰀂)Sn这里1R(󰀂)=0cos(2󰀂)0sin(2󰀂)cos(2󰀂)000

12、m33(󰀁)m34(󰀁)00000m34(󰀁)m33(󰀁)10(2)0-sin(2󰀂)00m11(󰀁)m12(󰀁)m11(󰀁)00(3)M(󰀁)是小球的Mueller矩阵9M(󰀁)=m12(󰀁)00(4)由将(3)、(4)式代入(2)式可以得到Stokes的每个分量的表达式,其第一个分量I即光强In+1=m11(󰀁)In+m12(󰀁)Qncos(2󰀂)+Unsin(2

13、83042;)比,所以可以定义󰀁和󰀂的概率密度函数P(󰀁,󰀂)=通过条件概率可得P()/(I(󰀁,󰀂)d!4(5)由(5)式可以知道散射后的光强在空间的分布,易知光子出现在某个方向的概率与这个方向的光强成正1101200C=(6)52光学仪器第27卷所以对󰀁和󰀂的采样可分为两个步骤,先采样󰀁,再在给定采样󰀁的情况下采样󰀂。󰀁和󰀂得概率密度函数分别为P󰀁(󰀁)

14、=P(󰀁,󰀂)d󰀂=22 m11(󰀁)I0(8)(9)P󰀂/󰀁(󰀂)=m11(󰀁)I0+m12(󰀁)Q0cos(2󰀂)+U0sin(2󰀂)即得󰀁和󰀂的概率密度函数f(󰀁)和g󰀁(󰀂)󰀁󰀁f(󰀁)=C󰀁,C=󰀂󰀁g(󰀂)=,Cb

15、3042;C󰀂P(󰀁)d󰀁=P(󰀂)d󰀂󰀁 0󰀂󰀁这里省去了(6)式中的归一化常数,实际计算时不用计算该常数,只要分别对(8)、(9)两式进行归一化处理(10)(11)由(10)式可见f(󰀁)与散射前光的Stokes量无关,只与m11有关,而m11是由介质的性质决定的,所以f(󰀁)中所有󰀁的概率都可以提前算好,而g󰀁(󰀂)则与散射前光的Stokes量有关,所以每次散射都必须重新计算。要对给

16、定概率分布进行采样即解方程p()d=a󰀁0x#,#(0,1)(12)由(12)式可得要采样󰀁和󰀂,分别要解以下两个方程f(x)dx=#g(x)dx=#󰀁󰀂0󰀁(13)(14)󰀂(13)式可以提前计算出#(0,1)和󰀁之间一一对应的表格,对󰀁的采样即通过查表,差值即得,蒙特卡洛算法是一种非常耗时的算法,通过提前计算可以大大减少计算所用的时间。(14)式无法提前计算,每次散射的时候通过Newton-Raphson方法解(14)式可得󰀂。当得

17、到󰀁和󰀂后,即得到光波的散射方向,通过(2)式,可以计算出散射后的Stokes量,在新的方向上,准备下一次散射。3数值模拟与分析根据OCT的原理,使用宽带光源的目的是为了在样品光中筛选出和参考光走过相等光学长度的那部分光,从而达到层析成像的目的。在模拟中,由于通过记录光子走过的光程就可以达到相同的效果,所以没有必要模拟OCT中的宽带光源,直接使用光源的中心波长模拟结果。模拟中,光波波长取800nm,介质折射率1.329,小球折射率1.565,模拟光子数为10个光子。使用一台P42.8G的PC机进行模拟计算。3.1多层介质模拟建立一个26层多层高散射介质,每层厚度

18、为0.2cm,小球直径为1.2m,0.1m交替,即第一层1.2m,第二层0.1m,第三层1.2m以此类推总共26层。模拟结果如图3所示,在每层的边界上信号都有突变,随着深度的增加信号的峰值逐渐衰减,层与层之间的分辨力逐渐减弱,15层之后层与层之间就难以分辨了,这是由于散散次数少的光主要集中在表层,随着深度的增加,光子的散射次数也增加,这种随着深度增加,信号减弱,信号对比降低的效应称为雨帘效应。3.2空间滤波设计一个4层介质(详见表1)。在该模拟中在接受面上设置一数值孔径为0.3的空间滤波器,结果(图4)很明显,没有滤波过的信号层之间的分辨力远远低于滤波后的信号,导致根本无法分辨第三层。空间滤波

19、后分辨力显著增加,下降延,7第4期黄亚达等:基于米氏理论的蒙特卡洛光学相干层析成像系统模拟研究53多次散射的光子,以大角度返回,这部分光子往往降低信号的分辨力。OCT中往往使用光纤端面作为空间滤波器,从而提高信号的分辨力,增加其探测深度。图3一个26层高散射介质的模拟结果横轴表示扫描深度;纵轴表示归一化光强。图4空间滤波的模拟结果虚线表示滤波前的扫描信号;实线表示滤波后的扫描信号;横轴表示扫描深度,纵轴表示归一化光强。表14层介质参数小球直径(m)123412345表25层介质参数小球直径(m)3.3偏振光一个5层介质(详见表2),54光学仪器第27卷去其偏振的特性,所以通过偏振片滤波可以滤掉

20、多次散射的光,而这部分光往往会降低信号的分辨力(详见图5)。图5偏振片滤波的模拟结果虚线表示无偏振的扫描信号;实线表示有偏振的扫描信号,横轴表示扫描深度,纵轴表示归一化光强。4结论通过结合米氏理论的蒙特卡洛方法模拟计算了光在高散射介质中的行为。与以前的蒙特卡洛模拟相比,米氏理论的引入更加贴近生物组织的实际情况;在模拟中Stokes量的引入,可以方便地记录光波的偏振状态;通过对光子光程的分类统计,巧妙地实现了参考臂扫描。使用模型成功反演出了多层复杂介质,对空间滤波和偏振光的作用也与理论分析和现有实现结果非常相符。模型比较完美地模拟了光学相干层析成像(OCT)系统。5参考文献1HuangD.Opt

21、icalcoherencetomographyJ,Science.1991,254(5035),11781181.2SchmittJM.Opticalcoherencetomography(OCT):areviewJ,IEEEJournalofSelectedTopicsinQuantumElectronics,1999,5(4):12051215.3SchmittJM.KnuttelA.ModelofopticalcoherencetomographyofheterogeneoustissueJ.JOptSocAmA,1997,14(6):12311242.4WangLHJacquesSL,ZhengLQ,MCML-MonteCarlomodelingoflighttransportinmulti-layeredtissuesJ,ComputerMethodsandProgrammsinBiomedicine,1995,47(2):131146.5LuQMonteCarlomodelin