电离层TEC，均值插值法，球谐函数模型

Empirical Model of Polar Ionospheric TEC Based on Spherical Harmonic Function

摘要: 本文参照极区电离层TEC的UT变化特征，统计分析了极区网格划分方法。引入网格均值插值方法，结合球谐函数理论，构建了基于球谐函数的极区电离层TEC经验模型，降低了拟合出现负值的情况。通过拟合太阳活动低年2007~2010年的极区电离层TEC观测数据，表明极区平均电离层TEC值随UT呈“波浪”式变化。通过广义三角级数函数模型与球谐函数模型的对比，发现后者对极区电离层TEC建模具有更好的适用性。

Abstract: Based on the characteristics of ionospheric TEC UT variation in polar region, the grid division of polar region is analyzed statistically. An empirical model of ionospheric TEC in polar region is proposed by introducing grid mean interpolation method to improve the spherical harmonic function. By fitting the ionospheric TEC data of low solar activity years from 2007 to 2010 using the improved the spherical harmonic function model, it is shown that the average ionospheric TEC value of the polar region changes with UT in a “wave” pattern. The improved spherical harmonic function model is also compared with the generalized trigonometric series function model. This is better applicability for the spherical harmonic function model to the polar ionospheric TEC modeling.

1. 引言

电离层是地球大气上层的一个电离区域，范围大致在地表以上高度60 km至1000 km处，极区作为地球磁力线汇聚的区域，是当前电离层研究的热点 [1]。宏观上电离层呈现电中性状态，但微观上呈带电状态 [2]，电离层中的等离子体在太阳辐射和大气碰撞的作用下电离产生大量的自由电子和离子，它们影响着经过其中的无线电波传播 [3]。电离层总电子含量(Total Electron Content, TEC)是表征电离层特性的重要参数 [4]，其变化与多种因素 [5] [6] 有关：地方时、地理位置、世界时、太阳活动等。TEC的空间分布及时间变化反映了电离层的主要特性，因此通过探测与分析极区电离层TEC参量，可以研究电离层不同时空尺度的分布与变化特性 [7]，对全球卫星导航系统(Global Navigation Satellite System, GNSS)导航和定位、地震 [8] 预报以及通讯具有重要的科学研究价值。

相关学者提出了许多基于GPS观测数据的电离层模型，球谐函数模型 [7] 是目前常用的区域电离层经验模型。早在1998年Schaer通过球谐系数展开，构建了8阶全球电离层球谐函数模型，实现了2小时分辨率的全球电离层地图构建 [8]。Wilson等人比较了球谐函数和三角网格内插法所建立模型的差异 [9]，指出球谐函数模型具有较好的预测效果。Li等人提出了一种球谐函数加广义三角级数函数的电离层建模方法 [10]，提高了模型的稳定性，中国科学院(Chinese Academy of Sciences, CAS/IGG)采用广义三角级数函数建立单站的局部模型，并结合球谐函数计算了全球模型 [11]。欧洲定轨中心(Center of Orbit Determination in Europe, CODE)将全球空间分辨率划分为2.5˚纬度和5˚经度 [12]，开发了15阶球谐函数模型。然而实测数据不足以覆盖整个极区，因此对极区电离层TEC建模的网格划分需要同时兼顾网格分辨率、覆盖范围和网格观测频数。

电离层世界时(Universal Time, UT)反映了电离层参数在某一固定纬度不同经度的均值的日变化，是影响电离层特征的重要参数 [13]。本文参照极区电离层TEC的UT变化特征，统计分析了极区网格划分。结合网格均值插值法和球谐函数理论，构建了基于球谐函数的极区电离层TEC经验模型，降低了球谐函数拟合出现负值的情况。为分析该模型对极区电离层TEC的拟合效果，本文利用太阳活动低年2007~2010年极区电离层TEC观测数据，综合比较了广义三角级数函数模型和球谐函数模型在时空分布上的差异，结果表明后者对极区电离层TEC建模有较好的效果。

2. 数据处理

2.1. 观测数据

研究中采用的GPS数据由MIT Haystack Observatory Madrigal database提供，Madrigal是一个上层大气科学数据库，能够以各种格式提供来自各种高层科学仪器的档案和实时数据。我们处理分析了太阳活动低年2007~2010年的TEC观测数据，数据采样间隔为5分钟。通过网格划分，全球电离层TEC地图如图1所示。

Figure 1. Global ionospheric TEC map at 0 UT on 1 March 2007

图1. 2007年3月1日0 UT全球电离层TEC地图

图1表明，全球电离层TEC观测在陆地区域覆盖比较全面，观测空白主要集中在低纬度海洋区域。高纬度区域观测覆盖全面，TEC值大致在5~10 TECu之间，符合研究要求，其他时间全球电离层TEC地图分布与之基本相同，仅TEC值有些变化。

2.2. 网格划分

CODE将全球空间分辨率划分为2.5˚纬度和5˚经度，然而实测数据不足以覆盖整个极区。这里根据电离层UT变化特征，构建地磁坐标系(MLat-MLT)。根据MIT Haystack Observatory Madrigal database提供的数据进行时空划分 [14]：1) 时间划分：为排除穿刺点不足对研究的影响，以UT (UT−0.5, UT+0.5)为划分单位。2) 空间划分：考虑高分辨率网格有助于提高模型的适用范围和精度，按照如下统计方法，获得空间网格划分：a) 网格面积尽可能小，b) 网格观测分布空白尽可能少，c) 单位网格内有足够多的观测频数(不少于10个)。

通过时空划分方法确定网格划分，每个网格以中心点的坐标表示，划分结果如表1所示：

Table 1. Statistical table of polar region grid division

表1. 极区网格划分统计表

磁地方时(MLT)由磁经度(MLon)转化获得。0.33 MLT对应的磁经度约是5˚，0.67 MLT对应的磁经度约是10˚，1.00 MLT对应的磁经度约是15˚。划分结果网格观测分布空白少，网格面积较小，网格内有足够多的观测频数。

以2007年3月份09:00 UT网格划分为例，如图2所示，极区网格观测空白较少，观测覆盖面较全；网格内观测频数总体高于10个，观测频数最大超过300个。网格划分情况较好，符合研究要求，其他时间分布情况与之相同。

为避免极端值对研究的影响，我们统计计算了极区网格TEC中值，如图3所示，网格TEC中值最大超过10 TECu，最小值为1.1 TECu，12:00 MLT附近TEC较大，00:00 MLT附近TEC较小，其他UT时刻的TEC随地磁纬度和磁地方时的分布情况与之基本一致。

Figure 2. Ionospheric TEC observation frequency distribution in UT polar region at 09:00 March 2007

图2. 2007年3月份09:00 UT极区电离层TEC观测频数分布图

Figure 3. Ionospheric TEC distribution in polar region at 09:00 UT, March 2007

图3. 2007年3月份09:00 UT极区电离层TEC中值分布图

3. 基于球谐函数的极区电离层TEC经验模型

3.1. 球谐函数

球谐函数在建立全球电离层模型中得到广泛的应用。1998年Schaer [8] 通过球谐系数展开，构建8阶全球电离层球谐函数模型，实现了2小时分辨率的全球电离层地图构建(Global Ionospheric Map, GIM)，且具有较高的精度。CODE利用15阶球谐函数模型实现了每小时一次的GIM绘制。在实际建模中，球谐函数零阶项 [15] 表征的是区域内的平均电离层TEC值。球谐函数 [16] 表达式如下：

$\text{VTEC}={\displaystyle \underset{n=0}{\overset{n_{\max }}{\sum }}{\displaystyle \underset{m=0}{\overset{n}{\sum }}{\tilde{P}}_{nm}\left(\sin \theta \right)}}\left(C_{nm}\cos m\phi +D_{nm}\sin m\phi \right)$ (3-1)

式中 ${\tilde{p}}_{nm}\left(\sin \theta \right)$ 是勒让德多项式； $\theta$ 为穿刺点的磁纬度， $\phi$ 为电离层穿刺点处的经度与太阳直射点的经度之间的夹角； $n_{\max }$ 是球谐函数最大展开阶数； $C_m$ ， $D_m$ 为球谐函数模型的系数，即模型所求的参数。研究中 $n_{\max }=15$ ，这意味着每UT有256个球谐系数。

3.2. 网格均值插值法

在实际应用中极区部分区域出现了拟合值为负的问题，这是由于这些区域缺少观测值造成的，对此我们需要对球谐函数的输入进行处理。本文结合网格均值插值方法，利用缺失值所在网格的“九宫格”TEC平均值填充缺失值，网格均值差值公式如下所示：

${\mathrm{TE}{C}^{\prime }}_i=\frac{1}{N}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left({\text{TEC}}_i\right)}$ (3-2)

式中 ${\text{TEC}}_i$ 为观测空白网格的TEC值， ${\mathrm{TE}{C}^{\prime }}_i$ 为均值差值方法填充后的网格TEC值。在实际应用中，我们采用“九宫格”均值插值法，N为以该网格为中心的相邻8个网格中有观测值的网格数量。

4. 模型拟合结果

我们将填充后的极区电离层网格TEC中值代入球谐函数(3-1)，推算每UT 256个球谐系数，即可得到极区网格内的TEC拟合值。

4.1. 球谐函数模型拟合结果

拟合结果如图4所示，拟合情况与极区电离层TEC中值分布图(图3)接近，拟合情况较好，部分观测空白的区域也得到了较为平滑的拟合，其他UT拟合结果相似。

4.2. 极区平均电离层TEC值

极区平均电离层TEC值可由球谐系数的零阶项表示，图5(a)~(d)表示2007~2010年球谐函数零阶项分布情况。极区平均电离层TEC值在18:00~20:00 UT达到最大值，最大超过5 TECu，在06:00~08:00 UT达到最小值，最小接近2 TECu。极区平均电离层TEC值随UT变化特征显著，呈现“波浪”式变化。在00:00~07:00 UT和19:00~23:00 UT，极区平均电离层TEC值随UT的增大而减小；然而在07:00~19:00 UT，随UT的增大而增大。与武业文 [2] 提出的北极电离层平均电子含量(mPEC)随UT呈三角函数的“正弦”形状变化相似。

其次，极区平均电离层TEC值与月份也有显著的关系，我们发现在12~23 UT内，4、5、6月，即第二季度的极区平均电离层TEC值随UT的变化波动幅度较小；而10、11、12月，即第四季度的极区平均电离层TEC值随UT变化波动幅度显著。

Figure 4. Ionospheric TEC fitting distribution in polar region, 09:00 UT March 2007

图4. 2007年3月份09:00 UT极区电离层TEC拟合分布图

Figure 5. The average ionospheric TEC value in polar region from 2007 to 2010 varies with UT

图5. 2007~2010年极区平均电离层TEC值随UT变化

5. 模型精度分析

为准确评定极区电离层球谐函数模型的精度，我们用平均绝对误差(MAE)和均方根误差(RMSE)来验证构建的模型在拟合极区电离层TEC的效果，误差计算公式如下：

$\text{MAE}=\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}\left|{\text{TEC}}_j\left(obs.\right)-{\text{TEC}}_j\left(pre.\right)\right|}$ (5-1)

$\text{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \underset{j=1}{\overset{N}{\sum }}{\left[{\text{TEC}}_j\left(obs.\right)-{\text{TEC}}_j\left(pre.\right)\right]}^2}}$ (5-2)

式中 ${\text{TEC}}_j\left(obs.\right)$ 和 ${\text{TEC}}_j\left(pre.\right)$ 分别是某网格的观测值和拟合值。通过拟合太阳活动低年2007~2010年极区电离层TEC观测数据，统计分析球谐函数MAE和RMSE结果如下：

Figure 6. Error of spherical harmonic function model from 2007 to 2010

图6. 2007~2010年球谐函数模型误差结果

图6(a)~(d)表明2007~2010年球谐函数模型的平均绝对误差(蓝色虚线线条)以及均方根误差(红色实线线条)情况，横坐标为每月UT，左侧纵坐标为平均绝对误差，右侧纵坐标为均方根误差，单位均为TECu。2007~2010年平均绝对误差均小于0.9 TECu，平均绝对误差随UT变化波动总体较小，部分时间段波动显著，如：2007年2月0 UT~3月0 UT。2007年均方根误差小于2 TECu，2008~2010年均方根误差均小于1 TECu，且波动随UT变化总体较小，部分时间段均方根误差波动比较剧烈，如：2007年2月0 UT~3月0 UT、2010年2月0 UT~3月0UT有显著的变化，其他时间均方根误差波动幅度较小。误差分析结果表明球谐函数模型拟合值与极区电离层TEC观测值基本一致，且比较稳定。

我们对比了球谐函数模型和广义三角级数函数模型 [17] 对极区电离层TEC的拟合精度。其中广义三角级数函数模型 [18] [19] 中 $N_2=0$ ， $N_i=N_j=1$ ， $N_3=2$ ， $N_4=6$ ，通过求解模型系数，反推2007~2010年极区电离层TEC拟合值，得到模型误差对比如下：

如图7所示为两个模型的极区电离层TEC平均绝对误差折线图。通过精度指标变化的情况，可以发现球谐函数模型和广义三角级数函数模型有着相似的变化趋势。球谐函数模型平均绝对误差较小，总体平均绝对误差约为0.4 TECu，2007年2月0 UT~4月23 UT的平均绝对误差较大，最大超过为0.8 TECu。广义三角级数函数模型平均绝对误差较大，总体平均绝对误差约为0.9 TECu，2007年2月0 UT~7月23 UT和2010年3月0 UT~11月23 UT的平均绝对误差较大，最大接近1.3 TECu，且波动比较剧烈。由平均绝对误差可以直观的看到，球谐函数模型拟合效果优于广义三角级数函数模型。

如图8所示为两个模型的极区电离层TEC均方根误差折线图。球谐函数模型均方根误差较小，总体均方根误差约为0.53 TECu，2007年3月0 UT~23UT的均方根误差较大，最大超过为1.2 TECu。广义三角级数函数模型均方根误差较大，且波动剧烈。总体均方根误差约为1.25 TECu，2007年2月0 UT~7月23 UT和2010年3月0 UT~11月23 UT的均方根误差较大，最大接近1.8 TECu。结合平均绝对误差和均方根误差，可得球谐函数模型对极区电离层TEC建模具有优越性。

Figure 7. Mean absolute ionospheric TEC errors of the two models in the polar region

图7. 两个模型极区电离层TEC平均绝对误差折线图

Figure 8. Root mean square error of ionospheric TEC of the two models in polar region

图8. 两个模型极区电离层TEC均方根误差折线图

6. 结论

本文参照极区电离层TEC的UT变化特征，统计分析了极区网格划分。引入网格均值插值方法，构建了基于球谐函数的极区电离层TEC经验模型。通过拟合太阳活动低年2007~2010年的极区电离层TEC观测数据，利用平均绝对误差分析和均方根误差分析方法，验证了模型精度，得到如下结论：

1) 结合网格插值方法的球谐函数模型能够实现极区电离层TEC建模，且精度较高，平均绝对误差约为0.37 TECu，均方根误差为0.56 TECu。

2) 极区平均电离层TEC值随UT变化特征显著，呈现“波浪”式变化。在00:00~07:00 UT和19:00~ 23:00 UT，极区平均电离层TEC值随UT的增大而减小；而在07:00~19:00 UT，随UT的增大而增大。其次，第二季度的极区平均电离层TEC值波动最小，第四季度的极区平均电离层TEC值波动最大

3) 通过广义三角级数函数模型和球谐函数模型精度对比，后者对极区电离层TEC有更好的拟合效果。

对于更精细的模型比较和电离层分析，需要更长的时间尺度和更均匀分布的观测，还需进一步考虑太阳活动、行星际磁场条件下的电离层数据特征。

NOTES

^*通讯作者。