CN1346291A - Monte carlo法医学加速器光子输出相空间的相关方框图表示法 - Google Patents

Abstract

提供了一表示发自医疗加速器的光子的计算相空间的方法。该加速器用在光子远距离放射治疗中。该方法再现了光子的能量分布和轨迹,这些光子是产生自韧辐射靶的光子和由加速器头内的元件散射的光子。该方法从源到高达几厘米的辐射范围内再现出能量和方向信息,这样,该方法有望很好地适应由不同的厂商制造的加速器。该方法的计算快速而高效,对大多数场尺寸有80%总采样效率或更高。计算成本不依赖用在治疗方案中的射束数目。

Description

Monte Carlo法医学加速器 光子输出相空间的相关方框图表示法

本专利申请对提交于1999年3月8日、专利号为60/123,242的美国临时专利申请要求有优先权。

按照美国能源部与加州大学之间关于Lawrence Livermore国家实验室运作的合同No.W-7405-ENG-48，美国政府对此发明具有权利。

                   发明背景

发明领域

本发明涉及辐射治疗，更具体地说，涉及通过医学加速器MonteCarlo模拟获得的聚集光子能量和角分布的方法。

相关技术描述

Monte Carlo方法对计算辐射剂量很有用，这是因为它们对癌症辐射治疗中遇到的全部情况有精确预测剂量的能力。它们能正确地描述自线性加速器发出的电子束和光子束，还有从近距放射治疗源发出的光子、电子和α粒子的特征。进而它们可精确地模拟通过波束成形设备或对近距的治疗源屏蔽设备的粒子传输，然后考虑病人身体表面的不规则性和内部组织结构，包括任何金属移植和修复术，可类似地模拟穿过病人的粒子的传输来确定沉积在病人体内的沉积剂量。这个揭示提供了使用线性加速器进行光子治疗的方法。

虽然Monte Carlo方法是预测辐射剂量最精确的手段，用于日常临床应用的一般Monte Carlo软件包的使用现在不合实际，部分原因是由于计算机的限制，特别是由于获得收敛计算而需要的运行时间。一般Monte Carlo软件包的冗长的运行时间保证了特别为手头上的问题优化的算法的开发。这个问题可分为三个步骤：(1)确定从加速器头发出的光子辐射的特性；(2)通过波束调节器(如果有的话)和病人上方的空气把该辐射传给病人；以及(3)原始辐射以及在加速器头和病人之间产生的或分散的辐射在病人体内的实际传输。

                发明概要

本发明的目的之一是提供一种特别是用于从医学加速器MonteCarlo模拟获得的聚集光子能量和角分布的方法。该方法将输出表示为一系列相关矩形图因而非常适合用于Monte Carlo程序的光子源包，这样的Monte Carlo码用于确定光子治疗的剂量分布。

该方法描述了光子能量谱分布随着离辐射束中心轴距离的增加而发生的等中心面(isocenter-plane)的变化。辐射产生在韧致辐射靶及由加速器头内的各种治疗器部件散射而成。把用该算法算出的等中心能注量(energy fluence)与潜在的全物理Monte Carlo光子相空间的能注量比较，表明能注量误差小于开放场尺寸范围内最大能注量的1％。档板边缘半阴影(penumbrae)的比较显示光子的角分布被精确再现。算法的Monte Carlo采样效率(生成光子的通过准直器档板的部分)对一个10×10开放场大致为83％，对一个40×40开放场上升为大约96％。和潜在全物理相空间文件1GB大小比较，在一给定能量时，一个典型医疗加速器数据文件大小约为150KB。

               附图简要描述

图1表明示出一个高能电子束撞击一个薄，高Z靶。

图2表明示出光子如何向初始点回踪，表示加速器头的哪部分对输出能注量流量有贡献。

图3显示从Varian 6WV加速器发出的能注量分布。

图4A-D显示子源的直径和其相对挡板开口的位置导致在等中心面处不同的照射面积。

图5A-D表明了挡板的位置对来自主准直器和矫直过滤器的能注量的影响。

图6A-6F显示轴上光子能量分布强烈地依赖于产生光子的治疗头部件和这些子源在距加速器中心轴20cm处的等中心面上的能量分布。

图7A和7B显示靶光子概率分布的近视图。

图8A-8C显示源位置对光子能量的影响。

图9A-9E显示用于生成主准直器光子的起始位置的径向和角形分布的近视图。

图10A-10F显示用于生成矫直过滤器光子起始的径向和角形分布的近视图。

图11A-F为图8的修改形式，表明用于描述构成治疗器的子源的相关矩形图。

图12A-C显示几个光子能量的概率箱(probability bin)高度的径向变化。

图13显示示于图12的分段线性函数的“结点”是从原始数据直方图化过程中得到的箱中心和箱高度。

图14显示在一个40×40场位置由上挡板确定的场边缘处矫直过滤光子的光子能注量的陡降。

图15显示生成光子轨道的算法，光子轨道可与挡板边缘相交。

图16显示对方位各向异性的子源的大多数粒子导致中心外分布的最小化R_closest。

图17A和17B显示假采样技术使任何分布可被有效地采样。

               发明的详细描述

本发明的公开集中在背景部分中所描述的步骤1上，确定发自加速器头的光子辐射特性。该步骤称为全部计算中的“病人—独立”部分，这是因为加速器头的辐射输出不受下游存在的波束调节器的影响。在背景部分中描述的步骤2和3必须根据治疗方案而变化，可称为全部剂量计算中“病人—依赖”部分。步骤1有利于详细的治疗有两个原因。第一个原因为在每个新的治疗方案期间不必特意重新模拟加速器输出；时间被更好地用在问题中的特定病人部分。第二个原因为实用性：医学加速器头的Monte Carlo模拟需要难以接受的长时间运行(使用单片CPU机器要几十小时)，而且要存贮和操纵(按照日常临床)最终的千兆字节大小的数据文件是不切实际的。

问题就变为开发出一种方法产生光子，其能量和轨道代表加速器头的潜在的Monte Carlo模拟的能量和轨道。由此产生三个问题。第一，由这样算法产生的光子的能量分布和角形分布必须能精确地模拟在潜在的千兆字节大小数据文件中建立起来的能量和角形分布。第二，用于产生这些光子的算法必须在计算上高效，使整个剂量计算的运行时间不会受到影响。第三，如果该算法依赖于压缩千兆字节大小文件为较小的数据文件以便在运行时间内使用，该压缩过程应是自动的，这样用于大量的加速器数据文件可以很容易地生成。

几种参考文献建议给病人的剂量应精确到5％以内，并且计算的剂量应精确到给病人最大剂量的2％以内。所以，相空间信息至少也必须一样精确，而且用于来自状态空间的采样粒子的任何算法必须生成同样精确的数据。进一步讲，采样算法必须计算上高效而不会影响计算效率和整个运行时间。

模拟问题的源特性在图1中说明。一束电子撞击一块由高z材料构成的薄靶10。产生的韧致辐射光子被圆锥形准直仪12准直。光子束穿过矫直过滤器14，它的中心较厚，这样韧致辐射光子分布的中心部分相对于边缘部分就会受到衰减。期望的结果是在水影(waterphantom)特定深度处平坦的剂量分布。虽然剂量分布不均匀，光子能注量和能谱作为与离开中心射线距离的函数而变化。这有两个原因。第一原因为韧致辐射光子能量随散射角增大而减小。第二原因为矫直过滤器在波束中心比在边缘去除更多的低能量光子。

另外，不可忽略的辐射量从准直仪和矫直过滤器散射，穿过挡板16到达病人面18。合成辐射场特点在于具有几种韧致辐射和散射光子分布。必须理解这些分布以便了开发出一种用于输入到Monte Carlo剂量计算码中去的有用的源算法。

该方法使用一系列相关矩形图来再生成从加速器治疗头原型Monte Carlo模拟获得的相空间信息。当与潜在的全物理Monte Carlo模拟给出的误差比较时，等角点处的全部能注量误差小于照射区域最大能注量的1％。出自于任何一个硬部件，诸如，靶或矫直过滤器的光子能量流的误差在照射区保持在最大能注量的1％标准以下。

算法的计算效率有两个基础。第一个是概念的：算法设计成仅产生有穿过被档板限定的孔的高的概率的光子。产生不会到达病人的光子所花费时间极少。第二个基础是使用假采样(alias-sampling)技术，该技术允许使用精细放箱(finely-binned)分配而没有码运行时间的相应增加。

该算法使用的相关矩形图可以从产生于加速器头的Monte Carlo模拟的数据文件自动生成。(即：不需要用户介入)。

本公开涉及从加速器头发出的光子。注意，这些光子一旦生成，接着就必须通过“病人-依赖”波束调节器，及病人上方的空气柱被传输。

方法与材料

使用Monte Carlo辐射传输模拟码研究加速器头的输出已被很好地确立。图1示出模拟过程的概略图。具有在(或近)额定工作能量的动能的电子入射在加速器头顶部的靶材料上。生成的韧致辐射光子通过加速器头被跟踪。产生在加速器头内的二次光子和电子也被跟踪。到达加速器头底部的粒子被记录。“记录”包括记录粒子的相空间信息(位置、方向、能量和粒子类型)、它的Monte Carlo重量(用于说明各种使用的物理偏向设计)和关于粒子的其它期望信息。这可以包括：最后的粒子散射位置、散射作用类型、经几次散射作用等等。

这里描述的模拟使用基于加速器制造商提供的机器图和材料数据的Monte Carlo码BEAM96和MCNP4B来运行的。(为了本公开之目的，将补丁用于“发行版”的MCNP4B以含有“BEAM-型”锁存(LATCH)能力及使用BEAM型韧致辐射分裂)。大约5×10⁶到2×10⁷入射电子用于模拟，这取决于哪个加速器在被模仿。(至今我们已经模拟了8台Varian公司的加速器。它们是600C：4MV和6MV，600CD：6MV，2100C：6MV，8MV，10MV，15MV和18MV)。给定使用的方差缩减设计(例如，强制碰撞，粒子分裂)，产生的相空间文件含有用于约3×10⁷个光子信息(可变重量)，每个占有为约1GB到约3GB空间，这取决于用于生成文件的码。

除非另外特别标明，本公开中的图都基于一个用BEAM96计算的6MV Varian 2100C数据设置。

分析与算法的设计

治疗头的Monte Carlo模拟生成的相空间文件的初步分析。

光子相空间文件的初始检查包括画出文件中第一个约10⁵个光子的形成(或最后散射)位置(x_e，y_e，z_e)。这样用于一个格式化(虚拟的)加速器头的一个x_e对z_e散射图示于图2。该步骤在相空间分析中有两个目的—第一个是实用的，第二个是概念的。第一个是对程序的输入的有用检查，因为光子产生的位置应与加速器头的物理结构相关；第二，它显示了治疗仪器头的哪部分有助于仪器的输出。如示于图2中例子，产生自靶的光子(“靶光子”)来自于清晰限定的点。来自主准直仪的光子(“主准直光子”)数量上较少，而且它们倾向于来自准直仪的上边缘。这样，主准直仪的内表面不是均匀的光子源。进一步说，主准直仪更像一个“环”源。矫直过滤器也是一个光子源(“矫直过滤器光子”)。然而，不像主准直仪那样，在那里光子仅仅发生在准直仪很小部分，光子发生于矫直过滤器的全部体内。

在等中心面的能注量分配示于图3。大部分能量直接来自靶，百分之几级的能量来自矫直过滤器和主准直仪。这将治疗仪器头部件的相对重要性做了更加定量的描述。来自每个治疗仪器头部件的能量百分化和包括的物理机制在表I中给出。韧致辐射是靶光子的主要产物机制；非相关散射产生了绝大部分来自于主准直仪和矫直过滤器光子能量。

                      硬件

产生机制   靶+靶回反    主准直仪    矫直过滤器    总计

韧致辐射    0.9045      0.0004       0.0003      0.9052

非相干光    0.0238      0.0247       0.0376      0.0861

散射

相干        0.0030      0.0014       0.0034      0.0079

散射

正电子      0.0004      0.0002       0.0001      0.0007

湮灭

合计        0.9318      0.0266       0.0415       …

                         表I

等中心面的照射区域

下面确立病人体表的区域，病人的体表被来自电子靶，主准直仪和矫直过滤器的光子照射。

靶光子

图3显示靶光子的能注量在距中心轴大约25厘米处非常迅速地减少。靶光子能注量的空间延伸受限于主准直仪的开放角度。距中心轴29厘米以外，靶光子的能注量曲线减少至低于最大值的0.3％。根据上述给定的小于1％的精度要求，距中心轴29厘米是需要考虑的来自韧致辐射靶的光子“能达到”的最大半径。

主准直仪和矫直过滤器光子

图3显示了来自主准直仪和矫直过滤器的光子，在没有准直挡板情况下，照射等中心面非常大的区域。然而，正如下面表明的那样，挡板确定了被这些光子照射的最大区域。

图4A-4C表明，对于一块10厘米×10厘米的开放场，生成于加速器头不同部件(子源)的光子照射不同的但是清楚的等中心面区域。应该理解为可使用许多其它子源(subsource)，这取决于使用的加速器头的类型。如图2所示，靶光子(图4A)源自很小(约2毫米直径)区域，被很好地移除于挡板。主准直仪光子(图4B)源自稍大一点的源，(约1.5厘米直径)的区域与挡板稍微更接近。最后，矫直过滤器光子(图4C)源自一个更大(约3厘米直径)区域，它与挡板相当接近。图4D包括一个简图，该简图表明垂直位置的差异引起在等中心面上不同的但仍清晰的照射区域。这些区域的边缘位置归因于治疗头部件和可移动挡板的相对位置。这些边的位置随着每个挡板移动而移动。上部的挡板比下部挡板给出稍微不同的边位置，这是由于对主准直仪和矫直过滤器挡板位置如图5A-D变化时，边位置的垂直位置变化差异所致。(通过阅读来自该加速器的相空间文件的光子轨迹并看一下该轨迹是否和挡板交叉，便生成了这些图。档板假设是无限衰减的。)挡板的阻挡效果对每一档板位置清晰可见。

辐射跟踪截止

保持前述的将能注量误差保持在最大能注量的大约1％以下的目标时，可大胆断言一个子源仅仅照射那些提供多于最大内场能注量0.2％的区域。这一能量水准定义了场边缘，超出它，来自给定子源的Monte Carlo光子跟踪就不能执行。这个0.2％截止水准用图5A-D每个画面中的破折线表示。所得的场尺寸限制在表II中列出。

档板放置    主准直仪    矫直过滤器

        上     下      上      下

+20    22.6    21.7    30.5    26.6

+15    17.3    16.6    24.1    20.8

+10    12.0    11.4    17.7    15.0

+5     7.0     6.4     11.4    9.2

0      1.7     1.2     5.2     3.2

-5     -3.3    -3.9    -1.5    -2.7

-10    -8.7    …      -7.9        …

                表II

表II和图5显示当档板被移至它们最大范围时，相应一40×40cm²区域，在大于离开束轴的距离 $\sqrt{{22.6}^2+{21.7}^2}\≈31.3\mathrm{cm}$ 地方在等中心面上来自主准直仪的光子将不被跟踪。类似地，在这个同样场尺寸下，在等中心面上距波束中心轴大于 $\sqrt{{30.5}^2+{26.6}^2}\≈40.5\mathrm{cm}$ 的地方来自矫直过滤器的光子将不被跟踪。这定义了最大半径，在这些最大半径上这些部件的光子输出将被特性化。

对源算法应用

表II中列出的用表格表示的照射范围，在源算法中可用来在等中心面上采集x和y。对于一个10cm×10cm对称场的情况(即所有档板移至其各自+5厘米设定处)，主准直仪光子照射于±7.0×±6.4≈180cm²的面积。类似地，对此相同的10×10cm区域，来自矫直过滤器的光子照射一±11.4×±9.2≈420cm²的面积。这个相同的分析可应用于x或y(或二者都)不对称场。确定在一个给定治疗中使用的场大小决定了档板的位置，并且进而-通过图5和表II-这些位置确定用于等中心面处采样光子位置的x与y的最小值和最大值。另外，使用表II中的信息在等中心面上来选择一光子位置(x，y)确保了光子到达病人而不交切挡板的先验的高的概率。

能量分布

各种部件的光子能量分布，如在等角正面和在射束中心轴计数的那样，示于图6A-C。光子能量分布随产生它们的治疗头部件变化很大。来自靶的光子能量分布在初始入射电子能量到很低，但不为零的能量范围内。这与矫直过滤器低能光子的去除是一致的。来自主准直仪的能量分布反映了该过滤处理(在低能量方面)和光子是通过大角度被进行康普顿散射(这样在高能方设置了一个上限)。来自矫直过渡器的光子分布反映了缺少低能量过滤(因为矫直过滤器是光子透过的最后一个治疗仪器头部件)，也反映了来自靶的高能光子的小角散射(因此能量损失很小)的概率。

光子能量分布随距束中心轴的距离而变化。这种分布示于图6D-6F，它们表明了在等中心面和距束中心轴20cm处记录的分布。所有三个子源的分布显示了随着离中心轴距离的增加以低能量光子(这里大约为1Mev)的比例增加，这与这些光子经过的矫直过滤器的厚度的减小相关。

光子初始分布

多个子源，除了照射病人的不同大小面积外，还具有明显不同的源分布。靶光子源是一个平圆盘，主准直仪光子来自一个圆环样的源，矫直过滤器光子来自一个宽阔的中心有尖的源。

虽然这些分布差异很大，每个源可以由辐射分布和角分布精确描述，讨论如下：

靶光子

靶光子来自一个清晰界定的圆盘，圆盘半径与射入电子束半径相匹配。穿过该点的交叉部分示于图7A和7B。该源分布相对平坦。这些射向病人之前在靶内散射的光子在大于R处还有一低密度的尾。(电子束外散射的靶内光子能量低，严格地说，不应与由电子束限定的半径内产生的光子结合。给定尾部小的辐射延伸范围下，在实用中我们把靶光子当成来自一个源)。示于图7A的能注量图被等价地表达为示于图7B径向概率函数。靶光子的方位分布是均匀的。由于这个原因，就产生了均匀分布的(cosΦ，sinΦ)函数对，然后得出：

(X_start，Y_start)＝R_start(cosΦ，sinΦ)。    (1)

散射光子

散射光子以与来自靶的光子相似方式处理，即源位置是从(R_start，Φ_start)对得到的。这里不同的是，Φ_start也是选自一个分布，而不是一致地以间隔(0，2π)选择。

使用角形分布的动机在表1和图8A、8B中给出。在等中心面20上而非在仪器对称轴上任意点，对于发生在散射器“近”边和“远”边的现象散射角有明显的差异，无论是示于图中的矫直过滤器还是主准直仪。计算表明非相干散射对于在主准直仪和矫直过滤器中的光子是主要产生机制(表1)。因为散射光子能量与散射角相反地变化，较多的光子能量将来自治疗头部件的“近”部分。

这将我们引导到用于描述光子轨迹的起始位置的坐标系统上来。如图8C所示，一个极坐标系统可被定义在粒子源面的位置。该系统的角坐标起始点可由含有束轴和交叉在等中心面的光子位置的平面来确定。按照相对方位角_rel而言，这就产生了一个起始位置的简单分解，这里_rel范围从“近”边的0°值到散射器“远”边的180°(假定含有反映对称)。

主准直仪

源对小半径而言是相当各向同性的。该各向同性消失在轴外，伴着这样的分布，在最后板中的分布显示了准直仪角度以外的移动效果：“近”边被阻挡而且只有来自准直管后表面的光子可以抵达等中心面(图9A-9F)。

矫直过滤器

这一同样的效果也出现在矫直过滤器的情况下。这里源分布也倾向于矫直过滤器的“近”边。方位各向异性随着R_iso增加而增加(图10A-10F)。

采样算法

概要

现在话题转到Monte Carlo光子远距放射疗法码的源算法的设计上。“源问题”，如上所述，能够以既精确又高效的方式生成光子。在上述的讨论过程中，可看到出现两个主模式。第一个趋势是能谱、辐射能注量变化和起始分布对不同子源差异很大。第二，对一个给定子源，能量分布、辐射能注量和起始分布作为等中心面上光子辐射位置的函数全都平滑且相对缓慢地变化。

把全部光子相空间分为单个的相空间，每个由产生光子的治疗头部件确定。接着将每个子源的相空间细分为“片”：在等中心面上确定的等宽环状区域。在上面已表示出，给定准直挡板移动的最大范围——光子有相当大着陆概率的位置距中心轴的最大距离，对靶光子约29厘米，对主准直光子约31厘米，对矫直过滤器光子约41厘米。现在对每个子源十个片被使用。每个片有其自己的光能谱；子源的辐射缓慢变化的能谱表现为依赖于片的能量谱。除了它的能量谱之外，每个片有其自己的R_start和_rel分布以描述光子起始位置。

实现

用于描述光子相空间的相关矩形图示于图11A-11F。考虑到治疗问题对一部机器含有一个场的尺寸；对多重场尺寸的概括是直接的。在计算治疗方案中相关矩形图的使用除在运行时间执行的步骤之外，需要几个预处理步骤：

(1)预处理步骤；

(a)为计算中使用的机器读数据文件；

(b)把每个能量、角度和辐射概率矩形图转换为假采样要求的形式(见附录)；

(c)使用在治疗方案中确定的档板位置来确定每束的每个子源的(x_iso，y_iso)采样区间的极限；

(d)对每个子源，用(x_iso，y_iso)采样区间来计算传给等中心面的能量；并且

(e)用相对能量贡献作为每个子源的采样概率。

(2)运行期产生光子

(a)用随机采样来确定将产生哪种类型的光子(如：靶、主准直仪或矫直过滤器)。如表1所指出的，这些光子子源将按大约93∶3∶4的相对概率被选出。实际的相对概率依赖于特定治疗方案的档板设定。

(b)对该光子类型，给定档板设定，查找等中心面上照射的x和y极限(图4)。在该区域内产生一个随机的、均匀分布(x_iso，y_iso)的坐标对。

(c)给出该(x_iso，y_iso)对，计算R_iso。调整粒子的重量以说明该光子类型的慢变化的能注量(图3)。给出R_iso，确定片指数，计算Φ_iso。

(d)基于片指数，从适与这个R_iso(图6)的能量分布采样粒子能量。

(e)从适当的分布采样起始半径R_start；

(f)采样相对起始角_rel。对靶光子，使用方程(1)。对其它光子，采样适当的分布。

这样粒子的能量和重量就已确定了(步骤c和d)，要有两个点确定粒子的轨迹(步骤b、e和f)。确定轨迹的点确定了粒子的方向余弦，用于开始跟踪该粒子所需的相空间信息现在完整了。光子被反向跟踪至监视室到底，然后穿过波束调节仪(有的话)传到病人。

精度

下一个问题是上面描述算法的精度。这里定义的精度指我们的算法再现潜在的光子相空间的程度。这是一个双部分的问题。第一个问题处理在等中心面上光子的能量分布。第二个问题处理光子轨迹。

能量在等中心面上的分布

这些分布被精确地再现，因为(a)在该面上实际地完成能量分布的分箱(binning)，(b)由于该分布随轴距变化足够慢(图12)，以致分段线性分布能够跟踪这种变化(图13)。

这种变化随辐射距离的变化很慢。图12A-12C表明几个能量的箱概率(bin probability)曲线作为等中心面处半径的函数而绘制。这些曲线都缓慢变化，而且都精确地按分段线性函数表示。另外，仅有很小块数片(大约是10)就足够了。

因为图12中示出的分段线性函数的结点实际上是矩形图的箱中心(bin-center)，矩形图是通过把光子能量分箱(binning)到该范围上得到的(图13)，通过在等中心面上，在少数的环状区域内作出能量矩形图，在所关心的区域上可再现能量分布的辐射变化。然后示于图6的该形状光子能量分布就在每个片内产生。最大片半径由所关心的用于特定光子分布的最大半径确定。对来自韧致辐射靶的光子而言，最大半径约为29厘米。类似地，如上所提到的，对主准直器和矫直过滤器，所要求的最大片半径为33.6和40厘米。

轨迹信息

第二个问题处理光子的轨迹。我们的源算法通过“连接点”确定光子的轨迹，其中一点在等中心面上，其它点在韧致辐射靶位置上或几个厘米以内。源分布的辐射范围对靶光子来说大约是1毫米级，对从主准直器或矫直过滤器散射的光子约为1厘米至3厘米。这样期望的轨迹误差最大值预期不大于δθ≈arctan3/100≈2度。实际误差将更小，因为使用了起始R和Φ分布。

轨迹精度可通过对每个源将作为挡板位置函数的挡板边缘半阴影(penumbra)的位置和宽度和用潜在的相空间文件得到的位置和宽度比较而检查。(“刀刃”测试)。当随机选择等中心面上的光子位置时，在光子方向余弦中的“误差”将导致假想轨迹中的误差。例如，这个轨迹误差会影响，靠近治疗场边缘的光子是否被挡板边缘阻挡。对矫直过滤器光子，这种类型测试示于图14，同时上部挡板移动到距光束中心线最大距离处。能注量降落由本文中提出的源算法很好地匹配，表明使用相关R和分布精确地再现光子轨迹。图14也示出了流量误差，这是当被假设为均匀分布时引起的。

计算效率

该算法的计算效率有两个基础。第一个是概念结构，其中我们首先只选那些进入用于被特定挡板照射的等中心面部分的光子轨迹。只有这样我们才产生一个起始位置来确定“它们从哪来。”第二个基础是“假方法”(见下)的有效性，去采样用于描述光子相空间的各种概率分布。这种方法需要用两个随机数以从任意长度和复杂性的分布中选择一个箱(bin)。如果期望在箱内均匀采样的话，需要第三个随机数。

该算法的计算要点在表III中详列出。注意源算法确定了六个必需的粒子相空间信息(x，y，u，v，E，q，重量)，用了大约14个随机数和3个或4个表查找。(回想一下，给定问题的上下条件，z是等中心面的位置，

步骤1，使用假采样来确定选择机器的那个治疗头部件，使得该算法的运行时间独立于治疗方案中使用的射束数。步骤数#            操作                            所需RNS数#   采样的矩形数#1     选择哪一个光束的哪一个部件                       2           1^A2     在等中心面上选择(I_iso，y_iso)，调整光子重量      2           03     为能量分布确定名称记号                           1^B         04     确定能量                                         3^C         15     确定起始半径R_start                               3^C         16     定起始角Φ_start                                  2或4^D      0或1总计                                                     13或15      3或4

表III

步骤2和6对算法效率影响最大。通过步骤2、4和5计算出的许多光子轨迹将交于挡板的拐角，因此不会到病人处。这种行为是想要的，因为这就是半阴影是如何确定的。该效果对小场尺寸是最明显的，并且随着治疗头器件辐射范围增加该效果变得更明显。对很多方形场大小表示于图15。对大多数场(field)，“靶”光子轨迹通过档板的机会大于90％。来自主准直器和矫直过滤器的光子生存概率较低。然而由于大部分光子起源于靶内，在研究过的范围内大部分场总效率大于80％。

相空间文件的压缩

如上所述一个自动化程序对于把从加速器头的Monte Carlo模拟获得的千兆字节大小的数据文件转换为这里描述的相关矩形图是有用的。转换需要五遍通过1GB相空间文件，而且需要大约23分钟完成。(用的计算机是Digital AlphaServer 8400 Model 5/440，运行OSF1 V4.0软件)。

五遍是执行下述操作：

(1)该步骤计算示于图5的曲线。该步骤也计算表II中给出的数据。在这样做过程中就也计算了每个子源需要的最大片半径和由此得到的每个子源的每个片的辐射极限。

(a)对相空间文件中的每个粒子：

(i)读每个粒子并辩认在其中最后散射的硬件部件。

(ii)基于在计数平面(x_tally，y_tally，z_tally)的粒子位置和粒子的方向余弦(u，v，w)，确定粒子在何处和等中心面相交。同时，确定表II中的档板设置，在该设置粒子能碰不到挡板。把该粒子加入示于图5中的曲线。

(b)在此处，已计算出图5中的曲线。接着查找每一曲线来确定曲线位置，在该位置曲线高度下降低于在等中心面上最大能注量的0.2％(标称值)。

(c)该交叉点位置进入到表II中。

(d)根据挡板位置的最大允许量，为每个子源计算出R_max，这样就确定了每个子源的每个片的辐射范围。

(2)该步骤对每个子源的每个片的“源面”(示于图8)位置做了一个初始估计。回想一下，通过首先选择在等中心面的一个位置然后询问“粒子从哪来？”，源粒子算法计算粒子的轨迹。因为确定粒子起始位置时有误差而出现了轨迹误差。将辐射分布范围最小化可以反过来使运行过程中计算出的轨迹误差最小化。为此目的，确定粒子的z值，在该值上每个粒子最靠近束中心轴。

(a)对每个粒子：

(i)确定其子源和其片指数；

(ii)当粒子交于计数平面时，使用粒子位置和方向；向后跟踪到粒子最靠近束中心轴的路径 $Z_{\mathrm{closest}}=Z_{\mathrm{tally}}-\frac{w({\mathrm{ux}}_{\mathrm{tally}}+{\mathrm{uy}}_{\mathrm{tally}})}{1-w^2}$

(iii)根据z_closest值，对这子源的这个片把粒子加到N(z)矩形图适当的箱。

(b)对每个子源的每个片，设置z_guess等于N(z)矩形图有最大值的z值。

(3)此步为每个子源的每个片确定高端截止R_start值。

(a)对每个粒子：

(i)确定其子源和该片指数；

(ii)在上面确定的z＝z_guess的位置，确定其距束中心轴的距离R；

(iii)对适当的片和子源将粒子加给R_start矩形图的适当的箱。

(b)对每个子源的每个片，找出超出99％的粒子的R_start的最小R值。

(4)给每个子源的每个片微调z_start值。上述步骤2确定的z_guess值产生了一个好但非最优的在该值产生R_start和_rel分布的z值。一个方位上各向同性的起始分布，在图16中举例示出了主准直器的例子，当属于该分布的光子被向后跟踪到z_guess时出现了“偏心”。这一图案可用于对最大化(R，)表示的精度的ΔR值的寻找。

(a)对每个粒子：

(i)确定其子源、片指数和z_guess；

(ii)计算其在z_guess的(x，y)位置；

(iii)对适当的片和子源，将粒子加到二维矩形图的适当箱。

(b)对每个这样计算出的二维矩形图，沿着x轴并在每点：

(i)假设该点为起始点，产生一个上面计算出的(x，y)分布的(R，)表示；

(ii)找出x值，当其用作起点时，能将上面的(x，y)分布和从基于起始点的(R，)表示算出的分布的差别最小化。

(c)确定不需要任何偏移的z_start。

(5)此步骤产生了用于描述加速器输出的相关矩形图。对每个粒子：

(a)确定其子源、片指数和Φ_iso；

(b)确定该片在z_best的R_start和_rel，并更新适当的矩形图；

(c)更新该片的能量分布。

该假采样方法对包括很多数目的箱的采样概率分布是一有效方法，这些箱依次有大宽度变化范围的高度。并且，虽然假采样对采样分布是一种很好的方法，但它不是可得到的唯一方法(其它方法更易编码，但会更慢)。这里展示的分布一般同时有这两个特点。业已发现该假方法比其它采样分布方法更快，例如，搜索累计分布方法，同时避免与使用同等概率分布有关的精度损失问题。

假采样方法是一个两阶过程。第一阶段做为预处理步骤，且是将示于图17A的初始分布再映射到示于图17B的表格形式。第二阶段发生在Monte Carlo计算运行期间。产生了两个随机数ξ1和ξ2。ξ1用来选取图表的一栏。ξ2用于与该栏的截止级(cutoff-level)相比较，并且根据结果，或者低的或者高的箱指数被返回。这样，对显示的例子，一对数(ξ1，ξ2)＝(0.5，0.9)表明图17a中的“箱1”被选中。类似地，数对(ξ1，ξ2)＝(0.9，0.07)表明图17A的“箱5”应被选中。

如图17B所看到的，概率分布的图表形式需要三个阵列：一阵列用于箱边，一阵列容纳截止级，一列容纳“上部箱”值。后面两个列阵是在假方法的预处理阶段得到的。预处理阶段是一个迭代过程，在此过程中“高度”从“高”箱转移到“低”箱。步骤描述如下并概括在表IV中。

	箱内容					被选箱		截止级					上部箱指数
																		步骤	1	2	3	4	5	HIGH	LOW	1	2	3	4	5	1	2	3	4	5
0	2.5	1.25	0.75	0.4	0.1	1	3	…	…	…	…	…	…	…	…	…	…
																		1	2.25	1.25	1.0	0.4	0.1	1	4	…	…	0.75	…	…	…	…	1	…	…
2	1.65	1.25	1.0	1.0	0.1	1	5	…	…	0.75	0.4	…	…	…	1	1	…
																		3	0.75	1.25	1.0	1.0	1.0	2	1	…	…	0.75	0.4	0.1	…	…	1	1	1
4	1.0	1.0	1.0	1.0	1.0	…	…	0.75	1.0	0.75	0.4	0.1	2	…A	1	1	1

            表IV

1.找出最高的现有箱，和高度低于平均值的第一个箱。对示于图17A中的分布，对于第一个迭代它们分别为箱1和3。

2.加入来自“高”箱中足够的“高度”使“低”箱上升到平均水平。现已定义了用于“低”箱的反-假指数(anti-alias index)。

3.以在步骤2中使用的量减少“高”箱的“高度”。

4.到步骤1，重复直到所有箱要么被“上升”要么“降低”到平均“高度”，视情况所需。

尽管在该部分中使用的例子为均匀箱宽的矩形图分布，假采样方法已被扩展到有非均匀间隔横座标的分段线性函数。

为了说明和描述，已给出了本发明的前面描述，而无意穷举或限制本发明到公开了的精确形式。按上面的教导下许多调整和变化是可能的。选择和描述实施例来最佳解释本发明原理及其实际应用，这样就可使其他本领域技术人员在各实施方案下和以适合所考虑的满足特定使用的各种修改情况下最好地使用本发明。发明的范围将由权力要求书确定。

Claims (20)

1.一种计算加速器治疗头相空间信息的方法，包括：

获取特定辐射加速器头的输出数据文件；

将上述输出数据文件转换成相关矩形图，在剂量计算期间该矩形图可被有效地采样以产生精确的粒子相空间信息，其中所述矩形图至少描述一个子源；并且

对来自所述特定的辐射加速器头的具有通过由一套档板确定的开口的高的概率的粒子，从上述的相关矩形图，计算用于所述的特定辐射加速器头的所述相空间信息。

2.权利要求1的方法，其中所述输出数据文件是采用所述特定辐射加速器治疗头的辐射输出的Monte Carlo模拟得到的。

3.权利要求1的方法，其中将所述输出数据文件转换成相关矩形图这一步包括：

对每一个所述至少一个子源计算最大片(tile)半径；

对每一个所述至少一个子源的每一片计算半径极限；

对每一所述至少一子源的所述的片，估算源平面的位置；

对每一所述至少一子源的所述片，确定R_start的高端截止值；

对每一所述至少一子源的每一所述片，微调Z_start；

用所述最大片半径、所述半径极限、所述源面、所述R_start的高端截止值以及Z_start值，对每个所述至少一个子源计算所述相关矩形图。

4.权利要求3中方法，其中将所述输出数据文件转换成相关矩形图的这一步需要5遍通过所述输出数据文件。

5.权利要求1中方法，其中所述相关矩形图包括来自每一所述至少一个子源的粒子的能量、角度以及径向的概率信息，其中将输出数据文件转换成相关矩形图这一步包括将所述粒子能量，角度，径向的概率信息分箱到矩形图中的步骤。

6.权利要求1中方法，进一步包括确定来自每一所述子源的将到达等中心面的粒子的能量或数量。

7.权利要求6的方法，其中确定来自所述至少一子源的将到达等中心面的粒子的能量或数量这一步包括：

使用在治疗方案中规定的挡板位置，对于每个所述至少一子源确定(X_iso，Y_iso)采样区域的极限；

对每一子源由所述(X_iso，Y_iso)采样区域极限计算(x_iso，y_iso)采样面积；

从所述(x_iso，y_iso)采样面积，计算从每一所述至少一个子源传送至等中心面的粒子的相对能量或数量。

8.权利要求7的方法，进一步包括使用传送至等中心面的所述的相对能量，计算传送至等中心面的粒子的绝对能量或数量。

9.权利要求8的方法，进一步包括对多个子源重复计算传至等中心面的粒子的绝对能量或数目这一步骤。

10.权利要求8的方法，进一步包括对于尽可能识别的许多子源重复计算传送至等中心面的粒子的绝对能量或数量的步骤；

11.权利要求1的方法，其中所述至少一个子光源选自含有靶，准直器和矫直过滤器的组。

12.权利要求1的方法，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步包括，对于每一所述至少一个子源，用所述的(x_iso，y_iso)采样区域极限，产生粒子轨迹和等中心面交叉的(X₁Y)位置。

13.权利要求1的方法，其中从所述的相关矩形图为所述的特定幅射加速器头计算相空间信息这一步包括，从一个粒子的(x，y)位置，确定等中心面上的距离R_iso和角Φ_iso。

14.权利要求1的方法，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步包括调整粒子的Moute Carlo重量以解决基于R_iso的等中心面上非均衡能注量或子源的数目流(number fluence)。

15.权利要求1的方法，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步，包括使用R_iso选择一套相关矩形图，该套相关矩形图描述适合于来自特定子源的，以此特定R_iso入射到等中心面上的一个粒子的能量、角度和径向概率信息。

16.权利要求1的方法，其中所述相关的矩形图被采样以得到一个粒子的能量，R_start和_rel。

17.权利要求1的方法，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步包括使用R_start、_rel和Φ_iso以对一个粒子计算在位置Z_start的X_start和Y_start。

18.权利要求1的方法，其中一个“隼(PEREGRINE)”规定面是“隼(PEREGRINE)”开始追踪穿过空气射向病人的粒子的平面，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步包括使用用于特定加速器的X_start、Y_start、Z_start、x_iso、y_iso及已知的z_iso值以确定粒子的轨迹和所述的隼(PEREGRINE)规定面的交叉。

19.权利要求1的方法，其中从所述的相关矩形图为所述的特定辐射加速器头计算相空间信息这一步包括，使用用于特定加速器头的X_start、Y_start、Z_start、x_iso、y_iso及已知的z_iso值以在“隼”规定面上所述的粒子轨迹位置处确定粒子的轨迹的方向余弦。

20.权利要求1的方法，其中计算相空间信息这一步包括以下步骤：

为给定一套由治疗方案特定的档板位置，确定来自每个所述至少一个子源将到达等中心面的粒子的能量的大小或粒子的数量。

使用所述粒子的能量或数量，以确定一相对概率，用于采样每一所述至少一个子源描述所述辐射加速器头。

使用所述的相对概率从每一所述至少一个子源采样。

使用所述治疗方案中特定的档板位置，为每一所述至少一个子源确定(x_iso，y_iso)采样区域的极限；

为所述至少一个子源从所述的(x_iso，y_iso)采样区域的极限，计算在所述等中心面上的(x_iso，y_iso)位置；

从所述的等中心面从所述的(x_iso，y_iso)位置计算距离R_iso和Φ_iso；

基于R_iso调整粒子的Monte Carlo重量以解决在所述等中心面上每个至少一个子源的非均衡能注量或数目流；

使用特定的R_iso选择一套相关矩形图，该套相关矩形图描述适合于来自所述至少一个子源的特定子源，以特定R_iso入射到等中心面上的一个粒子的能量、角度或径向概率信息；

采样所述的一套相关矩形图以得到粒子的能量，R_start和_rel；

在位置Z_start处，用R_start、_rel和Φ_iso为所述粒子计算X_start和Y_start；

用所述的用于特定加速器的X_start、Y_start、Z_start、x_iso、y_iso及已知的z_iso值确定粒子轨迹与“隼”规定面的交点；

使用用于特定加速器的X_start、Y_start、Z_start、x_iso、y_iso及已知的z_iso值确定粒子轨迹在所述的“隼”规定面上的所述的粒子轨迹位置上的方向余弦。

Images