针对射程不确定性的质子重离子调强放疗鲁棒优化方法

摘要

本发明公开了一种针对射程不确定性的质子重离子调强放疗鲁棒优化方法，首先导入质子重离子调强放疗计划的基本信息和准备计算的数据，然后根据均值和方差的鲁棒性度量方法构建基于最优性和鲁棒性相结合的射程不确定性数学模型，最后采用基于滤波算法的自适应粒子群鲁棒优化得到鲁棒最优解。本发明不仅可以准确地建立鲁棒优化问题的射程不确定性模型，而且具有求解速度快、算法稳定、鲁棒性强等特点，使用户可以通过鲁棒最优解来减小射程误差对放射治疗计划带来的影响。

说明书

技术领域

本发明属于放射治疗的优化方法领域，涉及一种针对射程不确定性的质子重离子调强放疗鲁棒优化方法。

背景技术

质子重离子放疗技术作为最先进的放疗手段之一，在精准杀灭肿瘤细胞的同时，能够有效保护正常组织(较高的治疗增益比)，是肿瘤治疗的“利器”。质子重离子放疗技术到达肿瘤病灶时瞬间释放大量能量杀灭肿瘤细胞，随后能量迅速衰减，在剂量学上叫做“布拉格峰”，峰前剂量约为峰值的20％。质子重离子调强放射治疗通过优化质子重离子扫描束的能量和强度得到靶区和危及器官的预期剂量目标，是质子重离子放疗主要的模式，其精度受射程、呼吸运动等影响较大。质子重离子放疗犹如“立体定向爆破”，具有有效杀灭肿瘤而不产生明显毒副作用的特点。因此，与光子放疗相比，质子重离子放疗在治疗儿童肿瘤和临近重要器官的肿瘤(脊索瘤、软骨肉瘤等)时更具优势，能明显减少二次原发肿瘤发生机率。

质子重离子调强放疗的治疗精度高，因此误差管理非常重要。质子重离子调强放疗中多个环节都可能存在误差(即不确定性)，并且误差无法消除只能减小。如果质子重离子束布拉格峰深度低于或高于预期，那么肿瘤附近的危及器官可能会受到高剂量照射，引起放疗并发症，甚至危及生命。不同于光子放疗，质子重离子调强放疗难以通过肿瘤靶区(GTV)外扩形成计划靶区(PTV)的方式降低不确定性的影响。虽然质子重离子调强放疗可以产生比光子放疗更优的剂量分布，但如果不能有效控制不确定性，则很难进一步提高疗效。

研究普遍认为影响质子调强放疗精度的主要因素之一是射程不确定性。专家学者提出了质子调强放疗不确定性两种主要来源：一是由CT伪影、菲尔德单位转换、解剖结构变化等引起的射程不确定性；二是由于患者位置偏差或组织密度不均匀导致的患者建模不确定性。有专家认为质子射程不确定性有多种来源，主要源于患者几何解剖结构变化，另外治疗物理量和计算也会导致射程不确定性。

虽然质子重离子调强放疗射程不确定性无法避免，但仍可通过鲁棒优化射程不确定性模型降低其负面影响。目前质子重离子调强放疗射程不确定性模型主要包括两类：极大极小模型和最差剂量分布模型。极大极小模型公式如下：

min_x{max_A'∈R{f(g(x,A'))},x≥0}>

其中x∈Rⁿ表示扫描点的强度向量，A'∈R表示不确定性解集，R＝{A₁,A₂,…,A_n}是剂量贡献矩阵不确定性解集。函数g表示器官上所有体素点的剂量，g(x,A')＝A'x。f是一个罚函数，表示优化剂量和期望剂量之间的偏差。质子重离子束的射程被认为是随机变量。因此，照射剂量和目标函数值也是随机变量。通过引入不确定性解集概率p(A')来建立一般随机性模型，如：

min_x{max{∑_A'∈Rp(A')·f(g(x,A'))},x≥0}>

其中p(A')∈R^|R|，∑_A'∈Rp(A')＝1，根据经验来定义p(A')。从公式(6)中可以看出，极大极小模型是在射程不确定性解集中找出“最差”解，再通过优化使其最小化。

最差剂量分布模型中需要统计射程不确定区间中肿瘤靶区剂量的最小值和危及器官剂量的最大值。虽然最差剂量分布模型在物理学上不成立，但提供了一个保守的下限。模型公式如下：

其中A₀表示期望值的解集，和λ_A'分别表示期望函数和最差剂量罚函数的权重因子，g_normal表示正常情况下体素点的剂量，g_worst-case表示最差情况下体素点的剂量，f是一个罚函数。

上述的两类模型主要问题在于仅根据当前解在扰动区间内的最“差”情况得到鲁棒最优解是不合理的。

粒子群优化算法(PSO)是1995年提出的一种启发式算法，通过鸟群的集体协作和竞争使群体达到觅食目的。定义f(x)为在D维欧氏空间E^D的实函数。x_i＝(x_i1,x_i2,…,x_iD)^T表示第i个粒子的当前位置。p_i＝(p_i1,p_i2,…,p_iD)^T表示第i个粒子曾经历的最佳位置，p_g＝(p_g1,p_g2,…,p_gD)^T表示种群的最佳位置。v_i＝(v_i1,v_i2,…,v_iD)^T表示第i个粒子的速度，fitness_i＝f(x_i)表示第i个粒子的适应度值。粒子通过速度和位置更新的方式产生新一代个体，见式(8)和(9)，整个种群不停的迭代从而达到寻优目的。通常，粒子在第d(1≤d≤D)维的位置变化限制在区域内，速度变化限制在区域内。

v_id(t+1)＝ω·v_id(t)+c₁·r₁·(p_id(t)-x_id(t))+c₂·r₂·(p_gd(t)-x_id(t))>

x_id(t+1)＝x_id(t)+v_id(t+1)(9)

其中，i＝1,2,…,M，d＝1,2,…,D，M为种群数量，ω为惯性权重，c₁为认知参数，c₂为社会参数，r₁和r₂分别为两个相互独立的在[0,1]之间的随机数。

粒子群算法在求解鲁棒优化问题时有不错的表现，但是要解决其易陷入局部极值的问题。

发明内容

为了解决上述现有技术存在的不足之处，本发明提出一种针对射程不确定性的质子重离子调强放疗鲁棒优化方法，可以准确地建立鲁棒优化问题的不确定性模型，并采用基于滤波算法的自适应粒子群鲁棒优化方法求解该模型，求解速度快，算法稳定，鲁棒性强，用户可以得到鲁棒最优解。

本发明为解决技术问题采用如下技术方案：

针对射程不确定性的质子重离子调强放疗鲁棒优化方法，包括如下具体步骤：

步骤1：构建基于最优性和鲁棒性相结合的射程不确定性数学模型；模型中体现了射程不确定性引起的放疗计划最优性和鲁棒性，具体如下：

上述公式中，表示不同射程下的目标函数值，体现了最优性。体现了鲁棒性，即不同射程下目标函数值的偏差。w表示质子重离子束的笔形束的强度，z表示器官采样点在射程为下受到的实际剂量。表示质子重离子束集合，b(i)表示射束i的笔形束集合(i∈B)。第i个射束的第j个笔形束的射程用表示，表示正常射程，表示高出/低于正常射程的最大值，第i个射束的第j个笔形束的射程表示成对于任意的射束i∈B，采用表示笔形束j∈b(i)的不确定性离散后的可能射程集，R表示射程集中的元素个数。笔形束j和里的每个元素是一一对应的，所以任意射程由索引唯一确认。使用表示射程索引元素的向量。为了方便定义是正常射程的索引。采用分别表示靶区(Target)和危及器官(OAR)，v表示第v个体素。和分别表示靶区和危及器官的权重。采用w_j表示第j个笔形束的强度，θ_v表示第ν个体素的期望剂量值。

步骤2：采用基于滤波算法的自适应粒子群鲁棒优化方法优化得到鲁棒最优解。

步骤2.1：首先采用滤波算法对违反约束的解进行“过滤”处理。公式(1)中的解表示笔形束的权重，权重是非负数，如果在优化过程中得到的权重出现负数，则采用置零的滤波方式进行处理。

步骤2.2：为了避免粒子群算法陷入局部最优，采用一种有效的动态变异策略。引入多样性函数D和适应度方差函数σ²。

f＝max{1,max{|f(i)-f_avg|}}>

其中，M为种群数量，L为搜索半径。f(i)为第i个粒子的适应度函数值，f_avg为当前粒子群适应度函数值的平均值。

粒子速度的变异函数如下所示：

其中，v_i为第i个粒子的速度，t为迭代次数，G^max为迭代的总次数，θ为服从标准正态分布N(0,1)的随机数。

步骤2.3：采用粒子群算法优化公式(1)，优化过程中采用步骤2.1的方法对违反约束的解进行滤波处理，同时粒子采用公式(2)和公式(3)计算得到D和σ²，如果D和σ²超过预先设定的阈值，则采用公式(4)进行变异操作。满足停止准则后，最终得到针对射程不确定性的质子重离子调强放疗鲁棒最优计划。

该方法还包括在步骤1之前进行如下顺序操作：

1)导入放疗计划信息

放疗计划信息包括患者CT影像、射束参数、加速器硬件参数、器官约束参数、射程不确定性参数；采用标准DICOM RT格式作为输入输出接口，导入商用计划系统的放疗计划文件，完成放疗计划后导出质子重离子扫描点权重和剂量信息到DICOM文件，支持的标准文件主要有：CT、DICOM RT-Structure、DICOM RT-Plan、DICOM RT-Dose；

2)计算剂量贡献矩阵

计算不同能量不同位置质子重离子扫描点对所有器官采样点的剂量贡献，采用快速解析方法-笔形束剂量计算算法，将待计算位置处质子重离子的通量权重与对应水等效深度的IDD数据作卷积叠加，用公式可以表示为：

Dose(r，z)＝Σ_kΦ_k(r，z)·IDD(z)>

其中z代表沿束流的深度位置，r代表横向位置，k代表质子重离子笔形束斑的索引，Dose代表某一个位置对应的体素内的总剂量，Φ是质子重离子的通量权重，IDD是预先通过实际测量或者Monte Carlo方法计算好的质子重离子笔形束在水中沿深度方向的积分深度剂量数据。

顺序步骤2)具体步骤为：

首先，在射束方向观方向上根据质子重离子束斑半高宽进行布点；再将每个束斑划分成19个子射束；然后采用改进的Siddon算法追踪每一子射束，计算射线路径上每个体素对应的水等效深度，查表得到对应的IDD值；根据Fermi-Eyges理论，采用二维校准方式，计算射线穿过的体素所对应多库伦散射的高斯分布标准差σ_MCS；然后根据Soukup等提出的理论公式计算核散射对束流通量横向散射的影响，得到核散射的权重ω_NS以及其高斯分布的标准差σ_NS；考虑到高斯分布的截断误差δ，得到总的通量权重公式为：

Φ(x,y,z)＝(1-ω_NS)·[G(x,y,σ_MCS)+δ_MCS]+ω_NS·[G(x,y,σ_NS)+δ_NS](11)

根据束流监测系统敏感因子S_m，将所有结果按权重叠加，得到最终的剂量dose(x,y,z)：

dose(x,y,z)＝∑_e∑_jS_m·ω_e,j·D_e,j(x,y,z)>

其中x,y,z代表位置坐标，e代表能量层索引，j代表笔形束斑索引，S_m是与束流能量有关的监测系统敏感因子，ω_e,j是能量层e第j个笔束的权重，D_e,j(x,y,z)是能量层e第j个笔束对位置(x,y,z)的剂量。

与现有技术相比，本发明的有益效果在于：

1、本发明构建了一种新型的质子重离子调强放疗射程不确定性模型，克服了之前模型根据当前解在扰动区间内的最“差”情况进行优化的缺点；另外，根据质子重离子调强放疗射程不确定性模型的特点，提出了一种基于滤波算法的自适应粒子群鲁棒优化方法，提高了鲁棒优化的效果，改善了质子重离子放疗计划的治疗效果。

2、本发明构建的新型射程不确定性模型更加合理。前人研究的极大极小模型和最差剂量分布模型可在一定程度上减小质子调强放疗射程不确定性的负面影响，但是由于只考虑了当前解的最“差”情况，未考虑解在不确定性区间中的分布，导致求出的解过于保守，难以得到鲁棒最优解；本专利的方法建立了同时体现最优性和鲁棒性的射程不确定性模型，提高了模型的合理性和实用性。

3、本发明采用的基于滤波算法的自适应粒子群鲁棒优化方法求解射程不确定性模型更加有效、快速；虽然前期国内外研究也涉及到鲁棒优化，但基本上是基于常规优化方法的改良如L-BFGS-B代码和CPLEX工具包；常规优化方法通用性强，但求解不确定性问题能力一般，导致优化时间长(最长可能数十个小时)，优化效果待提高，制约了该方法在临床的应用；因此，本发明具有极其重要的医学应用价值、而且克服了尚未解决的技术难题。

附图说明

为了便于本领域技术人员理解，下面结合附图对本发明作进一步的说明。

图1为本发明的整体流程图；

图2为本发明的基准几何解剖模型。

具体实施方式

本实施例首先导入质子重离子调强放疗计划信息和准备计算数据，然后根据均值和方差的鲁棒性度量方法构建基于最优性和鲁棒性相结合射程不确定性模型，最后采用基于滤波算法的自适应粒子群鲁棒优化得到鲁棒最优解。本发明不仅可以准确地建立鲁棒优化问题的不确定性模型，而且求解速度快，算法稳定，鲁棒性强，用户可以得到鲁棒最优解。

如图1所示，本发明一种针对射程不确定性的质子重离子调强放疗鲁棒优化方法，具体实现步骤如下：

(1)数据准备

1)导入放疗计划信息

放疗计划信息主要包括患者CT影像、射束参数、加速器硬件参数、器官约束参数、射程不确定性参数等。本项目拟采用标准DICOM RT格式作为输入输出接口，可以导入商用计划系统的放疗计划文件，完成放疗计划后导出质子重离子扫描点权重和剂量信息到DICOM文件，支持的标准文件主要有：CT、DICOM RT-Structure、DICOM RT-Plan、DICOM RT-Dose。

2)计算剂量贡献矩阵

需要计算不同能量不同位置质子重离子扫描点对所有器官采样点的剂量贡献，计算量大，采用相对较成熟的快速解析方法-笔形束剂量计算算法。其核心思想是卷积叠加方法，具体而言就是将待计算位置处质子重离子的通量权重与对应水等效深度的IDD数据作卷积叠加。用公式可以表示为：

Dose(r,z)＝∑_kΦ_k(r,z)·IDD(z)(10)

其中z代表沿束流的深度位置，r代表横向位置，k代表质子重离子笔形束斑的索引，Dose代表某一个位置对应的体素内的总剂量，Φ是质子重离子的通量权重，IDD是预先通过实际测量或者Monte Carlo方法计算好的质子重离子笔形束在水中沿深度方向的积分深度剂量数据。

首先，在射束方向观(BEV)方向上根据质子重离子束斑半高宽进行布点；再将每个束斑划分成19个子射束；然后采用改进的Siddon算法追踪每一子射束，计算射线路径上每个体素对应的水等效深度，查表得到对应的IDD值；根据Fermi-Eyges理论，采用二维校准方式，计算射线穿过的体素所对应多库伦散射的高斯分布标准差σ_MCS；然后根据Soukup等提出的理论公式计算核散射对束流通量横向散射的影响，得到核散射的权重ω_NS以及其高斯分布的标准差σ_NS；考虑到高斯分布的截断误差δ，得到总的通量权重公式为：

Φ(x,y,z)＝(1-ω_NS)·[G(x,y,σ_MCS)+δ_MCS]+ω_NS·[G(x,y,σ_NS)+δ_NS](11)

根据束流监测系统敏感因子S_m，将所有结果按权重叠加，得到最终的剂量dose(x,y,z)：

dose(x,y,z)＝∑_e∑_jS_m·ω_e,j·D_e,j(x,y,z)>

其中x,y,z代表位置坐标，e代表能量层索引，j代表笔形束斑索引，S_m是与束流能量有关的监测系统敏感因子，ω_e,j是能量层e第j个笔束的权重，D_e,j(x,y,z)是能量层e第j个笔束对位置(x,y,z)的剂量。

(2)构建质子重离子调强放疗射程不确定性模型

该实施方式中，记载的(10)、(11)、(13)、(14)、(15)、(16)、(17)表示的是不确定性模型公式和优化算法用到的公式，其中，公式(10)和(11)推导(12)，(13)和(14)推导(15)，(16)和(17)是独立的；表示质子重离子束集合，b(i)表示质子重离子束i的笔形束集合那么所有笔形束数量第i个质子重离子束的第j个笔形束的射程用表示，表示正常射程，表示高出/低于正常射程的最大值。第i个射束的第j个笔形束的射程表示成对于任意的射束采用表示笔形束j∈b(i)的不确定性离散后的可能射程集。笔形束j和里的每个元素是一一对应的，所以任意射程由索引唯一确认。使用表示射程索引元素的向量。为了方便定义是正常射程的索引。

采用分别表示靶区(Target)和危及器官(OAR)的体素。采用w_j表示第j个笔形束的强度，θ_v表示第ν个体素的期望剂量值。是第j个笔形束(j∈b(i))在射程为的情况下对第ν个体素的剂量贡献。第ν个体素受到的总剂量为

可以写成

每个笔形束射程引入不确定性后，目标函数和约束必须依赖可能射程集。可能射程集可以看成是均匀分布，不确定性模型如下：

其中，R为在可能射程集中的射程采样数量，和为变量w_j在可能射程集中的第k和第p个射程，和分别为射程和对应的目标函数值。

(3)基于滤波算法的自适应粒子群鲁棒优化方法

本步骤采用常规粒子群优化流程，即初始化粒子、计算适应度、位置更新、变异。对其中的关键环节采用新的约束处理方法和动态变异策略。采用滤波算法对违反约束的解进行“过滤”处理。此外，为了避免粒子群算法在变量维数高和形式复杂的优化问题中陷入局部最优，采用一种有效的动态变异策略。粒子的速度随着迭代次数动态变化，而粒子的位置随着速度的改变而改变。下面引入多样性函数D和适应度方差函数σ²。

f＝max{1,max{|f(i)-f_avg|}}1≤i≤M

其中，M为种群数量，L为搜索半径。f(i)为第i个粒子的适应度函数值，f_avg为当前粒子群适应度函数值的平均值。

对种群适应度方差函数和种群多样性函数分别设定阈值，当种群在迭代过程中，一旦连续出现种群适应度方差函数值或者种群多样性函致值低于阈值时，则对种群个体粒子的速度进行扰动。但在算法迭代后期，若变异范围过大，则会对粒子进行较大的扰动，影响其搜索。所以需要不断减少粒子速度扰动范围，扰动公式如下所示：

其中，v_i为第i个粒子的速度，t为迭代次数，G^max为迭代的总次数，θ为服从标准正态分布N(0,1)的随机数。

(4)系统性测试

通过基准例题和头颈部肿瘤及胸腹部肿瘤质子调强放疗临床例题各5例，以验证算法的有效性。采用肿瘤放射治疗小组基准几何解剖模型(图2所示)，这是一个用于测试不同治疗方案疗效的标准模型。OAR是半径为1.5厘米的圆，在一个半径为10cm圆形模体的中心。肿瘤(Tumor)的形状像一个“马蹄”，包裹OAR的一半，内径为1.8厘米，外径为4厘米。体素分辨率是1×1毫米，共有2186个Tumor体素，441个OAR体素。该模型是等效水的均匀模体。质子束从三个方向进行照射，从左到右分别是315°、0°和45°。对于每个质子束来说，束斑在轴向(平行)和横向(垂直)方向的间隔均是4毫米，穿过整个肿瘤及其外扩1厘米(各向同性)的体素。三个质子束总共有877个束斑位置(不同能量质子束可以有相同的位置)。

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。