利用香农信息熵的方法验证空间辐射屏蔽后各向同性

Evaluation of Space-Radiation-Spectrum Isotropy after Shielding Using Shannon Entropy Method

摘要: 空间辐射是制约载人深空飞行的重要环境因素之一，而人体器官所受的辐射剂量是评估辐射风险和研究防护方法的重要参数。但是，航天员深部器官的剂量只能通过计算得到，无法直接测量。器官剂量计算有直接和分步两种方法，后者可极大减少计算耗时。而经过飞行器屏蔽后舱内辐射能谱是否仍为各向同性，是决定分步计算方法准确性的关键。本文提出了一种利用香农信息熵理论对各向同性进行验证的方法。结果表明虽然屏蔽结构和统计位置会影响粒子分布，但出射粒子仍近似满足各向同性条件，并适用于分步计算的要求。

Abstract: Space radiation is one of the major factors that limit manned spaceflight. Organ doses are the important parameters for evaluating the health risks encountered by astronauts due to radiation, and are useful parameters in space radiation protection research. However, these organ doses cannot be directly measured in the bodies of astronauts, and can only be calculated. The direct and step-by-step calculation methods are used for calculating organ doses, and the latter method can reduce the calculation time considerably. But the validation of the particle isotropic distribution after shielding by spacecraft is critical in the step-by-step calculation method. In this paper, a method for testing the isotropy of the space radiation spectrum using Shannon entropy is proposed. The results demonstrate that although the shielding structure and statistical position can affect the isotropic distribution of particles, the spectrum of the emergent particles still satisfies an isotropic distribution approximately, and the requirements of the simulation calculation.

1. 引言

空间辐射是航天员进入太空不可避免遇到的危险因素，高能辐射粒子不仅足以造成DNA链断裂，也会使人体水分子发生电离，引起氧化等效应。空间辐射生物损伤程度，目前主要通过器官接受的辐射剂量进行评估 [1] [2] [3] [4]。由于航天员器官剂量无法直接测量，一般借助数字人体模型，通过蒙特卡罗方法对空间复杂的混合辐射场在器官的能量沉积进行计算 [5] [6]，这样的计算更加准确，通常有直接模拟和分步计算两种方法。

直接计算需要建立飞行器屏蔽模型和数字人体模型，通过计算得到经飞行器屏蔽后，进入舱内的辐射粒子在人体模型微小体元的能量沉积，统计属于同一器官体元的能量沉积后得到器官平均剂量。由于人体模型体元数量通常在百万个以上，且统计抽样次数较大，所以通常计算耗时较长，但是计算更为准确。

分布计算采用这样的思路：首先获得粒子注量–剂量转换数据，得到单位立体角单能单粒子造成的不同器官的剂量，该转换数据计算同样需要借助人体数字模型；然后计算得到空间辐射粒子经飞行器屏蔽进入舱内的能谱；最后，通过粒子注量–剂量转换数据和屏蔽后辐射能谱的加权求和，可快速得到器官剂量。不同的人体数字模型计算得到的粒子注量–剂量转换数据会有差异 [7]，对于空间辐射一般使用ICRP 123号报告的数据 [8]。

尽管直接计算方法的实现思路更为直观，但当飞行器模型有变化，如屏蔽厚度增减或几何结构改变，整个计算需要重新进行，因此花费时间较多。而分布计算方法，只需要对屏蔽后粒子谱进行重新计算，相较于直接计算节省了统计在粒子人体模型传输和能量沉积的计算步骤，大大节省了计算时间。但是，由于分布计算方法建立在辐射场各向同性的条件上，空间辐射谱(即入射能谱)通常认为是各向同性的，因此对屏蔽后辐射场的各向同性情况需要分析。本文借助香农信息熵的理论，提出了一种空间辐射场经屏蔽后各向同性的验证方法。

2. 分析方法

分步计算方法器官剂量通过下式得到：

$H_{T,R}={\displaystyle \underset{R}{\sum }{\displaystyle {\int }_E{h}_{T,Q,R}\left(E\right)\frac{\text{d}{\varphi }_E}{\text{d}E}\text{d}E}}={\displaystyle \underset{R}{\sum }{\displaystyle {\int }_E{Q}_{T,R}\left(E\right)d_{T,R}\left(E\right)\frac{\text{d}{\varphi }_E}{\text{d}E}\text{d}E}}$ (1)

H_T,R是特定器官组织特定辐射粒子的剂量当量，H_T,Q,R是能量为E的粒子造成的剂量当量，E是粒子能量，R指代粒子类型，T指代器官组织，Q_T,R(E)是品质因数，d_T,R(E)是注量–剂量转换系数， $\text{d}{\varphi }_E$ 是粒子通量。ICRP-123号报告发布了不同粒子在人体不同器官的注量–剂量转换系数，在计算时可以直接使用。

初始空间环境辐射场近似为各向同性，因此品质因数和注量–剂量转换系数均在各向同性条件下计算 [8]。当进行器官剂量计算时，必须验证屏蔽后辐射场各向同性情况。

2.1. 计算条件设置

2.1.1. 空间辐射环境能谱

本文选用1972年8月典型太阳粒子事件质子能谱作为环境输入 [9]。通过测量数据的分析拟合，质子能谱遵循下式幂指数变化规律：

$\frac{\text{d}\varphi }{\text{d}E}=2.2\times {10}^7{\text{e}}^{-\left(\frac{E-100}{30}\right)}$ (2)

由于质子能谱为各向同性，因此在使用蒙特卡罗方法进行粒子输运计算时，也需采用各向同性分布抽样，包括两个抽样步骤：入射粒子位置抽样和该位置的入射方向抽样，如图1所示。

Figure 1. Isotropic distribution sample steps

图1. 各向同性抽样步骤

图中A表示入射粒子能谱，B表示抽样球，其中中心的蓝色圆点表示统计位置。抽样步骤具体如下：

1) 设置一个半径为R的入射能谱抽样球，在球面上对入射粒子的空间位置均匀抽样，抽样球与屏蔽结构(简化起见，飞行器近似看作对称结构)同几何中心；

2) 在第一步确定的某个位置，对粒子入射方向进行抽样，为提高效率，入射方向需指向抽样球球心，并且遵照余弦分布被限定在特定角度α内，同时设定一个半径r，其定义为可包含屏蔽结构的最小球体的半径；

夹角α、r和R的关系如图1所示。最终由设定夹角α引起的影响可以通过归一化算法消除。设定X_r为真实值，X_s为仿真计算设定值，二者关系由下式表示：

$X_r=X_s\frac{N_r}{N_s}$ (3)

N_r是实际入射粒子个数，N_s是仿真计算设定值，通过统计误差确定；N_r/N_s定义为归一化系数P。在实际计算中，P通过如下经验公式获得：

$P=\frac{\Phi 4{\pi }^2{R}^2{\sin }^2\alpha }{N_s}$ (4)

φ是入射粒子积分谱通量，单位count⁻¹∙s⁻¹∙cm⁻²∙sr⁻¹，R是抽样球半径，α是入射角，N_s是抽样个数，本计算抽样个数均设定为600万次。

2.1.2. 屏蔽条件和粒子输运

本计算设定屏蔽为5 g/cm²的铝材料，这是载人飞行器典型质量屏蔽厚度。针对三种屏蔽结构进行了计算：球形、圆柱形和立方体。粒子输运计算使用Geant4软件，Geant4是一种面向对象的软件包，应用在空间、加速器、物理和生物物理等多领域 [10]。软件包包含7种物理作用：电磁相互作用、强相互作用、输运过程、衰变、光学过程、光轻子–强子过程、参数化过程。离子物理，包括电离激发和轫致辐射，是本计算种主要的物理过程。屏蔽后的出射粒子，统计在一个半径为1米的球面上，包括粒子位置和角度信息。

2.2. 信息熵及验证方法

2.2.1. 信息熵

熵的概念最早来源于热力学统计。1865年，R. Clausius提出了一个用于描述热传递程度的函数S。对于大多数热力学过程，这个函数可以如下表示：

$\Delta S=S_B-S_A\ge {\displaystyle {\int }_A^B\frac{\text{d}Q}{T}}$ (5)

dQ是表征系统从热源吸收热量的品质。

1872年，在Clausius的热动力学熵的基础上，Boltzmann提出了著名的Boltzmann方程：

$S=k\ln \Omega$ (6)

k是Boltzmann常数，Ω是微观状态个数。

由于熵的统计学特征，其概念如今已经从热力学，甚至是物理领域延伸开来，应用于如生命科学 [11]、信息学和社会科学 [12]。1948年，Glaude Shannon发表了文章《A Mathematical Theory of Communication》，首次提出了革命性的概念“信息熵” [13]。通过设定一组离散随机变量， $X=\left\{x_1,x_2,\cdots ,x_n\right\}$，信息熵可表示为如下形式：

$H\left(X\right)=-{\displaystyle \underset{i=1}{\overset{N}{\sum }}{P}_i\log \left(P_i\right)}$ (7)

N是总特征个数，P_i是X_i分布概率，因此 $0\le P_i\le 1$，且 ${\displaystyle \underset{i=1}{\overset{N}{\sum }}{P}_i=1}$。

以上公式可以看出，信息熵的值取决于变量的分布概率。根据定义，当每个变量的分布概率相同时，求和值最大。

$H\left(p_1,p_2,\cdots ,p_3\right)\le H\left(1/n,1/n,\cdots ,1/n\right)=\log n$ (8)

对于一组随机离散变量，当变量的分布概率相同时，信息熵的值最大。如果将统计球面上的粒子的位置信息(X, Y, Z)和角度信息(θ, φ)作为变量，则当屏蔽后粒子满足各向同性时，其求和值最大，因此可以通过该值判断各向同性情况 [14] [15] [16]。

2.2.2. 验证方法

通过上述分析，当各变量分布概率相等时，信息熵最大值 $H_{\max }=\log N$，具体到计算统计球面的变量时，N应当设置为比值的形式。统计位置信息时，以1 cm为间隔，则对于半径为1m的球，N_x，N_y，N_z设定为200，即将每个坐标轴分为了200个等间距的小格，统计屏蔽后出射粒子在该小格内的个数，再除以总抽样次数，得到该坐标值对应的分布概率P_i，然后利用公式(7)计算熵值；角度θ和φ进行同样的处理，θ和φ角度范围分别是0˚~180˚和−180˚~180˚，以1˚为间隔，则计算角度数据时，N_θ和N_φ设定为180和360。

3. 结果与讨论

3.1. 不同屏蔽结构各向同性分布的验证

图2所示为不同屏蔽结构条件，位置信息和角度信息算得的信息熵与最大值的比值随着抽样次数的变化规律，统计位置位于屏蔽结构几何中心。图中X，Y，Z表示位置变量与最大值的比值H/H_max，Theta和Phi表示角度变量与最大值的比值H/H_max。

(a) (b) (c)

Figure 2. Shannon entropy under different shielding structures. (a) Spherical shielding structure; (b) Cylindrical shielding structure; (c) Cubical shielding structure

图2. 不同屏蔽结构香农熵随抽样次数变化。(a) 球形屏蔽结构；(b) 圆柱体屏蔽结构；(c) 立方体屏蔽结构

由图2可以看出：

1) 球形结构对屏蔽后各向同性分布无影响，5个变量的H/H_max均超过99.99%；

2) 圆柱体结构对屏蔽后各向同性分布影响很小，除θ的H/H_max接近99.76%外，其他4个变量均超过99.99%；

3) 立方体结构对屏蔽后各向同性分布影响同样很小，除θ和φ的H/H_max接近99.79%和99.90%外，其他3个变量均超过99.99%。

由图2还可看出，屏蔽结构的改变不影响位置变量(X, Y, Z)的各向同性分布，但是会影响到角度变量(θ, φ)。这可能是因为入射能谱的抽样是在球面进行，且满足各向同性条件，当抽样粒子经过屏蔽结构时，因屏蔽结构和抽样球同几何中心，因此粒子的位置仍满足各向同性条件，但是其入射方向必然与屏蔽结构的几何形状有关，球形屏蔽结构自然满足各向同性条件，但圆柱体由于其结构在轴向方向同球形的差异，因此夹角θ的各向同性分布受到影响，同理，立方体结构在θ和φ两个方向均受到影响。但是根据分析结果，在抽样次数达到300万次时，大部分变量的H/H_max值均达到99.99%，抽样次数600万次时，不同几何结构屏蔽后的粒子能谱仍近似满足各向同性条件，引起的计算误差可以忽略不计。

3.2. 不同统计位置各向同性分布的验证

3.2.1. 球形屏蔽结构

图3所示为球形结构屏蔽后不同统计位置香农熵和最大值的比值随抽样次数变化规律。

(a) (b) (c)

Figure 3. Shannon entropy under spherical shielding conditions at different statistical positions. (a) Center; (b) Moved by 1 m along the X-axis; (c) Moved by 2 m along the X-axis

图3. 球形结构屏蔽后香农熵不同统计位置随抽样次数变化。(a) 中心位置；(b) 沿X轴移动1 m位置；(c) 沿X轴移动2 m位置

由图3可见，对于球形结构屏蔽，统计位置变化对各向同性分布无影响，5个变量的比值均超过99.99%。

3.2.2. 圆柱体屏蔽结构

图4所示为圆柱体结构屏蔽后不同统计位置香农熵和最大值比值随抽样次数变化规律。

(a) (b) (c) (d) (e)

Figure 4. Shannon entropy under a cylindrical shielding condition at different statistical positions. (a) Center; (b) Moved by 1 m along the X-axis; (c) Moved by 2 m along the X-axis; (d) Moved by 1 m along the Z-axis; (e) Moved by 2 m along the Z-axis

图4. 圆柱体结构屏蔽后香农熵不同统计位置随抽样次数变化。(a) 中心位置；(b) 沿X轴移动1 m位置；(c) 沿X轴移动2 m位置；(d) 沿Z轴移动1 m位置；(e) 沿Z轴移动2 m位置

由图4可见：

1) 在X和Z轴方向从几何中心移动1 m的位置进行统计，位置变量(X, Y, Z)和角度变量(θ, φ)的各向同性分布与在中心位置统计近似相同；角度变量θ的H/H_max可达到99.76%，其他4个变量的H/H_max超过99.99%。

2) 当统计位置沿X轴移动2 m时，X位置和φ的H/H_max不超过99.83%和99.85%；θ的H/H_max值降低到99.72%；其它两个变量的H/H_max值仍然可以超过99.99%。

3) 当统计位置沿Z轴移动2 m时，Z位置的H/H_max值只能达到99.87%，θ的H/H_max值降低到99.33%，其它3个变量的H/H_max值仍然可以超过99.99%。

以上结果可看出，对于圆柱体屏蔽结构，当统计位置位于几何中心，各向同性分布是最优的。当统计位置沿X轴移动时，X位置，θ和φ的各向同性分布变差；当统计位置沿Z轴移动时，Z位置和θ的各向同性分布变差。

3.2.3. 立方体屏蔽结构

图5所示为圆柱体结构屏蔽后不同统计位置的信息熵和最大值比值随抽样次数变化规律。

(a) (b) (c)

Figure 5. Shannon entropy under a cubical shielding condition at different statistical positions. (a) Center; (b) Moved by 1 m along the X-axis; (c) Moved by 2 m along the X-axis

图5. 立方体结构屏蔽后香农熵不同统计位置随抽样次数变化。(a) 中心位置；(b) 沿X轴移动1 m位置；(c) 沿X轴移动2 m位置

由图5可见：

1) 当统计位置沿X轴移动1 m时，θ和φ的H/H_max值降低到99.77%和99.86%；其它3个变量的比值仍然可以超过99.99%。

2) 当统计位置沿X轴移动2 m时，X的H/H_max值无法超过99.92%；θ和φ的H/H_max值降低到99.75%和99.64%；其它2个变量的比值仍然可以超过99.99%。

以上结果可看出，对于立方体屏蔽结构，当统计位置位于几何中心时，各向同性分布是最优的。当统计位置沿X轴移动时，X，θ和φ的各向同性分布变差。

4. 结论

本研究利用信息熵的方法对经屏蔽后的空间辐射能谱的各向同性进行了验证。尽管几何结构和统计位置对各向同性分布有影响，但是当抽样次数超过300万次时，位置(X, Y, Z)和角度(θ, φ)的H/H_max值均可以超过99%，并且当抽样次数继续增加时，H/H_max值并不改变。因此当抽样次数足够大时，经飞行器屏蔽后的舱内辐射能谱可近似认为是各向同性分布的，满足器官剂量快速计算的要求。所以，本文提出的验证方法可以用于分步剂量计算方法。

基金项目

本研究获得试验技术研究青年基金项目2017SY54C0401的资助。