CN106354896A - 确定经过射束调整装置的粒子强度分布的方法和装置 - Google Patents

Abstract

本发明提供了一种确定经过射束调整装置的粒子强度分布的方法和装置。该方法包括以下步骤：计算粒子经过射束调整装置后的强度分布，包括：获取粒子穿过射束调整装置前的状态，包括粒子在射束坐标系中的权重和/或能量、位置与方向；使用该射束调整装置的各限束部件的交线长度矩阵和粒子在射束坐标系中的状态来计算粒子经过射束调整装置的交线长度，其中所述交线长度矩阵包括不同入射位置的粒子与不同入射方向的粒子与该射束调整装置的对应限束部件的交线长度；使用粒子穿过射束调整装置前的权重和/或能量以及粒子穿过射束调整装置的交线长度来计算粒子经过射束调整装置后的权重和/或能量改变；以及将粒子的位置沿着粒子的方向由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

Description

确定经过射束调整装置的粒子强度分布的方法和装置

技术领域

本发明主要涉及粒子输运领域，尤其涉及一种确定经过射束调整装置的粒子强度分布的方法和装置。

背景技术

粒子经过具有限束功能的结构之后，其强度如何分布，是粒子输运领域关注的课题。粒子输运可应用在核物理实验、反应堆设计和辐射屏蔽等场合。在医疗领域中，粒子输运问题是当前放射治疗的一个研究热点，基于此可以模拟模体内的剂量分布。

具体来说，在放射治疗中常用到多叶准直器(Multi-Leaf Collimator,MLC，又称多叶光栅)，在不同机架角度形成一个或多个开口，并且多个开口在病人体内产生的剂量叠加，从而使得不同角度的开口形状与病人靶区形状贴合，且其强度分布满足靶区不同位置处不同剂量大小的要求，从而达到更好的治疗效果。

剂量计算算法根据上述信息以及其他信息，计算出经过MLC及其他装置后的射束粒子的能量、位置、权重与方向分布，并且利用病人或模体CT数据，计算出病人或模体接受到的剂量分布。算法的计算速度和计算精度具有重要的临床意义。

对于粒子输运来说，射束调整装置的形状决定了强度分布计算的速度和精度。在放射治疗设备中，由于MLC形状复杂，叶片数目繁多，对剂量计算的速度和精度造成很大影响。目前主要的经过MLC后的强度分布算法如下：

第一种方法是，首先得到所有叶片的端面(叶片闭合时靠近射束中心轴的一端)位置，将所有叶片端面所在的坐标点相连形成一个闭合的轮廓线，将位于轮廓线内的射束强度赋一个较大值，轮廓线外的射束强度赋一个较小值，而轮廓线附近的射束强度置为上述较大值和较小值之间的数。该方法得到强度分布速度较快，但不能准确考虑端面效应、叶间漏射与凹凸效应。

第二种方法为，将MLC的叶片位置与材料信息及其他部件信息建模，通过算法抽取一个具有能量、位置、方向的粒子，模拟该粒子在MLC及其他部件内可能发生的各种物理反应，从而改变方向、位置与能量。从多叶准直器穿出的粒子的能量、位置、方向即形成强度分布，用于进一步模拟粒子在病人或模体中的剂量分布。该方法可以准确地考虑端面效应、叶间漏射与凹凸效应，但由于每个粒子都要实时模拟各种物理反应，因此耗时非常长，不能用于快速剂量计算。

第三种方法为，将每个叶片用方程表示，在剂量计算过程中，通过解方程或者不等式的方法求得每个粒子与每个叶片的交线长度，利用一定的公式对粒子状态直接作出改变。该方法可以准确地考虑端面效应、叶间漏射与凹凸效应，但是在计算过程中，每个粒子都需要进行解方程或解不等式的运算，耗时较长，远远不能满足快速剂量计算的要求。

发明内容

本发明要解决的技术问题是提供确定经过射束调整装置的粒子强度分布的方法和设备，能够提高强度计算的精度和速度。

为解决上述技术问题，本发明提供了一种确定经过射束调整装置的粒子强度分布的方法，包括以下步骤：计算粒子经过射束调整装置后的强度分布，包括：获取粒子穿过该射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向；使用该射束调整装置的各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度，其中所述交线长度矩阵包括不同入射位置的粒子与不同入射方向的粒子与对应的限束部件的交线长度；使用粒子穿过该射束调整装置前的权重和/或能量以及粒子穿过该射束调整装置的交线长度来计算粒子经过该射束调整装置后的权重和/或能量改变；以及将粒子的位置沿着粒子的方向由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

可选地，在获取粒子穿过该射束调整装置前的状态之前，还包括预处理，所述预处理包括产生不同入射位置的粒子与不同入射方向的粒子与该射束调整装置的各限束部件的交线长度矩阵。

可选地，使用该射束调整装置的各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度包括：首先计算粒子与射束调整装置的各限束部件的交线长度，再求和获得粒子经过该射束调整装置的交线长度。

可选地，使用交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度包括：利用粒子在该射束坐标系中的位置与方向，及各限束部件坐标系相对于该射束坐标系的平移或平移与旋转，求得粒子在各限束部件坐标系中的位置与方向；利用各限束部件对应的交线长度矩阵，获取与粒子在各限束部件坐标系中的位置与方向所对应的交线长度，分别作为粒子与各限束部件的交线长度。

可选地，该预处理还包括根据各限束部件配置信息计算出各限束部件坐标系相对于该射束坐标系的平移或平移与旋转。

可选地，在使用交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度的步骤中，在求得粒子在各限束部件坐标系中的位置与方向之前还包括：判断粒子是否穿过且只穿过限束部件部分高度，如果是则在使用该交线长度矩阵的条件下求得粒子在各限束部件中的交线长度，否则在不使用该交线长度矩阵的条件下求得粒子在限束部件中的交线长度。

可选地，该预处理还包括产生该射束调整装置的各限束部件的散射强度分布。

可选地，该计算粒子经过射束调整装置后的强度分布还包括利用各限束部件坐标系相对该射束坐标系的平移或平移与旋转，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

可选地，产生该射束调整装置的各限束部件的散射强度分布包括：将该射束调整装置的各限束部件的形状和材料在蒙特卡罗程序中进行建模，将经过该射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序，模拟得到在每个限束部件中产生的次级粒子的能量、位置、权重与方向，作为各限束部件的散射强度分布。

可选的，在将经过该射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序之前，对于无权重的粒子，定义其权重为任意正的常数且各粒子的权重相等，对于无能量的粒子，定义其能量为任意正的常数且各粒子的能量相等。

可选地，该计算粒子经过射束调整装置后的强度分布还包括根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子。

可选地，该根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子包括：设置阈值W0和存活概率P0，若粒子权重小于W0，则在0到1之间抽样得到一个随机数，若该随机数小于P0，则判定该粒子存活，将其权重除以P0；若该随机数大于或等于P0，则将该粒子丢弃。

可选地，该射束调整装置为多叶准直器，该限束部件为多叶准直器的叶片。

本发明还提出一种确定经过射束调整装置的粒子强度分布的装置，包括预处理模块，配置为产生不同入射位置与方向的粒子与该射束调整装置的各限束部件的交线长度矩阵，其中各限束部件确定限定该射束调整装置的形状；

强度计算模块，耦接该预处理模块，配置为按照下述步骤计算粒子经过该射束调整装置后的强度分布：获取粒子穿过该射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向；使用该射束调整装置的各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子与各限束部件的交线长度，并求和获得粒子经过该射束调整装置的交线长度，其中所述交线长度矩阵包括不同入射位置的粒子与不同入射方向的粒子与对应的各限束部件的交线长度；使用粒子穿过该射束调整装置前的权重和/或能量以及粒子穿过该射束调整装置的交线长度来计算粒子经过该射束调整装置后的权重和/或能量改变；以及将粒子的位置由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

可选地，该预处理模块还根据各限束部件配置信息计算出各限束部件坐标系相对于该射束坐标系的平移或平移与旋转。

可选地，该预处理模块还产生各限束部件的散射强度分布。

可选地，该强度计算模块还利用各限束部件坐标系相对该射束坐标系的平移或平移与旋转，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

可选地，该强度计算模块还根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子。

与现有技术相比，本发明的上述技术方案具有以下优点：

(1)直接获取粒子与限束部件交线长度矩阵，只需要进行坐标变换和叠加就可以求得粒子穿过射束调整装置之后的强度分布，相比实时计算交线长度的方法，本发明为强度计算节省大量时间。

(2)在本发明中，限束部件的端面及侧面形状不受限制，其对粒子的影响结果反映在交线长度矩阵中，因此可以在不增加存储空间和计算时间的情况下，充分考虑端面效应、叶间漏射以及凹凸效应，使得计算的强度分布更加精确。

附图说明

图1是本发明一实施例的确定经过射束调整装置的粒子强度分布的装置结构框图。

图2是本发明一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。

图3是本发明另一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。

图4是本发明另一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。

图5是根据本发明一实施例的多叶准直器与射束示意图。

图6是根据本发明一实施例的叶片坐标系、射束坐标系与粒子示意图。

具体实施方式

为让本发明的上述目的、特征和优点能更明显易懂，以下结合附图对本发明的具体实施方式作详细说明。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其它不同于在此描述的其它方式来实施，因此本发明不受下面公开的具体实施例的限制。

在本发明的上下文中，射束调整装置是通过对射束的衰减作用来实现对射束状态的调节。射束调整装置可包括一个或多个用于限定射束形状的限束部件。这些限束部件的组合形成一个或多个形状被设计的开口，开口处允许粒子通过。此处，所述粒子可以指实际的粒子，例如光子、正电子、负电子、中子、重离子或其它粒子中得至少一种，也可以指确定位置和方向的具有一定宽度或无限窄的射束，例如笔形束，但是笔形束的位置和方向可能是一个很窄的范围，而不是固定值，在这种情况下，笔形束的位置和方向由射束中心轴决定。不通过开口而是穿过限束部件的粒子，其状态(例如能量，权重)会发生改变，从而影响经过射束调整装置后的强度分布。本发明的实施例描述确定经过射束调整装置的粒子强度的方法和装置。

图1是本发明一实施例的确定经过射束调整装置的粒子强度分布的装置结构框图。如图1所示，装置100可以例如使用通用计算机设备作为基础硬件。在一实施例中，可以通过执行安装在计算机设备中的处理器上的运算程序来实现所需的处理和计算。这时，可以通过预先将运算程序安装到计算机设备中或者通过将运算程序记录在磁盘、磁光盘、光盘、半导体存储器等中或在网络上发布运算程序并将运算程序安装到计算机设备中来实现装置100的各模块。在另一实施例中，可以通过诸如逻辑电路的硬件来实现装置100的各模块。在又一实施例中，也可以通过组合硬件和软件来实现装置100的各模块。

装置100可包括预处理模块110和强度分布计算模块120。预处理模块110输入射束调整装置的材料、形状等信息。预处理模块110配置为生成强度分布计算所需的数据。例如，预处理模块110可配置为求得以下数据中的一项或多项：交线长度矩阵、限束部件散射分布和限束部件在射束坐标系中的位置。在各个实施例中，将根据具体方法的不同选择产生上述的一项或多项数据。强度分布计算模块120耦接预处理模块110以输入前述的数据。强度分布计算模块120配置为进行强度计算，输出粒子经过射束调整装置对应的强度分布。例如，强度分布计算模块120可配置为计算原粒子强度和可选的散射强度分布，其中原粒子是相对次级粒子而言，是指在经过射束调整装置之前就已经存在的粒子。

图2是本发明一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。参考图2所示，方法可包括以下步骤：

步骤210，预处理。

这一步骤可由图1所示的预处理模块110完成。在此实施例中，这一步骤可包括：

步骤211，产生不同入射位置的粒子与不同入射方向的粒子与射束调整装置的各限束部件的交线长度矩阵。在此，不限定入射到射束调整装置的粒子的产生方式。此外，形状相同或相似的限束部件可以共用一个交线长度矩阵，而不必重复计算。

步骤212，根据各限束部件配置信息计算出各限束部件坐标系相对于射束坐标系的平移或平移与旋转。

在此实施方式中，需要计算各限束部件坐标系相对于射束坐标系的位置变化。通常情况下，为了形成开口，各限束部件会进行平移移动，因此，需要计算出各限束部件坐标系相对于射束坐标系的平移；有时，由于各限束部件的安装方位不同或其它原因，有的限束部件相对于射束坐标系还会存在一定的旋转，此时需要计算出各限束部件坐标系相对于射束坐标系的平移与旋转。

预处理步骤可以在强度计算之前进行从而不占用强度计算时间，预处理产生的数据可先行储存，储存的数据只需占用少量存储空间。在之后的强度计算时，只要射束调整装置没有变化，则预处理的步骤不需要重复进行，直接调用储存的数据即可。

步骤220，计算粒子经过射束调整装置后的强度分布。

在此实施例中，计算原粒子强度分布作为强度分布。具体步骤包括：

步骤221，获取粒子穿过射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向。

步骤222，使用各限束部件的交线长度矩阵和粒子在射束坐标系中的位置与方向来计算粒子与各限束部件的交线长度，并求和获得粒子经过射束调整装置的交线长度。

步骤223，使用粒子穿过射束调整装置前的权重和/或能量以及粒子穿过射束调整装置的交线长度来计算粒子经过射束调整装置后的权重和/或能量改变。

步骤224，将粒子的位置由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

在上述的步骤222中，可以直接使用预处理中预先获得的各限束部件的交线长度矩阵来计算交线长度，因此有利于缩短计算时间。另外，获得粒子经过射束调整装置的交线长度的方式并不限于前述的通过各限束部件的交线长度矩阵来求和计算，也可以是其它计算方式。

进一步，步骤222计算粒子与各限束部件的交线长度的过程可包括：首先利用粒子在射束坐标系中的位置与方向，及各限束部件坐标系相对于射束坐标系的平移和旋转，求得粒子在各限束部件坐标系中的位置与方向；然后利用各限束部件对应的交线长度矩阵，获取与粒子相对于各限束部件坐标系的位置与方向对应的交线长度，分别作为粒子与各限束部件的交线长度。

在一实施例中，可以通过插值方式，利用各限束部件对应的交线长度矩阵来获得粒子在限束部件中的交线长度。

在另一实施例中，也可以先对粒子做判断，若粒子穿过且只穿过限束部件部分高度，则按照前述方式获得粒子在限束部件中的交线长度，否则(即完全不穿过限束部件或穿过限束部件的完整高度)可以按照前述方式通过交线长度矩阵获得粒子在限束部件中的交线长度，也可以直接或经过简单计算得到粒子在限束部件中的交线长度，而不需要做前述的坐标变换和从交线长度矩阵中取值。例如，若粒子完全不穿过限束部件，则该粒子与限束部件的交线长度显然为0；若粒子穿过限束部件的完整高度，则该粒子与限束部件的交线长度为h/cos(f)，其中h为限束部件沿Z轴方向的高度，f为粒子的方向矢量V相对于Z轴正方向旋转的角度。

在上述的步骤223中，粒子状态的改变示例性计算方法如下：

若粒子为光子，则其状态的改变在于粒子的权重下降，下降量与穿过的材料厚度及在材料中的线衰减系数有关。穿过射束调整装置后的粒子权重表示为：

$W2=W1\×e^{-{\mathrm{\Σ}}_{i=1}^N\left({\mathrm{\μ}}_i{\mathrm{\Σ}}_{j=1}^M{l}_{ij}\right)}---\left(1\right)$

其中W2为粒子改变后的权重，W1为粒子穿过射束调整装置之前的权重，N为射束调整装置个数，M为第i个射束调整装置中的限束部件数目，μ_i为粒子在第i个射束调整装置的平均线衰减系数，l_ij为粒子在第i个射束调整装置第j个限束部件中的交线长度。

若粒子为电子或其他粒子，则其状态的改变在于粒子的能量下降，下降量与穿过的材料厚度及在材料中的阻止本领有关。穿过射束调整装置后的粒子能量表示为：

$E2=E1-{\mathrm{\&Sigma;}}_{i=1}^N\left(S_i{\mathrm{\&Sigma;}}_{j=1}^M{l}_{ij}\right),$ 若E2<0则E2＝0 (2)

其中E2为粒子改变后的能量，E1为粒子穿过射束调整装置之前的能量，N为射束调整装置个数，M为第i个射束调整装置中的限束部件数目，S_i为粒子在第i个射束调整装置的平均线阻止本领，l_ij为粒子在第i个射束调整装置第j个限束部件中的交线长度。

可以理解，由于射束调整装置可具有一个或多个开口。对于射束调整装置的每个开口均需要执行上述的步骤221至224。对于不同开口来说，同一个限束部件在不同开口中具有不同的位置，也就是说在射束坐标系下的一个粒子，变换到两个开口情况下的同一个限束部件中时具有不同的权重和/或能量、位置与方向，可以分别求得各自的权重和/或能量、位置与方向对应的交线长度即可。

对于多开口的射束调整装置。以光子为例，可利用下列公式得到合成的粒子权重。

$W={\mathrm{\&Sigma;}}_{k=1}^L({\mathrm{MU}}_k\&CenterDot;W{2}_k)---\left(3\right)$

其中W为多个开口合成的粒子穿过射束调整装置后的权重，L为开口数目，MU_k为第k个开口的机器跳数，W2_k为第k个开口求得的粒子的权重。

对于电子，分别求得每个开口形成的强度分布。

在上述的步骤224中，形状及材料相同或相似的叶片可以共用一个原粒子强度分布，而不必重复计算。

对于只修改权重的粒子如光子而言，对于多开口射束调整装置来说有利的是，每增加一个开口，只需增加一次坐标变换，因此强度分布的求取效率很高，且不显著增加粒子数目。

在上述过程中，对于已知在射束坐标系中位置和方向的粒子，根据限束部件在射束坐标系中的位置，求得粒子在限束部件坐标系中的位置和方向，然后求得当前粒子与限束部件的交线长度。不管粒子的入射位置与方向如何，总能在交线长度矩阵中插值求得对应的交线长度，也就是说交线长度矩阵是可以被无限次重复利用的。而且坐标转换与插值计算占用的计算机时间与实时计算相比短的多。因此本实施例有助于降低实时计算时间。并且，限束部件的端面及侧面形状不受限制，其对粒子的影响结果反映在交线长度矩阵中，因此可以在不增加存储空间和计算时间的情况下，充分考虑端面效应、漏射以及凹凸效应，使得计算的强度分布更加精确。

图3是本发明另一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。参考图3所示，方法可包括以下步骤：

步骤310，预处理。

这一步骤可由图1所示的预处理模块110完成。在此实施例中，这一步骤可包括：

步骤311，产生不同入射位置的粒子与不同入射方向的粒子与射束调整装置的各限束部件的交线长度矩阵。在此，不限定入射到射束调整装置的粒子的产生方式。此外，形状相同或相似的限束部件可以共用一个交线长度矩阵，而不必重复计算。

步骤312，根据各限束部件配置信息计算出各限束部件坐标系相对于射束坐标系的平移或平移与旋转。

在此实施方式中，需要计算各限束部件坐标系相对于射束坐标系的位置变化。通常情况下，为了形成开口，各限束部件会进行平移移动，因此，需要计算出各限束部件坐标系相对于射束坐标系的平移；有时，由于各限束部件的安装方位不同或其它原因，有的限束部件相对于射束坐标系还会存在一定的旋转，此时需要计算出各限束部件坐标系相对于射束坐标系的平移与旋转。

步骤313，产生各限束部件的散射强度分布。

预处理步骤可以在强度计算之前进行从而不占用强度计算时间，处理产生的数据可先行储存，储存的数据只需占用少量存储空间。在之后的强度计算时，只要射束调整装置没有变化，则预处理的步骤不需要重复进行，直接调用储存的数据即可。

步骤320，计算粒子经过射束调整装置后的强度分布。

在此实施例中，计算原粒子强度分布和散射强度分布，其集合作为强度分布。具体步骤包括：

步骤321，获取粒子穿过射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向。

步骤322，使用各限束部件的交线长度矩阵和粒子在射束坐标系中的位置与方向来计算粒子与各限束部件的交线长度，并求和获得粒子经过射束调整装置的交线长度。

步骤323，使用粒子穿过射束调整装置前的权重和/或能量以及粒子穿过射束调整装置的交线长度来计算粒子经过射束调整装置后的权重和/或能量改变。

步骤324，将粒子的位置由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

步骤325，利用各限束部件在射束坐标系中的坐标，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

上述的步骤313的示例性做法是，将射束调整装置的各限束部件的形状和材料在蒙特卡罗程序(常见的如EGSNRC，MCNP)中进行建模，将经过射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序，模拟得到其中每个限束部件产生的次级粒子的能量、位置、权重与方向。由于有些粒子没有权重或者没有能量，因此在将经过该射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序之前，对于无权重的粒子，定义其权重为任意正的常数且各粒子的权重相等，对于无能量的粒子，定义其能量为任意正的常数且各粒子的能量相等。对于形状及材料相同或相似的限束部件，不需要利用蒙特卡罗程序重复模拟，可以共用一个散射强度分布。

步骤321-324的进一步细节可以参考图2所示实施例，在此不再展开。

在上述过程中，对于已知在射束坐标系中位置和方向的粒子，根据限束部件在射束坐标系中的位置，求得粒子在限束部件坐标系中的位置和方向，然后求得当前粒子与限束部件的交线长度。不管粒子的入射位置与方向如何，总能在交线长度矩阵中插值求得对应的交线长度，也就是说交线长度矩阵是可以被无限次重复利用的。而且坐标转换与插值计算占用的计算机时间与实时计算相比短的多。因此本实施例有助于降低实时计算时间。并且，限束部件的端面及侧面形状不受限制，其对粒子的影响结果反映在交线长度矩阵中，因此可以在不增加存储空间和计算时间的情况下，充分考虑端面效应、漏射以及凹凸效应，使得强度分布更加精确。

图4是本发明又一实施例的确定经过射束调整装置的粒子强度分布的方法流程图。参考图4所示，方法可包括以下步骤：

步骤410，预处理。

这一步骤可由图1所示的预处理模块110完成。在此实施例中，这一步骤可包括：

步骤411，产生不同入射位置的粒子与不同入射方向的粒子与射束调整装置的各限束部件的交线长度矩阵。在此，不限定入射到射束调整装置的粒子的产生方式。此外，形状相同或相似的限束部件可以共用一个交线长度矩阵，而不必重复计算。

步骤412，根据各限束部件配置信息计算出各限束部件坐标系相对于射束坐标系的平移或平移与旋转；

在此实施方式中，需要计算各限束部件坐标系相对于射束坐标系的位置变化。通常情况下，为了形成开口，各限束部件会进行平移移动，因此，需要计算出各限束部件坐标系相对于射束坐标系的平移；有时，由于各限束部件的安装方位不同或其它原因，有的限束部件相对于射束坐标系还会存在一定的旋转，此时需要计算出各限束部件坐标系相对于射束坐标系的平移与旋转。

步骤413，产生各限束部件的散射强度分布；

预处理步骤可以在强度计算之前进行从而不占用强度计算时间，处理产生的数据可先行储存，储存的数据只需占用少量存储空间。在之后的强度计算时，只要射束调整装置没有变化，则预处理的步骤不需要重复进行，直接调用储存的数据即可。

步骤420，计算粒子经过射束调整装置后的强度分布。

在此实施例中，计算原粒子强度分布和散射强度分布，其集合作为强度分布。具体步骤包括：

步骤421，获取粒子穿过射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向；

步骤422，使用各限束部件的交线长度矩阵和粒子在射束坐标系中的位置与方向来计算粒子与各限束部件的交线长度，并求和获得粒子经过射束调整装置的交线长度；

步骤423，使用粒子穿过射束调整装置前的权重和/或能量以及粒子穿过射束调整装置的交线长度来计算粒子经过射束调整装置后的权重和/或能量改变；

步骤424，将粒子的位置由射束调整装置前变换到射束调整装置后，得到原粒子强度分布；

步骤425，利用各限束部件在射束坐标系中的坐标，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

步骤426，根据设置的权重阈值和存活概率筛选穿过射束调整装置的粒子。

步骤426使用称为俄罗斯轮盘赌的技巧。具体地说，可以设置阈值W0和存活概率P0，若粒子权重小于W0，则在0到1之间抽样得到一个随机数，若该随机数小于P0，则判定该粒子存活，将其权重除以P0；若该随机数大于或等于P0，则将该粒子丢弃。

存活的粒子参与接下来的模体内的剂量计算，丢弃的粒子将不参与接下来的模体内的剂量计算。剂量计算是指根据一定的粒子强度分布和模体材料、密度等分布，计算出粒子在模体内造成的单位质量的能量沉积。

上述的步骤421-424的进一步细节可以参考图2所示实施例，在此不再展开。上述的步骤425的进一步细节可以参考图3所示实施例，在此不再展开。

在上述过程中，对于已知在射束坐标系中位置和方向的粒子，根据限束部件在射束坐标系中的位置，求得粒子在限束部件坐标系中的位置和方向，然后求得当前粒子与限束部件的交线长度。不管粒子的入射位置与方向如何，总能在交线长度矩阵中插值求得对应的交线长度，也就是说交线长度矩阵是可以被无限次重复利用的。而且坐标转换与插值计算占用的计算机时间与实时计算相比短的多。因此本实施例有助于降低实时计算时间。并且，限束部件的端面及侧面形状不受限制，其对粒子的影响结果反映在交线长度矩阵中，因此可以在不增加存储空间和计算时间的情况下，充分考虑端面效应、漏射以及凹凸效应，使得强度分布更加精确。再者，使用俄罗斯轮盘赌技巧，可以大大减少穿过射束调整装置后的粒子数目，从而减少模体内的剂量计算时间。

上述的实施例可以在医疗领域中实施，用于确定经过射束调整装置的粒子的强度分布，其中射束调整装置可以是多叶准直器(MLC)、钨门等医疗领域中的准直器。以MLC为例，图5是根据本发明一实施例的MLC与射束示意图。参考图5所示，多叶准直器500具有许多叶片501，这些叶片501的组合决定了多叶准直器的开口形状。粒子射束510从MLC 500的入射面(在图5中为上表面)方向入射，经过MLC 500的各个叶片501后，从MLC500的出射面(在图5中为下表面)射出。确定粒子经过MLC后的强度分布是剂量计算的前提。可以使用图2-图4所示的任一实施例的方法来确定粒子强度分布。

下面以MLC为例描述上述方法的进一步细节。

图6是根据本发明一实施例的叶片坐标系、射束坐标系与粒子示意图。参考图6所示，按照各步骤211，311和411，粒子与叶片的交线长度矩阵的产生过程示例如下：

建立叶片坐标系。如图6所示，L为MLC中一个叶片，OXYZ为叶片坐标系，其选取方式为：选取叶片端面最外侧的点为原点，以距叶片上、下底面距离相等、由叶片内部指向叶片端面的方向为X轴正方向，以垂直于叶片上、下底面、从上指向下为Z轴正方向，以平行于叶片上、下底面的直线且从内指向外为Y轴正方向。

叶片的表示。在数学上，一个叶片由多个曲面(包括平面)围成，而每个曲面都可以用一个方程表示，曲面所隔开的两个空间用该方程对应的不等式表示，因此一个叶片可以用多个不等式表示。同时满足这些不等式的所有点的集合即为该叶片。若存在多种不同形状的叶片，则对于每一种不同形状的叶片，都有一个不等式集合。

区间划分。如图6所示，将粒子Q沿着方向矢量V投影到叶片坐标系OXY平面并交于点P，此时其坐标z值变为0。将此时x坐标，y坐标，方位角t和天顶角f分别取一个范围并分成多个区间，每个区间的中点所对应的x坐标，y坐标，方位角t和天顶角f构成一个粒子的可能状态。所述方位角t为粒子的方向矢量V在X-Y平面上的投影V'相对于X轴正方向旋转的角度，天顶角f为粒子的方向矢量V相对于Z轴正方向旋转的角度。

获得交线长度。将x坐标，y坐标，方位角和天顶角的值分别取某个区间的中点，构成一个粒子的可能状态。求得该粒子在叶片上底面所在平面的交点M与叶片下底面所在平面的交点N，将这两个交点构成的线段MN两端做少量延长构成扩展的线段M'N'，将扩展的线段分成多个区间。对每个区间做下面的循环：判断该区间中点是否满足所有不等式，若是，则将交线长度值加上该区间长度，然后判断下一个区间，若否，则直接判断下一个区间。将扩展的线段分成的所有区间全部操作后，得到的交线长度存在内存的一处。

获得交线长度矩阵。遍历所有的x，y，方位角t和天顶角f的组合，对于每个组合求具有该x，y，方位角t和天顶角f的粒子与叶片的交线长度，并按顺序进行存储。

具体的说，为交线长度矩阵文件设置文件头，文件头用于确定x，y，方位角和天顶角的具体值，交线长度矩阵内的值和文件头内的x，y，方位角和天顶角的值对应。

文件头的格式有多种，文件头内数据的顺序可以自行设置，在本发明中示例如下，但不限于此。

设置文件头(min(x)，resolution(x)，num(x))，(min(y)，resolution(y)，num(y))，(min(t)，resolution(t)，num(t))，(min(f)，resolution(f)，num(f))，其中min(x)表示x坐标的最小值，resolution(x)表示x坐标的分辨率，即相邻区间的中点之间的距离，num(x)表示x坐标的数量；类似的，min(y)、min(t)、min(f)分别表示y坐标、方位角、天顶角的最小值，resolution(y)、resolution(t)、resolution(f)表示y坐标、方位角、天顶角的分辨率，即各自相邻区间的中点之间的距离，num(y)、num(t)、num(f)分别表示y坐标、方位角、天顶角的数量；从x，y，方位角和天顶角的最小值开始，即从min(x)，min(y)，min(t)和min(f)的组合开始计算对应粒子与叶片的交线长度并存储，然后保持min(y)，min(t)和min(f)不变，以x坐标的分辨率为单位，逐渐增大x坐标，求取每个x，y，方位角t和天顶角f的组合对应的粒子与叶片的交线长度并按顺序进行存储；再以min(y)为基础，以y坐标的分辨率为单位，增大一次y坐标，再保持增大之后的y坐标、min(t)和min(f)不变，然后从min(x)开始，以x坐标的分辨率为单位，逐渐增大x坐标，求取每个x，y，方位角t和天顶角f的组合对应的粒子与叶片的交线长度并按顺序进行存储；依次类推，遍历所有的x，y，方位角t和天顶角f的组合，对于每个组合求具有该x，y，方位角t和天顶角f的粒子与叶片的交线长度，并按顺序进行存储，从而得到交线长度矩阵，并且该交线长度矩阵里的数值和文件头内的x，y，方位角t和天顶角f的值一一对应。依照上述方法，对所有形状的叶片都求得交线长度矩阵，并存于文件中。上述方法只是示出其中一种遍历的方式，本领域技术人员可以选用其它的遍历方式，只要得到所有x，y，方位角t和天顶角f组合对应的粒子与所有形状的叶片的交线长度即可。

然而需要指出，交线长度矩阵的产生方式并不限于上述方式，也可以通过求解粒子位置、方向所在直线与叶片曲面的方程获得。

按照各步骤212，312和412计算坐标变换时，由叶片配置信息计算出每个叶片坐标系相对于射束坐标系的平移、旋转。平移的原因有：为了形成开口而进行的叶片相对于X方向的平行移动，叶片安装位置不在射束坐标系原点，不同的叶片y坐标不同。由于各叶片的安装方位不同或其它原因，有的叶片相对于射束坐标系还会存在一定的旋转。

按照上述的各步骤222，322和422计算粒子与叶片交线长度时，对于一个MLC中的每个叶片，利用粒子在射束坐标系中的位置与方向x1,y1,t1,f1，及叶片相对于射束坐标系的平移与旋转，求得粒子相对于该叶片坐标系的位置与方向x2,y2,t2,f2，此处粒子的方向用方位角和天顶角表示。将该叶片对应的形状的交线长度矩阵从文件中读出，根据x2,y2,t2,f2查找到对应的交线长度并读取，如果查找不到对应x2,y2,t2,f2组合的交线长度，则利用相邻的组合进行插值求取x2,y2,t2,f2对应的交线长度，即为该粒子与该叶片的交线长度。插值为本领域常用的技术手段，在此不再赘述。进一步，求粒子与所有叶片的交线长度并求和，即为粒子穿过MLC的交线长度。若存在多个MLC，则分别求得粒子穿过每个MLC的交线长度。

以上得到的原粒子强度分布和散射强度分布，合并之后即为强度分布。此处“合并”是指本领域技术人员根据原粒子强度分布的形式和散射强度分布的形式采用对应的方法得到总的强度分布，总的强度分布的计算方法为本领域的现有技术，在此不再赘述。可以选择性地对粒子的状态进行统计，得到矩阵形式的强度分布。强度分布可以直接用于接下来的模体内的剂量计算。

本发明的实施例在节约计算时间方面的优势是明显的。原因在于，本发明实施例在预处理阶段求取不同入射位置的粒子与不同入射方向的粒子与限束部件的交线长度矩阵，实时计算的任务仅仅主要包括：对于已知在射束坐标系中位置和方向的粒子，根据限束部件在射束坐标系中的位置，求得粒子在限束部件坐标系中的位置和方向，然后求得当前粒子与限束部件的交线长度，进而利用一个指数衰减公式即可求得该粒子穿过射束调整装置后的状态。不管粒子的入射位置与方向如何，总能在交线长度矩阵中求得对应的交线长度，也就是说交线长度矩阵是可以被无限次重复利用的。相比之下，如果按照现有技术第二种和第三种方法实时计算粒子与限束部件交线长度，则粒子的入射位置与方向的组合是无穷多种的。对于每种组合，都要单独地计算交线长度，随着粒子数的增加，占用的计算机时间呈线性增长。而在本发明中，坐标转换与插值计算占用的计算机时间与实时计算相比短的多。

以下三种方法的计算时间在如下计算机配置的情况下适用：CPU速度不小于1.8GHz，内存不小于1GB RAM，硬盘不小于40GB。

以一个完整的大量粒子求强度分布的问题为例，第二种方法由于实时模拟粒子在MLC中的反应，其时间为数小时。第三种方法实时计算每个粒子在每个叶片中的交线长度，耗时相当长，以至于不能实际地用于临床治疗(Yan Chen在2000年提出的“Calculation of x-ray transmission through a multileafcollimator”第1721页D节指出“It may not practical to use full ray tracingcalculation for clinical applications because it is too slow.”)。本发明在预处理阶段计算叶片交线矩阵的方法耗时可以控制在1分钟之内，而实时强度分布计算时间控制在数分钟之内。

综上所述，本发明上述实施例的方法相比已知的方法具有如下的一个或多个优势：

(1)直接获取粒子与限束部件交线长度矩阵及限束部件散射强度分布，只需要进行坐标变换和叠加就可以求得粒子穿过射束调整装置之后的强度分布，相比实时计算交线长度的方法，本发明为强度计算节省大量时间。(2)在本发明中，限束部件的端面及侧面形状不受限制，其对粒子的影响结果反映在交线长度矩阵中，因此可以在不增加存储空间和计算时间的情况下，充分考虑端面效应、叶间漏射以及凹凸效应，使得计算的强度分布更加精确。

(3)使用俄罗斯轮盘赌技巧，大大减少穿过射束调整装置及后续射束调整装置的粒子数目，从而减少总的强度计算时间和模体内的剂量计算时间。

(4)每增加一个开口，只需增加一次坐标变换，因此强度分布的求取效率很高(平均一个开口的计算时间短)；且对于只修改权重的粒子来说不显著增加粒子数目，因此随后的粒子在病人或模体中的剂量计算时间效率也很高。

虽然本发明已参照当前的具体实施例来描述，但是本技术领域中的普通技术人员应当认识到，以上的实施例仅是用来说明本发明，在没有脱离本发明精神的情况下还可作出各种等效的变化或替换，因此，只要在本发明的实质精神范围内对上述实施例的变化、变型都将落在本申请的权利要求书的范围内。

Claims (18)

1.一种确定经过射束调整装置的粒子强度分布的方法，包括以下步骤：计算粒子经过射束调整装置后的强度分布，包括：

获取粒子穿过该射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置、与方向；

使用该射束调整装置的各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度，其中所述交线长度矩阵包括不同入射位置的粒子与不同入射方向的粒子与对应的限束部件的交线长度；

使用粒子穿过该射束调整装置前的权重和/或能量以及粒子穿过该射束调整装置的交线长度来计算粒子经过该射束调整装置后的权重和/或能量改变；以及

将粒子的位置沿着粒子的方向由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

2.根据权利要求1所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，在计算粒子经过射束调整装置后的强度分布之前，还包括预处理，所述预处理包括产生不同入射位置的粒子与不同入射方向的粒子与该射束调整装置的各限束部件的交线长度矩阵。

3.根据权利要求2所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，使用交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度包括：首先计算粒子与射束调整装置的各限束部件的交线长度，再求和获得粒子经过该射束调整装置的交线长度。

4.根据权利要求3所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，使用该射束调整装置的各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度包括：

利用粒子在该射束坐标系中的位置与方向，及各限束部件坐标系相对于该射束坐标系的平移或平移与旋转，求得粒子在各限束部件坐标系中的位置与方向；

利用各限束部件对应的交线长度矩阵，获取与粒子在各限束部件坐标系中的位置与方向所对应的交线长度，分别作为粒子与各限束部件的交线长度。

5.根据权利要求4所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，预处理还包括根据各限束部件配置信息计算出各限束部件坐标系相对于该射束坐标系的平移或平移与旋转。

6.根据权利要求4所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，在使用交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子经过该射束调整装置的交线长度的步骤中，在求得粒子在各限束部件坐标系中的位置与方向之前还包括：

判断粒子是否穿过且只穿过限束部件部分高度，如果是则在使用该交线长度矩阵的条件下求得粒子在各限束部件中的交线长度，否则在不使用该交线长度矩阵的条件下求得粒子在限束部件中的交线长度。

7.根据权利要求2所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，该预处理还包括产生该射束调整装置的各限束部件的散射强度分布。

8.根据权利要求7所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，计算粒子经过射束调整装置后的强度分布还包括利用各限束部件坐标系相对该射束坐标系的平移或平移与旋转，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

9.根据权利要求7所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，产生该射束调整装置的各限束部件的散射强度分布包括：将该射束调整装置的各限束部件的形状和材料在蒙特卡罗程序中进行建模，将经过该射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序，模拟得到在每个限束部件中产生的次级粒子的能量、位置、权重与方向，作为各限束部件的散射强度分布。

10.根据权利要求9所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，在将经过该射束调整装置之前的粒子的能量、位置、权重与方向，输入到蒙特卡罗程序之前，对于无权重的粒子，定义其权重为任意正的常数且各粒子的权重相等，对于无能量的粒子，定义其能量为任意正的常数且各粒子的能量相等。

11.根据权利要求1所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，计算粒子经过射束调整装置后的强度分布还包括根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子。

12.根据权利要求11所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，该根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子包括：设置阈值W0和存活概率P0，若粒子权重小于W0，则在0到1之间抽样得到一个随机数，若该随机数小于P0，则判定该粒子存活，将其权重除以P0；若该随机数大于或等于P0，则将该粒子丢弃。

13.根据权利要求1所述的确定经过射束调整装置的粒子强度分布的方法，其特征在于，该射束调整装置为多叶准直器，该限束部件为多叶准直器的叶片。

14.一种确定经过射束调整装置的粒子强度分布的装置，包括：

预处理模块，配置为产生不同入射位置与方向的粒子与该射束调整装置的各限束部件的交线长度矩阵，其中各限束部件确定限定该射束调整装置的形状；

强度计算模块，耦接该预处理模块，配置为按照下述步骤计算粒子经过该射束调整装置后的强度分布：

获取粒子穿过该射束调整装置前的状态，该状态包括粒子在射束坐标系中的权重和/或能量、位置与方向；

使用各限束部件的交线长度矩阵和该粒子在射束坐标系中的位置与方向来计算粒子与各限束部件的交线长度，并求和获得粒子经过该射束调整装置的交线长度，其中所述交线长度矩阵包括不同入射位置的粒子与不同入射方向的粒子与对应的限束部件的交线长度；

使用粒子穿过该射束调整装置前的权重和/或能量以及粒子穿过该射束调整装置的交线长度来计算粒子经过该射束调整装置后的权重和/或能量改变；以及

将粒子的位置由射束调整装置前变换到射束调整装置后，得到原粒子强度分布。

15.根据权利要求14所述的确定经过射束调整装置的粒子强度分布的装置，其特征在于，该预处理模块还根据各限束部件配置信息计算出各限束部件坐标系相对于该射束坐标系的平移或平移与旋转。

16.根据权利要求14所述的确定经过射束调整装置的粒子强度分布的装置，其特征在于，该预处理模块还产生各限束部件的散射强度分布。

17.根据权利要求15所述的确定经过射束调整装置的粒子强度分布的装置，其特征在于，该强度计算模块还利用各限束部件坐标系相对该射束坐标系的平移或平移与旋转，将各限束部件的散射强度分布中的位置进行变换，计算新的散射强度分布。

18.根据权利要求14所述的确定经过射束调整装置的粒子强度分布的装置，其特征在于，该强度计算模块还根据设置的权重阈值和存活概率筛选穿过该射束调整装置的粒子。