当前位置： >>
放射治疗剂量计算的蒙特卡罗方法研究
摘甓摘要放射治疗是治疗肿瘤的主要方法之一，剂量计算是放射治疗计划系统的重要组成部分，其准确性对于放射治疗的效果影响很大。长期以来临床上使用的 荆量计算模型在精度上还未能完全令人满意。蒙特｝罗方法擅长处理粒子输运问题，大量实验证明蒙特卡罗方法是计算精度
最高的剂量计算方法，但存在计算时间长的问题。本文分蔓个部分对蒙特卡罗剂量计算方法进行了系统研究。第一部分主要研究适用于放射治疗计划系统的蒙特卡罗剂最计算模型。介 绍了电子和光子耦合输运的蒙特卡罗模拟原理；详细研究了模拟电子输运的浓缩历史方法和多次散射理论；讨论了基于直线加速器的辐射源模拟和多源模型的使用；详细研究了根据ＣＴ值推导人体组织参数的多种方法及存在问题；给出了基于二次蓝面描述射线束修正设备的方法；阐述了人体内部蒙特常罗剂量 分布计算的特点。 第二部分主要研究ｒ蒙特卡罗剂量计算的快速算法，这是～个国际上的研究热点。在对前人成果深入研究的基础上，总结出提高蒙特卡罗剂量计算速度的四个基本原则，即空间换时间原则、降低计算方差原则、物理模型简化原则 和最大化利用计算资源原则。然后以此为基础，对电子输运、光子输运、多源模型和射线束修正设备模拟等四个方面的快速模拟进行了系统研究。在这里， 我们提出了可以提高模拟效率的描述复杂射线柬修正设备的离散网格方法，给 出了快速模拟光子和电子通过射线柬修正设备的方法以及对人体周围空气的简化模拟方法。第三部分重点研究了向量蒙特卡罗模拟方法，提出了一种基于ＳＳＥ技术的 电子输运向量化模拟算法，并将其推广到光子输运模拟；我们用这一方法实现 了快速蒙特卡罗剂量计算程序ＤＰＭ的向量化模拟，大量对比计算显示向量化 后的ＤＰＭ剂量计算速度提高了约１．５倍，计算偏差小于Ｏ．６％。本文同时给出 了基于向量化模拟的蒙特卡罗剂量计算快速算法模型，该模型以Ｋａｗｒａｋｏｗ等人提出的电子历史记录复用方法为基础并进行了修正，改进了ＤＰＭ的ＧＳ模型东南大学博七学位论文快速算法和随机链式算法，结合虚拟截面、光了分裂和俄罗斯轮擞赌等技巧， 实现了多个电子输运模拟的高度欠量化。在ＩＣＣＲ基准测试计算中，新的向量化模型取得了更好的结果。关键词：放射治疗、治疗计划系统、剂晟计算、蒙特卡罗模拟、电子输运、 光子输运、方差降低技巧、快速算法、多源模型、向量化模拟、ＳＳＥ技术。ＩｌＡｂｓｔｒａｃｔＡＢＳＴＲＡＣＴＲａｄｉｏｔｈｅｒａｐｙ ｉｓ ｃａｌｃｕｌａｔｉｏｎ ｐｌａｙｓａｌｌｏｎｅｏｆ ｔｈｅ ｍｏｓｔ ｉｍｐｏｒｔａｎｔ ｍｅｔｈｏｄｓ ｆｏｒ ｔｕｍｏｒ ｔｒｅａｔｍｅｎｔ．Ｄｏｓｅａｉｍｐｏｒｔａｎｔ ｒｏｌｅ ｉｎｒａｄｉｏｔｈｅｒａｐｙ ｔｒｅａｔｍｅｎｔ ｐｌａｎｎｉｎｇ（ＲＴＰ）ｓｙｓｔｅｍ．Ｔｈｅ ｐｒｅｃｉｓｉｏｎ ｏｆ ｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎ ｈａｓ ｇｒｅａｔ ｉｎｆｌｕｅｎｃｅＥｖｅｎ ｎｏｗ，ｔｈｅ ｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎ ｍｏｄｅｌｓ ｕｓｅｄ ｉｎ ｐｒａｃｔｉｃｅａｒｃ ｎｏｔｕｐｏｎｒａｄｉｏｔｈｅｒａｐｙ．ｐｒｅｃｉｓｅ ｅｎｏｕｇｈ ＴｈｅＬｏｔｓ ａｌｌＭｏｎｔｅ Ｃａｒｌｏ（ＭＣ）ｍｅｔｈｏｄ ｉｓ ｇｏｏｄｉｎｖｅｓｔｉｇａｔｉｏｎｓ ｓｕｇｇｅｓｔｅｄ ｔｈａｔ ｉｔ ｉｓａｔｓｉｍｕｌａｔｉｎｇ ｐａｒｔｉｃｌｅ ｔｒａｎｓｐｏｒｔ ｐｒｅｃｉｓｅ ｍｅｔｈｏｄ ａｍｏｎｇｏｆｔｈｅ ｍｏｓｔｄｏｓｅｃａｌｃｕｌａｔｉｏｎ ｍｅｔｈｏｄｓＴｈｅ ｐｒｉｎｃｉｐａｌ ｄｉｓａｄｖａｎｔａｇｅ ｏｆ ａｐｐｌｙｉｎｇ ｔｈｅ Ｍｏｎｔｅ Ｃａｒｌｏａｍｅｔｈｏｄｔｏｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎ ｉｓ ｃｏｎｓｕｍｉｎｇ ｔｏｏ ｍｕｃｈ ｔｉｍｅ．Ｔｈｉｓ ｔｈｅｓｉｓ ｉｎｔｅｎｄｓ ｔｏ ｇｉｖｅｏｎｓｙｓｔｅｍａｔｉｃ ｉｎｖｅｓｔｉｇａｔｉｏｎ ｃｏｍｐｏｓｅｄ ｏｆｔｈｒｅｅ ｐａｒｔｓ．Ｐａｒｔｏｎｅｔｈｅ Ｍｏｎｔｅ Ｃａｒｌｏｍｅｔｈｏｄｆｏｒ ｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎａｎｄｉｔ ｉｓｍａｉｎｌｙ ｄｅａｌｓ ｗｉｔｈ ｔｈｅ Ｍｏｎｔｅ Ｃａｒｌｏ ｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎ ｍｏｄｅｌ ｕＳｅｄ ｉｎｒａｄｉｏｔｈｅｒａｐｙ ｔｒｅａｔｍｅｎｔ ｐｌａｎｎｉｎｇ ｓｙｓｔｅｍｓ．Ｉｎ ｔｈｉｓ ｐａｒｔ，ｗｅＣａｒｌｏ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｅｌｅｃｔｒｏｎ－ｐｈｏｔｏｎ ｃｏｕｐｌｅｄｗｉｌｌｉｎ仃ｏｄｕｅｅ ｔｈｅＭｏｎｔｅｔｒａｎｓｐｏｒｔ；ｓｔｕｄｙｔｈｅ ｃｏｎｄｅｎｓｅｄ ｈｉｓｔｏｒｙｍｅｔｈｏｄ ａｎｄ ｍｕｌｔｉｐｌｅ ｓｃａｔｔｅｒｉｎｇ ｔｈｅｏｒｙ；ｄｉｓｃｕｓｓ ｔｈｅ ｄｅｓｃｒｉｐｔｉｏｎ ｍｅｔｈｏｄ ｏｆ ｂｅａｍｓｆｒｏｍ ａｃｃｅｌｅｒａｔｏｒｓ ａｎｄ ｉｎｔｒｏｄｕｃｅ ｔｈｅ ｍｕｌｔｉｐｌｅ－ｓｏｕｒｃｅｍｏｄｅｌ；ｉｎｖｅｓｔｉｇａｔｅｉｎ ｄｅｔｍｌ ｔｈｅｍｅｔｈｏｄｓ ｏｆ ｄｅｄｕｃｉｎｇ ｔｉｓｓｕｅ ｐａｒａｍｅｔｅｒｓ ｆｒｏｍ ＣＴ ｎｕｍｂｅｒｓ；ｐｒｅｓｅｎｔ ｔｈｅ ｄｅｓｃｒｉｂｉｎｇ ｂｅａｍ ｍｏｄｉｆｉｅｒｓ ｂａｓｅｄｏｎｍｅｔｈｏｄｏｆｑｕａｄｒｉｃｓｕｒｆａｃｅｓ；ａｎｄ矛Ｖｅ ｔｈｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓ ｏｆｃａｌｃｕｌａｔｉｎｇ ｐａｔｉｅｎｔ ｄｏｓｅ ｄｉｓｔｒｉｂｕｔｉｏｎｓ．Ｐａｒｔ ｔｗｏ ｍａｉｎｌｙ ｄｅａｌｓ ｗｉｔｈ ｔｈｅ ｒａｐｉｄ ａｌｇｏｒｉｔｈｍｏｆＭｏｎｔｅ Ｃａｒｌｏ ｄｏｓｅ ｃａｌｃｕｌｍｉｏｎ．ｏｎＴｈｅ ｔｏｐｉｃ ｈａｓ ｅｘｔｅｎｓｉｖｅｌｙ ｂｅｅｎ ｉｎｖｅｓｔｉｇａｔｅｄ ｂｙ ｍａｎｙ ｐｅｏｐｌｅ．Ｂａｓｅｄｔｈｅｉｒ ｅｆｆｏｒｔｓ，Ｗｅ ｓｕｍｍａｒｉｚｅ ｆｏｕｒ ｂａｓｉｃ ｌａｗｓ ｏｆ ｒａｐｉｄ ａｌｇｏｒｉｔｈｍ，ｉ．ｅ．ｓａｖｉｎｇ ｔｉｍｅ ｂｙ ｓａｃｒｉｆｉｃｉｎｇｓｐａｃｅ，ｒｅｄｕｃｉｎｇ ｖａｒｉａｎｃｅ，ｓｉｍｐｌｉｆｙｉｎｇ ｔｈｅ ｐｈｙｓｉｃｓｍｏｄｅｌａａｎｄ ｕｔｉｌｉｚｉｎｇ ｃｏｍｐｕｔｅｒ ｒａｐｉｄ ａｌｇｏｒｉｔｈｍ，ｉ．ｅ．ｆａｓｔｒｅｓｏｕｒｃｅ ａｄｅｑｕａｔｅｌｙ．Ｔｈｅｎ ｗｅ ｉｎｖｅｓｔｉｇａｔｅ ｔｈｅ ｆｏｕｒ ｐａｒｔｓ ｏｆｓｉｍｕｌａｔｉｏｎ ｏｆ ｅｌｅｃｔｒｏｎ ｔｒａｎｓｐｏｒｔ，ｆａｓｔ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｐｈｏｔｏｎｔｒａｎｓｐｏｒｔ，ｅｆｆｅｃｔｉｖｅｕｓｅｍｕｌｔｉｐｌｅ－ｓｏｕｒｃｅ ｍｏｄｅｌ，ａｎｄ ｆａｓｔ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｂｅａｍ ｍｏｄｉｆｉｅｒｓ．Ｗｅｍｅｓｈｔｏｔｒａｎｓｐｏｒｔｄｅｓｃｒｉｂｅｃｏｍｐｌｅｘ ｂｅａｍｍｏｄｉｆｉｅｒｓ；ｄｉｓｃｕｓｓａａｓｉｍｐｌｉｆｉｅｄ ｍｅｔｈｏｄ ｏｆｓｉｍｕｌａｔｉｎｇ ｔｈｅ ａｉｒ ａｒｏｕｎｄｐａｔｉｅｎｔ；ｐｒｅｓｅｎｔｆａｓｔｍｅｔｈｏｄｏｆ ｓｉｍｕｌａｔｉｎｇ ｔｈｅ ｔｒａｎｓｐｏｒｔｏｆｐｈｏｔｏｎｓ ａｎｄ ｅｌｅｃｔｒｏｎｓ ｔｈｒｏｕｇｈｂｅａｍ ｍｏｄｉｆｉｅｒｓ．Ｉｌｌ东南大学博士学位论文Ｉｎ ｐａｒｔ ｔｈｒｅｅ，ｔｈｅ ｖｅｃｔｏｒｉｚｅｄＭｏｎｔｅ Ｃａｒｌｏ ｍｅｔｈｏｄ ｉｓｏｎｅｏｆ ｔｈｅ ｋｅｙ ｉｓｓｕｅｓ Ｔｈｅｍｅｔｈｏｄ ｏｆ ａｃｈｉｅｖｉｎｇ ｖｅｃｔｏｒｉｚｅｄ Ｍｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｅｌｅｃｔｒｏｎ ｔｒａｎｓｐｏｒｔ ｂａｓｅｄｏｎＳＳＥ ｔｅｃｈｎｏｌｏｇｙ ｉｓ ｐｒｅｓｅｎｔｅｄ ａｎｄ ｉｓ ｅｘｔｅｎｄｅｄ ｔｏ ｐｈｏｔｏｎｔｒａｎｓｐｏｒｔＡ ｆａｓｔＭｏｎｔｅｉｎｔｏａＣａｒｌｏｃｏｄｅ，ｄｏｓｅ＃ａｎｎｉｎｇ ｍｅｔｈｏｄ（ＤＰＭ），ｗａｓｅｘｔｅｎｓｉｖｅｌｙ ｍｏｄｉｆｉｅｄｖｅｃｔｏｒｉｚｅｄｃｏｄｅ（Ｖ－ＤＰＭ）．Ｃｏｍｐａｒａｔｉｖｅ ｓｉｍｕｌａｔｉｏｎｓ ｗｅｒｅ ｃｏｎｄｕｃｔｅｄ ｆｏｒ ｓｏｍｅｃａｓｅｓｔｙｐｉｃａｌ ｓｉｍｕｌａｔｉｏｎｅｖｅｒｙ ｏｆＣａＳｅｉｎ ｂｏｔｈ ＤＰＭ ａｎｄｒｕｎｓＶ－ＤＰＭ ｃｏｄｅｓ．１１１ｅ ｒｅｓｕｌｔｓ ｓｈｏｗ ｔｈａｔ ｉｎ ＤＰＭ ｃｏｄｅ ｗｉｔｈ ｖａｒｉａｎｃｅｔｈｅ Ｖ－ＤＰＭ ｃｏｄｅ１ ５ ｔｉｍｅｓ ｆａｓｔｅｒ ｔｈａｎ ｔｈｅａ０．６％．Ｆｕｒｔｈｅｒｍｏｒｅ．Ｗｅｏｎｐｍｓｅｎｔｆａｓｔ ｍｏｄｅｌ ｏｆ Ｍｏｎｔｅ Ｃａｒｌｏ ｄｏｓｅ ｃａｌｃｕｌａｔｉｏｎｏｎｂａｓｅｄｖｅｃｔｏｒｉｚｅｄＭｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｉＯｎ．Ｔｈｅ ｎｅｗ ｍｏｄｅｌ ｉｓ ｂａｓｅｄｔｈｅｍｏｄｉｆｉｅｄ ｅｌｅｃｔｒｏｎ ｈｉｓｔｏｒｙ ｒｅｐｅｔｉｔｉｏｎ ｍｅｔｈｏｄ ｐｒｅｓｅｎｔｅｄ ｂｙ Ｋａｗｒａｋｏｗ ｅｔ ａ１．Ｔ１１ｅ ｆａｓｔｃｏｍｐｕｔｉｎｇ ｍｅｔｈｏｄ ｏｆ ＧＳ ｍｕｌｔｉｐｌｅ ｓｃａｔｔｅｒｉｎｇ ｄｉｓｔｒｉｂｕｔｉｏｎ ａｒｅ ｒｅｆｏｒｍｅｄ ｉｎ ｔｈｅ ｎｅｗａｎｄｒａｎｄｏｍ ｈｉｎｇｅ ｍｏｄｅｌｍｏｄｅｌ Ｗｉｔｈ ｔｈｅｈｅｌｐ ｏｆ ｗｏｏｄｃｏｃｋ ｔｒａｃｉｎｇ，ｐｈｏｔｏｎ ｓｐｌｉｔｔｉｎｇｔｅｃｈｎｉｑｕｅ ａｎｄ Ｒｕｓｓｉａｎ Ｒｏｕｌｅ廿ｅ ｔｅｃｈｎｉｑｕｅ，ｔｈｅ ｎｅｗ ｍｏｄｅｌ ｉｎｃｒｅａｓｅ ｔｈｅ ｐｒｏｐｏｒｔｉｏｎ ｏｆｖｅｃｔｏｒｉｚｅｄ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｅｌｅｃｔｒｏｎｓ ｔｒａｎｓｐｏｒｔ．Ｏｎ ｔｈｅ ＩＣＣＲ ｂｅｎｃｈｍａｒｋ ｃａｌｃｕｌａｔｉｏｎｓ，ｇｏｏｄ ｒｅｓｕｌｔｓ ａｒｅ ｏｂｔａｉｎｅｄ ｂｙ ｔｈｅＤＰＭ ｃｏｄｅ ｗｉｔｈ ｔｈｅ ｎｅｗ ｍｏｄｅｌＫｅｙＷｏｒｄｓ：Ｒａｄｉｏｔｈｅｒａｐｙ，ＴｒｅａｔｍｅｎｔＰｌａｎｎｉｎｇ Ｓｙｓｔｅｍ，Ｄｏｓｅ Ｃａｌｃｕｌａｔｉｏｎ，Ｍｏｎｔｅ Ｃａｒｌｏ Ｓｉｍｕｌａｔｉｏｎ，Ｅｌｅｃｔｒｏｎ Ｔｒａｎｓｐｏｒｔ，Ｐｈｏｔｏｎ Ｔｒａｎｓｐｏｒｔ，Ｖａｒｉａｎｃｅ ＲｅｄｕｃｔｉｏｎＴｅｃｈｎｉｑｕｅ，Ｒａｐｉｄ Ｔｅｃｈｎｏｌｏｇｙ．Ａｌｇｏｒｉｔｈｍ，Ｍｕｌｔｉｐｌｅ―Ｓｏｕｒｃｅ Ｍｏｄｅｌ，Ｖｅｃｔｏｒｉｚｅｄ Ｓｉｍｕｌａｔｉｏｎ，ＳＳＥ东南大学学位论文独创性声明水人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得¨７川允成求。尽我所知，除了文中特别加以标注和致谢的地方外，论文。ｐ不包含儿他人已经发表或撰写过的研究成果，也不包含为获得东南大学或其它教育机构。ｊ０’引讧或证书而使用过的材料。与我‘同工作的同志对本研究所做的任何贡献均Ｌ加沦文。ｐ作了明确的说明并表示了谢意。研究生签名阿期：２Ｑ塑生５旦东南大学学位论文使用授权声明柏何大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位ｅ１＿＝＝圳ｑ复印件和电ｆ文档，可以采用影印、缩印或其他复制手段保存论文。本人ｍ』ｉ艾档的内容和纸质沦文的内容相一致。除在保密期内的保密论文外，允许论 史破禽阅和借阅，可以公布（包括刊登）论文的全部或部分内容。论文的公布（包 川ＨｆＨ＃）授权东南大学研究生院办理。研究生签名导师签名Ｈ 期第一章绪论第一章绪论§１．１课题研究背景和意义恶性肿瘤是人类五大死因（心脑血管疾病、恶性肿瘤、呼吸器官疾病、车 祸和自杀）中的第二大死因，治愈肿瘤一直是人类追求的目标。肿瘤治疗的主要疗法有外科手术、放射治疗、化学治疗等。放射治疗主要用于治疗恶性肿瘤，据统计大约７０％的恶性肿瘤患者会接受放射治疗㈣ｆ¨引。放射治疗使用历史不过一百年时问，但发展迅速。随着计算机的使用，放射治疗开始向精确化、智 能化的方向大步前进。 放射治疗利用电离辐射对人体生物组织细胞的破坏作用进行治疗。生物细 胞经过一定剂量的射线照射后可发生一系列放射生物学变化，通过分子的电离 和激发，引起生物大分子损伤。肿瘤细胞通常比正常细胞更容易受到辐射的伤 害，临床上通过选择合适的照射剂量和照射方式可实现选择性地杀死肿瘤细胞 而保护周围的正常组织，这也正是放射治疗计划系统的主要功能。但是约有近 半数患者在肿瘤没有扩散的情况下经过放射治疗后仍然死亡，从理论上说，他 们是可以被治愈的。放射治疗采用的照射剂量不准确是死亡的一个主要原因。 如没有给肿瘤细胞足以致死的剂量，将造成肿瘤的复发和转移；如使重要器官 受到了它难以承受的辐射剂量，也将引起严重的并发症或直接置人于死地。根 据Ｂｏｙｅｒ和Ｓｅｈｕｌｔｈｅｉｓｓ两人的研究结果［４ＩＪ，照射粼量的准确性提高ｌ％，治愈率可提高２％。定义１．１剂量计算：就是计算光子或带电粒子射线照射到人体（或体模）中能量沉积的空间分布。 照射剂量的误差和放射治疗的全过程都有关系，大体可分病情诊断、计划 制定和治疗实施三个方面１７６１。在过去的几十年里，有关病情诊断和放射治疗实 施的技术发展很快，ＣＴ和ＭＲＩ的出现实现了对病变组织的精确定位，随着医 用加速器的进步，各种新的治疗方式不断涌现，放射治疗的精度正在逐步提高。 相对而言，放射治疗计划系统的发展较缓慢，其核心部分剂量计算的精度仍然 难以令人满意，剂量计算的误差问题已成为制约放射治疗疗效的一个瓶颈。东南大学博士学位论文§１．１．１剂量计算要求剂最计算是放射治疗计划系统的核心内容，给出计划设计的剂量分布是放 射治疗计划系统晟重要的功能之一。剂鼍分布是评价计划是否合理的主要标准，是医生决策的关键，因此剂量分布的精度显得尤为重要。剂量分布的精度主要取决于计划系统中采用的剂量计算模型，一个成功的计划系统在剂量计算方面应满足如下要求：１、精度关于肿瘤照射剂量的误差问题，存在各种不同的标准１４１１：５％（ＩＣＲＵ１９７６）、３．５％（Ｍｉｊｎｈｅｅｒ等１９８７）、３％（Ｂｒａｈｍｅ １９８８）。放射治疗是一个相当复杂的过程，有多种因素都会影响照射剂量的精度，具体包括：处方剂量的误差、肿瘤 位置和形状确定误差、治疗机器校准精度、治疗中病人体位移动、病人呼吸影 响、器官运动影响、照射野偏移等。要达到３％的照射剂量精度在目前条件下 还不大可能，一般把ＩＣＲＵ ２４号报告推荐的５％作为照射剂量误差的允许值 １］Ｏｊ［４”。要达到这个标准，计划系统中剂量计算的误差必须在３％以内‘３５１１４”。ＩＣＲＵ ４２号报告则要求计算的剂量值和测量值的偏差须在２％以内或者半影、组织界 面处等剂量线位置偏移必须在２ｍｍ以内婶９１。２、时间利用计划系统制定治疗计划的过程一般包括计划制定、交互式调整、精计 算或优化计算剂量分布等几部分，不同阶段对剂量计算的耗时要求是不一样的： （１）交互式剂量计算：在制定计划过程中，往往需要实时计算某些点的剂量，计算时间应在１分钟以内。（２）体积剂量计算：在无需人为干预的情况下，精确计算某一体积内所有点的剂量分布时间应在１小时以内。（３）剂量优化计算：由于往往需要进行５０到１００次的迭代剂量计算， 每次剂量计算的时间也应在１分钟以内。 概括起来，一个具备临床实用性的剂量计算模型应在１分钟以内完成单野、 低精度的剂量计算；在１小时内完成多野、高精度或优化剂量计算。●２第摩绪论３、真正三维真正三维剂量计算模型，必须在剂量计算过程中考虑Ｆ述影响因了－【●患者体外轮廓的；维形状；●三维电子密度（由ＣＴ值转换）及其对原射线的影响； ●射野或放射源的ｊ维位置和彤状： ●射野的二维扩散度： ●射野的二维平坦度、对称性； ●楔形板、挡块、补偿器等线束修止装置的二维散射的影响； ●不均匀组织的二维散射影响； 只有将这些因素都考虑进去，剂量计算的精度才能得到保证。§１．１．２当前算法分析基于汁算机的放射治疗计划系统（ＲａｄｉｏｔｈｅｒａｐｙＴｒｅａｃｍｅｎ ｃ ＰｌａｎｎｉｎｇＳｙｓｔｅｍ，ＲＴＰＳ）自７０年代初出现以来发展迅速，已基本从二维发展到了兰维。 作为其核心内容的剂量计算部分也随之不断发展，已基本完成ｒ从早期的基于修正（Ｃｏｒｒｅｃｔｉｏｎ―ｂａｓｅｄ）的剂量计算模式向基于模型（Ｍｏｄｅｌ－ｂａｓｅｄ）的剂量计算模式的转变ｏ”】。后一种模式一般以建立新的剂量计算模型为着眼点，全面考 虑放射线与人体的相互作刚，从而使基于理论计算的剂量分布与基于实验测量 的剂量分布更吻合。 １、光子剂量计算模型Ａｂｎｅｒ６等人的综述对放射治疗计划系统中光子剂量计算模型进行了完整的描述…１。根据对次级电予在介质中输运和能量沉积的处理方式可将光子剂量 计算模型分成两大类…】： （Ｉ）能量局部沉积类算法：认为光子与物质作用产生的次级电子不离开作用点位置，将全部能＿星＝沉积在作用点。属于此类的方法有：线性 衰减、ＲＴＡＲ（ＥＰＬ、有效源皮距、等剂量曲线平移等）、Ｂａｔｈｏ幂指数法、ＥＴＡＲ、ＤＳＡＲ、６体积法。 （２）能量非局部沉积类算法：认为次级电子有一定能量，具有～定的射 程，可将能量传递给作用点以外的组织。属于此类的算法有各种基东南大学博士学位论文于核的卷积／叠加（Ｃｏｎｖｏｌｕｔｉｏｎ／ｓｕｐｅｒｐｏｓｉｔｉｏｎ，ｃ／ｓ）模型的方法，如ＤＳＡ、变形ＤＳＡＲ、微分笔形束、筒串卷积、笔形束卷积等；本文研究的蒙特卡罗模拟法也属于这一类。 在原射线光子与散射线光子的能量很低时，能量局部沉积类算法的假设是 可以成立的；但当射线能量较高时，次级电子的射程可长达数厘米，这个假设显然将引起误差，特别是在非电子平衡区如剂量建成区和不同组织界丽处的误 差将更大。所以计划系统更倾向于采用能量非局部沉积类算法，其中基于Ｃ／Ｓ模型的算法已是当今的主流剂量计算方法。 在Ｃ／Ｓ模型中剂量计算是通过将输入的能量注量分布与能量沉积核相卷积 或叠加实现的。能量沉积核有两种类型：点核（Ｐｏｉｍ Ｋｅｒｎｅｌ）与笔形柬核（ＰｅｎｃｉｌＢｅａｍ Ｋｅｒｎｅｌ）。●点核：如图１．１所示，点核是描述无限大介质中一原射线作用点周围的能量 沉积模式，Ｍａｃｋｉｅ等将其称为剂量扩展阵列（Ｄｏｓｅ其称为微分笔形束（Ｄｉｆｆｅｒｅｎｔｉａｌ 数（ＰｏｉｎｔＳｐｒｅａｄ Ｆｕｎｃｔｉｏｎ）【４”。 Ｐｅｎｃｉｌ Ｓｐｒｅａｄ Ａｒｒａｙ），Ｍｏｈａｎ等将Ｂｅａｍ），而ＡｈｎｅｓｊＯ等将其称为点扩散函利用点核计算体模内点ｒ处剂量的公式为：Ｄ（ｒ）＝￡ｆｆｃ％（ｒＩ）＾（Ｅ，ｒ－ｒ＇）ｄ３ｒ＇ｄＥ（１．１）式中ｒ是比释总能Ｔｅｒｍａ，即原射线与物质相互作 用在作用点处原射线释放给单位质量介质的总能 量，ｈ为预先计算好的点核。在均匀介质中，对于 一定的能量，点核ｈ的分布与空间作用点无关，核 在空间的取值只依赖于原射线作用点，，与剂量沉积 点，间的相对几何位置，而与它们在空间的绝对位 置无关。此时式１．１中的三维卷积可用快速傅立叶≤夕图１ １点核变换（ＦＦＴ）来完成。在非均匀介质中，散射线和次级电子的输运都会受到非 均匀组织的影响，点核可能随位置的变化而变化，不再是空间不变量。此时可 用下面的叠加方法来计算剂量分布：４第一章绪论Ｄ（，）＝￡Ⅲ驰ｆ）祭２（一ｍ瞰ＶⅢ～ｆ）】办’把（１２）式中ｃ（，，ｒ’）为非均匀组织密度系数，，）（ｒ）是，处密度，ｈ是密度为ｐｎ的均匀介质中的点核。式１．２中的点核不再只依赖于两点间的相对几何位置，而且还依赖 于两点间的组织密度。 利用点核进行剂量分布计算的运算量很大。在均匀组织中用式１．１计算Ｎ３ 个点的剂量需要Ｎ６次乘法和加法运算，如在不均匀组织中用式１．２计算则需要 Ｎ７次操作。运算量大势必导致计算时间过长，限制了该方法的临床应用。为此 有些学者提出了快速算法，比如：Ｂｏｙｅｒ和Ｚｈｕ等提出了沿用空间不变点核计 算非均匀组织剂量分布的方法【５７】【１４９１，这样可用快速傅立叶变换进行离散卷积。 ＡｈｅｓｊＯ提出了筒串（Ｃｏｌｌａｐｓｅｄ Ｃｏｎｅ）卷积技术ｍＪ，将散射剂量的积分划分成 具有一定立体角的共轴锥筒串，由位于简轴上散射体积单元释放的、并落入具 有相同立体角、位于筒轴上筒串中的能量，被筒轴上的体积单元线性传输、衰 减和吸收，从而将式１．２中积分变成二维筒串积分在４蔸范围内的叠加。 ●笔形束核：如图１．２所示笔形束核描述的是一束纤细射线入射到无限大均匀 介质中的能量扩散模式，笔形束核可以看成对点核在入射方向上预先进行卷积。 笔形束核也可用蒙特卡罗法模拟或用实验测量推导得到【３３］。利用笔形束核计算 体模内某点，处剂量的公式为以ｒ）＝Ｌ旺垂。ｏ）Ｈ（Ｅ，ｒ，ｓ）ｄＥｄ２ｓ（１｜３）式中中是入射粒子通量，ｎ是预先计算的笔形束核。在不同入射点ｓ和不同的ｒ处笔形束核应该是不同 的。为了计算方便，式１．３常用空间不变的笔形束核进行计算。基于笔形束核的模型可以很好地进行 非规则野剂量计算ｍ】，但是对于人体组织的不均匀 性和体表不规则轮廓的处理仍须采用基于宽束的校图１．２笔形束核正措施，精度不高。式１．３是两维卷积，运算量大为减少。 速度快是笔形柬模 型的突出特点，已广泛用于优化算法中。东南大学博士学讯论文２、电子剂量计算模型 常用电子剂量计算模型有：经验模犁、陈化扩散方程模型、“原射”和散射分量分别计算模型、笔形束模型、双群模型、蒙特卡罗模拟。前三种模型的计算精度较低，无法很好处理介质的不均匀性。基于Ｆｅｒｍｉ．Ｅｙｇｅｓ理论的笔形束 模型在计算均匀和非均匀组织的剂量分布时都比较成功‘”Ｉ，并已为多数治疗计划系统所采用。但由于Ｆｅｒｍｉ―Ｅｙｇｅｓ理论是小角度多级散射理论．此模型还存在不少问题，比如：（１）全然没有考虑原射电子的反向散射问题；（２）忽略了骨．软组织界面处由骨组织反射回来的电子对剂量的贡献；（３）对电子斜入射校 正尚不完善；（４）由于电子剂量分布受电子准直器的影响更大，不规则野与标准野的数据转换还存在问题；（５）高能次级电子的处理在不均匀组织中误差较 犬。这些问题导致计算非均匀介质剂量分布时准确性较差，当治疗体内包含空 气腔或骨头时，该算法的结果和测量值或蒙特卡罗模拟结果相差可达１０％以上１５９１［ｓｑＩｌ０６１ １１１０】【１３０１。带电粒子输运双群（Ｂｉｐａｒｔｉｔｉｏｎ）模型由我国罗正明教授首次提出【１ ８Ｊ【１００］。 该模型将带电粒子在介质中的输运按角分布特征分为两群：向前运动的小角电 子和向四方运动的扩散电子。利用从遭受散射的向前运动电子中扣出大角散射电子的方法，得到一组计算电子能量损失分布的方程。该模型计算速度快，在 均匀介质中计算结果和蒙特卡罗方法吻合得很好【…１，但该模型在计划系统中应用不多。综合以上分析可知，目前计划系统中剂量计算模型在计算均匀介质剂量分 布时，精度和速度都可满足临床剂量计算要求‘１３ ７１。对于人体这样的非均匀介质， 在速度上满足临床要求的计算模型在精度上却尚未达到临床要求，特别是不同 组织界面处的计算精度较低。§１．２蒙特卡罗剂量计算方法放射线照射到人体后，光子和电子在人体组织中传输并不断沉积能量，剂 量计算求的正是沉积能量的分布，实际上这是一个粒子辐射输运问题。蒙特昔罗（Ｍｏｎｔｅ Ｃａｒｌｏ，ＭＣ）方法擅长求解粒子输运问题１６ｌ【７Ｊｆｌ２ｉ１２０】｜２９１，在很多复杂条件下它的求解精度是最高的，放射治疗剂量计算就是一例。６第‘章绪论放射治疗过程中商能光子和电子与物质之间的相互作用是一随机过程，剂 量计算所求物理量剂量分布是一个统计值，人体组织为结构复杂的非均匀介质，在此情况卜‘进行解析计算难度相当大，各种近似和大量的简化处理是必不可少 的，所得解的精度自然无法保证。蒙特卡罗方法是一种随机算法，通过对粒予与物质相互作用进行随机模拟 来获得粒子在人体组织中沉积能量的分布。如果完全模拟粒子发生的每一次作 用，那么蒙特膏罗模拟的精度将取决于它所采用的截而数据，只要这些截面数据的精度得到保证，那么蒙特＃罗模拟的精度也就得到保证。有关蒙特卡罗模拟光子和电子输运的理论研究已达到相当高度，在不考虑计算时间的前提下， 众多通用蒙特卡罗模拟程序具备＿『进行剂量计算的能力，由于它们采用的近似 和简化处理很少，因此蒙特卡罗剂量计算方法是当前所有剂量计算方法中最精 确的，在均匀和非均匀介质中都能满足临床剂量计算精度的要求。但是蒙特卡 罗剂量计算方法也有自身的不足，收敛速度慢、计算时间长是它的主要弱点， 故长期以来一直无法为临床放射治疗所接受，只能作为其它算法的验证基准或 是用来计算能量沉积核。如何在保持高精度特点的前提下，加快计算速度是蒙特＃罗剂量计算方法面临的主要课题。§１．３国内外的研究情况国外自：十世纪八十年代开始研究蒙特卡罗剂量计算方法，在当时条件Ｆ 蒙特卡罗方法计算时闯长得无法接受，故当时蒙特卡罗方法只能间接地应用于 剂量计算，比如计算能量沉积核、作为其它算法精度的验证工具等【ｌ ０１１。进入九 卜年代后，随着计算机速度的提高，蒙特卡罗方法本身也得到不断改进和发展。 目前，国际｝＝：蒙特一｝罗剂量计算研究水平较高的国家主要有美国、加拿大、德 国、英国等。蒙特卡罗剂量计算的研究工作主要集中在两个大方面，即如何保 证蒙特忙罗方法的高精度和如何提高蒙特卡罗方法的计算效率，具体内容为： （１）通过模拟粒子在直线加速器中的传输过程，得到加速器出射粒子束 的相空间信息：其结果不仅可用于蒙特卡罗剂量计算，同样可用于 现有剂量计算算法，提高它们的计算精度； （２）为避免直接使用辐射源相空间数据带来的效率损失，通过建立多源东南大学博士学位论文模型来重建入射粒子的相空间信息： （３）根据病人的ＣＴ扫描数据推导人体组织参数（物理密度和化学组成）：（４）可提高蒙特卡罗剂量计算速度的新方法和新技巧； （５）蒙特号罗方法在适形放射治疗、立体定向放射治疗、逆向计划系统 等方面的应用。在此期间，国际Ｌ陆续出现了一批蒙特卡罗剂量计算的快速算法和代码， 它们大多以通用蒙特卡罗模拟程序ＥＧＳ４１”５１的算法为基础，比较重要的有ＳＭＣｌ９２１１９”、ＭＭＣｔｌｌ们［１１７ｌ、ＶＭＣｔ“１１１１１】【１”１、ＤＰＭｌｌ２ ８１、ＭＣＤＯＳＥ［９８１等。这些蒙特卡罗剂量计算程序在速度上相对于通用模拟程序（如ＥＧＳ４、ＭＣＮＰ等）有 了较大提高，电子剂量计算效率增加约１０一３０倍，光子剂量计算效率增加约 ２０一１００倍。如§１．１．１所述临床上对剂量计算的速度要求是很高的，这些快速 算法虽然目前还无法完全满足要求，但已经向前跨出了一大步。很多研究者都 相信，在不远的将来蒙特卡罗方法将成为临床剂量计算的主要方法 １４６１［４７㈣】【１０６】【１０７１１１”Ｉ。国内进行蒙特旨罗剂量计算方法的研究始于九卜年代初，相关研究论文较少ｆＪｆ【Ⅻ１７Ｊｆ∞１，进展较慢，尚未发现关于蒙特卡罗剂量计算快速算法的研究，与国外研究水平差距较大。§Ｉ．４本文工作本文二ｒ：作可分为两个部分，一是研究适用于放射治疗计划系统的蒙特。号罗剂量计算模型，本文在前人工作的基础上进行了改进和补充；二是研究蒙特卡罗剂量计算的快速算法，重点研究了向量化蒙特卡罗模拟法，并建立了基于向 量化模拟的蒙特卡罗剂量计算快速算法模型。 本文主要研究光子和电子的蒙特卡罗剂量计算方法，其它使用范围较小的辐射源如中子和质子等剂量计算不考虑。本文具体内容安排如下：第一章为绪论部分，分析了各种剂量计算方法的优缺点，介绍了蒙特卡罗剂量计算方法的研究背景、现状及本文的目标。第一章绪论第二章的研究内容主要为蒙特砖罗方法的基本介绍及粒子输运过程的蒙特忙罗模拟。重点在于光了一电于耦合输运过程的模拟；详细讨论了浓缩历史方法和多次散射理论；给出ｒ糟二Ｆ适用于粒子输运的方差降低技巧；最后简要介 绍ｒ一些通用蒙特卡罗模拟程序。 第二三章捕述ｒ适用于放射治疗计划系统的蒙特｝罗剂量计算模型。主要研究了基于直线加速器的辐射源的蒙特卡罗模拟以及多源模型的使用；讨论了根据ＣＴ值推导人体组织参数的多种方法及存在问题；介绍了基于：次曲面描述 射线束修正设备的方法；阐述了人体内部蒙特卡罗剂量分布计算的特点。 第四章主要研究了蒙特＃罗剂量计算的快速算法。总结出提高蒙特卡罗剂量计算速度的四个基本原则；研究了提高多源模型抽样效率的方法；介绍了光 子和电子通过刺线束修正设备的快速模拟方法；给出了人体周围空气的简化模 拟方法：详细研究了电子和光子输运的快速模拟方法。 第五章主要研究了向量化蒙特卡罗模拟方法。提出了一种基于ＳＳＥ技术的 电子输运向量化模拟算法；研究ｒ该算法对于现有快速算法模拟效率的提升； 给出了基于向量化模拟的蒙特卡罗剂量计算快速模型，并对新模型进行了验证。 第六章对本文工作进行了总结，并展望了今后的研究方向。第二章蒙特卡罗方法及对粒子输运的模拟第二章蒙特卡罗方法及对粒子输运的模拟§２．１蒙特卡罗方法基础蒙特卡罗方法是一种应用广泛的计算方法，用蒙特｝罗这一世界闻名的赌城来命名反映了它的本质一一随机性和统计性，又称随机模拟（Ｒａｎｄｏｍｓｉｍｕｌａｔｉｏｎ）方法、随机抽样（Ｒａｎｄｏｍ ｓａｍｐｌｉｎｇ）技巧或统计实验（Ｓｔａｔｉｓｔｉｃａｌｔｅｓｔｉｎｇ）方法。二卜世纪四十年代中期Ｖｏｎ Ｎｅｕｍａｎｎ等人将其作为一种独立的方法提出，并在核武器的研制中得到了应用Ⅲ２０】【３ｌ】。蒙特卡罗方法的基本思想 并不新颖，可简述为：当所求问题的解是某事件出现的概率，或是某个随机变 量的数学期望时，可以通过某种试验方法，得到该事件出现的频率，或者该随 机变量若下个具体观察值的算术平均值，并用它们作为问题的解。这一点人们 很早以前就已发现并加以利用，如十八世纪法国学者蒲丰（Ｂｕｆｆｏｎ）采用随机 投针试验来计算兀值。 蒙特卡罗方法是独具风格的数值计算方法，它首先建立求解问题的概率模 型，使模型的参数等于问题的解，然后通过对模型的抽样试验来计算参数的统 计特征，最后给出所求解的近似值。对那些用纯数学方法处理太复杂或太困难 的问题，蒙特卡罗方法往往是最简便易行的近似计算方法。蒙特卡罗方法的理 论基础是概率统计理论，但它并不进行真实的实验，而是根据问题特征，用数 学方法对实际过程加以模拟。比如，蒙特卡罗方法计算积分，将所要计算的积 分看作服从某种分布密度函数坂力的随机变量烈，）的数学期望＜ｇ＞＝括（，）／（，．）毋通过试验，得到Ⅳ个观察值，－，，２（２．１）ｒＮ，将相应的Ⅳ个随机变量的值甙ｒ１），ｇ（，２）…烈憎）的算术平均值ｇ”２专善ｇ（‘）＠‘２’作为积分的估计值（近似值）。为获得一定精确度，Ⅳ值必须很大，通过人工方 法作如此大量的试验相对困难，甚至是不可能，计算机的出现使得人们可以把东南大学博士学位论史数目巨大的随机试验交给计算机完成，用计算机模拟替代复杂昂贵或无法进行 的实验，这使得蒙特矗罗方法发展成为一种应用面很广的重要的数值模拟方法， 图２．１姓示ｒ蒙特母岁方法解决实际问题的基本原理。蒙特卡罗方法的应用范 陶包括：粒子输运问题、统计物理、核物理、典型数学问题、交通流量控制、 医学、生物、矿产探测等，粒子（中子、光子和带电粒子）输运问题中，蒙特卡罗方法更是成为最主要的数值计算手段㈦【１３】【１５】【２５】１２７】【２９］１３８心ｊ。物理实验问题Ｉ鼍粼“Ｉ实际问题的近似解 具体物理量的估计物理规律数学定理’ｆ大数定理随机变量的数学 期望和方差估计ｔｌ中心极限定＿＿獗孚物。矗覆垂一－Ｉ｛ｌ概率统计方法｜随机数学屠１ｌｌ概率统计模型ｌ夸ｉ鬻慕机－…。．．．．―ｉｉ．一；：＝：＿一 随机变量的抽样，样本分布图２１蒙特譬罗方法原理蒙特卡罗方法的优缺点概括如下ｐ”，优点：能够比较逼真地描述具有随机性质的事物特点及物理实验过程，受几何条件限制小，收敛速度与问题的维数 无关，具有同时计算多个方案与多个未知量的能力，误差容易确定，程序结构清晰，分块性强：缺点：收敛速度慢，误差具有概率性，在处理粒子输运问题时，计算结果与系统大小有关。 蒙特号罗方法解决问题可分为三个步骤：构造或描述概率过程；从已知概 率分布抽样；建立各种估计量。１、构造或描述概率过程：对于本身就具有随机性质的问题，主要是正确 描述和模拟这个概率过程：对于不是随机性质的确定性问题，必须事先构造一 个人为的概率过程，它的某些参量正好是所求问题的解。概率模型的复杂程度决定ｒ蒙特卡罗计算的复杂程度。 ２、从已知分布抽样：即由已知分布的总体中抽取简单子样，是使用严格第二章蒙特卡罗方法及对粒子输运的模拟的数学方法，借助于随机数产生的。只要随机数序列满足均匀性及相互独寺的要求，由它产生的任何分布的简单＿『样严格满足具有相同的总体分布且相互独立。３、建立各种估计量：建立各种估计景，相当于对模拟试验的结果进行考 察和登记，从中得到问题的解。构造了概率模型并能从中抽样后，需要确定一个随机变量，作为所求问题解的估计量。如果这个随机变量的期望值正好是所 求问题的解，称为无偏估计，蒙特＃罗方法中使用最多的就足无偏估计１２”ｊ。§２．１．１收敛性和误差估计蒙特卡罗方法常以随机变量ｚ的简单子样ｚ１，Ｚ２…ＺＮ的算术平均值乙＝专善ｚ布，且期望值有限，则其均值以概率１收敛到ｚ的期望值，即（２．３）作为所求解的近似值。根据柯尔莫哥洛夫强大数定律，如ＺＩ，ｚ２…ＺＮ独立同分以ＺⅣ―币：著＿Ｅ（ｚ））＝１根据中心极限定理，当Ｎ充分大时（２ ４）Ｐ（ＩｚＮ一翰Ｉ＜丽ｘｏ＂）＊面１驴ｘ出＝１一“（２．５）式中ａ是ｚ的均方差，仅是置信度，１一ｎ为置信水平。在１一盘置信水平下，近似值与真值的误差为二等，根据问题的要求确定出置信水平后，查标准Ⅱ！态 ＿ｈ分布表，可确定工，得到误差。显然这里所说的误差为概率误差，均方差一是 未知量，实际使用的是其估计值盯。氍１Ｎ丽２１Ｎ２（２ ６）显然，一有限且不为零时，蒙特卡罗方法的收敛速度为Ｏ（Ｎ ２）。东南大学博士学位论文§２．１．２随机抽样方法由已知分布的随机抽样指的是由已知分布的总体中抽取简单子样，抽样的 ｆｉ｛ｆ提是获得均匀同相互独立的随机数序列，由此产生的简单了样将严格具有相 同总体分布且相互独立。衡最一种随机抽样方法的标准是其在计算机上的运行 效率，而不管它的实现如何复杂。多变量的随机抽样可以归结为单变量的随机 抽样，所以这里我们只考虑单变量的随机抽样方法。１．随机数发生器 在（０，１）范围内均匀分布的随机数是进行随机抽样的基础，通常计算机上生成的随机数序列是根据一个确定算法计算出来，并非真正的随机数，称“伪随机数”更合适。但只要算法选择得当，使用伪随机数对于蒙特卡罗方法的计 算精度几乎没有影响。一个高质量的伪随机数发生器应具有如下特性： （１）生成的伪随机数序列的相关性越小越好。 （２）生成的伪随机数序列的均匀性要好，即同一范围内的任何一个随机数 出现的概率相同。 （３）生成的伪随机数序列的周期越长越好。（４）算法的计算效率越高越好。 目前常用的伪随机数发生器有：线性同余发生器、多重递归发生器、斐波那契发生器、移位寄存发生器、反转同余发生器和复合线性同余发生器等。其中线性乘同余发生器最为简单，如Ｒ＝７５Ｒ。（ｍｏｄ２”一１），￡＝Ｒ／（２”一１）（２ ７）给一个种子Ｒ（＜２”一１），上述公式将产生一系列在０和ｌ之间均匀分布的随机数，序列的周期是１０９，这个周期对于蒙特卡罗方法来说是比较小的。Ｊａｍｅｓ的 方法可产生周期达１０１８的（Ｏ，１）之间均匀分布随机数序列【７＝Ｉ】，对实际蒙特卡 罗方法运算而言，此序列可以当成是无限长序列，本文随后进行的所有蒙特卡罗运算采用的都是该伪随机序列。１４第二二章蒙特乍岁方法及对粒子输运的模拟２、直接抽样法 连续分布的随机变量Ｘ的分布函数，’（Ｊ）＝ｒｐ（一）ｄｘ’定义一个新随机变量喜＝尸０），其概率密度函数为：（２ ８）是相对于Ｊ的单调递增函数，设Ｊ取值范围为［ａ，ｂ］，则ｐ（ａ）＝Ｏ，Ｐ（６）＝ｌ。心）叫工）ｆ譬卜１ Ｌ出／（２ ９）显然，ｆ在（ｏ，１）范围内均匀分布。如图２．２所示，如果＜是一个随机数，则ｘ＝Ｐ。（国在（ａ，ｂ）中以概率密度函数ｐ０）随机分布，Ｊ的随机性由ｆ的随机性来保证。这时对Ｘ的抽样就变成对下面方程求根掌＝ｅｐ（一）ｄｘ’性极高。（２．１０）当卢Ｏ）为简单解析表达式时，式２ １０可以得到解析解，此时直接抽样法的准确图２．２利用直接抽样法从分布ｐ０）中进行随机抽样即使厦ｘ）太复杂而无法进行解析抽样或连续分布的础）是用数值形式表示的，直接抽样法也可以进行有效的抽样。这时，取样函数ｐ（ｘ）＝善可通过反向插值解出，即根据表（萎，ｔ）进行插值得到：喜＝Ｐ（ｘ，）。当然必须注意进行数东南大学博士学位论文值插值和数值积分时的误差问题，网格点不能分得太粗。 直接抽样法同样适用于离散分布Ⅳｐ（』）＝∑Ｐ，占（工一ｆ）ｊｔｌ（２．１１）抽样公式为：ｚ＝１如果善－＜ＰｌＪ＝劢口果，，Ｉ＜ｆ≤Ｐｌ＋ｐ２：（２．１２）工＝／如果∑：Ｂ＜善≤∑：。ｐ，３．舍选抽样方法 在待抽样分布函数的反函数无法给出或其反函数太复杂、运算量太大的情 况下，直接抽样法是不合适采用的。Ｖｏｎ Ｎｅｕｍａｎｎ提出的舍选（Ｒｅｊｅｃｔｉｏｎ）抽 样法可以较好地克服这些困难。舍选法的一般形式为： ｐ（膏）＝Ｏｒ（ｘ）ｒ０）（２１３）式中巩ｘ）是一个适合直接抽样法进行抽样的概率密度函数，ｃ是一个正常数，函数，０）满足条件０＜如珏１。舍选法的抽样步骤为：（１）根据分布础）抽样得到一个随机值Ｘ。（２）生成一个随机数善。（３）如果善＞如），回到（１），否则接受工。图２．３利用舍选法从分布ｍ）中进行随机抽样１６第二章蒙特卡罗方法及对粒子输运的模拟从图２．３中可以看出，前两步得到的点（。∥）在满足Ｃ万（＿）ｆ）≥ｐ（ｘ）的区域Ａ中是均匀分布的，（３）步是拒绝那些户烈ｘ）的点，（４）步是接受肭）的点，显然，接受的点在ｚ轴和曲线，印０）围成区域中均匀分布，故ｘ值服从ｐ０）分布。舍选抽样法的抽样效率为接受一个生成ｘ值的概率，即ｓ＝ｊ：；ｍ）ｍ）出＝去（２ １４）只有当ｘ（ｘ）＝ｐ（砷时，ｃ＝１。显然舍选法是用牺牲一些抽样效率为代价，用 容易抽样的ｘ（ｘ）替代了难于抽样的ｐ（ｘ）。在使用舍选法时要注意以下两个问题：选取的石似）要容易抽样， ４．复合抽样方法 复合分布是实际问题中常见的分布，即一个随机变量服从的分布与一个参 数有关，而该参数也是～个服从确定分布的随机变量。下面为一复合分布 Ｃ值尽可能小。／（ｘ）＝ｆ正（ｊｆ／ｙ）ｃ蜗（ｙ）（２ １５）式Ａ（ｘ／ｙ）表示与参数ｙ有关的条件分布密度函数；Ｆｌ∽表示分布函数。 复合分布的抽样方法为：首先由分布函数Ｆ?Ｏ，）（或分布密度函数＾∽）中抽样玢，（或期），然后再由分布密度函数五㈨啪）中抽样确定Ｘ，＝ｘ＾＠／炸１）（２．１６）§２．２粒子输运的蒙特卡罗模拟在核武器研究、反应堆设计、辐射测量、屏蔽以及放射医学等很多领域都 涉及粒子的迁移输运问题，放射治疗中粒子迁移传输也扮演了重要角色‘”１。 定义２．１输运理论：研究微观粒子（光子、电子、中子等）在介质中的迁移统计规律的数学理论。定义２．２碰撞定律：用以描述粒子与物质碰撞之后发生的情况，常由下述函数表达，‘（Ｅ‘，而’斗Ｅ，壳）刎翻晓＝在第，类核介质中，当莱种ｘ型的相互作ｊ７东南大学博士学位论文用发生后，粒子的能量由Ｅ’和运动方向由壳’变为能量在Ｅ处ｄＥ范围内，运动方向在Ｑ处ｍ范围内的几率。定义２．３截面（ｃｒｏｓ ｓ ｓｅｃｔｉｏｎ）：是描述粒子与物质相互作用概率的物理 量，定义为一个入射粒子与单位面积上一个靶粒子发生相互作用的概率。如果 一个入射粒子与物质的相互作用有多种相互独立的作用方式，则相互作用总截 面等于各种作用截面之和，截面具有面积的量纲。 定义２．４微分截面（ｄｉｆｆｅｒｅｎｔｉａｌｃｒｏｓｓｓｅｃｔｉｏｎ，ＤＣＳ）：是描述相互作用后某种粒子的角分布特征，含义为一个八射粒子与单位面积上一个靶粒子发生相互 作用，并且作用后粒子飞行在某方向上单位立体角内的概率。 定义２．５散射：当电子通过物质时，由于与原子的相互作用，不仅因电离、激发和韧致辐射而损失能量，而且运动方向也会发生很大改变，这种入射粒子在通过物质时运动方向发生改变的现象就称为散射。光子的静Ｉ七质量为零，不算是粒子型物质，但从量子力学观点看它的行为却完全象是粒子【“。输运理论研究的不是个别粒子的运动，而是研究大量粒子运动所表现的非平衡统计运动规律。高能粒子（如光子、电子和正电子等）在介质中迁移时会与介质发生大量相互作用，并把能量转移给物质中的原子和分 子生成次级粒子，通过和介质的不断作用，一个高能粒子会引发一个粒子簇射 流（Ｃａｓｃａｄｅ），常被称为一个簇射（Ｓｈｏｗｅｒ）。粒子发生一次作用以后，它的能 量要减少，更多粒子可能被生成，簇射过程可以看成有效能量减少的过程。随 着时间推移，初始能量将不断地沉积到介质中，其它能量则在新生成粒子上。 粒子输运的基本规律，可以用如Ｆ稳态玻尔滋曼输运方程加以描述１２０】亏①（尹，Ｅ，奶△＋∥（Ｅ）中（尹，Ｅ，国＝盯ｍ（芦，∥，盎’）＂仃（盎７―＇盎，Ｅ’―争Ｅ）ｄＥ＇ｄ晓’＋ｓ（ｅ，Ｅ，Ｑ）（２．１７）式中通量①（尹，Ｅ，Ｑ）的物理意义为在与Ｑ方向垂直的平面上单位时间、单位能最、单位立体角内，在单位面积上通过的粒子数目。ｓ（Ｆ，Ｅ，ｆｉ）是初始粒子的源分布：单位时间内在尹处体积元中，能量在（Ｅ，Ｅ＋积）内，方向在Ｑ附近的ｌ Ｒ第二章蒙特卡罗方法及对粒子输远的模拟班２内生成的粒子数目。单位长度．Ｉ二发生任一类型作用的概率为“（』！１１）。Ｆ处单何体积内的散射巾一ＬＬ，＞ｂ ｎ，口（盎’＿叠，￡７斗Ｆ）是发生散射（从方向盎’到方向壶， 能茕从Ｅ’到，ｒ）｝精微分截嘶。这个线性玻尔滋曼微积分方程的解将产，卜粒了通 量的完整相位空间捕述，可以进一步得到吸收剂嚣等和辐射测量相关的量。 玻尔滋曼方程的解析计算难度相当大，必须采用大量的近似处理且只能应 用于一些简单几何体。输运方程的解析解、半解析解已从一维问题逐步发展到 ＿厂精确求解二维问题。计算机诞生以后，数值求解粒子输运方程的方法得到迅 速发展和广泛应用，首先应用的是球谐函数法，随后是蒙特卡罗模拟方法，近年来有限元方法也得到了较快发展。一般数值计算方法和解析法在处理截面与 能节相关、散射各向异性、非均匀介质、复杂几何、时间相关等方面碰到了较 大困难。而对于蒙特卡罗方法，只不过表示比通常情况Ｆ，维数增多了，分布 密度变化了，几何条件复杂了，根据蒙特Ｆ罗方法的特点，这不会产生什么本 质性困难。蒙特卡罗方法比其它近似计算方法更接近于真实情况，近似较少。 如果完全采用洋细模拟（按照发生的先后顺序模拟一个粒子所经历的所有作 用），它产生了和传输方程精确解同样的结果（不考虑统计误差）。粒子在介质中输运过程具有随机性质。一个由源发出的粒子，在它的运动方向上，在哪一点碰撞是偶然的，但有一定的几率分布；与原子核（或原子）发 生碰撞，又有各种几率不尽相同的类型；碰撞后的能量和方向，又遵从一定的 概率分布。粒子可能被介质吸收或从系统中逃脱，这时粒子运动过程就告结束， 否则，继续下一次类似的运动过程。一个粒子在介质中运动的情况，可由它所经历的碰撞反映出来。下一次碰撞位置、碰撞后的能量与方向的决定，只与这次碰撞情况有关，而与粒子以前的碰撞情况无关。显然，粒子输运过程是一种马尔可夫过程。只要粒子发生碰撞的规律在物理上清楚，那么粒子输运过程完全能够用蒙特卡罗方法正确地模拟，从而得到所需要的解‘２们。蒙特卡罗模拟方 法的弱点在于误差的概率性质，收敛速度慢，特别对于分布量计算问题。 一般而言，粒子输运过程的蒙特卡罗模拟可概括为三个步骤【３ｌＪ： （１）对源分布的抽样过程：源分布抽样目的是产生粒子的初始状态Ｓ。＝（芘，Ｅ，Ｑ。），具体包括源粒子位置、能量、运动方向的抽样。东南大学博士学位论文（２）空间、能量和运动方向的随机游动过程：一个粒子由状态岛到状态品＋，（ｍ＝ｏ，１，２…）后，需确定粒’Ｆ空间’位置ｉ∥能量Ｅｍ∥和运动方向壳。。当该粒子的能茕低于传输闽值或离开系统后，终止粒子的游动过程，开 始新粒子的随机游动。（３）记录贡献与分析结果：在粒子的随机游动过程中，同时记下每个粒子对所求量的贞献，对于同一个所求量，又随所用的蒙特卡罗技巧不同而不同。当所有粒子都模拟以后，得到所求量的估计值和方差的估计值，进而得到误差。在计算机上进行粒子输运蒙特卡罗模拟可看成是以上三个过程不断重复的过 程，在进行模拟前需要先解决以Ｆ问题： （１）粒子输运全部物理过程和相应物理数据及作用模型：（２）随机数生成和随机抽样方法； （３）所求量的记录及统计误差； （４）系统的几何描述。粒子的种类有很多，本文考虑的是放射治疗过程中使用最广的两种光子和 电子，通常放射治疗可看成是光子和电子耦合输运问题。所谓耦合输运即光子在输运过程中可产生电子，而电子又可以产生光子，这两种粒子互相产生的过程称为两种粒子的耦合输运过程。§２．３光子输运的蒙特卡罗模拟§２．３．１直接模拟方法光予输运模拟可采取直接模拟法，即直接从物理问题出发，模拟粒子的真 实物理过程ｐ”。图２．４是通用蒙特卡罗模拟程序ＰＥＮＥＬＯＰＥ模拟光子输运过程的流程图。（１）计算光子到下一个作用点的距离Ｄｓ光子和介质发生的作用主要有四种：相干散射（Ｃｏｈｅｒｅｎｔ ｓｃａｔｔｅｒｉｎｇ）、非相第二章蒙特卡罗方法及对粒子输逗的模拟ｌ散射（Ｉｎｃｏｈｅｒｅｎｔ ｓｃａｔｔｅｒｉｎｇ）、光电效戍（Ｐｈｏｔｏｅｌｅｃｔｒｉｃ ｅｆｆｅｃｔ）和电子对生成（Ｅｌｅｃｔｒｏｎ．ｐｏｓｉｔｒｏｎ ｐａｉｒｐｒｏｄｕｃｔｉｏｎ）。存当前介质中每种作朋发生的截面可由各自微分截面ＤＣＳ积分得到：盯（Ｆ）＝ｆｏＥｄＷＪ０［＂２ｒｅｓｉｎｏｄａ㈣２ｃｒ（Ｆ；∥，臼）（２ １８）式中Ｗ为能量损失，０和ｐ为光子散射极角’ｊ方位角。由于光子作用后方向偏 转是沿入射方向轴对称的，故ＤＣＳ与方位角无关。光子和介质发生作用的总截 面为各种作用截面之和： 口，（Ｅ）＝盯∞（Ｆ）＋仃ｍ（Ｅ）＋盯砷（Ｆ）十盯即（Ｅ）（２ １９）等式右边分别是相干散射、康普顿散射、光电效应和电子对生成的截面。光子 运行单位长度发生作用的概率为：∑＝Ⅳ盯，（２．２０）Ⅳ＿虬毒传输ｒ距离ｒ以后还没有发生任何作用的概率为：Ｑ?２Ｄ式中Ⅳ。是阿佛伽得罗常数，Ｐ是介质的质量密度，ＡＭ是介质的分子量。光孑刖＝Ｚｅ－∑‘构造一个分布函数：（２ ２２）Ｆ（ｒ）＝肛∑ｅ一∑。＝１百∑‘求出它的反函数：（２ ２３）卜古ｌｏｇ【１－砸）】在０和１之间取一个均匀分布的数‘，光子到下一次作用的距离就为：（２?２４）忙一古ｂ卅Ｄ（２ ２５’东南大学博士学位论文如１一ｆ在０和ｌ之间均匀分布，那么ｆ在０和１之间也是均匀分布，故可取ｆ：一士１０９（掌） ｋ―Ｆ吲孝’式中ｆ即所求的Ｄｓ。广Ｌ（２ ２６）ｌ从辐射源中抽样一初始光子ｊ一■■～―ｊ一一一―――――，――．．．．，――．―――Ｊ图２ ４光子输运的蒙特卡罗模拟（２）光子输运 计算光予沿当前方向（“，ｖ，ｗ）运行到当前区域边界的距离／）ｎｅａｒ，如Ｄｓ小于Ｄｎｅａｒ，则光子沿当前方向运行距离Ｄｓ，光子到达作用点位置：ｚ’＝工＋Ｉ，／×Ｄｓ （２．２７） （２ ２８）ｙ’＝ｙ＋Ｖ×胁第二章蒙特卡罗方法及对粒子输运的模拟：’＝：＋Ｗ×胍如Ｄｓ大于Ｄｎｅａｒ，光子沿当前方向运行距离／）ｎｅａｒ。 （３）决定作用类型（２．２９）光子到达作用点后，具体发生何种作用需要随机抽样。假设９０％可能发生 光电效应，１０％发生瑞利散射，选择一随机数如小于０．９，发生光电作用，否则发生瑞利散射。接下来对所发生作用进行模拟，计算光子能最改变和偏转角度 ０和口以及是否生成新粒子等等。（４）光子越过区域边界 当Ｄｓ大于／）ｎｅａｒ时，光子将运行到当前区域边晃，此时有两种模拟方式：一种是ＰＥＮＥＬＯＰＥ采用的方式，由于粒子输运为一马尔可夫过程，允许模拟过程的任意一点终止，以后冉从此点接着模拟‘７＂。故光子进入新区域后，立即重 新计算Ｄｓ；另一种是ＥＧＳ４中采用的方式，将Ｄｓ计算公式分解成如下三式：Ｄｓ＝以ｍ（２．３０）以：三１（２．３１）∑ｍ＝～ｌｏｇ（ｆ、（２．３２）式中ｈ是光子作用问的平均自由程（ＭｅａｎＦｒｅｅＰａｔｈ，ＭＦＰ），ｍ相当于光子将运行的ＭＦＰ数日。进入新区域后，如介质类型改变，ｈ需重新计算；而ｍ和介质类型无关，只减去／）ｎｅａｒ等效的ＩｖｌＦＰ数即可用式２．３０计算出剩余的Ｄｓ。。 （５）光子方向改变光子发生作用后运行方向将发生改变，通常模拟过程将计算出如图２．５所 示的偏转极角０和方位角妒。 如光子是沿着ｚ轴方向（０，０，Ｉ）运动，偏转角度秽和伊后，新方向为：（ｓｉｎＯｃｏｓ《ｏ，ｓｉｎＯｓｉｎｆａ，ｃｏｓＯ）（２．３３）如光子是沿任意方向（“，ｖ，ｗ）入射的，由旋转矩阵Ｒ可以得出光子的方向偏转：东南大学博士学位论文Ｃ｛宝一ｓｉｎ ｄｐ妒缈ｃｏｓ?ｐＲ＝Ｃ｛拿一ｓｉＩｌｐｃｏｓ妒１ ｓｉｎｏｓｉｍｐ Ｉｃｏｓ６（２＿３４）／，． ．。Ｌ口口．一§．１｜蚶０Ｊ光予的新方向为拈…ｓ臼＋考ｈ唧堋却】ｗ ｖ’＝ｖｅ：０ｓＯ＋蔫［ｗｃｏｓｑ，＋ｕｓｉｎ妒１、，１一。（２．３５）（２．３６）、『１一‘∥ｗＩ：ｗｃｏｓ０一厢ｓｉｎｏｃｏｓ驴§２．３．２光子作用机制介质的作用机制。 １．相干散射（２．３７）下面结合通用蒙特卡罗模拟程序ＰＥＮＥＬＯＰＥ［”１【７３ｌ的模拟过程，介绍光予和相干散射又称瑞利散射（Ｒａｙｌｅｉ曲Ｓｃａｔｔｅｒｉｎｇ），用电磁辐射的波动原理可解 释如图２．６所示的相干散射过程：当电磁波从电子附件通过时，会使电子产生 振荡，振荡电子按入射电磁波的频率以光子形式发射能量。散射光子的波长与 入射束相同，没有任何能量转换成电子动能，也没有任何能量在介质中被吸收， 唯一影响是使光子发生小角度散射。相干散射在高原子序数材料和入射光子能 量较低时容易发生。２４第二章蒙特卡罗方法及对粒子｛禽运的模拟图２６瑞利散射原理相干散射光子可用Ｂｏｍ微分截面加上一个原子形式因子的解析表达来模拟。每个单位立体角．Ｉ二的原子ＤＣＳ是：等＝ｒＪ １＋ｃｏｚｓ２（０）［Ｆ（弘）】２ＰＥＮＥＬＯＰＥ使用的是Ｈｕｂｂｅｌｌ等人编成表格形式的非相对论形式因子。２．非相干散射（２３８）这里托是经典电子半径，ｑ是传递的动量幅值，而Ｆ（ｑ，ｚ）是原予形式因子。非相干散射又称康普顿散射（Ｃｏｍｐｔｏｎ Ｓｃａｔｔｅｆｉｎｇ），其作用过程为：当入 射光子和原子内一个轨道电子发生相互作用时，光予损失一部分能量，并改变 运动方向，电子获得能量而脱离原子。损失能量后的光子称为散射光予，获得能量的电子称反冲电子。考虑到相对康普顿散射占优势的光子能量范围，轨道电子的结合能很小，因此往往可忽略结合能的作用，把康普顿散射看成是光子 和处于静止状态的自由电子之间的散射。图２．７康普顿散射原理ＰＥＮＥＬＯＰＥ模拟光子非相干散射时采用的微分截面为东南大学博七学位论丈赢ｄ２ｃｒｓ，＝量２㈨。．Ｅ Ｊ睁－ｓｉｎ２ Ｚｐｚ，㈦，。，式中Ｅ’和Ｑ分别是散射光予的能量和方｜｛ｌＪ，Ｅｃ是被静止电子散射后位于方向ｙ 上的光子的能量。函数跹ＥⅣ）解释了束缚作用和多普勒展宽，这些作用只对低 能光子起作用，所以当光子能量大－ＢＪＬＭｅＶ时，并不需要乘以函数预ＥⅣ）。 ３．光电效应能量为加的光子与物质原子的轨道电子发生相互作用，把全部能量传递给对方，光子消失，获得能量的电子挣脱原子束缚成为自由电子（称为光电子）； 原子的电子轨道出现一个空位而处于激发态，它将通过发射特征Ｘ射线或俄歇 电子的形式很快回到基态，对于高能光子和高原子序数材料，特征辐射是高能Ｘ射线。图２．８光电效应原理根据能量守恒定律，入射光子能量Ｅ和光电子巴的动能满足如下关系Ｅ＝Ｅ。＋Ｕ（２４０）式中（＾为原子｛层电子的结合能，与原子序数和壳层数有关。Ｋ层和Ｌ层电子发生光电效应的概率最大，如入射光子能量大干Ｋ层电子结合能，则Ｋ层电子光电效应截面占原子总截面的８０％以上。ＰＥＮＥＬＯＰＥ中光电效应的总截面由Ｂｅ唱ｅｒ和Ｈｕｂｂｅｌｌ编制的程序ＸＣＯＭ 生成。对化合物和混合物，ＰＥＮＥＬＯＰＥ为了减小模拟时间和内存需求，使用衰 减系数的倒数来替代组成元素的原子截面。如果光电效应发生了，也需要根据 原子截面来决定哪一个元素被电离。第二章蒙特卡罗方法及对粒予输运的模拟光电效应截面的解析表达式为：驴｛Ｇｐｊ驭Ｅ，ｚ、ｌｆ臣、＞Ｅ，兰５（ｚ＋１５）妇旷旷ｌｋｅＶ＜Ｅ＜Ｅｃ（２．４１）ｅｘｐ（呜一Ｂｓｙ＋ＣｓＹ―Ｄｓｙ。式中伉是一个依赖于Ｚ的参数，ｙ＝ｌｎＦ，系数Ａｓ到Ｄｓ足描述每个元素的参数。函数盯美（Ｅ，ｚ）是一个适用于高能范围的经验性公式，它可以正确描述出光电效应截面在高能区的变化趋势（盯“＊Ｅ。）。 ４．电子对生成 当光子从原予核旁经过时，在原子核库仑场的作用下形成一对正负电子，此过程称为电子对生成。由于原子核质量大，它获得的能量可以忽略，因此可认为光予能量的一部分转变为正负电子的静止能量２Ｍ。ｃ２，另一部分作为正负 电子的动能反和量。Ｅ＝Ｅ＋＋占一＋２镌ｃ２（２．４２）显然只有当入射光子能量大于２执ｃ２＝１．０２ ＭｅＶ时，电子对生成作用才可能发牛。图２９电子对生成原理ＰＥＮＥＬＯＰＥ中电子对生成的总截面同样由ＸＣＯＭ程序生成。生成粒子的 初始动能根据Ｂｅｔｈｅ．Ｈｅｉｔｌｅｒ微分截面来取样。该微分截面已考虑了库仑和低能校正，其表达式为：ｄ如Ｏ＇ｐｐ＝芬ｃｚｚ【ｚ＋叩】；『２（吾一占）２氟（￡）十九（ｓ）］（２．４３）这里Ｅ是入射光予能量，ｇ＝１／ｅ（Ｅ＋ｍｃ２）而参数Ｂ是电子动能。％和ａ分别东南大学博土学位论文是经典电子半径和精细结构常量。量＿，７解释了原予电子范围内的对生成，函数 曲１和妒２描述了电子的屏蔽。 ５．各种作用的相对重要性 瑞利散射只对于极低光子能餐和高原子序数材料来说才有意义，在放射治疗过程中可忽略不汁。低能光子与高原子序数介质作用冈光电效应占优势而导致质量衰减系数很大，衰减系数随能量的增大而急剧降低，直到光子能量远远超过电子结合能，康普顿成为占优势的相互作用为止。在康普顿能量范围，高原子序数物质与低原子序数物质的衰减系数变化不大，但是它们随能量的增大 而减小，直到电子对生成变得有意义为止。在能量远大于１．０２ ＭｅＶ时，电子 对生成占据优势。表２．１给出了光子在水中发生各类相互作用的相对重要性［８１１３５１。表２．１不同能量光予在水中发生光电效应、康普顿散射和电子对生成的相对重要性光子能量（ＭｅＶ）ＯＯｌ０相对重要性光电效应（ｕＰｒ２＃）９５５０ ７ ０Ｏ康普顿效应（№跏）５ ５０ ９３电子对生成（ⅣＰＰ加）０００２６００６０ ０】５０Ｏ０ ０ ６】００９４４．００１０００ ２４００ １００ ｏｏ０ ０ Ｏ７７ ５０ １６２３５０ ８４§２．４电子输运的蒙特卡罗模拟§２．４．１电子作用机制电子（正电子）和介质发生的相互作用类型主要有：弹性散射、非弹性碰 撞、韧致辐射和正电子湮灭四种。 １、弹性散射 当电子靠近原子核通过时，受原子核库仑场作用而改变运动方向，既不发第二章蒙特仁罗方法及对粒子输运的模拟射ｘ射线，也不激发原了。这种不改变作用体系总动能的过程称为弹性散射。由于原予核一般比电子重得多，原子核基本不发生改变，结果只是入射电子的运动方向发生偏转，电子能量的变化甚微，通常可忽略。电子因质量微小，散 射现象特别严重。不仅被原于核散射，也可被核外电子散射，但其贞献与原了核的相比相当小。图２１０弹性散射原理如图２．１０所示，电子发生弹性散射后改变的只是运行方向，需要确定的是 散射极角０和方位角妒。模拟弹性散射最准确方法是采用分波法微分截面的直接蒙特卡罗模拟法，但工作量巨大，不具实用性。ＰＥＮＥＬＯＰＥ采用了Ｗ２Ｄ模型，可解析表达的微分截面在细节上和分波法ＤＣＳ＿；；＿ｌｉ完全一致，但两者产生的 弹性平均自由程、一阶和二阶传输平均自由程相同。Ｗ２Ｄ单散射ＤＣＳ计算公式如下：警＝瓦１“＝ｍ础）汜ａｎ，舫础Ｈ１＿Ｂ）描邶占∽训ｆ归（１－ｃｏｓ０）值可根据以下两个公式：（２４５）（２．４６） 、厶式中参数Ａ、Ｂ和∥ｏ的取值范围为Ａ＞０，Ｏ≤Ｂ＜ｌ，畦／“ｏ＜１。具体确定它们的ｆ坤一似坳＝（∥）＝轰‘２４７）ｊ：１∥２ｍ础）舡＝（∥２）＝舞一丧旺。ｓ，东南大学博士学位论文这样Ｗ２Ｄ的ＤＣＳ和分波法ＤＣＳ将产生同样的平均自由程和一阶、二阶传输平均自由程，两者产生的多重散射分布将十分接近。有ｒ ｐｗ２Ｉ舡），可以通过直接抽样法抽样得到∥。由于散射是相对于电于入射方向轴对称的，方位角伊可直接在（Ｏ，２兀）范围内随机抽样得到。２．非弹－眭碰撞 带电粒子与原子、分子或其它束缚态的电子发生库仑相互作用，使电子跃迁到较高的能级，这个过程叫激发。如果相互作用的结果使电子离开了束缚体 系而成为自由电子，则称这个过程为电离。在电离和激发过程中带电粒子的一部分动能转变为柬缚体系的位能，这种相互作用称之为非弹性碰撞‘２引。图２．１１非弹性散射原理ＰＥＮＥＬＯＰＥ描述电子非弹性碰撞使用的是根据Ｓｔｅｒｎｈｅｉｍｅｒ－Ｌｉｌｊｅｑｕｉｓｔ的通 用振荡强度（ＧＯＳ）模型计算出的Ｂｏｒｎ微分截面。ＧＯＳ模型考虑了束缚能、密度效应和相对论作用，可区分电子和正电子之间的差别。电子在介质中输运时，通过碰撞损失和辐射损失两种方式损失能量，每种 方式有各自ＤＣＳ来描述，这里考虑的只是碰撞损失方式。 碰撞阻止本领和连续慢化近似（ＣｏｎｔｉｎｕｏｕｓＳｌｏｗｉｎｇ Ｄｏｗｎ Ａｐｐｒｏｘｉｍａｔｉｏｎ，ＣＳＤＡ）是和非弹性碰撞相关的两个概念。一个介质的碰撞阻止本领定义为：ｓ。＝一警砘Ｊ缈斋ｄ∥ｄＷ眨。。，式中―ｄｏ―＇ｔｎ是每个（原子或电子的）散射中心发生非弹性碰撞的ＤＣＳ，矿是能 量损失，Ｍ是单位体积中的散射中心数日。根据定义，５：。ｆ是单位路径长度上第二章蒙特卡罗方法及对粒子输运的模拟的、Ｆ均能量损失。碰撞阻止本领的解析表达式包括两个重要的量：（１）平均激 发能Ｉ，（２）密度效应修正项艿。平均激发能足～个常量，它是介质激发能的 几何平均，加权值为相应的振荡强度，它表征介质的阻止本领。密度效应修正项解释了由于电子极化作用所导致刚止本领的减小。电子在传输过程中，其能量损失和运行路径长度都是波动的，故宅子单位 路径长度上的能量损失也是在均值附近波动的。ＣＳＤＡ忽略这种波动行，假设电子按照阻止本领连续损失能量。ＣＳＤＡ的弱点在于假定电子通过无穷小量逐 渐损失能量，而电子实际上是通过有限个重大事件（Ｃａｔａｓｔｒｏｐｈｉｃ）损失能量的。Ｂｅｒｇｅｒ指出在水中每个电离事件导致的平均能量损失为８０ｅＶ（电子能量为 １ＭｅＶ），４５ｅＶ（电子能量为ｌｋｅＶ），３３ｅＶ（电子能量为０．１ｋｅＶ）。这样低能量 电子就会在不多的重大事件中损失掉所有能量，此时ＣＳＤＡ不再适用。在ＰＥＮＥＬＯＰＥ中碰撞损失的ＤＣＳ表达式为：―ｄ２ａ―ｃｏｌ：一２，，７ｅ２上―ａｆ（Ｑ―，Ｗ） ａＷｄＱ ｍＰ‘ＷＱｄＷ（２．５０）式中Ｗ 荡强度。Ｆａ………ｕ，，．：，Ｑ是反冲能（即靶在碰撞后的动能），生鬻婴是通用振＇７１ＷＷ－Ｑ３］Ｚ［［ｔｌ－上绘制的望篆笋曲面称为Ｂｅｔｈｅ，该曲面包含了电子和原子或分子发生非弹性碰撞的所有相关信息。ｆｔ满足如下相加法则：Ｊｏｆ。生！翌！兰！ｄ∥：Ｚｄ∥（２ ５１）式中Ｚ是原子数。平均激发能可通过下式进行计算：ｚＩｎ』：广Ｉｎｗａｆ（Ｏ＝ｏ，ｗ）ａｗｄ０ｄＷ（２．５２）式２ ５ｌ和２ ５２决定了Ｂｅｔｈｅ阻止本领。在高Ｑ值区域，Ｂｅｔｈｅ面除缈４Ｑ以外都为零。这意味着在高反冲能时，靶电子可作为静止的自由电子，此时Ｂｅｔｈｅ面退化为沿胆Ｑ的Ｂｅｔ＿ｈｅ脊。这样碰撞微分截面可分为两个部分，一是沿着Ｂｅｔｈｅ脊发生的激发（相当于自由两体碰撞），称为近碰撞，另一部分发生的碰 撞称为远碰撞（Ｑ＜矿）。碰撞微分截面计算公式如下：东南大学博士学位论文等＝烈努＋万ｄＯ＇ｄ１）｝ｉ式准确地描述了非弹性碰捕的平均性质。３．韧致辐射眩ｓｓ，当带电粒子在原子核附近穿过时，入射粒子在原子核电场中产生加速运动。按经典物理学观点，带电粒子将以正比于其加速度平方的几率辐射电磁波，这就是韧致辐射。图２１２韧致辐射原理ＰＥＮＥＬＯＰＥ对电子和正电子的韧致辐射用Ｂｅｔｈｅ．Ｈｅｉｔｌｅｒ微分截面模拟，该微分截面考虑了指数屏蔽和库仑校正，并根据实际试验结果进行了修改以提高 它在中低能区的可靠性。电子的角偏转假定已被包括在弹性散射ＤＣＳ中，故此处不模拟电子偏转方向的改变。电子Ｂｅｔｈｅ―Ｈｅｉｔｌｅｒ微分截面可表示为：！芝星：＃口ｓｚ（ｚ一叩）上【（１＋（１－８）：ｘｍｌ―４ｚ）一要（１一ｇ）（中２－－４五）】ｎ￡ ￡３（２．５４）这里占；Ｗ／（Ｅ＋ｍｃ２）是发射光子能量占电子能量比值，形是电子损失的能量，口“１／１３７是精细结构常量，，ｅ是经典电子半径，工是高能库仑校正项。此ＤＣＳ 仅适用于电予能量在光子发射前后都远大于其静止能量（小ｃ２＝５１１ ｋｅＶ）的情况。 （ｐｌ和中２是屏蔽函数，都包含原子形式因子以ｇ），其表达式为：Ⅷ５丽１‘２?５５’这里口是传递的动量而Ｒ是屏蔽半径，它作为一个可以调节的参数而存在。电了的韧致辐射微分截面的最后形式为：第二章蒙特卡罗方｝！圭及埘粒子输运的模拟警劬５犯＋Ｖ牝，＋知ｓ））盯Ｓ ＼Ｆ（２．５６）／仍（ｓ）和妒：（ｓ）中都包含Ｊ，低能校正因子来校ＪＥＢｏｒｎ近似在低能时低估ＤＣＳ的 错误，低能校正网子由对Ｓｅｌｔｚｅｒ和Ｂｅｒｇｅｒ的Ｅ＞５０ ｋｅＶ辐射阻止本领拟合得到。在高能条件Ｆ，】Ｅ电子ＤＣＳ可简化为电子的ＤＣＳ：但在ｌ－ＩＪ问能区，它小于电子的ＤＣＳ，可通过对电子ＤＣＳ加入校正因子来计算：豢叫ｚ固毕ａｏｅ式中的校正因子等于正电子和电子的辐射阻ｌＬ本领的比值。发射光子角分布为：眩ｓ，，发射光子的角分布可根据一个半经典参数推导，在高能区（ｙ＜＜ｌ，ｘ“１）础，兰蓦高籍ｚ式中ｒ是极角。旺ｓ鼬产而１小嘉４．正电子湮灭（２５９）图２１１３正电子湮灭原理一个带电粒子与其相应的反粒子相互作用，两者的静止质量转化为两个光东南大学博士学位论文子的能量，此为湮灭过程。正电子在介质中运行时会和介质中电子发生湮灭作用，称为Ｉ卜电予湮灭。通常假设介质中的电子是自由和静止的，则电子的束缚作用可忽略。ＰＥＮＥＬＯＰＥ使用Ｎｅｌｓｏｎ等计算的微分截面描述正电子湮灭。在实验室坐 标系中，该ＤＣＳ表达式为：：等＝赢阢删刊】蹦＝一㈣）２＋（ｙ２＋４７＂＋１）了１―７１ｃｏｓ（ａ＿）＝（ｙ２～１）。”（ｙ＋１―１／‘－） ｃｏｓ（Ｏ＋）＝（ｙ２－１）。”Ｉｒ＋ｌ－Ｉ／Ｏ－ｆ）ｌ眩ｓ。，（２ ６１）式中正电子的总能量（以其静止能量为单位）为，＝１＋Ｅ／（ｍｃ２），另外， ｆ＝ｔＬ／（Ｅ＋２ｍｃ２），下标“－”表示最低能量的光子。根据能量和动量的守恒定律，两个发射光子的方向极角为（２．６２）（２．６３）§２．４．２浓缩历史方法在§２．３中介绍的光子输运模拟方法称为直接模拟法或详细模拟法，即对粒 子输运过程中所发生的每次相互作用都采用单独的作用模型和截面进行模拟， 这是模拟粒子输运问题最精确的蒙特卡罗方法。但对于电子输运，此法碰到了计算方面的网难，在现有条件下几乎不具实用性。原因有三：（１）相同能景的电子截面比光予要高５到６个数量级，即在运行同样距离 过程中，电子和介质发生的相互作用次数比光子多５－－６个数量级。通常，具有一定能量的电子在停止输运以前和周围介质发生的相互作用次数在１０５～１０６之间１８６Ｉ。（２）电子每次作用损失的能量也比光子小５个数量级，意味着电子损失掉第二章蒙特卡罗方法及埘粒子输远的模拟所有能量需要发生更多次数的相互作用。在单次作用中，电了能量和运动方向 的改变都很小，电子轨迹形状是依靠大量的作用累积面成。（３）虽然电子在每次作用中能量和偏转角度的改变都很小，但它们却很重 要，否则于累积效应会造成较大的误差。所以不能采用简单的分布（如高斯分布）来近似描述电子偏转角分布。显然采用直接模拟法模拟电子输运过程的运算量过于巨大。为此Ｂｅｒｇｅｒ提出厂浓缩历史方法㈨蚓（Ｃｏｎｄｅｎｓｅｄｈｉｓｔｏｒｙｍｅｔｈｏｄ），该方法的基本思想是：把真实的物理ｌ二的随机游动划分为若干历史阶段，每一历史阶段包括好多次游动；也就是说，人为地把多次相互作用合并为一次作用，作为一步随机游动处理，这一步称为电子多次散射步或电子步，其能量和飞行方向的转移概率由近似的多次散射理论给出。应用浓缩历史法模拟电子输运能否得出理想结果，主要取 决于对电子历史阶段的划分是否合理，也取决于在各历史阶段中所采用的统计理论是否准确。浓缩历史法极大地提高了模拟效率，使得电子输运的蒙特卡罗模拟具备了 使用价值，在蒙特卡罗模拟程序中已获得Ｊ一泛应用。Ｌａｒｓｅｎ证明丁浓缩历史法的结果是波尔滋曼方程在小步长限制下的解１８６］。但在使用中必须注意到浓缩历史方法是以如下假设为基础：（１）假设电子在每个电子步中沿直线前进，而实际上电子连续不断地和介质 发生相互作用，每次作用都会导致电子发生极其轻微的偏转。显然，电子步长取得越大，这种假设带来的误差也越大；同样，在高原子序数的介质中此假设的误差也大于低原子序数的介质。（２）浓缩历史法用多次散射分布表示电子运行一步后运动方向的改变，多次 散射理论是～种近似理论，有可能引入系统误差。 （３）假设电子以～定的比率连续损失能量，即阻止本领（ｄＥ／ｄｘ）。实际上电 子每次损失的能量是有波动的，服从一个分布，为了方便大家常常用高 斯分布来近似表达，阻止本领只是此分布的均值。显然电子步长越小，此假设引起的误差也越小。 根据实现方式的不同，浓缩历史法又可以分为两种类型：’Ｓ东南大学博士学位论文第一类方法是将电子历史轨迹直接划分成若干电子步，在每个电子步中考 虑所有的相互作用，电子步＆预先设定，采用此类浓缩历史方法的蒙特乍．罗模 拟程序有ＥＴＲＡＮ和ＩＴＳ等。 第二类方法相当于混合处理方法，即对电子影响较大的相互作用（硬事件）采用直接模拟法，剩下的影响较小的相互作用（软事件）由多次散射步来模拟。 通常定义一个能量闽值或偏转角度闽值，将能量损失或偏转角度大于闽值的作用归为硬事件，将能最损失或偏转角度小于阙值的作用归为软事件。采用此类浓缩历史方法的蒙特卡罗模拟程序有ＥＧＳ４和ＰＥＮＥＬＯＰＥ等。无法准确模拟大能最损失的相互作用是第一种方法的主要问题，第二类方法则可对生成的次级光子和电子进行准确模拟，更接近于真实情况。故本文的 研究只涉及第二类浓缩历史方法。§２．４．３多次散射理论电子在介质中运行长度ｓ的过程中，会与介质发生大量的相互作用，因此方向在不断改变。多次散射理论就是求取电子运行距离Ｓ后其偏转角的分布。 电子多次散射理论从Ｗｅｎｔｚｅｌ工作算起已有八十余年的历史，早期工作研 究对象是一束单能电子穿过金属薄层后的角分布问题，主要的理论有Ｗｉｌｌｉａｍｓ的小角理论、ＭｏｌｉＳｒｅ的傅立叶一贝塞尔变换法‘“５１和Ｇｏｕｄｓｍｉｔ和Ｓａｕｎｄｅｒｓｏｎ的勒让得展开法‘７孙。它们都采用了小角近似和忽略电子能量损失，理论结果和 实验相比较也很符合。Ｌｅｗｉｓ从连续慢化的电子迁移方程出发，利用球谐函数展开得到空问矩方程‘７３１，在小角近似下统一地处理了前面作者的结果，Ｌｅｗｉｓ理论对大角散射同样是正确的。 目前在通用蒙特卡罗模拟程序中使用较广的两种多次散射理论是ＭｏｌｉＳｒｅ多次散射理论和Ｇｏｕｄｓｍｉｔ－Ｓａｕｎｄｅｒｓｏｎ（ＧＳ）多次散射理论，前者用于ＥＧＳ４等采用第二类浓缩历史法的程序：后者则多用于ＥＴＲＡＮ、ＭＣＮＰ等采用第一类浓缩历史法的程序中。ｌ、Ｍｏｌｉｅｒｅ多次散射理论电子运行距离ｓ后，偏转极角ｐ分布为：第二章蒙特卡罗方法及对粒子输运的模拟粤：旦压（２ ６４）％…叫∥＋窄＋铲扣?］∥ｃ舻ｉｌ Ｊ０删舭批＿ｙｑ２可由下面公式得到：（２ ６５）Ｉ彳２?ｎ等］“眩ｓｓ，式中乾和Ｂ是粒子能量与步长ｌ的函数，ｄｏ是零阶贝塞尔函数。Ｍｏｌｉ６ｒｅ参数Ｂ曰一ｔｎＢ＝?ｎ（薏］２＋ｔ―ｚｙ式中，为欧拉常数，７＝０．５７７２１５…，屏蔽角‰由下式得到ｃｚ?ｓ，，《＝嚅［¨３ｎ，６（南）２］光速。，参数ｘ；足和一定步长中作用概率相关的量，可计算如下＠ｓｓ，式中ｆ是电子动能，单位为电子静止质量，Ｚ是原子数，卢是电子速度，单位为Ｊ。２＝０．６００９２４２一Ｉ ｒ＋ｌ巧卜式中，Ａ是原子质量，ｓ是步长。眨ｓ∞２、ＧＳ多次散射理论 根据ＧＳ理论，电子越过距离ｓ后，偏转极角０的分布为一勒让得级数：吆（刚）如＝喜（，＋ｅｘｐ（一ｓ鲁）只＠№Ｊ神，０（２－７。）式中ｇＯ＝ｃｏｓＯ，Ｐｌ（∞）是，阶勒让得多项式，＾是弹性散射平均自由程（ＭＦＰ），盛是传输系数，它对整个分布起到了决定性作用，其定义如下：ｇ，：Ｉ―Ｃ，（，国只（彩）ｐ（国）（２．７１）东南大学博士学位论文式中“∞）是概率密度函数（ＰＤＦ），由单散射事件微分截面盯（∞）计算得到ｐ（ｃｏ）＝丁堕盟一Ｉ，ｄ（ｏ‘口（∥）３、两种理论比较（２．７２）Ｍｏｌｉｅｒｅ理论实现简单，用两个变景表示ｒ所有的参数（能量、原子序数、 偏转角度、运行距离），编码简单、极角抽样速度快，运行效率较高，非常适合根据任意大小步长抽样多次散射角。但存在步艮的上下限，步长下限（大于２０个散射事件）将降低不同介质交界处模拟精度。在低能范围，由于采用了小角 度近似，即ｓｉｎＯ＝０，限制该理论在多次散射角小于约２０度时有效，导致了步 长上限存在，限制了整个模拟效率的提高，当电子在高原子序数材料中传输时情况尤为严重。在高能范围Ｍｏｌｉ＆ｅ模型的不足在于没有考虑量子自旋和相对 论效应。ＧＳ理论比Ｍｏｌｉ＆ｅ理论更精确。可计算比Ｍｏｌｉ６ｒｅ理论更小的步长，对任 意大小的多次散射角都有效；由于未采用小角近似，步长上限可取得较大：可以和任何单次弹性散射截面联合使用，可以包括自旋和相对论效应；可以容易地分别处理电予和正电子。但由于Ｇｓ理论采用了勒让得级数，当步长较小时， 需要高阶保证勒让得级数的精度‘ｍ１。所以ＧＳ理论实现困难，编码复杂，计算 效率低。ＭＣＮＰ等属于第一类浓缩历史法的程序采用了如下快速算法：由于电 子步长是预先确定的，故可事先计算好固定步长的ＧＳ分布并保存起来，模拟 电子输运时只需要进行插值运算。 两个模型在高能和低原子序数材料小散射角条件下表现都很好，这也正是 放射治疗一般碰到的情况。但随着ＧＳ理论快速算法的出现，Ｇｓ理论更有生命 力。Ｍｏｔｉ６ｒｅ理论一直是ＥＧＳ４无法进一步提高模拟精度的主要原因，ＥＧＳ系列最新版本ＥＧＳｎｒｃ已改用了ＧＳ理论。§２．４．４电子步传输算法采用浓缩历史法模拟电子输运时，电子运行一个电子步之后，依据步长和 多次散射角来计算电子空间位置的算法叫电子步算法‘５７１（Ｅｌｅｃｔｒｏｎ－Ｓｔｅｐ第一璋蒙特幸罗方法及对粒子输运的模拟Ａｌｇｏｒｉｔｈｍ）ｎＰＥＮＥＬＯＰＥ采用的电：子步算法称为随机链式算法ｍ１（ＲａｎｄｏｍＨｉｎｇｅＭｅｔｈｏｄ）。假设电子步的步长为Ｓ，偏转极角为０，如图２ １４所示，随机链式算 法将此电子步分成两个子步，具体模拟过程为： （１）在（０，１）范围内取一随机数６（２）沿电子初始方向运行第一予步，子步长度为Ｉ：－＝－ＳＸ毒；（３）根据多次散射分布计算出偏转极角０，在（Ｏ，２ｎ）范围内随机抽样方位角妒，改变电子的运行方向。注意偏转极角０是根据整个电子步长Ｓ计算出来的。 （４）电子在新方向上运行第二个子步，子步步长为ｓ―Ｔ。阁２１４随机链式算法示意图§２．４．５电子步长和边界越过第二类浓缩历史方法只对硬事件进行详细模拟，各个硬事件之间发生的软 事件都采用多次散射步进行模拟。ＰＥＮＥＬＯＰＥ中硬事件的平均自由程为：上胪Ｉｌ上掣十上础 上豫＋＋忐（２．７３）等式右边分别为弹性散射、非弹性碰撞、韧致辐射和正电子湮灭四种作用的平 均自由程（大于阈值部分）的倒数。两个连续硬事件问的距离ｔ服从指数分布：肌）＝嘉ｅ醑万ｔ）（２．７４）每次只有一种类型的硬事件发生，每种类型发生的概率为：东南大学博士学位论文列６）ｐ一２万具体发生何种类型的作用，随机抽样得到。（２．７５）电子在硬事件之问发生的所有软事件采用多次散射步模拟，ＰＥＮＥＬＯＰＥ采用ＧＳ理论计算多次散射分布，采用随机链式模型传输电子步。如果模拟系统由非均匀介质组成，则电子步算法中需要处理边界越过情况。①吆Ｉ，’一．｜。－多矿‘ｒ’ｒＩ斗，ｃｄ―ｒ蛞ｄｆ图２１５电子边界越过处理如图２．１５所示，假设电子步开始时电子位于介质１，经过步长ｔ以后电子位于介质２，即电子将在此步中越过边界。ＰＥＮＥＬＯＰＥ中处理方法如下：（１）随机抽样得到第一个子步的步长Ｔ，电子沿当前方向运行ｒ后如位 于介质１中，则计算多次散射偏转角，电子方向发生偏转：电子沿新方向运行到边界结束当前电子步：（２）如电子越过边界后仍来运行完距离ｒ，则结束当前电子步，不再计 算多次散射偏转角，在介质２中重新开始一个新的多次散射不。 （３）只要电子越过边界后。都要重新取样到下个硬事件的距离ｔ。 注意，当电子进入新体元后必须再向前运行一个极其微小的距离，比如ｌＡ， 以避免由于截断误差导致电子在边界上发生振荡现象。因为ｆ＜＜磷’，且ｌ＿（ｃｏｓｚ）扣’反映了多次散射偏转角的均值，故有：１－ｃｏｓｚ）ｐ’＿ｌ－ｅｘｐ（一方）ａ方位７６）第二章蒙特卡罗方法及列粒子输运的模拟输入各种参数一模一一拟 一初 一始 ～化Ｌ弓一一入． 几 何 结 构一７丌蛳一７ｒ帮ｌ叫－ｆｆｊ朋ｔ旦Ｔ侵１％，’设置初始粒子的能量Ｅ、位置ｒ、运动｜方向ｄ，所在的体序号和材料序号。’，１，１初始化二次堆栈１ １０．；１设置下一事件为软事件 ，’――●－‘。。。。。｝．】●ｌ｜＿确定下一事件步长Ｄｓｒｋ判别径迹是｛ §跨越界面？ｌ是ｉ将Ｄｓ改为到交界点距离ｉＩｒ＝ｒ十Ｄｓ奉ｄ； ｛ｌ、ｌ１５日ｒ＝ｒ＋ＤＳ＊ｄｌ ●｝１５是否速逸？Ｉ●。士｛悬．是Ｉ模拟下一事件转入新的体Ｅ≤Ｅ＋ａｂｓ？ｉ从堆栈中取一个次级粒子并设置其初始状态Ｉ●是堆栈中次级粒子数＞０７－～●－－＿－是Ｉ模拟的初始粒晕数＜设定值？统计结果并输出ＰＥＮＥＬＯＰＥ模拟电子输运的流程图图２．１６４东南大学博士学位论文义由于电子发生偏转的位置是在两个硬事件之间均匀取样的，它发生的概率为ｓ／，。这样电子到达界面时，其平均偏转角为：，一（ｃｏｓＺ）＝∞ｃｏｓｚ）㈣】兰嘉第二类浓缩历史法模拟电子输运的结构。弦?，，陶２．１６给出了ＰＥＮＥＬＯＰＥ模拟电子输运的流程图，这是一个完整的采用§２．５模拟参数粒子输运的蒙特卡罗模拟是一个非常复杂的过程，影响模拟精度和速度的 因素有很多，所以需要使用一些参数来控制。ＰＥＮＥＬＯＰＥ中的相关参数有：ｃ１、 Ｃ２、Ｗｃｃ、Ｗｃｒ、Ｅａｂｓ、Ｅａｂｓｐｈ等六个，现分别介绍如Ｆ： ＰＥＮＥＬＯＰＥ中，弹性散射的硬事件平均自由程由下式确定：棚…｛２ｅｌ（Ｅ），ｍｉｎ［Ｇ弛ｃ：劫眩㈣该平均自由程减小意味着有模拟更多的弹性散射硬事件，精度增加，效率下降； 反之，则意味着将更多的弹性散射作为软事件用多次散射步进行模拟，精度将 卜降，效率上升。如图２．１７所示，在低能区，‘“Ｉ’取ｋｌ，在中能区，当钿趋近一个固定值时，Ｃ１决定了露１的大小，此时：ｌ―ＣＯＳ２＂）＝ｌ―ｅ砸一争）“ｃｌ ＾（２．７９）即ｃ】表示两个硬事件之间电子由于多次弹性散射导致角偏转的平均值。由于Ａ。随着能量增加而迅速变大，特别是高能区。此时露’将主要由Ｃ２决定，它表示两个连续硬弹性碰撞之间平均能量损失的最大值。４２第二章蒙特卡罗方法及对粒子输运的模拟一。＿【薯ｕ＼∞－Ｉ；【ｐ№｝ｕ口一ｘ口图２１７弹性散射硬事件平均自由程和弹性散射平均自由程、弹性散射一阶传输 自由程以及单位长度能量损失之间的关系。此处Ｃ１－－Ｃ２＝０ ０５。Ｗｃｃ是决定非弹性碰撞硬事件的能量阈值，即能量损失大于此阈值的非弹 性碰撞将作为硬事件模拟，反之作为软事件模拟。同样Ｗｃｒ是决定硬韧致辐射 硬事件的能量闽值。此两个阚值增大，则精度下降，效率增加；反之，精度变大，效率下降。Ｅａｂｓ是电子传输截止能量，即电子能量如小于Ｅａｂｓ，则此电子不再进行模 拟而是就地沉积。Ｅａｂｓｐｈ是光子传输截止能量，当光予能量小于Ｅａｂｓｐｈ时，光子就地沉积。这两个参数设置越大，模拟效率越高，但引入的误差也随之增 加。这两个参数的取值必须和能量沉积的空间分辨率相适应。§２．６方差降低技巧在模拟粒子输运问题方面，蒙特卡罗方法的系统误差很小，而统计误差对 计算精度和速度影响较大［９０１。由式２．５可知统计误差和计算结果的均方差成正 比和模拟粒子轨迹数目的平方根成反比，增大模拟粒子轨迹数Ⅳ和降低均方差盯是减小统计误差的两条途径。但前者效率不高，在均方差固定条件下，要提东南大学博士学位论文高精度一位数字，就要增加１００倍的模拟量，因此单纯增加Ⅳ１ｉ是有效的办法；另一方面，均方差减小一半，统计误差也就降低了～半，相当于Ⅳ增大４倍的收益，因此在有关蒙特卡罗模拟粒子输运问题的研究中，有关方差降低技巧的内容是很多的【２１１３ｎ１４ｊ【２２峭Ｊｉ＝｝６Ⅱ”１１”１。在实际应用中，人们利用方差降低技巧的目的在于减少模拟的粒子轨迹数（保持同样统计误差），从而减少模拟时间、提高 模拟的效率。因此，Ｂｉｅｌａｊｅｗ和Ｒｏｇｅｒｓ等人将所有可减少模拟时间（保持同样统计误差）、提高模拟效率的技巧都称为方著降低技巧‘５”。 方差降低技巧的主要问题在于通用性差，需要具体问题具体分析，对某些 问题有效的技巧对别的问题却可能用处不大（“甜。下面是一些适用于粒子输运问题的方差降低技巧。１、减小几何边界检查次数 蒙特卡罗方法模拟粒子在非均匀介质中输运时，粒子每前进一步都要检查粒子是否离开了当前模拟单元，因为粒子离开当前模拟单元后介质类型可能发生改变，所以这步通常是必不可少的。每次检查所需时间与模拟系统几何形状 的复杂程度有关，如几何形状复杂，化在这方面的时间就很可观了。电子的平均自由程很短，即使采用了浓缩历史法后，电子的多次散射步长也往往小于模拟系统的几何分辨率。此时电子往往要运行多个电子步后才会离 开当前模拟单元，这样就无须每次都进行几何边界检查。比如ＥＧＳ４就采用了 减小ＪＬ何边界检奁次数技巧，计算出电子到当前单元边界的最短距离／）ｎｅａｒ， 每运行一步就将此距离减去步长。当距离Ｄｎｅａｒ为正数时，就不需要进行几何 边界检查；只有当其为负数后再开始几何边晃检查。 这个技巧对计算精度没有影响，其效率提升来自于计算一次Ｄｎｅａｒ替代多 次进行的ＪＬ何边界检查，节约的时间应该是两部分的差。具体提升情况取决于 Ｄｎｅａｒ的计算难度、电子步长与几何分辨率的比等因素，通常该技巧的效率提升值在１ ３４左右１５”。２、区域抛弃和距离拒绝 粒子输运问题求解的结果往往是沉积能量的分布。每个模拟单元是最小的 模拟单位，其内部能量沉积的变化不需考虑。和光子相比电子的射程较短，低第＿二章蒙特卡罗方法及对粒子输运的模拟能时更为严重。如果菜电子的能最低到不足以使其离开当前模拟单元，显然就 没有必要继续进行模拟，将其能量就地沉积和继续模拟的结果是一样的，这就是区域抛弃的核心思想。具体做法是计算出当前能量下电子所能运行的ＣＳＤＡ距离，并将其和电子至模拟单元边界展短距离Ｄｎｅａｒ进行比较，如小于Ｄｎｅａｒ则就地沉积而不再模拟。 通常我们需要计算的是某些部分的能量分布，这些部分称为感兴趣区域。 此时可采用距离拒绝技巧加速模拟速度，如果电子的能量无法保证其运行到感兴趣区域，该电子对于计算结果就几乎没有贡献，故可立即抛弃不再模拟。 以上两个技巧实际上引入，极小的误差，因为电子会发生韧致辐射，生成 的次级光子射程很长，它是有可能离开当前模拟单元或到达感兴趣区域的。不 过发生此种现象的概率相当小，一般小于Ｏ．１％。除非专门研究韧致辐射光子的剂量分布，此技巧的误差可以忽略。如果整个模拟系统中存在多种介质，电子的ＣＳＤＡ距离不易准确计算，因 而此时距离拒绝技巧用处不大；而单一介质时，此技巧可获得较高的效率提升。３、强迫作用光子发生相互作用的概率分布为：ｐ（Ａ）ｄ名＝ｅ－，Ｚｄ２（２．８０）式中Ａ是光子运行的距离（单位是平均自由程数），０≤五＜ｍ。＾可由下式抽样得到：五＝一Ｉｎ（１一善）（２．８１）式中ｆ是０和１之间均匀分布}