CN103949015A - 一种伽玛刀放射治疗逆向计划设计方法 - Google Patents

Abstract

本发明公开一种伽玛刀放射治疗逆向计划设计方法，首先，进行重建靶区和危及器官感兴趣区域三维模型；建立优化模型；确定优化参数与遗传算法参数；整个过程应用图形处理器(GPU)并行计算进行加速，随后根据改进的遗传算法进行求解，包括：基因编码初始化族群；计算评估个体适应度选择繁殖个体；个体繁殖变异更新族群；调整权重，最终得到可行的治疗计划。本发明的优点为：优化模型具有较好的适应性，可以通过调整目标函数中的权重适应不同的病症；针对伽玛刀治疗计划设计的特点，改进了遗传算法的基因编码，种群初始化，繁殖和变异规则，并根据优化模型提出了一种灵活有效的个体适应度评价函数。

Description

一种伽玛刀放射治疗逆向计划设计方法

技术领域

本发明涉及计算机技术领域，具体来说，是一种伽玛刀放射治疗逆向计划设计方法。 

背景技术

放射治疗是治疗癌症最主要的手段之一。伽玛刀放射治疗使用伽玛射线照射癌变的肿瘤组织，在尽可能减少损害正常组织的前提下，抑制癌变细胞的生长、繁殖和扩散。为了达到治疗的目的，一般在对病人进行治疗之前，通过治疗计划系统(TPS)制定治疗计划，使靶区内的剂量具有较好的覆盖度(Coverage)和适形度(Conformity)；危及器官(OAR)承受尽可能低的剂量(Avoidance)；同时兼顾治疗计划的可行性。传统的治疗计划正向设计是一个不断试错修改(trial and error)的过程：物理师根据临床经验设计出初步的治疗计划，然后基于靶区和危及器官内的剂量分布和治疗计划的可行度进行调整，重复这个步骤直到得出一个可行的治疗计划。然而，计划设计通常需要综合权衡考虑多个目标，设计的过程繁重且耗时较长，而且在很大程度上依赖于物理师的临床经验。不同于正向计划设计，逆向计划设计根据预期的剂量分布定义优化目标和约束条件，建立优化模型，从而将计划设计转化为优化问题。借助计算机强大的运算能力，使用高效的优化算法，逆向计划设计可以快速的制定出治疗计划，还可以提高治疗计划的可行度，同时降低了对临床经验的依赖程度。 

现有的优化方法主要针对头部伽玛刀提出，主要包括非线性优化方法和基于几何特征的优化方法。其中，针对头部伽玛刀的非线性优化方法，通过查表的方法计算剂量，进而建立优化模型并求解。对于体部伽玛刀，由于体部由于肋骨等骨骼的遮挡，不同位置靶点剂量分布相差较大，不能通过查表法计算剂量分布。此外，非线性优化方法求解时间较长，不适用于靶区较大的病例。由于头部伽玛刀使用直径较小的准直器，作用范围较小，不同靶点之间剂量分布相互影响较小。基于几何特征的优化方法近似的认为单个靶点的剂量分布为球形，不同直径的准直器对应大小不同的球体，在各个球体不相互交叠的前提下，利用尽可能少的球 体填充形状体积已知的靶区，使靶区内的覆盖度尽可能的高。但是对于体部伽玛刀伽玛刀，准直器直径较大，距离较近的靶点之间剂量分布相互影响较大，因此，基于几何特征的优化方法也不适用于体部伽玛刀治疗计划设计。 

发明内容

针对上述问题，本发明提出了一种适用于伽玛刀放射治疗逆向计划设计方法，建立了一种基于剂量分布的优化模型，并应用改进的遗传算法(Genetic Algorithm)进行求解。 

一种伽玛刀放射治疗逆向计划设计方法，具体通过下述步骤完成： 

步骤1：重建感兴趣区域的三维模型，并对靶区部分进行膨胀处理生成靶区边缘轮廓。 

步骤2：建立优化模型； 

优化目标为：1、保证靶区剂量的覆盖度；2、提高靶区剂量的适形度；3、降低危及器官承受的剂量。 

步骤3：确定步骤2建立的优化模型参数。 

步骤4：确定遗传算法参数。 

步骤5：基因编码及初始化族群。 

步骤6：对每个靶点序列进行剂量计算。 

步骤7：评价族群中各个个体，并通过计算适应度函数得到一个适应度数值。 

步骤8：选择繁殖个体； 

将族群中的个体按适应度数值由高到底排序，选取前m个个体两两任意配对，作为父代和母代个体繁殖下一代。m为族群中个体数与繁殖概率的乘积。 

步骤9：父代和母代个体繁殖； 

A、确定父代和母代个体的重心，得到重心位移方向向量d； 

父代和母代个体的重心位置H为： 

$H=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\φ}}_i\×S_i\×P_i$

其中，φ_i为靶点序列中第i个靶点对应的靶点准直器的半径；S_i为靶点序列中第i个靶点对应的靶点的使能状态；P_i为靶点序列中第i个靶点的位置坐标；n为靶点序列中靶点数量； 

B、以d为方向向量，随机的取正向或者反向，整体移动父代和母代中所有靶点，得到子代个体。 

步骤10：子代个体变异； 

令子代个体中的靶点i变异的平移方向向量Im_i，则： 

${\mathrm{Im}}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{(P_i-P_j)}\×S_j\×{\mathrm{\φ}}_j\×R_j,j\≠i$

其中，R随机取±1，j为子代个体中靶点i外的其余靶点。 

步骤11：将族群中进行繁殖的父代母代个体由其繁殖的子代个体代替，而未进行繁殖的个体保持不变。 

步骤12：重复步骤5～11进行迭代运算，直至满足终止条件后，进行步骤9。 

步骤13：调整权重 

对更新后的族群中最优个体的每个靶点，增加权重，随后计算个体的适应度；若适应度减小，继续增加该靶点权重直至适应度值不再减小；若适应度增大，则需减小该靶点的权重以提高个体适应度；因此将靶点权重重置，并减小权重，直至适应度不再减小。 

本发明的优点为： 

1、本发明伽玛刀放射治疗逆向计划设计方法，针对伽玛刀治疗计划设计的特点，提出了一种基于剂量分布的优化模型；首先，优化模型将优化目标转化为能够使用优化方法进行求解的数学描述；同时，优化模型具有较好的适应性，可以通过调整目标函数中的权重适应不同的病症； 

2、本发明伽玛刀放射治疗逆向计划设计方法，针对伽玛刀治疗计划设计的特点，改进了遗传算法的基因编码，种群初始化，繁殖和变异规则，并根据优化模型提出了一种灵活有效的个体适应度评价函数； 

3、本发明伽玛刀放射治疗逆向计划设计方法，优化过程应用图形处理器(GPU)并行计算进行加速，能够快速高效的得出可行的治疗计划。 

附图说明

图1为本发明伽玛刀放射治疗逆向计划设计方法流程图； 

图2为本发明伽玛刀放射治疗逆向计划设计方法中权重调节流程图。 

具体实施方式

下面结合附图对本发明做进一步说明。 

本发明为伽玛刀放射治疗逆向计划设计方法，如图1所示，具体通过下述步骤完成： 

步骤1：重建感兴趣区域的三维模型； 

重建包含靶区和危及器官区域部分的感兴趣区域(ROI，Region of  Interesting)的三维模型；并对ROI模型中的靶区部分进行膨胀处理生成靶区边缘轮廓(Margin)。为了提高计算速度，本发明中还对ROI三维模型进行采样处理。 

步骤2：建立优化模型； 

优化目标为： 

1、保证靶区剂量的覆盖度； 

2、提高靶区剂量的适形度； 

3、降低危及器官承受的剂量； 

其中，对于优化目标1来说，通过对靶区剂量的覆盖度进行约束实现，靶区剂量的覆盖度的约束条件为： 

$\frac{\underset{(x,y,z)}{\mathrm{\Σ}}C(x,y,z)}{T_N}\≥C_T---\left(1\right)$

式(1)中，C_T为预设的靶区覆盖度；D(x,y,z)为靶区中对应坐标为(x,y,z)处的剂量；T_N为靶点体素个数；C(x,y,z)为为靶区剂量是否达到当前50％最大剂量的计数点： 

$C(x,y,z)=\left\{\begin{array}{cc}1& \frac{D(x,y,z)}{D_{\max }}\\ 0& \mathrm{else}\end{array}\right.(x,y,z)\∈T---\left(2\right)$

式(2)中，D_max为当前靶区最大剂量值；T为Target，(x，y，z)为靶区的点。 

对于优化目标2和3来说，通过在靶区剂量覆盖度达到C_T的前提下，最小化靶区边缘轮廓和危及器官区域的剂量来实现，则有： 

$\mathrm{Obj}=\underset{(x,y,z)\∈M}{\mathrm{\Σ}}{w}_1\×D(x,y,z)+\underset{(x,y,z)\∈\mathrm{OAR}}{\mathrm{\Σ}}{w}_2\×D(x,y,z)---\left(3\right)$

其中，Obj为目标函数；w₁，w₂分别为靶区边缘轮廓和危及器官区域剂量值的权重，通过调节权重，以适用不同的病例。 

步骤3：确定步骤2建立的优化模型参数； 

优化模型参数包括：靶区最小覆盖度、靶区边缘轮廓和危及器官区域剂量权重、靶点序列个数； 

步骤4：确定遗传算法参数，包括：族群个体数、不同类型准直器变异出现的概率(单个靶点使用某种直径准直器的概率)、遗传算法终止条件(即：最小迭代次数)、以及繁殖概率。 

步骤5：基因编码及初始化族群； 

将靶点序列编码为个体，靶点编码为基因，基因参数包括靶点位置、使能状态和准直器类型(由准直器的半径进行区分)。通过启发性方法确定靶点序列中靶点的初始位置完成族群初始化，具体为： 

a、在靶区中随机确定第一个靶点的位置坐标； 

b、根据第一靶点位置坐标，计算得到剂量场； 

c、搜索当前剂量场中最大剂量值D_max； 

d、对剂量场进行归一化处理，得到第一个靶点的剂量百分比，为： 

$\frac{D_1(i_1,j_1,k_1)}{D_{\max }},(i_1,j_1,k_1)\∈T---\left(4\right)$

e、令第二个靶点在靶区中坐标为(i₂,j₂,k₂)，确定第二个靶点在靶区中出现的概率； 

出现的概率P₂(i₂,j₂,k₂)为当前第二个靶点的剂量百分比与100％的差值： 

$P_2(i_2,j_2,k_2)=1-\frac{D_2(i_2,j_2,k_2)}{D_{\max }},(i_2,j_2,k_2)\∈T---\left(5\right)$

f、根据第二个靶点的在靶区出现的概率，确定第二个靶点在靶区中的具体位置。 

g、根据第二靶点位置，更新剂量场； 

h、令第三个靶点的在靶区中坐标为(i₃,j₃,k₃)的一点，则根据更新后的剂量场，返回执行步骤e～g，得到第三个靶点的位置，并再次更新剂量场；由此，依次遍历全部靶点，即可得到靶点序列中各靶点的位置。 

步骤6：对每个靶点序列进行剂量计算； 

搜索当前场量场得到中最大剂量值D_max，并对剂量场进行归一化处理，统计靶区、靶区边缘轮廓与危及器官区域的剂量。 

步骤7：评价族群中各个个体，并通过计算适应度函数得到一个适应度数值； 

上述适应度函数为： 

$F=\left\{\begin{array}{cc}{w}_1+w_2& \mathrm{if}(\mathrm{Cov}\≤C_T)\\ w_1\×\stackrel{\‾}{D_M}+w_2\×\stackrel{\‾}{D_{\mathrm{OAR}}}& \mathrm{if}(\mathrm{Cov}\≥C_T)\end{array}\right.---\left(7\right)$

式(7)中，和为以当前剂量场中最大剂量归一化后靶区边缘轮廓和危及器官区域的平均剂量值，Cov为靶区覆盖度。和为个体中个个基因对应的剂量场中最大剂量归一化后靶区边缘轮廓和危及器官区域的平均剂量值。 

步骤8：选择繁殖个体； 

将族群中的个体按适应度数值由高到底排序，选取前m个个体两两任意配对，作为父代和母代个体繁殖下一代。m为族群中个体数与繁殖概率的乘积。 

步骤9：父代和母代个体繁殖； 

现有的繁殖规则(One-Point，two-point,cycle和order)的基本原则是按照一定的规律交换个体的部分基因。然而对于伽玛刀计划设计，个体的适应度由个体中所有基因共同决定，交换两个计划的部分基金的繁殖方式显然不再适用。 

因此，本发明提出了一种基于父代个体剂量分布，整体移动个体中全部基因的繁殖方式，具体为： 

A、确定父代和母代个体的重心，得到重心位移方向向量d； 

根据父代和母代个体中靶点对当前剂量分布的贡献值，得到父代和母代个体的重心位置； 

父代和母代个体的重心位置H为： 

$H=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\φ}}_i\×S_i\×P_i---\left(6\right)$

其中，φ_i为靶点序列中第i个靶点对应的靶点准直器的半径；S_i为靶点序列中第i个靶点对应的靶点的使能状态，使能为1，否则为0；P_i为靶点序列中第i个靶点的位置坐标；n为靶点序列中靶点数量。 

B、以d为方向向量，随机的取正向或者反向，整体移动父代和母代中所有靶点，得到子代个体。 

步骤10：子代个体变异； 

本发明提出了一种基于靶点之间相互作用，变异靶点位置的方法。靶点之间的相互影响由靶点的使能状态、相对位置和准直器类型决定；例：靶区中有四个靶点，靶点1、2对应的准直器半径相同；3、4对应的准直器半径相同；1、3、 

4均为使能状态；2为非使能状态。考虑靶点1、2、3对靶点4的作用。使能状态决定是否对其他靶点产生影响，靶点2为非使能状态，所以对靶点4没有影响。准直器的半径越大，靶点对当前剂量分布贡献值越大，同样对其他靶点的作用也越明显；因此，1对4的影响要大于3对4的影响。靶点之间的相互作用与相对距离成反比，距离越近影响越大，反之影响越小。 

令子代个体中的靶点i变异的平移方向向量Im_i，则： 

${\mathrm{Im}}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{(P_i-P_j)}\×S_j\×{\mathrm{\φ}}_j\×R_j,j\≠i---\left(8\right)$

式(8)中，R随机取±1，j为子代个体中靶点i外的其余靶点。 

按照这种基于靶点间相互影响的变异方法，子代个体中靶点的移动方向有2^n-1种可能。对于靶点个数较少的情况，变异平移方向较少，可能导致算法陷入局部最优的死循环。因此，本发明按照的变异概率随机地由子代个体中选取变异个体，保证了基因的多样性。 

步骤7：更新族群； 

将族群中进行繁殖的父代母代个体由其繁殖的子代个体代替，而未进行繁殖的个体则保持不变。 

步骤8：重复步骤5～7进行迭代运算，直至满足终止条件后，进行步骤9； 

步骤9：调整权重 

由于遗传算法没有优化靶点的权重，因此遗传算法得到的最优个体中的每个靶点的权重都为1。对于更新后的族群中最优个体(族群中适应度最高的个体)的每个靶点，本发明均采用以0.1为步长增加权重，随后计算个体的适应度F；若适应度减小，继续增加该靶点权重直至适应度值不再减小。若适应度增大，则需减小该靶点的权重以提高个体适应度；因此将靶点权重重置为1，并以0.1为步长减小权重，直至适应度不再减小。通过上述方法，遍历最有个体中所有靶点后得到的靶点序列，即为治疗计划。 

本发明伽玛刀放射治疗逆向计划设计方法中在进行计算个体中个个靶点剂量场的剂量分布，搜索当前剂量场中最大剂量值，进行归一化处理，以及统计靶区、靶区边缘轮廓和危及器官区域的剂量值评估个体的适应度时，这个过程计算量非常大，采用传统的CPU单线程计算方式耗时严重。因此本发明采用图形处理单元(GPU)对上述过程进行处理，提高了优化速度。 

Claims (8)

1.一种伽玛刀放射治疗逆向计划设计方法，其特征在于：具体通过下述步骤完成：

步骤1：重建感兴趣区域的三维模型，并对靶区部分进行膨胀处理生成靶区边缘轮廓；

步骤2：建立优化模型；

优化目标为：1、保证靶区剂量的覆盖度；2、提高靶区剂量的适形度；3、降低危及器官承受的剂量；

步骤3：确定步骤2建立的优化模型参数；

步骤4：确定遗传算法参数；

步骤5：基因编码及初始化族群；

步骤6：对每个靶点序列进行剂量计算；

步骤7：评价族群中各个个体，并通过计算适应度函数得到一个适应度数值；

步骤8：选择繁殖个体；

将族群中的个体按适应度数值由高到底排序，选取前m个个体两两任意配对，作为父代和母代个体繁殖下一代；m为族群中个体数与繁殖概率的乘积；

步骤9：父代和母代个体繁殖；

A、确定父代和母代个体的重心，得到重心位移方向向量d；

父代和母代个体的重心位置H为：

$H=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\φ}}_i\×S_i\×P_i---\left(1\right)$

其中，φ_i为靶点序列中第i个靶点对应的靶点准直器的半径；S_i为靶点序列中第i个靶点对应的靶点的使能状态；P_i为靶点序列中第i个靶点的位置坐标；n为靶点序列中靶点数量；

B、以d为方向向量，随机的取正向或者反向，整体移动父代和母代中所有靶点，得到子代个体。

步骤10：子代个体变异；

令子代个体中的靶点i变异的平移方向向量Im_i，则：

${\mathrm{Im}}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{1}{(P_i-P_j)}\×S_j\×{\mathrm{\φ}}_j\×R_j,j\≠i---\left(2\right)$

其中，R随机取±1，j为子代个体中靶点i外的其余靶点；

步骤11：将族群中进行繁殖的父代母代个体由其繁殖的子代个体代替，而未进行繁殖的个体保持不变；

步骤12：重复步骤5～11进行迭代运算，直至满足终止条件后，进行步骤9；

步骤13：调整权重

对更新后的族群中最优个体的每个靶点，增加权重，随后计算个体的适应度；若适应度减小，继续增加该靶点权重直至适应度值不再减小；若适应度增大，则需减小该靶点的权重以提高个体适应度；因此将靶点权重重置，并减小权重，直至适应度不再减小。

2.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：对感兴趣区域的三维模型进行采样处理。

3.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：所述优化模型为：

$\frac{\underset{(x,y,z)}{\mathrm{\Σ}}C(x,y,z)}{T_N}\≥C_T---\left(3\right)$

$\mathrm{Obj}=\underset{(x,y,z)\∈M}{\mathrm{\Σ}}{w}_1\×D(x,y,z)+\underset{(x,y,z)\∈\mathrm{OAR}}{\mathrm{\Σ}}{w}_2\×D(x,y,z)---\left(4\right)$

其中，C_T为预设的靶区覆盖度；D(x,y,z)为靶区中对应坐标为(x,y,z)处的剂量；T_N为靶点体素个数；C(x,y,z)为靶区剂量是否达到当前50％最大剂量的计数点：

$C(x,y,z)=\left\{\begin{array}{cc}1& \frac{D(x,y,z)}{D_{\max }}\≥50\%\\ 0& \mathrm{else}\end{array}\right.(x.y.z)\∈T---\left(5\right)$

D_max为当前靶区最大剂量值；T为Target，(x，y，z)为靶区的点；Obj为目标函数；w₁，w₂分别为靶区边缘轮廓和危及器官区域剂量值的权重，通过调节权重，以适用不同的病例。

4.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：所述步骤5中，通过启发性方法确定靶点序列中靶点的初始位置完成族群初始化，具体为：

a、在靶区中随机确定第一个靶点的位置坐标；

b、根据第一靶点位置坐标，计算得到剂量场；

c、搜索当前剂量场中最大剂量值D_max；

d、对剂量场进行归一化处理，得到第一个靶点的剂量百分比，为：

$\frac{D_1(i_1,j_1,k_1)}{D_{\max }},(i_1,j_1,k_1)\∈T---\left(6\right)$

e、令第二个靶点在靶区中坐标为(i₂,j₂,k₂)，确定第二个靶点在靶区中出现的概率；

出现的概率P₂(i₂,j₂,k₂)为当前第二个靶点的剂量百分比与100％的差值：

$P_2(i_2,j_2,k_2)=1-\frac{D_2(i_2,j_2,k_2)}{D_{\max }},(i_2,j_2,k_2)\∈T---\left(7\right)$

f、根据第二个靶点的在靶区出现的概率，确定第二个靶点在靶区中的具体位置。

g、根据第二靶点位置，更新剂量场；

h、令第三个靶点的在靶区中坐标为(i₃,j₃,k₃)的一点，则根据更新后的剂量场，返回执行步骤e～g，得到第三个靶点的位置，并再次更新剂量场；由此，依次遍历全部靶点，得到靶点序列中各靶点的位置。

5.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：步骤7中，适应度函数为：

$F=\left\{\begin{array}{cc}{w}_1+w_2& \mathrm{if}(\mathrm{Cov}\≤C_T)\\ w_1\×\stackrel{\‾}{D_M}+w_2\stackrel{\‾}{D_{\mathrm{OAR}}}& \mathrm{if}(\mathrm{Cov}\≥C_T)\end{array}\right.---\left(8\right)$

其中，和为个体中个个基因对应的剂量场中最大剂量归一化后靶区边缘轮廓和危及器官区域的平均剂量值，Cov为靶区覆盖度。

6.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：步骤10中，变异概率取

7.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：步骤13中，靶点的权重按0.1步长增加或递减。

8.如权利要求1所述一种伽玛刀放射治疗逆向计划设计方法，其特征在于：所述个体中个个靶点剂量场的剂量分布，搜索当前剂量场中最大剂量值，进行归一化处理，以及统计靶区、靶区边缘轮廓和危及器官区域的剂量值评估个体的适应度，均通过图形处理单元计算完成。