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

ACTA OPTICA SINICA

June,2011

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

林 林1 高应俊2 金重星2 林惠珍2 罗云瀚

2

1

广东医学院生物医学工程教研室,广东东莞523808

2

暨南大学理工学院光电工程系,广东广州510632

摘要 建立了混浊物质中低相干干涉系统的蒙特卡罗模型,结合相干面光强的高斯分布和低相干光源的相干条件,对干涉信号进行了数值分析。提出在蒙特卡罗计算的基线数据基础上快速模拟其他参数条件下干涉信号的方法。快速模拟结果与蒙特卡罗计算结果进行了比较,结果显示在满足前向散射条件下,两种方法所得模拟信号强度基本吻合。当各向异性因子g =0.9时,平均相对误差为2.37%,在g =0.5时,平均相对误差为4.69%。在深度相关干涉信号的模拟中,两种方法对浅表层干涉信号的模拟取得了一致的结果。实验中使用低相干干涉系统对不同浓度的脂肪悬乳注射剂与印度墨水混合液中深度为0~1mm 的漫反射光信号进行了干涉测量,并利用基线数据和快速蒙特卡罗模拟方法成功对实验数据进行了数值拟合。关键词 相干光学;混浊物质;蒙特卡罗法;高斯光束

中图分类号 O43 文献标识码 A doi :10.3788/AOS 201131.0626002

Fast Simulation and Inte rfe re nce Prope rties of Photons

Transmission in Turbid Mate rial

Lin Lin 1 Gao Yingjun 2 Jin Chongxing 2 Lin Huizhen 2 Luo Yunhan 2

1

Depa r t m en t of Biological Engineer ing ,Gu an gdon g Medica l College ,Donggua n ,Gua ngdong 523808,Chin a

2

Depa r tm en t of Optoelectr onic En gineer in g ,College of Scien ce an d En gin eer in g ,Jin an Un iver sit y ,

Gua ngzhou ,Gu an gdong 510632,China

Abstract I n order to analyze low coherent signal characteristics of turbid material,Monte Carlo(MC)numeric al m odel is established.The Gaussian distribution of optical intensity and c oherent gate of optical source are introduced into numerical simulation.An acc elerated method is introduced into simulation of coherent intensity of turbid m edium.Under condition of forword sc attering,anisotropy fac tor g =0.9,the results of scaling method match with those of MC simulation with mean relative error of 2.37%.When g decrea ses to 0.5,and mean relative error increases to 4.69%.I n calculation of coherent signals of different depths,similar data in surface layer are gotten by the t wo methods.I n experiment,optical signals from 0to 1mm depth in the m ixt ure liqid of Intralipid TM and I ndia ink are measured by low coherent system and fitted by accelerated Monte Carlo method succesfully.Key wo rds coherence optic s;turbid m aterial;Monte Carlo method;Gaussian beam OCIS co des 260.0260;260.3160;030.0030;290.0290

收稿日期:2010 12 10;收到修改稿日期:2011 01 18

基金项目:国家自然科学基金(61008057)资助课题。

作者简介:林 林(1979 ),男,博士研究生,主要从事生物医学光学和光学传感等方面的研究。E mail:ly nw indsent@yaho o.com.cn

导师简介:高应俊(1946 ),男,教授,主要从事光纤传感和生物医学光学成像相关等方面的研究。E mail:tyjgao @jnu.edu.cn

1 引 言

随着光学技术的快速发展,激光及相关技术在生物医学领域的应用得到了广泛的关注和发展,例

如共焦显微成像[1]、光学相干层析成像[2]

、荧光成

像[3]、漫反射成像[4]以及激光动力医疗[5]等技术手段已经得到了临床应用。光学干涉测量技术可以精

光 学 学 报

确定位目标、提供高分辨成像结果,因此得到了广泛的关注。在光的辐射传播过程中,生物组织对光子呈现出了吸收和散射并存的混浊特性,称为混浊物质。光子在混浊物质中受到的散射和吸收作用使得光子传输和干涉呈现出与透明介质不同的特征。在干涉测量系统中,由于光子在混浊物质中的传输受到不同尺寸粒子散射作用的影响呈现出随机性,随着探测深度的增加,多次散射作用将决定后向传输光子的空间分布和能量大小,从而决定了信号的强度与探测深度的依赖关系

[6]

。对光子在混浊物质中

的传播过程和相干性的研究一般采用基于辐射传输理论扩散近似法和蒙特卡罗(M onte Carlo )数值模拟。前者只能在特定边界条件下得到解析解[7]

。相较而言,蒙特卡罗数值模拟方法没有引入近似条件,适合各种组织结构,因此在对组织光学的研究中一直得到了广泛的应用[8],其模拟结果也被用来验证各种理论模型的正确性。但是,蒙特卡罗法也存在耗时长[9],对低概率事件模拟的统计偏差较大等问题。特别是在光学干涉测量系统中,经过相干门和探测光路筛选后的反射光子数量过少会导致统计偏差增加,这就了数值模拟在光学干涉系统中的应用。

本文以迈克耳孙光学低相干干涉系统作为研究对象,在高斯聚焦光束的入射条件下,对半无穷均匀混浊物质的光学干涉信号进行了蒙特卡罗数值模拟。根据光子在相干平面上的位置和时延确定其对干涉信号的贡献值。在蒙特卡罗数值模拟所得基线数据的基础上,对不同散射和吸收系数样品实现了快速数值模拟。实验中对不同浓度的脂肪悬乳注射剂(Intralipid TM )和印度墨水(India ink)混合液进行了干涉测量,并与数值分析结果进行了拟合对比。

2 光子传输特性的理论分析

图1干涉系统结构示意图

F ig.1Schematic o f optical interference system

低相干干涉测量系统(如图1所示)包括低相干

光源、光纤迈克耳孙干涉结构、参考臂的光学扫描平台、样品臂的聚焦光路、光电探测器以及数据采集和处理系统。其中,样品臂的聚焦结构如图2所示,样品臂输出光经过物镜会聚在样品上,来自样品的后向散射光经过物镜返回到样品臂并且最终达到光电

探测器与参考光叠加并产生干涉。

图2干涉系统样品臂结构示意图

F ig.2Schematic of sample ar m in the inter ference

sy st em

由干涉光信号形成的电流信号的均方值可以表示为[10]

i 2

(z ) =2 2

g( )2

Re

S

(p 1

,p 2

;z ) R

(p 1

,p 2

)d p 1

d p

2

,(1)

式中 =q /h 是光电探测器的转换系数,q 是单位电荷的电量, 是光电探测器的量子效率,h 是普朗克常数, 是光子频率。p 1,p 2是干涉平面上的空间坐标,|g( )|代表低相干光源时间相干函数的归一化模, 代表参考光与样品光的时差,当 =0时,|g( )|=1。被积函数为

R (p 1,p 2)=U R (p 1)U *R (p 2),(2) S (p 1,p 2;z )= U S (p 1;z )U *S (p 2;z ) ,

(3)

(2),(3)式分别代表了参考光和信号光在干涉平面上

的互相关函数, 表示对时间取平均。在高散射介质中和傍轴近似的条件下,近红外光的散射以前向散射为主,散射方向角余弦统计平均值,即各向异性因子g 接近于1。一般的,参考光和信号光均为高斯光束,其传输方向上横截面光场U R 和U S 分别为

U R (p ,t)=

P R 20

ex p -p 2

21 20+i k f exp i a t+ (t),

(4)U S (p ,t)=

P S 20ex p -p 221 20+i k f exp (i b t),

(5)

式中P R和P S分别是参考臂和样品臂的入射光功

率, 0是高斯光束的束腰半径,k=2 / , a和 b分

别代表参考光和样品光的角频率, (t)为两干涉臂

光场相位差,f为样品臂聚焦透镜的焦距,p代表径

向位置,k为传播常数,U*R和U*S为光场函数的共

轭。由(1)~(5)式可见,从样品反射回来的信号光

对干涉信号的贡献取决于两方面:1)信号光与参考

光进行叠加的位置与干涉平面中心点的相对距离;

2)信号光与参考光的延时差。当采用大量光子模拟

干涉光强特性时,(1)式的积分运算就转化为对光子

的求和。

3 光子传输特性的蒙特卡罗模拟

光子传输特性的蒙特卡罗模拟算法参考了多层

介质中光子传输的蒙特卡罗光传输模型(MCM L)

标准C程序[11]。M CML程序并不能描述聚焦光束

中光子的初始状态,也无法区分不同区域的散射光

子以及光子的延迟特性。结合低相干干涉系统的光

学结构和时间相干约束条件开发出用于光子传输特

性模拟的OCMC C程序。数值模拟过程主要包括:

1)光子初始状态的确定;2)光子在混浊物质中的传

输追迹;3)来自目标层有效光子的标记;4)光子在干

涉平面上与参考光的叠加。其中,只有2)采用了

MCM L中的相关子程序。当光子能量小于阈值能

量时或者累积光程超出相干长度时,对该光子的模

拟过程即结束。

3.1 入射光子的初始状态

在干涉系统中,样品光一般为聚焦状态的高斯

光束,其在传输方向的任意横截面上光强分布为[12]

I(p)=

20

2(z)ex p-

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=co s 0=z

x20+y20K2+z2

,(9)式中 1和 2为在[0,1]区间内不相关的随机数。r= 21+ 22,K= z/ 20,z代表高斯光束在自由空间传输时束腰的纵向坐标, z0为入射光子的初始方向角 0的余弦值。

3.2 样品光与参考光的叠加

由于参与蒙特卡罗模拟的单个光子的传输过程和能量状态不能用来表示光束的波动特性,因此根据数值孔径或入射角度来确定光子是否对干涉信号有贡献仍然存在较大争议。当入射光在干涉面上呈高斯分布时,来自样品的信号光对干涉信号的贡献与其在干涉平面上与中心的相对位置相关,干涉信号强度的均方值可表示为

I2(z) coher ent M

m=0

I r(p m)w m,(10)式中I r(p m)为参考光在干涉面上p m位置的强度密度,w m代表满足相干条件的第m个光子的能量, 表示结果是统计平均值。OCM C程序会统计每个后向散射光子的累积光程,并根据其延时情况和光源时间相干函数g( )确定其对干涉信号的贡献,则(10)式可以表示为

I2(z) coherent= M

m=0

g m( )I r(p m)w m.(11)

3.3 快速蒙特卡罗法

在对相干光信号的数值模拟中,不同深度的干涉信号都需要一次的模拟过程。而光子传输过程与物质的光学参数密切相关,其出射权重与散射系数和吸收系数满足一定得比例关系[13]。光子传输了距离s i后与物质发生相互作用的概率为

p i= t ex p(- t s i),(12)如果光子与物质发生了N次相互作用,则总概率为P= N i=1 t exp(- t s i)= N t ex p(- t S),(13)式中 t= s+ a为衰减系数,S= s i代表光子的累积传输程长。光子的最终权重为w0( s/ t)N,w0为初始权重。当散射系数由 s变为 *s,吸收系数由 a变为 *a时,由(13)式可得光子的出射权重W与W*的关系为[14,15]

W*=W

*s

s

N

exp-( *t- t)S.(14)利用(14)式所确定的运算关系和已经获得蒙特卡罗模拟数据可以快速推算具有不同散射和吸收系数的物质相同深度的干涉信号强度。该方法只需要运行一次完整的蒙特卡罗模拟,并保存参与干涉的光子信息:出射权重、被散射次数和累积传输距离,所存储的数据称为基线数据。当样品散射和吸收系数发生变化时,通过(14)式和基线数据确定出射光子

光 学 学 报

的新权重和对干涉信号的贡献。由于(14)式是在各向异性因子g 值不变的情况下得出的,因此,在对基线数据进行比例缩放的计算过程中要求g 为常数。

4 数值模拟结果和实验验证

以均匀半无穷混浊介质为对象进行了数值模拟。样品光路参数为物镜焦距f =16mm ,焦平面

上光斑半径 0=30 m,样品折射率n = 1.35,入射光聚焦在样品表面。为了验证快速蒙特卡罗法的有效范围,首先选择干涉光子来自样品内0.2mm 深处。在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.2mm 深处干涉信号进行快速蒙特卡罗模拟运算,所得结果如图3(a)和(b)中符号所示。可以看到两种方法所得结果基本相吻合。图3(c)和

(d)给出了快速运算结果与蒙特卡罗模拟结果之间的相对误差统计。可以看到当g =0.9时最大和平均相对误差分别是10.82%和2.37%。当g =0.5时,模拟结果的最大和平均相对误差分别是27.46%和4.69%,相对误差有所放大,特别是目标参数为高吸收和高散射系数时相对误差出现明显放大。造成误差的主要原因是基线数据参数与目标参数距离过

大,例如使用第一组基线数据( a , s )=(1cm -1

,30cm -1)计算( ai , si )=(15cm -1,150cm -1)条件

下干涉信号时,所得误差在两种g 值条件下都是最大的。

表1蒙特卡罗模拟中的光学参数

T able 1Optical par amet ers in M onte Carlo simulation G ro up a /cm -1

s /cm -1

g M 11300.5

4 1080.92 10822400.54 1080.92 10833500.54 1080.92 1084

4

80

0.54 1080.9

2 10

8

图3蒙特卡罗与快速模拟运算结果比较(a),(b)和相对误差统计(c),(d)

F ig.3Compar ism betw een the r esults o f M onte Carlo and accelerat ed simulation (a),(b)and

relat ive er ror statistics (c),(d)

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

图4给出了四组参数样品的不同深度(z r =0.1~1.0mm )干涉信号的快速运算结果与蒙特卡罗模拟结果的比较,其中各向异性因子为g =0.9,

0.5。基线数据所用参数为 a =3cm -1

, s =50cm -1。可以看到两种数值模拟结果在浅表深度基本吻合,随深度增加,相对误差会增加。这主要是蒙特卡罗模拟在深度增加时参与干涉的光子数量减

少导致统计偏差变大的结果,这在 a =6cm -1

, s =80cm -1对应曲线中尤为明显。在本次计算中,利用传统蒙特卡罗方法完成10个不同深度干涉信号模拟平均需要10.4h,而利用基线数据进行快速运算,平均只需要4.98s 。运算所使用的计算机配置为:配备Intel T81002.1GH z Co re TM 2Duo CPU 的IBM X61型计算机。

实验中,低相干光源采用中心波长1550nm 自激发辐射光源(ASE)光源,光谱半峰全宽为80nm 。样品臂光束通过f =16m m 的显微物镜聚焦在样品表面,焦面光斑半径为 0=30 m 。参考臂采用水平扫描结构,通过步进电机驱动反射镜匀速水平运动以实现纵向扫描(A Scan),扫描速度2mm/s 。

干涉信号经过差分放大器和带通滤波器后通过NI (National Instr ument Co t.)公司数据采集卡读入计算机,之后进行了希尔伯特变换解调并提取了包络以及归一化处理,最终得到了不同深度的干涉信号。样品为不同浓度脂肪悬乳注射剂和印度墨水的混合液。利用快速蒙特卡罗方法对实验数据进行了拟合,并通过最小二乘法获得如图5所示结果。拟合

使用的基线数据来自 a =3cm -1, s =50cm -1

的蒙

特卡罗模拟数据。根据文献[16],在1550nm 波段脂肪悬乳注射剂的g 0.36,并且可以认为加入印度墨水或改变脂肪悬乳注射剂浓度对g 值影响很小。图5(a)中样品为体积分数10%脂肪悬乳注射剂与体积分数0.5%印度墨水混合液,其拟合曲线

对应的参数为 a =8.12cm -1

, s =42.36cm -1

;

图5(b)中样品为15%印度墨水与0.3%印度墨水混合液,其拟合曲线对应的参数为 a =7.55cm -1

, s =47.25cm -1。拟合结果说明基于比例缩放的快

速蒙特卡罗方法可以对散射和吸收系数分别进行提

取或重构计算。

图4不同深度干涉信号的蒙特卡罗模拟结果(符号)与快速模拟计算结果(曲线)对比

F ig.4Compariso n of interfer ence sig nals betw een M o nt e Car lo (symbol)and scaling method (line)in different

depth

图5实验干涉信号(虚线)与快速蒙特卡罗拟合结果(实线)

F ig.5Dat a f rom inter ference ex periment (dott ed line)and fit by fast M onte Ca rlo fitting (so lid line)

5 结 论

对混浊物质中的低相干干涉信号进行了数值模

拟,实现了干涉信号的快速数值计算。在各向异性因子不变的条件下,可以在蒙特卡罗模拟数据的基

础上快速地获得其他光学参数条件下的干涉信号强度。通过与传统蒙特卡罗模拟结果进行对比,运算速度较之有大幅提高。在实验中对脂肪悬乳注射剂和印度墨水的混合液的干涉测量数据成功进行了拟合,并获得了样品的光学参数,验证了该方法在数据拟合和参数重构中的可行性。

参考文献

1S.G.Daniel, A.S an jee,R.M ilin d.Line s can ning reflectance confocal microscopy of hu man skin:comparison of fu ll pupil and divided pu pil con figurations[J].Op t.L ett.,2009,34(20): 3235~3237

2Li Pen g,H uang Ru n,Gao Wanrong.Experiment research on optical coheren ce tomography of h uman skin[J].Chinese J.

L asers,2009,36(10):2498~2502

李 鹏,黄 润,高万荣.光学相干层析术在人体皮肤成像方面的实验研究[J].中国激光,2009,36(10):2498~2502

3Li Zuangfan g,H uang Zufang,Chen Rong et al..T w o photon flu or escence im aging of thyroid tissu e[J].Chinese J.Laser s, 2009,36(3):765~768

李钻芳,黄主芳,陈 荣等.甲状腺组织的双光子荧光成像[J].

中国激光,2009,36(3):765~768

4Li Lanquan,W ei Hu ajiang,H e Bohua et al..Gastric cancer detection using diffu se reflectance spectral ratio R540/R575of oxyg enated hemoglobin band[J].Ac ta Op tic a S inica,2009, 29(10):2860~2865

李兰权,魏华江,何博华等.采用氧合血红蛋白吸收带的漫反射光谱比率R540/R575检测人胃癌[J].光学学报,2009,29(10): 2860~2865

5F.L.Primo,M.V.L. B.Bentley, A. C.T edesco.

Ph otophys ical studies an d in vitr o s kin per meation/retention of foscan/nanoemulsion(NE)applicable to photodyn amic therapy s kin cancer treatment[J].J.Nanosc ienc e N anotechn olog y, 2008,8(1):340~347

6J.M.S chmitt, A.Kn ttel.M od el of optical coherence tomography of heterogen eous tissue[J].J.Op t.Soc.Am.A, 1997,14(6):1231~1242

7Y.Deng,M.Ignor.Latest progres s of M onte Carlo sim ulation of ligh t tran smis sion in tis sues[J].Acta Ph ysica S inica,2010, 59(2):1398~1401

邓 勇,M.Ignor.M on te Carlo方法模拟光在生物组织中传播的新进展[J].物理学报,2010,59(2):1398~1401

8A.T ycho,T.M.J rgensen,H.T.Yu ra e t al..Derivation of a M onte Carlo method for modeling heterodyne detection in optical coheren ce tomograp hy s ystem s[J].Ap p l.Opt.,2002,41(31): 6676~6691

9Q.Fan g, D. A.Boas.M onte Carlo simu lation of photon migration in3D turbid m edia accelerated by graphics processin g units[J].O p t.Exp re ss,2009,17(22):20178~20190

10J.E.Bender,K.Vishw anath,L.K.M oore et al..A robu st M onte Carlo model for the extraction of biological abs orption and scatterin g in vivo[J].IE EE T rans.Biomed.Eng ng.,2009, 56(4):960~968

11L.H.Wang,S.J.Jacques,L.Q.Zheng.M CM L M onte Carlo modeling of photon tran sport in m ulti layered tis sue[J].Comp ut.

Me thod s P rogr ams B iomed.,1995,47(2):131~146

12L.Lin,M.Zh ang,H.L.Zhan g.Num erical sim ulation of focused optical beam in tissue samples[J].Chin.J.Med.

Imag ing T ec hnol.,2010,26(12):2375~2378

林 林,张 梅,张怀岺.聚焦光束在生物组织中传输的数值模拟[J].中国医学影像技术,2010,26(12):2375~2378

13Q.Liu,N.Ramanujam.Scaling meth od for fast M onte Carlo simu lation of diffu se reflectance s pectra from multilayered turb id media[J].J.O p t.Soc.Am.A,2007,24(4):1011~1025

14C.K.Hayakaw a,J.Spanier, F.Bevilacqua et al..Perturbation M onte Carlo methods to solve inverse photon migration problems in heterogeneous tissues[J].Op t.Lett.,2001,26(17): 1335~1337

15A.Sass aroli, C.Blumetti, F.M artelli e t al..M onte Carlo procedure for inves tigating light propagation and imaging of highly scattering media[J].Ap pl.Op t.,1998,37(31): 7392~7400

16C.Chen,J.Q.Lu,H.Ding e t al..A prim ary m ethod for determination of optical param eter s of turb id samples an d application to intralipid b etw een550an d1630n m[J].Opt.

Exp re ss,2006,14(6):7421~7435