空间设定差异对蒙特卡罗方法计算比吸收分数精度和效率的影响

引用本文

张卫媛, 卓维海, 陈波, 季华钧. 空间设定差异对蒙特卡罗方法计算比吸收分数精度和效率的影响[J]. 中国辐射卫生, 2022, 31(5): 597-600, 605. DOI: 10.13491/j.issn.1004-714X.2022.05.014.

ZHANG Weiyuan, ZHUO Weihai, CHEN Bo, JI Huajun. Effects of grid methods on precision and efficiency of specific absorbed fraction calculation in Monte Carlo simulation[J]. Chinese Journal of Radiological Health, 2022, 31(5): 597-600, 605. DOI: 10.13491/j.issn.1004-714X.2022.05.014.

基金项目

上海市卫生和计划生育委员会青年项目（2017Y0140）

通讯作者

卓维海，E-mail：whzhuo@ fudan.edu.cn

文章历史

收稿日期：2022-03-26

Contents Abstract Full text Figures/Tables PDF

空间设定差异对蒙特卡罗方法计算比吸收分数精度和效率的影响

张卫媛 , 卓维海 , 陈波 , 季华钧

复旦大学放射医学研究所，上海　200032

收稿日期：2022-03-26

基金项目：上海市卫生和计划生育委员会青年项目（2017Y0140）

作者简介：张卫媛（1986—），女，河南周口人，实验师，从事放射卫生与辐射防护工作。E-mail：wyzhang@ fudan.edu.cn.

通讯作者：卓维海，E-mail：whzhuo@ fudan.edu.cn.

摘要：目的比较2类栅元法计算参考人比吸收分数的计算结果精确度和计算时间代价。方法基于国际放射防护委员会提供的体素化成年女性参考人模体，以肝脏为源器官发射0.5 MeV单能光子为例，利用MCNP6.0^®的F4、F6和*F8卡，采用网格法和重复结构晶格法计算各个靶器官/组织的比吸收分数（Specific Absorbed Fraction, SAF），以计算结果和计算时间为评估参数，研究27个器官/组织的SAF值。结果与已知计算结果比较，除眼晶状体和皮肤外，当前研究结果得到25个器官/组织SAF值的相对偏差绝对值均小于5%。基于重复结构晶格的*F8卡计算结果各器官相对偏差均最小，但计算时间也最长；基于网格法F4卡计算时间比基于重复结构晶格的F6卡的计算时间稍长，F6晶格法的计算时间最短。各卡计算结果的标准误差在±5%内。结论在计算结果上，同时模拟计算初级和次级粒子输运的*F8卡的精度最高；在同类计数卡中，使用网格法所需的计算时间较使用晶格法长。

关键词：比吸收分数蒙特卡洛方法网格法重复结构晶格法

Effects of grid methods on precision and efficiency of specific absorbed fraction calculation in Monte Carlo simulation

ZHANG Weiyuan , ZHUO Weihai , CHEN Bo , JI Huajun

Institute of Radiation Medicine, Fudan University, Shanghai 200032 China

Corresponding author: ZHUO Weihai, E-mail: whzhuo@fudan.edu.cn.

Abstract: Objective To compare the precision and efficiency of computing the specific absorbed fraction (SAF) of a reference human with two grid methods in MCNP6.0. Methods Based on the adult female reference voxel phantom provided by the International Commission on Radiological Protection, assuming the liver as the source organ emitting single-energy photons (0.5 MeV), the SAF of each target organ/tissue was calculated by using the mesh method and repeated structure lattice method with the F4, F6, and *F8 tally cards in MCNP6.0. We compared the methods by assessing the relative deviation of SAF and computing time for 27 organs/tissues. Results Compared with reported data, the absolute values of relative deviations of SAF values for all the organs/tissues were less than 5%, except for the eye lens and skin. By using the repeated structure lattice-based *F8 tally, the relative deviations of SAF values of the organs/tissues were all smallest, but with the longest computing time. The computing time of the mesh-based F4 tally was slightly longer than that of the repeated structure lattice-based F6 tally, which was shortest. Conclusion The *F8 tally simultaneously simulating primary and secondary particle transport showed the highest precision. The mesh tally requireda longer computing time than the lattice tally when using the same tally card.

Key words: Specific absorbed fraction Monte Carlo method Mesh tally Repeated structure lattice tally

美国核医学会内照射剂量委员会（Medical Internal Radiation Dose Committee of Society of Nuclear Medicine, MIRD）推荐的MIRD模型及方法^[1]，至今仍是内照射剂量估算的常用方法。该方法将对于某一单能光子照射下的一对源器官和靶器官，沉积在靶器官的能量占源器官发射光子的能量比例定义为吸收分数（absorbed fraction, AF），而单位质量靶器官的吸收分数（specific absorbed fraction, SAF）定义为比吸收分数，SAF可方便地将源活度积分等物理量直接换算出内照射剂量。除实测外，多利用蒙特卡罗（MC）方法，基于MCNP^®、Geant^®、Fluka^®、PHITS^®等标准代码计算获得SAF系数（库），并针对不同性别、年龄和种族的参考人模型开展了大量的SAF研究工作^[2-7]。MC方法因其固有特点，在计算源距大、器官体积小等参数时往往存在计算时间长，结果精度低等缺点，研究者多采用并行计算和减方差计算技术克服上述问题^[8-9]。

本文从计算对象空间分割方法入手研究并获得了若干解决方法。具体是基于MCNP^®软件，利用不同剂量计数卡中的网格栅元法和重复结构晶格栅元法计算SAF并评估系列计算结果和所需时间并给出优化方案。

1 材料与方法 1.1 体素模体

本研究采用了ICRP 110^[10]推荐的基于“Laura”^[11]的体素化成年女性参考人计算模体，模体源于真实的医学图像数据，符合ICRP 89^[12]号报告中给出的女性成年人解剖生理学参数，主要参数见表1。该模体内含141个器官/组织和每一个器官/组织的元素组成和密度，单体素尺寸为1.775 mm×1.775 mm×4.84 mm。

表 1 成年女性参考人模体的主要参数 Table 1 Main parameters of adult female reference phantom

1.2 计算方法与模型参数 1.2.1 计算方法

本研究使用美国洛斯阿拉莫斯国家实验室开发的MCNP6.0^[13]软件开展MC模拟计算。利用MCNP6.0开展粒子输运模拟计算辐射剂量时，可选用的剂量卡有F4、F6和*F8卡；对计算问题的空间划分，可选择网格栅元方法和重复结构晶格栅元方法。网格栅元有矩形网格、圆柱形网格和球形网格共3种，这3种方法均可通过TMESH和FMESH来实现，网格粗细可自由定义，对问题几何的结构没有特殊要求。重复结构晶格栅元法只能用于有重复结构的问题几何，当使用填充卡并将重复结构填充到单个栅元内时，MCNP6.0中可直接用*F8和F6+SD卡得到各体素的能量沉积和吸收剂量。

1.2.2 模型参数设置

本研究假设肝脏为单一源器官，发射能量为0.5 MeV的单能光子。在MCNP6.0中采用光子和电子联合输运模式计算，采用默认的光子和电子的能量截止值；使用F6和*F8晶格栅元法以及F4和F6网格栅元法计算上述女性模体内的器官/组织的比吸收分数。为实现结果比较，网格栅元的大小和坐标位置与体模体素完全一致，即网格栅元对应于体模的每一体素。

1.2.3 计算条件和数据处理

模拟计算在40核80线程的双路服务器上进行，设定光子数为1.0×10⁷个，计算结果标准偏差在±5%内。通过Python^®自编程，经过坐标变换将计算结果归一化到器官/组织水平的比吸收分数。在转换过程中为简化数据的集成，以特定的逗号分隔格式和特定的命名约定记录SAF和相应的误差。

2 结　果 2.1 SAF计算结果比较

表2汇总了本研究基于不同剂量计数卡和不同栅元方法计算得到的参考人模体内27个器官/组织的比吸收分数（SAF）值，其中最右一列给出Maxime等^[14]研究结果用于可信度比对，其源项条件、模体参数设置与本研究相同。

表 2 27个器官/组织的SAF计算结果比较（kg⁻¹） Table 2 Comparison of specific absorbed fractions of 27 organs/tissues (kg⁻¹)

器官	F4网格法	F6网格法	F6晶格法	*F8晶格法	Maxime等
肌肉	1.019 × 10⁻²	1.014 × 10⁻²	1.014 × 10⁻²	1.016 × 10⁻²	1.019 × 10⁻²
皮肤	9.919 × 10⁻³	9.869 × 10⁻³	9.869 × 10⁻³	9.232 × 10⁻³	9.133 × 10⁻³
肝脏	1.093 × 10⁻¹	1.089 × 10⁻¹	1.089 × 10⁻¹	1.065 × 10⁻¹	1.063 × 10⁻¹
大脑	2.685 × 10⁻⁴	2.661 × 10⁻⁴	2.661 × 10⁻⁴	2.703 × 10⁻⁴	2.739 × 10⁻⁴
肺	5.898 × 10⁻²	5.887 × 10⁻²	5.887 × 10⁻²	5.900 × 10⁻²	5.890 × 10⁻²
红骨髓	8.578 × 10⁻²	8.513 × 10⁻²	8.513 × 10⁻²	8.580 × 10⁻²	8.580 × 10⁻²
小肠壁	9.409 × 10⁻³	9.363 × 10⁻³	9.363 × 10⁻³	9.317 × 10⁻³	9.314 × 10⁻³
乳腺	3.692 × 10⁻²	3.673 × 10⁻²	3.673 × 10⁻²	3.705 × 10⁻²	3.694 × 10⁻²
骨表面	1.037 × 10⁻¹	1.029 × 10⁻¹	1.029 × 10⁻¹	1.039 × 10⁻¹	1.040 × 10⁻¹
结肠	2.432 × 10⁻²	2.417 × 10⁻²	2.417 × 10⁻²	2.429 × 10⁻²	2.423 × 10⁻²
肾脏	1.616 × 10⁻¹	1.608 × 10⁻¹	1.608 × 10⁻¹	1.601 × 10⁻¹	1.611 × 10⁻¹
心脏	1.805 × 10⁻²	1.789 × 10⁻²	1.789 × 10⁻²	1.789 × 10⁻²	1.792 × 10⁻²
胃	2.857 × 10⁻²	2.857 × 10⁻²	2.857 × 10⁻²	2.772 × 10⁻²	2.772 × 10⁻²
脾脏	1.229 × 10⁻²	1.223 × 10⁻²	1.223 × 10⁻²	1.235 × 10⁻²	1.231 × 10⁻²
胰腺	3.450 × 10⁻²	3.434 × 10⁻²	3.434 × 10⁻²	3.450 × 10⁻²	3.430 × 10⁻²
子宫	7.112 × 10⁻⁴	7.039 × 10⁻⁴	7.039 × 10⁻⁴	7.257 × 10⁻⁴	7.166 × 10⁻⁴
淋巴结	9.148 × 10⁻³	9.092 × 10⁻³	9.092 × 10⁻³	9.104 × 10⁻³	9.114 × 10⁻³
唾液腺	1.669 × 10⁻³	1.655 × 10⁻³	1.655 × 10⁻³	1.714 × 10⁻³	1.689 × 10⁻³
膀胱	7.240 × 10⁻⁴	7.170 × 10⁻⁴	7.170 × 10⁻⁴	7.198 × 10⁻⁴	7.209 × 10⁻⁴
食管	1.520 × 10⁻²	1.515 × 10⁻²	1.515 × 10⁻²	1.512 × 10⁻²	1.506 × 10⁻²
口腔黏膜	2.042 × 10⁻³	2.031 × 10⁻³	2.031 × 10⁻³	2.037 × 10⁻³	2.084 × 10⁻³
胸腺	4.990 × 10⁻³	4.961 × 10⁻³	4.961 × 10⁻³	5.044 × 10⁻³	5.004 × 10⁻³
甲状腺	2.419 × 10⁻³	2.407 × 10⁻³	2.407 × 10⁻³	2.558 × 10⁻³	2.408 × 10⁻³
肾上腺	1.105 × 10⁻¹	1.101 × 10⁻¹	1.101 × 10⁻¹	1.066 × 10⁻¹	1.072 × 10⁻¹
卵巢	1.465 × 10⁻³	1.458 × 10⁻³	1.458 × 10⁻³	1.485 × 10⁻³	1.476 × 10⁻³
胆囊	8.549 × 10⁻²	8.513 × 10⁻²	8.513 × 10⁻²	8.271 × 10⁻²	8.326 × 10⁻²
眼晶体	1.455 × 10⁻³	1.444 × 10⁻³	1.444 × 10⁻³	1.225 × 10⁻³	1.269 × 10⁻³

表 2 27个器官/组织的SAF计算结果比较（kg⁻¹） Table 2 Comparison of specific absorbed fractions of 27 organs/tissues (kg⁻¹)

为进一步分析不同计算方法获得结果的差异，表3汇总了SAF值相对偏差值大于1%的器官或组织相关参数。

表 3 SAF值相对偏差大于1%的器官/组织 Table 3 Organs/tissues with specific absorbed fraction’s relative deviations greater than 1%

2.2 计算时间比较

在控制计算标准偏差小于±5%的前提下，表4汇总了本研究使用不同栅格方法的CPU计算时间。

表 4 不同剂量计数卡2类栅元方法的计算时间比较 Table 4 Comparison of computing time with different grid methods and tallies

3 讨　论

从表2可知，当肝脏为单一源器官并发射0.5 MeV光子时，SAF值较大的主要器官/组织有肾脏、肾上腺、肝脏、骨表面和红骨髓等，其SAF值均接近或超过每千克0.1。

从表3可知，部分器官/组织SAF值相对偏差大于1%主要由于质量较小、位置远离肝脏以及形状欠规则所致；其中偏差较大的皮肤和眼晶状体是因皮肤全身分布，对应栅元空间大且薄，眼晶状体体积和质量很小且距肝脏（源器官）较远，在同样抽样次数前提下，计算误差相比更大，此类远距离、小体积器官/组织的SAF值精度提升十分困难。同时考虑初级粒子（入射光子）和次级电子行为的*F8卡计算结果相对偏差最小，适用于当前需求。

从表4可知，网格法*F8卡计算时间最长，网格法F6卡较短，但两者远大于网格法F4卡和晶格法F6卡。这是因为*F8卡增加计算了次级电子的输运，晶格法F6卡可直接记录各体素内单位体积的初级光子沉积能量，所需计算时间最短。网格法F6卡记录栅元并不假定与模体体素一致^[13]，需要对粒子位置进行对偶计算，需要追踪每个历史步骤中粒子在体素单元中的位置（即使网格栅元与体素的大小和位置是否完全一致），重复跟踪会使网格计数变慢。F4卡只计算栅元中入射粒子的平均通量，结合通量和能量转换系数表可直接快速给出栅元内的能量沉积，但结果准确性依赖于系数转换表。F6卡晶格法计算速度快和定义简单，它还为剂量计算结果和其他统计信息的收敛提供了统计测试^[15]进而定性增加了剂量计算的置信度，晶格计数的更高计算效率有助于常规剂量模拟。

上述研究初步表明，使用不同剂量卡和栅格法均可获得较为准确的SAF值，同时考虑初级和次级粒子输运计算的*F8卡计算精度最高，计算时间会增加很多。同等条件下采用F6卡，网格栅元法所需计算时间远多于晶格栅元方法。对于体素不够精细的体素化模体，晶格栅元法对小器官SAF值的计算精度会受较大影响，建议在开展类似计算前，要综合考虑体素的精度和所要计算器官或组织的质量，选取合适的网格栅元方法。

在实际工作中，计算SAF值使用的源项不尽相同，受研究时间限制，本研究尚未系统完成。今后拟开展多能量下SAF计算结果的对比，以期获得更系统的研究结果。

参考文献

[1]	Bolch WE, Eckerman KF, Sgouros G, et al. MIRD Pamphlet No. 21: a generalized schema for radiopharmaceutical dosimetry-standardization of nomenclature[J]. J Nucl Med, 2009, 50(3): 477-484. DOI:10.2967/jnumed.108.056036
[2]	Abrefah RG, Essel PAA, Odoi HC. Estimation of the dose rate of nuclear fuel of Ghana Research Reactor-1 (GHARR-1) using ORIGEN-S and MCNP 6[J]. Progress Nucl Energy, 2018, 105: 309-317. DOI:10.1016/j.pnucene.2018.02.002
[3]	梅雪松, 张慧, 马玉刚, 等. 用Monte Carlo模拟法和理论计算法研究¹²⁵I粒子源的剂量分布[J]. 中国辐射卫生, 2011, 20(2): 151-153. DOI: 10.13491/j.cnki.issn.1004-714X.2011.02.013. Mei XS, Zhang H, Ma YG, et al. Study of dose distribution of ¹²⁵I seed using Monte Carlo simulation and the theoretical calculation[J]. Chin J Radiol Health, 2011, 20(2): 151-153.
[4]	Abo Shdeed T, Nahili M, Abo Issa N, et al. Study of absorbed dose in important organs during helical CT chest scan using MCNP code and MIRD phantom[J]. Egyp J Radiol Nucl Med, 2016, 47(3): 1649-1633. DOI:10.1016/j.ejrnm.2016.09.004
[5]	Leal-Acevedo B, Gamboa-deBuen I. Dose distribution calculation with MCNP code in a research irradiator[J]. Radiat Phys Chem, 2020, 167: 108320. DOI:10.1016/j.radphyschem.2019.05.010
[6]	洪天祺, 戴瑜. MCNP5计算¹⁸F药物诊断的吸收剂量率 [J]. 中国辐射卫生, 2015, 24(2): 125-126. Hong TQ, Dai Y. Absorbed dose ratecalculation of ¹⁸F drug diagnosis with MCNP5 [J]. Chin J Radiol Health, 2015, 24(2): 125-126. DOI:10.13491/j.cnki.issn.1004-714X.2015.02.010
[7]	Zubair M, Ahmed E, Hartanto D. Comparison of different glass materials to protect the operators from gamma-rays in the PET using MCNP code[J]. Radiat Phys Chem, 2022, 190: 109818. DOI:10.1016/j.radphyschem.2021.109818
[8]	Van Wijk A,Van den Eynde G, Hoogenboom JE. An easy to implement global variance reduction procedure for MCNP[J]. Ann Nucl Energy, 2011, 38(11): 2496-2503. DOI:10.1016/j.anucene.2011.07.037
[9]	Sweezy J, Brown F, Booth T, et al. Automated variance reduction for MCNP using deterministic methods[J]. Radiat Prot Dosimetry, 2005, 116(1-4): 508-512. DOI:10.1093/rpd/nci257
[10]	Zankl M. Adult male and female reference computational phantoms (ICRP Publication 110)[J]. Jpn J Heal Phys, 2010, 45(4): 357-369. DOI:10.5453/jhps.45.357
[11]	Zankl M, Becker J, Fill U, et al. GSF male and female adult voxel models representing ICRP Reference Man–the present status[C]. The Monte Carlo Method: Versatility Unbounded in a Dynamic Computing World, 2005, 1721(1).
[12]	Valentin J. Basic anatomical and physiological data for use in radiological protection: reference values: ICRP Publication 89[J]. Ann ICRP, 2002, 32(3/4): 1-277. DOI:10.1016/S0146-6453(03)00002-2
[13]	Goorley T, James M, Booth T, et al. Features of MCNP6[J]. Ann Nucl Energy, 2016, 87: 772-783. DOI:10.1016/j.anucene.2015.02.020
[14]	Chauvin Maxime, Borys D, Botta F, et al. Open Dose: Open-Access Resource for Nuclear Medicine Dosimetry[J]. J Nucl Med, 2020, 61(10): 1514-1519. DOI:10.2967/jnumed.119.240366
[15]	Zhang Rui, Fontenot JD, Mirkovic D, et al. Advantages of MCNPX-based lattice tally over mesh tally in high-speed Monte Carlo dose reconstruction for proton radiotherapy[J]. Nucl Technol, 2013, 183(1): 101-106.