2. 中国科学院地球科学研究院, 北京 100029;
3. 中国科学院大学地球与行星科学学院, 北京 100049
2. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
3. College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
我们把60~1000 km左右的含有大量自由电子和离子的高空大气层称为电离层.由于电离层中的自由电子的数量足以对无线电波的传播产生极其重要的影响,所以电离层观测对于航天、通讯、导航等都具有非常重要的应用价值.从20世纪20年代开始,电离层垂直探测一直是探测电离层的主要手段,全球数百个垂测站点为我们积累了大量的观测数据,使得我们能够更直观地了解和认识电离层(林兆祥等, 2014; Bibl and Reinisch, 1978).
数字测高仪只能观测F2层峰值以下的电子密度(Reinisch and Huang, 1982,1983; Huang and Reinisch, 1982),为了研究顶部电离层电子密度,科研工作者们做了大量的工作(Huang and Reinisch, 2001; Reinisch and Huang, 2001; Sibanda and McKinnell, 2009; Stankov et al., 2003).目前国际上普遍采用Reinisch和Huang(2001)提出的基于Chapman方程外推顶部电离层电子密度的方法(Reinisch et al., 2004; Bilitza et al., 2006).该方法做了两次假设,首先假设电子密度在高度上的分布符合Chapman理论,然后又假设Chapman方程中的标高是一个不随高度变化的常量.但是事实证明电离层的具体物理过程远比Chapman理论要复杂的多,而且标高H是一个随着高度、地方时、季节变化的变量,随着高度的增加标高也进一步增大(Zhang et al., 2006;Liu et al., 2015; Bilitza et al., 2004).所以由Chapman方程外推的顶部电子密度必然存在很大的误差.
基于这种现状,本文提出了一种使用经验正交函数(EOF)估算顶部电离层电子密度剖面的方法.该方法使用收敛速度极快的EOF经验正交分解方法,利用IRI 2012模型输出的电子密度作为正交分解的背景场,结合峰值以下的测高仪数据、MIT TEC数据建立一种估算顶部电离层电子密度的模式.下面将依次介绍IRI 2012模型、测高仪观测数据、MIT TEC数据、EOF方法和峰值以上电子密度的建模过程,然后展现基于该模式估算美国Millstone Hill (42.6°N、288.5°E)台站顶部电子密度的结果,并将该结果与该处非相干散射雷达ISR观测数据作对比.
1 数据介绍本文选用了目前国际上发展最为完善的电离层经验模型--国际参考电离层模型(International Reference Ionosphere, IRI)作为研究背景模型.该模型由全球测高仪台站数据、国际电离层研究卫星和加拿大科学研究卫星(ISIS)以及加拿大的百灵鸟科学研究卫星(ALOUETTE)的顶部观测数据、一些卫星和火箭的就地测量数据,以及非相干散射雷达观测数据建立的标准电离层经验模型组成.当我们给定确定的时间、地点、F10.7等参数,就可以得到60~2000 km高度范围内的电子密度、电子和离子温度、离子成分以及某个高度区间内的总电子含量等电离层参数(Bilitza et al., 2011).
测高仪是目前地面垂直探测电离层的主要常规设备,它的主要原理是在地面向电离层垂直发射脉冲调制的多频、高频无线电波,根据电离层反射回波到达接收机的时间得到反射点处等离子体的密度(熊年禄等, 1999).本文使用了由数字测高仪实测的峰值以下的电子密度数据;本文还使用了总电子含量TEC数据,该TEC数据来自于麻省理工学院Madrigal镜像数据库(MIT TEC).在数据处理中,由于GNSS卫星运行高度在20200 km以上而且90 km以下的电离层电子密度较低,所以文中设定的电子密度的高度区域范围是90~20200 km.考虑到电子密度在高度上大体呈Chapman函数形式的指数变化,所以在高度方向采用了对数坐标.为了保证测高仪数据的真实性,数据处理时间方面直接保留了测高仪电子密度数据的原始时间坐标.即对IRI电子密度数据、MIT TEC数据在时间上根据测高仪数据时间进行了一次样条插值,在高度上对IRI电子密度数据、测高仪数据在90~20200 km范围内根据对数坐标也做了一次样条插值,实现了对所有数据在时间和高度上的统一.
2 经验正交函数(EOF)电离层的物理过程极为复杂,受到光化学、动力学、电动力学等多个过程的复合影响.EOF是一种能够通过对探测数据分析将多个物理过程的贡献分离出来的正交分解方法.这种方法可以剔除原始数据冗余信息,它最大的特点是可以由观测资料自身分解出快速收敛的正交基函数,该基函数能反映数据自身特点,有助于我们抓住主要变化形态.关于经验正交函数(EOF)的分解方法可参考Von Storch和Zwiers(1999)专著.本文使用EOF方法分解了IRI 2012模型Millstone Hill (42.6°N、288.5°E)台站位置处2005-2015年11年间的电子密度Ne(h),得到基函数Ej(h),如式(1)所示.
(1) |
其中h是插值后网格点对应的高度,Ne(h)是某一时刻电子密度在高度上的一维分布,Ej(h)是正交分解出来的第j阶基函数,Aj是第j阶基函数所对应的系数.
EOF正交分解最大的特点是基函数由大量的数据在计算过程中自己生成,它更能反映数据固有特性.另外,它还具有收敛速度快、基函数完全正交等特点.下面我们讨论一下EOF收敛性及前三阶基函数的基本特征.
EOF基函数的实质就是观测数据协方差矩阵的特征向量,归一化之后的特征值代表了对应基函数占全部基函数的权重.实践证明归一化后特征值会随着阶数的增大而减小.本文工作中第一阶特征值为λ1=0.9722,其代表了第一阶基函数在全部基函数中占据了极大的比重,第二至五阶的特征值分别是λ2=0.0195,λ3=0.0074,λ4=0.007,λ5=0.0001,λ5之后的特征值逐渐减小,权重之和不超过0.00008.由此可见,前五阶特征函数就可以直接反应观测数据的基本特征,所以本文工作统一截取了前5阶基函数做为研究对象.
图 1a是IRI模型输出的Milstone Hill台站位置处2005-2015年间中午12:00的平均电子密度值Ne(h).图 1b-1d是对2005-2015年中午12:00处所有电子密度采用不去均值的EOF方法分解出的前一、二、三阶基函数.从图 1b可看出第一阶基函数E1恒为正且它的形状与平均电子密度形状相似,峰值高度与平均电子密度峰值都在270 km左右,由此可以判定电离层电子密度强度主要由第一阶EOF级数A1E1(h)控制.观察图 1c可以发现350 km以下为负,350 km以上为正,若A2为正则E2在350 km以下会使电子密度减小,350 km以上调节电子密度增大,总体看来EOF第二阶级数A2E2(h)可以控制电离层高度变化.第三阶基函数E3在190~300 km之间为负,大于300 km和小于190 km电子密度为正值,也就是说当A3是个正数时,第三阶级数会降低峰值附近的电子浓度而使得顶部电子浓度和底部电子浓度增加,也就导致电子浓度剖面的厚度增加.也就是说A3E3(h)主要影响的是电子浓度剖面的厚度即板厚.由于A4E4(h)及之后的级数贡献较小且形状不规则,这里就不再讨论了.
为了展示方法的特性,我们分两种情况进行顶部剖面重建,即通过测高仪观测峰下剖面估算和基于测高仪和GNSS TEC观测估算.
3.1 通过测高仪底部电子密度估算顶部电离层首先介绍根据峰值以下电子密度估算顶部电离层电子密度过程.通过公式(1)对IRI 2012模型输出的Millstone Hill台站位置处的2005-2015年电子密度做不去均值的EOF分解,得到与时间无关的基函数Ej(h),它基本上描述了电子密度在空间上的分布.因为Ej(h)是一组收敛很快并且完全正交的基函数,完全可以建立测高仪观测电子密度Ne(h1)与基函数Ej(h)之间的观测方程,如式(2).其中Ne(h1)代表的是测高仪台站实测峰值以下电子密度,h1指的是峰值以下的高度坐标,Ej(h1)是和Ne(h1)的高度对应保持一致的基函数,ε代表误差,AjB是测高仪底部观测电子密度基于最小二乘法向上述部分高度基函数Ej(h1)做投影后求得的不包含顶部信息的基函数系数:
(2) |
求得系数AjB,又已知基函数Ej(h),根据公式(1)就可以估算出相应时间点的全高度的电子密度Ne(h).
已知估算方法,根据此方法估算美国Millstone Hill台站电子密度Ne(h)并对估算结果中的几个重要电离层参量进行讨论.把根据峰下数据重建的电子密度称为峰下估算电子密度Bot Model Ne,相对应的电离层参量叫做峰下估算电离层参量.
经实践证明使用以上方法估算的峰下电子密度与测高仪实测数据对应高度的电子密度高度吻合,此处不再赘述.为了研究峰上电子密度估算精度,本小节重点讨论了临界频率、总电子含量等电离层参数.本文利用Millstone Hill台站2015-2016年测高仪实测峰下电子密度重建了该位置处相同时间段的全高度电子密度,并提取了对应时间点的临界频率并与测高仪测得的临界频率做统计对比,如图 2所示.其中,散点图展示了两种频率对比效果,实线是两种频率的等值分界线.在等值分界线之上测高仪实测临界频率DGS foF2大于峰下估算临界频率Bot Model foF2,而在等值线之下Bot Model foF2大于DGS foF2.如图 2所示,散点图主要分布在等值线之下,也就是Bot Model foF2大于DGS foF2.可以考虑是因为顶部约束的缺少导致峰值电子密度相应的发散变大.
图 3展示了峰下估算电子总含量Bot Model TEC与实测电子总含量MIT TEC的差距.本文将估算电子密度(Bot Model Ne)在高度90~20200 km做积分得到估算总电子含量Bot Model TEC,然后将Bot Model TEC与实测MIT TEC作对比.图 3中横坐标是峰下估算总电子含量Bot Model TEC,纵坐标是GNSS实测的电子总含量MIT TEC.散点图展示了两种TEC的对比效果,直线是两种TEC的等值线.通过观察发现在总电子含量低值区域MIT TEC大于Bot Model TEC,而在高值区域Bot Model TEC大于MIT TEC.总体看来Bot Model TEC并不等于MIT TEC,且分布零散.
根据以上两种对比结果基本可以认为底部电子密度仅包含底部信息,如果只根据底部电子密度外推顶部电子密度会使数据发散,估算结果不准确.为了提高估算精度,希望在此工作中能够融入顶部信息.
3.2 联合底部电子密度、垂直TEC两种数据估算顶部电离层众所周知,总电子含量TEC是电子密度沿高度的积分,它不仅包含底部信息同样也包含顶部信息,理论上如果在上小节工作中引入总电子含量会有助于提高估算精度.为了验证这一设想,本节将Millstone Hill台站处GNSS垂直TEC加入估算模型中.首先将垂直TEC观测路径沿高度进行划分,如图 4所示.其中O代表了地心,Re是地球半径,K是地面接收机位置,G是天顶位置,观测路径与各个薄层相交于Pi点.这就等同于电离层在高度方向上划分成了多个薄层,每个薄层之间的距离是ΔZi,并且假设薄层附近的电子密度主要集中在这个薄层上并且是一个常数.这种模型就是电离层多层模型,电离层多层模型更好地考虑了电离层电子密度在高度上存在较强的高度梯度,比经典电离层薄层模型更加实用(She et al., 2017).
已知垂直TEC是电子密度沿高度的积分:
(3) |
其中ε0代表的是拟合误差,EjT是电子密度第j阶基函数Ej(hi)在高度上的积分结果,如式(4)所示:
(4) |
因此,根据已知的电子密度基函数Ej(hi)就可以求出对应的积分EjT,这样就建立起垂直总电子含量TEC与基函数Ej(hi)的观测方程如式(3).因为TEC是电子密度在整个高度上的积分量,所以TEC实际上也包含了顶部电离层电子密度的信息,也就是说根据式(3)求出的系数Aj包含了顶部信息.
于是,已知IRI电子密度,根据式(1)求出一组完全正交的基函数Ej(h),把Ej(h)代入式(4)可进一步求出基函数的高度积分EjT.已知实测峰下电子密度Ne(h1)、GNSS TEC、基函数Ej(h)、基函数高度积分EjT,根据式(5)、(6),利用最小二乘法求出公共系数Ajc.因为最小二乘法的本质是通过最小化误差的平方和寻找数据的最佳函数匹配.公式(6)中的TEC数量级是1017,而公式(5)中电子密度数量级仅仅为1011左右,为了减少TEC数据在最小二乘法中产生的权重,本文工作在求解公共解之前首先将TEC、Ne(h1)分别归一化至0~1之间.若求出基函数Ej(h)、公共解Ajc,通过式(7)可以就估算出该时间点的全高度电子密度Ne(h).
(5) |
(6) |
(7) |
上述步骤提供了一种联合底部电子密度Ne(h1)、垂直TEC两种数据估算顶部电离层电子密度Ne(h)的方法,接下来就以使用该方法估算的2015-2016年美国Millstone Hill台站处电子密度Ne(h)为研究对象,对估算结果中的几个重要电离层参量进行讨论.为了区别于上小节提出的峰下估算电子密度等概念,本节把联合峰下电子密度Ne(h1)、总电子含量TEC重建的电子密度称为估算电子密度,相对应的电离层参量叫做估算电离层参量.
图 5展示了Millstone Hill台站处2015-2016年所有时间点测高仪临界等离子体频率DGS foF2与估算等离子体频率Model foF2的统计对比,实点代表了两种频率对比的散点图,直线是对散点图的一个拟合,直线的斜率是1.01,截距是-0.36,截距太小基本可忽略不计.根据这个统计结果基本可认定估算临界频率(Model foF2)等于测高仪临界频率(DGS foF2).
图 6的展示了Millstone Hill台站处估算峰值高度Model hmF2减去测高仪峰值高度DGS hmF2得到的峰高绝对偏差(Model hmF2-DGS hmF2)直方图,时间跨度为2015-2016年.横坐标代表了误差大小,纵坐标代表落在误差区间的数据量占总数据量的百分比.很明显,Model hmF2-DGS hmF2直方图主要分布在零点右侧,也就是说Model hmF2基本全部大于DGS hmF2.为了找出原因,我们作出了Millstone Hill台站处2015-2016年间IRI模型输出的峰值高度IRI hmF2与测高仪峰值高度DGS hmF2的绝对偏差值(IRI hmF2-DGS hmF2)的直方图,如6图所示.由图可见两个绝对偏差值的直方图形状、大小基本一致,也就是说IRI hmF2峰值高度也比DGS hmF2要大,并且二者误差变化趋势也大致相同,由此可以判定估算峰值高度大于测高仪峰值高度是一种系统误差,是背景场IRI电子密度峰值高度不准确导致的.
将估算的Millstone Hill台站处2015-2016年间的电子密度(Model Ne)在高度90~20200 km做积分得到估算总电子含量Model TEC.再次引入Model TEC、MIT TEC两种TEC数据之间的对比,如图 7所示.其中散点图代表了估算总电子含量Model TEC与实测总电子含量MIT TEC的对比,对散点图做拟合后得到散点图的拟合曲线,在图 7中还给出了两种TEC的等值分界线.观察分析认为拟合曲线基本与等值线吻合,散点基本分布在等值线附近,所以Model TEC与MIT TEC大致相等.即相较于仅根据底部电子密度估算顶部电离层电子密度,引入TEC数据后可以提高数据模型的估算精准度.
本工作下载了美国Millstone Hill台站处的非相干散射雷达(ISR)的电子浓度剖面观测数据,并将该数据与所建模式估算数据作比较.ISR是一种探测电离层强有力的手段,它能够同时对电离层多种重要参数进行测量,测量范围可覆盖电离层D层到两千公里左右的高度.受到非相干散射雷达的雷达脉冲长度的影响,ISR观测数据的高度分辨率在不同时刻完全不同.我们将高度分辨率小于5 km的观测数据称为AC(Alternate Code)数据,其余的则称为SP(Single Pule)数据.AC数据在200 km以下比较精确,而SP数据在顶部电离层表现更为精确.但是Millstone Hill非相干散射雷达AC数据量少而且坏点较多(梅冰和万卫星, 2007).由于本文主要目的是估算顶部电离层电子密度,所以本文选取了SP数据进行研究.
观察估算电离层高度剖面发现在峰值以下估算电子密度(Model Ne)与测高仪电子密度(DGS Ne)高度吻合,而在峰值以上随着高度的增加DGS Ne逐渐小于非相干散射雷达电子密度(ISR Ne),而Model Ne表现的更接近ISR Ne.为了验证这一结果,本文将2015-2016年400 km以上的Model Ne、DGS Ne分别减去对应高度的ISR Ne,得到顶部电子密度误差.经计算,400 km以上Model Ne的绝对误差均值是-4.13×1010el/m3,而DGS Ne绝对误差的均值是-8.88×1010el/m3,明显400 km以上DGS Ne对比Model Ne误差更大.图 8是两种电子密度Model Ne、DGS Ne相对于ISR Ne的绝对偏差直方图.横坐标代表了误差大小,纵坐标代表落在误差区间的数据量占总数据量的百分比.可以直观地看到Model Ne与ISR Ne的差值Model Ne- ISR Ne主要分布在零点附近,而DGS Ne与ISR Ne的差值DGS Ne-ISR Ne基本上分布在零点左侧.所以,估算电子密度Model Ne基本等于ISR Ne,而DGS Ne小于ISR Ne.
本文使用EOF方法,以IRI 2012模型输出的电子密度作为背景场,结合测高仪实测峰值以下电子密度建立模型估算顶部电离层电子密度,但是结果证明仅仅根据峰下电子密度建模估算电子密度会因为顶部信息的缺失导致估算顶部电离层电子密度有系统偏差.为此本文工作引入了包含顶部信息的GNSS TEC数据,在结合实测底部电离层电子密度、GNSS TEC数据的基础上估算出全剖面电离层电子密度.将最终估算结果与测高仪实测电离层参数foF2、hmF2、MIT TEC以及ISR测得的顶部电离层电子密度对比后发现估算临界频率、峰值高度与测高仪实测数据基本一致,400 km以上估算电子密度相较于非相干散射雷达实测对应高度电子密度的绝对误差平均值仅是测高仪400 km以上电子密度绝对误差平均值的一半左右.可以清晰地得出一个结论:使用EOF方法分解IRI模型,联合实测底部电子密度、TEC数据可以更精确地估算顶部电离层电子密度.
致谢 所使用的MIT TEC、ISR数据从madrigal数据库(http://madrigal.haystack.mit.edu/madrigal/)获取,Millstone Hill台站测高仪数据从DIDBase (http://umlcar.uml.edu/SAO-X/DIDB-connect.html)数据库获取, IRI模式从http://irimodel.org/网站下载.
Bibl K, Reinisch B W. 1978. The universal digital ionosonde. Radio Science, 13(3): 519-530. DOI:10.1029/RS013i003p00519 |
Bilitza D, Huang X Q, Reinisch B W, et al. 2004. Topside Ionogram Scaler With True Height Algorithm (TOPIST):Automated processing of ISIS topside ionograms. Radio Science, 39(1): RS1S27. DOI:10.1029/2002RS002840 |
Bilitza D, Reinisch B W, Radicella S M, et al. 2006. Improvements of the International Reference Ionosphere model for the topside electron density profile. Radio Science, 41(5): RS5S15. DOI:10.1029/2005RS003370 |
Bilitza D, Mckinnell L A, Reinisch B, et al. 2011. The international reference ionosphere today and in the future. Journal of Geodesy, 85(12): 909-920. DOI:10.1007/s00190-010-0427-x |
Huang X Q, Reinisch B W. 1982. Automatic calculation of electron density profiles from digital ionograms:2. True height inversion of topside ionograms with the profile-fitting method. Radio Science, 17(4): 837-844. |
Huang X Q, Reinisch B W. 2001. Vertical electron content from ionograms in real time. Radio Science, 36(2): 335-342. DOI:10.1029/1999RS002409 |
Lin Z X, Zhang Y T, Wu Q, et al. 2014. Ionospheric data assimilation and inversion experiment with using the single ionosonde and GPS data. Journal of South-Central University for Nationalities (Natural Science Edition) (in Chinese), 33(2): 85-88. |
Liu L B, Huang H, Chen Y D, et al. 2015. Deriving the effective scale height in the topside ionosphere based on ionosonde and satellite in situ observations. Journal of Geophysical Research:Space Physics, 119(10): 8472-8482. DOI:10.1002/2014JA020505 |
Mei B, Wan W X. 2007. An empirical orthogonal function (EOF) analysis of ionospheric electron density profiles based on the observation of incoherent scatter radars at Millstone Hill. Chinese Journal of Geophysics (in Chinese), 51(1): 10-16. |
Reinisch B W, Huang X Q. 1982. Automatic calculation of electron density profiles from digital ionograms:1. Automatic O and X trace identification for topside ionograms. Radio Science, 17(2): 421-434. |
Reinisch B W, Huang X Q. 1983. Automatic calculation of electron density profiles from digital ionograms:3. Processing of bottomside ionograms. Radio Science, 18(3): 477-492. |
Reinisch B W, Huang X Q. 2001. Deducing topside profiles and total electron content from bottomside ionograms. Advances in Space Research, 27(1): 23-30. DOI:10.1016/S0273-1177(00)00136-8 |
Reinisch B W, Huang X Q, Belehaki A, et al. 2004. Modeling the IRI topside profile using scale heights from ground-based ionosonde measurements. Advances in Space Research, 34(9): 2026-2031. DOI:10.1016/j.asr.2004.06.012 |
She C L, Wan W X, Yue X A, et al. 2017. Global ionospheric electron density estimation based on multisource TEC data assimilation. GPS Solutions, 21(3): 1125-1137. DOI:10.1007/s10291-016-0580-7 |
Sibanda P, McKinnell L A. 2009. Evaluating the IRI topside model for the South African region:An overview of the modelling techniques. Advances in Space Research, 44(6): 707-714. DOI:10.1016/j.asr.2009.04.032 |
Stankov S M, Jakowski N, Heise S, et al. 2003. A new method for reconstruction of the vertical electron density distribution in the upper ionosphere and plasmasphere. Journal of Geophysical Research:Space Physics, 108(A5): 1164. DOI:10.1029/2002JA009570 |
Von Storch H, Zwiers F W. 1999. Statistical Analysis in Climate Research. Cambridge: Cambridge University Press.
|
Xiong N L, Tang C C, Li X J. 1999. Introduction to Ionospheric Physic (in Chinese). Wuhan: Wuhan University Press.
|
Zhang M L, Reinisch B W, Shi J K, et al. 2006. Diurnal and seasonal variation of the ionogram-derived scale height at the F2 peak. Advances in Space Research, 37(5): 967-971. DOI:10.1016/j.asr.2006.02.004 |
林兆祥, 张雨田, 吴祺, 等. 2014. 利用单站电离层测高仪与GPS数据的同化反演试验. 中南民族大学学报(自然科学版), 33(2): 85-88. |
梅冰, 万卫星. 2007. 基于Millstone Hill非相干散射雷达观测的电离层电子浓度剖面的经验正交函数分析. 地球物理学报, 51(1): 10-16. DOI:10.3321/j.issn:0001-5733.2007.01.002 |
熊年禄, 唐存琛, 李行健. 1999. 电离层物理概论. 武汉: 武汉大学出版社.
|