近地空间质子环境严重影响着人类航天活动，给地球近地轨道运行的航天器带来危险，使卫星载荷仪器上的敏感元件受到损伤，影响仪器乃至整颗卫星的使用寿命^[1-3]。影响航天器的高能带电粒子除了被地磁场捕获的辐射带粒子外，主要的还有来自太阳的太阳宇宙线和来自银河系的银河宇宙线。因此，准确描述空间高能质子的分布情况，获取卫星轨道宇宙线质子能谱，对卫星的合理使用和我国的深空探测事业十分重要。通过对宇宙线几十年的持续观测和研究，已经有一些模式对宇宙线的通量和空间分布情况进行描述^[4-5]。如20世纪90年代初，Badhwar和O’Neill基于1951年至1992年气球和卫星探测数据提出了Badhwar-O’Neill银河宇宙线模型^[6]；1998年，莫斯科州大学经过对宇宙线长达四十年的直接观测后提出了新的银河与太阳宇宙线模式，该模式描述了带电荷数为1～28的粒子在1 MeV～1 GeV能量范围的通量情况^[7]；2004年P. M. O’Neill利用要素/同位素成分高级探测器（ACE）的数据重新修正的Badhwar-O’Neill 银河宇宙线模型等^[8]。比较接近实用且目前仍被广泛采用的模式是Adams模式，该模式总结了大量地面、气球和卫星的探测数据，非常全面地描述了不同能量的粒子在不同的太阳活动条件时的通量水平^[9-11]。Nicola Tomassetti等进行了数据驱动的分析，在太阳周期的银河宇宙射线（GCR）通量的时间依赖性。利用AMS-02和PAMELA的空间GCR观测数据，对GCR扩散平均自由程进行了全球统计推断，确定了GCR扩散平均自由程的刚性和时间依赖性，给出了GCR质子扩散平均自由程的结果，讨论了它们在粒子输运基本过程中的解释及其与日球等离子体动力学的关系^[12]。虽然对于宇宙线的探测取得了很大的成功，但是目前应用最广泛的仍然是质子能谱的经验公式，其计算速度慢，平均误差在200%以内，局部区域误差可达450%^[13]，部分结果仍需要结合实际卫星的探测结果进行优化，无法满足航天工程需求，所以本文基于机器学习算法以卫星实测数据为驱动进行宇宙线能谱反演研究，达到对宇宙线质子能谱的快速响应。

1 原理分析及数据集建立

从宇宙线研究早期，人们就开始广泛地研究宇宙线的时间变化。作为初级近似，海平面的宇宙线通量通常被当作常量，但多年观测资料的分析发现宇宙线的强度存在着涨落。这些涨落部分是周期性的，其余则是非周期性的。已经观测到的周期性变化有几小时变化、周日变化、27天变化、年变化和11年变化等周期性特征，这些全部都是太阳活动的调制过程所造成的。非周期性变化也有起源于太阳的，如Forbush下降^[14]。这是一种短期宇宙线通量全球性剧烈下降，下降幅度可达10%～50%，通常与太阳耀斑关系密切。宇宙线随地磁场的变化主要表现有纬度效应、经度效应、东西效应及辐射带的形成。初级宇宙线以及它在大气层中产生的大多数次级粒子都是带电的，宇宙线带电粒子穿越地磁场到达近地轨道的过程中会受到地磁场的作用，因而轨道会发生偏转；同时当宇宙线进入大气层后，会与大气中各种分子碰撞，发生弹性和非弹性散射、康普顿散射、光电效应、轫致辐射等过程，并产生μ子、中微子、电子、正电子、中子、质子等次级粒子。这些粒子能够到达地球表面，从而被地面探测器监测到^[15-16]。基于此，本文建立以下三个数据集：地面中子数据集、太阳活动观测数据集、宇宙线探测数据集。

地面中子探测数据来自中子探测数据库（NMDB），该网站提供了来自世界各地站点的中子探测器探测结果的访问。本文共选取了59个站点的数据，时间跨度由1957年至今。所有数据都经过了大气压力修正及各个站点的修正（建筑物的不同、探测器数量、探测器类型等）。图1给出了所有站点的磁纬度及其截止刚度的分布，可以看出站点主要分布在高磁纬度地区，且截止刚度与磁纬度的绝对值之间的相关系数 ${r}^{2}=-0.908\;7$，证明两者具有较强的负相关性。

图 1. 所有站点的磁纬度与其截止刚度的分布

Fig. 1. Distribution of magnetic latitude and cutoffrigidity of all stations

下载图片 查看所有图片

太阳黑子数据来自国际数据中心（WDC）从全世界约85个观测点收集的太阳黑子数据并进行平均，时间跨度由1818年至今；地磁活动指数K_p数据来自德国地球科学研究中心（ GFZ），时间跨度由1932年至今。

宇宙线中高能质子探测数据来自于美国国家环境信息中心提供的地球静止环境业务卫星（GOES）自1974年7月以来的探测数据。该卫星位于地球静止（GEO）轨道，环境中的高能质子主要来自于磁层外，其搭载的空间环境探测器中的高能粒子传感器能够对质子、α粒子及电子作出响应。本文根据探测数据的特点，首先对于数据集 $ \{\left({x}_{1},{y}_{1}\right),\cdots ,\left({x}_{t},{y}_{t}\right)\} $进行归一化处理，即

$ {x}_{i}'=\dfrac{{x}_{i}-{x}_{\mathrm{m}\mathrm{i}\mathrm{n}}}{{x}_{\rm{max}}-{x}_{\mathrm{m}\mathrm{i}\mathrm{n}}} $ (1)

式中： $ {x}_{i} $为待归一化的数据， $ {x}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{a}\mathrm{x}\{{x}_{1},\cdots ,{x}_{t}\} $， ${x}_{\rm{min}}=\mathrm{m}\mathrm{i}\mathrm{n}\{{x}_{1},\cdots ,{x}_{t}\}$。归一化后的 $ {x}_{i}' $分布在[0,1]区间，由于粒子通量的变化范围跨越多个数量级，对通量数值取对数，形成训练所需训练集。

因此，本文构建的模型输入分别为：（1）地面中子计数率；（2）太阳黑子数：太阳黑子数观测历史长、数据量多，可以准确反映太阳活动变化。（3）地面探测站的磁纬度：磁纬度会影响不同能量的宇宙线粒子在不同地区的分布。（4）地磁活动指数K_p：由各种情况导致的地磁扰动同样会对宇宙线粒子的地区分布产生影响。

在航天工程领域，大多数计算与应用都是基于完整的粒子能谱输出，但是XGBoost模型的输出只能是离散值，所以本文选择粒子能量作为模型的一个额外输入，粒子通量作为模型输出。通过不断改变输入能量，从而得到连续的输出通量，最终得到连续的粒子能谱。基于以上输入输出，本文的数据集由地面中子观测数据、太阳活动观测数据和宇宙线探测数据三个部分构成，数据集的总体情况如表1所示。

2 模型构建

本节使用遗传算法来获得XGBoost的初始参数最优组合，具体步骤如下：

① 设置初始种群的规模为30，最大进化代数为40，初始化每个个体的参数，组成父代子群 $ {C}_{1} $。

② 计算每个个体的适应度，使用转盘式选择法挑选出15个优势个体，组成精英子群 $ {P}_{1} $。

③ 在 $ {P}_{1} $中随机两两结合，经过交叉和变异操作产生15个后代，组成后代子群 $ {Q}_{1} $。

④ 精英子群与后代子群混合为新的父代子群 $ {C}_{2} $，重复步骤②～③。

⑤ 迭代结束后，选取种群中适应度最高的个体，该个体所代表的方案即为最佳方案。

训练过程中，种群的最佳个体（即适应度最高）的模型参数方案的均方误差(MSE)随进化代数的变化以及进化结束后的最佳方案参数如图2所示。

图 2. 模型的MSE随迭代次数的变化

Fig. 2. MSE of the model changes with the number of iterations

下载图片 查看所有图片

本文使用MSE、均方根误差（RMSE）和平均相对误差（MRE）作为反演精度的总体衡量标准，计算公式为

$ {E}_{\mathrm{M}\mathrm{S}\mathrm{E}}=\dfrac{1}{N}\displaystyle\sum _{i=1}^{N}{\left(\hat {{y}_{i}}-{y}_{i}\right)}^{2} $ (2)

$ {E}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt{\dfrac{1}{N}\displaystyle\sum _{i=1}^{N}{\left(\hat {{y}_{i}}-{y}_{i}\right)}^{2}} $ (3)

$ {E}_{\mathrm{M}\mathrm{R}\mathrm{E}}=\dfrac{1}{N}\displaystyle\sum _{i=1}^{N}\dfrac{\left|\hat {{y}_{i}}-{y}_{i}\right|}{{y}_{i}}\times 100{\text{%}} $ (4)

式中：N为样本总数。上文提到通量数据在进行归一化处理时取10为底的对数，则以上衡量标准也是针对取对数后的输出。

为了测试遗传算法对模型训练结果的改进，针对使用XGBoost默认参数、进化开始时最佳个体参数、进化结束后最佳个体参数这三种情况，分别使用上式计算MSE、RMSE、MRE，结果如表2所示。默认参数方案的MSE为0.605，种群进化结束时最佳个体的MSE为0.499，降低了17.6%，决定系数 $ {r}^{2}=0.837 $，在测试集上的平均相对误差为26.9%。作为对比，传统的辐射环境模型AE8和AP8的相对误差通常在200%以内，某些区域甚至会达到450%，表明模型在反演宇宙线环境上具有较好的精确度。

3 结果及分析

为了全面地比较本文所使用模型的反演效果，使用支持向量回归（SVR）、误差反向传播算法（BP）神经网络和长短期记忆（LSTM）神经网络同时建立对宇宙线的反演模型，并与XGBoost模型进行对比^[17]。本文对SVR超参数设置为：核函数为高斯径向基函数，惩罚系数为1×10³， $ \gamma =0.1 $；对BP神经网络的结构设置为：隐含层共3层，神经元个数分别为64、32、16，激励函数为ReLU，训练次数为1000次；对LSTM网络结构设置为：地面中子探测数据的输入变更为过去24 h内的中子计数率（间隔为1 h），隐含层结构与其余参数与BP神经网络相同。

表3给出了计算机电脑在处理器为12th Gen Intel(R) Core(TM) i9-12900KF，3187 MHz的条件下，这些模型的各项指标的对比，从表中可以看出，本文所使用的GA-XGBoost模型的各项指标综合表现最优。原因有以下几个方面：（1）其平均相对误差MRE仅为12.5%，在四种模型中表现为佳，说明该模型预测的数据集稳定性最高，与测试集的差距最小；（2）均方根误差RMSE为0.706、均方误差MSE0.499，比SVR支持向量机模型和BP神经网络模型精度高，但次于深度学习网络LSTM模型的0.361/0.131，说明该模型的拟合程度较好，预测准确度也较高；（3）训练时间为135 s，远低于其他模型的训练时间，说明该模型占用最少的时间和资源成本，能更加灵活地应对突发事件，更加适用于航天工程应用。

为了更直观地分析各种算法的适用性，分别在太阳活动极小年和极大年选取两个时间点——2010年7月21日与2015年7月21日，进行GEO轨道质子能谱反演，图3和图4给出了实际质子通量与模型反演质子通量的对比。从图3中可以看出，在太阳活动极小年，质子微分通量随着能量的增加先下降后上升，SVR未能捕捉到这个特征，只表现出指数形式下降的整体趋势；BP捕获到了部分特征，但是与测量值偏差较大，精度低；LSTM和XGBoost模型能完整地捕获能谱变化特征，都具有较高的反演精度。从图4中可以看出，在太阳活动极大年，质子微分通量也是随着能量的增加先下降后上升，但是下降点右移，SVR仍然无法捕捉到这个特征，只表现出指数形式下降的整体趋势；BP捕获到了部分特征，但是整体向右偏移，与测量值偏差较大，精度低；LSTM和XGBoost模型能完整地捕获能谱变化特征，具有较高的反演精度，且XGBoost模型与实测值符合更好。

图 3. 太阳活动极小年模型的预测值与探测值及其他算法的对比

Fig. 3. Fluxes calculated by our model and SVR, BP, LSTM models, in comparison with flux data detected by GOES10 detector in solar minimum

下载图片 查看所有图片

图 4. 太阳活动极大年模型的预测值与探测值及其他算法的对比

Fig. 4. Fluxes calculated by our model and SVR, BP, LSTM models, in comparison with flux data detected by GOES13 detector in solar maximum

下载图片 查看所有图片

图5和图6给出了模型的预测通量与宇宙线模型CREME96的对比，从图中可以看出，模型预测值与宇宙线模型CREME96在高能部分符合得较好，在低能部分与AP8模型符合较好。因为在GEO轨道，卫星监测到的辐射质子通量不仅包括宇宙线高能质子通量，还包括地球外辐射带中的低能质子通量，所以预测值/探测值都接近宇宙线质子通量与外辐射带质子通量之和。

图 5. 太阳活动极小年模型的预测值与CREME96模型、AP8模型的对比

Fig. 5. Calculated fluxes comparison among our model, the CREME96 model, and AP8 model in solar minimum

下载图片 查看所有图片

图 6. 太阳活动极大年模型的预测值与CREME96模型、AP8模型的对比

Fig. 6. Calculated fluxes comparison among our model, the CREME96 model, and AP8 model in solar maximum

下载图片 查看所有图片

综上所述，本文建立的GA-XBGoost算法能够准确地绘制出GEO轨道的质子微分能谱，更加有助于设计和优化航天器的辐射防护系统，提前预知可能的辐射事件，并采取措施确保卫星正常运行的稳定性和可靠性。

4 结　论

在大数据时代，针对海量数据的机器学习是国内外研究的重点。本文针对主要的辐射环境研究目标——宇宙线环境，分析了影响宇宙线的各个因素，包括地磁场、太阳活动等。本文使用极端梯度提升决策树（XGBoost）作为反演模型，使用遗传算法对XGBoost的超参数进行了优化，基于GOSE探测数据、全球地面中子探测器的探测数据、太阳黑子数与地磁活动指数建立了训练集。构建的模型能较好地实现对平静时期同步轨道宇宙线环境的反演，MSE为0.499。

然而，受限于探测数据，模型的反演结果无法准确延伸到GeV量级。随着辐射探测技术的发展，如果能够获得更多高能粒子的探测数据，并应用在模型的构建中，则能够扩充模型的结果，提高模型的精度。对于宇宙线的各个组成部分，考虑到现有的公开探测数据的体量和范围，本文仅选择了宇宙线质子作为研究对象，但是重离子在辐射环境对人体造成的辐射剂量中占据了较大的比例，如果能够收集到关于宇宙线重离子的探测数据，就能扩大模型的应用范围，使模型更好地为航天员的生命安全服务。