2015, Vol. 58
2. 中国电波传播研究所青岛分所, 青岛 266107
2. China Research Institute of Radiowave Propagation, Qingdao 266107, China
The advantages of function-based model and voxel-based model are combined in this new algorithm which will effectively reduce the dependence of the ionospheric tomography on the background ionosphere and reconstruct ionospheric disturbance structure of small scale as well. The absolute mean error, RMS and relative error of NmF2 derived from simulation is 0.89×1011el·m-3, 1.00×1011el·m-3 and 7.8% respectively compared with 1.6×1011 el·m-3, 2.3×1011 el·m-3 and 16.2% derived from CT result of Low-latitude Ionospheric Tomography Network (LITN) observations. The inversion results of numerical simulation and actual LEO beacon observations from LITN validate the feasibility and reliability of this new algorithm.
电子密度分布是表征电离层状态变化的一个重要参数，研究其时空变化规律和特征对卫星通信、卫星导航、空间天气等领域具有重要的理论意义和应用价值. 为了获得电离层电子密度时空分布的精细结构，美国伊利诺大学的Austen等人在1986年首 次提出了电离层层析成像(Computerized Tomography，CT)的设想，即通过LEO卫星发射信标信号结合沿子午面分布的地面台链接收的方式实现对电离层电子密度的二维层析成像反演(Austen et al.，1986，1988). 从20世纪80年代末起，美国、俄罗斯、西北欧、东亚等相继建立了基于LEO卫星信标的电离层CT观测台链并获得了有关电离层中纬槽、电离层行进式扰动等很多具有较大影响意义的探测结果(Afraimovich et al.，1992; Andreeva et al.，1992; Bust et al.，1994; Kunitsyn et al.，1995; Huang et al.，1997; Nygrén et al.，1997). 20世纪90年代初，武汉大学等与美国伊利诺大学以及中国台湾“中央”大学和中山大学合作，建立了北 起上海南至马尼拉 的国际上第一个基于LEO卫星信标信号的低纬电 离层CT探测台链(Low-latitude Ionospheric Tomography Network，LITN)并实现了东亚磁赤道异常区电离层CT成像(Huang et al.，1997; 徐继生等，2000). 此外，自2007年起，日本京都大学也相继在日本及东南亚部分国家建立了无线电信标接收站网(GNU Radio Beacon Receiver Network，GRBRN)并实现了对低纬区域的电离层扰动结构的探测(Ram et al.，2012).
对于LEO卫星信标测量而言，要实现高精度的电离层电子密度CT反演需要重点解决两方面的问题：第一，传统的电离层CT方法主要采用绝对TEC数据进行反演，而卫星信标的直接测量量为相对TEC数据. 为此，有学者研究了基于两个或多个台站的相对TEC数据估计绝对TEC的方法，如双站法(Leitinger et al.，1975)，多站法(Huang et al.，1997)等. 但由于电离层存在较大的水平不均匀性，这些方法重构得到的绝对TEC往往存在较大的误差(Watthanasangmechai et al.，2014)，从而影响到电离层CT的反演精度; 第二，LEO卫星相对地面运动的角速度较快，一般15~20 min以内即可实现一次电离层CT成像，主要用于获取电离层内部的快速变化特征，对背景电离层的敏感度较高.传统CT方法通常直接将经验电离层模型的输出设定为背景初值，由于经验模型与真实电离层间存在的较大偏差，这同样也会导致CT的反演结果存在较大的误差(Al-Fanek，2013)，无法较为“忠实”地反演出电离层的真实结构.
目前主流的电离层CT方法主要分为两大类:像素基模型法和函数基模型法(闻德保，2013). 像素基模型法将待反演的电离层区域离散化为一系列小的像素，然后在选择的参考框架和反演时间内，假定每个像素内的电离层电子密度为一常量并在此基础上进行电离层电子密度反演(Raymund et al.，1993; Pryse et al.，1998; Hernández-Pajares et al.，2000; 邹玉华，2004; Ma et al.，2005; Wen et al.，2007; Das and Shukla，2011). 像素基模型法通常使用线性代数类迭代算法实现电离层CT反演 矩阵的求解，如代数重构算法(Algebraic Reconstruction Technique，ART)、 联合迭代重构算法(Simultaneous Iterative Reconstruction Technique，SIRT)、乘法 代数重构算法(Multiplicative Algebraic Reconstruction Technique，MART)及相应的改进算法等(Wen et al，2007; 姚宜斌等，2014; 闻德保等，2014). 在给定精度较高的迭代电子密度初值情况下，像素基模型法的优点是能够很好地重建出电离层中小尺度结构，但缺点是方法本身对迭代初值的精度要求较高(Raymund et al.，1993; Ossama，2013). 为克服像素基模型电离层CT方法存在的不足，函数基模型法开始被很多学者使用(Fremouw et al.，1992; Hansen et al. 1997; Howe et al.，1998; Bhuyan et al.，2002 ; Mitchell et al.，2003; Amerian et al.，2010; Al-Fanek，2013). 函数基模型法的特点是利用一组函数来描述电离层电子密度的空间分布特性，再通过求解基函数的权重系数来重构电离层电子密度的区域变化. 函数基模型法常用非迭代类方法如随机反演(Stochastic Inversion)、广义奇异值分解(Generalized Singular Value Decomposition，GSVD)、截断奇异值分解(Truncated Singular Value Decomposition，TSVD)等方法来求解反演方程，而无需给定背景初值(Fremouw et al.，1992; Nygrén et al.，1997; Thampi et al.，2004; 欧明等，2014). 相比像素基方法，函数基模型法需要求解的未知数较少，反演结果稳定性高. 但该类方法缺点是反演结果过于平滑，有时会“掩盖”一些电离层的小尺度扰动结构特征(Nygrén et al.，1997).
因此，本文针对LEO卫星信标的特点及现有电离层CT算法存在的不足，提出了一种函数基模型与像素基模型组合的电离层CT新算法. 算法首先采用函数基模型法反演获取平滑但精度较高的电子密度，再利用像素基模型法对函数基模型法的反演结果进行二次迭代重构，该算法一方面能有效降低电离层CT方法对背景电离层模型的依赖，同时也能满足电离层小尺度扰动结构的反演要求.
2 基于LEO卫星信标的电离层层析成像新算法 2.1 电离层CT原理电离层TEC可以表示为沿信号传播路径上电子密度的积分(Austen et al.，1988)，有
由于LEO卫星信标测量得到的是相对TEC数据，若想利用绝对TEC数据实现电离层CT，首先需要估算δ0，然后再利用式(2)进行CT反演. 由于估算δ0一般存在较大的误差，这势必会影响电离层CT的精度. 为消除δ0对电离层CT的影响，本文采用差分相对TEC的方法实现电离层CT成像.
对于某个接收机的一段连续测量的数据而言，可先选择一个参考值rtec0，再将其他观测值依次对该参考值相减，得到
像素基模型电离层CT法需要求解每个网格内的电子密度值，未知数较多，加上LEO卫星信标的观测数据较为稀疏，反演结果将严重依赖于给定的背景初值的精度. 为降低背景初值对电离层CT结果的影响，本文选择函数基与像素基组合的方法进行电离层CT，具体分为如下两步：
第1步 ：先利用球谐分析(Spherical Harmonic Analysis，SHA)和经验正交函数(Empirical Orthogonal Functions，EOF)组合的函数基模型法进行电离层CT：
为连带勒让德函数(Associated Legendre Function)；(φ，λ，h)分别表示网格点的纬度、经度和高度值；fj表示第j阶的EOF函数、J为EOF函数的最大阶数、N为SHA的最大阶数、Sjnm和Cjnm为待求解的权重系数. 通过级数展开，式(5)可以表示为
相较于像素基模型法而言，函数基模型法的未知数会有较大减少，但矩阵 A 通常情况下依然是病态的. 为了获得式(9)的稳定解，本文采用TSVD方法进行求解(Hansen，1990)，计算过程如下：
第2步 ：将函数基模型的反演结果作为迭代初值，采用像素基模型的ART算法(Gordon et al.，1970)对式(4)进行第二次迭代重构，以获取更小尺度的电离层扰动结构，计算方法如下：
为利用LEO卫星信标实现电离层CT探测，2003年起，台湾“中央”大学在中国台湾、菲律宾、印尼等地各设立了LEO卫星信标接收站，建立了新的 低纬电离层层析成像网(new Low-latitude Ionosphere Tomography Network，简称LITN). 通过接收COSMIC、COSMOS 2414、 OSCAR、RADCAL、DMSP F15等卫星发射的信标信号(赵运超等，2014)，LITN可对北半球120°E子午面附近部分中低纬度区域的电离层进行CT探测. 通过对各站接收数据的连续性和数据质量的分析，本文选择LITN网的中坜(NCU)、 嘉义(Jiayi)、车城(Checheng)三个站的数据进行电离层CT反演，具体台站地理位置信息如表 1所示.
 | 表 1 LITN接收站地理位置信息表 Table 1 Geographical information of LITN stations |
首先采用数值仿真的方法验证本文电离层CT 新算法的有效性. 选择经度120°E、纬度10°N—50°N、 高度100~800 km的区域作为CT成像区域，其中 沿纬度和高度方向上网格的间隔分别设定为0.5°和25 km. 每个卫星信标地面接收站对卫星的观测截止仰角设定为15°，理想情况下，每个台站一次卫星过境期间的可观测时间约15 min. 仿真时，观测所需的电子密度“真值”由国际参考电离层模型IRI-2012给出.
利用式(1)的积分方程计算出绝对TEC，再将每个台站的绝对TEC值减去该站TEC的最小值，得到卫星信标测量的相对TEC数据. 由于观测噪声和离散误差的存在，仿真过程中在相对TEC中加入了随机噪声. 电离层CT过程中函数基模型采用了15阶SHA和2阶EOF函数. 通过分析可知，此次仿真若采用像素基模型法需要求解的未知数为 81×29=2349，而采用函数基模型法仅为(15+1)2×2=512个，不足像素基模型法的1/4，这有助于降低反演矩阵的条件数从而增强CT反演的稳定性. 在计算量上，新算法的计算量主要集中在反演矩阵的构建、截断奇异值分解及ART算法的迭代上面，在普通的个人计算机上(Intel Core I5处理器，4G内存)，3个台站的观测数据利用MATLAB软件完成一次CT成像计算只需2~3 min.
为了验证电离层CT算法对在不同时间和空间的变化条件下的电子密度重构能力，选取2012年5月1日00 ∶ 00 LT—21 ∶ 00 LT作为仿真时间，每间隔3 h进行一次成像. 图 1给出了IRI模型计算的“真实”电子密度分布，从图中可以看出电离层随着当地时(Local Time，LT)和纬度的变化特征，其主要表现为低纬度电子密度高于中纬度区域的，白天高于夜间的，其中夜间00 ∶ 00LT前后电子密度值最小，白天12 ∶ 00前后电子密度达到最大.
 | 图 1 IRI-2012模型给出的“真实”电离层电子密度分布Fig. 1 “Real” distribution of ionospheric electron density derived from IRI-2012 model |
图 2为电离层CT的反演结果，由于台站位置的限制，部分区域没有观测射线穿越，因此图中用空白显示. 从图 2可以看出，本文的电离层CT较好地重构出了电子密度随时间和空间的变化规律，反演结果随LT与纬度的变化规律与真实值较为一致，但反演的峰值电子密度高度与真实值之间存在一定的偏差，造成该偏差的主要原因分析是由于卫 星与地面接收机间缺乏水平射线，地基电离层CT的垂直分辨率较低造成的(Zhao et al.，2010).
 | 图 2 电离层CT反演得到的电离层电子密度分布Fig. 2 Distribution of ionospheric electron density derived from ionospheric CT |
电离层F2层峰值电子密度(NmF2)反演精度是验证电离层CT算法反演性能的重要指标.图 3给出了电离层CT反演的NmF2与真实值之间的比较. 其中图 3a所示为真实的NmF2分布，从图中可以看出，与电子密度的变化规律相同，真实的NmF2呈现出随纬度的增加而逐渐减小的趋势，NmF2在12 ∶ 00 LT时最大，在00 ∶ 00LT时最小. 图 3b为电离层CT反演结果，从图中可以看出，CT反演的NmF2与真实值之间的一致性较好. 图 3c为电离层CT反演相比真实值的绝对误差. 从图中可以发现20°N—24°N之间的区域在LT06 ∶ 00和LT18 ∶ 00的反演误差要高于其他时段，最大误差约2.0×1011el·m-3，且反演结果呈现出了类似的“双峰”结构，由于此区域处于磁赤道驼峰异常区北 侧，且正好处于晨昏交替时刻，电离层一般存在较大的梯度变化，因此会引起较大的反演误差. 而在28°N—32°N 之间的区域，则是正午(LT12 ∶ 00)前后反演误差最大，最大NmF2误差约3.0×1011el·m-3，分析引起该误差的原因与该区域接收机与卫星间的射线普遍仰 角相对较低，射线对反演的影响权重过小有关. 对 于整个仿真场景，电离层CT反演NmF2的绝对平均误差为0.89×1011el·m-3，RMS为1.00×1011el·m-3，相对百分比误差为7.8%. 仿真结果验证了本文的电离层CT方法的可行性.
 | 图 3 真实的电离层NmF2与电离层CT反演结果比较Fig. 3 Comparison of ionospheric NmF2 between the true value and CT results |
应该指出的是，图 2中给出的IRI模型电子密度分布主要验证本文新CT算法对电离层时变特征的重构能力，由于IRI模型是一个经验电离层模型(采用URSI作为模型的默认输入系数)，只能反映电离层的平均变化状态，因此与实际情况(一般情况NmF2 LT14 ∶ 00最大，LT06 ∶ 00最小)会存在一定的偏差.此外，由于几何因素的限制，卫星信标地面观测站与LEO卫星间缺乏水平射线，基于LEO卫星信标的电离层CT的垂直分辨率是有限的.从仿真结果来看，本文新算法反演的峰值高度比真实值要系统地低约100 km，该系统偏差与观测视角、算法及仿真所用的电离层模型等很多影响因素有关，并非恒定值，因此无法在利用实测数据进行电离层CT时对该偏差进行修正.
受太阳、地磁、大气等多种影响因素的共同作用，电离层扰动状况时有发生.为验证本文算法对不同尺度电离层扰动结构的成像能力，选定一个有扰动的槽状电子密度分布模型进行CT成像反演，该扰动模型在Chapman函数的基础上通过添加高斯型电子密度耗尽因子和扰动因子来模拟电离层大尺度的“槽”结构和小尺度的增强扰动结构(吴雄斌，1999).模型最大电子密度值、峰值高度和标高分别设定为1×1012el·m-3、300 km和60 km，电子密度剖面见图 4a所示. 由于像素基模型电离层CT算法需要在反演时先给定背景电子密度分布，为分析像素基算法对背景电离层模型的依赖性，给定两种不同的背景电离层电子密度值作为像素基模型法的迭代初值(以下分别简称为CIT-1和CIT-2). 其中，“CIT-1”给定的背景电子密度参数为最大电子 密度值、峰值高度和标高分别 设定为8×1011el·m-3、 280 km和60 km；而“CIT-2”给定的背景电子密度 参数则相应调整为8×1011el·m-3、250 km和40 km. 像素基模型法在两种不同条件下的CT反演结果分别如图 4(b、c)所示；函数基模型算法的CT反演结果(简称CIT-3)如图 4d所示；本文算法的CT反演结果(简称CIT-4)如图 4e所示.
 | 图 4 不同算法电离层CT结果 (a) 真实电离层电子密度分布； (b) 像素基算法CIT-1结果； (c) 像素基算法CIT-2结果； (d) 函数基算法CIT-3结果； (e) 本文算法CIT-4结果. Fig. 4 Ionospheric CT results calculated by different algorithms (a) “Real” electron density; (b) Result of voxel-based algorithm CIT-1; (c) Result of voxel-based algorithm CIT-2; (d) Result of function-based algorithm CIT-3; (e) Result of proposed algorithm CIT-4. |
从图 4中的电离层CT结果来看，对于存在扰动变化的电离层而言，像素基、函数基及本文算法获取的四组电离层CT结果均较好地重构出了纬度27°N—33°N附近的大尺度电离层“槽”状结构，验证了电离层CT对于大尺度电离层结构的反演能力，但在CT反演精度及对纬度30°N、高度400 km处的电离层小尺度的增强扰动结构的重构方面，各种算法的表现差别较大. 从图 4b、4c的对比可以看出，即使迭代算法相同，但在输入不同的电离层背景电子密度的条件下，像素基模型法的反演结果也存在很大区别. 图 4b “CIT-1”中，由于设定的背景电离层模型参数与“真实”电离层结构较为接近，因此其反演得到的峰值电子密度(1.06×1012el·m-3)与峰值高度(280 km)均与结构“真实”分布更为接近，且在30°N附近区域也重构出了电离层的部分扰动特征；而“CIT-2”由于设定的背景电离层模型参数与“真实”电离层结构偏差较大，其CT反演效果则明显要差于“CIT-1”的表现.从图 4c中可以看 出，“CIT-2”反演的峰值电子密度(~1.4×1012el·m-3)相比真实值过高，而峰值高度(~250 km)则过低.对比结果表明，像素基模型法的反演精度在很大程度上依赖于背景电子密度值的准确性，其反演精度存在较大的不确定性；当背景电子密度值与实际的电离层状态存在较大偏差时，即使进行多轮次的ART迭代重构，CT成像结果依然非常不理想；而在利用实测数据进行电离层CT时，由于无法事先获知电离层的结构特征，通常只能使用经验电离层模型(如IRI)的输出值作为背景模型值，这势必会影响电离层CT的精度，这也正是像素基模型法的不足之处. 另外，从图 4d可以看出，函数基模型法能够较为真实地对电离层分布及较大尺度的“槽”状结构特征进行成像，重构的电离层峰值电子密度(0.95×1012el·m-3)和峰值高度(~280 km)与“真实”值也非常接近. 函数基模型法的优点在于无需设定背景电子密度分布即可实现CT反演，因此算法的稳定性较强.但图 4d的结果同时也印证了函数基模型法重构的电子密度剖面较为“平滑”的特点，其CT结果没有重构出电离层的小尺度增强扰动结构；从图 4e可以看出，本文新算法充分结合了像素基和函数算法的优点，不管是在大尺度“槽”状结构还是在小尺度增强结构上，反演结果均与图 4a中的“真实”场景非常一致，仿真结果验证了新算法相比其他算法在重构电离层不均匀体结构方面的优越性.
3.3 实测数据验证利用LITN网NCU、Jiayi、Checheng等三个台站的实测数据进行电离层CT成像，CT算法采用与3.2节数值仿真时相同的参数设置. 由于LEO卫星轨道高度和卫星数目的限制，LITN网暂时不能进行全天候连续的电离层CT，因此本文主要对卫星过境期间的夜间黎明前时段和白天午后时段的CT成像结果进行讨论.
图 5所示为5月9日的白天下午时段(08 ∶ 27UT 转换为当地时即16 ∶ 27LT)电离层CT成像结果，图中绿色三角形表示参与CT反演的台站所在的地理纬度，白线代表用于CT精度验证时采用的台湾花莲(坐标23.9°N，121.6°E)的动态式电离层探测仪(dynasonde)的位置. 从图中可以看出，CT算法较为清晰地重构出了白天磁赤道北驼峰区域的特征，重构的峰值区域分布在16°N—19°N之间. 随着纬度的增加，电子密度表现出先增加再逐渐下降的变化趋势. 其中，驼峰区最大电子密度值约2.0×1012el·m-3，而驼峰区南侧及北侧的最大电子密度约为1.5×1012~1.6×1012el·m-3，相比峰值区有20%的下降.另外，峰值电子密度高度在驼峰区南侧与北侧均相比驼峰区也出现了不同程度的降低.
 | 图 5 2012年5月9日08 ∶ 27 UT (16 ∶ 27 LT) 的电离层CT结果Fig. 5 Ionospheric CT result at 08 ∶ 27 UT (16 ∶ 27 LT) in May 9, 2012 (16 ∶ 27 LT) |
图 6所示为5月9日的夜间黎明前时段(20 ∶ 02UT，即5月10日04 ∶ 02LT)电离层CT成像结果. 从图 6可以看出，与白天不同，在夜间较高纬区域的电离层电子密度要高于较低纬区域，其峰值出现在33°N—35°N附近，电子密度值约为3×1011el·m-3左右，当纬度高于33°N时，电子密度开始逐渐减小. 通过图 5及图 6的比较可以看出，随着时间的变化，峰值电子密度呈现出从白天的低纬区域逐渐向夜间的更高纬度区域移动的趋势. 以上电离层CT结果与电离层其他观测手段观测到的变化规律 是基本吻合的(Wen et al.，2007; 闻德保等，2014).
 | 图 6 2012年5月9日 20:02UT (04 ∶ 02LT)的电离层CT结果Fig. 6 Ionospheric CT result at 20:02 UT (04 ∶ 02 LT) in May 9, 2012 |
由于结合了函数基模型与像素基模型的优点，本文算法能够较好地反演出电离层扰动结构.如图 7所示为2012年5月8日19 ∶ 36 UT的卫星信标地面接收站获取的相对TEC测量值及电离层CT结果. 其中图 7a所示为NCU、Jiayi、Checheng三个台站测量的电离层相对TEC值，图中横坐标表示卫星星下点的纬度值；图 7b为电离层CT反演结果.从图 7a可看出，在整个卫星过顶期间，当卫星处于较低纬度时，电离层TEC的相对变化较小，以NCU站的测量结果为例，卫星纬度由5°N变化到20°N期间，电离层TEC的相对变化仅为不到2.5TECU，而在20°N到45°N期间，相对TEC的变化则增加到约12.5TECU，从电离层TEC变化的形态分析来 看，卫星与地面接收机信号穿越区的电离层电子密度应该满足低纬较小而高纬较高的变化特征，对比图 7b可以看出，电离层CT的反演结果与TEC推断的电子密度变化趋势是一致的.另外值得注意的是，当卫星位于41°N和46°N附近时，NCU站测量到了卫星信号传播路径上TEC的两处小“凸起”结构(图中圆圈处位置)，由于斜TEC反映的是信号穿越区域电子密度的积分情况，因此，只要该区域的电离层出现了较为明显的增强情况，对应的三个台站的电离层TEC一般会出现一一对应的扰动特征.由此可以推断NCU站接收的信号路径经过了两处电离层电子密度的扰动增强区. 通过计算NCU站在300 km左右的电离层穿刺点，得到这两处电子密度增强的位置约为30.8°N和34.3°N，而图 7b给出的CT结果非常准确地重构出了在31°N和34°N两处电离层电子密度增强的特征.反演结果验证了本文算法重构电离层小尺度扰动特征的可行性.
 | 图 7 2012年5月8日19 ∶ 36 UT (03 ∶ 36 LT)卫星信标TEC测量值及电离层CT结果 (a) LITN地面卫星信标接收机相对TEC测量值； (b) 电离层CT结果. Fig. 7 TEC from beacon observation and ionospheric CT result at 19 ∶ 36 UT (03 ∶ 36 LT) in May 8, 2012 (a) Relative TEC observations from ground beacon receivers of LITN; (b) Ionospheric CT result. |
在对图 7的分析过程中，选择NCU站的数据为例分析斜TEC的增强与反演Ne小尺度增强间的对应关系，主要原因在于NCU站观测的电离层斜TEC的采样点数最多，观测时间最长.而Checheng 站点卫星星下点41°N附近观测的斜TEC数据同样出现了类似的增强，计算其300 km左右的电离层穿刺点，反推得到这处电子密度增强的位置也为30.8°N，这与NCU站的分析结果是非常一致的.然而，由于Checheng站并未获取卫星星下点45°N更 往北后的电离层斜TEC数据，因此只通过Checheng 站的数据无法进一步验证NCU站反演的第二个扰动增强结构.另外对Jiayi站而言，由于其观测数据覆盖的区域太小(13°N—25°N)，其接收的LEO信标信号并未穿越对应的两处增强区域，因此无法通过其观测的电离层斜TEC分析得到电离层的增强结构.
此外值得注意的是，在图 7a第二个虚圈所示的TEC凸起结构之后约5°同样也出现了一个扰动结构，但图 7b的CT结果中并没有反映出该扰动特征，分析其原因主要在于NCU站该TEC凸起结构的星下点位置约50°N左右，计算该扰动结构对应350 km高度的对应的信号穿刺点位置应为纬度41°，从图中可以看出，在350 km高度处41°已处于电离层CT反演区域的边缘位置，由于该区域仅有少量NCU的观测射线通过，而Checheng和Jiayi站均无观测射线通过这部分边缘区域，电离层CT的视角非常有限，这极大限制了CT算法对该扰动结构的成像能力，导致CT算法无法对这一扰动结构进行清晰的成像.
为进一步验证算法的可靠性，本文对LITN网2012年5月1日至5月31日共计66组数据进行了电离层CT反演，同时利用台湾花莲站动态式电离层探测仪探测获得的NmF2数据对反演结果进行验证，结果如图 8所示，图中横轴对应时刻均为UT时. 从图中可以看出，本文算法的反演结果与实测结果较为 一致，计算NmF2的最大绝对误差6.1×1011 el·m-3，绝对平均误差为1.6×1011 el·m-3，RMS为2.3×1011 el·m-3，相对百分比误差约为16.2%.
 | 图 8 花莲站dynasonde探测的电离层NmF2与电离层CT反演结果比较Fig. 8 Comparison of ionospheric NmF2 between dynasonde observation at Hualien and ionospheric CT result |
LEO卫星信标作为一种电离层探测的有效手段，可用于电离层电子密度的快速成像，从而满足空间科学研究及空间信息系统对精细化电离层电子密度参量的需求. 本文针对LEO卫星信标的特点，提出了一种新的电离层CT方法. 该算法的特点主要表现为：(1)无需估算电离层绝对TEC值，直接利用差分相对TEC数据进行电离层CT反演;(2)利用基于SHA及EOF的函数基模型与像素基模型组合的方法进行电离层CT，同时将TSVD和ART算法综合用于反演矩阵的求解，降低了传统电离层CT方法对背景电离层模型的依赖，增强了反演的稳定性.
仿真结果表明：电离层CT结果与“真实”电子密度分布较为吻合，整体上CT得到NmF2的绝对平均误差为0.89×1011el·m-3，RMS为1.00×1011el·m-3，相对误差为7.8%.利用LITN网的实测数据对东经120°子午面部分中低纬度区域进行了电离层CT反演：其中白天数据重构得到了清晰的电离层磁赤道异常区北驼峰的结构，且发现夜间电离层峰值有从低纬向更高纬度移动的趋势；通过与实测的TEC结果的比较分析，验证了本文算法对电离层的小尺度扰动结构的重构能力.对2012年5月期间电离层CT反演结果的统计表明：本文算法反演的NmF2最大绝对误差6.1×1011 el·m-3，绝对平均误差为1.6×1011 el·m-3，RMS为2.3×1011 el·m-3，相对误差约为16.2%. 仿真实验和实测数据的CT验证了本文算法的可靠性.
必须指出的是，由于地面观测站与卫星间缺乏水平射线，基于LEO卫星信标的电离层CT的垂直分辨率较为有限(Zhao et al.，2010)，本文峰值高度(hmF2)的反演精度相比经验电离层模型(如IRI模型)并没有明显的提高.融合地基垂测仪、天基无线电掩星等数据进行多手段联合CT是提高成像垂直分辨率的重要途径(Hajj et al.，2000; 邹玉华，2004; Bust et al.，2004; Li et al.，2012)，这也将是下一步研究的方向.
致谢 本文使用的LITN卫星信标观测数据及花莲dynasonde数据均由台湾“中央”大学太空及遥测研究中心蔡龙治教授提供下载. 另外，两位匿名审稿人为论文提供了宝贵的修改意见，作者一并在此表示感谢.| [1] | Afraimovich E L, Pirog O M, Terekhov A I. 1992. Diagnostics of large-scale structures of the high-latitude ionosphere based on tomographic treatment of navigation-satellite signals and of data from ionospheric stations. J. Atmos. Terr. Phys., 54(10): 1265-1273. |
| [2] | Al-Fanek O J S. 2013. Ionospheric imaging for Canadian Polar regions. Calgary: University of Calgary. |
| [3] | Amerian Y, Hossainali M M, Voosoghi B, et al. 2010. Tomographic reconstruction of the ionospheric electron density in terms of wavelets. Journal of Aerospace Science and Technology, 7(1): 19-29. |
| [4] | Andreeva E S, Kunitsyn V E, Tereshchenko E D. 1992. Phase-difference radio tomography of the ionosphere. Ann. Geophys., 10: 849-855. |
| [5] | Austen J R, Franke S J, Liu C W, et al. 1986. Application of computerized tomography techniques to ionospheric research.//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. |
| [6] | Austen J R, Franke S J, Liu C H. 1988. Ionospheric imaging using computerized tomography. Radio Sci., 23(3): 299-307. |
| [7] | Bhuyan K, Singh S B, Bhuyan P K. 2002. Tomographic reconstruction of the ionosphere using generalized singular value decomposition. Current Science, 83(9): 1117-1120. |
| [8] | Bust G S, Cook J A, Kronschnabl G R, et al. 1994. Application of ionospheric tomography to single-site location range estimation. International Journal of Imaging Systems and Technology, 5(2): 160-168. |
| [9] | Bust G S, Garner T W, Gaussiran T L. 2004. Ionospheric Data Assimilation Three-Dimensional (IDA3D): A global, multisensor, electron density specification algorithm. J. Geophys. Res., 109: A11312, doi: 10.1029/2003JA010234. |
| [10] | Das S K, Shukla A K. 2011. Two-dimensional ionospheric tomography over the low-latitude Indian region: An inter-comparison of ART and MART algorithms. Radio Sci., 46, doi: 10.1029/2010RS004350. |
| [11] | Fremouw E J, Secan J A, Howe B M. 1992. Application of stochastic inverse theory to ionospheric tomography. Radio Sci., 27(5): 721-732. |
| [12] | Gordon R, Bender R, Herman G T. 1970. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography. Journal of Theoretical Biology, 29(3): 471-481, doi: 10.1016/0022-5193(70)90109-8. |
| [13] | Hajj G A, Lee L C, Pi X Q, et al. 2000. COSMIC GPS ionospheric sensing and space weather. Terrestrial Atmospheric and Oceanic Sciences, 11(1): 235-272. |
| [14] | Hansen A J, Walter T, Enge P. 1997. Ionospheric correction using tomography.//Proceeding of Institute of Navigation ION GPS-97. Kasas City, Missouri, USA, 249-257. |
| [15] | Hansen P C. 1990. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. SIAM J. Sci. Stat. Comput., 11(3): 503-518. |
| [16] | Hernández-Pajares M, Juan J M, Sanz J, et al. 2000. Application of ionospheric tomography to real-time GPS carrier-phase ambiguities resolution, at scales of 400-1000 km, and with high geomagnetic activity. Geophysical Research Letters, 27(13): 2009-2012. |
| [17] | Howe B M, Runciman K, Secan J A. 1998. Tomography of ionosphere: Four-dimensional simulations. Radio Sci., 33(1): 109-128. |
| [18] | Huang C R, Liu C H, Yeh H C, et al. 1997. The low-latitude ionospheric tomography network (LITN)-initial results. J. Atmos. Sol.-Terr. Phys., 59(13): 1553-1567. |
| [19] | Kunitsyn V E, Andreeva E S, Tereshchenko E D, et al. 1995. Investigations of the ionosphere by satellite radiotomography. Ann. Geophys., 13(12): 1263-1276. |
| [20] | Leitinger R, Schmidt G, Tauriainen A. 1975. An evaluation method combining the differential Doppler measurements from two stations that enables the calculation of the electron content of the ionosphere. Zeitschrift fur Geophysik, 41: 201-213. |
| [21] | Li H, Yuan Y B, Li Z S, et al. 2012. Ionospheric electron concentration imaging using combination of LEO satellite data with ground-based GPS observations over China. IEEE Transactions on Geoscience and Remote Sensing, 50(5): 1728-1734. |
| [22] | Ma X F, Maruyama T, Ma G Y, et al. 2005. Three-dimensional ionospheric tomography using observation data of GPS ground receivers and ionosonde by neural network. J. Geophys. Res., 110: doi 10.1029/2004JA010797. |
| [23] | Mitchell C N, Spencer P S. 2003. A three-dimensional time-dependent algorithm for ionosphere imaging using GPS. Ann. Geophys., 46(4): 687-696. |
| [24] | Nygrén T, Markkanen M, Lehtinen M, et al. 1997. Stochastic inversion in ionospheric radiotomography. Radio Science, 32(6): 2359-2372. |
| [25] | Ou M, Zhen W M, Yu X, et al. 2014. A computerized ionospheric tomography algorithm based on TSVD regularization. Chinese Journal of Radio Science (in Chinese), 29(2): 345-352, doi: 10.3433/cjors.2013052401. |
| [26] | Pryse S E, Kersley L, Williams M J. 1998. Electron density structures in the polar cap imaged by ionospheric tomography. Adv. Space Res., 22(9): 1385-1389. |
| [27] | Ram S T, Yamamoto S M, Tsunoda R T, et al. 2012. On the application of differential phase measurements to study the zonal large scale wave structure (LSWS) in the ionospheric electron content. Radio Sci., 47(2): RS2001, doi: 10.1029/2011RS004870. |
| [28] | Raymund T D, Pryse S E, Kersley L, et al. 1993. Tomographic reconstruction of ionospheric electron density with European incoherent scatter radar verification. Radio Sci., 28(5): 811-817. |
| [29] | Thampi S V, Pant T K, Ravindran S, et al. 2004. Simulation studies on the tomographic reconstruction of the equatorial and low-latitude ionosphere in the context of the Indian tomography experiment: CRABEX. Ann. Geophys., 22(10): 3445-3460. |
| [30] | Watthanasangmechai K, Yamamoto M, Saito A, et al. 2014. Latitudinal GRBR-TEC estimation in Southeast Asia region based on the two-station method. Radio Sci., 49(10): 910-920, doi: 10.1002/2013RS005347. |
| [31] | Wen D B, Yuan Y B, Ou J K. 2007. Monitoring the three-dimensional ionospheric electron density distribution using GPS observations over China. J. Earth Syst., 116(3): 235-244. |
| [32] | Wen D B. 2013. GNSS-based Ionospheric Tomographic Algorithms and Applications (in Chinese). Beijing: Surveying and Mapping Press. |
| [33] | Wen D B, Lü H Z, Zhang X. 2014. A new method of ionospheric tomographic reconstruction. Chinese J. Geophys. (in Chinese), 57(11): 3611-3616, doi: 10.6038/cjg20141114. |
| [34] | Wu X B. 1999. A study of algorithm and experiments for computerized ionospheric tomography in the equatorial anomaly region (in Chinese). Wuhan: Wuhan University. |
| [35] | Xu J S, Ma S Y, Wu X B, et al. 2000. Low-latitudinal ionospheric effects during a moderate storm by tomographic imaging. Chinese Journal of Geophysics (in Chinese), 43(2): 145-151. |
| [36] | Yao Y B, Tang J, Zhang I, et al. 2014. An adaptive simultaneous iteration reconstruction technique for three-dimensional ionospheric tomography. Chinese J. Geophys. (in Chinese), 57(2): 345-353, doi: 10.6038/cjg.20140101. |
| [37] | Zhao H S, Xu Z W, Wu J, et al. 2010. Ionospheric tomography by combining vertical and oblique sounding data with TEC retrieved from a tri-band beacon. J. Geophys. Res., 115: A1303, doi: 10.1029/2010JA015285. |
| [38] | Zhao Y C, Gui X C, Hong Z J, et al. 2014. The Kalman Filter imaging studies of ionosphere TEC. Chinese J. Geophys. (in Chinese), 57(11): 3617-3624, doi: 10.6038/cjg20141115. |
| [39] | Zou Y H. 2004. A study of time-dependent 3-D ionospheric tomography with ground-based GPS network and occultation observations (in Chinese).Wuhan: Wuhan University. |
| [40] | 欧明, 甄卫民, 徐继生等. 2014. 一种基于截断奇异值分解正则化的电离层层析成像算法. 电波科学学报, 29(2): 345-352, doi: 10.3433/cjors.2013052401. |
| [41] | 闻德保. 2013. 基于GNSS的电离层层析算法及其应用. 北京: 测绘出版社. |
| [42] | 闻德保, 吕慧珠, 张啸. 2014. 电离层层析重构的一种新算法. 地球物理学报, 57(11): 3611-3616, doi: 10.6038/cjg20141114. |
| [43] | 吴雄斌. 1999. 低纬电离层CT实验与算法[博士论文]. 武汉: 武汉大学. |
| [44] | 徐继生, 马淑英, 吴雄斌等. 2000. 一次中强磁暴期间低纬电离层响应的CT成像. 地球物理学报, 43(2): 145-151. |
| [45] | 姚宜斌, 汤俊, 张良等. 2014. 电离层三维层析成像的自适应联合迭代重构算法. 地球物理学报, 57(2): 345-353, doi: 10.6038/cjg.20140101. |
| [46] | 赵运超, 桂晓纯, 洪振杰等. 2014. 电离层TEC卡尔曼滤波成像研究. 地球物理学报, 57(11): 3617-3624, doi: 10.6038/cjg20141115. |
| [47] | 邹玉华. 2004. GPS地面台网和掩星观测结合的时变三维电离层层析[博士论文]. 武汉: 武汉大学. |