地球周围的磁层空间中充满着大量的高能电子，这些高能电子主要被地球磁场俘获于距离地心（1.1～8）R_E（R_E为地球半径）的内磁层里。辐射带是内磁层中高能带电粒子组成的环绕地球的结构，在地球磁场的约束下，带电粒子以回旋、弹跳、漂移的形式在内磁层做周期性运动。辐射带主要包括内外两层可以长期稳定存在的结构。内辐射带处在距离地心（1.1～2.5）R_E范围内，其中既含有高能电子也含有高能质子，电子和质子的分布相对稳定；外辐射带处在（3～7）R_E处，该范围内高能质子的回旋半径较大，难以被磁场稳定捕获，因此外带粒子以高能电子为主，且与内辐射带相比外辐射带的电子动态特性更强。两层辐射带之间存在电子通量相对较低的槽区，该区域的电子容易与哨声波发生波粒相互作用，通过投掷角散射进入损失锥，在弹跳过程中注入大气发生损失。

高空强爆炸产生的裂变衰变β电子被地磁场俘获注入空间自然辐射带中，可以形成比自然辐射带电子强度高出几个量级的增强人工电子辐射带。辐射带增强效应是高空强爆炸的一种重要地球物理效应^[1-3]。人工辐射带持续时间长，影响范围广，会对穿越其中的空间飞行器的性能和寿命带来影响^[4-6]，因此，研究人工辐射带电子的长期扩散损失规律具有重要的应用价值。

目前已开展了许多基于数值模拟方法的辐射带研究。CCMC（Community Coordinated Modeling Center）开发了VERB（Versatile Electron Radiation Belt）模型，用于三维模拟天然辐射带中的相对论电子动力学，考虑了电子径向扩散、能量扩散和倾角扩散的耦合效应^[7]。美国戈达德太空飞行中心开发了CIMI（Comprehensive Inner Magnetosphere-Ionosphere）模型^[8]，该模型综合考虑了内磁层和电离层的耦合效应，能够模拟环电流、辐射带、等离子体层等多个物理过程的相互作用。这两个模型都主要应用于天然辐射带的相关研究。对于人工辐射带，西北核技术研究所对高空强爆炸产生的辐射环境进行了数值模拟^[9-10]，给出了包括核辐射、X射线、高空电磁脉冲和光辐射等瞬时环境的模拟方法，同时也对碎片云、附加电离层和人工辐射带等长期环境进行了数值模拟研究，得到爆后初期赤道面电子注量分布。武汉大学基于Stömer理论和高空强爆炸烟云模型^[6]，研究了高空强爆炸形成的人工辐射带的电子密度通量，分析了爆点纬度、高度和当量对人工辐射带的影响。西北核技术研究所牛胜利等人利用描述倾角扩散的福克-普朗克方程研究了高空强爆炸注入辐射带电子的大气扩散损失^[11-12]，建立了低L壳（这里L是赤道面上的空间映射点距地心距离与地球半径的比值）辐射带电子通量分布函数的数值计算方法，展示了大气散射对辐射带电子损失的影响。美国洛斯阿拉莫斯实验室针对Starfish高空强爆炸试验开展数值模拟研究，结果表明磁壳分裂与库仑碰撞耦合的新经典扩散是低L壳层人工辐射带电子长寿命的主导机制，并解决了Starfish人工辐射带衰变异常问题^[13]。

综上，虽然目前对人工辐射带已有一些研究报道，但大多关注的是爆炸初期或低纬度低L壳层人工辐射带，且文献[12]只对投掷角一个参量进行了研究。但高纬度高L壳层的人工辐射带更容易受到地磁活动的影响，故本文使用VERB3D模型，从径向、能量和投掷角三个参量对高纬度强爆炸产生的高L壳层人工辐射带电子的扩散过程开展模拟研究。首先我们模拟分析了2012年10月1日～11月31日的天然辐射带的演化，并与范艾伦探测器相对论电子质子望远镜（REPT）仪器的原位测量及A. Y. Drozdov等人在2014年的模拟结果^[14]进行了对比。结果表明，VERB程序的模拟结果可以较好地展现相对论粒子在地球磁层中的扩散过程。在此基础上，本文基于VERB程序构建了人工辐射带计算模型，对不同地磁水平下，高纬度地区强爆炸产生的人工辐射带的扩散过程进行了模拟研究。并分析了地磁活动水平对高L壳层人工辐射带电子扩散过程的影响规律。

1.   电子的扩散理论及计算模型

1.1   福克-普朗克（Fokker-Planck）方程

福克-普朗克方程是一个偏微分方程，用于描述粒子在相空间中的分布随时间的演化。在地球辐射带的背景下，该方程可用于描述高能电子在辐射带中的扩散过程^[15]，即

式中：f代表相空间密度（Phase Space Density, PSD），表示在相空间中粒子的分布；${D_{{I_i}{I_j}}}$为扩散系数，i和j是微分因子编号。

方程（1）可以展开为由电子的壳层L、动量p和投掷角α 来表示的方程形式^[7,16]，即

式中：μ和J分别是第一和第二绝热不变量，径向距离L与第三绝热不变量ϕ成反比，因此以L代替ϕ，p 是粒子的相对论动量，α 是投掷角。参数τ是粒子完成一次完整反弹运动的时间。${D_{LL}}$、${D_{pp}}$、${D_{\alpha \alpha }}$、${D_{p\alpha }}$是扩散系数，分别表示径向扩散、能量扩散、投掷角扩散及混合扩散的速率。

方程（2）中右侧的第一项描述了粒子在径向位置 L上的扩散。径向扩散通常由超低频（ULF）波引起，导致粒子在磁场中的径向位置发生变化。带电粒子的速度方向与地磁场方向有夹角，垂直于磁场的速度分量使带电粒子受洛伦兹力的影响做绕磁力线的回旋运动，平行于磁场的速度分量引导粒子沿磁力线运动。由于偶极场在两极附近磁力线汇聚，汇聚磁场产生的力使得粒子运动到场强高的区域时发生反射，该力足够使平行于磁场的速度减小为零。粒子随后沿磁力线作反向运动，粒子被反射回赤道区，穿过赤道向另一半球运动，在另一半球发生同样的反射，粒子将在两反射点间来回振荡，反射点又称磁镜点，粒子在两磁镜点之间四维运动被称作反弹运动。除了回旋运动和反弹运动外，还会围绕地球做漂移运动。径向扩散不改变粒子的第一和第二绝热不变量 μ 和 J，但会导致粒子的能量变化。第二项描述了粒子在投掷角 α 上的扩散。T（α）是与粒子反弹时间相关的函数。投掷角扩散由甚低频（VLF）波（如合声波和嘶嘶波）引起，导致粒子在投掷角上的分布发生变化。粒子在两磁镜点之间的反弹运动中，投掷角的变化会影响粒子的运动轨迹。投掷角扩散可以将粒子散射到损失锥中，导致粒子损失到大气层中。第三项描述了粒子在动量 p 上的扩散。能量扩散同样由VLF波引起，导致粒子的能量发生变化。粒子在磁场中的运动会导致能量的变化，特别是在辐射带的恢复阶段，能量扩散是粒子加速到相对论能量的重要机制。方程（2）中的${D_{p\alpha }}$项表示能量和投掷角扩散之间的耦合效应。混合扩散在高能粒子的动力学中也起着重要作用，特别是在低投掷角区域。最后一项f/τ项表示电子损失率。

1.2   VERB数值模型

VERB 是用于模拟地球辐射带中高能电子行为的数值模型。通过求解福克-普朗克方程，VERB 可以模拟辐射带在太阳活动周期和地磁活动变化下的长期演化。实际空间环境中存在的多种波粒相互作用可能会破坏三个绝热不变量，并导致径向输运、局部加速和电子损失。随Kp指数（Kp 指数是全球地磁活动指数）变化的径向扩散系数取自Brautigam和Albert在2000年提出的模型^[16]。合声波的参数取自Subbotin等人在2011年提出的模型^[17]。对于嘶嘶声波，采用了频率范围在100～2000 Hz、峰值频率550 Hz的真实数据^[18]。

为了简化计算并提高数值稳定性，VERB使用（V, K, L）坐标代替（μ，J，L）坐标^[19]。K是与投掷角相关的变量，V是一个与能量相关的变量。这种坐标变换可以将复杂的扩散方程简化为更易处理的形式，尤其是在处理混合扩散项时，传统的（μ，J，L）坐标在处理这些混合项时可能会遇到数值困难。（V, K, L）坐标则能够更自然地描述能量和投掷角之间的耦合关系，更好地处理混合扩散项，同时也降低了计算需要的网格数量。

2.   人工辐射带数值计算模型的构建

原始的VERB程序用来模拟天然辐射带的扩散过程，但高空强爆炸产生的人工辐射带比天然辐射带的电子通量密度有明显提高，所以对VERB中的初始辐射带电子密度分布进行调整和修改，从而建立人工辐射带的计算模型。

2.1   模型验证

首先验证VERB模型在进行电子扩散演化计算的正确性，本文模拟了2012年10月1日至11月30日期间的电子通量演化，并与REPT (Relativistic Electron-Proton Telescope) 仪器测量的数据及A. Y. Drozdov等人在2014年用VERB模型得到的模拟结果进行对比。

图1和图2分别是能量为2.0 MeV、3.6 MeV投掷角为85°的电子的通量变化，伪彩图中的白色虚线指的是等离子层顶的位置。图1和图2是本研究VERB的模拟结果，上面两张曲线图分别为地磁活动水平（Kp指数）及径向扩散系数（DLL scale）和合声波水平（chorus scale）。模拟结果分别与REPT仪器测量的数据和A. Y. Drozdov等人在2014年用VERB模型得到的模拟结果进行对比验证，即与文献[14]中的图7（a, b）和图（e, f）进行对比。图3为REPT卫星数据、A. Y. Drozdov等人的模拟和本研究模拟得出的等离子层顶的位置对比。

图  1  2012年10月1日至11月30日的能量为2.0 MeV投掷角为85°的电子通量

Figure  1.  The electron flux with an energy of 2.0 MeV and an injection angle of 85° from October 1, 2012, to November 30, 2012

下载: 全尺寸图片 幻灯片

图  2  2012年10月1日至11月30日的能量为3.6 MeV投掷角为85°的电子通量

Figure  2.  The electron flux with an energy of 2.0 MeV and an injection angle of 85° from October 1, 2012, to November 30, 2012

下载: 全尺寸图片 幻灯片

图  3  等离子体层顶变化对比

Figure  3.  Comparison of plasmapause variations

下载: 全尺寸图片 幻灯片

在2012年10月1日～11月30日期间一共发生了三次较为强烈的地磁暴，地磁暴会使等离子层顶向内收缩，同时使电子向内磁层扩散，在地磁暴发生时电子的通量密度显著提高。对比显示，模拟结果与REPT中的等离子体层顶的形状一致，可以验证模型中的等离子体模型的正确性；三次强地磁暴发生时，本研究进行的VERB模拟结果与探测器和A. Y. Drozdov等人的模拟结果中一致，都会产生电子通量显著增强、粒子向内辐射带扩散的现象。表明本次研究进行的VERB模拟结果是可靠的。

2.2   初始条件设置

数值模拟结果表明，在一定的条件下0.1～1 Mt TNT 当量的高空强爆炸预计在地球周围可形成电子通量密度比自然辐射带高3～4个量级的人工辐射带^[20]。人工辐射带的位置与爆点的磁纬度密切相关，爆点的磁纬度明显影响高空强爆炸产生的人工辐射带的中心位置；随着爆点磁纬度的增加，人工辐射带的中心位置逐渐向外推移，对应的L值也就越高。当爆点的纬度在45°左右时，强爆炸产生的人工辐射带位于内外辐射带之间的槽区；当爆点纬度进一步提高后，强爆炸产生的人工辐射带位于外辐射带。与内辐射带相比，外辐射带电子更具动态性，更易受到地磁活动的影响，故本研究主要针对高纬地区强爆炸形成的人工辐射带在地磁活动作用下的变化开展数值模拟研究。

VERB模型的计算网格分别有62×61×41个能量点、投掷角点和径向扩散点，能量和投掷角网格点按对数分布。初始条件在网格点上被计算为径向扩散方程的稳态解。方程（2）也需要设置六个边界条件，可分为三个变量对应电子PSD的上下限：下能量边界（10 keV）的PSD导数值设置为0，上能量边界（10 MeV）设置PSD为0；上投掷角边界（85°）的PSD导数值设置为0，下投掷角边界（0°）的PSD为0；下径向边界（L=1）处的PSD设置为0，表示大气中的损失参考Subbotin的方法^[16]，通过地球同步轨道上长期测量的PSD谱的平均形状计算得出上径向边界条件（L=5.5）所需的PSD。

以上是VERB3D中对天然辐射带电子初始分布的设定，本研究对天然辐射带进行的模拟沿用了这些设定。为了模拟分析高空强爆炸后产生人工辐射带的情况，在天然辐射带的基础上，本研究在L=5.5处加入了大量电子，以此来模拟强爆炸后大量电子进入辐射带的情况，下文对于人工辐射带的模拟都遵循此方法，具体设定如下：

在L=5.5处搜索能量在2.5～3 MeV的电子对应的V、K参量值，将PSD数组中这部分电子对应的PSD提高10个量级，电子通量相应提高了3～4个量级，最终注入部分的电子可近似为高斯分布，通量密度峰值在2.75 MeV处为10¹³（cm²·s·MeV·sr）⁻¹，注入电子的能量范围为2.5～3 MeV。除了L=5.5的其他位置与天然辐射带的设定一致。

3.   计算结果与分析

3.1   地磁活动水平对径向扩散过程的影响

按照2.2节的初始条件设置方式，我们模拟了人工辐射带电子的演化过程。图4和图5分别为不同地磁活动水平下，人工辐射带粒子在L=2.6～6.6范围内的扩散情况，这里分别选取了具有不同地磁活动特征的2012年10月7日至10月17日和10月21日至10月31日两个时期进行对比分析。图中最上面的曲线为Kp指数， 10月7日至10月17日期间Kp指数波动更为剧烈，整体活动水平较高，表明此期间发生了多次强地磁暴，10月21日至10月31日期间Kp指数虽然也有频繁波动，但整体水平相对较低。

图  4  地磁活动强烈时期的人工辐射带扩散情况

Figure  4.  The diffusion of artificial radiation belts during periods of intense geomagnetic activity

下载: 全尺寸图片 幻灯片

图  5  地磁活动平静时期的人工辐射带扩散情况

Figure  5.  The diffusion of artificial radiation belts during periods of quiet geomagnetic activity

下载: 全尺寸图片 幻灯片

在图4和图5两种场景下，人工辐射带电子均从L=5.5附近的径向位置注入。之后受Kp指数的影响，图中白色虚线代表的等离子层顶和ULF波都会发生变化，在波粒互作用下，注入的电子发生径向扩散，分布在更大范围的L范围。扩散分布的边界均与等离子体层顶位置密切相关，原因是等离子体层是内外辐射带的分隔区，外辐射带的电子很难穿过等离子体层进入内带，而等离子体层受地磁活动水平影响，因此图5中层顶位置不会发生显著变化，变化范围为L=4.5～6，而图4中的层顶甚至会到达L=3以下，而辐射带电子不会向内扩散穿过等离子层顶至槽区，但向外扩散至磁层顶，因此相较于地磁平静时期，地磁活动更强烈时，人工辐射带的扩散范围更大。除了空间分布，通量密度图5降低的更慢，可知地磁活动会加快人工辐射带电子的扩散过程，减小辐射带寿命。

图6和图7为人工辐射带中，投掷角为85°，磁地方时（MLT）为6时3 MeV电子的径向分布情况，图中给出了L=4.5、5.0、5.5、6.0、6.5五处的电子通量随时间的变化。

图  6  地磁活动强烈时期人工辐射带的径向扩散情况

Figure  6.  The radial diffusion of artificial radiation belts during periods of intense geomagnetic activity

下载: 全尺寸图片 幻灯片

图  7  地磁活动平静时期人工辐射带的径向扩散情况

Figure  7.  The radial diffusion of artificial radiation belts during periods of quiet geomagnetic activity

下载: 全尺寸图片 幻灯片

图6中随着10月8日第一次强磁暴发生，人工辐射带注入的粒子迅速扩散，在径向上呈现梯度分布，L值与电子通量分布负相关，此时电子径向分布较为稳定。此后在10月13日发生了又一次的强磁暴事件，各个径向位置上的电子通量有小幅度的减少，总体分布无明显变化。在图7所示地磁平静期的（2012年10月21日至31日），电子同样是由高密度处向低密度处扩散，但是扩散速率比地磁活动水平剧烈时的扩散速率慢很多。经过十天的扩散，仍未达到稳定分布状态。

通过以上分析不难看出，当电子在径向上的通量密度差距较大时，电子分布不稳定，强地磁活动加快扩散过程，会使其迅速达到稳定分布状态。当电子的通量密度稳定分布时，强地磁活动仅会使其通量密度小幅度降低，并不会显著改变其通量分布特征。

3.2   地磁活动水平对能量扩散和投掷角扩散过程的影响

上述分析表明，Kp指数对径向扩散的影响主导了电子在L壳间的输运过程，但其能量和投掷角的演化还需进一步考虑局地波粒相互作用的影响。与径向扩散不同，能量和投掷角扩散由甚低频（VLF）波（如合唱波、等离子体层嘶嘶波）直接驱动，其扩散系数${D_{pp}}$和${D_{\alpha \alpha }}$的时空变化与背景等离子体参数（如等离子体密度、波功率谱）紧密相关。下文将详细讨论这些局地扩散过程的物理机制及其数值实现方法。

辐射带电子在地磁场中运动时，三个绝热不变量（μ，J，ϕ）在缓变场中近似守恒。径向扩散要求粒子沿恒定的μ和J的路径移动，导致能量E和投掷角α随L壳变化。若直接在（E，α，L）的正交网格上计算径向扩散，会破坏μ和J的守恒性，导致非物理的数值误差。因此VERB引入非正交的能量和投掷角网格，通过分离径向扩散与局部扩散的网格体系，确保物理过程的精确描述。如图8所示，局部扩散网格为在每一L壳上构建(E，α)的非正交网格。

图  8  非正交的能量-投掷角网格

Figure  8.  Non-orthogonal Energy-Pitch Angle Grid

下载: 全尺寸图片 幻灯片

图9、图 10分别为2012年10月7日至17日（地磁活动强烈期）和2012年10月21日至31日（地磁活动平静期），电子在L=5.5处的能量和投掷角扩散情况。注入电子会同时向高能和低能方向扩散，也会向高投掷角和低投掷角方向扩散。然而，由于不同区域的损失机制存在显著差异，这种扩散过程呈现出明显的非对称性。注入电子的由于损失锥的存在，低投掷角的电子更容易进入损失锥，高投掷角的电子在辐射带中会更加稳定；高能电子的动量较大，在绕地球做漂移运动的过程中会相对更快地进入高L壳层，而进入高L壳层的电子更容易发生磁层顶逃逸，因此低能电子会更稳定。这种选择性损失机制导致辐射带电子总体上呈现出向低能区和高投掷角区扩散的特征。

图  9  2012年10月7日至17日期间电子的能量和投掷角演化(L=5.5)

Figure  9.  The evolution of electron energy and pitch angle distributions during October 7-17, 2012 (L=5.5)

下载: 全尺寸图片 幻灯片

Kp指数作为表征全球地磁活动的关键参数，通过改变波粒相互作用的背景环境，间接影响能量扩散和投掷角扩散的速度与空间分布。Kp指数升高直接改变激发VLF波（合唱波、嘶嘶波）的波功率和空间分布，进而影响扩散系数。对比两图可知，地磁活动会显著增强电子的能量和投掷角扩散，图9中的电子在经历10月8日的强磁暴后，电子迅速完成向低能区和高投掷角区扩散，分布趋于稳定；而在10月21日至31日，经过十天的扩散，电子能量仍然集中在2.0～3.5 MeV，能量和投掷角均是缓慢扩散。因此，地磁活动水平会通过大幅度增强波粒互作用来显著加快粒子的能量和投掷角扩散。

图  10  2012年10月21日至31日期间电子的能量和投掷角演化(L=5.5)

Figure  10.  The evolution of electron energy and pitch angle distributions during October 21-31, 2012(L=5.5)

下载: 全尺寸图片 幻灯片

4.   结　论

本文采用VERB3D数值模型模拟了人工辐射带的演化过程，并分析了不同空间天气水平下人工辐射带的扩散情况和磁层中等离子体层顶的变化。通过分别对比地磁活动强烈期和平静期进行的模拟结果对比，我们发现等离子体层顶在地磁活动平静期相对稳定，变化范围在L=4.5～5.5；在强地磁活动背景下，等离子体层会受磁暴影响向内收缩，最大可收缩至L=3处。这种等离子体层的动态变化对人工辐射带的径向演化产生重要影响：一方面，等离子体层的向内收缩导致人工辐射带的径向分布范围扩大；另一方面，增强的波粒相互作用显著加速了人工辐射带电子的能量和投掷角扩散过程，使电子迅速向低能区和高投掷角区扩散。本文的研究结果有助于对人工辐射带的演化进行预测，并分析人工辐射带对空间飞行器的长期影响。