基于卡方分布函数的辐射带电子蝴蝶状投掷角分布的优化判别方法

罗琼，倪彬彬,2*，曹兴，顾旭东，易娟，郭英杰，郭德宇，周若贤，龙敏义

1 武汉大学电子信息学院空间物理系，武汉 430072 2 中国科学院比较行星学卓越创新中心，合肥 230026

0 引言

地球辐射带电子受多种物理过程影响，其空间分布对不同物理过程的响应各不相同.电子的投掷角分布(Pitch Angle Distribution，PAD)及其演化作为电子空间分布的基本特征可为研究这些物理过程提供重要信息.地球辐射带电子投掷角分布主要分为以下三种：①90°峰值(90°-peak)分布，表现为通量在90°附近达到峰值，并随着投掷角向0°和180°变化而迅速减小；②平顶状(flattop)分布，通量在以90°为中心的较宽投掷角范围内变化不大；③蝴蝶状(butterfly)分布，通量在90°附近有极小值、在较低投掷角达到峰值(Fritz et al.，2003；Gannon et al.，2007；Thorne et al.，2010；Gu et al.，2011；Chen et al.，2014；Zhao et al.，2014；Ni et al.，2015).

不同类型的电子投掷角分布往往与不同物理过程相关联，其中，蝴蝶状分布的形成机制相对复杂且更引人关注.多项研究表明，与地球磁层的地方时不对称性相关联的磁层顶阴影效应和漂移壳分裂效应可能使电子形成蝴蝶状分布(Roederer，1970；Roederer et al.，2014；Sibeck et al.，1987；Selesnick and Blake，2002；Turner et al.，2012)，其原因在于赤道附近弹跳的电子更易在日侧漂移更远从而穿过磁层顶离开地球磁层.波粒相互作用也是蝴蝶状分布的可能成因之一，研究发现磁声波可以通过朗道共振机制使电子形成蝴蝶状分布(Horne et al.，2007；Mourenas et al.，2013；Xiao et al.，2015；Li et al.，2014，2016；Yang et al.，2017；Fu et al.，2019；Zhou et al.，2020；Ni et al.，2020)，而合声波则可能导致蝴蝶状分布的迅速平坦与恢复(Yang et al.，2016).此外，非绝热散射如磁力线曲率散射(Artemyev et al.，2015)或向外的绝热输运(Su et al.，2010)等都可能形成蝴蝶状分布.

上文列举的判别方法在原则上都是对特定投掷角的通量比值进行限定，这样挑选的电子分布尽管满足数值标准，但其完整“通量-投掷角”剖面不一定符合蝴蝶状特征，这为准确研究电子蝴蝶状分布的现象学规律及其背后物理机制带来了一定困难.因此，本文建立了一个基于卡方分布函数(Chi-square Distribution Function，简称2分布)的模型，对符合传统方法挑选标准的事例进行二次筛选，以实现对电子蝴蝶状分布的优化判别.1.1节以Ni等(2016)中判别方法为例，通过具体电子投掷角分布事例阐述传统方法的不足.1.2节详细介绍了建模过程，2.1和2.2节则分别从具体事例、统计结果上验证了传统方法结合该模型后对辐射带电子蝴蝶状投掷角分布的判别效果的提升.

1 蝴蝶状投掷角分布判别方法

1.1 传统方法判别电子蝴蝶状分布的不足

本文使用了范艾伦A星上相对论电子质子探测器(Relativistic Electron Proton Telescope，REPT)(Baker et al.，2013)测量的2013年至2015年期间磁壳(L-shell)大于3空间区域内2.1 MeV电子通量随投掷角变化的数据，时间精度约为11 s.根据研究目标，本文对Ni等(2016)中的方法进行细微调整后作为判别电子蝴蝶状分布的传统方法，相较于调整前，筛选可用数据的标准更加严格，而从可用数据中判别蝴蝶状分布的标准不变.使用传统方法对REPT提供的17个投掷角电子通量数据进行筛选，从符合该方法标准的所有电子投掷角分布事例中随机选取8376个事例，组成传统方法判别的电子蝴蝶状投掷角分布样本集，该方法具体标准如下：

①通量等于零视作无效通量，筛选第2～16个投掷角的通量均有效(>0)的数据.REPT提供固定17个投掷角的通量数据，第2～16个投掷角分别为15.9°、26.5°、…、153.5°和164.1°.其中，场向(一般为20°以下或160°以上)通量、中间投掷角通量、90°通量是判断电子分布类型的重要信息.

②j(90°)<β×javg(90°-α:90°+α)，其中，j(90°)表示90°的电子通量；javg(a∶b)为投掷角a～b之间通量的平均值，α在10°～45°之间变化，最终取使javg(90°-α:90°+α)最大的α值；β为“蝴蝶状分布指标”，表征90°通量小于其他投掷角通量的程度，本文中β=0.95.

③90°投掷角以下(上)的通量峰值应介于20°～80°(100°～160°)之间，以此排除场向分布或90°峰值分布.

本文选取了3个符合以上标准的2.1 MeV电子通量随投掷角分布事例，来说明传统方法判别蝴蝶状分布的不足，如图1所示.图1a展示了卫星2015年3月22日于磁壳(L-shell)4.9，磁地方时(MLT)0.2 h，磁纬(MLAT)-2°观测到的一个清晰、标准的电子蝴蝶状分布事例，图中，90°通量为极小值，在47.6°和132.4°处达到峰值，随后通量降低，通量剖面关于90°基本对称，整体呈清晰、标准的“蝴蝶翅膀”形；尽管图1b、1c在数值上均满足挑选标准，但图1b的剖面仅在0°～90°之间较好地符合“蝴蝶状”特征，而在90°～180°之间无明显起伏，倾向于各向同性分布；图1c展示了一个存在明显90°处通量“凹陷”和90°以外通量峰值但通量剖面更倾向于各向异性分布的事例，正是它起伏不定的通量使之满足判别标准，造成了误判.

图1 3个REPT测量的符合传统蝴蝶状分布判别标准的2.1 MeV电子投掷角分布事例(a)清晰的蝴蝶状投掷角分布；(b)混合型投掷角分布；(c)各向异性投掷角分布.Fig.1 Three examples of 2.1 MeV electron PADs measured by REPT that meet the traditional identification criteria of butterfly PAD(a)A clear butterfly PAD.(b)A mixed PAD.(c)An anisotropic PAD.

以上结果表明，传统方法判别的电子分布虽满足数值标准，但其通量剖面并不一定符合蝴蝶状特征，且这种“非蝴蝶状分布”事例占比不小，在本文样本集中约占42.9%.人工筛选判别后的电子通量剖面是解决该问题的方法之一，但工作量大幅增加.基于此，本文建立了一个模型，对符合传统方法挑选标准的事例进行二次筛选，挑选出通量剖面真正具有蝴蝶状特征的电子分布事例.

1.2 基于卡方分布函数的判别模型

模型的建立主要分为以下2步：

(1)选择数学函数，利用函数曲线描述蝴蝶状投掷角分布的“双峰”特征.

理想蝴蝶状分布的“通量-投掷角”剖面关于90°对称，下文仅以0°～90°为例进行分析，其通量峰值处投掷角(记为αpeak)一般可取20°～80°之间的任意值，峰值两侧通量的下降速率各不相同，因此选择极大值两侧曲线斜率不等的卡方分布函数来描述其“双峰”特征.

图2展示了使用卡方分布函数建模，描述“双峰”特征的具体过程.若随机变量X服从自由度为n的卡方分布，则其概率密度函数为(Lancaster and Seneta，2005)

(1)

函数曲线如图2a所示，随着n的增大，曲线逐渐出现极大值点(记为xpeak)，函数斜率也随之变化，因而，改变自由度n即可得到不同坡度的曲线，以此对应不同通量峰的理想蝴蝶状剖面.以n=9为例(如图2b)，xpeak≈7，取一段包含xpeak且起点为xstart=3.5终点为xend=14的曲线，将函数f(x)在[xstart,xend]之间的值一一映射到投掷角[0°,90°]之间的通量上.由于卡方分布函数极大值左侧的曲线必定比右侧更倾斜，而理想蝴蝶状分布通量峰值两侧的斜率并没有固定大小之分，故采用两种映射方式，此处引入参数isFlip来表征映射方式：①isFlip=0(如图2c)，保持曲线不变，xstart和xend分别对应0°和90°投掷角，相应αpeak≈30°，通量剖面关于90°对称；②isFlip=1(如图2d)，将f(x)在[xstart,xend]之间的曲线翻转，xstart对应90°，xend对应0°，相应αpeak≈60°，且剖面关于90°对称.从图2c、2d可以看出，对于任意一段包含xpeak的卡方曲线，“翻转”与“不翻转”映射的剖面完全不同.

图2 基于卡方分布函数建立判别模型的流程图(a)不同自由度(n)的卡方分布函数曲线示意图；(b)自由度为9的卡方分布函数曲线示例，xstart和xend为截取的曲线段的两个端点；(c)通过将所截取的位于xstart和xend之间的曲线段直接映射到0°到90°的投掷角范围内或(d)先翻转曲线段再进行映射得到的两个不同的蝴蝶状剖面；αpeak代表通量峰值处的投掷角；n、xstart、xend和isFlip组成该判别模型的变量.Fig.2 Flow chart of the process of establishing the identification model based on the chi-square distribution function(a)Function curve graphs of the chi-square distribution with different degrees of freedom (n);(b)The function curve of chi-square distribution with 9 degrees of freedom.xstart and xend are associated with the two endpoints of the selected curve segment.(c)and (d)Two different butterfly profiles obtained by directly mapping the selected curve segment between xstart and xend to the pitch angle range of 0°～90° or by first flipping the curve segment over and then mapping.αpeak represents the pitch angle corresponding to the flux peak.n,xstart,xend and isFlip constitute the model variables.

由此，得到基于卡方分布函数建立模型的四个参数(n,xstart,xend,isFlip)，其中，n为卡方分布的自由度，xstart和xend分别为卡方曲线的起点和终点，isFlip表征从卡方曲线到“通量-投掷角”剖面的映射方式，值可取0(不翻转)或1(翻转).改变这四个参数，即可映射不同的通量剖面.

(2)遍历模型参数(n,xstart,xend,isFlip)，尽可能全地描述多种理想蝴蝶状投掷角分布.

不同参数(n,xstart,xend,isFlip)表征的卡方曲线映射出的“通量-投掷角”剖面各不相同，根据以下条件，遍历出各参数的合理取值组成参数集，参数集越完备，模型所描述的蝴蝶状分布越全面.

①3≤n≤20，n以步长1遍历3～20.当n<3时，卡方分布函数没有极大值，不可取，而20自由度卡方分布对于本模型参数集的建立已经足够.

②xpeak

③当n<15时，0≤xstart

⑤对任一卡方曲线映射的通量剖面，应有j′(0°)≤0.5.一个0°附近通量相对整个通量剖面普遍较高(尽管没有达到通量峰值)的电子分布更倾向于场向分布，故限制j′(0°)≤0.5.

根据以上条件，共遍历出9869组(n,xstart,xend,isFlip)，每组分别映射一个理想蝴蝶状“通量-投掷角”剖面，根据通量值计算剖面间的相关系数(Correlation Coefficient，记作CC)，仅保留相似剖面(CC>0.95)中的一组，最终精简得到一个包含320组参数的基于卡方分布函数的判别模型(模型下载网址http:∥doi.org/10.6084/m9.figshare.14616738).

图3是该模型所描述的320组理想蝴蝶状“通量-投掷角”剖面的简单示意图，不同颜色代表不同组别.每组αpeak各不相同，为了更好的展示细节，将剖面按αpeak从小到大变化划分为3个子图：20°<αpeak<40°(图3a)；40°<αpeak<60°(图3b)；60°<αpeak<80°(图3c).以图3b为例，随着线条颜色由绿向橘变化，曲线极大值点逐渐从40°向60°移动，且即使极大值位于同一投掷角，极大值两侧的斜率也各不相同，与之相比，图3a中剖面的αpeak相对较小，0°～αpeak之间曲线斜率的变化相对简单，αpeak～90°之间的相对复杂，曲线覆盖的“面积”也更多，图3c遵循类似规律.综上，这320组参数组成的模型能基本完备地描述各种理想状态下的电子蝴蝶状投掷角分布.

图3 基于卡方分布函数的判别模型所模拟的320组理想蝴蝶状投掷角分布的剖面图不同颜色代表不同组别，按照峰值通量处的投掷角将这320组剖面划分为3个部分展示：(a)20°～40°，(b)40°～60°，(c)60°～80°.Fig.3 320 groups of color-coded profiles of ideal butterfly pitch angle distributions simulated by the identification model based on the chi-square distribution functionDividing these 320 profiles into 3 parts according to the peak-flux pitch angle:(a)20°～40°,(b)40°～60°,(c)60°～80°.

2 结合传统方法与基于卡方分布函数的模型对电子蝴蝶状分布的判别效果

为了更准确地判别电子蝴蝶状投掷角分布，本文在传统方法的基础上，创新性地引入了基于卡方分布函数的数学模型，本节将分别从具体事例和统计结果上验证并展示传统方法结合卡方分布模型后对电子蝴蝶状分布的判别效果.

2.1 具体事例分析

模型使用方法如下：将电子观测投掷角输入模型(只输入通量有效的投掷角点)，返回320组模型通量值，计算归一化观测通量与每一组归一化模型通量的相关系数，其中最大值记为CCmax，这表征了电子观测“通量-投掷角”剖面能在多大程度上符合模型描述的某一理想蝴蝶状“通量-投掷角”剖面，这也是模型判别蝴蝶状分布的根本依据.

对本文1.1节中3个符合传统方法挑选标准的电子分布事例应用该模型，结果如图4，图中电子分布事例和布局均与图1一致，纵坐标为最小-最大归一化通量，蓝色点虚线为观测通量，红色为CCmax对应的那一组模型通量，该组模型参数与CCmax值均展示在子图标题中.图4a中CCmax=0.997，表明电子观测通量剖面与模型(n,xstart,xend,isFlip)=(12,1,18,0)所映射的理想蝴蝶状通量剖面高度相似，图中90°处“凹陷”、通量峰值点均匹配很好；相比之下，图4b中CCmax仅为0.711，虽然该电子观测通量只在0°～90°之间符合蝴蝶状分布，模型仍较好地模拟出了～60°处的峰值特征；图4c的电子分布更倾向于各向异性型，CCmax极小，仅为0.391，说明该事例几乎不与任何理想蝴蝶状剖面相似.综上，该模型能在一定程度上模拟出观测通量的特征，如通量峰值点、90°通量点、0°通量点等，能依据CCmax对电子投掷角分布进行量化判别.

图4 3个REPT测量的符合传统蝴蝶状分布判别标准的电子投掷角分布事例图上展示了归一化观测通量(蓝色点虚线)和与CCmax相关联的由卡方分布模型模拟的归一化模型通量(红色点虚线)，CCmax定义为观测通量与320组模型通量的相关系数的最大值，n、xstart、xend和isFlip是模型变量，图中3个投掷角分布事例与图1中一致.Fig.4 Three examples of electron PADs measured by REPT that meet the traditional identification criteria of butterfly PADIncluding:the normalized observation fluxes (blue dot-dashed line),and the normalized model fluxes (red dot-dashed line)associated with CCmax and simulated by the chi-square distribution model.CCmax represents the maximum of the correlation coefficients between observation fluxes and 320 groups of model fluxes.n,xstart,xend and isFlip are the model variables.These three PAD events are consistent with Fig.1.

2.2 统计分析结果

将CCmax≥0.8定为模型判别标准，基于1.1节中的传统方法随机选取的包含8376个电子投掷角分布事例的数据集来验证引入卡方分布模型后对电子蝴蝶状分布判别效果的提升程度.通过人工筛选，将这8376个符合传统方法判别标准的样本分为4785个“真蝴蝶状分布”和3591个“非蝴蝶状分布”，应用该模型，计算每个电子分布事例的对应CCmax，并按不同样本集画出柱状图，如图5所示.需要说明的是，对这8376个传统方法判别的事例进行人工再筛选时，并不要求通量剖面严格符合理想蝴蝶状分布，即允许0°～90°(或90°～180°)之间除通量峰值外存在其他细微的通量极大值，但要求通量整体符合蝴蝶状变化趋势，要求0°(或180°)附近通量相对整个通量剖面处于中等或较低通量值水平，要求通量剖面关于90°具有一定对称性.

图5a展示了符合传统方法挑选标准且人工判断为蝴蝶状分布事例的对应CCmax柱状图，横坐标为CCmax，纵坐标为相应CCmax区间内事例数与该样本集总事例数的比(记为Ratio)，具体事例数展示在柱状图上方.图中，Ratio随CCmax的增大而增大，于CCmax≥0.95处达到峰值，约为65.7%，表明该样本集中大部分电子分布与模型描述的理想蝴蝶状剖面具有较高相似度.类似地，图5b为符合传统方法但人工判断不是蝴蝶状分布的对应CCmax分布图，Ratio在0.6～0.7和0.9～0.95处有两个峰值，整体变化趋势不明显，且这些“非蝴蝶状分布”事例中仍有一部分CCmax≥0.8，原因主要有以下2点：

(1)CCmax计算的是观测通量与模型通量的整体相似度，故无法准确排除在某个或某几个投掷角处通量激增或锐减的事例，而这种可能由仪器观测误差引起的通量突变事件是可以人工筛选的.

(2)该模型可能无法判别某些高度类似蝴蝶状分布的事例.由于平顶状分布(或场向分布)和蝴蝶状分布在定义上并非完全互斥，甚至某些平顶状分布可以视作90°峰值分布与蝴蝶状分布的过渡阶段(Gannon et al.，2007)，故实际存在一些具备多种分布结构特征的“复合事例”，这些事例并不属于蝴蝶状分布，但它们的CCmax可能很大.具体类型有：①某些20°～160°之间较符合蝴蝶状分布、但0°(或180°)附近通量处于相对较高水平(虽未达到峰值)，故更倾向于人工判定为场向分布的通量剖面(参考附图1a中90°以下通量部分)；②某些在0°、90°和180°符合蝴蝶状特征但在整体上倾向于平顶状分布的通量剖面(如附图1b)；③某些在90°以下和90°以上分别符合场向分布和蝴蝶状分布特征的通量剖面(如附图1a)，或分别符合平顶状分布和蝴蝶状分布特征的通量剖面(如附图1c).这些高度类似蝴蝶状分布的事例的CCmax可能很大，进而造成模型误判，在现有观测通量投掷角精度条件下，这种误判可以通过增加其他标准得到一定改善，具体在本文的第3节结论中阐述.

附图1 3个符合传统判别方法且具有较大CCmax值但人工判定为“非蝴蝶状分布”的电子投掷角分布事例图片布局与图3一致.Appendix Fig.1 Three examples of electron PADs that meet the traditional identification method,have relatively large values of CCmax but are manually identified not as butterfly PADsThe picture layout is consistent with Fig.3.

REPT除提供具有固定17个投掷角的通量数据外，同时也提供具有变化36个投掷角的数据，且17个投掷角的数据是根据36个投掷角的数据处理得到的.为验证投掷角分辨率对模型判别效果的影响，对这8376个事例同时刻的36个投掷角通量数据做类似处理，结果如图5d—f.由于17个投掷角与36个投掷角的通量剖面并不相同，人工筛选的结果也不会完全一样.与图5b不同，图5e中柱状图峰值集中于0.4～0.6，之后Ratio随CCmax的增大而减小，CCmax≥0.95时Ratio接近于0，可见，提高投掷角精度后，CCmax表征观测通量剖面与理想蝴蝶状剖面相似性的准确度大大提高，定量上看，误判率从仅使用传统方法时的51.2%变为结合卡方分布模型后的23.7%，下降了27.5%，而丢失率仅上升了0.7%.对比图5c和图5f可知，提高投掷角精度能显著提高卡方分布模型对电子蝴蝶状分布的优化判别效果.

图5 不同数据集中电子投掷角分布事例的CCmax柱状图Ratio为CCmax网格内电子投掷角分布事例数占该数据集总事例数的比值.数据集包括：(a)和(d)符合传统判别方法且人工判定为“真蝴蝶状分布”的电子投掷角分布事例；(b)和(e)符合传统判别方法但人工判定为“非蝴蝶状分布”的电子投掷角分布事例；(c)和(f)以CCmax≥0.8作为模型判别标准时，“单独使用传统判别方法”与“结合传统方法及卡方判别模型”的误判率(RF)和丢失率(RM)对比.(a)—(c)与(d)—(f)对应于相同8376个投掷角分布事例但具有不同投掷角分辨率的数据集.Fig.5 The distribution of electron PADs as a function of CCmax for different data setsY-axis is the ratio between the number of PAD events in each CCmax bin and the total event number of that data set.The data sets are:(a)and (d)the electron PAD events that meet the traditional identification method and are manually identified as butterfly PADs,(b)and (e)the electron PAD events that meet the traditional identification method but are manually identified not as butterfly PADs.(c)and (f)Comparisons of false identification rate (RF)and missing rate (RM)when using different identification methods,i.e.,“using the traditional identification method alone”and “combining the traditional identification method and the chi-square identification model”.Setting CCmax≥0.8 as the model identification criterion.(a)—(c)and (d)—(f)correspond to the electron PAD datasets of the same 8376 PAD events but with different pitch-angle resolutions.

3 总结与讨论

传统的判别电子蝴蝶状分布的方法主要是对特定投掷角的通量比值进行限定，这样挑选的事例其通量剖面并不一定符合蝴蝶状特征.为了改善这些不足，本文建立了一个包含320组参数的基于卡方分布函数的模型，对传统方法挑选的事例进行二次筛选，并使用范艾伦卫星上REPT仪器提供的L-shell>3空间区域内两种投掷角分辨率的电子通量数据对该模型方法进行验证，主要结论如下：

(1)本文建立了一个包含320组较完备的理想蝴蝶状分布模拟参数的模型，并通过计算模型描述的蝴蝶状“通量-投掷角”剖面与实际观测的电子“通量-投掷角”剖面的最大相关系数(CCmax)，创新性地量化了电子投掷角分布符合蝴蝶状分布的程度.

(2)使用卡方分布模型对传统方法的判别结果进行二次筛选，可以更准确地挑选出真正符合蝴蝶状分布特征的通量剖面，进而优化对电子蝴蝶状投掷角分布的判别效果.

(3)统计结果表明，卫星观测数据的投掷角精度越高，卡方分布模型对判别蝴蝶状投掷角分布的优化效果越好.相较于仅使用传统方法的结果，结合卡方分布模型后在使用17个投掷角和36个投掷角通量数据时误判率分别下降了12.6%和27.5%，而丢失率基本不变，接近于0.

引入卡方分布模型在一定程度上优化了对电子蝴蝶状投掷角分布的判别效果，但仍有一些改进空间，具体包括以下几方面：

(1)改进建模过程.本文基于理想蝴蝶状分布在通量峰值两侧斜率不等的特征，经过一定对比后选择卡方分布函数作为建模基础，并规定模型描述的蝴蝶状剖面关于90°对称，若要改进模型，可尝试使用其他函数(如两个正态分布函数相加)、以更精细的步长和标准遍历模型参数、允许模型通量关于90°有一定范围的不对称等.

(2)增加其他判别标准.如2.2节对图5b的描述所示，某些高度类似蝴蝶状分布的场向分布和平顶状分布可能较难被该模型识别，但可通过增加标准进行一定优化，如：①通过减小β值来加深90°处通量“凹陷”程度进而降低平顶状分布被误判的概率；②对0°(或180°)通量与峰值通量、90°通量的相对大小进行限定，进一步剔除场向分布的干扰.

本文的研究结果是基于REPT数据得到的，不同行星环境中电子投掷角分布的发展趋势会随磁场结构的改变而有所不同，如地球辐射带中稳定、闭合的磁场结构更利于电子向双边损失锥分布发展，而近火空间中场向分布更为普遍，故不同行星环境中卡方分布模型的应用效果也会不同.此外，CCmax≥0.8为本文经验标准，可根据具体研究需求提高该模型判别阈值，在有选择性地牺牲丢失率的前提下降低对电子蝴蝶状分布的误判率.

致谢感谢范艾伦卫星科研组提供数据.范艾伦卫星电子通量数据来源于(https:∥www.rbsp-ect.lanl.gov/data_pub/).本文建立的卡方分布判别模型可通过网址(http:∥doi.org/10.6084/m9.figshare.14616738)下载.

附录A