学习笔记 | 中国空间站模拟剂量评估

发表于2022-07-24更新于2022-07-24字数统计8.5k阅读次数

本文仅为某次课程大作业，有诸多不严谨、遗漏之处，全文仅供参考！

2021年发射中国空间站天和核心舱，标志着中国空间站全面建设正式开启。随后神舟十二号乘组和神舟十三号乘组先后在轨工作了3个月与6个月，并且在将来将会实现空间站的长期在轨驻留和轮换(1)。在此背景下，空间站舱内外辐射环境，以及航天员在轨工作期间关键器官的辐射剂量与效应研究非常重要。

由于空间任务的特殊性，难以实地进行实验以获取剂量数据，因此主要方式是通过探测器获取轨道上辐射环境情况，随后建立空间站模型与人体模型，在传统的载人任务中，常使用BRYNTRN、HZETRN等基于Boltzmann输运方程的程序进行模拟，而Geant4等基于Monte Carlo方法的模拟工具在近些年来也备受人们关注。

Sun(2)等人在2012年建立了Visible Chinese Human Adult Female Astronaut（VCH-FA）体模，用于模拟女性航天员在神舟飞船类的剂量。石苗(3)根据中国成年男性数字化人体模型，结合航天员的实际情况，建立了航天员关键器官的数学参数化模型，借助Mulassis一维多层辐射屏蔽仿真程序和Geant4模拟程序计算了在各类情况下的辐射粒子在航天员关键器官内的能量沉积、吸收剂量等参数。

欧洲航空局（European Space Agency，ESA）的T. Ersmark等人曾进行了Dose Estimation by Simulation of the ISS Radiation Environment（DESIRE）项目用于模拟国际空间站（ISS）哥伦布舱，（Columbus module）的辐射粒子通量和剂量(4-7)。

该项目可以划分为4个部分：

  1. 测试Geant4用于模拟空间任务的准确性，与NASA的HZETRN、俄罗斯的SHIELD等模拟程序以及实验数据进行比较(4)。
  2. 建立哥伦布舱的不同详细程度的模型，进行模拟比较[5]。
  3. 使用SPENVIS等工具，获取ISS轨道上外部辐射场的模型[6]。
  4. 结合人体模型和粒子能谱，使用Geant4工具计算辐射剂量[7]。

本文以T. Ersmark等人的工作作为参考，结合石苗文中提供的部分中国航天的信息，建立了空间站天和核心舱的大致模型，通过SPENVIS模拟了2022年5月空间站轨道上的辐射环境，并使用MIRD体模进行了模拟。本文将分为5个部分，第1部分为轨道和空间辐射环境情况，第2部分为Geant4模型的建立，第3部分为Geant4模拟流程，第4部分为结果分析，第5部分为总结。

空间站轨道与空间辐射环境

轨道参数

航天器的轨道是进行空间任务模拟的必要条件。用于确定轨道的各类参数被称为轨道参数（轨道根数），一般情况下为了确定一个开普勒轨道（椭圆轨道）需要6个轨道参数，根据选取的测量装置不同，有若干种方式可以定义轨道参数。

两行轨道参数（Two Line Elements，TLE）是一种常用的表示轨道的方法，最初由NASA和北美防空司令部为部分较大的地球轨道航天器生成并公布，如今众多航天机构也在采用这一格式。

TLE格式的数据主体有两行，另外还可能添加表示航天器名称的第0行。随后两行提供了若干信息。中国载人航天工程在官方网站[1]公布了中国空间站的TLE信息，查询得到2022年5月11日的轨道参数如图 1所示，从中可以获得表 1的信息。

图1 中国空间站TLE参数
图1 中国空间站TLE参数
表1 提取的轨道信息
时间倾角(deg)升交点赤经(deg)偏心率近地点辐角(deg)平近点角(deg)真近点角(deg)
22.5.11 00:0041.469684.98020.001152308.145644.415944.4629

为了便于后续的分析，我们需要进一步转化，得到常用的六轨道根数。平近点角MM与偏近点角$E满足：

M=E−e⋅sin⁡EM=E-e\cdot\sin{E}

真近点角TT与偏近点角EE满足：

tan⁡T2=1+e1−etan⁡E2\tan{\frac{T}{2}}=\sqrt{\frac{1+e}{1-e}}\tan{\frac{E}{2}}

由此可以计算得到轨道的真近点角。除此之外也可借助轨道模型，由上述数据计算得到轨道近地点、远地点，过程较为复杂，本文不做描述。最终结果如表 2所示。通过表2中参数，我们可以使用SPENVIS[2]的Coordinate generators生成空间站的预计轨道，总计30天。

表2 处理后的轨道信息
倾角(deg)升交点赤经(deg)偏心率近地点辐角(deg)近地点(km)远地点(km)真近点角(deg)
41.469684.98020.001152308.1456379.2393.844.4629

空间辐射环境

空间辐射环境主要包括地球辐射带（Van Allen radiation belt）、太阳粒子事件（Solar Particle Event，SPE）和银河宇宙射线（Galactic cosmic rays，GCR）。辐射的强度不是恒定的，受太阳活动影响很大，需要考虑太阳活动极小时或极大时等极端情况，但在本文中只考虑模拟阶段30天的情况。

银河宇宙射线

银河宇宙射线是从系外进入太阳系的高能带电粒子，其通量与太阳活动负相关，主要由质子、电子和完全电离的原子核组成，其中核成分约占98%，分别为H（≈87%）、He（≈12%）和少量较重的原子核（≈1%）(8)。

GCR一般各向同性，通量约为1-10 cm-2s-1，以高能粒子为主（可达到1011 GeV），贡献绝大部分有效剂量，一般在太阳活动极小时最大，除此之外还取决太阳磁性(8)。

GCR中还包含有所谓异常成分（anomalous component），星际间的中性粒子进入太阳系后，太阳辐射使其失去电子，离子太阳风作用和碰撞加速，其穿透地磁场的能力比其他成分更大，又可分为部分电离（singly-ionized anomalous component）和完全电离（fully-ionized anomalous component）。

GCR常使用的模型有ISO-15390模型、CRÈME[3]模型 和Nymmik模型等。ISO-15390模型不包含异常成分，而异常成分是有效剂量的主要来源；Nymmik模型与ISO-15390模型基本相同，但是在能量低于10 MeV/u时，通量随着能量降低而增加；而CRÈME模型是美国海军研究实验室负责开发的，基于Nymmik模型，但额外考虑了GCR中的异常成分。

借助CRÈME86模型，我们可以得到其在M=3情况下的GCR通量，该情况表示在任意情况下，只有10%的概率实际通量会超过这一值。在SPENVIS中使用CRÈME86模型可以生成1号~28号元素的通量，其中H和He的通量如图 2所示。

图2 GCR中H和He的通量
图2 GCR中H和He的通量

太阳辐射

太阳辐射的主要来源是日冕抛射和太阳耀斑，太阳辐射以质子为主（90~95%），通量可达10 9 cm-2s-1；其次为α粒子（5~8%），另外还有少量电子和Z>2的离子(9)。太阳活动具有明显的周期性，主要的周期有11年、22年等，当太阳活动较强时，质子的比例会上升(9)。

主要考虑太阳质子和α粒子事件，其模型主要有King模型、JPL模型、ESP模型和SAPPHIRE模型等。SAPPHIRE是一个较新的模型，参考了其他模型太阳活动周期的规定，由于SAPPHIRE模型能够提供任务期间累计太阳质子注量，因此我们使用这一模型生成任务期间总通量，实际情况只有5%的可能大于这一估计值，结果如图 3。

图3 太阳质子和α粒子通量
图3 太阳质子和α粒子通量

可以发现整个任务期间的总通量不超过10 3 cm-2，远远小于银河宇宙射线和地球辐射带的通量，因此在后续模拟中忽略这一因素的剂量。

地球辐射带

近地地球辐射带被称为Van Allen带，是部分宇宙射线、宇宙射线与地球磁场和大气层作用产生的带电粒子被地磁场俘获形成的辐射带，分为内外两个辐射带(8)。

内带的高度一般在600 km以上，但是在高纬度地区和南大西洋地区，由于地磁异常，内带高度会下降至约300 km，因此只有经过这些区域时才会有明显的俘获带电粒子通量(8)（见图 4），为方便计算，将通量平均到整个轨道上。

图4 空间站轨道上辐射带质子（左）和电子（右）的分布
图4 空间站轨道上辐射带质子（左）和电子（右）的分布

近地地球辐射带的常用模型是AE/AP系列模型，这是由NASA综合20多颗不同卫星的轨道探测数据得到的，目前最新版本是AE9/AP9模型。使用该模型计算轨道上平均通量，结果如图 5所示。

图5 AE9/AP9得到的质子（左）和电子（右）微分通量
图5 AE9/AP9得到的质子（左）和电子（右）微分通量

辐射模型的拟合

上述模型得到的结果是分立的，Geant4中可以使用General Particle Source（GPS）进行直方图的拟合，也可以自行编写函数，使用Particle Gun进行。考虑到本文所述部分模型较为复杂，同时便于获取发射粒子的能谱、位置、方向等信息，因此采用自行编写的方法。

银河宇宙射线

定义粒子刚度（particle rigidity）R=pcZeR=\frac{pc}{Ze}，其中pp为粒子动量，cc为光速，ZeZe为粒子电荷量，则有：(10)

F(p,t)=CiβαiR(p)γi[R(p)R(p)+(0.37+3×10−4W(t)1.45)]b⋅W(t)+cAiZi1βF\left(p,t\right)=\frac{C_i\beta^{\alpha_i}}{R\left(p\right)^{\gamma_i}}\left[\frac{R\left(p\right)}{R\left(p\right)+\left(0.37+3\times{10}^{-4}W\left(t\right)^{1.45}\right)}\right]^{b\cdot W\left(t\right)+c}\frac{A_i}{Z_i}\frac{1}{\beta}

其中pc=2m0c2E+E2pc=\sqrt{2m_0c^2E+E^2}，EE为粒子动能；W(t)W\left(t\right)为表征太阳活动的量，本文由于时长（30天）远小于太阳周期，因此可认为是常数；CiC_i，γi\gamma_i，αi\alpha_i为与粒子有关的常数，AiA_i，ZiZ_i分别为质量数与电荷数，β=vc=pcm02c4+E2\beta=\frac{v}{c}=\frac{pc}{\sqrt{m_0^2c^4+E^2}}，则可简化为：

F(pc)=C1βα−1pc(pcpc+C2)C3F\left(pc\right)=\frac{C_1\beta^{\alpha-1}}{pc}\left(\frac{pc}{pc+C_2}\right)^{C_3}

同时考虑双指数模型：

F(E)=C1e−C2E(1−e−C3E+C4)F\left(E\right)=C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)

两个模型同时进行拟合，随后进行贝叶斯信息准则测试（BIC），最终发现H使用简化模型效果较好，其余元素均使用双指数模型。

地球辐射带电子

定义LL为轨道高度（含地球半径）与地球半径的比值。当L>2.5L>2.5时，有64%的能谱呈现指数分布，有32%的能谱呈现bump‐on‐tail（BOT）分布，剩余4%呈现幂律分布(11)。

中国空间站轨道较低（L≈1.06L\approx1.06），不考虑较为复杂的BOT分布，使用指数分布（j=j0e−EE0j=j_0e^{-\frac{E}{E_0}}）和幂律分布（j=j0E−αj=j_0E^{-\alpha}）同时进行拟合，最终选择采用BIC值较低的指数分布模型。

地球辐射带质子

由于未搜索到低轨道质子能谱的相关资料，因此采用与电子类似的方法，考虑增加平移项将幂律分布模型修改为j=j0(E−β)−αj=j_0\left(E-\beta\right)^{-\alpha}，最终选用BIC值较低的幂律分布模型。

Geant4模型

空间站

空间站核心舱模型包括几何模型与材料填充，由于核心舱公开资料少，因此参考网络资料、ISS舱段以及石苗文中提供的信息，进行核心舱的建模。

几何尺寸

核心舱分区如图 6所示。DESIRE项目通过建立三个不同详细程度的哥伦布舱模型，研究了几何模型对模拟结果的影响(5)，结果显示舱段细节对剂量率模拟基本没有影响，尤其是边角处的部件重要性有限，因此将天和核心舱简化为5个部分，几何模型与尺寸见表 3。

图6 核心舱分区
图6 核心舱分区
表3 核心舱模型
部位几何模型尺寸(m)
大柱段圆柱形Φ4.2×L 6.72\Phi4.2\times L\ 6.72
小柱段圆柱形Φ2.8×L 5.18\Phi2.8\times L\ 5.18
节点舱球形Φ2.8\Phi2.8
过渡段1圆台形Φ2.8−2.13×L 0.815\Phi2.8-2.13\times L\ 0.815
过渡段2圆台形Φ4.2−2.8×L 0.815\Phi4.2-2.8\times L\ 0.815

材料

空间站主要可分为外壳、填充材料与中间区域3个部分，中间区域主要为空气，填充区域主要为各类金属与有机材料，外壳较为复杂。ISS舱段外壳如图 7所示，均由3层组成，内外为Al合金，中间为Kevlar等有机材料。

国内Kevlar类产品主要有Taparan（泰普龙）等，其分子式为[C14H10N2O2]n\left[C_{14}H_{10}N_2O_2\right]_n，密度为1.44 g/cm3。外壳主要为5系铝合金，在本文中给直接使用纯Al，最终外壳为3层，为2 mm铝与10 mm Taparan + 5 mm铝。模型总重为22.1吨，与实际质量相当。

图7 ISS舱段外壳
图7 ISS舱段外壳

辐射源

通量

考虑各类粒子均为各向同性入射，核心舱全长16.6 m，因此设置粒子均匀分布在半径为12 m的球面上，并保证粒子入射方向为球内，则每秒钟入射的粒子数目为：

F=A∫EbEaϕ(E)d E, A=4π R2F=A\int_{E_b}^{E_a}\phi\left(E\right)\mathrm{d}\ E,\ A=4\pi\ R^2

经测试，外壳能够屏蔽50 MeV以下的质子和100 keV以下的电子，结合CRÈME86等模型的能量区间，对于不同粒子设定各自的能区。

由于各类粒子通量相差极大，因此将其通量进行放缩，得到缩放因子（ZF），并保证每一类粒子的模拟数目不小于106。模拟得到剂量率后，乘以缩放系数和任务时长，即可得到任务期间的剂量。

表4 各类粒子的模型、能区与缩放因子（ZF）
粒子模型R2R^2EbE_b(MeV)EaE_a(MeV)模拟数目ZF
e−e^-j=j0e−EE0j=j_0e^{-\frac{E}{E_0}}0.981070.110923242099999.98096
pj=j0(E−β)−αj=j_0\left(E-\beta\right)^{-\alpha}0.99883501000100000042.85399415
HC1βα−1pc(pcpc+C2)C3\frac{C_1\beta^{\alpha-1}}{pc}\left(\frac{pc}{pc+C_2}\right)^{C_3}0.98915220100000100000012.85431092
HeC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9893622010000010000001.720188128
LiC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9975522010000010000000.01298897
BeC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9970822010000010000000.006691177
BC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9956922010000010000000.017105837
CC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9893622010000010000000.052286554
NC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9889222010000010000000.013576554
OC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9893622010000010000000.048846635
FC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9932322010000010000000.001123492
NeC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9901422010000010000000.008071585
NaC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9926322010000010000000.00186687
MgC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9897222010000010000000.01093646
AlC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9921722010000010000000.001940381
SiC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9887322010000010000000.008269257
PC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9918322010000010000000.00041389
SC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9887322010000010000000.00166099
ClC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9908722010000010000000.000372027
ArC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9943222010000010000000.000751798
KC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9896322010000010000000.000466712
CaC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9893522010000010000000.001161901
ScC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9925122010000010000000.000231554
TiC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.993422010000010000000.000829281
VC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9942422010000010000000.000404172
CrC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9931322010000010000000.00078409
MnC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9938522010000010000000.000570973
FeC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9935122010000010000000.00602806
CoC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9942622010000010000002.09E-05
NiC1e−C2E(1−e−C3E+C4)C_1e^{-C_2E}\left(1-e^{-C_3E+C_4}\right)0.9923122010000010000000.000292996

能谱与位置

对于一个在[a,b]\left[a,b\right]的分布，设其概率密度函数(PDF)为F(x)F\left(x\right)，要生成符合该分布的随机数，主要有两种方法：

  1. Inverse Transform Method（ITM）
    对PDF积分可获得累积分布函数(CDF)，记为y=C(x)=∫axF(u)d uy=C\left(x\right)=\int_{a}^{x}F\left(u\right)\mathrm{d}\ u，其反函数为C−1(y)C^{-1}\left(y\right)。设yy为[0,1)\left[0,1\right)上均匀分布的随机数，则x=C−1(y)x=C^{-1}\left(y\right)满足F(x)F\left(x\right)的分布。

  2. Acceptance-Rejection Method（ARM）
    生成一个均匀分布随机数X∼U(xmin,xmax)X\sim U\left(x_{min},x_{max}\right)，随后独立生成另一个均匀分布随机数Y∼U(ymin,ymax)Y\sim U\left(y_{min},y_{max}\right)，如果Y≤F(X)Y\le F\left(X\right)，则保留XX，保留下来的XX满足F(x)F\left(x\right)的分布。

ITM也被称为反演法，关键是需要获取累积分布函数的逆函数，效率较高；ARM本质上是一种模拟算法，效率很低，但是适应性更广，特别是对于一些复杂分布函数。

位置

发射位置要求在球面上均匀分布，因此可以采用ITM方法生成位置，可以得到θ\theta的反函数为C−1(z)=cos⁡−1(z)C^{-1}\left(z\right)=\cos^{-1}{\left(z\right)}，而ϕ\phi在[0,2π)\left[0,2\pi\right)上均匀分布。

方向 (12)

在Geant4中使用的坐标系是以中心为原点的三维坐标系（简称为W系），而我们随机生成的符合均匀分布的指向是在以粒子的位置矢量r⃗\vec{r}反方向作为z轴正方向建立的右手坐标系（简称为P系）中。

为了保证发射方向d⃗\vec{d}朝向球内，因此需要限制θ\theta的最大角度为90∘90^\circ。随后需要将d⃗\vec{d}转化到W系中，已知两个坐标系有两个公共点，在W系中表示为原点OO和发射点r⃗=(x0,y0,z0)\vec{r}=\left(x_0,y_0,z_0\right)，在P系中则为(0,0,ρ)\left(0,0,\rho\right)和原点O′O^\prime,则有：

[xyz]=λR[uvw]+[x0y0z0]⋯(1)\left[\begin{matrix}x\\y\\z\\\end{matrix}\right]=\lambda R\left[\begin{matrix}u\\v\\w\\\end{matrix}\right]+\left[\begin{matrix}x_0\\y_0\\z_0\\\end{matrix}\right]\quad\cdots\quad(1)

其中λ\lambda为尺度比例因子，设置为1；RR被称为罗德里格矩阵，定义反对称矩阵：

S=[0,−c,−bc,0,−ab,a,0]S=\left[\begin{matrix}0,&-c,&-b\\c,&0,&-a\\b,&a,&0\\\end{matrix}\right]

则R可表示为R=(I+S)(I−S)−1R=\left(I+S\right)\left(I-S\right)^{-1}，展开后则有：

R=11+a2+b2+c2[1+a2−b2−c2,−2c−2ab,−2b+2ac2c−2ab,1−a2+b2−c2,−2a−2bc2b+2ac,2a−2bc,1−a2−b2+c2]⋯(2)R=\frac{1}{1+a^2+b^2+c^2}\left[\begin{matrix}1+a^2-b^2-c^2,-2c-2ab,-2b+2ac\\2c-2ab,1-a^2+b^2-c^2,-2a-2bc\\2b+2ac,2a-2bc,1-a^2-b^2+c^2\\\end{matrix}\right]\quad\cdots\quad(2)

将两个公共点代入(1)，相减则可消去平移项，可得到：

[x2−x1y2−y1z2−z1]=R[u2−u1v2−v1w2−w1]\left[\begin{matrix}x_2-x_1\\y_2-y_1\\z_2-z_1\\\end{matrix}\right]=R\left[\begin{matrix}u_2-u_1\\v_2-v_1\\w_2-w_1\\\end{matrix}\right]

随后代入(2)中，可得到（u21=u2−u1u_{21}=u_2-u_1）：

[0,w21+z21,v21+y21w21+z21,0,−u21−x21−v21−y21,−u21−x21,0][abc]=[u21−x21v21−y21w21−z21]\left[\begin{matrix}0,&w_{21}+z_{21},&v_{21}+y_{21}\\w_{21}+z_{21},&0,&-u_{21}-x_{21}\\-v_{21}-y_{21},&-u_{21}-x_{21},&0\\\end{matrix}\right]\left[\begin{matrix}a\\b\\c\\\end{matrix}\right]=\left[\begin{matrix}u_{21}-x_{21}\\v_{21}-y_{21}\\w_{21}-z_{21}\\\end{matrix}\right]

代入本文中则有：

[0,−ρ+z0,y0−ρ+z0,0,−x0−y0,−x0,0][abc]=[−x0−y0−ρ−z0]\left[\begin{matrix}0,&-\rho+z_0,&y_0\\-\rho+z_0,&0,&-x_0\\-y_0,&-x_0,&0\\\end{matrix}\right]\left[\begin{matrix}a\\b\\c\\\end{matrix}\right]=\left[\begin{matrix}-x_0\\-y_0\\-\rho-z_0\\\end{matrix}\right]

令a=1a=1，则可得：

{b=z0+ρ−y0x0c=z0−ρ+y0x0\begin{cases} b=z_0+\rho-y_0x_0 \\ c=z_0-\rho+y_0x_0 \end{cases}

计算得到RR，继而根据(1)计算得到d⃗\vec{d}在W系中的表达式。

能谱

地球辐射带电子和质子的能谱较为简单，采用ITM方法，其CDF逆函数分别为：

{−E0ln⁡(C−yj0E0), C=j0E0e−a/E0(y+C)(1−α)j01−α+β, C=j0(a−β)1−α1−α\begin{cases} -E_0\ln{\left(\frac{C-y}{j_0E_0}\right)},\ C=j_0E_0e^{-a/E_0} \\ \sqrt[1-\alpha]{\frac{\left(y+C\right)\left(1-\alpha\right)}{j_0}}+\beta,\ C=\frac{j_0\left(a-\beta\right)^{1-\alpha}}{1-\alpha} \end{cases}

其中系数j0j_0均已归一化。GCR中各类粒子的函数较为复杂，因此采用ARM方法。粒子能谱的随机结果见图 9。

图9 各类粒子能谱随机结果
图9 各类粒子能谱随机结果

人体模型

模拟中可以使用ICRU球进行剂量的大概估算，也可以使用人体模型。辐射防护用虚拟人体模型经历了从简单的数学模型到复杂的面元模型的发展。

在Geant4的示例中，附带有MIRD体模，效果见图 8。MIRD模型是由橡树岭国家实验室（ORNL）基于白种人参考模型建立的数学模型，在外照射问题上的转换系数相当精确(13)，符合本文情况。

图8 Geant4中MIRD体模
图8 Geant4中MIRD体模

由于工作量的问题，在本文中只是用MIRD男性全身体模进行模拟，身体各器官的体积与质量见表 5。MIRD体模中，除肺之外的各类器官使用ICRU规定的软组织材料进行填充；肺使用专用材料；骨使用ICRU规定的骨质，而包裹在外的头部、躯干和四肢则使用脂肪、水等材料，使得人体总重约为86.8 kg。

表5 MIRD男性体模各器官信息
器官质量(g)体积(cm3\mathrm{cm^3})器官质量(g)体积(cm3\mathrm{cm^3})
头颅4595.484656.48骨盆897.554603.925
头骨1256.16845.219胃396.856402.124
大脑14511470.27上大肠429.223434.92
躯干43406.143982.3下大肠339.834344.345
(左/右)腿10252.110388.2脾脏173.625175.929
(左/右)手臂骨1217.78819.39胰腺60.241461.0411
(左/右)腿骨2080.291399.74肝脏1440.351459.47
上脊椎187.947126.462左肾脏142.102143.988
左肩胛骨152.24102.436右肾脏142.091143.977
右肩胛骨154.486103.947膀胱45.440946.044
(左/右)肾上腺7.751097.85398心脏362.795367.611
胸腺24.803525.1327左肺499.7431689.46
(左/右)锁骨20.309913.6657右肺499.8141689.7
小肠1006.851020.22甲状腺12.742112.9113
胸腔1021.01686.992男性生殖器226.245229.248
中低脊椎1120.57753.982(左/右)睾丸18.540618.7867

Geant4模拟流程

使用Geant4实现前文所述核心舱建模、人体建模与粒子定义，最终模型如图 10所示，图中从前往后以此为节点舱、过渡段1、小柱段、过渡段2与大柱段。人体模型站立在大柱段与过渡段2的交界处，如图 11所示。

对30种粒子依次模拟，通过Sensitive Detector记录进入人体模型的粒子，最终在User Run Action中将结果输出至文件中。

图10 核心舱建模效果
图10 核心舱建模效果
图11 人体模型位置
图11 人体模型位置

结果与分析

使用Python对数据进行处理，得到按粒子、器官分类的沉积能量信息，同时加入缩放因子和任务时长的因素，得到30天内航天员各器官沉积能量结果，按器官分类的结果如图 12所示，按粒子分类的结果如图 13所示。

图12 按器官分类的沉积能量结果
图12 按器官分类的沉积能量结果
图13 按粒子分类的沉积能量结果
图13 按粒子分类的沉积能量结果

吸收剂量

通过表 5中各器官的质量，可以计算得到各类器官的吸收剂量，按器官分类的结果如图 14所示，按粒子分类的结果如图 15所示。

图14 按器官分类的吸收剂量结果
图14 按器官分类的吸收剂量结果
图15 按粒子分类的吸收剂量结果
图15 按粒子分类的吸收剂量结果

吸收剂量受通量的影响很大。由于地球辐射带电子的通量远大于其余粒子，因此吸收剂量中地球辐射带电子占据的份额最大。

当量剂量

考虑粒子的辐射权重因子ωR\omega_R，由此可以计算得到各器官的当量剂量。按器官分类的结果如图 16所示，按粒子分类的结果如图 17所示。

由于电子、质子的ωR\omega_R较小，而重离子的ωR=20\omega_R=20，因此当量剂量中重离子占据了相当大的份额，电子的贡献显著减小。

图16 按器官分类的当量剂量结果
图16 按器官分类的当量剂量结果
图17 按粒子分类的当量剂量结果
图17 按粒子分类的当量剂量结果

有效剂量

结合不同器官的组织权重因子ωT\omega_T，可以得到人体的有效剂量。按器官分类的结果如图 18所示，按粒子分类的结果如图 19所示。

在计算中，男性生殖器的ωT\omega_T使用性腺的ωT\omega_T为0.08，躯干、头颅、四肢的ωT\omega_T使用皮肤的ωT\omega_T为0.01，其余器官的ωT\omega_T来自于ICRP 103报告。

图18 按器官分类的有效剂量结果
图18 按器官分类的有效剂量结果
图19 按粒子分类的有效剂量结果
图19 按粒子分类的有效剂量结果

分析

对结果的分类分析中可以发现一些现象：

左右不对称

由于粒子入射是各向同性，而核心舱模型也是各向同性的，理论上左右对称的器官应当有大致相同的能量沉积，但在结果中却发现左右腿等部分的沉积能量存在较大差异[4]，如表 6所示。

表6 左右部位吸收剂量差异值
部位左右差异
腿0.0259290550.03098487519.50
手臂骨0.0003338580.00024384836.91
腿骨0.0005267970.00075606443.52
肩胛骨0.0001273650.000095932.81
肾上腺00.0000463\
锁骨0.0000000770\
睾丸0.000001120\

对于如腿、手、肩胛骨等靠外侧的部位，可能是由于单次模拟的随机性造成的；而对于肾上腺，由于其位于人体内，由于肝脏、胰腺等器官并不左右对称，因此可能会导致左右不等。通过多次独立的模拟取均值的方式，可以验证这一问题的产生原因。

男性生殖器

在图 16和图 18中可以发现，男性生殖器的剂量是所有器官中最高的，特别是贡献了近2/3的有效剂量。这一结果可能是由于单次模拟的随机性造成的，同样需要多次独立模拟的验证。

总体而言，在此次30天的任务期间，航天员全程站立于大柱段与过渡段2的交界处所受到的全身当量剂量为287.6 mSv，而有效剂量为8.57 mSv。假设其每5年平均在轨工作12个月，则年平均剂量为20.568 a-1，略高于国家对于职业照射的剂量限值要求。

总结

本文使用SPENVIS建立了空间站天和核心舱轨道辐射环境的模型，并在Geant4中建立了核心舱简易几何模型与MIRD人体模型进行航天员辐射剂量的模拟。

模拟结果显示，在所有种类的射线中，地球辐射带的电子由于其通量大，贡献了绝大部分的沉积能量，但是由于电子难以穿透皮肤，因此主要贡献在表层。相反，质子、α粒子等重离子通量可观，能量高，因此在体内的器官中贡献了相当大一部分的有效剂量。

总体而言，任务期间航天员所受有效剂量基本符合国家规定的基本剂量限值，能够保障航天员的健康安全。当然本文使用的方法也存在局限性与诸多不足，后续可考虑改进的方向有：

  1. 建模
    1. 更精细化的人体模型
    2. 核心舱的细致材料
    3. 航天员卧室的建模
  2. 辐射环境
    1. 太阳活动极大/极小时
    2. 更长的任务时长
    3. 各向异性的环境
  3. 模拟
    1. 使用屏蔽计算工具配合Geant4
    2. 加入更多的入射事件
    3. 多次模拟对比
    4. 体模在不同位置的模拟（居留因子）

附录

参考文献

  1. 中国空间站建造进展情况新闻发布会召开_中国载人航天官方网站[EB/OL]. [2022-04-17]. http://www.cmse.gov.cn/xwzx/202204/t20220418_49553.html.
  2. SUN W, JIA X, XIE T, 等. Construction of boundary-surface-based Chinese female astronaut computational phantom and proton dose estimation[J/OL]. Journal of Radiation Research, 2013, 54(2): 383-397. DOI:10.1093/jrr/rrs100.
  3. 石苗. 空间站航天员关键器官辐射剂量学研究[D]. 南京航空航天大学, 2013.
  4. ERSMARK T, CARLSON P, DALY E, 等. Status of the DESIRE project: Geant4 physics validation studies and first results from Columbus/ISS radiation Simulations[J/OL]. IEEE Transactions on Nuclear Science, 2004, 51(4): 1378-1384. DOI:10.1109/TNS.2004.832648.
  5. ERSMARK T, CARLSON P, DALY E, 等. Influence of geometry model approximations on Geant4 simulation results of the Columbus/ISS radiation environment[J/OL]. Radiation Measurements, 2007, 42(8): 1342-1350. DOI:10.1016/j.radmeas.2007.06.001.
  6. ERSMARK T, CARLSON P, DALY E, 等. Geant4 Monte Carlo Simulations of the Belt Proton Radiation Environment On Board the International Space Station/Columbus[J/OL]. IEEE Transactions on Nuclear Science, 2007, 54(4): 1444-1453. DOI:10.1109/TNS.2007.896344.
  7. ERSMARK T, CARLSON P, DALY E, 等. Geant4 Monte Carlo Simulations of the Galactic Cosmic Ray Radiation Environment On-Board the International Space Station/Columbus[J/OL]. IEEE Transactions on Nuclear Science, 2007, 54(5): 1854-1862. DOI:10.1109/TNS.2007.906276.
  8. 程彭超, 闵锐. 近地空间辐射环境与防护方法概述[J]. 辐射防护通讯, 2017, 37(01): 14-21.
  9. BOURDARIE Sé, XAPSOS M. The Near-Earth Space Radiation Environment[J/OL]. IEEE Transactions on Nuclear Science, 2008, 55(4): 1810-1832. DOI:10.1109/TNS.2008.2001409.
  10. MATTHIÄ D, BERGER T, MRIGAKSHI A I, 等. A ready-to-use galactic cosmic ray model[J/OL]. Advances in Space Research, 2013, 51(3): 329-338. DOI:10.1016/j.asr.2012.09.022.
  11. ZHAO H, JOHNSTON W R, BAKER D N, 等. Characterization and Evolution of Radiation Belt Electron Energy Spectra Based on the Van Allen Probes Measurements[J/OL]. Journal of Geophysical Research: Space Physics, 2019, 124(6): 4217-4232. DOI:10.1029/2019JA026697.
  12. 韩梦泽, 李克昭. 基于罗德里格矩阵的空间坐标转换[J/OL]. 测绘工程, 2016, 25(04): 25-27. DOI:10.19349/j.cnki.issn1006-7949.2016.04.006.
  13. 徐玉海, 李桃生. 辐射防护用虚拟人体模型的发展研究[Z]. 中国辐射卫生: 卷 21. 2012: 119-120.

部分缩写列表

  1. DESCSS : Dose Estimation by Simulation of China Space Station
  2. DESIRE : Dose Estimation by Simulation of the ISS Radiation Environment
  3. ISS : International Space Station 国际空间站
  4. TLE : Two Line Elements 两行轨道参数
  5. SPE : Solar Particle Event 太阳粒子事件
  6. GCR : Galactic cosmic rays 银河宇宙射线
  7. GPS : General Particle Source 通用粒子源
  8. BIC : Bayesian information criterion 贝叶斯信息量准则
  9. ZF : Zoom Factor 缩放因子
  10. ITM : Inverse Transform Method
  11. ARM : Acceptance-Rejection Method

  1. http://www.cmse.gov.cn/gfgg/zgkjzgdcs/ ↩︎

  2. https://www.spenvis.oma.be/ ↩︎

  3. https://creme.isde.vanderbilt.edu/ ↩︎

  4. 差异η=max−minmin\eta=\frac{max-min}{min} ↩︎