（光学工程专业论文）光学相干层析成像系统模拟研究.pdf

摘 要 光学相干层析成像 ( o p t i c a l c o h e r e n c e t o m o g r a p h y , o c t ) 是近年来快速发展 的一种生物组织高分辨率实时成像技术。 由于生物组织的复杂性， 一直没有完美 的o c t理论模型。本文希望通过对 o c t理论模型的探索，选择一种合适的模 型能解释光在生物组织中传播的各种现象。 本文共分6 章，各章的主要内容如下: 第1 章 绪论。综合的介绍和评述光学相干层析成像的理论和实践发展的近 况，提出本文的研究目的。 第2 章 o c t原理。 介绍了o c t的提出、 发展、以及o c t的传统物理模型。 用经典的光学模型分析 o c t . 第3 章电磁波在高散射介质中的传播。电磁波理论是成熟的经典理论，它 在现实的 各个方面都有涉及， 本章主要介绍电 磁理论中电 磁波在高散射介质中 传 播的理论和方法。 第4 章 蒙特卡罗模拟方法。本章介绍传统蒙特卡罗方法模拟光在高散射介 质中传播的模型和计算机模拟的经典的实现， 并在此基础上提出了基于mie 散射 的蒙特卡洛模拟。 第5 章 数值模拟。运用第四章发展的理论，实际编写程序模拟简单的 o c t 探测例子。同时比 较了模拟结果和现有实验结果。 第6 章 总结和展望。最终总结了本文的工作内容和得到的结果，并在此基 础上 提出了以后进一步研究的方向。 本文利用结合米氏理论( m i e t h e o r y ) 的蒙特卡洛算法研究光在生物组织( 高散 射介质)内的行为，并以此为基础建立了一种 o c t理论模型。然后运用该理论 编写程序模拟偏振光在单层高散射介质中的传播并与前人的实验结果比较， 再将 该程序运用于重建多层介质的图像，分析了空间滤波和偏振光对图像质量的影 浙江大学工程硕十学位论文 响，都获得了较理想的结果。 关键词:光学相干层析成像， 蒙特卡洛， 米氏 理论， 高散射介质 abs tract o p t i c a l c o h e r e n c e t o m o g r a p 勿 ( o c t ) i s a f as t d e v e l o p i n g t e c h n o l o g y f o r m e d i c a l d i a g n o s i s i n v i v o . b u t t i l l n o w , t h e r e s t i l l i s n t e x i s t w e l l - a c c e p t e d p 勿s i c a l m o d e l s o f o c t . t h i s p a p e r u s e s mo n t e c a r l o s i m u l a t i o n a n d mi e t h e o r y t o b u i l d u p t h e o c t m o d e l . o u r r e s e a r c h i s f o c u s o n f in d i n g a n e w m o d e l w h i c h c a n e x p l a i n t h e b e h a v i o r s o f l i g h t p r o p a g a t i o n i n h i g h s c a t t e r i n g m e d iu m s u c h a s t i s s u e . t h e p r e s e n t t h e s i s c o n s i s t s o f s i x c h a p t e r s . t h e c o n t e n t o f e a c h c h a p t e r i s l i s t e d b e l o w: c h a p t e r 1 : t h e h i s t o r y a n d t h e l a t e s t t r e n d o f t h e o p t i c a l c o h e r e n c e i m a g e d e v e l o p m e n t , a n d t h e a i m o f t h i s r e s e a r c h p r o j e c t . c h a p t e r 2 : b a s i c t h e o r y a n d p r in c i p a l o f o p t i c a l c o h e r e n c e t o m o g r a p 冲. c h a p t e r 3 : d i s c u s s s o m e w e l l - d e v e l o p e d e l e c t r o m a g n e t i c w a v e t h e o r i e s w h i c h a r e u s e d i n d e s c r i b i n g l i g h t p r o p a g a t i o n i n h i g h s c a t t e r i n g m e d i u m . l i g h t s c a tt e r i n g t h e o r y w h i c h i s w i d e l y a c c e p t e d t o s o lv e t h e p r o b l e m a b o u t p h o t o n s s c a t t e r i n g i n t i s s u e . w e a l s o n e e d t o r e v i e w t h e c l a s s i c e l e c t r o m a g n e t ic w a v e t h e o r y : mi e t h e o r y , wh i c h wi l l b e u s e d i n o u r s o l u t i o n . c h a p t e r 4 : mo n t e c a r l o s im u l a t i o n . a n i n t r o d u c t i o n a b o u t mo n t e c a r l o s i m u l a t i o n . a n d d i s c u s s a n d e x p l a i n h o w m o n t e c a r l o m o d e l i n g o f l i g h t t r a n s p o r t i n m u l t i - l a y e r e d t i s s u e s ( mc ml ) w o r k s . t h e n w e e x t e n d e d mc ml u s i n g mi e t h e o r y t o s i m u l a t e c o h e r e n c e l i g h t t r a n s p o r t i n m u l t i - l a y e r e d t i s s u e s . c h a p t e r 5 : n u m e r i c a l s i m u l a t i o n . i n t h i s c h a p t e r w e c o d e d o u r mo n t e c a r lo s imu l a t i o n b a s e d o n mi e t h e o r y i n c + + , a n d t h e n w e u s e d t h e p r o g r a m t o s i m u la t i n g s o me c a s e s . t h e s i mu l a t i o n r e s u l t s we r e d i s c u s s e d c h a p t e r 6 : c o n c l u s i o n a n d f o r e s e e . wh a t w e g e t a t l a s t a n d w h a t s h o u l d w e d o i n t h e f u t u r e . t h i s t h e s i s d e v e l o p e d a n e w m o d e l t o s i m u l a t e d c o h e r e n c e l i g h t t r a n s p o r t i n h i g h s c a t t e r i n g m e d i u m . o u r s i m u l a t i o n b a s e d o n t h i s m o d e l fi t t h e e x p e r i m e n t a l d a t a v e r y i i i 浙汀大学工程硕士学位论文 w e l l . t h i s m e t h o d p r o v i d e s u s a g o o d a n d r e l i a b l e w a y t o s t u d y t h e o p t i c a l c o h e r e n c e i ma g ing . k e y w o r d s : o p t i c a l c o h e r e n c e t o m o g r a p h y , mo n t e c a r l o s i m u l a t i o n , mi e t h e o r y , h i g h s c a tt e r i n g m a t e r i a l n 第 1 章 o c t原理 第， 章 绪论 1 . 1 引言 2 1世纪己经被认为是生物医学技术的世纪，新兴的生物基因技术、先进的 医疗探钡 诊断技术等生物医学技术已经开始对人们的生活产生重大的影响。 随着 对学科的交叉发展， 生医、 物理、 信息学等越来越多的专家， 科研技术人员开始 对生物医学技术领域投以更多关注。 借助现代电子技术和高性能电 子计算机技术 的发展，结合现代光学、电磁学、声学等理论，许多医学监测技术如:c t ，核 磁共振，p e t和o c t等先进的检测手段已经在实际的医疗诊断广泛的运用。为 了提高医生的诊断能力， 需要不同的检测手段从不同的侧面提供反映人体组织特 定参数的医学图像来描述组织生物变化的多层面相互补充的医学信息， 本文的工 作正是基于这样一种背景中产生的， 利用电 磁波在高散射介质中的散射吸收理论 结合现有的光学相干层析成像技术从而达到多层生物组织的高精度探测。 1 . 2光学相干层析成像的发展概况 1 9 9 1 年哈佛大学医学院的d a v i d h u a n g 博士等人在s c ie n c e 上发表了o p t i c a l c o h e r e n c e t o m o g r a p h y 1 一 文， ( o c t光学相 干 层析成 像技术) 宣布了 他们的 实 验结果。 他们利用一种类似麦克尔逊干涉仪的装置， 加上横向纵向扫描， 成功的 实现了眼角膜和眼前室的活体动态层析成像。由于o c t系统具有原理简单，高 分辨率 ( 0 1,) = r sr , p 。 二 一 (a 1pa k )z c o s(2 4 lok o) 4 ( 2 . 8 ) o c d r 或 者 o c t 系 统 的 设 计 用 于 高 精 度 测 量 群 延 迟 1k 当 i、 的 大 于 六 时，由 ( 2 .8 ) 可得， 最后的相干信号以非常快的 速度衰减。 可以 看出与用单色光检 测不同， 在轴向深度方向， 不会有两个不同的深度得到相同的结果， 所以我们可 以通过调整麦克尔逊两臂的长短， 影响两臂的群延迟， 最终达到纵向深度扫描的 效果。 一 一 1 f w h p _ _一 一/ g r o u p d e la y m is m a t c h a l . 圈 2 - 3纵向 分拼率 由 图2 - 3 我 们 可 以 看 到 通 常 情 况 下 信 号 的 分 布 ， 当 几 等 于0 的 时 候 信 号 到 达 峰 值 ， 调 整 两 臂 的 长 度 差 ， 达 到 扫 描 ig 的 效 果 】 可 以 看 到 当 1g 偏 离0 时 ， 信号迅速衰减至0 ，衰减的速度越快，获得信号就越代表局部的特性。一般，人 们定义信号衰减到峰值的一半 ( f wh m f u l l - w i d t h - h a l f - m a x i m u m)时，我们称这 个宽 度为 纵向 的 最小分辨长 度 今 、 ， 一 般情况下在8 1 0 - 8 3 0 的 红 外宽 带激光器 第2 章o c t 原理 的，。的达1 0 微米左右。可见o c t 是一种很高分辨率的探测器，这在生物组 织的监测和诊断中具有很大的优势。 2 2 2o c t 组织散射模型 o c t 一经提出，就引起很多的关注，它由于利用光的相干特性，可以提供 很高的分辨率。但是在实际的测量中，人们发现o c t 虽然具有很高的测量精度， 但是o c t 的探测深度很有限，在高散射的介质如皮肤组织探测深度只有2 毫米， 在透明组织如眼镜探测深度也只能达到2 厘米，这极大地限制了o c t 的适用范 围。究竟是什么原因导致o c t 的探测深度只有这么浅呢? 在上面的推导中，我们有几个假设：我们假设组织是多层结构的理想反射面， 光波在组织中传播，碰到一个反射介面即有一部分光被原路反射回来，其只影响 该光波最终和参考臂光波的时间延迟，对光波的其他特性的影响我们并没有考虑 在内。在现实而言，这个假设一般不能成立。 尽管把多层组织当作理想的多层镜面【2 4 】【2 5 反射极大的方便了对o c t 的理 论分析，但是大多数组织是光稠密介质，并不能很好的符合这个模型。软组织是 由胶原质和弹性蛋白纤维混合了许多细胞，血管，神经等其他生物组织的胶凝状 物质。这些组织的大小通常在1 0 0 纳米到几个毫米都有。当光被聚焦到组织上是， 组织的折射率在各项是异性的，就这样导致光向各个角度散射。各种光在总散射 光中的比例如何，它们对最后探测的结果又分别起到什么作用的，它们是如何影 响返回光信号的时间相干性和空间相干性的。由此一些研究小组开始研究和解决 这些问题，发展出一些更具普遍性的光学相干成像理论。 2 3 1 1 3 3 世界各地的研究小组发展了很多物理模型用于解释o c t 的理论，其中的一 些也成功地解释了某些物理现象，但还是没有一种物理模型能够比较完美的模拟 o c t 的所有物理现象。不过几年的研究结果，我们还是得到一些十分有用的结 论：a 组织的主要信息是由单次散射光携带( 图2 4 中的标号为1 的光) ；b 小 浙江大学t 程硕士学位论文 角度散射光( 图2 4 中的标号为4 的光) 是影响信号对比度的主要因素，而大角 度散射光( 图2 - 4 中的标号为2 的光) 却对结果没有很明显的影响；c 最大探 2 圈2 _ 4 光波在组织中的散射分4 种类型 测深度主要受组织的单次散射系数和平均散射角度的制约。一般类似皮肤组织等 高散射组织的最大探测深度在1 0 1 5 毫米( 约为1 3 0 0 n m 的光波在组织中的5 8 个平均散射自由程) ：d 光源的相干长度或者说相干时间极大地影响o c t 成像的 分辨率。上面这些结论在后面o c t 的设计和制造方面无疑具有很高的指导意义， 今年来越来越多的高性能o c t 已经丌始被运用于i 临床，在眼科等需要高精度探 测的临场诊断中起着重要的作用。 o c t 理论研究的难点就在于由于其探测对象生物组织的复杂性。生物组织 不同于普通的电磁波传播介质，简单的运用光学，电磁学中的折射率，电导率， 磁导率等物理参数远远不足以描述生物组织这种复杂的物质。企图运用传统光 第2 章o c t 原理 学，或者传统的电磁场理论m a x w e l l 方程从一开始就有很大的局限性。人们 开始寻求其它手段，其它更好的物理模型，用于描述光波在组织中的传播。 2 3o c t 实验 根据o c t 的原理图2 - 1 ) 我们搭建了一个简单的o c t 实验平台图2 5 ) 。 圈2 - 5 0 c t 实验装置照片 是实验装置的实物照片。实验使用一个1 5 3 0 1 5 7 0 纳米的宽带光源，光源频谱如 图2 - 6 所示，一个5 0 5 0 的2 2 光线耦合器，上图中央部分是一个一维高精度平 移台用作深度扫描，平移台上是参考臂的反射镜，右下部是样品臂和用于对准的 h e n e 激光器。这个装置暂时没有引入o c t 的二维横向扫描功能，现在可以实现 浙江大学1 ：程硕士学位论文 - 2 5 3 0 芒3 5 量 ： l “ 丑 一5 5 0 1 5 2 0 1 5 3 01 5 4 0 1 5 5 01 5 6 01 5 7 01 5 8 0 w a v ei e n g t h f n m l 圈2 - 6 竟带光源囊潜 最简单的1 维层析成像。相干信号传入上图左侧的光电探测器，将信号输入计算 机，用l a b v i e w 控制计算机上的a d d 卡对信号进行模数转换及相应的滤波处理 最后在屏幕上输出探测到的波形。 围2 - 7 简单多层样本 图2 7 是我们实验中使用的一个简单多层样本，样本的第一层折射率1 5 ， 第二层是空气，最后。个反射面是全反射镜。通过这个简单多层样本的纵向扫描， 我们可以得到它的纵向扫插图样图2 - 8 ) 。从图样中我们可以清楚的看到三个界 而的信号。由于实验中的扫描是由深至浅的，所以实验所得的信号中第个极 1 4 第2 章o c t 原理 uz4581 u1 z o p t i c a ll e n g t h 【c m 】 圈2 - 8 深度扫描圈祥 大值信号是由实验样本的反射镜反射的信号，其能量最强，另外两个略小的界面 信号则是介质到空气、空气到介质的两个约为4 的表面反射。而图2 - 9 是峰值 信号的放大信号，可以看到相干叠加的图样，相应的相干长度在介质中为1 0 0 微 米。 6 0 0 l e n g t h u m 】 图2 - 9 相干值信号图样 9 8 7 6 5 4 3 2 z，。口曼2 7 b 5 e、ooel石 浙江大学工程硕 七 学位论文 第3 章 电磁波在高散射介质中的传播 3 . ， 光散射理论 早在 1 8 0 2年 r i c h t e : 就曾观察到光束通过金溶胶时有光径的光散射现象。 1 8 6 9 年t y n d a l l 详细的 研究了 这种现象， 并使光线通过一个盛有亚硝酸丁酷和盐 酸混合热气的管子， 制成了模拟的蓝色天空。 他还发现在白光照射 卜 ， 散射光是 部分偏 振的。 1 9 7 1 年r a y l e i g h 详细的 研究了t y n d a l l 现象后指出 ， 散射光强 度与 入射光强度、 单位体积内的粒子数及微粒的体积平方成正比， 而与 光频率的四次 方成反比， 这就是 著名的r a y l e i g h 散 射定 律。 r a y l e i g h 散 射的 特点是散射光的 频 率不变，即光子的能量不变， 所以 又称为弹性散射。1 9 1 4 年b r i l l o u i n 开始对散 射光的频谱进行研究， 计算出散射体内存在声波引起介质密度涨落时散色射光中 的频率分布， 发现在入射光频率( r a y l e i g h 谱线) 附 近对称的 分布着分离的 谱线， c r o s 。 与 1 9 3 0年首先在液体中观察到b r i l l o u i n双线，这种由声波引起的光散射 称为b r i l l o u i n 散 射。1 9 2 8 年r a m a n y o n g 汞灯照射胶体苯时， 发现散射光谱中 有 比汞灯本身更多的谱线，这种由分子现晶格振动所引起的散射，称为 r a m a n散 射。 b r i l l o u i n散射和r a m a n散射对于入射光频率都发生了 频移， 所以 称为非弹 j性散射。 r a y l e i g h 散射中没有考虑吸收 和散射粒子的空间尺寸， 当 这种尺寸变得可以 和入射辐射的波长相比 较时， r a y l e i g h散射不再适用， 将会被 m i e散射取代。 m ie 散射理论相当复杂，与r a y l e i g h 散射相比 表现出 对波长更弱的 依赖性， m i e 散 射的 相函 数具有前h 性 ， 而r a y l e i g h 散 射的 前向 和后向 散 射强 度相同。 现在对于波在介质中的传播和散射问题的研究主要分为两种理论进行讨论， 一种是 “ 输运理论” 或称之为 “ 辐射传输理论” 3 4 - 3 7 ; 另一种是“ 解析理论” 或称之为 “ 多重散射理论” 。输运理论不是从波动方程出发，是直接讨论能量强 第 3 章 电磁波在高散射介质中的传播 度的传播。 输运理论中基本的微分方程称之为输运方程或者传输方程， 它等价十 气体分子动力学理论和中子输运理论中的b o l t z m a n n 方程。 该理论的研究比较直 观， 但它缺少解析理论那种数学上的严谨， 尽管在描述单个粒子的散射和吸收特 性使他考虑了衍射和干涉效应， 但输运理论本身并不包括衍射效应。 在输运理论 中假定辐射场之间不存在相关性， 因此用到的是功率的叠加而不是场的叠加。 尽 管如此， 它己经成功地被用于研究大气和水下能见度、 海洋生物学、 感光乳胶以 及辐射能量在行星、 恒星和银河系大气中的传播、 光学在生物医学中的应用等问 题。解析理论从诸如m a x we l l方程或者波动方程这种基本微分方程出发、引 进粒子的散射和吸收特性， 并求出方差和相关函数这些统计量的适当的微分方程 或者积分方程，原则上， 这种理论考虑了多次散射、 衍射和干涉效应。 在这个意 义上说，它在数学上是严格的。 但是， 实际上它不可能得到完全包括这些效应的 解， 产生有用解的各种理论也都是近似的， 每种理论只在一定的参数范围内有用。 光波在复杂介质中的散射理论尚在发展之中， 还没有一个统一的理论可以解 决所有的问题， 对于特定的问题需要灵活的采用特定的方法来解决。 如上面提及 的输运理论， 在非相干光成像中己经广泛的运用， 实践证明输运理论能够很好的 解释光子在高散射介质中的能量传播。但是输运理论只能用于描述光子的能量， 无法描述光场的一些诸如电场， 磁场等矢量信息， 所以很难用于解决相干光成像 的问题。 解析理论由于生物组织的复杂性， 无法使用常用的物理参数描述， 所以 也难以对相应的微分方程或者积分方程求解。 两个方法都具有优点和缺点， 有没有可能结合这两个方法的优点创造一种新 的办法来解决这个复杂的问题呢? 研究人员作出了很多尝试， 发展出许多新的理 论和方法， 诸如: 单次 散 射理论 3 9 、 线 性系统理 论2 3 . 蒙 特卡洛方法 4 0 、 扩 展的h u y g e n s - f r e s n e l 2 2 4 1 解析 模型 等等。 扩展的h u y g e n s - f r e s n e l 解析 模型 在 传统的生物组织散射模型中， 运用光学的h u y g e n s - f r e s n e l 原理，即当 前波的 波 前是一个小的二次光源， 这些光源发出的光波相互叠加形成下一个波的波前， 这 是一个经典的解析理论， 这样光子与光子之间不再是互不相关的个体， 它们之间 浙江大学工程硕十学位论文 会有相互作用，成功地在此模型内描述光的相千特性。 我们知道h u y g e n s - f r e s n e l 也是电 磁波理论的远场近似，如果结合经典的电 磁波理论， 用其更加精确的ma x we l l方程解析解， 不是可以得到更好的结果? 而且现阶段o c t 发展迅速， 扩展的h u y g e n s - f r e s n e l 无法描述光的 偏振， 相位等 特性，很容易想到应该寻找更加本质的物理模型。 互 3 . 2米氏 ( mi e )散射理论 现在我们把样品看作是一个悬浮着很多小球的均匀介质， 这种物理模型与大 多数悬浊液以 及o c t仿真实验中 经常采用的聚乙烯小球水溶液等高散射介质完 全一致， 在生物组织中这个假设也被证明是合理的。 做这个假设是为了满足mi e 散射理论 3 8 的要求。 mi e 散射理论是计算经平面波小球散射后的电磁场的解析解的理论。 平面波 延k 方向 入 射， 经 过 散 射后 原 坐标 系 绕k 转 过必 角， 再 绕e i 转 过b 角， 入 射 光 和 散射光坐标系之间的关系如图 3 - 1 所示。根据 mi e散射理论得到散射光和入射 光的关系 e k ( - , ) - i k r oz s, i ce , i ( 3 . 1 ) -工 哪习 /rleses、 其中 s , = e s 2 一 l r 2 n+l n ( n + 1 ) 2 n+1 n ( n + 1 ) ( a n n n 十 b n z n ) ( 3 .2 ) ( a n r + b o ffi n ) a n 和氏 是m ie 散 射的 散 射 场系 数，;7 n 和t o 是m ie 散 射的 角 度函 数 第3 章 电 磁波在高散射介质中的传播 圈 3 - 1 傲射前后局部坐标系关系 a = h= m z j n ( m x ) i x j ( x ) 一 p j ,( x ) l m x j ( m x ) l m j n ( m x ) x h , ( x ) 一 a h n . .( x ) m x j ( 二) n , j ( m x ) i x j . ( x ) 一 j ., ( x ) m x j ( m x ) a m-) l 戒 ( x ) l 一 彭 p ( x ) m x j . ( m x ) ( 3 . 3 ) m是小球相对小球周围媒质的相对折射率: x = k a 是尺寸参数， a 是小球的半径， k = 2? c 是 波 数 元 是光波 在介 质中的 波长;a , 是小 球的 相 对磁介电 常数， 一般 情 况f 等 于1 ; j ( z ) 和叼 ) ( x ) = j ( z ) + ry ( z ) 是n 阶 球 谐b e s s e l 函 数。 其中 的 两 个导数表达式分别为 zj ( z ) t = zj - , ( z ) 一 n j ( z ) , z h ,0 ( z ) = z h !, ( z ) 一 n h ) ( z ) ( 3 .4 ) 所以 为了 求a n 和瓦需 要求得大 ( z ) 和y ( z ) ， 而人 ( z ) 和y ( z ) 满足 递推关系 浙江大学工程硕士学位论文 人 - 1(z ) 十 f - , ( z ) = 2 n + l f n ( z ), 其中 f 可以 是 i n 或vn ( 3 . 5 ) 当n 等于 0 和 1 的时候有如下表达式 、 s i n ( 习_ 、 i n k z ) = , hk z l =s i n ( z ) z zz c o s ( z ) y o ( z ) = z c o s ( z ) ， y , ( z ) = 一 c o s ( z ) s in ( z ) ( 3 . 6 ) 2 ，乙 艺 至此当我们知道小球的相对折射率m，小球的直径a ， 和光波在介质中的波长a 后， 根 据以 上 ( 3 .3 ) ( 3 . 6 ) 的 一系列关系式， 我 们求解 得到不同n 对 应a 。 和b o 、 。 和: 。 是6 的 函 数 所以 称 之为 角 度函 数， 也 是由 两 个 递 推 关 系 2 n - 1 n一t c o s 口 气- 1 - n - 了 r n - 2 1 t n = n c o s b - ; r n 一 ( n + 1 ) i r n - , 其中 ;t u = 0 , )r , = 1 ， 二 ， = 3 c o s b ; ( 3 . 7 ) t o = 0 , t , = c o s b , t 2 = 3 c o s 2 b 所以在一开始得到散射系统基本参数m, a , 兄的基础上， 每次散射的发生， 我们都可以算出其在各个散射角b 下得电场分布。在实际的计算中， n 可以取前 儿项。 在实际的使用中， 我们也可以用s t o k e s 量来表示。 s t o k e s 是能个在实验中被 测定的物理量，用 s t o k e s量可以很方便的表示偏振光，所以在实际运用中，比 直接操作无法测量的电场矢量要容易得多。 ee , + ee , ee ;, 一 e , e , e e l + ee ,* i ( e e l 一 e , e ii ) ( 3 . 8 ) 声一、 -一 、. 了quv 一、 一- s 而入射光和散射光的s t o k e s 量之间同样满足一个关系式 第3 章 电磁波在高散射介质中的传播 ( 3 . 9 ) 去9鱿犷冲 一 一 oom3毯 nun 舰脚 m3 3 -m3 今 门o 叭mloo 1一洲 - 159认代 m il = i (is 2 i + is i), u2一 合 (is l 一 is 12) m 33 = 合 ( a + s 2s i. )l m 34s 2= 告 ( as 2 一 、 s ), ( 3 , 1 0 ) 在小球模型的基础下，mi e 散射理论是ma x we l l方程在此条件下的精确 解析解。 在后面的章节里我们可以看到， 当在组织或者介质颗粒半径接近光波波 长时， 我们不能再认为光的散射场是各向同性的了， 通过mi e 散射理论， 我们可 以精确的得知散射场在空间的分布 4 2 1 ，从而可以 在理论上更贴近实际情况。 浙江大学工程硕士学位论文 第4 章 蒙特卡罗模拟方法 4 . 1 蒙特卡罗方法 蒙特卡罗( mo n t e c a r l o ) 方法的定名和系统发展始于二十世纪四十年代中期。 约在1 9 4 6 年， 物理学家v o n n e n m a n n 等在电子计算机上用随机抽样的方法模拟 了 种子连锁反应， 并把这种方法称为m o n t e c a r l 。 方法。 m o n t e c a r l o 方法以 概率 统计理论为主要理论基础，以随机抽样为其主要手段。 它的基本思想是， 为了求 解数学、 物理、 工程技术以及生产管理等方便的问 题， 首先建立一个概率模型或 随机过程， 使他的参数等于问题的解; 然后对模型或者过程的观察或者抽样实验 来计算所求参数的统计特征， 最后给出所求的近似解， 二阶的精确度可用估计值 的标准误差来表示。 假设所要求的 量x 是随 机变量咨 的 数学期 望e ( 约， 那么 近似确定x 的 方法是 对咨 进行n次重复抽样，产生相互独立的睿 质的序列氢 ， 鑫 ， ， ， . ，氛，并计算算数平 均值s n =1 nn y c ， 根 据 柯 尔 莫 哥 罗 夫 加 强 大 数 定 理 有 p ( 息 n 一 x ) 一 1 ， 因 此 当 我们取足够多的n后，所求得算数平均值以极大的概率等于我们所要求得量x . mo n t e c a r l o 方法可以解决各种类型的问 题，但总的来说，视其是否涉及随 机过程的性态和结果，用 mo n t e c a r l o 方法处理的问题可以分成两类: 第一类是确定性的数学问题，用mo n t e c a r l o 方法求解这里问题的方法是， 首先建立一个与所求解相关的概率模型， 是所求的解就是我们所建立模型的概率 分布或者数学期望; 然后对这个模型进行随机抽样观察，即产生随机变量; 最后 用其算数平均值作为所求解的近似估计值。 计算多重积分、 求逆矩阵、 解线性代 数方程组、 解积分方程、 求某些偏微分方程边值问题和计算微分算子特征值等都 属十这一类。 第4 章 蒙特卡罗模拟方法 第二类是随机性问题，例如中子在介质中的扩散等问题就属于随机性问题， 这是因为中子在介质内部不仅受到某些确定性的影响， 而且更多的是受到随机性 的影响。 对于这类问题， 虽然有时可表示为多重积分或某些函数方程， 并进而可 考虑用随机抽样方法求解， 然而一般情况下都不采用这种间接模拟方法， 而是采 用直接模拟方法， 即根据实际物理情况的概率法则， 用电子计算机进行抽样实验。 原子核物理问题、 随机服务系统中的排队问题， 动物的生态竞争、 传染病的蔓延 还有输运问题都属于这一类。 我们研究的问题基本上就是随机性问题， 通过模拟大量的光子在高散射介质 中的随机过程，来求得最后接受面的光场分布。mo n t e c a r l o 方法使用起来很简 单，只要能对一个光子的随机过程完整的描述，通过简单的随机重复那个过程， 即能求得问题的解，所以mo n t e c a r l o 方法对空间要求是很小的。但由于随机过 程的 收 敛 是 与习 气 万同 阶 的， 要 把结 果 的 精 度 提高1 0 倍， 模 拟次 数 要 增 加1 0 0 倍， 所以m o n t e c a r l o 方法是一个非常耗时的模拟方法， 对计算机的运算能力要求比 较高。不过m o n t e c a r l o 模拟的随机过程与过程之间具有很好的并行性，随着现 在计算机性能和计算机阵列的发展， mo n t e c a r l o 模拟算法可以 很方便的移植成 为并行程序运行在计算机阵列上进行高速并行运算。 4 .2 多 层组织的欲特卡罗 模型 m c m l ) 1 9 8 3 年， wi l s o n 和a d a m首次使用m o n t e c a r l o 模拟方法模拟激光和组织的 相互作用。随后mo n t e c a r l o 模拟方法被广泛的运用与模拟光和组织相互作用。 1 9 9 5 年， l i h o n g w a n g 用c语言实现m o n t e c a r l o 模拟光子在多 层组织中 传输， 并 在i n t e r n e t 上公开了 源代码。 由 于c 语言良 好的 跨平台 移植性能， l i h o n g w a n g 的m c ml程序己经成为多层组织m o n t e c a r l o 模拟的 标准模版，人们纷纷下载 他的源程序，并在mc ml的基础上加以修改来满足自己的特殊需求。 用mo n t e c a r l o 方法模拟光子在组织中的传输具有灵活性和鲁棒性。它能够 浙江大学工程硕士学位论文 同时获得多个物理量的模拟结果。 在这里光子在组织中的传输模型十分简单， 山 z 个随机过程组成，一 个描述了光子在两次散射过程中走过的步长，另一个过程 解决光子在散射过程中偏转的角度。 这种物理模型基于一些假设: 虽然组织复杂 多变的混合物，但是在宏观上其表现出统一的光学特性;在mc ml里光子被认 为是经典的粒子， 光的偏振和波的特性被忽略在外; 最后， 假设组织是各项同性 的。 4 . 2 . 1 m c m l里的坐标系 m c m l是多层组织最简单的 模型， 一束无穷细的光线垂直于入射在一 个多 层结构的组织表面。 每一层组织在x y方向上是无穷大的， 至少每层的横向尺度 远远大于它的厚度。 每一层我们用以下几个物理量描述: 厚度， 折射率，吸收系 数p a ( c m - , ) ， 散 射 系 数从 ( c m - , ) ， 和各 项 异 性因 子s 。 我 们 还需 要 周围 介 质 的 折射 率，比 如空 气 和基 板的 折 射率。 吸 收系 数p o 代表的 单 位长 度光 子 被吸 收 的概率，散射系数,u s ( c m - ) 单位长度光子被散射的概率，总作用系数 从= p a + p l ， 代 表单 位长 度 光子 参 加 相 互 作 用的 概率。 各 项 异 性因 子9 是 散 射 角 度余弦的平均值。 m c m l里面同时使用3 个坐标系，一个笛卡尔坐标系，用于追踪光子的 位 置，坐标系的原点即垂直光束和组织表面的交点，z轴和光束的入射方向一致， 垂直组织表面指向组织内部，而x y平面落在了组织表面上 ( 如图 4 - 1 ) . 因为该模型是无穷细的光束垂直组织表面入射，所以这个问题是圆柱对称 第4 章 蒙特卡 罗模拟方法 厄甲叭u p ,一压.i_l一 一 1仁 l la-ayer n n 1一 j . 口 胡一一 . . - 月一 三 二 ， 一 . . . 月， 二- 圈 4 - 1 多层组织结构的 全局笛卡尔坐标系， 丫 轴垂直纸面向 外 的。 因此我们建立起一个柱坐标系， 可以获得组织内部光子被吸收的分布， 该分 布是深度z 和离轴距离: 的函数， 柱坐标系和笛卡尔坐标系共享的z轴方向和同 一个原点，r轴在x y平面内。 另外还需要一个移动的球坐标系， 其z轴总是与光子的前进方向一致， 用于 对散射角度随机采样( 如图 3 - 1 所示的入射光和散射光的球坐标系之间的关系) 。 在这个球坐标系中， 偏转角b 和 方 位角0 在散 射过 程中 被随 机 采样， 根据图 3 - 1 所示，使用坐标变换可以得到新的光子前进的方向。 最后物理量的获得，通过按在接受表面， 或者组织内部划分的网格，分别统 计得到物理量在某个维度上的变化趋势。 互 4 .2 . 2 模拟光子传翰 浙江大学工程硕士学位论文 图 4 - 2 典型光子在介质中 传播的过程 在这节我们提出mo n t e c a r l o 模拟中光子在传输遵循的规则。我们的模型主 要是针对光子在多层介质的， 当然多层玻璃或者其它非高散射复杂多层介质时最 极端的例子。图 4 - 2 展示的是模拟过程中光子运输的典型过程，垂直表面入射， 经过一系列的散射后最终从表面出射。图 4 - 3 则是mc ml程序的流程图。下面 我就来看一下mc ml的模拟流程。 第4 章 蒙特卡罗模拟方法 图 4 - 3 m c m l 流程图 ( $ 是光 子每步 要走的 步长 z 7 浙江大学工程硕十学位论文 a . 舫次萝奋动买 # m o n t e c a r l o 算法的核心是对已知概率密度分布的随机变量进行采样。对于 随机变量x 在区间( a , b ) 上的概率密度分布是p ( x ) ，这个变量可以 是m c m l中 的每步光子走的步长， 也可以是每次散射光子偏转的角度， 或者是其他随机过程 中的问题。 计算机可以提供给我们一个伪随机数， 这是产生一个随机数最廉价的 方法。 虽然产生的是伪随机数， 但是在翻转周期足够长的情况下，我们可以认为 在有限次的伪随机数采样过程中，采样的前一次和后一次的结果是相互独立的， 在这种条件下， 我们利用计算机的伪随机数产生代替随机数的产生。 计算机可以 产生的随机数通常是在( 0 , 1 ) 得范围内 均匀分布。下而我们就利用这个( 0 , 1 ) 的 均 匀分布来产生( a , b ) 区间l 的ax ) 分布。 由 于ax ) 是x 在( a , b ) 上概率 密度函 数， 所以 f p (x )d x 一 ( 4 . 1 ) 定义p ( x ) 的 分布函 数 f , ( x ) 0 , k x p ( ic ) d ir , 1 , xa a _ 2 b ( 4 . 2 ) 厂leseseses少、weweesesse 工- 假设聋 是( 0 , 1 ) 上均匀分布的随机变量，显然它的分布函 数 乓( ,; ) 否 0 0 1 ( 4 . 3 ) 0，么1， 11、 工一 由 于 分 布函 数的 特点 ， 它 是 一 个 单 调 递 增 函 数， 所以 可以 证 明 对 于凡( 劝e 0 , 1 和 乓( 幻e 0 , 1 ， 对于 任意随 机参 数x , e ( a , b ) ， 一 定 可以 找到 一 个自 ( 0 , 1 ) ， 使 得 f ( x 卜乓( 5 j 成立， 而且x和4 . 是一一 对 应的， 将( 4 .2 ) ( 4 .3 ) 代入可得 ( 4 . 4 ) 第4 章 蒙特卡罗模拟方法 f p ( x )d : 一 ， 对 任 意 么 。 (。 ， : ) 都 成 立 ， 整 理 表 达 式 可 得 f p ( x )d : 一 ; ，; e ( 0 , 1) ( 4 . 5 ) 图 4 - 4由 均匀分布的艘机变，产生任宜的随机变f 由 式( 4 .5 ) 可以 看到一 个( 0 , 1 ) 上均匀分布的随 机变量咨 和在( a , b ) 上按p ( x ) 分布 的随机变量x 用一个积分方程联系 在了 一起， 只要求解这个积分方程我们就可以 得到所需要的任意分布的随机变量，图 4 - 4 用图示的方法展示了映射的过程。 b . 芳子含内璐逻 数据结构在程序中起了很重要的作用， 在 c语言中 把有逻辑上联系得参数 用结构 s t r u c t u r e )来表示有利于程序的可读性和维护性。那我们就来看一下 m c ml的源代码里使用怎么样的结构来表示一个光子包的。一 卜 面是程序中一段 定义光子包结构体的源代码: 浙江大学工程硕士学位论文 ty p e d e f s t r u c t d o u b l e x , y , z ; d o u b l e