漫射光层析成像的正问题理论模型(一)—— 蒙特卡洛模型
组织光学领域的根本问题是光在生物组织中的传输问题。目前有三种理论来研究光在介质中的传输与散射:一种是基于随机数思想的离散统计理论,也就是Monte Carlo方法,它主要通过描述单个光子与组织的相互作用来模拟光在生物组织中的吸收和散射,然后通过发射大量的光子来得到宏观的光强、反射率、透射率等物理量从而描述光在介质中的传输行为。第二种是通过求解辐射传输方程(Radiance TransportEquation,RTE)来描述光在介质中传输的辐射传输理论。第三种是通过对RTE做一定的近似处理从而使时间效率得到大幅提高的漫射方程理论,它将解析理论中对RTE的求解通过做漫射近似(DA,Diffusion Approximation)转换为对漫射方程(DE,Diffusion Equation)的求解从而减少计算量加快计算速度。通过这三种方法,我们可以推导出在特定光照条件下确定光学参数分布的生物组织中的漫射光强分布及其他物理量的理论值,从而描述光在生物组织中的传输行为。
Monte Carlo(中文译为蒙特卡洛)方法是一种非确定性方法,在数值模拟领域,我们通常用Monte Carlo方法来表示随机实验。Monte Carlo是地中海小国摩纳哥的第一大城市,以Monte Carlo赌场闻名于世。20世纪40年代,John Von Neumann,Stanislaw Ulam和Nicolas Metropolis在洛斯阿拉莫斯国家实验室为核武器计划工作时,用随机抽样的方法在计算机上模拟了中子连锁反应,因为Ulam的叔叔经常在Monte Carlo赌场输钱,他们把这种方法叫做Monte Carlo方法。
Monte Carlo方法是一种基于统计学的数值计算方法,以随机抽样为主要手段,其核心思想是使用随机数或更为常见的伪随机数来解决一些计算问题,很适合用在一些复杂系统的仿真中,例如空气动力学计算或者宏观经济学中。与那些通过求解一组不同的方程得到确定的解来解决问题的确定性方法不同,MonteCarlo方法是通过发射大量的相互独立的随机事件,并对其结果进行分析从而得到所求解问题的概率分布,再用统计方法得到这种概率分布的数字特征,从而进一步得到实际问题的数值解的一种方法。由于与其他的一些确定性方法相比,Monte Carlo方法在某些条件下的求解更为方便、快捷和灵活,所以他已经被广泛应用在粒子输运计算、量子热力学计算、金融工程学计算、生产管理以及生物医学等领域中。
在使用MC方法来解决相关问题时,首先要对所求解的问题建立一个便于实现且具有一定意义的概率统计模型,这个模型的概率分布能够反映我们要求解的问题;然后确定随机变量的产生方法,也就是如何根据这个模型对随机数进行抽样,包括随机数或伪随机数的产生以及根据概率统计模型的特点建立起产生随机变量的抽样方法;接着按照实验的规则,进行成千上万次相互独立的随机抽样实验,得到实验结果,分析其数学特性,再用统计的方法求出所求解的数学期望、方差标准差等值,从而对问题中的各个变量的关系进行研究。
模拟光子在混浊介质中的传输也是Monte Carlo方法被证明适合使用的一个领域,尤其是在生物成像领域的应用已经有多篇文章发表在国际知名的期刊和杂志上,例如人脑功能成像、乳腺癌的检测、以及小动物的成像等。凭借其在实际应用中的准确性,灵活性和通用性,Monte Carlo方法已经成为解决这些问题的标准方法,尤其是在低散介质(散射系数较小)中,其误差相比于其他方法更小,甚至可以用来作为评价其他模型的标准。光在不均匀的人体组织中的传输的模拟需要一个准确而高效的模型,Monte Carlo方法能够满足这样的要求。
使用MC方法来模拟光在介质中传输的原理是用随机数来描述每次光在组织中传输时的步长的变化以及由于散射引起的散射角的变化,从而将确定性问题转换为随机问题。在向组织发射大量的光子,通过模拟每个光子在组织内部的按照与其宏观物理行为相一致的规律和组织随机发生相互作用,通过对这些大量的随机事件的追踪和统计来获得我们所需的宏观的物理量如光强等。MC光传输模型是忽略光的波动性如干涉和衍射效应的。
如图1所示,MC仿真的流程为:首先生成光子并为其分配随机数,这个随机数称为这个光子的种子,发射这个光子,由这个随机数可以得到光子的步长即下一次散射会走的距离,由光子现在所处的位置以及它的步长可以知道光子的这次散射是否会到达样品的边界,也就是说光子的下一个位置是否还在样品内部。如果没有到达边界,则使光子移动相应的步长到下一次发生散射的位置,确定光子由于吸收而引起的强度衰减,并确定下一次散射所引起的方向变化。如果到达边界则先将光子移动到边界上,计算由这段吸收造成的强度衰减,再确定是反射还是透射,如果反射则进行下一次的散射行为,如果透射则停止对这个光子的追迹。如果光子的权重小于阈值或者光子已经出射,则对这个光子的追踪结束,开始追踪下一个光子。当所有的光子都被发射并追踪完毕后,模拟结束。
整个过程可以总结为:光子的产生,轨迹的产生,吸收,消亡和检测。
图1:MC仿真流程图
Monte Carlo模型的优点思想极为简单,无需建立辐射传输方程RTE或输运方程DE,也无需为求解方程而煞费苦心地需找各种近似方程或数值解,只要知道组织的基本光学参数及其分布,Monte Carlo模型就可以根据追踪各个光子在组织中的传输来计算出漫射光强等各种物理量理论值的时间空间分布,从而描述光子的输运过程,但是它有一个致命的缺点,就是计算时间太长。Monte Carlo方法以概率论为基础,实质上是使用计算机来模拟很大数量(如
)的光子随机移动,然后用统计的方法来得到结果,其仿真输出的误差正比于随机变量的标准差,因而
反比于取样次数(模拟中发射的总光子数)。
光在介质中传输时其强度会随着传输距离的增大而指数衰减,所以为了达到同样的相对误差,光源与探测器的距离越大,所需的计算时间t也越长。也就是说,在相同的预设精度条件下,介质的尺寸越大,MC仿真所需要的时间也越长并且是一个指数增加的关系,这才是MC仿真最为耗时的地方。所以,提高MC仿真的时间效率,加快它在预设精度下的收敛,对于MC仿真的实用意义重大。
近年来,为了改善MC方法的性能,科学家提出了许多技术手段或者理论模型来加快MC仿真:例如在硬件上对MC仿真进行加速的基于GPU(GraphicProcessorUnit,图形处理单元)的MC,它利用近年来快速发展的多处理器内核以及多线程并行计算理论在硬件上对MC仿真进行加速,通常来说,基于GPU的MC方法可以比基于CPU的MC快几百倍到上千倍之多。还有在算法上对MC方法进行改进的白蒙特卡(WMC,White,Monte Carlo)和扰动蒙特卡洛(PMC,Perturbation Monte Carlo)方法,它利用光子种子的概念和对种子的一些处理来对MC进行算法上的改进,经过改进后的MC方法相比于传统的MC仿真加速比也可以达到几百倍到上千倍。有了这两种手段,MC方法的计算速度得到了有效的提高,向实用化的方向又迈出了一大步。
虽然这些手段使得MC方法的时间效率得到了大幅度的提高,计算时间大大减少,但是与其他方法如求解RTE或DE相比,它仍然是一种较慢的方法,如对 一定尺寸大小的人脑成像,用计算机通过MC来仿真可能需要一个小时或更多的时间,但求解DE通常只需要几秒钟的时间,而求解RTE一般来说也只需要几分钟。所以MC方法仍然需要继续加速才能够更加适用于日常的临床诊断和治疗。
还有一些其他的比较通用的用来解决光在非均匀介质中传输的方法,例如求解辐射传输方程RTE,虽然它在数学上是严格的,但是实际应用时它的求解过于复杂,计算规模巨大,且只有少数情况下才能够得到精确解,其他情况无法直接求解,所以一般不用它来作为解决光在混浊介质中传输的正问题模型。而对RT做了输运近似的漫射近似理论虽然从效率上相比RTE有了很大的提高,但是它只适用于吸收相对较弱而散射较强的介质
,当吸收增强如以
时漫射近似理论的误差会显著变大。因此它的时间效率是以求解精度为代价的,只有在高散射介质中它的误差才能与其他方法相当。但是实际使用时有很多生物样品并不能满足高散射低吸收的这个条件,例如人脑的脑脊液层,乳腺的囊肿区等,因此它并不适用于这些领域的研究。所以综合效率以及精度的考量,Monte Carlo方法是一种性能良好的可以用于描述光在介质中传输的正问题模型。
郭则飞. 蒙特卡洛仿真在漫射光层析成像中的应用[D].