顶部电离层和等离子体层电子密度分布——基于GRACE星载GPS信标测量的CT反演

刘裔文，徐继生＊，徐 良，尹 凡，马淑英，H.Lühr

1 武汉大学电子信息学院，武汉 430072

2 Germany Research Center for Geoscience，Telegrafenberg，14473Potsdam，Germany

1 引 言

1986年，美国伊利诺大学波传播实验室首次提出电离层层析（简称 CT）的构想［1－2］.随后，国际上多个研究小组先后开展电离层CT探测与研究，取得了一系列有意义的成果，表明电离层CT是一种强有力的电离层探测技术［3－6］.叶公节等通过理论分析讨论了电离层CT的潜力及其局限性［7］.

1994年，24颗GPS卫星星座布设完成，利用遍布全球的GPS接收台网，开展地基电离层CT探测成为可能.20世纪90年代中后期，已有一些研究者提出利用GPS台网测量数据开展电离层CT成像［8－9］.为将电离层CT从准静态二维推广到时变三维，还有一些小组先后研究了时变三维电离层CT的原理和算法，并给出了初步的研究结果［10－13］.不过，仅用地基GPS接收台网的观测数据反演电离层电子密度分布存在严重的局限性，一个主要的问题是数据严重不完备，特别是缺少近水平向的射线，导致重建的电子密度分布，其垂直分辨率较低.为改善反演质量，提高垂直分辨率，一些研究者先后把GPS掩星数据和电离层垂测数据加入到电离层CT反演中，以补充近水平向射线的缺失［14－15］.

对于地基CT反演，沿传播路径的电子密度质心高度通常在350km和450km之间，射线在1000km以上的部分对总电子含量（简称TEC）的贡献很小，所以地基CT探测区域的高度范围通常取在100～1000km.利用电离层电子密度分布模型，通过正演可估计出，对于轨道高度在500km上下的低轨卫星，沿传播路径的电子密度质心高度通常在1000km上下，因此，利用低轨卫星载GPS双频信标测量解算的TEC数据与合适的CT反演算法，得到低轨卫星高度以上顶部电离层和等离子体层的电子密度分布的信息，是一条有效的技术途径.最近，Yizengaw等曾利用FedSat卫星GPS信标观测数据，开展天基CT反演，重建顶部电离层和等离子体层电子密度分布，得到了一些有意义的结果［16－17］.

本文将利用两颗跟飞的GRACE卫星载GPS双频信标测量解算的TEC数据，借助差分相对TEC层析算法，开展全球范围的顶部电离层和等离子体层（450～5000km）天基CT成像.对不同地磁活动条件下的CT反演结果表明，等离子体层电子密度随纬度的分布是不均匀的；在低纬赤道带，经常出现近似垂直于磁力线的电子密度柱状增强结构，从顶部电离层向上延伸直到等离子体层.在下一节，将首先介绍基于差分相对TEC的层析算法，然后介绍所使用的数据，在第4节将给出磁静日和不同强度磁扰动条件下的CT反演结果，最后给出总结和简要的讨论.

2 原理和算法

CT技术的基础是所谓“投影”（projections），即介质对波作用的积分效应，它可以表述为，已知待测量穿过探测区域的线积分，重建待测参量的空间分布.在天基电离层－等离子体层CT反演问题中，线积分是沿GPS卫星信标发射机至低轨卫星载信标接收机路径的TEC，而待测量为探测区域内的电子密度分布.观测方程为

式中Si表示GPS卫星至低轨卫星载信标接收机的第i条射线传播路径，NTi表示沿Si的TEC，Ne（r）表示沿路径的电子密度分布，r是位置矢量.

借助计算机数值求解，需要将方程（1）离散化，把探测区域沿垂直和子午方向划分为N个网格，观测方程变为

式中di表示沿第i条射线测量的TEC，aij表示第i条射线穿过第j个网格的截距，bj是基函数，Nej表示第j个网格中的电子密度.

对于不同的接收机，存在不同的硬件偏差和整周模糊度，导致在同一时刻沿同一路径测量的TEC的量值不一致.精确消除这种硬件造成的不一致性很困难.一些作者先后提出，用差分相对TEC进行反演，可以回避硬件造成的不一致性这一难题［18－19］.计算相对TEC的差要求对每一段连续的数据进行，求第k条射线和第1条射线之差后，由（2）式可以得到

若有M个形如（3）式的差分表达式，则离散化的观测方程可写成矩阵形式，

式中D是一个列向量，表示沿M条路径的差分相对TEC，A表示M×N维系数矩阵，X也是一个列向量，表示待测的各网格内的电子密度.

于是，CT问题转化为线性代数方程组（4）的求解问题.为求解（4）式，本文利用代数重建算法（ART）来重建电子密度二维分布.ART是一种迭代算法，在每一次迭代过程中计算观测值与上一轮迭代结果之差，然后根据这个差值对待测量反复修正，最终得到一个稳定的收敛解.在数据不完备的条件下，选择迭代过程所需要的合适的初始分布非常重要.

3 数 据

2002年，德国航空中心与美国国家航空航天局合作，发射升空了GRACE卫星.GRACE卫星星座含两颗跟飞的小卫星，即GRACE－A和GRACE－B，两颗卫星在同一轨道面内，前后相距约200km，都为近极圆轨道，轨道倾角为89°，轨道设计高度为500 km，运行周期约1.57h.两颗GRACE卫星都载有GPS双频信标接收机.本文所使用的GRACE卫星GPS信标观测数据采样间隔为10s，在相邻两次采样时间内，卫星飞行了约70km.为了获得较好的射线覆盖，同时又能近似作为一个二维问题进行反演，本文定义满足以下条件的射线为有效射线：从GRACE卫星轨道高度至5000km高度，射线与GRACE卫星轨道面的夹角都在±15°之内.图1给出GRACE－GPS有效射线覆盖的一个示例.

从图1中可以看出，两个半球绝大部分区域都有很好的射线覆盖，且由于每颗GPS卫星同时被GRACE－A和GRACE－B接收，因此覆盖区域的射线交叉情况也较好，有利于CT反演.不过，由于GPS卫星轨道面倾角只有55°，在高纬地区较高高度上，没有射线覆盖.

图1 GRACE－GPS射线覆盖示意图Fig.1 A sketch map of the coverage for the GRACE－GPS effective ray paths

图2a和图2b分别给出了地基GPS台站观测到的TEC和GRACE卫星观测到的TEC及相应的接收机对GPS卫星仰角随时间的变化.从图2a可以看出，当仰角较低时，地基GPS台站观测到的相对TEC曲线出现周跳和明显的快速随机起伏，这主要是多路径效应所引起.这种现象经常出现，因此，在地基CT中，只有较高仰角的射线可用，缺少接近水平向的射线，数据严重不完备.这导致反演得到的电子密度分布垂直分辨率较低.从图2b可以看出，GRACE卫星观测到的相对TEC，即使在仰角低至0°时依然很平滑，没有周跳和明显抖动.所以，对于天基CT，仰角接近0°的射线依然可以使用，有利于改善反演结果的垂直分辨率，这是天基CT与地基CT相比的优势之一.此外，地表大部分区域为海洋所覆盖，在海洋覆盖的区域，无法获取地基测量数据进行CT反演，对于天基CT，则不受此限制，这是天基CT与地基CT相比的另一优势.

4 反演结果

图2 地基观测台站测量的斜TEC（a）和GRACE－A卫星载接收机测量的斜TEC（b）比较.实线代表斜TEC，虚线代表仰角Fig.2 A comparison of slant TEC measured by ground－based GPS receiver （a）and slant TEC from GRACE－A borne receiver（b）.The solid and the dash lines represent variations of TEC and elevation with universal time，respectively

本文设定CT反演区域的高度范围从450km至5000km，纬度范围从70°S至70°N.在南北方向均匀划分网格，网格尺度为1°（～110km）.在垂直方向，取不等间距网格，1000km以下高度，网格尺度为50km；1000至2000km高度，网格尺度为100km；高度在2000km以上，网格尺度为200km.

图3是利用上述算法得到的一个反演实例，以及反演结果与CHAMP卫星实测电子密度和DMSPF13卫星实测离子密度的比较.

图3a中在约800km以下，赤道电离异常结构清晰可见.从图3b可以看出，CHAMP卫星载探针实地测量得到的电子密度与本文反演得到的电子密度随纬度分布的形状基本一致.图3c是850km高度上CT反演的电子密度与同一高度上DMSP－F13卫星实地测量的离子密度随纬度变化的比较，表明两者无论量值还是随纬度的变化趋势均基本一致.图3表明，本文CT反演结果是可靠的.

4.1 磁静日的CT反演结果

图3 2003年4月4日CT反演得到的电子密度剖面（a）；以及反演结果分别与CHAMP卫星实测结果（b）和DMSP－F13卫星实测结果（c）对比图Fig.3 Reconstructed 2－D electron density profile on April 04，2003（a）；and comparison of reconstructed electron density at 500km and electron density measured by CHAMP at about 400km （b）；and comparison of reconstructed electron density at 850km and plasma density measured by DMSP－F13at about 850km （c）

图4 2004年6月27日CT反演得到的全球范围顶部电离层－等离子体层电子密度二维剖面Fig.4 Reconstructed 2Dplasma density profile of the topside ionosphere and plasmasphere on June 27，2004

图4为2004年6月27日CT反演得到的全球电子密度二维剖面，图中白色弧线给出了偶极磁场的磁力线分布，箭头所指表示磁赤道的位置，上半图对应15UT（10LT）前后，磁赤道所在地理纬度约为11°S，下半图对应16UT（22LT）前后，磁赤道所在地理纬度约为10°N.在得到图4及以下的反演结果时，迭代所需的初始电子密度分布均由NeQuick模型计算得到.

2004年6月27日地磁活动非常平静，全天∑Kp＝6，平均Ap指数为3.25.从图4中可以看出，在磁静日，在大约1000km以下的顶部电离层，电子密度随高度升高下降很快，从450km至1000km，电子密度下降了一个量级左右；昼夜子午面内电子密度随纬度的分布均呈现出南北不对称，夜间半球南北不对称性更为显著，北半球电子密度值较高；此外，在等离子体层高度上，电子密度值较低，昼夜子午面内电子密度随纬度变化都呈现一定程度的不均匀性，但没有明显的扰动或大尺度的异常结构出现.

4.2 磁扰日的CT反演结果

图5a给出了2003年4月10日07∶10UT（19 LT）前后的CT反演结果.图5b是NeQuick模型得到的初值分布.图中白色弧线给出了偶极磁场的磁力线分布，磁赤道所在地理纬度约为5°N，箭头所指表示磁赤道的位置.这一天存在弱的地磁扰动，Kp指数最高为4，全天∑Kp＝29＋，Ap指数最高为32，平均约28.

由图5可以看出，175°E子午面正处于日落时段（19LT前后），此时电子密度量级依然处于较高水平，在顶部电离层和等离子体层的较低高度上，电子密度分布基本上关于磁赤道对称.图中一个值得注意的现象是，在等离子体层中，电子密度随纬度的分布是不均匀的，特别是在磁赤道及邻近区域上空，电子密度随高度下降的速度比其他纬度更慢，这就使得磁赤道附近顶部电离层存在一个电子密度显著增强的区域，向上一直延伸到等离子体层.在约4000km以上，磁赤道以南5°左右存在一个相对孤立的电离增强云团.

为了定量地分析图5a所示的电子密度分布的特征，我们选取了若干高度，给出在这些高度上反演得到的电子密度随纬度的变化，如图6所示.

从图6中可以看出：（1）在所有高度上，低纬地区电子密度都明显高于中、高纬地区；（2）在450km高度峰值电子密度最大，接近9×1011m－3，至1000km高度，峰值电子密度降低至1011m－3左右，2000km高度，峰值电子密度降为1000km处的50%左右，3000km高度，峰值电子密度降为2000km的40%左右.若电子密度分布处于扩散平衡状态，设n＝noexp（－h／H），利用电子密度随高度的衰减速率，可粗略估计出，800km至1000km，平均等离子体标高H约345km；2003年4月10日06∶00UT前后，DMSP F13号卫星的轨道面位于180°E附近（约18∶00LT），非常接近图5CT反演所在的子午面，利用DMSP F13号卫星实测离子成分及离子和电子温度数据以及标高的定义，可以近似估算出，DMSP卫星高度（约850km）的等离子体标高约为340km，与依据CT反演结果估计的等离子体标高基本一致；（3）电子密度随纬度变化呈现不均匀分布，在某些纬度，存在电子密度耗尽，而在另一些纬度，出现电子密度增强，且在等离子体层中，随着高度增加，电子密度随纬度的相对起伏更加明显；（4）在低纬赤道区，3500km至5000km高度上，电子密度随高度出现异常变化，其值不是随高度增高而减小，而是出现弱的增强；（5）电子密度的极大值在1000km高度位于磁赤道偏北上空，随高度增高，电子密度增强区的中心跨过磁赤道，逐渐向南移动，至5000km高度，电子密度增强区的中心偏移至磁赤道以南约7°左右.

图6 2003年4月10日175°E子午面CT反演得到的不同高度上电子密度随纬度的变化.图中竖虚线代表磁赤道Fig.6 Electron density variations with latitude in different altitudes obtained by CT reconstruction on Apr.10，2003，175°E. The vertical dash lines indicate the geomagnetic equator

图5 （a）2003年4月10日，CT反演得到的顶部电离层－等离子体层电子密度分布；（b）2003年4月10日，NeQuick模型得到的初始电子密度分布Fig.5 （a）Reconstructed 2－D electron density profile of the topside ionosphere and plasmasphere on April 10，2003；and （b）initial electron density profile derived from NeQuick model on April 10，2003

2003年10月29日至11月1日发生了一次多主相型的超强磁暴，磁暴急始发生在29日06∶12UT，第3个主相出现在10月30日，SYM-H指数在当天23∶36UT达到极小值－403nT.图7a为2003年10月30日CT反演结果.反演结果对应的地方时在16LT前后，即该超强磁暴第3个主相环电流开始增强的阶段.2004年7月22至10月28日发生了一个多主相型的强磁暴，磁暴急始发生在22日10∶39UT，第2个主相出现在7月25日，SYM-H指数在18∶30UT达到极小值－168nT.图7b为2004年7月25日CT反演结果.该图给出的电子密度分布对应的地方时在20LT前后，即磁暴第2个主相环电流达到极大期间.从图7a和图7b可以看出，与图5类似，在低纬地区也出现近似垂直于磁力线分布的电子密度增强结构.图7a中的电子密度增强结构出现在磁赤道以南，从顶部电离层开始一直延伸到5000km，形成一个较强的柱状电离增强结构，该结构随着高度升高逐渐向南偏移，至5000km处到达地理纬度20°S（磁纬约11°S）左右，增强结构形态下粗上细，随高度增高，其覆盖的纬度范围逐渐变小，在等离子体层中，增强结构的中心最大电子密度比背景电子密度高2至5倍.Yizengaw等［16］曾利用FedSat卫星GPS信标测量解算的TEC数据进行了天基CT反演，结果表明，在10月30日21∶50UT前后，160°E子午面（07∶50LT），南半球中纬度地区也曾有与图7a形态类似的近似垂直于磁力线的电子密度增强结构出现，这表明等离子体层中，这一类电子密度增强结构可以重现，具有一定的普遍性.图7b表明，在低纬赤道区顶部电离层直至2500km，较大的纬度范围内出现的电子密度增强，也许是超级赤道喷泉效应所致.在等离子体层约3500km以上，地理纬度约15°N至30°N（磁纬5°－20°N）之间，出现一个相对孤立的电离增强云团.在约2500km至3500km之间，低纬赤道区存在一些不连续的弱电离增强.

5 总结与讨论

图7 2003年10月30日（a）和2004年7月25日（b）CT反演的顶部电离层－等离子体层电子密度分布Fig.7 Reconstructed electron density profile of the topside ionosphere and plasmasphere on（a）Oct.30，2003and（b）July.25，2004

本文利用两颗跟飞的GRACE卫星载GPS信标TEC测量数据和基于差分相对TEC的CT反演算法，实现了全球范围顶部电离层和等离子体层的CT成像.研究结果表明，基于低轨卫星的星载GPS双频信标测量解算的TEC数据进行天基CT反演，是获得顶部电离层－等离子体层电子密度二维分布的一种有效的技术手段.

天基CT技术与地基CT或天地基结合的CT技术相比主要有下列不同点.首先，主要由于多路径效应的影响，开展地基CT反演，只有较高仰角的射线可用，缺少接近水平向的射线，而在数百千米高度飞行的低轨道卫星，其星载GPS信标接收机的电磁环境很好，仰角接近0°的射线依然可以使用，有利于改善反演结果的垂直分辨率.其次，在海洋覆盖的区域，无法获取地基测量数据，而天基CT则不受此限制，这有利于在广阔的海洋上空进行CT反演，如本文图5所示，175°E子午面电离层几乎都位于海洋之上，只有通过天基CT反演，才能得到该子午面的电离层－等离子体层电子密度分布.此外，对于地基CT，射线的电离层穿刺点一般位于350至450km高度，有效反演区域仅能覆盖电离层，而对天基CT，穿刺点一般在1000km或更高的高度，可以同时反演电离层－等离子体层电子密度分布.不过，天基CT只能得到轨道高度以上的电离层－等离子体层电子密度分布，不能得到轨道高度以下电离层的电子密度分布，这是它的局限性.

静日的CT反演结果表明，昼夜子午面内电子密度随纬度的分布均呈现出南北不对称，夜间半球南北不对称性更为显著；在等离子体层高度上，电子密度值较低，没有强的大尺度扰动或异常结构出现.在不同强度的磁暴期间，几个CT反演示例表明，在低纬地区都出现近似垂直于磁力线分布的电子密度柱状增强结构和等离子体层中孤立的电离增强云团，电子密度柱状增强结构下宽上窄，随高度增高，其覆盖的纬度范围逐渐变小，在等离子体层中，增强结构的中心最大电子密度比背景电子密度高2至5倍.文献［16］和［20］的CT反演结果，也显示有类似的垂直于磁力线的电子密度增强结构出现，说明这种现象具有一定的普遍性.

本文反演结果中出现的近似垂直于磁力线从顶部电离层向上延伸到等离子体层的电离柱状增强以及等离子体层中局地的电离增强云团，是近年来CT反演发现的一种新的现象.研究表明，磁暴期间，跨极盖电位的快速变化，有利于极区对流电场向中低纬电离层穿透［21］，强的东向即时穿透电场，通过低纬赤道区离子E×B漂移机制，引起较长时间持续的强离子上行，是形成近似垂直于磁力线从顶部电离层向上延伸到等离子体层的电离柱状增强结构的可能原因之一；此外，源于磁层的热粒子向下注入，也是形成等离子体层中局地的电离增强云团的一种可能的过程.如前所述，这种现象多次反复出现，具有一定程度的普遍性，其形成机制是一个有待解决的问题，值得进一步深入研究.

（References）

［1］Austen J R，Franke S J，Liu C H，et al.Application of computerized tomography technique to ionospheric research.／／ Tauriainen A ed.URSI and COSPAR International Beacon Satellite Symposium on Radio Beacon Contribution to the Study of Ionization and Dynamics of the Ionosphere and to Corrections to Geodesy and Technical Workshop，Oulu，Finland.Proc.Part I，25.University of Oulu，ISBN 951－42－2256－3，1986.

［2］Austen J R，Franke S J，Liu C H.Ionospheric imaging using computerized tomography.Radio Sci.，1988，23（3）：299－307.

［3］Raymund T D，Austen J R，Franke S J，et al.Application of computerized tomography to the investigation of ionospheric structures.Radio Sci.，1990，25（5）：771－789.

［4］Pryse S E，Kersley L，Rice D L，et al.Tomographic imaging of the ionospheric mid－latitude trough.Ann.Geophys.，1993，11（2－3）：144－149.

［5］Pryse S E，Mitchell C N，Heaton J A T，et al.Travelling ionospheric disturbances imaged by tomographic techniques.Ann.Geophys.，1995，12（12）：1325－1330.

［6］徐继生，马淑英，阳其罕等.东亚赤道异常区电离层CT诊断－实验及初步结果.地球物理学报，1995，38（5）：553－563.Xu J S，Ma S Y，Yang Q H，et al.Tomographic diagnosis of the ionospheric equatorial anomaly in the region of East－Asiaexperiment and early results.Chinese J.Geophys.（in Chinese），1995，38（5）：553－563.

［7］Yeh K C，Raymund T D.Limitations of ionospheric imaging by tomography.Radio Sci.，1991，26（6）：1361－1380.

［8］Hajj G A，Ibañez－Meier R，Kursinski E R，et al.Imaging the ionosphere with the Global Positioning System.Int.J.Imaging System Tech.，1994，5（2）：174－184.

［9］Kunitsyn V E，Andreeva E S，Razinkov O G.Possibilities of the near－space environment radio tomography.Radio Sci.，1997，32（5）：1953－1963.

［10］Mitchell C N，Spencer P S J.A three－dimensional timedependent algorithm for ionospheric imaging using GPS.Ann.Geophys.，2003，46（4）：687－696.

［11］徐继生，邹玉华.时变三维电离层层析成像重建公式.地球物理学报，2003，46（4）：438－445.Xu J S，Zou Y H.Reconstruction formula of time－dependent 3－D computerized ionospheric tomography. Chinese J.Geophys.（in Chinese），2003，46（4）：438－445.

［12］徐继生，邹玉华，马淑英.GPS地面台网和掩星观测结合的时变三维电离层层析.地球物理学报，2005，48（4）：759－767.Xu J S，Zou Y H， Ma S Y. Time－dependent 3－D computerized ionospheric tomography with ground－based GPS network and occultation observations.Chinese J.Geophys.（in Chinese），2005，48（4）：759－767.

［13］Ma X F，Maruyama T，Ma G，et al.Three－dimensional ionospheric tomography using observation data of GPS ground receivers and ionosonde by neural network.J.Geophys.Res.，2005，110（A5），A05308，doi：10.1029／2004JA010797.

［14］Rius A，Ruffini G，Cucurull L.Improving the vertical resolution of ionospheric tomography with GPS occultations.Geophys.Res.Lett.，1997，24（18）：2291－2294.

［15］Hernández－Pajares M，Juan J M，Sanz J，et al.Global observation of the ionospheric electronic response to solar events using ground and LEO GPS data.J.Geophys.Res.，1998，103（A9）：20789－20796.

［16］Yizengaw E，Dyson P L，Essex E A.A study of the spatial density distribution in the topside ionosphere and plasmasphere using the FedSat GPS receiver.Advances in Space Research，2006，38（11）：2318－2323.

［17］Yizengaw E，Moldwin M B，Dyson P L，et al.First tomographic image of ionospheric outflows.Geophys.Res.Lett.，2006，33（20）：L20102，doi：10.1029／2006GL027698.

［18］Kunitsyn V E，Preobrazhensky N G，Tereshchenko E D.Reconstruction of ionospheric irregularities structure based on the radioprobing data.Dokl.Akad.Nauk SSSR，1989，306：575－579.

［19］Kunitsyn V E，Tereshchenko E D，Andreeva E S，et al.Radio tomography of Global Ionospheric Structure.Preprint Polar Geophys.Inst.，90－10－78.Apatity：Kola Science Center AS USSR，1990：1－30.

［20］Xiao R，Xu J S，Ma S Y，et al.Abnormal distribution of ionospheric electron density during November 2004superstorm by 3DCT reconstructions from IGS and LEO／GPS observations.Sci.China （Ser.E－Tech.Sci.），2012，55（5）：1230－1239.

［21］Xiong W，Xu J S，Wang H，et al.Effect of temporal variation rate of cross polar cap potential on the equatorial ionospheric vertical drift：A statistical study.Sci.China（Ser.E－Tech Sci.），2012，55（5）：1217－1223.