单电荷高能离子与较冷中性原子之间发生共振电荷交换而生成高能中性原子(ENA，Energetic Neutral Atoms)是太阳系空间等离子体中一种比较普遍的现象(Fahr et al.,2007).高能中性原子卫星成像测量可广泛应用于地球和其他行星的磁层与电离层，以及太阳日冕与太阳风等离子体的探测(Gruntman，1997).对于地球内磁层，环电流中的高能离子通过与中性冷地冕氢间的共振电荷交换生成高能中性原子，它们不再受磁场与电场的约束，以原有能量和动量几乎沿直线轨迹从磁层逃逸；这些中性原子携带着发生电荷交换处的能量离子性态(成份、能量、方位等)的消息，是环电流能量离子的远程信息使者；利用专门的探测仪器可在卫星上对环电流ENA进行全球遥感成像.20世纪90年代后期，瑞典发射了首个ENA成像卫星Astrid，用专门设计的成像仪，从1000 km的低轨道上对磁层进行了短期的扫描成像(Barabash et al.,1997; C:son Brandt et al.,1997; C:son Brandt et al.,2001).21世纪初，美国实施了首个对磁层进行多手段全面成像的IMAGE卫星计划(Burch，2000，2003)，其中包括多能段ENA成像，取得有价值的数据，深化了人们对磁暴和亚暴全球过程的认识(Perez et al.,2001; C:son Brandt et al.,2002; Mitchell et al.,2003; Pollock et al.,2003)，揭示了暴时非对称环电流的最强离子通量会出现在午夜后，而不是传统认为的总是出现在午夜前等新现象(Fok et al.,2003).中国科学家在21世纪初即开展了卫星ENA成像的模拟研究(Shen and Liu，2002)，中国双星Double Star(Liu et al.,2005)的极轨道星TC-2上携载了高能段ENA成像仪，并由TC-2 ENA成像数据反演获得赤道面环电流离子分布(Lu et al.,2008，2010).

虽然ENA成像是空间等离子体的一种重要遥测诊断工具，但是从ENA二维图像定量反演其生成离子源的三维空间(加上投掷角和能量共五维)分布是一个具有挑战性的难题.IMAGE和TC-2卫星在极区上空几个地球半径以外对内磁层进行ENA全球成像，都是单颗星观测，一次扫描所生成的二维图像的观察视角有限，这使得人们难以由此提取得到可靠的离子三维空间分布.自从1987年Roelof(Roelof，1987)利用ISEE-1卫星观测到的ENA(它们是被原本用来测量带电粒子的仪器观测到的)，采用环电流赤道离子分布零阶模型等简化假设，首次反演得到粗略的环电流离子分布以来，已发展了多种反演方法以获取ENA的离子源分布，包括使用参数化离子分布模型求解模型参数的方法(Roelof and Skinner，2000; Perez et al.,2000，2012)和仿照“有约束的线性反演”(Twomey，1977)直接解得赤道离子通量强度(C:son Brandt et al.,2002; DeMajistre et al.,2004; Lu et al.,2010)，以及利用环电流理论模型的数据融合(data assimulation)方法(Nakano et al.,2008)，等等.但至今所进行的反演都是对赤道面二维环电流离子分布的重建，为了使求解得以进行，大都采用绝热不变假设，借此约束赤道面与其他纬度上离子分布之间的定量关系而将三维分布的求解降为二维.然而磁暴主相/亚暴期间的环电流离子分布往往会打破绝热不变.美国TWINS(Two Wide-angle Imaging Neutral-atom Spectro-meters)是世界上首个对地球磁层进行立体成像的星座(McComas et al.,2009；Goldstein & McComas，2013)，它利用相互交叉莫尔尼亚轨道上的两颗卫星，在北极上空有利位置，同时从不同角度对磁层ENA进行全球扫描成像；这为由多幅二维图像利用体像素(Voxel)CT技术重建环电流离子能谱的三维空间分布提供了有利条件.

本文以下将介绍我们发展的磁层ENA发射CT反演方法，以及在此基础上，为克服多卫星ENA成像存在的仪器不一致性造成反演误差增大甚至迭代发散无解的问题，而引入的差分ENA体像素CT反演方法；然后利用离子分布模型和TWINS卫星实测成像几何构形，对差分反演方法进行数值模拟检验；并利用该方法由卫星实测ENA数据重建2012年7月15日磁暴期间的环电流离子三维分布与能谱，将其与THEMIS卫星当地测量结果进行比较，验证此方法的可靠性.

CT探测与反演的基本思想可以简要地概括为：从物理上，对于一个有限的空间区域，如果有某种波(包括辐射)从所有各种不同方向穿过该区域，波与该区域介质相互作用的路径积分效应能够被观测到，波与介质相互作用的物理规律已知，那么可以由所有这些路径积分效应，确定该区域介质(或辐射源)的分布；从数学上来说，CT探测与反演分别对应获得Fredholm型积分方程组和求解该方程组，理论基础是Radon变换与求逆(Radon，1917；Deans，1984)，及其三维推广(Chiu et al.,1980；Tuy，1983；徐继生和邹玉华，2003).路径积分值在CT领域称为投影，投影与待重建的目标分布之间，由Furiuer切片定理或称投影定理(Kak and Slaney，1988；徐继生和邹玉华，2003)直接联系起来.

磁层卫星ENA成像测量到的高能中性原子强度(计数)C，是到达仪器的ENA微分通量J_ena被仪器响应函数A加权的多重积分：

式中A_i是仪器因子，表示t时刻，第i个像素对能量为E，以仰角η和方位角β入射的ENA粒子的响应.

ENA微分通量J_ena含有两个部分，一是较高高度上能量离子(对于我们所研究的情况，主要是环电流离子)与稀薄地冕氢发生电荷交换生成的ENA，称为高高度发射(High Altitude Emission)，记作J_HAE；另一部分是到达外逸层底较浓密大气中的高能离子与氧原子发生电荷交换生成的，称为低高度发射(Low Altitude Emission)，记作J_LAE，即

其中

式中σ_O表示高能离子与中性氧原子共振电荷交换截面，σ_strip表示ENA被碰撞剥离又生成离子的交换截面.这部分涉及复杂的多重碰撞和多种类型的电荷交换，需要单独专门处理.本文主要关心较高高度上环电流ENA发射，其微分通量可近似表示为

其中

上两式中n_H是冷的背景氢中性原子数密度，σ_H表示高能离子与中性氢原子共振电荷交换截面，其上标表示离子的种类，J_ion是沿视线方向的高能离子微分通量，s是从卫星出发的视线，沿着由仰角η和方位角β决定的方向，积分的上限是视线首次和外逸层底相交的地方，如果没有相交则理论上为+∞，但针对实际的探测源区，积分上限是有限的.

由(4a)式可见，到达仪器的ENA微分通量J_ena是沿观测视线上低能中性原子数密度与高能离子强度(通量)及其电荷交换截面的乘积的路径积分.这样，如果电荷交换截面和地冕氢分布已知，且假定高能离子源存在于有限的空间，背景地冕氢为光学薄的稀疏介质，并进一步假定中性原子成像测量在待测区域之外的多点从不同角度同时进行，那末从原则上就可以采用CT反演技术重建高能离子源的分布.此情况可归属于文献Bates(Bates et al.,1983)中所述遥感CT，在某种意义上类似医学中的单光子发射CT(Knoll，1983).

由于卫星观测的投影数据不完备，不适合采取直接依据投影定理导出的富里叶积分反演算法，需采用离散的有限级数展开法，在处理问题的最开始就进行离散化(Censor，1983).将反演区域划分为三维网格即体像素，式(4a)离散化，省略左端的上标HAE，代之以观测视线的序号i，得到

其中i、j分别是视线序号和网格序号，M为视线数目，N为网格数目，s_ij是第i条视线在第j个网格内的截距.

上述离散化ENA微分通量方程组可表示为以下矩阵形式：

其中

称为投影矩阵，它依赖于反演区域网格与观测视线族之间的几何关系.于是，反演问题转化为对以J_ion为未知数的代数方程组(6)的求解.本文利用代数重建法(Algebraic Reconstruction Technique)求解方程组，以得到反演区域离子通量的三维分布.ART是一种迭代算法，在每一次迭代过程中计算观测值与上一轮迭代结果之差，然后根据这个差值对反演迭代结果进行反复修正，最终得到一个稳定的收敛解.ART迭代公式如下：

式中k表示第k次迭代，μ_k为松弛因子，取值在0到2之间，对含测量误差的数据，适当选取松弛因子的值可以改善重建图像的质量和迭代效率.一般情况下，迭代过程中松弛因子可以取一个定值1.由于ART算法并不能保证迭代结果为正，因而从实际的物理意义考虑，在ART算法中需要加上离子通量的非负约束.

随着时间的推移，星载探测仪器性能会有所降低，从而探测结果可能存在因仪器而异的偏差，造成多卫星测量的不一致性，这将导致相互矛盾的ENA测量代数方程组，使得反演误差增大，甚至迭代反演不收敛.为了解决这个问题，本文引入差分体像素 CT方法，此方法的思想与电离层/等离子体层无线电射线CT中为消除多站点接收不同GPS信号所得TEC的不一致性而使用的差分TEC反演方法(Xiao et al.,2012)是相似的.我们假定这个偏差是加性误差，并假定同一个探测仪器的偏差为一个常数，在这样的简化假设下，(5)式可表示为

式中ε^q为卫星q的仪器偏差.对于某一卫星，选择一条合适的观测视线m作为参照，用沿任意第l条视线测量到的ENA通量，减去沿参考视线m的ENA通量，得到第l条视线差分ENA通量，它不再含有这个偏差，从而消除了任意射线l上的仪器偏差，可用公式表示为

令

δJ_ena^l=J_ena^l－J_ena^m，δG_lj=G_lj－G_mj，

由(9)式得

此为差分体像素 CT所依据的公式，其中用来进行反演的ENA通量变为差分通量，投影矩阵亦作相应差分转换.

对于本文所用的磁层ENA成像卫星观测区域，ENA离子源区主要包括环电流离子源区和极光带沉降离子源区.我们假定，环电流离子源位于偶极地磁场中2 ＜L＜7.5，覆盖0-24所有磁地方时，磁纬介于±55°之间的区域，这里我们忽略更远处磁尾等离子体片离子源区，其影响将另文讨论.对于极光带离子源区，我们考虑外逸层底高度以上600~9000 km，磁纬55°～78°区域.南半球极区不能被卫星探测到，因此不包括在反演区域内.

考虑到环电流粒子被束缚在磁壳上，并依赖于磁纬和磁地方时，我们采用地磁偶极子场模型，将上述反演区域依L值、偶极磁纬λ、以及磁地方时MLT划分为三维网格(即体像素)，采用太阳磁坐标系(Solar Magnetic Coordinates)，径向地心距离r(以地球半径R_E为单位)与L值和偶极磁纬λ之间由下式联系：

根据卫星锥束扫描的视线分辨率，以及视线与反演区域的相对位置，合理确定体像素的大小.鉴于本文考虑的卫星位于北半球极区上空，角度分辨率4°×4°，锥束扫描视线在北半球较密，在南半球较疏，因此环电流区的纬度网格采取不等间隔划分.低纬赤道区域与较高纬度区不同，南北两半球的较高纬度区又各不相同.在北纬30°以上较高纬度区，间隔为12.5°，而在南半球30°以上较高纬度区，间隔则为20°，南北纬30°以下赤道低纬区域间隔为15°. 又鉴于纬度越高，相邻L壳之间的空间距离越小，对于北纬30°以上环电流区，L值间隔取为1；对于 其他纬度的环电流区，L值间隔为0.5.所有反演区 域的磁地方时MLT的网格间隔均为1小时(经度间隔15°).

对于磁纬55°~78°的北半球极区，不再对纬度进行划分，磁地方时0-24 h划分间隔与环电流区相同，为1小时，高度以5000 km为界划分为高低两个区域.

这样划分网格的结果，反演区域体像素总数为1608个，略少于有效视线(即观测方程)数；对于所关注的环电流反演区，两颗卫星所处的位置比较有利时，每个体像素内至少有数条，多至数十条不同方向的视线穿过.

环电流离子通量严格来说是各向异性的，即随投掷角而变，这使CT反演变得复杂；另一方面，卫星当地测量表明，对于本文所考虑的几keV至几十keV的中等能量离子，离子通量投掷角分布的各向异性相对较弱(Fok et al.,1996)，而且在扰动期间特别是大磁暴期间，离子分布各向异性有弱化的趋势(Chen et al.,1998)，因此本文采用离子通量投掷角分布为各向同性的假设，这也是以往由ENA反演环电流的经典研究采用的假设(Roelof，1987；C:son Brandt et al.,2001，2002).各向异性投掷角分布情况下的处理将另文讨论.

地球外逸层内中性氢通常占主要部分，电荷交换主要发生在高能离子和氢原子之间.对于外层大气氢原子密度分布，本文参考Rairden等(1986)和C:son Brandt等(2002)，采用一个修正的Chamberlain模型，如下式：

式中n_H表示氢原子数密度，单位cm^-3 ；r是地心距离，以地球半径为单位；φ是地方时角，正午为零；θ是SM坐标系下从z轴出发的极角.其中

H的单位为地球半径，此量给出氢原子密度分布的昼夜不对称性.

环电流区域能量离子主要有H⁺，O⁺，其他成 分(如N⁺，He⁺，He⁺⁺等)含量很低，均可以忽略.环电流区域高能H⁺，O⁺与地冕中性H发生电荷交换产生ENA，其电荷交换截面随能量E变化，如式(15)和图 1所示(Lindsay and Stebbings，2005)：

图 1

Fig. 1

图 1 H+和O+与H的电荷交换作用截面 Fig. 1 Cross sections for charge transfer of ions H+

如式(4b)所示，忽略次要离子成分后，电荷交换截面加权的能量离子通量包含氢离子和氧离子两种成分，精确的反演应先将氢ENA和氧ENA加以分离，然后分别反演得到氢离子和氧离子的分布.中等能量不同成分ENA的分离需要特殊的技术来实现，我们对这一技术的有关工作将另文介绍.本文通过使用一个氧和氢离子通量强度比率 δ=J_O/J_H，将(4a)式简化为

取δ=0，即假定环电流能量离子全部是氢离子，本文在数值模拟反演中取δ=0；在实测数据反演中，针对本文所考虑的大磁暴事件，采用δ=1/4.

在本文这一节，采用环电流离子分布理论模型和实际卫星观测几何构形，对差分ENA体像素CT 反演进行数值模拟，并就存在仪器偏差与不存在偏差情况，对差分与非差分CT反演效果进行比较，以检验差分CT方法的可靠性和优势.

图 2

Fig. 2

图 2 数值模拟过程框图 Fig. 2 Block diagram of numerical simulation

数值模拟的框图如图 2所示.

本文模拟采用Shen和Liu(2002)的环电流离子通量模型，离子通量J_ion随能量E、投掷角α、L值、磁纬λ和地方时角φ的变化如下式所示(Lu et al.,2008)：

其中κ是kappa指数(Christon et al.,1991)，κ=5.5，e=(1+1/κ)^κ+1≈2.962 ；L₀=7.3是环电流离子注入区域外边界，E₀是注入离子特征能量，对于本文模拟所考虑的氢离子，此特征能量约为7 keV，J₀是离子注入最大微分通量，取J₀=2.0×10⁵ (s·sr·cm²·keV)^-1(Shen and Liu，2002)；因子h(φ)反映随地方时的分布：

φ_s是环电流离子通量最大值所在的地方时角，本文模拟中取φ_s=180°，即环电流离子通量在午夜最大，ξ反映地方时分布不对称性的程度，取为0.75.我们 假定通量只依赖于L值而不随纬度变化；还假定通量随投掷角分布为各向同性.取α_eq=90°，(11)式化简为

考虑能量E=20 keV的环电流离子，L取值范围2~7.5，磁纬 -55°＜λ＜55°，地方时角φ取值0°~360°.

我们选取2012年7月15日磁暴主相期间16 ∶ 45-17 ∶ 00UT，美国TWINS卫星(McComas et al.,2009)的位置和姿态及相应扫描视线用于本文的模拟.表 1给出SM坐标系中两颗卫星在上述时间段大约中间时刻的位置，可见两颗星都在北半球较高的中纬区，距地心5.3 R_E以远，均处于本文考虑的反演区域以外；所处磁地方时相差约14个小时，比较接近同一个磁子午面；两颗卫星进行扇形锥束扫描的中轴线接近正交，有比较大的视角跨度.

表 1 (Table 1)

表 1 TWINS两颗卫星的位置(2012-07-15，16 ∶ 45-17 ∶ 00 UT) Table 1 Positions of TWINS satellites(2012-07-15，16 ∶ 45-17 ∶ 00 UT) r(RE) λ(Mlat) MLT S1 5.48 54.20°N 17 ∶ 02 S2 5.33 42.93°N 02 ∶ 50

表 1 TWINS两颗卫星的位置(2012-07-15，16 ∶ 45-17 ∶ 00 UT) Table 1 Positions of TWINS satellites(2012-07-15，16 ∶ 45-17 ∶ 00 UT)

图 3a和图 3b分别给出在所考虑的时间段的中点左右TWINS两颗卫星在磁子午面和赤道面内的观测视线覆盖示意图.图 3a中，绿色和蓝色射线分别代表S1、S2的探测视线，红色星号是探测视线与反演网格边界的的交点.图 3b是两颗卫星探测视线在纬度-15°~15°网格内的截距在赤道面上的投影，红色为卫星S1，蓝色为卫星S2.

图 3

Fig. 3

图 3 TWINS视线覆盖示意图.(a)子午面；(b)赤道面 Fig. 3 Illustration of projection of satellite line-of-sight on(a)Meridian plane and (b)Equatorial plane

可以看出，由于两颗卫星所在位置相隔较远，在几乎相对的地方时扇区，因此覆盖区域的射线交叉情况较好，有比较大的角度覆盖，有利于CT反演.在昏侧半球的环电流区域视线覆盖情况好于晨侧半球.由于两颗卫星均处于北半球，所以南半球高纬区域的视线覆盖不理想.

图 4给出在上文所述条件下，对TWINS两颗卫星观测分别进行模拟计算得到的ENA成像.在模拟ENA成像中，能量取E=20 keV，因此由(15)式得到电荷交换截面σ_H=5×10^-16 cm².图像由鱼眼视角给出，图中的颜色值代表ENA微分通量大小，两幅图使用相同的色标.图像中间的白色实线圆圈代表地球所在位置，红色和黄色磁力线分别代表正午和黄昏.图像外圈所标示的角度代表仪器探测的方位角β，取值-90°~270°，径向标示的角度代表仪器仰角η，中心仰角为90°.

图 4

Fig. 4

图 4 鱼眼视角下TWINS两颗卫星模拟得到的ENA图像.红色磁力线为正午，黄色为黄昏 Fig. 4 Simulated ENA images in view of ‘fish-eye’ from the two satellites of TWINS，in which the dipole magnetic field lines are shown，with red line for those at noon and yellow at dusk

利用本文第2节所述的普通CT和差分CT方法，分别对3.3.1节正演模拟得到的ENA微分通量进行反演，ART迭代初始分布取为0，使用非负约束.

图 5给出普通ENA CT重建得到的反演区域离子微分通量分布(中间一列)和差分CT重建得到的离子分布(右边一列)；为了比较，图中左边一列给出了正演模拟所采用的离子微分通量分布模型.普通CT和差分CT反演迭代次数相同.图中上面一行是赤道区(-15°＜λ＜15°)离子微分通量随L(2~7.5)和地方时的分布，左侧为午夜，右侧为正午；下面一行是子午面00 ∶ 00-12 ∶ 00MLT内(左侧为午夜)离子微分通量随L(2~7.5)和纬度(-55°＜λ＜55°)的分布.所有图像均使用相同的色标.

图 5

Fig. 5

图 5 (a，b)模拟的离子微分通量分布模型;(c，d)利用普通CT反演得到的离子微分通量分布；(e，f)利用差分CT反演结果.上面一行是赤道面内分布，下面一行是00-12 MLT子午面内分布 Fig. 5 Comparison of the ion flux distribution from(a，b)model，(c，d)classic CT retrieval and (e，f)differential CT retrieval.The top raw is for equatorial plane，the bottom raw is for meridian plane. Red parts of the earth indicate the dayside

由图 5所示结果可见，在没有加入仪器偏差的情况下，差分CT和普通CT都能很好地重建离子微分通量分布.其中赤道面的反演结果(图 5c和图 5e)在04 ∶ 00-08 ∶ 00 MLT扇区的外圈(L大于3)有明显的误差.结合图 3b可以发现，这个区域只能被一颗卫星探测到，从而没有足够的交叉视线，与CT反演需要的条件相差较远.子午面的反演结果中(图 5d和图 5f)，南半球区域也存在观测视线覆盖不足情况，因此这些区域也存在较大误差.迭代终止时两种反演方法的总体相对均方根误差都降到0.14左右(如图 6).可以说，没有仪器偏差的情况下，两种反演方法都可以有效地重建待反演区域的离子微分通量分布.

图 6

Fig. 6

图 6 无仪器偏差情况下反演的相对均方根误差.(a)普通CT反演；(b)差分CT反演 Fig. 6 Relative root-mean-square error for reconstruction without instrument bias by (a)classic CT method and (b)differential CT method

随着卫星运行时间增加，探测仪器的性能会有所改变，测量参数出现偏移，比如仪器响应因子的增大或减小，且难以及时进行在轨标定，从而探测结果可能存在仪器偏差，特别是两颗卫星偏差不一致.这种偏差和不一致性造成反演涉及的代数方程组出现相互矛盾的方程，导致反演迭代误差增大甚至迭代不收敛.我们对仪器偏差作一简单的加性误差模拟，在3.3.1节得到的两颗卫星的模拟ENA分布(图 4)中，分别加入以下偏差值：在TWINS-1的模拟ENA中减去其均值的十分之一；在TWINS-2中加上其均值的五分之一.使用相同的迭代条件，分别用普通CT和差分CT两种方法进行重建，得到的反演结果如图 7所示.图中左边一列仍为离子分布模型，中间一列为利用普通CT方法重建得到的离子微分通量分布，右边一列为利用差分CT方法得到的通量分布.

图 7

Fig. 7

图 7 对于具有仪器偏差的模拟ENA图像，反演得到的离子微分通量分布. (a，b)模拟的离子微分通量分布模型 ;(c，d)利用普通CT反演得到的离子微分通量分布；(e，f)利用差分CT反演结果.上面一行是赤道面内分布，下面一行是00-12 MLT子午面内分布 Fig. 7 The same as Fig. 5 but the instrument biases are applied in ENA images which are further used to reconstruct the RC ion flux distribution

由图 7可以发现，加入仪器偏差后，普通CT(图 7(c，d))反演得到的结果杂乱无章，无法重建反演区域离子微分通量分布，而差分CT(图 7(e，f))可以有效地消除仪器偏差的影响，得到接近原始模 型(图 7(a，b))的结果，其结果也与图 5相近.通过比较两者的相对均方根误差随着迭代次数的变化(图 8)，可以发现在加入仪器偏差后，普通CT迭代反演过程中，均方根误差越来越大，迭代无法收敛；而差分CT反演结果在10次迭代内快速收敛，其相对均方根误差最后降到0.14左右.

图 8

Fig. 8

图 8 加入仪器偏差情况下反演的相对均方根误差.(a)普通CT反演；(b)差分CT反演 Fig. 8 Relative root-mean-square error of inversion by(a)classic CT method and (b)differential CT method from images with instrument bias included

由此可见，在两颗卫星ENA图像具有不一致的固定偏差时，普通CT反演难以重建反演区域离子微分通量分布，而差分CT方法能消除这个偏差带来的影响，很好的重建得到反演区域离子微分通量分布，表明差分CT方法处理此类问题的优越性.

在这一节，我们由TWINS卫星实测ENA数据，利用差分CT方法反演环电流离子微分通量强度的三维分布与能谱，并与THEMSIS卫星的离子能谱当地测量数据进行对比.

选取2012年7月15日磁暴主相期间，TWINS两颗卫星在大约15 min里(16 ∶ 45~17 ∶ 00 UT)约12轮扫描平均的ENA图像，由其重建环电流离子分布.单次扫描可能会因干扰而不完整或者有较大噪音，所以选取多轮扫描的累计平均以提高图像的质量.如果采用每颗卫星的12~13幅图像来进行反演，虽然观测视线的数目明显增多，但每颗卫星观测视线的角度变化很小(总计约2°以内)，对于改进观测视角覆盖范围的作用不大，同时又损失了用来进行反演的ENA图像的质量.二者权衡，采用多次扫描平均为宜.考虑能量为4，8，12，16，20，25，30，50 keV的离子.磁暴期间Dst指数如图 9所示，其中灰色竖线表示进行环电流反演的时间段，处在临近主相极大时.

图 9

Fig. 9

图 9 磁暴期间Dst变化，灰色竖线是本文研究的时间段 Fig. 9 Dst index during the storm on July 14-16，2012.The gray bar is for concerned time period

TWINS卫星的位置如3.2节中的表 1所示；这里取15 min时间段大约中点时刻卫星所在的位置.在本文分析的时间段内，TWINS两颗卫星的地心距离变化约0.2R_E，磁纬变化约2°，磁地方时改变约12 min；此变化相对于本文反演区域及其网格大小来说很小，不足以影响反演结果，因此忽略卫星的位置变化，以TWINS卫星在中间时刻的位置来代表.此间卫星扫描视线对反演区域的覆盖情况如3.2节的图 3所示.所用原始数据是由美国西南研究所提供的经过平滑的ENA微分数通量二维图像数据，取自cdaweb网站.参照McComas等(2012)给出的方法，对图像进行消除背景，得到用来进行反演的来自两颗卫星的一对ENA图像，图 10给出能量为20 keV时的一个示例，在仪器坐标系中给出，径向为仰角η，取值30°~90°；角向为方位角β，取值-90°~270°.在其他能量上做相同处理，得到分别在4，8，12，16，20，25，30，50 keV上的TWINS两颗卫星的用来进行反演的图像.

图 10

Fig. 10

图 10 消除干扰背景后的TWINS ENA图像 Fig. 10 TWINS ENA images with background interferences removed

利用消除背景等预处理后的ENA数据，反演重建环电流离子微分通量的三维分布.如2.4节所述，在反演中假设氧离子和氢离子的通量强度比为常数0.25，离子通量分布为各向同性；依然采用ART迭代算法，初始分布取作零.

图 11a给出在几个不同能量(4，8，16，30 keV)上，赤道区(-15°-15°N)H⁺离子微分通量随L值和MLT的分布.图 11b给出在01-13 MLT子午面内不同能量的离子通量随L值和纬度的分布.图中的空缺部分(白色网格)是由于缺少卫星射线覆盖以及数据不可用；所有各幅图像都使用相同的色标.

图 11

Fig. 11

图 11 图 11 差分CT重建的环电流区离子分布.(a)赤道区H+微分通量强度随L值和MLT的分布；(b)01 ∶ 00 MLT子午面内H+离子微分通量强度随L值和纬度的分布.图中左侧为夜晚 Fig. 11 Reconstructed RC ion distribution at different energy levels by differential voxel CT method.(a)The distribution of H+ differential flux versus L-value and MLT in equatorial plane；(b)The distribution of H+ differential flux versus L-value and Mlat in meridian plane. Red parts of the earth indicate the dayside

从图 11可以看出，反演所得离子通量呈现出显著不对称的部分环电流特征，强的离子通量主要分布在磁午夜前后至黎明前数小时L值约3.5~6.5的区域，对于较低的能量(如4，8 keV)，在午夜前和黄昏扇区也有较强的通量；从子午面内的分布可以看出，强的通量主要出现在低纬-赤道区，与环电流的一般特征相符.

在本文研究的磁暴期间，THEMIS-D卫星刚好处于环电流区域内：r=5.31 R_E，MLT=01 ∶ 45，Mlat=15.07°S附近，该卫星携载的ESA(ElectroStatic Analyzer)和SST(Solid State Telescope)两种仪器联合，能够探测6 eV到300 keV的离子能通量谱(其间20~30 keV有一小段能量空缺).图 12给出THEMIS-D在2012年7月15日16 ∶ 00到17 ∶ 00 UT期间ESA和SST测量到的离子能通量谱，横坐标是时间，纵坐标是能量，颜色代表能通量强度.其 中ESA数据的能量范围是6 eV到20 keV，时间分 辨率为97 s；SST数据的能量范围是30 keV到300 keV，时间分辨率为3 s.

图 12

Fig. 12

图 12 THEMIS-D卫星实测离子能谱随时间的变化. 上半部为SST 测量(30 keV～300 keV)，下半部为ESA测量(6 eV~20 keV) Fig. 12 Ion energy spectra measured in situ onboard THEMIS-D satellite during 16 ∶ 00-17 ∶ 00 UT on July 15，2012. The top one is measured by SST(at 30 keV~300 keV) and the bottom is by ESA(6 eV~2 0 keV)

我们将THEMIS-D实地探测结果与利用ENA 图像反演得到的离子能通量谱作比较，以验证反演结果的可靠程度.首先将 16 ∶ 45-17 ∶ 00 UT期间反演得到的在5.0 ＜L＜6.5，-15°＜Mlat＜0°，00 ∶ 00＜MLT＜03 ∶ 00区域内的离子微分数通量求空间平均，将不同能量上所得平均数通量对相应的能量段求积，得到能通量谱.图 13中较粗的红色点实线是由16 ∶ 45-17 ∶ 00 UT期间TWINS ENA数据，利用差分ENA体像素 CT方法反演得到的离子能通量谱；其他颜色实线是THEMIS-D卫星在16 ∶ 45-17 ∶ 00 UT期间，在相继的不同时刻测量得到的离子能通量谱.可以看出，二者的变化趋势一致，彼此符合得很好，说明本文的差分ENA体像素CT反演方法是可靠的.

图 13

Fig. 13

图 13 TWINS ENA反演得到的离子能通量谱(红色粗点线)与THEMIS-D实测结果的比较 Fig. 13 Comparison of the ion energy spectra between retrieved(red line)from TWINS ENA data and measured by THEMIS-D

本文发展了由卫星ENA二维图像重建磁层环电流能量离子通量三维分布的CT反演方法，综合考虑电荷交换物理机制与卫星扫描视线几何特征，比较合理地确定了反演区域及其体像素划分方案；针对卫星ENA成像可能存在的仪器偏差问题进一步发展了差分ENA 体像素CT方法，并通过正/反演数值模拟，证明差分CT可以比较有效地降低仪器整体偏差造成的大的反演误差，避免迭代反演发散，比普通CT具有更宽的适用性和优越性.将此差分ENA体像素CT方法应用于2012年7月15日磁暴事件中TWINS两颗卫星的ENA数据，反演环电流离子微分通量的三维分布，得到了合理的结果，这是在不受绝热不变条件限制下反演得到的环电流区离子三维分布，至今国内外已有反演大都使用绝热不变假设重建赤道面内二维分布(C:son Brandt et al.,2002; DeMajistre et al.,2004; Lu et al.,2010; Perez et al.,2012).所得离子分布呈现显著昼夜不对称，夜侧远大于白天，这与磁暴主相期间的环电流为部分环电流的特征相呼应；对于所分析的磁暴事件，较强的离子通量主要出现在低纬-赤道区磁午夜前后至黎明前数小时，这与传统认为的强的环电流离子通量只出现在黄昏至午夜前不同，但与IMAGE卫星揭示的现象(Fok et al.,2003)相符；而且我们的反演结果还显示，环电流强度的地方时分布依赖于离子的能量，在此次事件中较低能量(10 keV以下)的离子在午夜前和黄昏扇区也有较强的通量.将反演所得磁午夜后赤道区环电流离子通量能谱，与THEMIS卫星当地测量得到的此区域同时的能谱作比较，两者符合得很好，证明本文发展的差分ENA体像素CT是由多卫星ENA二维图像重建暴时环电流离子分布的有效方法.

应该指出，由于磁层成像卫星探测到的ENA的离子源并非单一的被地磁场捕获的环电流离子，还有近地等离子体片中的离子，以及到达较低高度大气层的沉降离子(来自环电流和等离子体片)；与高能离子发生电荷交换的背景冷中性气体既有光学薄的地冕氢，也有浓密的光学厚大气层中的氧原子，后者涉及多重碰撞、往复的不同类型电荷交换与能量损失过程.因而下一步工作还需要更加合理细致地区分和确定反演区域以及相应区域电荷交换过程的物理性质与理论公式，以便对沉降离子源引起的低高度光学厚介质中ENA发射做出合理修正，对近地等离子体片ENA发射源的影响作出评估.高空磁场与地冕氢分布模型对反演具有一定的影响，采用更接近实际的非偶极子磁场模型，以及利用TWINS卫星的Ly-α成像反演得到的地冕氢分布，将会有助于改进反演结果.另外，尽管TWINS两颗卫星从不同角度扫描给出的两幅ENA二维图像(即层析技术的所谓“投影”)，提供了比单颗卫星更为丰富的离子源分布信息，有可能进行CT 反演，但是投影数据还是很不完备的，还需要更为充分的合理先验知识的约束以求得更为可靠的物理上合理的解，反演算法还有进一步优化的空间；在条件允许情况下，发展ENA卫星编队CT遥测，应是从根本上解决反演数据不完备的最好办法.

致谢 感谢美国TWINS卫星计划团队与PIs通过互联网提供TWINS ENA数据；感谢J. Goldstein博士对本文工作的支持与有益建议.感谢美国国家宇航局通过SPDF提供THEMIS卫星数据.