CN102711913B - 放射治疗法的粒子束优化方法及其设备 - Google Patents

Abstract

通过提供标准形式的粒子束集的模型和目标剂量来优化放射剂量。根据模型确定格拉姆矩阵。对目标剂量进行二次抽样，以确定粒子束集的初始强度值。之后，在收敛之前对下述步骤进行迭代。向各强度值添加很小的正值0<ε<<1，以确保所述强度值大于零。将各强度值乘以格拉姆矩阵，以确定积，该积逐元素被划分成标准形式的目标剂量，以确定对应的比。在误差宽容限内如果这些比均接近于1，则输出粒子束集的强度值。否则，在下次迭代之前将强度值乘以该比。

Description

放射治疗法的粒子束优化方法及其设备

技术领域

本发明大体上涉及放射疗法，更具体地说，涉及对优化由粒子束集传送到组织的放射剂量进行优化的方法。

背景技术

在放射剂量优化的最简形式中，剂量矩阵表示m个粒子束集传送到n个体素中的放射目标剂量。此处，体素代表组织中的小立方体(mm³)。非负粒子束权重矢量与该粒子束线性组合，以产生矢量

的累积体素放射剂量。通常，剂量矩阵A相对较大（n>>m），并且稀疏。问题在于在给出目标剂量图案时如何确定确保无剂量不足并且确保过剂量最小的粒子束权重。此为线性规划（LP）：

$\underset{x}{\min }{\left|\right|\mathrm{Ax}-b\left|\right|}_1$ 其中Ax≥b，x≥0        （1.1）

在一台2GHz的奔腾4处理器上，一个小问题(m～10³个粒子束，n～2m个体素）通常可以在大致一秒的时间内被求解。然而，LP的求解时间大体上与问题规模（problem size）的立方成正比，因此目前该方法对于求解更大规模的实际临床问题并不现实。

LP与二阶锥规划（SOCP）方法因它们适合于使剂量传送误差最小的公式而具有吸引力。然而，问题仅限于几千个变量和约束，并且需要若干个小时来求解。因此，问题优选目标函数中的L₂规范（即，在不准许剂量不足的同时使误差平方和最小化）。此为带有不等式的最小平方问题（LSI）：

$\underset{x}{\min }{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2$ 其中Ax≥b，x≥0        （1.2）

能够证明LSI可以优化LP的上限。使过剂量平方和最小给出了稍光滑些的解，但是因过剂量更大而不安全。此LSI是二次规划（QP）的具体的典型实例，其假设矩阵A具有行满秩，并且QP是严格凸的。此处，LSI和QP被互换使用。

快速逼近法

在现有技术中，经常使用避免有关不足剂量的约束并且简单地寻求减少误差平方和

${\left|\right|\mathrm{Ax}-b\left|\right|}_2^2$

或加权误差平方和

${\left|\right|W(\mathrm{Ax}-b)\left|\right|}_2^2$ 的方法。

根据问题的规模，通过梯度、共轭梯度、或者准牛顿法，或者优选的有限记忆布洛伊登-弗莱彻-戈德法布-香农（L-BFGS）法而迭代计算出解。为避免物理上不可能的解，这些方法通常被改进为将任何负粒子束权重x_i设置为零。这样并不确保能找到最优解或避免不足剂量，但似乎对临床应用足够有效。

更原理性的方法（例如梯度投影法（gradient projection））能够确保非负性和最优性，但是会增加处理时间从而变慢。因此，有时使用基于图形处理器单元（GPU）的并行解法来在几秒钟内求解小问题。

发明内容

本发明的实施方式提供了一种优化由粒子束集传送的放射剂量的方法。所述方法使用一组快速乘法更新，以求解非负最小二乘（NNLS）问题以及一类特殊的二次规划（QP）。通过问题转化，所述迭代更新能够求解最小距离问题（LDP）、最小平方不等式问题（LSI）以及一般的二次规划问题（QP）。

意在的应用为放射剂量优化，所述放射剂量优化是一个很大的线性规划（LP），以至于在实践中，所述约束矩阵不能被储存在存储器中。

所述方法利用L₂规范限定LP的L₁目标的上限来得到所述QP，并且通过乘法迭代产生不同的精确解和近似解。

附图说明

图1是根据本发明的实施方式的方法流程图。

具体实施方式

本发明的实施方式提出了一种相当快速的方法，该方法确保平方误差最优性以及放射疗法粒子束的权重非负。

对于“笔形射束（pencil beam）”疗法，大致具有n～10⁶个体素，以及m～10⁴个粒子束。因此，计算时间主要由体素的个数确定，并且理想的是以将基于体素的计算最小化的方式来修订该问题。该方法的步骤能够在本领域中已知的、连接到存储器与输入/输出接口的处理器中进行。并且，与这些步骤对应的装置也能够在处理器中实现。因此根据本发明的设备由上述装置组成。

如图1所示，该方法的输入包括：粒子束模型101、目标剂量102、以及敏感器官（OAR）的位置103。

第一步110确定经标准化的模型和经标准化的目标剂量。

预先确定粒子束模型的格拉姆矩阵其中^T为转置运算符。之后，将相同的转换应用到目标剂量，以得到内积空间中的矢量集的格拉姆矩阵为内积的厄密矩阵。

对目标剂量进行二次抽样（120），以确定粒子束集的初始强度值。

现在，忽略了不等式约束。目标函数被改写为无约束最小二乘（LS）问题

该式比原始LSI小两个数量级。注意在标准化形式中数值精度有损耗，但剂量矩阵A的稀疏性减轻了该损耗。

无约束最小二乘问题以非严格平方误差约束Ax≥b取代了严格不等式约束Ax≈b。由于此问题非负(A,b≥0)并且过约束(n>>m)，所以该LS最优解具有两处不足：一些LS优化粒子束权重的将是负值（x_i<0）；以及一些体素会不足剂量(A_ix<b_i)。

根据本发明的方法完全消除了第一个问题（负粒子束权重）。给出了非负最小二乘（NNLS）问题

$\underset{x}{\min }{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2$ 其中x≥0,        (2.2)

令：

乘法更新130，

是这样的缩并式（contraction）：对于任何正值的h和正定的Q（positive definite Q），该式收敛于根据任何正权重矢量x>0的最优化。该方程将x的每个元素乘以与梯度相关的比。实践中，当全部乘数（在数值误差宽容限内）等于1时，实现收敛（150）。在此情况下，如果收敛，则输出粒子束集的强度值（160）。否则，在下次迭代前将强度值乘以该相乘结果（170）。

在剂量优化的语境中，

${\&ForAll;}_{\mathrm{ij}},Q_{\mathrm{ij}}\&GreaterEqual;0,h_i>0,$

因此该迭代可以被简化为单矩阵矢量积，然后逐粒子束权重进行步骤130的相乘和步骤140的划分，使得适于精细粒子化的平行处理。迭代的代价由格拉姆矩阵Q的非零元（entry）的个数控制，所述非零元的个数为具有无消失交叠性质的粒子束对的个数。

第二个问题（剂量不足）因最小二乘解使误差光滑散布而出现，为了使对邻接的、与健康组织对应的体素过剂量最小化，而造成要治疗的组织容量边界上的体素剂量不足的结果。对于等式（2.2）中的NNLS问题，避免剂量不足仅为非严格约束。有若干加强该约束的、不会从根本上使得问题求解更困难的推导法：

（1）通过向Q添加而显式地惩罚到要被治疗的体素的剂量的方差，其中S是选择这些体素的二元矩阵，I为单位矩阵（identity matrix），1为全部是1的矩阵，以及t为这些体素的个数，以及τ为平衡剂量误差和方差的权重；

（2）将平方误差和重新加权，使得边界体素的剂量误差比其他位置被大大地惩罚，并同样地对治疗区域外的体素减权。换言之，存在体素误差权重矢量w，并令

Q=A^T(diag(w)A以及h=A^T(diag(w)b.

（3）Q=A(diag(w)A以及h=A(diag(w)b；以及

（4）修改被“一裙（skirt）”正值围绕着的目标组织容量，该正值的衰减与粒子束的衰减剖面（profile）近似匹配。

精确解法

如果将避免不足剂量设为严格约束，那么因为需要在体素剂量的大很大的空间内设置运算，并且有亏秩矩阵参与运算，从而使得问题变得更加困难；这些条件都减慢迭代运算器的收敛速度。另一方面，不需要使用推导法。

本节考虑确保无不足剂量的初始LSI/QP的重新构建，并通过乘法更新提供解法。

最小化总剂量

理想的是将服从无不足剂量约束的总剂量||Ax||₁最小化。非负的A在没有一般性的衰减的情况下具有单位列和，从而||Ax||₁=||x||₁。由于||x||₂≥||x||₁，总剂量的上限能够在下述最小距离问题（LDP）中被最小化：

$\underset{x}{\min }{\left|\right|x\left|\right|}_2^2$ 其中Ax≥b，x≥0

已知其具有等效的NNLS形式

其中u，v≥0，其可以通过NNLS迭代而被求解。LDP和NNLS优化因（假定原始约束是一致的，否则等式右边为零）而相关联。注意由于目标函数||x||+2倾向于更加光滑的解，所以LDP允许比等式（1.2）中原始QP更大的过量。并且，在dim(u)+dim(v)=#voxels+#beams为未知的情况下LDP很大。然而，这不允许为追求速度而利用矩阵A的稀疏性。

将LSI旋转为LDS形式

原始LSI能够通过如下地将变量变为z，将其转换成稍更复杂的LDP而被精确地解出：

${\left|\right|\mathrm{Ax}-b\left|\right|}_2^2={\left|\right|z\left|\right|}_2^2.$

之后，LSI能够被改写为等效的LDP。A的“小”QR分解为UR=A，其中U是酉的，而R是在对角线上非负的上三角矩阵。利用变量

的L₂规范保留变换改写QP，来得到最小值

$\underset{z}{\min }{\left|\right|z\left|\right|}_2$ 其中Cz≥d

其中表示拟逆。

对应的NNLS为

其中u≥0，

其被然后相关联地优化。

对应的乘法更新解决了剂量问题，无剂量不足，无负权重，并有使过量平方和最小。然而，由于先前的LDP/NNLS是稀疏的，该NNLS是密集的，所以在每次迭代中需要更多的运算。并且由于该问题被实际转换了两次，因而存在精度的降低。

求解对偶QP问题

给出最初原始形式的严格凸二次规划，

其中Ax≥b，    （3.1）

二次规划的拉格朗日对偶采用以下形式

其中λ≥0    （3.2）

并且h＝b-AP^-1q

该最初优化和对偶优化被x^*=P^-1(q+A^Tλ)相关联，并且该二次优化能够通过NNLS迭代法解出。注意除非x≥0被显式地编辑为约束Ax≥b，否则并不确保x≥0。为了将其限定于等式（1.2）的剂量优化LSI/QP，设

以及

$h\stackrel{\&CenterDot;}{=}{b}^{\&prime;}-0b^{\&prime;},$ 解出，

此处，

以及

以包含约束x≥0。能够通过使用上述的卷积分解法而有效地计算逆矩阵(A^TA)^-1。这里需要注意，因为(A^TA)很可能是条件不足的。

通过迭代法求解经简化的最初QP的方法

通过向这些比的分母添加辅助矩阵矢量积而将放射剂量中的误差缩小到目标剂量的矩形子区，该辅助矩阵矢量积代表针对该子区中的剂量值缺失的额外惩罚。

特别地，考虑到对最初QP的约束集进行限制而只禁止在k≤2m体素处的剂量不足（例如，在边界上的体素，以及离粒子束中心最远的体素）。对于有非负Q、h、A、b以及对于全部i，Q_ii>0,h_i>0的约束的QP，提供辅助矩阵矢量积

λ←diag(A″x)^-1diag(b″)λ

其中b″为对应k个约束体素的子集b，而A″为对应于剂量矩阵A的行的子集。

辅助矩阵矢量积包含对偶变量（dual variable），也被称为影子价格，其在约束为有效时，在解中为正值（否则为零）。在约束优化中，如在本文所述，影子价格代表通过将约束放宽一个单位而得到的优化问题的优化解的目标值的变化。

本方法也可以被用于在任何体素集施以最大或最小累积剂量的严格约束。λ的更新能够沿针对λ的拉格朗日梯度被改变为小步长。

迭代的计算细节

乘法迭代程序具有线性复杂性，并从正锥（positive cone）的某处以线速度收敛。在各次迭代中精度一点点稳定提高。乘法更新具有所有非负约束被隐性增强的优点。然而，可能的是，数值误差会将很小的强度值圆整到零，之后该强度值会“粘”在零边界处。

因此，在步骤130中，在做乘法前，将一个很小的正值，0<ε<<1添加到各强度值上，以便产生正强度值。这确保迭代脱离零边界。逐变量更新能够如信息传递过程那样以矢量平行方式进行，或者在Q_ij≠0时粒子束权重x_i的更新得益于权重x_j的更新的情况下，该更新可以顺序地进行。这使得可以更快一些地收敛。

本应用推荐使用提前终止，因为误差为±2.5%的治疗剂量并不罕见，所以权重x的7比特之后的精度是多余的。

剂量优化的计算细节

可分离卷积

矩阵Q=A^TA的形成通常是这些方案中最耗时的操作。通常情况下，其可以完全被避免。关键点是Q_ij为两个粒子束的内积，因而其可以由卷积来计算。这点可以在最简单的可能的粒子束模型中被利用。该粒子束是标准偏差为σ的xy平面高斯函数和z中布拉格曲线的外积。因此，由p、q、r偏移的两个粒子束的3D卷积被分解为1D卷积的外积

在上述方程中，仅为的逆。因此，整个矩阵Q的元为自卷积的布拉格曲线和两个1D高斯方程的外积。

$Q=Q_B\&CircleTimes;Q_G\&CircleTimes;Q_G.$

该结构可以被应用于通过一个改型的等式

$(C\&CircleTimes;D)x=\mathrm{vec}\left(D^{|}{\mathrm{vec}}_{\mathrm{rows}\left(D\right)}\left(x\right)C\right).$

快速确定Q对x的行为。

对每个外积使用一次该等式，标量乘法的总数就被减少了因子(B_G)²这么多，（(B_G)²即平面中束点的个数）。作为一个优点，对于x+y+z个体素来说，仅使用了矢量Q_B、Q_G、Q_G中的x+y+z个元，而不是格拉姆矩阵Q中的(xyz)²个元。

同理，矢量h=A^tb可以通过沿x轴和y轴将3D目标剂量b与1D高斯方程卷积，然后沿z轴与布拉格曲线卷积而有效地形成。

有界域中的权重误差

注意高斯卷积贯穿整个无限xy平面，这暗示了各处1均为最小平方误差。由于目标剂量几乎在各处均约等为零，无限的跨度有助于目标边界处的滚降。该效果似乎很小，但在某些情况下，积分被限制在非常小的容量

其中Φ()为误差方程。再次重复，此为非负外积。

加权误差

通过将代表额外惩罚来忽略子区中的剂量值的、不在加权误差区域的第二特定公式中的辅助矩阵矢量积加到所述比的分母，而将放射剂量中的误差缩小到目标剂量的矩形子区域。

格拉姆矩阵Q的边界形式在将加权平方误差||diag(w)Ax-diag(w)b||₂最小化时变得有用。为了与标准化的形式一起工作，加权积

并且

不幸的是，w的大部分选择在先使用了外积分解来构建Q。但是在某些实用情况下，该外积可以通过将加权的Q分解为两个矩阵来取得

例如，为了利用w_i=c>1均匀地使（例如，围绕敏感器官的）某个边界框过权重，第二加数简单地为边界积分Q矩阵，如上所述，其与c²-1相乘。为了使少数单独的像素过权重，设

$w_i^{\&prime;}=\sqrt{w_i^2-1},$

并且设W′包含对角（w′）的非零元，从而第二加数可以被写为

A^TW′l的列可以与计算h相同的方法被预先计算，并且其对x的作用被有效地计算为

注意使用外积乘法时，两个加数对x的作用以及结果的总和被分别确定。

发明效果

本发明的实施方式使用二次规划（QP）来求解很大并且复杂的放射治疗的优化粒子束问题。

QP在网络优化、机器学习、物流、控制以及许多其他应用中普遍存在。其在估计量时十分有用，经常以焦耳、时间或钱数计值，其不会变为负值。

QP的解法自WW2起被集中地发展，其中至少有一项诺贝尔奖项与此领域相关（哈利马克维茨-1990）。现今有大量有效的方法，但是问题的规模已经逐渐超出可控范围。

本发明提出了一种新的并且非常简单的固定点（fixpoint），可在线性时间、线性空间中以线性收敛速度求解QP。这变得对优化粒子束非常有用，在优化粒子束的情况下，QP能够具有>10⁴个变量以及>10⁶个约束。

由于它的简单性，该固定点十分快速并可并行处理。之前的QP需要数小时来求解问题。本文的所述方法能够以视频帧的速度求解问题。

虽然本发明已经通过优选实施方式的举例的方式被描述，但应该被理解的是可以在本发明的精神和范围中进行各种其他的变化和修改。因此，所附权利要求旨在涵盖落入本发明的真实精神和范围内的所有变型和修改。

Claims (10)

1.一种优化由粒子束集传送的放射剂量的方法，其中各粒子束为放射治疗粒子束，所述方法包括以下步骤：

提供所述粒子束集的模型、以及目标剂量；

确定包含格拉姆矩阵Q的标准形式的模型和目标剂量h，所述格拉姆矩阵源自所述粒子束集的模型；

对所述目标剂量进行二次抽样，以确定所述粒子束集的初始强度值，并进行下述步骤的迭代，直至收敛为在下述数值误差容限内下述比都接近于1的状态：

向各强度值添加很小的正值0<ε<<1，以确保所述强度值大于零；

将强度值矢量x乘以所述格拉姆矩阵Q，以确定积Qx；

逐元素地，根据式子将所述积Qx与所述标准形式的目标剂量h彼此相加，根据所述积Qx进行划分，以确定与梯度相关的比；以及

确定在数值误差容限内所述比是否都接近于1，如果接近1，则输出所述粒子束集的强度值；如果不接近1，则在下次迭代之前将所述强度值乘以所述比以更新所述强度值，

其中，所述格拉姆矩阵为内积的厄密矩阵，所述格拉姆矩阵的各元素针对各粒子束对，表示由该粒子束对传送的放射剂量的内积。

2.根据权利要求1所述的方法，所述方法还包括以下步骤：

在可能时将所述格拉姆矩阵分解为外积，并且将所述外积的因子直接应用到所述粒子束集的强度值。

3.根据权利要求1所述的方法，所述方法还包括以下步骤：

通过将辅助矩阵矢量积加到所述比的分母，将所述放射剂量的误差缩小到所述目标剂量的矩形子区域中，所述辅助矩阵矢量积表示对于所述子区域内的剂量值缺失的额外惩罚。

4.根据权利要求1所述的方法，所述方法还包括以下步骤：

通过将辅助矩阵矢量积加到所述比的分子，对所述放射剂量的任意子区域加以严格放射限制，所述辅助矩阵矢量积表示确保严格放射限制的影子价格。

5.根据权利要求1所述的方法，所述方法还包括以下步骤：

将所述模型转换成对偶形式，并且对该对偶模型应用二次抽样步骤，来计算针对所述目标剂量规定的各约束的影子价格。

6.一种优化由粒子束集传送的放射剂量的设备，其中各粒子束为放射治疗粒子束，所述设备包括：

用于提供所述粒子束集的模型以及目标剂量的装置；

确定装置，用于确定包含格拉姆矩阵Q的标准形式的模型和目标剂量h，所述格拉姆矩阵源自所述粒子束集的模型；

用于对所述目标剂量进行二次抽样，以确定所述粒子束集的初始强度值，并且进行下述步骤的迭代，直至收敛为在下述数值误差容限内下述比都接近于1的状态的装置：

向各强度值添加很小的正值0<ε<<1，以确保所述强度值大于零；

将强度值矢量x乘以所述格拉姆矩阵Q，以确定积Qx；

逐元素地，根据式子将所述积Qx与所述标准形式的目标剂量h彼此相加，根据所述积Qx进行划分，以确定与梯度相关的比；以及

确定在数值误差容限内所述比是否都接近于1，如果都接近1，则输出所述粒子束集的强度值；如果不是都接近1，则在下次迭代之前将所述强度值乘以所述比以更新所述强度值，

其中，所述格拉姆矩阵为内积的厄密矩阵，所述格拉姆矩阵的各元素针对各粒子束对，表示由该粒子束对传送的放射剂量的内积。

7.根据权利要求6所述的设备，所述设备还包括：

用于在可能时将所述格拉姆矩阵分解为外积，并且将所述外积的因子直接应用到所述粒子束集的强度值的装置。

8.根据权利要求6所述的设备，所述设备还包括：

用于通过将辅助矩阵矢量积加到所述比的分母而将所述放射剂量的误差缩小到所述目标剂量的矩形子区域中的装置，所述辅助矩阵矢量积表示对所述子区域内的剂量值缺失的额外惩罚。

9.根据权利要求6所述的设备，所述设备还包括：

用于通过将辅助矩阵矢量积加到所述比的分子，对所述放射剂量的任意子区域加以严格放射限制的装置，所述辅助矩阵矢量积表示确保严格放射限制的影子价格。

10.根据权利要求6所述的设备，所述设备还包括：

用于将所述模型转换成对偶形式，并且对该对偶模型应用二次抽样步骤，来计算针对所述目标剂量规定的各约束的影子价格的装置。