程凡杰，李开健，张 巍,*，杨历军，李 岩
(1.中国原子能科学研究院，北京 102413；
2.北京中科天马信息技术有限公司，北京 100089)
硼中子俘获治疗(BNCT)是一种二元治疗癌症的方法，经过几十年的研究，被认为是目前治疗癌症的有效方法之一。其原理是将含硼药物注入人体后，含硼药物会在肿瘤区域富集，然后利用中子束照射肿瘤部位，肿瘤中的硼原子与中子反应后产生高传能线密度粒子α和7Li，这两种粒子在细胞中的射程分别为9 μm和5 μm，粒子产生的电离辐射损伤作用仅限于含10B的细胞，并且可以引起肿瘤细胞DNA不可修复的损伤，从而达到选择性地杀死癌细胞的目的。
BNCT治疗一般采用超热中子进行治疗，超热中子束穿透能力强，可以治疗深部肿瘤，是各国科学家研究的重点[1]。IAEA-TECDOC-1223报告[2]给出了BNCT装置出口束流参数，常被作为束流整形装置(BSA)优化的参考标准。
早期BNCT使用反应堆做中子源，反应堆中子源具有中子注量高、稳定性好的特点，但由于造价高和核安全要求严等问题，难以在医院普及使用。相比于反应堆中子源，加速器中子源具有安装简单、维护费用低、相对安全的优点，更容易在医院中普及。近年来，世界各国都在致力于研制基于加速器的硼中子俘获治疗(AB-BNCT)设备。
加速器打靶产生的中子能量较高，而且含有大量γ无法直接用于治疗，必须经过特定的系统进行处理，才能满足BNCT要求，这一系统即为BSA。BSA可能的设计方案多，设计变量主要有材料种类、材料顺序以及几何尺寸。
传统的BSA设计方法需要在设计之前比较不同种类的材料在慢化、准直、反射以及屏蔽γ射线等方面的能力[3-9]，以找到不同区域合适的材料，选出材料后还需要进行材料的组合和厚度的调整。但是由于材料的散射、吸收等各种反应截面随能量是变化的，当中子穿过材料时散射、吸收同时都在发生，不同尺寸材料效果也不同。这导致优化BSA时各参数会相互影响，如在降低快中子成分时，往往需要快中子截面大的材料，结果也会导致超热中子注量率下降。因此，传统方法在不同区域推荐的材料也只能是定性的结论，而且这些结论往往只对某一能量区间的中子有效，对于优化设计作用较小。由于材料种类较多、尺寸可调整的范围巨大，同时调整尺寸和材料得到的方案几乎是无限的，这个过程将花费大量的时间。得到初步方案后，再通过理论计算验证设计方案的可行性，如果不满足要求，则针对性地提出优化方案。
传统优化方法本质上是从少数个体迭代寻求最优解，这容易陷入局部最优解，难以兼顾多个参数同时优化。利用传统方法进行BSA设计工作量大、设计周期长，优化效果并不好。
为解决AB-BNCT BSA优化设计问题，将BSA设计转化为确定各栅元的材料种类、顺序以及几何尺寸的问题，优化目标是达到IAEA推荐的4个参数的要求。BSA设计转化为多个变量、多个目标参数的优化问题。像这种问题可通过智能算法来解决，遗传算法就是其中的一种。
遗传算法是一种模仿生物进化机制发展起来的随机全局搜索和优化方法，它借鉴了遗传、突变、自然选择以及杂交等进化生物学中的一些思想，是解决搜索问题的一种通用算法。该算法对于被求解问题要求低，只需要合适的适应度函数。该算法从初始解的1个集合开始搜索，同时处理群体中的多个个体，按照多种路线进行平行搜索，并且采用轮盘赌的方法进行选择。遗传算法的优势是搜索范围大，有利于全局择优。遗传算法还具有鲁棒性、可并行、可移植性的特点，遗传算法在核电厂燃料管理[10-12]、屏蔽设计[13-15]以及反应堆BNCT BSA设计[16]等方面均有应用。
本文基于中国原子能科学研究院14 MeV质子回旋加速器生成中子的特性，研制遗传算法与蒙特卡罗程序相耦合的设计程序，对BSA的材料种类、材料顺序、几何尺寸进行优化设计，以得到满足AB-BNCT治疗要求的超热中子束。
BSA设计的本质是一个带有约束条件的多目标寻优过程。一种优化方法是将慢化区进行精细划分为若干个区，然后用多种材料进填充。这种方法从原理上可行，但是会导致可能的解过多，计算量巨大，难以得到最优解。另外一种方法是先将慢化区初步划分为较少数量的区域，通过优化算法选出最可能的材料；
然后使用可能的材料进行材料顺序的优化；
最后调整每个栅元的尺寸。这种BSA设计方法转化为材料种类、材料排列顺序、材料尺寸的优化，利用遗传算法可较快得到最优解，提高设计效率。
1.1 程序架构
遗传算法主要通过编码和解码、交叉、变异、适应度评估来进行寻优，根据遗传算法原理编制了遗传算法与蒙特卡罗方法耦合的设计程序，程序流程图如图1所示。该程序主要包括4个模块：遗传算法模块、适应度计算模块、输入文件生成模块和输出文件处理模块。
图1 程序流程图Fig.1 Program flow chart
遗传算法模块主要功能是设置参数，参数包括变量数目，变量上限、下限，种群规模，交叉概率，变异概率，初始种群，保留到下一代的精英数量，终止条件等。返回的是最优个体的基因及其适应度。
适应度计算模块主要功能是调用输入文件生成模块，执行蒙特卡罗程序以及调用输出文件处理模块，得到适应度。该模块在产生输入文件之前判断是否存在该文件输入文件、输出文件，如果存在就跳过输入文件生成和蒙特卡罗程序计算两个步骤，直接处理输出文件得到适应度。由于蒙特卡罗程序计算耗费了绝大多数时间，这样做的优势是可以避免重复计算、节约时间。
输入文件生成模块主要功能是将个体基因带入到模板文件中产生蒙特卡罗程序计算的输入文件。参数个数不受限制，可以通过合适的编码实现材料尺寸和材料种类同时优化，但这样会导致可能的解更多，优化时间更长。实际优化时按照材料种类、顺序和尺寸分步进行优化，每一步变量数目较少，这样可以提高效率。
输出文件处理模块主要功能是处理蒙特卡罗程序结果文件，得到对应计数卡的结果，根据适应度公式计算适应度并返回给适应度计算模块，同时产生记录文件。在输出文件处理模块中，通过调整每个参数的权重可以针对某个参数进行优化。记录文件中的信息包括蒙特卡罗程序输出文件中的原始数据、权重因子和适应度等参数。
该程序中断运行后可以继续运行，运行结束后的蒙特卡罗程序输出文件和记录文件保留在文件夹中。初始种群可以根据记录文件个体信息选择加快收敛速度。
1.2 适应度函数
遗传算法中适应度函数是决定结果的最关键因素，是描述个体优劣的主要指标，遗传算法根据个体适应度大小进行优胜劣汰[16]。根据优化的问题，定义适应度函数为：
(1)
式中：fitness为目标函数；
wi为权重因子；
S为解空间；
k为优化的参数个数；
Ti为每个目标函数的目标值，计算时采用IAEA的推荐值；
x为向量，进行材料优化时，x为代表材料顺序的向量，进行尺寸优化时，x为代表材料尺寸的向量；
fi(x)代表方案一定情况下利用蒙特卡罗程序计算得到的各参数值。
IAEA的推荐值如表1所列，其中：En为超热中子能量范围；
df为快中子产生的剂量；
dr为γ产生的剂量；
J为中子流量；
φepi为超热中子注量率，可以看出各目标参数值有数量级的差异，该适应度函数可消除参数之间数量级的差异，实现多个目标参数同时优化。
表1 IAEA推荐的参数Table 1 IAEA recommended parameter
结合14 MeV回旋加速器生成中子的特性，按照装置辐射防护要求，确定整型装置内慢化材料为含硼聚乙烯，屏蔽材料为铅。然后将慢化区域分成9个区域，准直区分为2个区域，慢化区中每个区厚度约为8 cm，优化前相关部分结构如图2所示。然后，利用程序依次进行材料种类、材料顺序、几何尺寸的优化。
图2 优化前相关结构Fig.2 Related structure before optimization
2.1 材料种类优化
如前文所述，必须使用适当的慢化材料将产生中子慢化到超热范围。这些材料必须在超热能段具有低散射截面，但在较高的能量段具有高散射截面。某些具有较高吸收截面的慢化材料在中子到达患者病灶之前会吸收较多中子，发生(n，γ)反应，导致γ射线污染。但是，传统方法很难给出满足BSA全部要求的材料组合。
BSA常用的材料有铝、Fluental、氟化铝、镍、铋、氟化锂、氟化镁、硫、铅、含硼聚乙烯、铍、聚乙烯、不锈钢、石墨、氟化钙等[3-9]，共计15种。这些材料多数是反应堆或加速器BNCT慢化材料，也有屏蔽材料，这里不进行人为区分，全部作为备选材料。优化时对超热中子注量率、快中子成分、γ成分、流量注量率比值4个参数的权重wi分别进行设定。
利用程序进行优化，得到每代个体适应度的变化，如图3所示。
图3 适应度随每一代的变化Fig.3 Fitness change with each generation
根据优化结果，得到最优的10种方案均含有铝、氟化铝、石墨、铅、Fluental，其中最优方案的出口束流的参数列于表2。由表2可看出，各项参数较靶后的参数均有很大提升，但此方案的超热中子注量率较低，γ成分偏高，不能满足要求，需要继续优化。
表2 材料种类优化后结果Table 2 Result after material type optimization
此次优化得到了合适的材料，与最初材料相比可看出含氢和硼的材料均被舍弃。通过分析记录文件可看出，方案中含有这两种元素会大幅度降低超热中子注量率，在优化过程中含硼聚乙烯、聚乙烯等材料很容易被淘汰，这与实事相符。
2.2 材料顺序优化
根据初步优化选出的这些材料进行材料顺序的优化。由于超热中子注量率低，快中子和γ成分高，需要提高两参数的绝对权重。进行优化时，每一代适应度的变化如图4所示。由图4可看出，在200代左右优化停止，得到最优的材料排列顺序，BSA第2次优化结果列于表3。
图4 材料顺序优化适应度随每一代的变化Fig.4 Fitness change with each generation for material sequence optimization
表3 材料顺序优化后的结果Table 3 Result after material sequence optimization
Fluental的主要成分是27Al和19F，第2次优化后BSA主体元素排序基本上是208Pb、27Al和19F，在出口处用的主要是208Pb和12C。经过分析，208Pb、27Al、19F 3种核素共振峰所在的能量区间逐渐降低，其弹性散射截面也逐渐下降，这种排序可以使能量高于10 keV的中子都能得到慢化，大幅度降低了快中子成分，而在出口处采用208Pb可以减少γ成分。
2.3 几何尺寸优化
进行材料种类和材料顺序优化时，每个区域几何尺寸是固定的，两次优化后快中子成分仍然偏高，其他参数满足IAEA要求。然后优化各层之间的尺寸，直到满足参数要求。
优化后相关结构如图5所示。调整后的方案与之前方案相比减少了氟化铝厚度，增加了铅和Fluental的厚度，减少了出口铅的厚度。这是由于铅和Fluental在相应的能量范围慢化性能优于氟化铝，这样调整可进一步减少快中子份额。由于屏蔽γ所需要的铅比较少，减少铅厚度可以减少超热中子的损失。几何尺寸优化后的出口束流参数列于表4。
图5 优化后相关结构Fig.5 Related structure after optimization
表4 几何尺寸优化后的结果Table 4 Result after geometry optimization
另外模拟了加速器打靶产生的中子能谱和BSA出口处的中子能谱，如图6所示。优化后在0.5 eV～10 keV能量范围内，超热中子份额约为84.2%，如图7所示。对比图6可看出，利用该方法优化BSA设计可以大幅度提高超热中子份额，减小了快中子份额。
图6 加速器打靶产生的中子能谱Fig.6 Neutron spectrum behind cyclotron target
图7 BSA优化后的中子能谱Fig.7 Neutron spectrum after BSA optimization
本文将遗传算法与蒙特卡罗程序进行耦合，并利用该程序开展了AB-BNCT BSA优化设计，包括材料种类、材料顺序和几何尺寸。优化前只需要修改生成蒙特卡罗计算输入文件的模板和各变量的上下限，优化过程自动完成。计算结果表明，利用该程序优化得到的BSA模型满足IAEA的推荐参数要求，该方法为AB-BNCT BSA设计提供了一种有效的技术方案。该程序通用性强，也可应用于反应堆设计、核临界安全和屏蔽设计中。
猜你喜欢热中子蒙特卡罗中子(70～100)MeV准单能中子参考辐射场设计宇航计测技术(2019年3期)2019-10-293D打印抗中子辐照钢研究取得新进展表面工程与再制造(2019年3期)2019-09-18利用蒙特卡罗方法求解二重积分智富时代(2019年6期)2019-07-24利用蒙特卡罗方法求解二重积分智富时代(2019年6期)2019-07-24单晶硅受照热中子注量率的双箔活化法测量研究原子能科学技术(2019年4期)2019-05-13基于PLC控制的中子束窗更换维护系统开发与研究工业设计(2016年8期)2016-04-16DORT 程序进行RPV 中子注量率计算的可靠性验证核科学与工程(2016年3期)2016-01-03脉冲中子－裂变中子探测铀黄饼的MCNP模拟应用科技(2015年3期)2015-05-15探讨蒙特卡罗方法在解微分方程边值问题中的应用梧州学院学报(2015年3期)2015-02-28复合型种子源125I-103Pd剂量场分布的蒙特卡罗模拟与实验测定同位素(2014年2期)2014-04-16