具体实施方式

下面结合附图详细说明本发明的优选实施方式。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为了达到本发明的目的，一种放疗计划强度分布优化方法、装置及放射治疗设备的其中一些实施例中，如图1所示，该放疗计划强度分布优化方法，包括以下步骤：

S1：通过神经网络，将PB剂量映射到MC剂量上；

S2：求解MC剂量相对于PB剂量的梯度，得到基于MC算法的沉积矩阵；

S3：将基于MC算法的沉积矩阵带入FMO中，实现放疗计划强度分布优化。

值得注意的是，S1中，PB剂量和MC剂量可以但不限于通过计划系统TPS计算得到。

其中PB剂量的计算方法可根据公式或试验结果计算主射线和散射线在组织内的剂量沉积，MC剂量的计算方法是模拟数十亿个粒子在组织内的剂量沉积。

为了进一步地优化本发明的实施效果，在另外一些实施方式中，其余特征技术相同，不同之处在于，S1具体包括以下步骤：

S1.1：数据准备，选取N1个病人CT，每个病人添加N2个射野，分别计算每个射野的PB剂量和MC剂量，形成N1*N2组数据，每组数据均包括对应的CT、PB剂量和MC剂量；

S1.2：数据划分，将N1*N2组数据按比例划分为训练集测试集；

S1.3：创建卷积神经网络；

S1.4：采用训练集数据训练创建好的卷积神经网络；

S1.5：在训练好的卷积神经网络上测试测试集数据，得到预测剂量，计算预测剂量和MC剂量的gamma通过率。

进一步，卷积神经网络包括：输入层、融合层、卷积层、激活层、池化层、上采样层、输出层，如图6所示。

在具体地实施例中，PB剂量到MC剂量的映射方法，具体可以分为如下步骤，如图5所示；

S1.1：数据准备：

选取50例病人，病人CT涵盖头颈、胸腹、盆腔等部位；每个病人添加72个射野(间隔0.5°)，射野的等中心位置尽可能分布在不同 CT层，使得射野可以覆盖更多的CT范围；

分别计算每个射野的PB剂量和MC剂量，值得注意的是，这里剂量网络保持统一，可以采用3mm网格；

整体来看，共有3600组CT、PB剂量、MC剂量的数据，如图3 所示；

由于CT数据的采样间距、大小与剂量网格不一致，因此需要将 CT数据插值到与剂量网格一致的大小和间隔。

S1.2：数据划分：将3600组数据按7：2：1的比例划分为训练集、验证集和测试集。

S1.3：创建卷积神经网络：卷积神经网络由输入层、融合层、卷积层、激活层、池化层、上采样层、输出层组成。

对CT与PB剂量图像执行随机裁剪(crop)，从而提高神经网络的泛化能力，处理后的图像大小为128×128，两者经融合层形成一组输入数据(128×128×2)，而后经64个卷积核处理，形成128×128×64 的张量，继续卷积处理，形成128×128×64的张量，经最大池化处理后形成64×64×64的张量，以此类推，直至上采样、卷积处理后得到 128×128×64的张量，再次卷积后得到128×128的输出图像。

卷积神经网络的输出真值为MC剂量，损失函数定义为MC剂量和预测剂量之间的均方误差(Mean Square Error,MSE)，是目标变量和预测值之间的平方距离之和的均值：

其中：n为128×128。

S1.4：训练卷积神经网络；

网络训练的过程就是在不断迭代求解卷积核，是特征提取的过程；

训练的具体步骤为：

i.随机选取数据(batch size＝512)，对输入图像进行在线数据增广，方法包含平移、2D旋转，目的是增大数据量，提高卷积神经网络的泛化能力；

ii.使用随机梯度下降法反向传播得到网络所有的权重和偏移的梯度，使用学习率10^-5更新权重和偏移，在此基础上，得到新的预测值，与MC剂量对比，得到均方误差，重复此过程200次，完成一个epoch；

iii.重复步骤i),ii)，直至均方误差小于阈值或达到最大epoch(如： 5000)，迭代停止，卷积神经网络训练完毕；

上述计算过程在NVIDIARTX 2080Ti上执行。

在该实施例中，将数据划分为训练集、验证集和测试集。在训练的过程中，一般在训练时，几个epoch结束后跑一次验证集看下效果。从而可以及时发现模型或者参数的问题，及时终止训练，重新调参或者调整模型。此外，可以验证模型的泛化能力，如果在验证集上的效果比训练集上差很多，就应考虑模型是否过拟合了。

S1.5：测试效果；

在训练好的卷积神经网络上测试，输入测试集的CT和PB剂量数据，得到预测剂量，将所有的预测剂量数据保存在本地，使用剂量验证软件(如：DVS,LinaTech)计算预测剂量和MC剂量的gamma 通过率(3mm，3％，10％)，通过率在95％以上即为合格。

gamma通过率是衡量“测量”和“计算”间差异的定量方法，是吻合距离(DistanceTo Agreement,DTA)和剂量差异(Dose Difference,DD) 的组合：

r(r_m,r_c)＝|r_c-r_m|；

δ(r_m,r_c)＝D_c(r_c)-D_m(r_m)；

其中：

r(r_m,r_c)代表对于测量点r_m与计算点r_c间的距离；

δ(r_m,r_c)代表两者之间的剂量差；

Δd_M代表距离阈值(可以但不限于取3mm)；

ΔD_M代表剂量的阈值(可以但不限于取最大剂量的3％)；

当γ(r_m)≤1时，该点通过，否则不过。

若只考虑距离因素(δ(r_m,r_c)＝0)，那么测量点通过的依据是：在距离阈值(如：3mm)内，存在一个计算点的剂量值等于测量点；

若只考虑剂量因素(r(r_m,r_c)＝0)，那么测量点通过的依据是：在同一位置，计算剂量和测量剂量之间的差不大于剂量阈值(如：最大值的3％)。

如图7所示，其显示某一层PB剂量和neural剂量(基于神经网络的剂量，即为卷积神经网络的预测剂量)的gamma通过率，MC剂量为真值，可见neural剂量更接近于MC剂量，统计结果显示，neural 剂量的gamma通过率达到99％，而PB剂量的gamma通过率仅有62％。

为了进一步地优化本发明的实施效果，在另外一些实施方式中，其余特征技术相同，不同之处在于，获取基于MC算法的沉积矩阵，其原理为：

其中，A为基于MC算法的沉积矩阵；

D_PB与D_MC分别为基于PB算法的剂量与基于MC算法的剂量；

a_jm为基于MC算法的沉积矩阵的每一个值；

代表该体素点(第j个)的MC剂量相对于第k个PB 剂量的导数；

dD_PB,k/dx_m代表第k个PB剂量相对于第m个强度值的导数；

N_dose代表剂量个数；

由上式可得，通过对PB剂量到MC剂量映射的卷积神经网络的梯度求解，能够得到基于MC算法的沉积矩阵。

值得注意的是，上述公式将问题“使用MC算法计算沉积矩阵”转化为“求解MC剂量对PB剂量的导数”，其中：PB剂量对强度值x的导数 dD_PB,k/dx_m是原算法得到的沉积矩阵，即TPS软件里使用主射线和散射线叠加的方法计算出的矩阵。

而“求解MC剂量对PB剂量的导数”的本质就是“PB剂量到MC 剂量映射的卷积神经网络”的梯度求解，求解方法可参考反向传播算法。

1)求解各个类型网络的梯度，下面分别描述卷积层(带激活层)、池化层、上采样层的梯度求解。

i.带激活层的卷积层梯度：如图8所示，记输入为X，卷积层的权重和偏置分别为W和b，输出为Y，激活函数为那么有：

conv2(W,X)代表W与X的离散二维卷积，也可以写作：

因此有

其中：

k反映了卷积核大小(2*k+1＝kernel size)，当采用3×3卷积核时，k为1；

激活函数若为ReLU函数，那么它的导数为：

激活函数若为Sigmoid函数，那么它的导数为：

依照上述公式可以得到卷积层(带激活函数)的导数。

ii.池化层的梯度：以最大池化(Max Pooling)，2×2的卷积核为例，该层的输出等于对应输入的最大值，即 y(i,j)＝max{x(p,q),p＝i-1,i,q＝j-1,j}，那么梯度可表示为：

iii.上采样层的梯度：上采样(Up sampling)有三种方法：插值法、反卷积和反池化；

a)如图10所示，插值法：将输入图像放大指某一尺寸，使用插值法计算每一点的像素值，以双线性插值为例，上采样后的图像某一点像素值表示为下式，其中上采样图像i,j处对应了原图的I，J:，输出对输入的梯度记为插值权重w；

b)如图11所示，反卷积：卷积可以表示为两个矩阵的乘积，反卷积就是卷积核矩阵的转置与原矩阵相乘，是小尺寸到大尺寸的一种转换方法，本质也是一种卷积，它的梯度计算方法可以参考卷积的方法。

c)如图9所示，反池化：以最大池化为例，记录下max-pooling 在对应kernel中的坐标，在反池化过程中，将一个元素根据kernel进行放大，根据之前的坐标将元素填写进去，其它位置补0。

它的梯度计算类似于池化层的方法：最大值处的梯度直接传递，其他位置梯度为0。

2)依次反向计算梯度：考虑包含一个输入层A、隐藏层B、输出层C的网络，信息传递方向为A→B→C，C中每一个点C_ij相对于 A中每一个点A_st的梯度可通过B传递，写作:

考虑整个卷积神经网络，从最末尾的输出层出发，依次计算出相对于前一层、前两层、前三层的梯度，不断向前递归，直至计算得到相对于输入层的梯度。

为了进一步地优化本发明的实施效果，在另外一些实施方式中，其余特征技术相同，不同之处在于，FMO具体包括以下步骤：

S3.1：将射野分割成笔形束；

S3.2：引入基于MC算法的沉积矩阵；

S3.3：设置初始强度值，将其作为优化的初始值；

S3.4：利用基于MC算法的沉积矩阵求解目标函数的导数；

S3.5：求解当前梯度下的最优步长；

S3.6：更新强度值、剂量和目标函数；

S3.7：重复步骤S3.4-S3.6，直至满足停止条件后，迭代停止。

在具体的实施例中，强度分布优化是计划优化两步法的第一步，使用梯度下降法执行此优化的步骤见图4，包括：

S3.1：射野分割成笔形束(beamlets)：将设置的射野钨门范围划分成N×M个矩形网络，N代表光栅的宽度方向，M代表光栅运动方向，通常情况下N方向的网格按照光栅宽度划分，便于执行下一步的光栅子野分割，M方向按照等间距分割；

S3.2：引入基于MC算法的沉积矩阵A；

S3.3：设置初始强度值：作为优化的初始值；

S3.4：求解目标函数的导数：这里定义的优化目标函数为：

其中：

是目标函数，表示为靶区(PTV)与危及器官(OAR)的剂量二次型误差的加权；

是危机器官的剂量误差值，α是它的加权；

是危机器官的剂量误差值，β是它的加权；

是待优化变量,即为所有射野离散成的笔形束的强度通量图；为OAR的二次型误差，是所有器官中所有体素点计算剂量和理想剂量之差的二范数；

N_OAR为危及器官的个数；

N_i为第i个危机器官的体素点个数；

w_j为第j个体素点的权重；

δ_j为惩罚因子，当计算剂量满足预设的约束时，即为0，否则记为1；

为根据射野强度通量图计算的j点剂量值，是计算剂量；

p_j为理想剂量值；

同理，为靶区的二次型误差，其中N_PTV为靶区的体素点个数；

a_jm为第m个笔形束对第j个体素的剂量贡献，组成的矩阵就是基于MC算法的沉积矩阵A。

因此，目标函数的导数为：

S3.5：求解最优步长：求解当前梯度下的最优步长；

S3.6：更新强度值、剂量和目标函数；

S3.7：判断是否满足停止条件：当自变量(强度值x)或目标函数J_obj(x)的绝对变化值或相对变化值小于阈值ε时，迭代停止，否则重复步骤S3.4-S3.4。

为了进一步地优化本发明的实施效果，在另外一些实施方式中，其余特征技术相同，不同之处在于，该放疗计划强度分布优化方法还包括以下步骤：

S4：根据放疗计划强度分布优化方法计算人体内的剂量和各个器官的DVH，评估计划效果。

在具体的实施例中，直接使用gDPM计算所有笔形束(beamlets) 的沉积矩阵(总粒子数约为2×10⁹)，将其代入优化流程，得到最优强度分布，进一步地获得剂量分布，作为强度分布优化的基准(baseline)。

使用本发明公开的方法得到基于MC算法的沉积矩阵，加入优化流程中去，得到优化强度分布，使用MC算法计算人体内的剂量，并统计出器官的DVH。

对比两者的DVH，验证本发明的计划效果。

两者的结果图可见图12，图12中：neural dose代表基于神经网络的结果(即本发明的结果)，MC剂量代表基于随机计算的结果，可见两者结果几乎完全一致，验证了本发明的可行性。

另一方面，本发明实施例还公开一种放疗计划强度分布优化装置，包括：

神经网络映射模块，用于通过神经网络，将PB剂量映射到MC 剂量上；

沉积矩阵获取模块，用于求解MC剂量相对于PB剂量的梯度，得到基于MC算法的沉积矩阵；

强度分布优化模块，用于将基于MC算法的沉积矩阵带入FMO 中，实现放疗计划强度分布优化。

另一方面，本发明实施例还公开放射治疗设备，用于执行上述任一实施例公开的放疗计划强度分布优化方法，或，包括：上述实施例公开的放疗计划强度分布优化装置。

本发明一种放疗计划强度分布优化方法、装置及放射治疗设备具有以下有益效果：

本发明将MC算法应用在FMO中可以确保优化过程得到的强度分布是真实可靠的，使用这样的强度分布进行多叶光栅子野分割，得到的子野走位是更加合理的，保证了最终得到的计划剂量分布更接近于临床要求，计划效果有显著提升。

上述实施例只为说明本发明的技术构思及特点，其目的在于让本领域普通技术人员能够了解本发明的内容并加以实施，并不能以此限制本发明的保护范围，凡根据本发明精神实质所作的等效变化或修饰，都应涵盖在本发明的保护范围内。