引言

电离层层析(CIT)技术是利用星载无线电设备发射的无线电波从外部扫描电离层,通过接收路径上的总电子含量(TEC),计算电离层电子密度的空间分布.电离层电子含量和电离层不规则结构会对雷达信号造成延迟,影响高频通信和GPS定位的准确性.利用电离层层析成像技术能研究地震、台风、磁暴等现象的电离层变化特征.

Austen等(1986,1988)首次使用低轨卫星GPS/TEC数据,使用算数迭代法ART对电离层电子密度进行二维成像.为了提高算法精度,出现了很多衍生算法,例如乘法迭代法MATR(Raymund et al., 1993)、同时迭代法SIRT(Afraimovich et al., 1992).迭代法的缺陷是严重依赖初值并对噪音误差极其敏感.

由于接收机最小仰角的限制,使得缺少接近水平方向的射线,因此求解过程为病态问题,很小的噪音误差会导致结果很大的偏离.为了克服病态问题,一些研究加入先验信息作为约束项,使用地面接收机的数据进行区域三维CIT成像.这些求解方法包括基函数方法(Zhao et al., 2005; Lin et al., 2014; Chen et al., 2015)、神经网络法(Ma et al., 2005)、约束最小二乘法(Seemala et al., 2014)、正则化法(Fehmers et al., 1998; Lee et al., 2007; Wang et al., 2016).加入时间维度构成四维CIT成像(徐继生等,2005;Xiao et al., 2012;Farzaneh and Forootan, 2017).由于观测数据不完备造成的病态问题,仍然是CIT成像的关键问题(Yao et al., 2015).

天基CIT成像虽然只能得到二维电子密度信息,但是相对于地基三维CIT成像有以下优势.首先,主要由于多路径效应的影响,地基CIT成像只有较高仰角的射线可用,缺少接近水平向的射线,而低轨道卫星的TEC数据在仰角为0的时候,数据质量依然很好,有利于改善病态性(刘裔文等,2013).其次,在海洋覆盖的区域,无法获取地基测量数据,而天基CIT则不受此限制.

本文将利用两颗伴飞的Swarm卫星载GPS双频接收机测量的TEC数据,利用Tikhonov正则化算法进行天基二维电离层层析反演,同时在两个卫星轨道平面电离层电子密度的层析成像,实现对顶部电离层电子密度的三维观测.第二节详细描述了我们的反演算法,第三节给出了数值验证以及在不同地磁活动条件下的真实数据反演结果.

1、数据和方法

1.1卫星数据以及网格划分

GPS星座包括超过24颗卫星,它位于距地表20,200 km的上空,分布在6个轨道面上(每个轨道面4颗),轨道倾角为55°.GPS接收机能接收沿着接收机与卫星之间轨道路径的双频GPS信号(f1=1.57542 GHz,f2=1.22760 GHz).在全球任何地方、任何时间都可观测到4颗以上的卫星,可以使用GPS信号测量TEC.

Swarm星座包括3颗相同的卫星.它们于2013年11月22日发射升空,初始轨道高度为500 km.自2014年1月起3颗卫星逐渐分离,在2014年4月17日达到最终星座位置.至此Swarm A和C卫星在470 km的高度伴飞,经度间隔约为1.4°(～150 km),通过赤道的时间相差不超过10 s,轨道倾角均为87.5°.Swarm B卫星的轨道高度约为520 km,轨道倾角为88°.Swarm A和C卫星覆盖全部地方时需要约133天,Swarm B卫星需要约141天.Swarm卫星带有朗缪尔探针测量等离子体密度,时间分辨率为2 Hz.所有三颗卫星均搭载双频GPS接收机,这些GPS接收机有8个频道,最多同时接收8个GPS卫星信号.GPS-TEC数据的时间分辨率为1 Hz.本文联合对Swarm A和C卫星两个平面进行反演(适用于两颗卫星之间),便于观测经度相差百公里级别,等离子体密度的差异.由于GPS和SWARM卫星的连线主要集中在中低纬区域,因此本文方法在中低纬的结果较为可靠.

Swarm卫星的飞行速度约为7.6 km·s-1,穿过南北半球±70°N(-代表南半球,+代表北半球)区域时,约需要37 min.本文假定在此时间范围内,顶部电离层电子密度不随时间变化,或者其变化可忽略.将卫星通过赤道时间记为层析成像时间.Swarm卫星与GPS卫星连线为空间三维分布,本文选择连线与Swarm卫星轨道平面夹角小于30°的TEC数据,投影到卫星轨道平面的二维网格(刘裔文等,2013).网格范围为-70°～70°纬度,480～1000 km.经线方向的分辨率为1°,高度方向的分辨率为40 km.

1.2 CIT方法

在Swarm接收机所获得的TEC数据,为信号传播路径上电子密度的积分, TEC可表示为

TEC=∫ne(s)ds,         (1)

式中TEC为总电子含量,ne(s)为电子密度随路径的变化,s为地面接收机至卫星的视线路径.对式(1)进行离散化有

yi=∑j=1Naijxj+e,i=1,2,⋯,M         (2)

式中yi表示第i条路径的TEC,aij表示第i条路径穿过第j个网格中的网格的截距,xj表示第j个网格的平均电子密度,e为误差项,M为总路径数,N为总网格数.e包含卫星和接收机的硬件时延偏差、等离子体层电子密度的贡献以及网格划分造成的偏差.对于不同的接收机,存在不同的硬件偏差,导致在同一时刻沿同一路径测量的TEC不一致,我们使用差分TEC的方法来消除硬件延时偏差(Van De Kamp,2013;刘裔文等,2013).对于每一段连续数据弧(Swarm卫星与同一颗GPS卫星的连线)硬件误差相同,我们对第k条射线和第k-1条射线做差,这样就可以消除硬件造成的误差,由(2)式可以得到

yi−yi−1=∑j=1N(a)ij−ai−1j)xj+e,         (3)

将上式写为矩阵的形式

b=Ax+e,         (4)

b是一个(M-P)×1向量,表示差分相对TEC,P表示连续数据弧段数,A表示(M-P)×N维差分系数矩阵,x是N×1向量,表示待测的各网格内的电子密度.本文同时对Swarm A/C两个卫星轨道平面进行层析成像,式(4)为两个平面方程组的联立.

由于观测条件的限制,矩阵A通常是奇异的,式(4)通常为病态问题.电子密度x对误差项e极其敏感.为了克服解的不稳定性,针对病态问题Hansen(1997)做了大量的尝试,采用正则化的方法事先对结果做合适的约束,可提高结果的稳定性.为方便引入物理约束,本文使用Tikhonov正则化方法求解式(4),其方程定义如下:

minxJ(x)=∥b−Ax∥2+λ∥Lx∥2,         (5)

其中,λ为正则项参数,Lx为引入的约束项,用于将病态问题转化为适定问题.我们使用水平约束矩阵H和垂直约束矩阵V使得反演结果符合电子密度变化规律,再通过整体约束矩阵C来控制各约束矩阵的权重,定义为

Lx=∑j=1NCj(H+V)xj.         (6)

首先我们考虑水平方向的约束条件.相对于地基电离层层析,天基电离层层析有一大优势,即可以使用卫星实测电子密度数据作为约束条件,从而提高反演结果的准确性.之前的研究(刘裔文等,2013)通常利用模式电子密度作为初值,再根据卫星实测数据缩放整个高度电子密度初值.这个方法的一个后果是卫星高度的电子密度在整个反演高度都发挥着重要影响.为减少对其他高度的影响,我们只利用卫星实测数据来构造最低反演高度的水平矩阵H.

由于无外力驱动下,电离层水平方向电子密度应是平滑连续的.在其他高度上,我们采用与其他学者(Seemala et al., 2014)的相同处理方式:

H(xA1)=−xA1−1+3xA1−xA1+1−xC1,

其中,xA1-1、xA1、xA1+1为Swarm A卫星水平方向3个连续网格电子密度.xC1为与Swarm A卫星xA1相邻Swarm C卫星网格的电子密度.由于Swarm A和C卫星是伴飞的,除了引入不同纬度方向的平滑约束xA1-1、xA1+1,我们还引入不同经度方向的平滑约束xC1.这样我们可以使得两个轨道平面的层析结果建立联系、相互约束从而提高结果的准确性.

在垂直方向,通过IRI模式引入先验信息,很多学者使用EOF分解方法提取特征向量(Zhao et al., 2005; Lin et al., 2014),这样的好处是可以极大地减少未知数的数量,但在磁暴期间的表现并不理想.本文采用相比特征向量略微宽松的约束条件,使垂直方向相邻网格电子密度的变化关系尽量与模式一致:

V(xA2)=xA2−x0A2x0A2+1xA2+1,

其中,xA2、xA2+1为Swarm A卫星垂直方向2个连续网格的电子密度,x0A2、x0A2+1为与之对应IRI模式的电子密度.同时此方法有利于我们对不同的区域采用不同的约束强度.

最后我们考虑总约束矩阵C的选择.以往学者将约束条件作用于整个区域(Zhao et al., 2005;Lee et al., 2007; Lin et al., 2014),我们认为对不同区域使用不同的约束强度更为合理(Seemala et al., 2014).在密度较大的区域,实际电子密度与模式偏离的绝对量往往较大;在密度较小的区域,实际电子密度与模式偏离的绝对量往往较小.因此我们在不同区域引入不同的约束强度.在密度较大的区域,我们允许电子密度与模式的偏离更大,施加的约束较弱,相应的Cj值较小;在密度较小的区域,我们施加的约束较强,相应的Cj值较大.这样的约束有利于提高反演结果整体的相对偏差.

本文研究中,我们利用IRI模式密度确定矩阵C.图1a给出了2014年8月20日05∶00UT,IRI模式在10°N、40°N和60°N,142°E电子密度随高度的分布,图1b给出了与之对应Cj的取值.从图中我们可以看出Cj呈二维分布同时受纬度和高度的影响.电子密度随高度增加而减小,而Cj随高度增加而增加.在500～600 km高度,电子密度较大,Cj比较小,在这些区域电子密度的变化受IRI先验信息的影响较小.随着纬度的增加,电子密度逐渐减小,Cj逐渐增加,正则项‖Lx‖2的约束作用也逐渐增强.

正则项系数λ是另一个重要的参量,决定了残差项‖b-Ax‖2和正则项‖Lx‖2的权重.λ大代表正则项重要性大(强约束),λ越小代表残差项重要性大(弱约束).有三种常用的方法选择最佳的λ,分别为无偏风险估计预测法(Mallows,1973),L-曲线法(Hansen and O′Leary, 1993)和广义交叉验证法(Golub et al., 1979)等.其中无偏风险估计预测法需要知道噪音的方差,使用广义交叉验证法经常导致结果欠平滑,使用L-曲线法的结果较为理想.本文使用L-曲线法,即在残差项和正则项之间找到一个平衡点.

2、反演结果

2.1数值验证

为了验证反演方法的有效性,我们借助IRI模式数据进行数值验证.日侧使用2014年8月20日05∶00UT的142°E的IRI模式电子密度作为x0.并依据方程y0=Ax0计算矩阵y0,将y0视为不含误差的TEC值.使用2014年8月1日05∶00 UT的160°E的IRI模式电子密度作为先验信息,记为x1,用于构建矩阵C和V.与日侧的方法类似,夜侧使用2014年8月20日05∶42 UT的50°W的IRI模式电子密度作为模式值x0,使用2014年8月1日05∶42 UT的60°W的IRI模式电子密度作为先验信息.

图2a给出了电子密度模式值(x0)与先验信息(x1)偏差的绝对值,即|x0-x1|.左半图为日侧,右半图为夜侧.可以看出日侧在22°N 500 km附近,最大偏差可达1.2×1011el./m3,均方根误差(RMSE)约为1.64×1010 el./m3,最大相对偏差为15.8%,平均相对偏差约为7.7%.夜侧在18°N 500 km附近,最大偏差约为6×1010el./m3,RMSE约为9.4×109 el./m3,最大相对偏差约为33.0%,平均相对偏差约为10.3%.

图1 (a) 2014年8月20日05∶00 UT IRI模式在10°N,40°N,60°N,142°E电子密度剖面; (b)与之对应的总约束系数

图2 (a)模式值和先验信息的绝对偏差; (b)模式值和CIT结果(未考虑误差项e)的绝对偏差; (c)模式值和CIT结果(包含误差项e)的绝对偏差

使用不含偏差的TEC(y0),可求解出电子密度x2.图2b给出了电子密度模式值(x0)与反演结果(x2)偏差的绝对值,即|x0-x2|.可以看出两者相当接近,日侧最大偏差约为6.8×109 el./m3,RMSE约为1.5×109 el./m3,最大相对偏差约为6.8%,平均相对偏差约为0.9%,夜侧最大偏差约为1.5×109 el./m3,RMSE约为5.2×108 el./m3,最大相对偏差约为4.5%,平均相对偏差约为1.2%.由此看出,在不引入误差的情况下,反演结果与模式值x0极为接近.

为了进一步比较模式值和CIT结果(图2c),图3a显示了2014年8月20日05∶00UT在三个不同高度处(600、720和840 km)的模式值(x0)和CIT值(x2)电子密度分布.可以看到在三个高度上,反演结果和模式值均有很好的吻合.

由于式(4)中A往往为病态矩阵,我们还需要测试误差项e对反演结果的影响.Swarm卫星观测到的顶部TEC含有等离子体层的贡献,等离子体层对TEC的贡献约为1～3.3 TECU(Lee et al., 2007),在电离层反演过程中,我们将等离子体层的贡献计入测量误差.假定系统误差e为均值为2 TECU、方差为1 TECU2的正态分布.使用含偏差的TEC(y0+e),可求解出电子密度x3.图2c给出了考虑系统误差后,电子密度模式值(x0)与反演结果(x3)偏差的绝对值,即|x0-x3|.可以看出虽然偏差较图2b略有增大,但模式值和反演结果仍然很接近,日侧电子密度最大偏差约为1.0×1010el./m3,RMSE约为2.0×109 el./m3,最大相对偏差约为5.8%,平均相对偏差约为1.2%,夜侧电子密度最大偏差约为1.2×1010el./m3,RMSE约为2.9×109 el./m3,最大相对偏差约为7.8%,平均相对偏差约为2.8%.

表1列出在不同观测误差条件下,CIT值和模式值之间电子密度的偏差.可以看出,日侧的电子密度高于夜侧,日侧电子密度的绝对误差也高于夜侧.引入误差项e后,反演计算的电子密度的偏差均有所上升,但仍在可接受范围内.即使引入均值为2方差为3的误差项e后,最大相对不超过15%,平均相对偏差不超过10%.这表明本文算法对通常观测误差具有较强的包容性.

表1不同误差项下,电子密度最大偏差、RMSE、最大相对偏差和平均相对偏差

图3 2014年8月20日05∶00 UT东经142°模式值和CIT值电子密度分布,至上而下分别代表600,720,840 km高度.

(a)x0为IRI模式值;(b)x0为IRI模式值在10°N处人为引入一个的耗空

IRI模型主要描述电离层的大尺度空间结构,为了验证本文方法对电离层不规则结构的重现能力,我们在IRI模式预测值的基础上,在10°N处人为引入一个密度耗空.图3b显示了引入耗空后的模式值(x0)和CIT值(x2)电子密度分布.可以看到,我们的重建方法较好地捕捉到了百公里级的电离层场向不规则结构,佐证了我们引入约束的有效性.

2.2实测数据反演结果

在本节,我们将逐一给出利用Swarm卫星实测的顶部TEC观测数据进行电离层电子密度的反演结果.

2.2.1地磁平静期

图4a为2014年8月20日05∶42 UT Swarm-A卫星轨道平面约50°W的CIT电子密度剖面,Dst指数最小值为-9 nT.图4b为Swarm C卫星轨道平面约48°W的CIT电子密度剖面.A卫星和C卫星的轨道平面仅相差约150 km,两个卫星平面的CIT结果比较接近.我们可以看出,在25°N附近电子密度有一个明显的峰值.图5比较了CIT 480 km高度电子密度与Swarm A和C卫星约460 km高度实测电子密度.除了Swarm A卫星的反演结果与实测值在25°N附近的峰值绝对量略有差异,相对偏差为12%,其他位置无论量值还是随纬度的变化趋势均高度一致.Swarm卫星在60°N附近观测到的槽状结构在反演结果中也得到了较好地重现,这体现了反演计算过程中水平约束矩阵H的有效性.

为了验证不同高度下的反演结果,图6给出了2014年8月20日约05∶42 UT反演结果与COSMIC卫星掩星测量结果的比较.图6a为COSMIC掩星测量(蓝线)和Swarm C卫星轨道平面CIT(红线)在9°S附近电子密度剖面.图6b为与之对应的经度信息;图6c为COSMIC掩星测量与CIT电子密度误差;(d)为COSMIC掩星测量(蓝线)和Swarm A卫星轨道平面CIT(红线)约52°N电子密度剖面.(e)为与之对应的经度信息.(f)为COSMIC掩星测量与CIT电子密度误差.

图4 2014年8月20日05∶42 UT的CIT电子密度剖面; (a) Swarm A卫星轨道平面约50°W;

(b) Swarm C卫星轨道平面约48°W

图5 2014年8月20日05∶42 UT的CIT 480 km高度电子密度与Swarm A卫星约460 km高度实测电子密度.

(a) Swarm A卫星轨道平面约50°W; (b) Swarm C卫星轨道平面约48°W

图6 (a) COSMIC掩星测量(蓝线)和Swarm C卫星轨道平面CIT(红线)在9°S附近电子密度剖面; (b)与之对应的经度信息; (c) COSMIC掩星测量与CIT电子密度误差; (d) COSMIC掩星测量(蓝线)和Swarm A卫星轨道平面

CIT(红线)约52°N电子密度剖面; (e)与之对应的经度信息; (f) COSMIC掩星测量与CIT电子密度误差

图6a两组密度数据,时间相差约22 min,经度相差0～2°,最大电子密度误差约为1.98×1010el./m3,平均电子密度误差约为5.7×109 el./m3,最大相对误差约为20.1%,平均相对误差约为12.6%.图6d两组密度数据,时间相差约22 min,经度相差0～2°,最大电子密度误差约为1.38×1010el./m3,平均电子密度误差约为8.4×109el./m3,最大相对误差约为43.7%,平均相对误差约为30.1%.考虑到反演结果和COSMIC掩星测量数据在观测时间与空间位置上均不完全吻合,上述电子密度量值上的差异应该在合理范围内.

2.2.2磁暴期间反演结果

(a) 2015年3月17—18日磁暴

图7所示2015年3月17—18日发生了两次磁暴,磁暴的急始发生在17日06∶00 UT,第一个较小磁暴的主相极大发生在17日09∶37 UT,Dst最小值为-101 nT,AE为指数～1016 nT,Kp指数为5+.在17日～12∶00 UT,发生第二个磁暴,Dst最小值为-231 nT,AE达到指数～2145 nT,Kp指数达到8-.

图8a为2015年3月17日09∶20 UT Swarm A卫星轨道平面约154°E的CIT电子密度剖面,图8b为Swarm C卫星轨道平面约156°E的CIT电子密度剖面.该时段Kp值为5+,Dst约-100 nT.和地磁平静期不同,Swarm A和C卫星轨道平面的反演结果在10°S—35°N观测到两段等离子体密度耗空.Swarm A卫星轨道平面的等离子体密度耗空出现在8°S—2°N和28°N—32°N区域.Swarm C卫星轨道平面的等离子体密度耗空出现在1°N—16°N和24°N—28°N区域.Zhou等(2016)使用Swarm和C/NOFS卫星在09∶30 UT均观测到赤道附近的等离子体耗空,推测这些不规则结构是由东向的快速穿透电场引起的.

为了验证其他高度的反演结果,我们将其与DMSP实测结果进行对比.在2015年3月17日约09∶20 UT Swarm A卫星和DMSP F17卫星同时通过40°S—70°S 140°E—160°E区域.图9a显示DMSP F17卫星(蓝线)与Swarm A卫星(红线)的运行轨迹.两者相差7 min先后通过65°S,155°E位置.图9b F17卫星(蓝线)在约840 km高度测得的电子密度,Swarm A卫星(红线)在对应高度下的CIT值.两者同时在40°S—70°S观测到一个槽状结构,而且量值非常接近,相对偏差约10%,验证了在840 km高度上的反演结果是可信的.

(b)2017年9月7—8日磁暴

图10所示2017年9月7—8日发生了两次磁暴,磁暴的急始发生在7日00∶00 UT,第一个磁暴的主相极大发生在8日01∶00 UT,Dst最小值为-146 nT,AE达到指数～24∶00 nT,Kp指数达到8.在8日～13∶00 UT,发生第二个磁暴,Dst最小值为-115 nT,AE达到指数～2400 nT,Kp指数达到8+.图11a为2017年9月8日01∶28 UT Swarm A卫星轨道平面约48°W的CIT电子密度剖面,图11b为Swarm C卫星轨道平面约47°W的CIT电子密度剖面.该时段Kp值为8+,Dst约-140 nT.与地磁平静期不同,Swarm A和C卫星轨道平面的反演结果在40°S—20°N范围存在不规则结构.Swarm A卫星轨道平面的等离子体密度耗空出现在纬度0°和32°S附近.Swarm C卫星轨道平面的等离子体密度耗空出现在4°S和28°S附近.Aa等(2019)使用Swarm和DMSP F17卫星在01∶28 UT均观测到赤道附近和南半球中纬区域的等离子体耗空,认为赤道附近不规则结构是由东向的快速穿透电场引起的,南半球中纬区域的等离子体耗空可能与东向的快速穿透电场及行进式电离层扰动有关.

图7 2015年3月17—18日AE、Dst、Kp指数的变化

为了进一步比较不同高度的电子密度分布特征,图12显示Swarm A和C卫星轨道平面电子密度在四个高度(在460、520、640和840 km)的纬度分布.其中460 km为Swarm卫星实测值,其他高度为CIT值.

从反演结果不难发现:(1)低纬区域的电子密度随纬度发生剧烈扰动,表明低纬区域出现大尺度不规则结构,这与Aa等(2019)的观测结果相吻合;(2) Swarm A和C两个轨道平面电子密度存在明显的差异,反演结果在520 km高度上,Swarm A卫星平面电子密度耗空出现在32°S和0°附近,Swarm C卫星平面电子密度耗空出现在28°S和4°S附近.这种差异随着高度增加而逐渐减弱,这表明本文反演方法不仅能重现顶部电离层百公里级别的不规则结构,还能区分～150 km的两个卫星轨道平面电子密度的差异.

图8 2015年3月17日09∶20 UT的CIT电子密度剖面. (a) Swarm A卫星轨道平面约154°E; (b) Swarm C卫星轨道平面约156°E

图9 (a) 2015年3月17日约09∶20 UT DMSP F17卫星(蓝线)与Swarm A卫星(红线)的运行轨迹; (b) F17卫星

(蓝线)在约840 km高度测得的电子密度,Swarm A卫星(红线)在对应高度下的CIT值

3、总结

本文利用两颗伴飞的Swarm A和C卫星GPS/TEC数据进行电离层层析成像.在正则化求解过程中,我们引入了水平约束矩阵H和垂直约束矩阵V实现相对宽松的电子密度空间变化约束,引入整体约束矩阵C调节不同区域权重因子,通过L-曲线法确定正则化系数.

图10 2017年9月6—9日AE、Dst、Kp指数的变化

图11 2017年9月8日01∶28 UT的CIT电子密度剖面. (a) Swarm A卫星轨道平面约48°W; (b) Swarm C卫星轨道平面约47°W

图12 2017年9月8日01∶28 UT的Swarm A和C卫星轨道平面的电子密度剖面.460 km为卫星实测值,520,640,840 km为CIT值

我们通过数值方式验证了本文反演方法求解病态问题的能力.利用Swarm卫星实测数据,分别对比了地磁平静期和扰动期间反演结果与第三方测量数据间的异同.结果表明本文方法不但能够较好重现实测电子密度空间结构和分布,尤其是在强磁暴期间,本文反演方法较好的重现了顶部电离层的百公里级不规则结构,而且能构区分经度方向相隔～150 km的两个卫星轨道平面电子密度的差异.

未来的工作有两个方向可以提高反演结果的准确性.其一是进一步改进约束模型,以提供更准确的先验信息.其二结合地面观测数据作为约束条件,提高反演结果的准确性.

周康俊,蔡红涛,宋方雷,谷骏,潘建宏.基于Swarm卫星GPS/TEC观测的顶部电离层三维电子密度分布的层析法研究[J].地球物理学报,2020,63(09):3219-3230.

基金:国家重点研发计划(2018YFF013702);国家自然科学基金项目(41374159,41431073,41574148)资助