放射治疗的出射射束强度的计算方法和装置的制造方法
【技术领域】
[0001] 本发明涉及放射治疗领域，尤其是涉及放射治疗的出射射束强度的计算方法和装 置。
【背景技术】
[0002] 放射治疗是利用放射线治疗肿瘤的一种局部治疗方法。放射线包括放射性同位素 产生的α、β、γ射线和各类X射线治疗机或加速器产生的X射线、电子线、质子束及其他 粒子束等。
[0003] 在现代放射治疗计划中广泛使用二次规划问题模型。二次规划是非线性规划中 的一类特殊数学规划问题。在各种放射治疗类型，例如容积调强放射治疗（Volumetric ArcRadiationTherapy,VMAT)和调强放射治疗（IntensityModulatedRadiation Therapy，頂RT)中，二次规划问题可用于通量图优化、直接子野优化，来计算需要的出射射 束强度。
[0004] 为了通过求解大规模二次规划问题来计算需要的出射射束强度，一种主流方法是 Mehrotra预测-修正（predictor-corrector)格式的内点法，该方法涉及求解两次大规模 的线性方程组。直接法求解此类线性方程组，如高斯消去法（LU分解）、对称不定分解（LDLT 分解）和传统迭代法（Jacobi,Gauss-Seidel，S0R等）直接求解的代价是昂贵的。
[0005] 无等式约束情形下，线性方程组可以写成拟定系统。此时，可以通过Vanderbei对 拟定方程组求解的技巧，使用块Cholesky分解进行加速。不利的是，此方法在带等式约束 的问题中，无法保证子线性系统仍然是正定的。尤其是当迭代进行了若干步时，线性系统接 近奇异，此时的Cholesky分解由于数值不稳定性会导致失败。使得Vanderbei的拟定算法 无法进行。
【发明内容】
[0006] 本发明所要解决的技术问题是提供一种放射治疗的出射射束强度的计算方法和 装置，以缓解放射治疗的优化方法中二次规划问题求解速度慢以及容易失败的问题。
[0007] 为解决上述技术问题，本发明提供了一种放射治疗的出射射束强度的计算方法， 包括以下步骤：
[0008] a.输入二次规划问题，该问题表述为：
[0012] 其中，t表示每个三维体素上的剂量分布，向量d是平均处方剂量向量，向量dmax 是最大处方剂量向量，向量dmin是最小处方剂量向量，系数向量wd,wmax,wmin分别是平均处方 剂量权重、最大处方剂量权重、最小处方剂量权重，向量V是每一个控制点处的出射射束强 度，向量V_是每一个控制点处的出射射束的最大强度；其中向量d，d_，d_，wd，w_，w_的 向量长度是CT图像离散化之后三维体素的个数，向量V，7_的向量长度为计划控制点的个 数；矩阵D是剂量沉积矩阵，该矩阵的行数为三维体素的个数，其列数为控制点的个数；记 号Θ是Hadamard乘积；记号（·）+是取正部函数，即s+=max(s, 0);
[0013] b.对该二次规划问题进行预求解；
[0014] c.根据预求解结果判断该二次规划问题是否正常，如果是则进入步骤d，否则进 入步骤f;
[0015] d.对该二次规划问题进行Mehrotra预测修正迭代，其中优先使用Vanderbei策略 求解拟定系统，当分解遇到不正定的情形，则退回到使用对称不定分解求解增广系统；
[0016] e.后处理；
[0017] f.输出该二次规划问题的解，其中当该二次规划问题正常时，该二次规划问题的 解为该向量v的值，当该二次规划问题不正常时，该二次规划问题无解。
[0018] 在本发明的一实施例中，在该步骤d中，当退回到对称不定分解求解增广系统之 后，在后续迭代中不再返回到使用Vanderbei策略求解拟定系统。
[0019] 在本发明的一实施例中，该步骤d包括以下步骤：在步骤dl，根据该步骤b探测稠 密行与稠密列的结果对剂量沉积矩阵D进行排序；在步骤d2,生成增广系统的固定部分； 在步骤d3,初始化不定标记为正定；在步骤d4,判断不定标记，如果标记为不定，则进入步 骤d5,否则进入步骤d6 ;在步骤d5,使用对称不定分解增广系统，记录分解因子，跳到步骤 dlO;在步骤d6,使用Vanderbei策略求解拟定系统；在步骤d7,判断求解拟定系统分解是 否成功，如果是则进入步骤d8,否则进入步骤d9 ;在步骤d8,记录求解拟定系统的分解因 子；在步骤d9,更新不定标记为不定；在步骤dlO,根据不定标记，计算仿射步，得到关于解、 对偶变量、拉格朗日乘子的更新；在步骤dll，执行线搜索，找出仿射步长；在步骤dl2,根据 仿射步骤更新残量；在步骤dl3,判断不定标记，如果标记为不定，则进入步骤dl4,否则进 入步骤dl5 ;在步骤dl4,使用对称不定分解因子求解增广系统，进入步骤dl6 ;在步骤dl5， 使用Vanderbei策略求解拟定系统；在步骤dl6,清空分解因子内存；在步骤dl7,执行线搜 索，找出修正步长；在步骤dl8,更新当前解、对偶变量、拉格朗日乘子向量；在步骤dl9,判 断停机准则是否满足，如果是则进入步骤d20,否则返回步骤d04 ;在步骤d20,停机。
[0020] 在本发明的一实施例中，该步骤d中，使用Vanderbei策略求解拟定系统：
[0021]
其中矩阵AjpAD是不等式约束条件A和等式约束条件 A^i过稠密行和稠密列重新排列得到的矩阵，向量x、ys、yD分别是二次规划主问题的变量、 对偶问题稀疏列部分所对应的变量、对偶问题稠密列部分所对应的变量，向量分别 是三个方程的残量，且包括如下步骤：
[0022] 步骤d21，对Q_ASTASAS=LLT进行Cholesky尝试分解，如果分解不成功，则退出拟 定系统求解，并返回不定标记；
[0023] 步骤d22,求解三角方程形成辅助矩阵Z=L
[0024] 步骤d23，Cholesky尝试分解AD-ZTZ=KKT，如分解不成功，则同步骤d21，退出拟 定系统求解并返回不定标记；
[0025] 步骤d24，计算辅助变量》= & - 4 八
[0026] 步骤d25,计算辅助变量π;,,
[0027] 步骤d26,计算yD=ΚΤΚV;
[0028] 步骤d27,计算解向量χ=ΚΤΚ1 (u_ADyD);
[0029] 步骤d28,计算ys=Λ1 (rs_Asx)。
[0030] 本发明还提出一种放射治疗的出射射束强度的计算装置，包括：
[0031] 用于输入二次规划问题的模块，
[0032] 其中该二次规划问题表述为：
[0036] 其中，t表示每个三维体素上的剂量分布，向量d是平均处方剂量向量，向量dmax 是最大处方剂量向量，向量dmin是最小处方剂量向量，系数向量wd,wmax,wmin分别是平均处方 剂量权重、最大处方剂量权重、最小处方剂量权重，向量v是每一个控制点处的出射射束强 度，向量v_是每一个控制点处的出射射束的最大强度；其中向量d，d_，d_，wd，w_，w_的 向量长度是CT图像离散化之后三维体素的个数，向量V，7_的向量长度为计划控制点的个 数；矩阵D是剂量沉积矩阵，该矩阵的行数为三维体素的个数，其列数为控制点的个数；记 号Θ是Hadamard乘积；记号（·）+是取正部函数，即s+=max(s, 0);
[0037] 用于对该二次规划问题进行预求解的模块；
[0038] 用于根据预求解结果判断该二次规划问题是否正常的模块；
[0039] 用于该二次规划问题正常时对该二次规划问题进行Mehrotra预测修正迭代的模 块，其中优先使用Vanderbei策略求解拟定系统，当分解遇到不正定的情形，则退回到使用 对称不定分解求解增广系统；
[0040] 用于后处理的模块；
[0041] 用于输出该二次规划问题的解模块，其中当该二次规划问题正常时，该二次规划 问题的解为该向量v的值，当该二次规划问题不正常时，该二次规划问题无解。
[0042] 与现有技术相比，本发明将Vanderbei的策略作为优先策略，如果方程能够使用 拟定系统技巧求解时，就使用该方法，通过对稠密行与稠密列的重新排序，保证了对于良态 问题的速度优势；当分解遇到不正定的情形，本发明退回到使用对称不定分解求解增广系 统的策略。这使得求解器能够处理一些病态问题，并使得在求解器迭代到边界时，降低数值 不稳定性对求解器的影响。
【附图说明】
[0043] 图1是本发明一实施例的放射治疗的出射射束强度的计算方法流程图。
[0044] 图2A、2B是图1中的Mehrotra预测修正迭代步骤的流程图。
【具体实施方式】
[0045] 为让本发明的上述目的、特征和优点能更明显易懂，以下结合附图对本发明的具 体实施方式作详细说明。
[0046] 在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以 采用其它