磁暴事件对高精度TEC二维分布的影响

余侯芳,郭文兴,张玉强,冯健,甄卫民

(1. 中国电波传播研究所,山东 青岛 266107;2中兴通讯股份有限公司,广东 深圳 518057;3. 西安电子科技大学 物理与光电工程学院,陕西 西安 710071)

0 引 言

电离层总电子含量(TEC)是表征电离层特性的重要特征参量之一,研究其时空变化规律对卫星通信、导航定位、空间天气预报等领域有着极其重要的理论意义及应用价值[1-2]．

实验探测中,电离层TEC值主要由地面接收站接收的卫星信号获得．在信号由卫星向接收机传播时,由于电离层的色散性质,不同频率信号所引起的伪码时延和载波相位时延也不同,因此可以通过不同频率信号间的差分求得信号路径上的电离层总电子含量[3]．二十世纪七十年代以后,卫星导航定位系统的出现,极大地促进了电离层探测技术的发展,其对电离层TEC的研究有着极其重要的作用及影响．由于全球定位系统(GPS)发展时间较早,卫星星座多,接收站分布广泛,因此基于GPS的电离层研究工作较多[4-5]．

目前,国际上已有许多电离层研究机构实现了基于GPS的TEC地图重构数据产品的发布,包括欧洲定轨中心(CODE)、美国喷气动力实验室(JPL)以及欧洲空间局(ESA)等机构分别完成了全球电离层地图重构算法研究,并于事后提供了基于IGS的全球卫星导航系统(GNSS)观测站数据的全球电离层TEC地图产品——GIM[6-7]．此外,日本、欧洲、美国、澳大利亚等国家和地区都建立精度更高、实时性更好的地区性TEC预报系统[8-11]．

近年来,国内也有很多单位开展了GPS-TEC的研究与观测工作．其中,中科院地质与地球物理研究所(IGGCAS)开发了我国首个覆盖全境的TEC现报系统[6]．毛田等[12]用Kriging方法构建了中纬度区域的电离层TEC地图．雷宵龙等[13-14]利用COSMIC数据对全球电离层TEC地图进行了构建,证明了利用COMSIC掩星资料构建电离层TEC地图的可行性．

随着全球电离层TEC地图的建立,人们对电离层暴的理解与研究越来越深入,例如:Balan等[15-16]对中、低纬电离层暴的产生原因进行了深入的探讨;邓忠新等[17]提出了TEC扰动指数DI,借此分析了中国地区电离层TEC暴扰动的产生缘由和分布特性;吴佳姝等[18]通过对2001—2010年间的156次磁暴事件的统计分析了欧洲扇区赤道到极光带不同纬度的电离层暴特征．除此之外,国内外学者对特定磁暴期间不同纬度地区的电离层的变化情况进行了细致的分析[19-23],例如:Oluwaseyi等[19]重点研究了赤道地区电离层在磁暴期间的变化情况,刘海涛等[20]通过构建北美地区的二维分布图像,研究了一次强磁暴期间的电离层扰动的变化情形及其原因．

虽然在电离层TEC领域已经开展了大量的研究工作,也取得了不错的研究成果,但随着科技的进步,研究的发展,必然对电离层TEC的精度和分辨率提出更高的要求．本文提出了一种高精度TEC地图重构方法,利用亚大地区60个GPS台站的实测数据,基于Kriging内插方法实现了亚大区域高分辨率TEC地图重构,通过与实测数据进行对比并对成像精度进行了验证;在此基础上,详细分析了南北半球不同纬度的电离层TEC值在磁暴期间的不同表现．

1 Kriging插值方法

1.1 半变异函数(semi-variogram)

要使用Kriging法进行插值,必须先求出半变异函数,这是用于描述插值中所用数据之间空间相关性的函数,它之所以重要是由于它是Kriging算法的主要输入量．半变异函数可通过在固定距离dl±Δdl/2处对不同观测数据做方差求得:

(1)

式中:γ*为经验半变异函数;m(dl)为在距离dl处的不同观测对的数量,Zi和Zj是相应点xi、xj在距离|xi-xj|=dl处的观测值．

一旦求得经验半变异函数,下一步便是通过满足各种数学条件转化为可以应用在Kriging方程中．

1.2 Kriging方法

Kriging是一种属于最优线性无偏估计的线性内插方法,因此,Kriging方法的主要目的就是利用已知量的线性组合来估计未知变量Z*:

(2)

式中:ωi是权重因子,对应着每个已知量Zi．

(3)

式中:E[(Z*-Z)2]是Z*的方差,能够用半变异函数来表示:

(4)

减去方程(3),可得简化后的方程:

(5)

因此,为了求得权重因子ωi,必须先计算出下列参数:

Ω=Γ-1Γ0,

(6)

式中：Ω为包含权重因子ωi和拉格朗日乘子λ的向量；Γ为包含已知量和方位半变异的矩阵；Γ0为包含未知量和方位半变异的矩阵．

2 TEC的二维分布结果与分析

2.1 数据来源

本文使用的卫星数据全部来自于IGS的全球GPS台网,从ftp://cddis.gsfc.nasa.gov/pub/gps数据中心下载而得,选取观测数据的经度和纬度范围为:90°E～150°E、50°S～50°N．此区域涵盖了约60个观测台站,观测台站接收数据的时间分辨率均为30 s．

2.2 TEC的二维分布与精度验证

由于测量精度、仪器偏差等问题的存在,在解算TEC的过程中,不同台站获得的TEC值可能会存在一定的误差,这会对TEC二维分布重构产生较大的影响．为此需要对各台站的TEC数据进行预处理,具体方法为:首先把同一段时间内(本文取15 min)的垂直TEC值按经度进行划分,经度间隔取2.5°;然后画出相同经度区间内垂直TEC随纬度的变化曲线,采用数据线性拟合的方法对曲线中相比其它数据点偏差过大的点进行滤除,以保证所有的垂直TEC曲线都较为平滑．

对所有经度区间的数据点处理过后,即可得到所测电离层垂直TEC数据的分布图,如图1所示．得到垂直TEC值的分布以后便可进行插值处理,本文使用Kriging法进行插值处理．由于数据量过大,计算过程中会出现超出计算机内存的情况,需要对原始数据进行取中值处理:对于每一对卫星接收站而言,在15 min内大约可以获得30个位于电离层穿刺点处的垂直TEC值,按每十组数据对垂直TEC取中值的方法对卫星信号进行处理,从而每对卫星接收站之间可以获得的数据量可以减少到约3组,这样可以大大降低运算内存且能优化计算结果．经过插值后便可以得到TEC的二维分布图,如图2所示．可以看出:去掉测量过程中一些偏差较大的数据可明显提高TEC值的重构质量．

图1 数据处理前后观测值分布情况

图2 由Kriging内插得到的TEC二维分布图

在获得高分辨率的TEC二维分布图后,采用以下方法对精度进行分析:在60个接收站随机选取59个接收站的数据进行TEC重构,剩下1个站的数据作为真实值对我们的成像结果进行验证,以评估本文方法的TEC重构误差．

设重构得到的TEC值为TECsimulation,TEC真实值为TECstation,则误差为

图3是TEC重构误差分布图,选取不同日期不同时段的共611个TEC值进行对比分析,从图中我们可以看出,经过插值得到的TEC地图的误差主要集中在10%以下,精度还是比较理想的,验证了本文方法的有效性和可靠性．

图3 TEC成像结果的误差分布图

2.3 实测数据成像结果分析

图4示出了东经120°,不同纬度下TEC值在一天内随时间的变化情况,很明显地出现了一天中TEC值先升高后降低的变化趋势,而且纬度越低变化趋势越明显;在白天可以很清晰地看出TEC值从赤道向纬度升高的方向变小的趋势,而晚上各纬度TEC值相差不大;且从图中可以看出白天TEC值在20°纬度附近有一个反常变大的过程,这是由于赤道异常区的影响．一般来讲,赤道异常在日出后形成于地方时10:00并逐渐增强,在地方时14:00左右达到最大值,然后又逐渐衰退,在地方时04:00左右有极小值．一天之中,异常峰先在低纬度发展出现,并于增强时逐渐移向较高纬度,在其之后衰退时,又逐渐返回低纬度,最后消失[24]．图5所示为2014年02月21日00:00-24:00UT时刻的TEC二维分布重构结果．在图中可以很明显地看出TEC极大值区域在随着时间自东向西偏移,这是因为地球自西向东自转因素的存在;由于白天电子密度会升高,白天各地区的TEC值比晚上大很多．本文重构的TEC结果较好反映了电离层的日变化和区域变化特征．

图4 垂直TEC随纬度和时间的变化结果

图5 同一天内不同时刻的TEC二维分布图

在得到亚大地区的高精度TEC二维分布图后,便可以借此分析TEC随纬度、时间等的变化情况以及在各种突发事件中的变化趋势．相比于磁静日,磁暴期间电离层会发生显著的变化．如图6所示,给出了2014年2月25日-3月5日一次磁暴事件期间电离层TEC的变化结果．从图中可以看出Dst变化情况和东径120°不同纬度处TEC变化情况，本次磁暴发生在2月27日夜晚,最小Dst值为-94 nT;其中黑色实线表示实时的TEC值,红色虚线表示静日TEC平均日变化．本文以Dst指数作为静日的选取标准,期间所有Dst指数均不小于-20 nT,以此来选择磁暴前后共7天的TEC值进行平均以作为本文的静日参考．可以看出随着磁暴的发生不同纬度的TEC值都出现了不同程度的变化,在磁暴恢复相期间,各纬度TEC值出现不同程度的减小,呈现负暴效应,尤其在低纬赤道异常区的负暴效应十分明显,这应该是由于磁暴期间出现的扰动发电机电场阻碍了赤道异常驼峰的形成所造成的[16]．由于磁暴主相持续时间较短,因此主相期间TEC的变化情况我们以图7进行详细分析．

图6 一次磁暴期间的Dst指数以及不同纬度TEC的变化情况

图7 磁暴主相期间TEC与静日TEC对比图

图7示出了磁暴主相期间TEC值与宁静日里相对应地方时的TEC值之间的差异,本次磁暴急始发生于27日17:00UT时,并于27日23:00UT时达到主相极大．26日作为我们的静日参考,从中可以看出在主相期间,大部分地区TEC值都出现了明显的增大,整体呈现正暴相的特征．不过由于磁暴期间的物理机制极其复杂,不同区域内TEC值的变化情况也有所不同．例如南北半球TEC值的变化差异极大,TEC值在北半球东亚地区升高幅度较为明显,最高可达100%,而在澳大利亚南部区域则出现了下降的趋势,下降最大百分比接近60%,这种南北半球的差异化表现可能是由于磁暴带来的带电粒子由南向北的迁移现象造成的．

图8示出了磁暴恢复相期间TEC值与宁静日里相对应地方时的TEC值之间的差异,可以看出在恢复相期间,大部分地区TEC值都出现了不同程度的下降,呈现负暴的特征,而且南半球TEC的下降幅度要比北半球大很多,这也符合图6中显示的TEC变化趋势．

图8 磁暴恢复相期间与宁静日TEC对比

由于磁暴期间复杂的物理机制,不同经纬度地区的TEC值变化情况也有所不同．例如在北半球中低纬区域有着明显的先正暴后负暴的现象,而在澳大利亚南部区域却只有负暴产生,全程未见正暴的存在．这一南北半球的差异化表现对于深入了解电离层暴的特征及其动力学过程无疑是非常有益的．

图9示出了2015年与2016年的两次磁暴期间的电离层TEC变化情况．在第一次磁暴中,磁暴急始发生于12月20日04:00UT时左右,并于22:00UT时达到主相极大,Dst指数达-170 nT．从图中可以看出,本次磁暴期间大体呈现了先正暴后负暴的现象,而不同纬度的具体表现也有着细微的不同．图10(a)示出了本次磁暴期间TEC与宁静时的差值,在主相期间,各纬度的TEC值都有明显的增大,而且有着增大幅度从赤道向两极逐渐变强的趋势;而恢复相期间不同纬度的TEC变化情况有着较大的差异,赤道与中低纬地区的TEC值有着明显的降低,而高纬区的TEC值继续保持着增大的趋势．

图9 磁暴期间TEC的变化情况

第二次磁暴发生于2016年5月7日-2016年5月9日．在本次磁暴中,磁暴急始发生于8日01：00UT时左右,并于08：00UT时达到主相极大,Dst指数达到-92 nT．图10(b)示出了本次磁暴期间TEC与宁静时的差值,在本次磁暴期间不同纬度的TEC变化趋势有较大的差异,在主相期间,赤道地区TEC变化程度不太明显,而高纬地区出现了明显的正暴现象;而在恢复相期间,南北半球出现了较大的差异,北半球中低纬区域主要出现负暴,而南半球相应纬度出现了明显的正暴,这与2014年2月27日磁暴期间北半球TEC值升高,南半球TEC值减小的现象正好相反．

由这几次磁暴事件的对比可以看出,亚大区域在磁暴期间出现先正暴后负暴的现象比较普遍,不过不同纬度区域的具体表现情况有所不同:例如主相期间高纬TEC值的升高程度要普遍大于低纬地区,而且南北半球相同纬度的TEC变化趋势也有所不同,经常会出现相反的变化趋势．

(a) 第一次磁暴 (b) 第二次磁暴图10 磁暴期与宁静时的TEC差值

3 结 论

本文提出了一种高精度TEC二维分布重构方法,基于亚大地区60个GPS台站的实测数据,利用Kriging内插方法实现了亚大区域高分辨率TEC地图重构．基于高分辨率TEC地图重构结果,本文对亚大区域几次磁暴事件期间的电离层TEC变化进行了分析,研究结果如下:

1)通过与实测结果的对比可以看出,本文的TEC地图绘制结果精度较为理想,经计算误差普遍在10%以下,与实测值平均相差2.6 TECU左右,验证了本文方法的有效性．

2)基于TEC的二维分布分析了亚大地区在不同年份的几次磁暴事件的影响,分析了不同纬度地区在磁暴主相与恢复相期间的差异化表现,大部分地区会出现先正暴后负暴的现象,高纬地区的TEC变化幅度较之低纬地区要更为明显,而且可以看出南北半球相同纬度地区在磁暴期间的表现也存在极大不同,这对进一步了解磁暴的物理机制是极其有意义的．

致谢：本文所用实测数据均下载自IGS组织提供的全球观测资料,在此感谢他们对电离层观测研究工作做出的无私贡献．