基于蒙特卡罗光传输算法的激光消融生物组织温度场分布计算

王亚芬，白景峰

1 上海交通大学生物医学工程学院生物医学研究所，上海市，200030

2 上海交通大学Med-X研究院，上海市，200030

基于蒙特卡罗光传输算法的激光消融生物组织温度场分布计算

【作 者】王亚芬1,2，白景峰1,2

1 上海交通大学生物医学工程学院生物医学研究所，上海市，200030

2 上海交通大学Med-X研究院，上海市，200030

利用蒙特卡罗方法模拟了无限细光束在生物组织中的传输分布情况，通过卷积计算了有限宽光束光分布，在此基础上建立了生物组织传热模型，并根据Pennes方程，利用ANSYS求解得到了模型中生物组织温度场分布结果，通过分析激光功率及激光光束半径对生物组织热效应的影响，得到了对于临床激光消融有重要意义的结论。

蒙特卡罗；有限宽光束；卷积运算；生物传热

激光诱导间质热疗[1]利用激光照射肿瘤病变组织，使病变组织温度升高，最终组织凝结坏死，达到治疗目的。在实际治疗过程中，激光参数的选择对于治疗效果非常重要，包括激光功率、激光照射时间、激光光束半径等。为了得到有效治疗参数，对生物组织中的激光分布及因光分布引起的组织热分布进行数值计算非常必要[2-3]。李小霞等[4]假设激光在生物组织中呈指数衰减分布，依据Pennes生物传热方程，建立传热模型，得到生物组织光热分布，江世中等[5]在未考虑血液灌注项的前提下，利用蒙特卡罗方法，仿真计算了激光在生物组织中的光分布，并通过Pennes方程采用有限差分法求解了离体肝脏组织的热场分布。本文在蒙特卡罗方法计算无限细光束能量分布的基础上，利用有限元方法，计算了一定束宽的高斯光束照射牛肝脏组织后，组织的温度场分布情况。

1 理论模型

1.1 蒙特卡罗仿真光分布

描述激光的传输问题，较常用的理论模型是传输模型，即忽略光的波动特性，考虑光子在传输过程中，被散射或者被吸收。蒙特卡罗方法从概率统计的角度出发，跟踪光子的运动轨迹，通过大量样本的叠加统计，得到光在介质中的能量传输结果[6]。在模拟过程中，蒙特卡罗方法以传输介质的吸收系数(μa)和散射系数(μs)为根据，通过相函数控制光子碰撞后的散射方向。

蒙特卡罗模拟光子的运动主要分为以下几个步骤：首先，光子发射的初始位置设为坐标原点，光子垂直入射至生物组织，每个光子的初始权重为1（光子被吸收后，权重减少）；其次，光子在生物组织中被散射后的运动方向通过相函数产生，相函数与生物组织的各向异性因子有关，而运动步长则由随机数产生，光子散射改变光子的运动方向，光子吸收仅改变光子的权重大小；最后，光子的权重减少至0时，光子死亡，开始新的跟踪[7]。

利用以上方法模拟了生物组织的无限窄光束的脉冲响应，但是在利用激光治疗肿瘤时，用到的大多为高斯光束，所以还需要计算高斯光束光分布情况。由于用于治疗的生物组织可以看成一个线性系统，根据信号系统响应可知，在知道脉冲响应的前提下，可以通过求卷积得到任何输入的输出结果。

假设无限窄光束的光分布函数为G(x, y, z)，有限宽光束的光束形状为S(x, y)，则它在生物组织中的光束分布函数T(x, y, z)可以用式(1)表示。

高斯光束形状S(r)=S0exp(-2(r/ω0)2)，其中S0是光束在r=0 cm处的光强，ω0是激光强度下降到1/e2的距离，即光束束腰半径，由于高斯光束具有轴对称特性，因此将高斯光束表达式及式(2)代入式(1)，并转换为柱坐标系统，得到：

最终得到：

其中I0是贝塞尔函数，由下式给出：

因此，在用蒙特卡罗模拟得到无限窄光束的能量分布后，可以通过上式，利用数值算法求解得到高斯光束的分布。

1.2 生物组织传热方程

在求得生物组织光分布，即外加热源后，采用经典的Pennes方程[8]作为二维柱坐标生物传热理论模型：

其中，ρ是生物组织的密度，c是生物组织的比热容，T是生物组织温度，k是生物组织的热导率，Cb为血液的比热，Wb为血液的灌注率，Tv为血液的温度，等式右侧第二项构成生物组织中由于血管存在引起的热量对流项，Q为外加热源。本文模拟计算的对象是牛肝脏组织，在仿真计算过程中，未考虑血液灌注的影响。Pennes方程的数值求解方法有有限差分法和有限元法，本文采用有限元方法，利用ANSYS软件，通过APDL编程施加蒙特卡罗热源载荷，求解方程。

1.3 生物组织模型及求解参数

激光在生物组织中的能量分布是轴对称的，本文采用二维柱坐标系统进行仿真，模型如图1所示，网格为划分的有限元单元，r为激光照射时半径方向，z为轴向。计算范围为2 cm×2 cm，网格大小为0.005 cm×0.005 cm。分别计算了功率为1 W、2 W、 4 W时的温度场分布情况。计算时边界条件设置为z轴绝热，照射表面即z=0处为自然对流，对流系数为数6.2×10-4W/(cm·oC)。

图1 生物组织模型示意图Fig.1 Geometry model of tissue

组织模型对象为牛肝脏，激光波长为1 064 nm，激光照射时生物组织光热学参数参考文献[4]，吸收系数μa为0.5 cm-1，散射系数μS为80 cm-1，生物组织各项异性因子g为0.9，密度ρ为1.093 g/cm-3，比热C为3.37 J/(g·oC)，热导率K为4.5 mW/(cm·oC)。

2 数值分析结果

蒙特卡罗仿真高斯光束分布，跟踪记录了105个光子的运动轨迹，根据Pennes热传导方程，对激光热疗时的生物组织温度场进行建模仿真，模型求解工具ANSYS。初始温度为37oC，照射时间为300 s，计算时间600 s。

图2是P=2 W时，光束半径为0.05 cm，激光照射时生物组织在不同时刻的温度场分布仿真结果，(a)～(d)分别为t=100 s、300 s、400 s、600 s时组织温度场分布情况。

图2 P=2 W不同时刻生物组织温度场分布Fig.2 Temperature distribution at different time for 2 W power

图3是P=2 W时生物组织表面(r = 0 及 z = 0)沿径向及轴向的温度分布情况，从图2、图3可以看出，随着照射时间的增长，生物组织内部热传输范围扩大，照射300 s时，组织温度扩散范围在r、z方向均不超过1.5 cm，组织产生热损伤范围较小。

图3 P=2 W生物组织表面温度沿径向(a)及轴向(b)的变化Fig.3 Temperature distribution in r (a)and z (b) direction at tissue surface for 2 W power

为了分析功率大小对温度场分布结果的影响，分别仿真计算了激光功率为1 W、2 W、4 W时，高斯光束的束腰半径为0.05 cm时，生物组织的温度场分布，如图4、图5所示。图4为不同功率激光生物组织温度场变化情况，图4(a)、4(b)为激光功率是1 W时在t=300 s及t=600 s时组织温度等温线分布结果，图4(c)、4(d)为功率是2 W时的计算结果，图4(e)、4(f)为功率是4 W时的计算结果。

图4 不同功率下生物组织温度场分布Fig.4 Temperature distribution for different laser powers

图5为不同功率激光(1 W、2 W及4 W)照射后，照射原点(r = 0，z = 0)处组织温度升高大小随时间变化情况。在t=300 s时停止照射没有外加热源后，组织温度迅速下降。

图5 不同功率下生物组织在 r =0，z = 0处温度随时间分布Fig.5 Temperature distribution vary with time for different laser powers at r =0 and z =0

为了分析能量密度对温度场分布结果的影响，计算了激光功率2 W时，束腰半径分别为0.05 cm、0.08 cm及0.1 cm时组织温度分布情况，图6为激光照射300 s，在t =300 s时刻生物组织表面(r = 0，z = 0)温度升高大小沿径向和轴向的变化情况。

图6 t=300 s不同光束半径下生物组织表面温度沿径向及轴向的变化Fig.6 Temperature distribution in r (a) and z (b) direction at tissue surface for different laser radius at 300 second

3 结论

基于蒙特卡罗方法的光分布理论，建立了生物组织传热模型，并利用ANSYS有限元软件求解了Pennes生物传热方程，得到了生物组织的温度场分布结果，通过结果得到了以下结论：功率大小对组织的损伤区域影响低于其对温升的影响，功率越大，温度变化越迅速；高斯光束照射生物组织时，照射初期生物组织温度沿径向和轴向有差别，照射一定时间后，由于能量传递，径向和轴向温度分布差别较小；在功率相同的情况下，激光束腰半径对照射点温升有影响，半径越小，能量越集中，温升越高。

[1] Bown SG. Phototherapy of tumors [J]. World J Surg, 1983, 07(6): 700-709.

[2] Puccini S, Bar NK, Bublat M, et al. Simulations of thermal tissue coagulation and their value for the planning and monitoring of laser-induced interstitial thermotherapy (LITT)[J]. Magn Reson Med, 2003, 49(2): 351-362.

[3] 谢树森, 龚玮, 李晖. 生物组织的选择性光热解效应[J]. 激光与光电子学进展, 2004, 41(8): 48-51.

[4] 李小霞. 激光照射下生物组织热效应的数值分析与实验研究[D]. 天津大学, 2004.

[5] 江世臣, 张学学. 表面照射下激光与生物组织的光热作用分析[J]. 光电子·激光, 2005, 16(06): 752-756.

[6] Wang LH, Jacques SL, Zheng L. MCML-Monte Carlo modeling of light transport in multilayered tissues[J]. Comput Meth Prog Bio, 1995, 47(1): 131-146.

[7] Kienle A, Patterson MS. Determination of the optical properties of turbid media from a single Monte Carlo simulation[J]. Phys Med Biol, 1996, 41(12): 2221-2227.

[8] Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm[J]. J Appl Physiol, 1948, 1 (2): 93-122.

Temperature Distribution Based on Monte Carlo Method of Optical Transmission in Tissues of Laser Ablation

【 Writers 】Wang Yafen1,2, Bai Jingfeng1,2
1 Biomedical Instrument Institute, School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai, 200030
2 Med-X Research Institute, Shanghai Jiao Tong University, Shanghai, 200030

Monte Carlo, fi nite-diameter laser, convolution operation, bio-heat equation

R318.51

A

10.3969/j.issn.1671-7104.2013.04.005

1671-7104(2013)04-0252-03

2013-01-13

白景峰，E-mail: jfbai@sjtu.edu.cn

【 Abstract 】Monte Carlo method was used for calculation of finite-diameter laser distribution in tissues through convolution operation. Photo-thermal ablation model was set up on the basis of Pennes bioheat equation, and tissue temperature distribution was simulated by using finite element method by ANSYS through the model. The simulation result is helpful for clinical application of laser.