一种剂量分布确定和放射治疗计划优化的方法及设备

摘要

本公开涉及放射治疗计划中剂量分布的确定方法和装置，包括：将射束和感兴趣区域分别离散化为多个射束元和多个体素，确定每个射束元对每个体素的剂量贡献值；根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布；根据所述连续分布计算当前叶片位置处的剂量分布。本公开还涉及放射治疗计划优化方法和装置。通过本发明的方法或装置，可以将多叶光栅的叶片移动到任何最优解的位置，不受限于射束元的分辨率，并且本发明的方法或装置，计算量小，速度快。

说明书

技术领域

本申请一般涉及放射治疗领域，尤其涉及放射治疗计划优化方法及设备、以及放射治疗计划中剂量分布的确定方法及设备。

背景技术

多叶准直器(MLC)是用来产生适合形状的辐射野的设备，其广泛应用于医学领域。多叶准直器通常通过每个叶片的移动来达成射野动态或静态成形。放射治疗计划优化中实现多叶准直器的功能的关键包括了子野形状和强度的确定。

现有技术中，在确定子野形状和强度时存在以下挑战。

放射治疗优化问题中叶片的可连续移动性实现困难。在制定逆向优化计划时，需要给出每个子野形状中MLC叶片的移动位置，通常通过拼凑一些射束元(beamlet)剂量矩阵来实现。这时叶片只能定位到射束元的边缘处，叶片移动的最小距离为一个射束元，因此叶片的移动连续性被射束元的分辨率所限制；分辨率过低，导致叶片移动的最小距离过大，从而使得叶片移动位置的可调节性降低，达不到预期的优化效果；而分辨率过高，则会使射束元的剂量矩阵数据量急剧增大，并使计算时间显著延长。

另外，射束元剂量矩阵的数据量过大。在VMAT计划优化过程中，首先要计算射束元剂量矩阵，用来存放该射束元射束元附近采样点或称体素(voxel)的剂量沉积情况(cGy/MU)。在一些多器官病例中为了较好的控制各器官上的剂量分布，体素的数目选取比较多，这就使得即使在射束元的分辨率较低的情况下也会产生大量数据，导致消耗大量的内存资源，并占用计算时间。此外，计算速度慢也是现有技术的问题之一。在放射治疗计划的优化过程中，需要计算当前的体素上的剂量分布情况，并根据当前剂量分布情况，来确定下一步的叶片移动方向；每次更新剂量分布情况或者改变叶片位置都需要对子野剂量做累计处理，从而使计算量增大，降低优化速度。

例如，Jiang、Men和Jia于2011年在专利申请号为US13/695,015的美国专利申请中提到，可利用列生成(columngeneration)方法，在两个问题中分别求解子野形状和强度。该方法对于求解大规模优化问题具有一定效果，但该方法变量空间大，消耗较多计算机资源，并且在优化过程中将子野权重和形状分别优化，并不是在全局范围内寻找局部最优解。

又如，美国专利US2013/0077751A1通过构建网络模型，并将叶片移动约束加入其中，该模型更加注重叶片移动的连续性，而相对忽视优化的效果，且其子野形状和子野权重也是分别优化的。

在以上两篇专利中，叶片的可移动位置都是离散的，其可移动位置仍受限于射束元，故其精细调整能力将会受到很大程度上的限制。

发明内容

根据本公开的一方面，一种放射治疗计划中剂量分布的确定方法包括将射束和感兴趣区域分别离散化为多个射束元和多个体素，确定每个射束元对每个体素的剂量贡献值；根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布；以及根据所述连续分布计算当前叶片位置处的剂量分布。

根据该方面的一实施例，根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布进一步包括对每个体素，将所述每个射束元对所述体素的剂量贡献值按行存储，所述行与叶片的移动方向平行；对每行的射束元，计算每个射束元依次叠加后，所述体素上的累计剂量随射束元的分布；以及使用拟合函数对所述累计剂量随射束元的分布进行拟合从而得到所述体素上的累计剂量随射束元的连续分布。所述拟合函数可为一维函数，所述累计剂量随射束元的连续分布可为累计剂量随位置变化的曲线。替换地，所述拟合函数可为二维函数，所述累计剂量随射束元的连续分布可为累计剂量随位置和行数变化的曲面。

根据该方面进一步的实施例，在拟合得到所有体素上的累计剂量随射束元的连续分布之后，还包括对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正。根据又一实施例，可利用插值方法对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正。根据进一步的实施例，对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正进一步包括：建立放射治疗计划的优化模型，利用拟合得到的每个体素上的累计剂量随射束元的连续分布计算所述优化模型的初始解，根据所述初始解所在范围对每个体素上的累计剂量随射束元的连续分布进行修正。进一步，根据所述初始解所在范围对每个体素上的累计剂量随射束元的连续分布进行修正可包括对每个体素，根据所述初始解确定需要修正的射束元的位置范围，利用插值方法对需要修正的射束元的位置范围内的累计剂量进行修正。

根据本公开的另一方面，一种放射治疗计划优化方法包括将射束和感兴趣区域分别离散化为多个射束元和多个体素，确定每个射束元对每个体素的剂量贡献值；根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布；建立放射治疗计划的优化模型；以及利用每个体素上的累计剂量随射束元的连续分布计算所述优化模型的最优解，并作为放射治疗计划的优化结果。

根据该方面的一实施例，根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布进一步包括对每个体素，将所述每个射束元对所述体素的剂量贡献值按行存储，所述行与叶片的移动方向平行；对每行的射束元，计算每个射束元依次叠加后，所述体素上的累计剂量随射束元的分布；以及使用拟合函数对所述累计剂量随射束元的分布进行拟合从而得到所述体素上的累计剂量随射束元的连续分布。所述拟合函数可为一维函数，所述累计剂量随射束元的连续分布可为累计剂量随位置变化的曲线。替换地，所述拟合函数可为二维函数，所述累计剂量随射束元的连续分布可为累计剂量随位置和行数变化的曲面。

根据该方面的又一实施例，在拟合得到所有体素上的累计剂量随射束元的连续分布之后，还包括对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正。例如，可利用插值方法对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正。

根据进一步的实施例，对拟合得到的所述体素上的累计剂量随射束元的连续分布进行修正包括利用拟合得到的每个体素上的累计剂量随射束元的连续分布计算所述优化模型的初始解，根据所述初始解所在范围对每个体素上的累计剂量随射束元的连续分布进行修正。

根据该方面的又一实施例，根据所述初始解所在范围对每个体素上的累计剂量随射束元的连续分布进行修正包括对每个体素，根据所述初始解确定需要修正的射束元的位置范围，利用插值方法对需要修正的射束元的位置范围内的累计剂量进行修正。进一步的实施例还可包括根据修正后的所述体素上的累计剂量随射束元的连续分布计算最优解情况下的剂量分布。

本公开还包括与上述各方面相应的设备。

附图说明

图1是根据本公开一实施例的射束元剂量分布的示例。

图2是根据本公开一实施例的射束元剂量累积分布的示例。

图3是根据本公开一实施例的射束元矩阵单行拟合结果的示例。

图4是根据本公开一实施例的对累计射束元区间进行分段插值的结果示例。

具体实施方式

现在参照附图描述各个方面。在以下描述中，出于解释目的阐述了众多具体细节以提供对一个或多个方面的透彻理解。然而，明显的是，没有这些具体细节也可实践此种(类)方面。

本公开涉及放射治疗计划中剂量分布的确定方法及设备、放射治疗计划优化方法及设备。例如，本公开的各方面可以涉及但不限于以下一者或多者：通过曲线拟合的方法来构造叶片可连续移动的优化模型；用拟合所得的参数来代替射束元矩阵，从而大幅减小了数据规模；用拟合的方法来确定累计剂量的连续分布，从而增强了优化问题的平滑性，使得可以收敛到理想的解，不受限于射束元的分辨率；采用解析方法来快速寻找最优解，对于拟合误差采用例如插值法来修正，并可基于修正来重新优化，兼顾了计算速度和精度。

以下给出本公开的各种具体实施例。应当了解的是，本公开的具体实施例是以示例而非限定的方式给出的，本公开的保护范围不被限定于任何特定实施例。确切而言，本公开所请求保护的范围仅由权利要求书来限定。例如，本公开涉及但不限于多叶准直器叶片位置连续优化。如本领域普通技术人员可理解的，本公开同样可适应于其他相关应用。

根据本公开的方面，可以将射束和感兴趣区域分别离散化为多个射束元(beamlet)和多个体素(voxel)，确定每个射束元对每个体素的剂量贡献值；并且根据所述剂量贡献值确定每个体素上的累计剂量随射束元的连续分布。

例如，确定每个射束元对每个体素的剂量贡献值可包括计算射束元剂量矩阵。可进一步将计算出的射束元剂量矩阵按行存储。本领域技术人员可知，也可使用其他存储方式，诸如按列存储等，而不会脱离本公开的范围。根据一实施例，行与叶片的移动方向可为平行的。

图1根据一示例性实施例，示出了一个体素在一行中的各射束元的剂量分布情况，其中x轴可为该行处从左到右的射束元的位置，y轴可为射束元对于该体素的剂量贡献值(cGy/MU)，并可记第i个射束元对应的剂量贡献值为b_i。

根据进一步的示例性而非限定性的实施例，本公开的优化方法可包括从左到右对射束元的剂量贡献值进行累加，以得到该体素上的累计剂量分布。如本领域普通技术人员所能认识到的，本公开并不限于从左到右地进行累加，而是可以有其他累加方向和/或方式。例如，对每行的射束元，可以计算每个射束元的剂量贡献值依次叠加后，该体素上的累计剂量随射束元位置的分布。

根据一示例，该累加的结果可如图2所示，其中x轴为从左到右的射束元位置，y轴为对应位置处射束元的剂量贡献值累加之后的剂量累计值。例如，y_i可为第i个射束元的剂量贡献值累加之后的该体素上的剂量累计值，其可例如如下式(1)所示地来计算：

y_i＝y_i-1+b_i(1)

根据进一步的示例性而非限定性的实施例，根据所述连续分布计算当前叶片位置处的剂量分布可以进一步包括使用拟合函数对累计剂量随射束元的分布进行拟合从而得到所述体素上的累计剂量随射束元的连续分布，例如可包括选取拟合函数对该体素上的累计剂量分布进行拟合。根据非限定性示例，可选取S型拟合函数对该体素上的累计剂量分布进行拟合。该拟合函数的选取方法有很多，作为示例性而非限定性示例，比如有Logistics曲线、Gompertz曲线、Bertalanffy曲线、曲正切函数(tanh(x))、多项式函数等。

下面，为便于解说，结合Logistics函数为例对本公开的实施例进行说明。应用其他曲线来拟合的示例在本公开的基础上对于本领域普通技术人员是可现成的，故而被包括在本发明的范围中。

例如，典型的Logistics函数可以写为：

$>y=\Phi \left(x\right)=\frac{K}{1+e^{a+bx}}---\left(2\right)$>

其中a、b和K可为待定参数。

例如，为了利用此类Logistics函数对例如该体素上的累计剂量分布进行拟合，可通过对待定参数a、b和K中的任何一者或多者进行调节以使累计剂量分布与Logistics函数达到最佳拟合。例如，可以使用最小均方误差作为拟合函数的目标。例如，均方误差最小化可以写为：

>

其中N代表射束元的个数，由最优化理论可知，通过例如分别对a、b、K求偏导，则在以下方程组(4)成立时所对应的a、b、K值为该优化问题的最优解。

$>\left(\begin{array}{c}\frac{\partial H}{\partial K}=\underset{i=1}{\overset{N}{\Sigma }}(\frac{K}{1+e^{a+{\mathrm{bx}}_i}}-y_i)\frac{1}{1+e^{a+{\mathrm{bx}}_i}}=0\\ \frac{\partial H}{\partial b}=\underset{i=1}{\overset{N}{\Sigma }}(\frac{K}{1+e^{a+{\mathrm{bx}}_i}}-y_i)\frac{-K}{{(1+e^{a+{\mathrm{bx}}_i})}^2}\left(x_i{e}^{a+{\mathrm{bx}}_i}\right)=0\\ \frac{\partial H}{\partial a}=\underset{i=1}{\overset{N}{\Sigma }}(\frac{K}{1+e^{a+{\mathrm{bx}}_i}}-y_i)\frac{-K}{{(1+e^{a+{\mathrm{bx}}_i})}^2}\left(e^{a+{\mathrm{bx}}_i}\right)=0\end{array}\right)---\left(4\right)$>

本领域普通技术人员可知，若选用了其他拟合函数，则可相应地对其参数进行调节以使累计剂量分布函数与所选用的拟合函数达到最佳拟合。

根据一示例性实施例，对单行处射束元矩阵的拟合结果如图3所示。对于每个体素而言，原有每行上的射束元矩阵中的N个数据可以用拟合后的3个参数K、a、b来代替，这样就可以大幅度减少射束元矩阵所占有的数据量，减少存储和/或传输所占用的资源。不仅如此，还可以进一步利用解析方法来加快目标函数和/或变量导数的计算。

在上例中，拟合函数可为一维函数，并且累计剂量随射束元的连续分布可为累计剂量随位置变化的曲线。在另一实施例中，该拟合函数可为二维函数，累计剂量随射束元的连续分布可为累计剂量随位置和行数变化的曲面。

例如，在一示例性而非限定性的示例中，例如，仍以Logistics函数为例，二维的拟合函数可以写为：

$>y=\Phi (x,r)=\frac{K}{(1+e^{a_1+b_{1}x})(1+e^{a_2+b_{2}r})}---\left(10\right)$>

其中x代表叶片位置，r表示叶片所在的行。这样数据量可以得到进一步的压缩，其代价是可能会有相应的精度降低。

在确定了每个体素上的累计剂量随射束元的连续分布后，在一方面，可根据所述连续分布计算当前叶片位置处的剂量分布。

根据本公开的另一方面，在根据剂量贡献值确定每个体素上的累计剂量随射束元的连续分布后，可进一步建立放射治疗计划的优化模型，并且利用每个体素上的累计剂量随射束元的连续分布计算所述优化模型的最优解，并作为放射治疗计划的优化结果。

根据本公开的一方面，一种建立放射治疗计划的优化模型的解析方法可包括设优化的总目标函数为F，则F为关于子野叶片位置x和子野权重w的函数：

F＝F(x,w)＝F(d(x,w))(5)

其中d是d_j组成的向量，d_j为第j个体素(或称采样点)处的剂量值，其表达式可以为例如：

$>d_j(x,w)=\underset{k=1}{\overset{K}{\Sigma }}{w}_k\underset{r\in R}{\Sigma }[{\Phi }_r^{j,k}\left(x_{right}^{k,r}\right)-{\Phi }_r^{j,k}\left(x_{left}^{k,r}\right)]---\left(6\right)$>

K为子野的个数；w_k为第k个子野的权重，k的范围从1到K；R为MLC叶片的行数；r为行，其范围为R；和分别为第k个子野第r行的左右叶片的位置；为第k个子野第r行第j个体素的拟合函数。

于是计算当前叶片位置处的剂量分布情况只需代入上式即可，而无需重新对射束元剂量矩阵的值进行累加。由此，目标函数的计算可得以加快。用拟合的方法来确定累计剂量的连续分布可增强优化问题的平滑性，使得可以收敛到理想的解。

另一方面，使用基于梯度的优化方法进行放射治疗计划时需要确定每个变量的导数，利用上式也可以很方便的求得：

$>\frac{\partial F}{\partial x}=\frac{\partial F}{\partial d}·\frac{\partial d}{\partial x}$>

(7)

$>\frac{\partial F}{\partial w}=\frac{\partial F}{\partial d}·\frac{\partial d}{\partial w}$>

其中

$>\frac{\partial d_j}{\partial x^{k,r}}=w_k{\Phi }^{k,j^{\prime }}\left(x^{k,r}\right)·flag=w_k·flag·\frac{-bK·e^{a+{\mathrm{bx}}^{k,r}}}{{(1+e^{a+{\mathrm{bx}}^{,k,r}})}^2}---\left(8\right)$>

其中flag为标志，当x^k,r为左叶片是flag＝-1，x^k,r为右叶片时flag＝1。

$>\frac{\partial d_j}{\partial w_k}=\underset{r\in R}{\Sigma }[{\Phi }_r^{j,k}\left(x_{right}^{k,r}\right)-{\Phi }_r^{j,k}\left(x_{left}^{k,r}\right)]---\left(9\right)$>

由此，可以很容易获得当前变量处的导数，以用于例如但不限于基于梯度的优化方法等。

通过例如上述拟合来得到的剂量累计曲线与真实情况可能存在一定的误差。故在利用例如上述的解析方法得到初步解后，需要利用精确的射束元剂量进行修正。作为示例性而非限定性例子，可利用各种插值方法对拟合得到的体素上的累计剂量随射束元的连续分布进行修正。例如，可以使用多项式插值，诸如拉格朗日(Lagrange)插值、牛顿(Newton)插值；Hermit插值；分段插值；样条插值等等。

根据示例性实施例，在建立放射治疗计划的优化模型后，可以利用通过例如前述拟合得到的每个体素上的累计剂量随射束元的连续分布计算该优化模型的初始解，然后根据该初始解所在范围对每个体素上的累计剂量随射束元的连续分布进行修正。

根据一示例性实施例，例如，可首先根据当前叶片位置，来确定需要精确修正的射束元位置范围。例如对每个体素，可根据基于拟合所得到的初始解来确定需要修正的射束元位置范围。接着，可利用插值方法对需要修正的射束元位置范围内的累计剂量进行修正。根据一示例性实施例，可以以各射束元位置为节点，对该区间进行插值，从而可以得到一个经过各射束元累计剂量值的精确剂量变化曲线。

在一个示例性而非限定性的示例中，利用hermit插值来进行修正。图4中示例性地示出了以300为中心，对区间260-360利用分段3次hermit插值后的到的结果，其中x轴的单位可为毫米(mm)。

如从图4中所可见，其中‘*’处为真实射束元位置，实线为插值后的结果。利用该方法修正后，可以得到叶片在特定位置处的精确剂量。

根据一进一步的示例性实施例，可以基于此来重复前述解析方法的优化过程，从而可以对初步的结果加以修正。例如，可根据修正后的该体素上的累计剂量随射束元的连续分布，来计算最优解情况下的剂量分布。通过采用拟合法来快速寻找最优解，对于拟合误差采用例如插值法来修正，并基于修正来重新优化，兼顾了计算速度和精度。

根据本公开一示例性而非限定性的实施例的流程概括可包括：

(1)计算射束元，对每个体素将剂量沉积例如按行累加，得到该行上离散的累计剂量沉积点，对这些位置上的点利用例如‘S’型曲线进行拟合，记录拟合后的参数用于后续计算；

(2)采用例如基于梯度的最优化方法，求解优化模型的最优解，在求解剂量累计和导数时可以利用已有的拟合参数快速计算；

(3)以(2)中得到的叶片位置作为中心点，选取一个小的区间作为修正区域，利用插值方法对(2)中的优化结果进行局部修正，得到精确剂量分布情况，并修正局部最优解。

类似地，使用例如如前所述的二维拟合函数以解析方法进行优化后，同样可以类似地利用插值方法来进行修正。

根据本公开的方法和设备实现了各种优势，包括但不限于例如以下一者或多者：

1.减少了射束元的数据量。根据本公开的方法和设备用拟合参数来代替射束元剂量矩阵，使数据量明显降低，从而减少了存储/传输所需占用的资源量并大大提高了效率。

2.利用拟合参数，通过解析方法计算当前形状下的剂量分布情况和每个叶片所处位置的导数，与剂量叠加方法相比，大大加快了优化的迭代速度。

3.增加的优化模型的平滑性，使优化问题的收敛速度加快并使问题更加逼近最优解。

4.采用了两步优化策略，第一步利用拟合方法快速定位出局部最优解的粗略位置；第二步，利用插值法弥补拟合方法中的误差，精确修正局部最优解。在不丢失优化精度的前提下，大幅度提高优化速度。

本领域普通技术人员应理解，本公开的有益效果并非由任何单个实施例来全部实现。各种组合、修改和替换均为本领域普通技术人员在本公开的基础上所明了。

此外，术语“或”旨在表示包含性“或”而非排他性“或”。即，除非另外指明或从上下文能清楚地看出，否则短语“X采用A或B”旨在表示任何自然的可兼排列。即，短语“X采用A或B”藉由以下实例中任何实例得到满足：X采用A；X采用B；或X采用A和B两者。另外，本申请和所附权利要求书中所用的冠词“一”和“某”一般应当被理解成表示“一个或多个”，除非另外声明或者可从上下文中清楚看出是指单数形式。

各个方面或特征将以可包括数个设备、组件、模块、及类似物的系统的形式来呈现。应理解和领会，各种系统可包括附加设备、组件、模块等，和/或可以并不包括结合附图所讨论的全部设备、组件、模块等。也可以使用这些办法的组合。

结合本文所公开的实施例描述的各种说明性逻辑、逻辑块、模块、和电路可用通用处理器、数字信号处理器(DSP)、专用集成电路(ASIC)、现场可编程门阵列(FPGA)或其它可编程逻辑器件、分立的门或晶体管逻辑、分立的硬件组件、或其设计成执行本文所描述功能的任何组合来实现或执行。通用处理器可以是微处理器，但在替换方案中，处理器可以是任何常规的处理器、控制器、微控制器、或状态机。处理器还可以被实现为计算设备的组合，例如DSP与微处理器的组合、多个微处理器、与DSP核心协同的一个或多个微处理器、或任何其它此类配置。此外，至少一个处理器可包括可作用于执行以上描述的一个或多个步骤和/或动作的一个或多个模块。

此外，结合本文中所公开的方面描述的方法或算法的步骤和/或动作可直接在硬件中、在由处理器执行的软件模块中、或在这两者的组合中实施。

本公开中通篇描述的各种方面的要素为本领域普通技术人员当前或今后所知的所有结构上和功能上的等效方案通过引述被明确纳入于此，且意在被权利要求书所涵盖。此外，本文所公开的任何内容都并非旨在贡献给公众——无论这样的公开是否在权利要求书中被显式地叙述。