混浊物质中光子传输的快速模拟与干涉特性分析.pdf

第 31 卷?第 6 期光?学?学?报Vol.31,No.62011 年 6 月ACTA OPTICA SINICAJune,2011混浊物质中光子传输的快速模拟与干涉特性分析林?林1?高应俊2?金重星2?林惠珍2?罗云瀚21广东医学院生物医学工程教研室,广东 东莞 5238082暨南大学理工学院光电工程系,广东 广州 510632摘要?建立了混浊物质中低相干干涉系统的蒙特卡罗模型,结合相干面光强的高斯分布和低相干光源的相干条件,对干涉信号进行了数值分析。提出在蒙特卡罗计算的基线数据基础上快速模拟其他参数条件下干涉信号的方法。快速模拟结果与蒙特卡罗计算结果进行了比较,结果显示在满足前向散射条件下,两种方法所得模拟信号强度基本吻合。当各向异性因子 g=0.9 时,平均相对误差为 2.37%,在 g=0.5 时,平均相对误差为 4.69%。在深度相关干涉信号的模拟中,两种方法对浅表层干涉信号的模拟取得了一致的结果。实验中使用低相干干涉系统对不同浓度的脂肪悬乳注射剂与印度墨水混合液中深度为 0 1 mm 的漫反射光信号进行了干涉测量,并利用基线数据和快速蒙特卡罗模拟方法成功对实验数据进行了数值拟合。关键词?相干光学;混浊物质;蒙特卡罗法;高斯光束中图分类号?O43?文献标识码?A?doi:10.3788/AOS201131.0626002Fast Simulation and Interference Properties of PhotonsTransmission in Turbid MaterialLin Lin1?Gao Yingjun2?Jin Chongxing2?Lin Huizhen2?Luo Yunhan21Department of Biological Engineering,Guangdong Medical College,Dongguan,Guangdong 523808,China2Department of Optoelectr onic Engineering,College of Science and Engineering,Jinan University,Guangzhou,Guangdong 510632,ChinaAbstract?In order to analyze low?coherent signal characteristics of turbid material,Monte Carlo(MC)numericalmodel is established.The Gaussian distribution of optical intensity and coherent gate of optical source are introducedinto numerical simulation.An accelerated method is introduced into simulation of coherent intensity of turbidmedium.Under condition of forword scattering,anisotropy factor g=0.9,the results of scaling method match withthose of MC simulation with mean relative error of 2.37%.When g decreases to 0.5,and mean relative errorincreases to 4.69%.In calculation of coherent signals of different depths,similar datain surface layer are gotten bythe two methods.In experiment,optical signals from 0 to 1 mm depth in the mixture liqid of IntralipidTMand Indiaink are measured by low?coherent system and fitted by accelerated Monte?Carlo method succesfully.Key words?coherence optics;turbid material;Monte Carlo method;Gaussian beamOCIS codes?260.0260;260.3160;030.0030;290.0290?收稿日期:2010?12?10;收到修改稿日期:2011?01?18基金项目:国家自然科学基金(61008057)资助课题。作者简介:林?林(1979?),男,博士研究生,主要从事生物医学光学和光学传感等方面的研究。E?mail:lynwindsent 导师简介:高应俊(1946?),男,教授,主要从事光纤传感和生物医学光学成像相关等方面的研究。E?mail:tyjgao 1?引?言随着光学技术的快速发展,激光及相关技术在生物医学领域的应用得到了广泛的关注和发展,例如共焦显微成像 1、光学相干层析成像 2、荧光成像 3、漫反射成像 4以及激光动力医疗 5等技术手段已经得到了临床应用。光学干涉测量技术可以精0626002?1光?学?学?报确定位目标、提供高分辨成像结果,因此得到了广泛的关注。在光的辐射传播过程中,生物组织对光子呈现出了吸收和散射并存的混浊特性,称为混浊物质。光子在混浊物质中受到的散射和吸收作用使得光子传输和干涉呈现出与透明介质不同的特征。在干涉测量系统中,由于光子在混浊物质中的传输受到不同尺寸粒子散射作用的影响呈现出随机性,随着探测深度的增加,多次散射作用将决定后向传输光子的空间分布和能量大小,从而决定了信号的强度与探测深度的依赖关系 6。对光子在混浊物质中的传播过程和相干性的研究一般采用基于辐射传输理论扩散近似法和蒙特卡罗(Monte?Carlo)数值模拟。前者只能在特定边界条件下得到解析解 7。相较而言,蒙特卡罗数值模拟方法没有引入近似条件,适合各种组织结构,因此在对组织光学的研究中一直得到了广泛的应用 8,其模拟结果也被用来验证各种理论模型的正确性。但是,蒙特卡罗法也存在耗时长 9,对低概率事件模拟的统计偏差较大等问题。特别是在光学干涉测量系统中,经过相干门和探测光路筛选后的反射光子数量过少会导致统计偏差增加,这就限制了数值模拟在光学干涉系统中的应用。本文以迈克耳孙光学低相干干涉系统作为研究对象,在高斯聚焦光束的入射条件下,对半无穷均匀混浊物质的光学干涉信号进行了蒙特卡罗数值模拟。根据光子在相干平面上的位置和时延确定其对干涉信号的贡献值。在蒙特卡罗数值模拟所得基线数据的基础上,对不同散射和吸收系数样品实现了快速数值模拟。实验中对不同浓度的脂肪悬乳注射剂(IntralipidTM)和印度墨水(India ink)混合液进行了干涉测量,并与数值分析结果进行了拟合对比。2?光子传输特性的理论分析图 1 干涉系统结构示意图Fig.1 Schematic of optical interference system低相干干涉测量系统(如图 1 所示)包括低相干光源、光纤迈克耳孙干涉结构、参考臂的光学扫描平台、样品臂的聚焦光路、光电探测器以及数据采集和处理系统。其中,样品臂的聚焦结构如图 2 所示,样品臂输出光经过物镜会聚在样品上,来自样品的后向散射光经过物镜返回到样品臂并且最终达到光电探测器与参考光叠加并产生干涉。图 2 干涉系统样品臂结构示意图Fig.2 Schematic of sample arm in the interferencesystem由干涉光信号形成的电流信号的均方值可以表示为 10?i2(z)?=2?2g(?)2?Re?S(p1,p2;z)?R(p1,p2)dp1dp2,(1)式中?=q?/h?是光电探测器的转换系数,q 是单位电荷的电量,?是光电探测器的量子效率,h 是普朗克常数,?是光子频率。p1,p2是干涉平面上的空间坐标,|g(?)|代表低相干光源时间相干函数的归一化模,?代表参考光与样品光的时差,当?=0 时,|g(?)|=1。被积函数为?R(p1,p2)=UR(p1)U*R(p2),(2)?S(p1,p2;z)=?US(p1;z)U*S(p2;z)?,(3)(2),(3)式分别代表了参考光和信号光在干涉平面上的互相关函数,?表示对时间取平均。在高散射介质中和傍轴近似的条件下,近红外光的散射以前向散射为主,散射方向角余弦统计平均值,即各向异性因子 g 接近于 1。一般的,参考光和信号光均为高斯光束,其传输方向上横截面光场 UR和 US分别为UR(p,t)=PR?20exp-p221?20+ikf?exp i?at+?(t),(4)US(p,t)=PS?20exp-p221?20+ikf?exp(i?bt),(5)0626002?2林?林等:?混浊物质中光子传输的快速模拟与干涉特性分析式中 PR和 PS分别是参考臂和样品臂的入射光功率,?0是高斯光束的束腰半径,k=2?/?,?a和?b分别代表参考光和样品光的角频率,?(t)为两干涉臂光场相位差,f 为样品臂聚焦透镜的焦距,p 代表径向位置,k 为传播常数,U*R和 U*S为光场函数的共轭。由(1)(5)式可见,从样品反射回来的信号光对干涉信号的贡献取决于两方面:1)信号光与参考光进行叠加的位置与干涉平面中心点的相对距离;2)信号光与参考光的延时差。当采用大量光子模拟干涉光强特性时,(1)式的积分运算就转化为对光子的求和。3?光子传输特性的蒙特卡罗模拟光子传输特性的蒙特卡罗模拟算法参考了多层介质中光子传输的蒙特卡罗光传输模型(MCML)标准 C 程序 11。MCML 程序并不能描述聚焦光束中光子的初始状态,也无法区分不同区域的散射光子以及光子的延迟特性。结合低相干干涉系统的光学结构和时间相干约束条件开发出用于光子传输特性模拟的 OCMC C 程序。数值模拟过程主要包括:1)光子初始状态的确定;2)光子在混浊物质中的传输追迹;3)来自目标层有效光子的标记;4)光子在干涉平面上与参考光的叠加。其中,只有 2)采用了MCML 中的相关子程序。当光子能量小于阈值能量时或者累积光程超出相干长度时,对该光子的模拟过程即结束。3.1?入射光子的初始状态在干涉系统中,样品光一般为聚焦状态的高斯光束,其在传输方向的任意横截面上光强分布为 12I(p)=?20?2(z)exp-p2?2(z),(6)式中?0为高斯光束在焦面上的束腰半径,?(z)为距束腰纵向 z 处的高斯光斑半径。采用了单页双曲面的几何结构模拟高斯光束的纵向传输,光子密度在与传输方向垂直的任意横截面上服从高斯分布,且光子沿双曲面的任一直母线运动,在样品表面的初始位置和方向可表示为x0=0.5?0-2ln r/r(?1+K?2),(7)y0=0.5?0-2ln r/r(?2+K?1),(8)?z0=cos?0=zx20+y20K2+z2,(9)式中?1和?2为在 0,1 区间内不相关的随机数。r=?21+?22,K=?z/?20,z 代表高斯光束在自由空间传输时束腰的纵向坐标,?z0为入射光子的初始方向角?0的余弦值。3.2?样品光与参考光的叠加由于参与蒙特卡罗模拟的单个光子的传输过程和能量状态不能用来表示光束的波动特性,因此根据数值孔径或入射角度来确定光子是否对干涉信号有贡献仍然存在较大争议。当入射光在干涉面上呈高斯分布时,来自样品的信号光对干涉信号的贡献与其在干涉平面上与中心的相对位置相关,干涉信号强度的均方值可表示为?I2(z)?coherent?Mm=0Ir(pm)wm,(10)式中 Ir(pm)为参考光在干涉面上 pm位置的强度密度,wm代表满足相干条件的第m 个光子的能量,?表示结果是统计平均值。OCMC 程序会统计每个后向散射光子的累积光程,并根据其延时情况和光源时间相干函数 g(?)确定其对干涉信号的贡献,则(10)式可以表示为?I2(z)?coherent=?Mm=0gm(?)Ir(pm)wm.(11)3.3?快速蒙特卡罗法在对相干光信号的数值模拟中,不同深度的干涉信号都需要一次独立的模拟过程。而光子传输过程与物质的光学参数密切相关,其出射权重与散射系数和吸收系数满足一定得比例关系 13。光子传输了距离 si后与物质发生相互作用的概率为pi=?texp(-?tsi),(12)如果光子与物质发生了 N 次相互作用,则总概率为P=?Ni=1?texp(-?tsi)=?Ntexp(-?tS),(13)式中?t=?s+?a为衰减系数,S=?si代表光子的累积传输程长。光子的最终权重为 w0(?s/?t)N,w0为初始权重。当散射系数由?s变为?*s,吸收系数由?a变为?*a时,由(13)式可得光子的出射权重W 与W*的关系为 14,15W*=W?*s?sNexp-(?*t-?t)S.(14)利用(14)式所确定的运算关系和已经获得蒙特卡罗模拟数据可以快速推算具有不同散射和吸收系数的物质相同深度的干涉信号强度。该方法只需要运行一次完整的蒙特卡罗模拟,并保存参与干涉的光子信息:出射权重、被散射次数和累积传输距离,所存储的数据称为基线数据。当样品散射和吸收系数发生变化时,通过(14)式和基线数据确定出射光子0626002?3光?学?学?报的新权重和对干涉信号的贡献。由于(14)式是在各向异性因子 g 值不变的情况下得出的,因此,在对基线数据进行比例缩放的计算过程中要求 g 为常数。4?数值模拟结果和实验验证以均匀半无穷混浊介质为对象进行了数值模拟。样品光路参数为物镜焦距 f=16 mm,焦平面上光斑半径?0=30?m,样品折射率 n=1.35,入射光聚焦在样品表面。为了验证快速蒙特卡罗法的有效范围,首先选择干涉光子来自样品内 0.2 mm 深处。在g=0.5 和 0.9 两种散射条件下,选取同样的15 组散射和吸收系数组合进行蒙特卡罗模拟,即每组参数为(?ai,?si)=(i cm-1,10i cm-1),i=1,?,15。为了减少统计偏差,每组参数进行 5 次模拟后取平均值,结果如图 3(a)和(b)中曲线所示。之后利用表 1 中的四组参数进行蒙特卡罗模拟,以获得四组基线数据,表中 M 表示单次蒙特卡罗运算的光子数。利用这四组基线数据再次对以上 15 组参数样品在 0.2 mm 深处干涉信号进行快速蒙特卡罗模拟运算,所得结果如图 3(a)和(b)中符号所示。可以看到两种方法所得结果基本相吻合。图 3(c)和(d)给出了快速运算结果与蒙特卡罗模拟结果之间的相对误差统计。可以看到当 g=0.9 时最大和平均相对误差分别是 10.82%和 2.37%。当 g=0.5时,模拟结果的最大和平均相对误差分别是27.46%和4.69%,相对误差有所放大,特别是目标参数为高吸收和高散射系数时相对误差出现明显放大。造成误差的主要原因是基线数据参数与目标参数距离过大,例如使用第一组基线数据(?a,?s)=(1 cm-1,30 cm-1)计算(?ai,?si)=(15 cm-1,150 cm-1)条件下干涉信号时,所得误差在两种 g 值条件下都是最大的。表 1 蒙特卡罗模拟中的光学参数Table 1Optical parameters in Monte Carlo simulationGroup?a/cm-1?s/cm-1gM11300.54?1080.92?10822400.54?1080.92?10833500.54?1080.92?10844800.54?1080.92?108图 3 蒙特卡罗与快速模拟运算结果比较(a),(b)和相对误差统计(c),(d)Fig.3 Comparism between the results of Monte?Carlo and accelerated simulation(a),(b)andrelative error statistics(c),(d)0626002?4林?林等:?混浊物质中光子传输的快速模拟与干涉特性分析?图 4 给出了四组参数样品的不同深度(zr=0.1 1.0 mm)干涉信号的快速运算结果与蒙特卡罗模拟结果的比较,其中各向异性因子为 g=0.9,0.5。基线 数据所 用参 数为?a=3 cm-1,?s=50 cm-1。可以看到两种数值模拟结果在浅表深度基本吻合,随深度增加,相对误差会增加。这主要是蒙特卡罗模拟在深度增加时参与干涉的光子数量减少导致统计偏差变大的结果,这在?a=6 cm-1,?s=80 cm-1对应曲线中尤为明显。在本次计算中,利用传统蒙特卡罗方法完成 10 个不同深度干涉信号模拟平均需要 10.4 h,而利用基线数据进行快速运算,平均只需要 4.98 s。运算所使用的计算机配置为:配备 Intel T8100 2.1 GHz CoreTM2 Duo CPU的 IBM X61 型计算机。实验中,低相干光源采用中心波长 1550 nm 自激发辐射光源(ASE)光源,光谱半峰全宽为 80 nm。样品臂光束通过 f=16 mm 的显微物镜聚焦在样品表面,焦面光斑半径为?0=30?m。参考臂采用水平扫描结构,通过步进电机驱动反射镜匀速水平运动以实现纵向扫描(A?Scan),扫描速度 2 mm/s。干涉信号经过差分放大器和带通滤波器后通过 NI(National Instrument Cot.)公司数据采集卡读入计算机,之后进行了希尔伯特变换解调并提取了包络以及归一化处理,最终得到了不同深度的干涉信号。样品为不同浓度脂肪悬乳注射剂和印度墨水的混合液。利用快速蒙特卡罗方法对实验数据进行了拟合,并通过最小二乘法获得如图 5 所示结果。拟合使用的基线数据来自?a=3 cm-1,?s=50 cm-1的蒙特卡罗模拟数据。根据文献 16,在 1550 nm 波段脂肪悬乳注射剂的 g?0.36,并且可以认为加入印度墨水或改变脂肪悬乳注射剂浓度对 g 值影响很小。图 5(a)中样品为体积分数 10%脂肪悬乳注射剂与体积分数 0.5%印度墨水混合液,其拟合曲线对应的参数为?a=8.12 cm-1,?s=42.36 cm-1;图 5(b)中样品为 15%印度墨水与 0.3%印度墨水混合液,其拟合曲线对应的参数为?a=7.55 cm-1,?s=47.25 cm-1。拟合结果说明基于比例缩放的快速蒙特卡罗方法可以对散射和吸收系数分别进行提取或重构计算。图 4 不同深度干涉信号的蒙特卡罗模拟结果(符号)与快速模拟计算结果(曲线)对比Fig.4 Comparison of interference signals between Monte Carlo(symbol)and scaling method(line)in different depth图 5 实验干涉信号(虚线)与快速蒙特卡罗拟合结果(实线)Fig.5 Data from interference experiment(dotted line)and fit by fast Monte?Carlo fitting(solid line)5?结?论对混浊物质中的低相干干涉信号进行了数值模拟,实现了干涉信号的快速数值计算。在各向异性因子不变的条件下,可以在蒙特卡罗模拟数据的基0626002?5光?学?学?报础上快速地获得其他光学参数条件下的干涉信号强度。通过与传统蒙特卡罗模拟结果进行对比,运算速度较之有大幅提高。在实验中对脂肪悬乳注射剂和印度墨水的混合液的干涉测量数据成功进行了拟合,并获得了样品的光学参数,验证了该方法在数据拟合和参数重构中的可行性。参考文献1 S.G.Daniel,A.Sanjee,R.Milind.Line?scanning reflectanceconfocal microscopy of human skin:comparison of full?pupil anddivided?pupil configurations J.Op t.Lett.,2009,34(20):3235 32372 Li Peng,Huang Run,Gao Wanrong.Experiment research onoptical coherence tomography of human skin J .Chinese J.Lasers,2009,36(10):2498 2502?李?鹏,黄?润,高万荣.光学相干层析术在人体皮肤成像方面的实验研究 J.中国激光,2009,36(10):2498 25023 Li Zuangfang,Huang Zufang,Chen Rong et al.T wo photonfluorescence imaging of thyroid tissue J.Chinese J.Lasers,2009,36(3):765 768?李钻芳,黄主芳,陈?荣 等.甲状腺组织的双光子荧光成像 J.中国激光,2009,36(3):765 7684 Li Lanquan,Wei Huajiang,He Bohua et al.Gastric cancerdetection using diffuse reflectance spectral ratio R540/R575ofoxyg enated hemoglobin band J.Acta Op tica Sinica,2009,29(10):2860 2865?李兰权,魏华江,何博华 等.采用氧合血红蛋白吸收带的漫反射光谱比率 R540/R575 检测人胃癌 J.光学学报,2009,29(10):2860 28655 F.L.Primo,M.V.L.B.Bentley,A.C.T edesco.Photophysical studies and in vitro skin permeation/retention offoscan/nanoemulsion(NE)applicable to photodynamic therapyskin cancer treatment J.J.Nanoscience N anotechnology,2008,8(1):340 3476 J.M.Schmitt,A.Kn?ttel.M odelofopticalcoherencetomography of heterogeneous tissue J.J.Op t.Soc.Am.A,1997,14(6):1231 12427 Y.Deng,M.Ignor.Latest progress of Monte Carlo simulationof light transmission in tissues J.Acta Physica Sinica,2010,59(2):1398 1401?邓?勇,M.Ignor.M onte Carlo 方法模拟光在生物组织中传播的新进展 J.物理学报,2010,59(2):1398 14018 A.T ycho,T.M.J?rgensen,H.T.Yura et al.Derivation of aMonte Carlo method for modeling heterodyne detection in opticalcoherence tomography systems J.Ap p l.Opt.,2002,41(31):6676 66919 Q.Fang,D.A.Boas.Monte Carlo simulation of photonmigration in 3D turbid media accelerated by graphics processingunits J.Op t.Exp ress,2009,17(22):20178 2019010 J.E.Bender,K.Vishwanath,L.K.Moore et al.A robustMonte Carlo model for the extraction of biological absorption andscattering in vivo J.IEEE T rans.Biomed.Engng.,2009,56(4):960 96811 L.H.Wang,S.J.Jacques,L.Q.Zheng.MCML?Monte Carlomodeling of photon transport in multi?layered tissue J.Comp ut.Methods Programs Biomed.,1995,47(2):131 14612 L.Lin,M.Zhang,H.L.Zhang.Numerical simulation offocused optical beamin tissue samples J.Chin.J.Med.Imaging T echnol.,2010,26(12):2375 2378?林?林,张?梅,张怀岺.聚焦光束在生物组织中传输的数值模拟 J.中国医学影像技术,2010,26(12):2375 237813 Q.Liu,N.Ramanujam.Scaling method for fast M onte Carlosimulation of diffuse reflectance spectra from multilayered turbidmedia J.J.Op t.Soc.Am.A,2007,24(4):1011 102514 C.K.Hayakawa,J.Spanier,F.Bevilacqua et al.PerturbationMonte Carlo methods to solve inverse photon migration problemsin heterogeneous tissues J.Op t.Lett.,2001,26(17):1335 133715 A.Sassaroli,C.Blumetti,F.Martelli et al.Monte Carloprocedure forinvestigating light propagation and imaging ofhighly scatteringmedia J .Ap pl.Op t.,1998,37(31):7392 740016 C.Chen,J.Q.Lu,H.Ding et al.A primary method fordetermination ofopticalparameters ofturbid samplesandapplication to intralipid between 550 and 1630 nm J.Opt.Exp ress,2006,14(6):7421 74350626002?6