1 引言

双频信标电离层观测起源于20世纪60年代的卫星导航技术，70年代后期由于 GPS(Global Positioning System)技术的提出，使得信标导航卫星的主要功能淡化，信标技术主要转向电离层探测，Leitinger等(1975)曾使用双站法完成电离层TEC的测量.随着电离层探测技术的深入发展，双频信标技术精度瓶颈问题越来越突出，使其无法很好的完成一些卫星科研任务，于是三频信标技术应运而生.

1997年，中国台湾和美国联合制定了COSMIC(Constellation Observing System for Meteorology Ionosphere and Climate)计划，它是一项多任务科学卫星，致力于包括大气层、电离层、地球引力场等多学科领域的研究(Lei， et al.， 2007).2006年该卫星发射成功，COSMIC携带有美国海军研究实验室研制的三频信标发射机(Tri-B and Beacon，简称TBB).TBB发射三种频率(150 MHz，400 MHz和1067 MHz)，其发射的信号能被地面或其他低轨卫星上的接收机接收.根据不同频率之间的相位差，计算出发射机到接收机路径上的总电子含量.此外，利用三频信标接收机提供的稠密、精确和全球性的电子密度测量数据，有利于开展空间天气物理模型的发展研究.同时，COSMIC得到的大量高精度电离层观测数据将推动空间天气预报技术的研究(郭鹏等，2002).

三频卫星信标技术的发展，为高精度的电离层电子密度成像研究带来了新的希望.三频信标有三个频率载波，通过三个频率之间的两两差分并结合有关的数学处理，使得相位模糊系数的求解精度大幅提高.2012年，中国科学院上海天文台联合中国地震局地震预测研究所和温州大学分别在上海和温州建立了两个卫星三频信标接收站，用于中国中低纬区域电离层探测，其三频信标接收机为ITS33s(黄飞，2013).

国内的三频信标尚处于研制阶段，许多关键问题仍未解决，一些研究者(吴健，2000；徐继生等，2005；夏明清，2012；陈必焰，2012)在这方面取得了重要进展.本文利用三频信标接收机的三频数据结合双频信标探测原理，提出了一种高精度电离层TEC测量的新方法，并以电离层NeQuick模型(王军等，2007)为背景场，利用卡尔曼滤波方法对三频信标TEC数据进行反演，获得二维电子密度分布.为高精度的太空探测活动打开方便之门.

卡尔曼滤波作为一种最优状态估计方法，它是解决实时应用问题的一种重要手段.因此，通常情况下，使用卡尔曼滤波解决实时问题时要考虑“时间”这一变量.由于本实验是静态的二维平面反演，即基于某一天具体时刻的观测数据进行实验.因此，本实验中用的是不考虑时间变量的简化的卡尔曼滤波模型.模型是根据相关文献(崔锦泰和陈关荣，2012)推导而来，演算过程如下：

本次实验收集了COSMOS 2414卫星在UT时间2013年1月20日18 ∶ 00左右(北京时间为凌晨2 ∶ 00左右)运行穿过中国上空时的数据，其绕行轨道大约为120°E.实验选取了Shanghai、Wenzhou、Kim、Jayi、Ludao、Checheng 6个接收站的6份数据，所有数据资料的接收时间在30 min以内.因此，可假设电离层电子密度分布不随时间变化，将其放在一起作为反演的输入量.对于初始TEC的处理，本文采用的办法为减最小项差分TEC，即令同一台站的每一条射线TEC减去所有射线中最小TEC.

其中COSMOS 2414卫星的各参数(黄飞，2013)如表 1.

表 1(Table 1)

表 1 COSMOS 2414卫星特征 Table 1 Feature of COSMOS 2414 周期(min) 飞行高度(km) 倾角(°) 频率偏置(10-6) 105.0 966 83.0 -200

表 1 COSMOS 2414卫星特征 Table 1 Feature of COSMOS 2414

本次实验要用到的数据从Shanghai、Wenzhou、 Kim、Jayi、Ludao、Checheng这6个接收站中获得，该6个台站的经纬度与高度信息如表 2.

表 2(Table 2)

表 2 各接收站地理位置信息表 Table 2 Each station location information table 站点 北纬(°N) 东经(°E) 高度(m) Shanghai 31.100 121.200 28 Wenzhou 27.920 120.690 60 Kim 24.411 118.292 44 Jayi 23.566 120.466 114 Ludao 22.675 121.492 40 Checheng 22.050 120.698 32

表 2 各接收站地理位置信息表 Table 2 Each station location information table

实验选取反演区的纬度范围为17°~36°，纬度间隔为1°；高度范围为100~800 km，网格高度划分为：100~200 km以及400~500 km间隔为20 km，200~400 km间隔为10 km，500~800 km间隔为50 km.

纬度范围17°~36°的电离层估测电子密度(以下称为初始电子密度)可由以上6个接收站的接收数据得到.其流程如图 1.

图 1

Fig. 1

图 1 数据处理流程图 Fig. 1 The flow chart of data processing

图 1流程解析：通过观测数据剖离出两行元素和初始相对TEC，再由两行元素计算出卫星在某一固定时刻的地固坐标，再由地固坐标和卫星接收站的坐标计算卫星在该时刻的方位角.由于卫星发射的信号在被地面接收站接收之前，可能产生漂移，因此还需要对初始相对TEC结合之前计算得到的卫星方位角进行误差纠正，此时计算得到无偏相对TEC，利用无偏相对TEC减最小项差分法计算出估测TEC(用b表示)，这里，减最小项差分法可以帮助我们消去在计算TEC时的产生的积分常数.最后，使用N_eQuick模型得到初始电子密度(用aN_e表示).其中，N_eQuick模型中的太阳活动参数(用flx表示)取值待定.

在进行电离层电子密度反演的实验中，需要计算出卡尔曼滤波模型中的观测矩阵 A，背景误差协方差矩阵 P，以及观测误差协方差矩阵 R，A 、 P 、 R 的计算方式如下：

A 矩阵的计算：当卫星发射的信号射线穿过反演区网格时，利用计算射线在每个网格中的长度来推导 A 矩阵.其中，A 矩阵的行维数为所有接收站接收到的射线总条数，A 矩阵的列维数为反演区的网格数，A 矩阵的每个分量的值就是对应的射线在反演区网格中所截得的线段的大小.在计算 A 矩阵 之前，需先规定好步长.在本实验中，步长取为0.1 km. 通过图 2清楚表达 A 矩阵.

图 2

Fig. 2

图 2 A 矩阵示意图 Fig. 2 A matrix diagram

由图 2可知，射线X只穿过网格2、3、4、5网格，假设射线在这4个网格中的长度分别是a₁、 a₂、a₃、a₄，则射线X所确定的 A 矩阵为[0，a₁，a₂，a₃，a₄，0，0，0，0].

P 矩阵的计算： P 矩阵中的每个分量值使用以下公式(Bust et al., 2004)计算：

R 矩阵通过经验给出，本实验取 R = k E，其中，E 表示单位阵，其维数为所有卫星信号射线数量之和，k =0.1×dim(E).

利用卡尔曼滤波模型对上面得到估测电子密度aN_e进行过滤并反演重建，过程如下：

第一步：计算出卡尔曼增益矩阵 K _g，K _g= P × A ^T×(A × P × A ^T+ R)；

第二步：计算反演电子密度cN_e，cN_e=aN_e+ K _g×(b－ A ×aN_e)；

第三步：对第二步得到的cN_e进行电子密度的反演重建.

最后，将反演值与苏州的测高仪数据进行对比，比较之后得出相关结论.其中，各参数含义可参照前文所述.

A组实验：实验分为A1、A2、A3三组小实验.实验的思路是：先取定flx的值，然后让k_α的值变化，通过与测高仪值比较，找到最接近测高仪数据的k_α值.

A1 当取flx=95时，误差随k_α的变化如图 3.

图 3

Fig. 3

图 3 A1实验误差变化曲线 Fig. 3 Curves of A1 experimental error

由图 3可知，当k_α=0.37时，反演结果与测高 仪数据的误差最小.图 4即为k_α=0.37时的图像.

图 4a表示估测电子密度重建图像，作为模型初 始值的仿真图，以下统称为初始电子密度图像，简称初始图；图 4b表示经过卡尔曼滤波过滤之后的电子密度重建图像，以下统称为反演电子密度图像，简称反演图；图 4c表示前两幅图差的百分比，即(cN_e－aN_e)/aN_e，图 4c可以反映经过卡尔曼滤波过滤之后反演图相对初始图的改善，以下统称为差图(后同).

图 4

Fig. 4

图 4 A1电离层卡尔曼滤波反演图 Fig. 4 Ionosphere Kalman filtering inversion figure of A1 experiment

本次实验，估测电子密度aN_e在苏州测高仪地理位置的峰值电子密度为M_aNe=1.2900×10¹¹el·m^-3，误差为|M_aNe－N_m|÷N_m×100%=10.01%，这里N_m表示由苏州测高仪数据计算的峰值电子密度值.反演电子密度为M_cNe=1.3547×10¹¹el·m^-3，误差为5.50%.

为证明实验结果的可靠性，可利用反演残差作进一步验证.计算可得，初始值的反演残差为ΔTEC₀=202.84×10¹⁶el·m^-2，经卡尔曼滤波过 滤之后的反演残差为ΔTEC₁=103.91×10¹⁶el·m^-2. 由此可知，反演值较初始值有一定改善.

A2 当flx=85，k_α=0.44时，反演结果与测高仪数据的误差最小.初始电子密度和反演电子密度分别与测高仪数据对比，其误差分别为：16.99%和14.36%.初始值与反演值的反演残差分别为 ΔTEC₀=130.18×10¹⁶el·m^-2和ΔTEC₁=85.74×10¹⁶el·m^-2.

A3 当flx=75，k_α=2.5时，反演结果与测高仪数据的误差最小.初始电子密度和反演电子密度分别与测高仪数据对比，其误差分别为：23.26%和21.82%.初始值与反演值的反演残差分别为ΔTEC₀=79.26×10¹⁶el·m^-2和ΔTEC₁=48.63×10¹⁶el·m^-2.

根据以上A组实验结果，可以得到以下相关结论：

(1)反演精度受电子密度的初始场误差影响，当初始电子密度的误差较小时，反演结果较好.如A1组的反演图像对初始电子密度的改进程度均要好于A2、A3两组；

(2)由图 3—4的电离层卡尔曼滤波反演图的差图知道，数据丰富的区域经过卡尔曼滤波过滤之后的电离层改善更加明显.表 2数据说明，本次试验的接收站集中分布在20°—25°纬度之间，此区间的数据反演结果与初始值相比改变更为明显.

下面再取定k_α的值，观察不同的flx初值对反 演结果产生的影响.为此，进行B组实验，B组实验包括B1、B2、B3三组小实验.在B组实验中取定k_α=0.3.

B1 取flx=95，此时的实验结果见图 5.

图 5

Fig. 5

图 5 B1实验电离层卡尔曼滤波反演图 Fig. 5 Ionosphere Kalman filtering inversion figure of B1 experiment

初始电子密度和反演电子密度与测高仪数据对比之后的误差分别为：10.01%和5.54%.初始值与反演值的反演残差分别为ΔTEC₀=202.84×10¹⁶el·m^-2和ΔTEC₁=110.14×10¹⁶el·m^-2.

B2 取flx=85时，实验结果见图 6.

图 6

Fig. 6

图 6 B2实验电离层卡尔曼滤波反演图 Fig. 6 Ionospheric Kalman filter inversion figure of B2 experiment

初始电子密度和反演电子密度与测高仪数据对比之后的误差分别为：16.99%和14.43%.初始值与反演 值的反演残差分别为ΔTEC₀=130.18×10¹⁶el·m^-2和ΔTEC₁=91.93×10¹⁶el·m^-2.

B3 取flx=75时，实验结果见图 7.

图 7

Fig. 7

图 7 B3实验电离层卡尔曼滤波反演图 Fig. 7 Ionospheric Kalman filter inversion figure of B3 experiment

初始电子密度和反演电子密度与测高仪数据对 比之后的误差分别为：23.26%和22.77%.初始值与反 演值的反演残差分别为ΔTEC₀=79.26×10¹⁶el·m^-2和ΔTEC₁=70.90×10¹⁶el·m^-2.

从以上B组实验中可以知道：

(1)flx作为实验的输入参数，其大小决定了初始电子密度场的大小.在k_α不变的情况下，随着初始场的减小，反演值与初始值之间的差值也在不断减小.由B组实验图的差图可以看出，最大差值由实验B1的30%减小到实验B3的5%左右，即反演效果在减弱.

(2)B组三个实验得到的误差随着初值的减小在增大，并且反演误差与初始误差在慢慢靠近.

通过(1)、(2)两点的分析知道，B1实验的反演差值跨度最大，B3最小，但B1、B3都不是理想的反演结果.B1反演差值跨度太大，容易造成反演结果失真；B3反演差值跨度太小，效率较低.B2的反演结果相对比较可靠，其图像反演过度更趋平稳.因此，B2反演结果要比B1、B3的可靠性更高.

此次实验收集了COSMOS 2414卫星在UT时间2013年1月25日4 ∶ 00左右(北京时间为12 ∶ 00左右)飞过中国上空时的数据，其绕行轨道大约为120°E.实验选取了Shanghai、Wenzhou、Ludao、Checheng四个接收站的三频数据，所有数据资料的接收时间在30 min以内.卫星信息同实验一.

初始电子密度和反演电子密度分别与测高仪数据对比后的误差百分比为：9.92%和8.07%.初始值与反演值的反演残差分别为ΔTEC₀=350.53×10¹⁶el·m^-2和ΔTEC₁=297.03×10¹⁶el·m^-2.

因此，通过卡尔曼滤波模型过滤得到的反演电子密度误差更小更接近真实值.但是，这里改善的不很明显，只有不到2%，但这不能表示卡尔曼滤波法仿真效果不太好.由图 8的差图可以看到，在纬度区间18°—28°之间反演电子密度的改善最为明显，在5%以上，尤其是20°—22°，甚至可以达到或超过10%，而测高仪站点的纬度是在31°左右，此处的改善不很明显.

图 8

Fig. 8

图 8 实验二电离层卡尔曼滤波反演图 Fig. 8 Ionospheric Kalman filter inversion figure of experiment 2

三频信标技术是近几年新兴的空间探测技术，其探测技巧与手段在某些方面还有待完善，但是凭借着其相对双频信标探测的巨大优势，已经在空间探测领域占据越来越重要的位置.

本次实验在三频卫星信标技术的基础上，应用卡尔曼滤波模型对初始场的电子密度进行反演，得到误差小、可靠性高的反演电子密度.但是，由于客观条件的限制，在本次实验中只有苏州一处测高仪数据，很多结论还有待完善.但是，通过本次实验知道，只要对输入条件设置恰当，应用卡尔曼滤波模型可以很好的完成电子密度反演工作，降低误差.因此，我们有理由相信，卡尔曼滤波法在电离层研究领域必会大有作为.

致谢 中国科学院上海天文台的全体老师学生对本次研究工作带来的巨大帮助，上海佘山接收站、温州接收站及台湾各接收站为本次研究提供的三频数据，苏州接收站为本次实验提供测高仪数据，特别是中国科学院上海天文台的郭鹏博士为本次研究提供的宝贵意见，一并表示衷心感谢.