电离层是位于地面上空约60~1 000 km的大气区域，其中部分中性粒子被电离为自由电子和离子^[1]。因此，电离层是一种色散介质，能引起无线电信号的群延迟和相位超前，影响程度与信号频率及电离层电子密度有关。为了评估电离层对无线电信号传播的时延影响，研究人员提出了一些描述电离层电子密度的理论模型或经验模型^[2-3]。理论模型虽然能反映电离层变化的物理机制，但精度无法满足电离层延迟改正的要求。经验模型和真实电离层之间也存在无法忽略的偏差，只能改正部分电离层延迟。因此，通过观测对电离层进行建模已经成为一种必要手段。一些观测仪器例如测高仪和非相干散射雷达已经应用于探测电离层的电子密度，然而这些仪器只能得到测站上空的电离层剖面，无法进行大规模的电离层观测。随着全球导航卫星系统的运行，电离层观测可以通过全球导航卫星系统多频无线电信号进行。由于全球导航卫星系统是基于倾斜轨道卫星星座连续运行的，不同时空尺度的电离层变化都可以通过全球导航卫星系统进行探测。

本文利用全球导航卫星系统，电离层电子密度的积分测量值即总电子含量作为参数指标进行电离层反演。不同频率无线电信号通过卫星和接收机硬件通道时产生的硬件偏差作为系统误差必须剔除，全球导航卫星系统多站观测网通常用于电离层硬件偏差的估算和电离层参数的反演。尽管如此，单站电离层反演作为一种灵活简便又成本低的方法仍然具有研究价值。单站电离层反演的算法主要有网格法、卡尔曼滤波法和多项式法^[4]。多项式法可以是二次多项式、三角级数或球谐函数^[5]，一般来说，二次多项式和三角级数只适合较小范围的反演，球谐函数通过改变阶数和次数能更好地描述电离层在不同空间尺度的变化。

电离层反演的拟合方法通常有最小二乘法和卡尔曼滤波法。人工神经网络提出之后不久就应用于电离层反演^[6]。基于误差反馈的前向传播神经网络能够以任意精度逼近一个连续的实函数。这种拟合通过迭代的方式进行，各参数对应的权值可以动态调整。因此这种动态拟合适应于电离层的准实时反演。

基于人工神经网络，我们提出一种准实时的单站电离层反演方法。首先对电离层反演原理和提出的方法进行了推导。接下来利用位于北纬30.28°、东经109.46°的单站在2014年春分、夏至、秋分和冬至的全球导航卫星系统数据进行了实验，其中，位于东经110°子午线的18个站作为验证站。在另一个实验中，2016年10月13日电离层暴事件的数据用于电离层反演，并且将反演结果与前后磁静日的结果进行了比较。

1 电离层反演算法 1.1 球谐函数拟合

全球导航卫星系统以两个不同的L波段频率播发信号，双频接收机接收的伪距和载波相位观测值为

$ \begin{gathered} P_{i}=\rho+c\left(\delta_{\mathrm{t}}-\delta_{\mathrm{s}}\right)-I_{i}+T+\varepsilon, \end{gathered} $

(1)

$\frac{c}{f_{i}} L_{i}=\rho+c\left(\delta_{\mathrm{t}}-\delta_{\mathrm{s}}\right)+I_{i}+T+\frac{c}{f_{i}} N_{i}+\varepsilon, $

(2)

其中，i为1或2，分别为f₁或f₂频率；P_i和L_i分别为f_i频率的伪距和载波相位观测值；ρ为卫星和接收机间的真实距离；c为光速；δ_t和δ_s分别为接收机和卫星钟差，并且伪距包含群延迟，载波相位包含相位超前；I_i为f_i频率的电离层延迟；T为对流层延迟；N_i为L_i的整周模糊度；ε为观测值的残差项。

由不同频率观测值的差分，沿着传播路径的斜向总电子含量(Slant TEC, STEC)可以分别通过伪距和载波相位观测值求出：

$ \begin{gathered} S_{\mathrm{TEC}_{\mathrm{P}}}=\frac{\left(f_{1} f_{2}\right)^{2}}{40.3\left(f_{1}^{2}-f_{2}^{2}\right)}\left(P_{2}-P_{1}\right), \end{gathered} $

(3)

$ S_{\mathrm{TEC}_{\mathrm{L}}}=\frac{\left(f_{1} f_{2}\right)^{2}}{40.3\left(f_{1}^{2}-f_{2}^{2}\right)}\left(\frac{c L_{1}}{f_{1}}-\frac{c L_{2}}{f_{2}}+\frac{c N_{1}}{f_{1}}-\frac{c N_{2}}{f_{2}}\right) . $

(4)

考虑到伪距包含较大的观测噪声，且载波相位含有未知的整周模糊度，高精度斜向总电子含量可以通过载波相位平滑伪距求得^[7]

$ S_{\mathrm{TEC}}(t)=\omega_{\mathrm{t}} S_{\mathrm{TEC}_{\mathrm{P}}}(t)+\left(1-\omega_{\mathrm{t}}\right)\left[S_{\mathrm{TEC}}(t-1)+S_{\mathrm{TEC}_{\mathrm{L}}}(t)-S_{\mathrm{TEC}_{\mathrm{L}}}(t-1)\right] \quad(t>1), $

(5)

其中，t为历元；ω_t为权重，与卫星高度角有关。对于第1个历元，S_TEC与S_TECP相等。

电离层反演即使用斜向总电子含量得到反演区域内沿垂直方向的总电子含量。通常假设电离层的所有电子集中于一个无限薄的单层上，为此引入了如图 1的电离层单层模型^[8]，其中点S，R和O分别表示卫星、接收机和地心。卫星和接收机视线与电离层单层的交点定义为穿刺点(Ionospheric Piercing Point, IPP)，其经纬度等于对应的星下点(SubIonospheric Point, SIP)的经纬度。因此总电子含量可以根据图 1的几何关系得

$ \begin{aligned} T E C &=\left(S_{\mathrm{TEC}}-b_{\mathrm{s}}-b_{\mathrm{r}}\right) \cos \chi, \end{aligned} $

(6)

$ \chi =\arcsin \left(\frac{R_{\mathrm{e}} \cos E}{R_{\mathrm{e}}+H}\right), $

(7)

其中，b_s和b_r分别为卫星和接收机的硬件偏差；χ和E分别为天顶角和高度角；R_e为地球的平均半径；H为电离层单层的高度。

另一方面，总电子含量在模型中可以认为是一个球谐函数，参数为穿刺点对应的余纬和日地坐标系下的经度^[5]：

$ T E C=\sum\limits_{n=0}^{N} \sum\limits_{m=0}^{n}\left(A_{n m} \cos m \phi+B_{n m} \sin m \phi\right) P_{n m} \cos \theta, $

(8)

其中，θ和ϕ分别为余纬和日地坐标系下的经度；对于单站，通常取阶数和次数为5(n=m=5)；P_nm为归一化的勒让德函数；A_nm和B_nm为球谐函数的系数，可以通过最小二乘法或其他拟合方法求得。

1.2 人工神经网络

人工神经网络由许多神经元组成，这些神经元组成了一个输入层、一个输出层和若干隐藏层，通常一个隐藏层是足够的。一个基于误差反馈的多层前向传播神经网络能够以任意的精度趋近于一个连续实函数^[6]。图 2为用于准实时电离层反演的人工神经网络示意图，其中，k表示卫星和接收机的数目。由于在电离层反演时，卫星和接收机的硬件偏差需要同时估计并剔除，额外的网络用于硬件偏差的估计。输入参数为斜向总电子含量对应的余纬、日地坐标系下的经度以及高度角。输出参数为(8)式的球谐函数的系数以及硬件偏差，观测站上空的总电子含量可以通过球谐函数的系数求得。输入层的神经元除硬件偏差的初值外，还有穿刺点的余纬和日地坐标系下的经度以及天顶角。输出层包括硬件偏差以及球谐函数系数。隐藏层的节点数是一个小于200的随机值。输出层的值通过动态权值进行调整。权值通过梯度下降算法根据观测值与模型估计值的残差求出^[9]，

$ E=\frac{1}{2}\left[V_{\mathrm{TEC}} / \cos \chi-\left(S_{\mathrm{TEC}_{\mathrm{obs}}}-D C B\right)\right]^{2}, $

(9)

其中，V_TEC为电离层穿刺点处的垂直向总电子含量。

根据(9)式，以每15 min内各卫星对应的斜向总电子含量为训练数据集，在对神经网络进行训练的过程中不断调整权值和输出层节点，最终使(9)式收敛。由于电离层的变化具有连续性，下一个时段可用当前时段训练的权值和输出层节点作为初始值，并继续用新的数据集进行训练，以达到快速收敛的目的。对于第一时段，前一天的硬件偏差作为初值，之后每个时段都用前一个时段的硬件偏差作为初值。输入值经过非线性变换进入输入层再经过隐藏单元逐层处理，传向输出层，输出层判断估值与观测值的残差，大于阈值则转入反向传播，通过梯度下降算法修改各神经元的权值，直至输入满足条件。由于不同卫星对应的观测值包含非均一的观测误差^[10-11]，而且电离层受到的扰动也随着时间动态变化，因此，神经网络是一种良好的准实时估计算法。

我们基于图 2的神经网络架构使用MATLAB语言编写了准实时单站电离层反演程序。该程序每15 min通过实时接收的全球导航卫星系统数据进行一次电离层反演，得到一组球谐函数系数，并更新硬件偏差。使用3.4 GHz的中央处理器和16 GB内存的工作站测试后证实每一次电离层反演可以在1 min左右完成，因此可以实现准实时的单站电离层反演。

2 磁静日数据验证 2.1 测试观测网和数据处理方法

如图 3为全球导航卫星系统观测站的分布，其中位于北纬30.28°、东经109.46°的观测站为测试站，其余位于东经110°子午线附近的18个站为验证站。为了对提出的准实时神经网络法(Real-Time ANN, RTANN)进行详细的评估，我们使用测试站的全球定位系统(Global Positioning System, GPS)的双频观测数据分别用准实时神经网络法和最小二乘球谐函数法反演电离层总电子含量。同时，使用18个验证站的全球定位系统双频观测数据，我们利用1°×1°的自适应网格法估计硬件偏差^[12]，再利用2°×2°的网格法每15 min得到电离层总电子含量^[13]，并将它们作为参考值与准实时神经网络法和最小二乘球谐函数法得到的结果进行比较。自适应网格法已经在文^[12]中应用于电离层暴时的硬件偏差估计，因此本文使用该方法得到硬件偏差参考值。除此之外，尽管多站网的反演结果也含有误差，但在观测站分布合理的情况下，多站网的反演结果比单站有更高的精度和更好的稳定性，因此作为参考值。在实验中，我们选择2014年春分、夏至、秋分和冬至的观测数据，这4天均为磁静日。历元间隔为30 s，截止高度角选取20°，电离层单层高度设为400 km。

2.2 估计的硬件偏差

由于在电离层反演中，卫星和接收机的硬件偏差作为整体耦合在一起，我们采用零均值处理分离卫星和接收机硬件偏差，即对每一颗卫星硬件偏差减去它们的均值，将均值并入接收机硬件偏差。图 4为通过准实时神经网络法和最小二乘球谐函数法估计的卫星硬件偏差与参考值的差值，(a)~(d)分别为春分、夏至、秋分和冬至的结果。在这4天中，通过准实时神经网络法得到的卫星硬件偏差与参考值的差值在整体上比最小二乘球谐函数法更小，原因在于神经网络的引入使得硬件偏差可以动态调整，减小了观测噪声以及电离层扰动对反演结果的影响。另外有极少数卫星的硬件偏差与参考值的偏差较大，这是由零均值处理造成的。还有部分卫星在测试时处于不可用状态，因此其结果在图中没有显示。经过计算，通过准实时神经网络法估计的硬件偏差在4天的均方根误差分别为0.38 ns，0.25 ns，0.37 ns和0.31 ns。而最小二乘球谐函数法的结果为0.84 ns，0.61 ns，0.67 ns和0.68 ns。春分时精度提高最为明显，这与春分时电离层总电子含量较大，对硬件偏差估计的噪声干扰较大有关，这将在下文进一步分析。

2.3 反演的总电子含量

图 5为通过准实时神经网络法和最小二乘球谐函数法反演的北纬30.28°、东经109.46°单站的电离层总电子含量，(a)~(d)分别为春分、夏至、秋分和冬至的结果。参考值是通过18个站采用网格法得到的，并且每15 min独立反演。我们可以发现图 5(a)中春分时的总电子含量明显大于另外3天，因此, 通过(6)式与总电子含量同时估计的硬件偏差受噪声影响也更大。将单站反演的总电子含量与参考值对比后发现，通过准实时神经网络法反演的总电子含量比最小二乘球谐函数法更接近于参考值，尤其是正午和午夜。正午提高较明显是因为正午时总电子含量较大，而神经网络的加入使得反演受噪声影响更小。而午夜是因为此时总电子含量少，最小二乘球谐函数法反演的结果受硬件偏差的影响较大，而准实时神经网络法估计的硬件偏差更准确。另外，准实时神经网络反演的总电子含量也不平滑，我们推测这可能与电离层扰动有关，因为大部分波动趋势与参考值相对应。经计算，准实时神经网络法的均方根误差分别为2.76，1.00，0.99和1.73 TECU(TEC Unit)，最小二乘球谐函数法的为3.46，2.32，0.94和1.46 TECU。因此，通过准实时神经网络法反演的总电子含量比最小二乘球谐函数法更准确。

3 电离层暴数据测试

我们分别用准实时神经网络法和最小二乘球谐函数法通过单站反演了2016年10月13日电离层暴的总电子含量。图 6为0°子午线和40°N纬线的全球导航卫星系统观测站的分布，其中位于北纬40.52°、西经3.09°的观测站为测试站，其余位于0°子午线和40°N纬线附近的7个站为验证站。根据文^[14]，暴的起始发生于12日2100 UT，初相从2300 UT持续至13日的0600 UT，接下来主相持续了18 h，并且图 6的7个验证站观测到13日白天的正暴及夜间的负暴，结合测高仪观测到的电离层抬升证实了引起这个电离层暴的主要物理机制是扰动中性风。

如图 7为不同方法估计的2016年10月8日至17日的卫星硬件偏差相对于10天日均值的均方根误差。由于10月13日发生了电离层暴，估计的硬件偏差相对于磁静日变化较大。但实际上硬件偏差取决于仪器本身以及所处的外部环境，只是由于反演过程中与总电子含量耦合在一起才受到电离层变化的影响。而准实时神经网络法的均方根误差在这10天中均小于最小二乘球谐函数法，特别是在13日暴发生期间。这说明准实时神经网络法估计的硬件偏差受电离层变化的影响更小，因此优于最小二乘球谐函数法。

如图 8为不同方法反演2016年10月8日至17日北纬40.52°、西经3.09°的总电子含量。相比8~11日磁静日的总电子含量，13日白天的总电子含量表现为明显的正暴效应，总电子含量增大率超过了100%，而从13日夜间开始表现为明显的负暴效应，并持续至15日。另外，反演的总电子含量与文^[14]的结果相符，并且与沿赤道方向的扰动中性风作用于电离层电子密度的物理机制相符。暴时的赤道向中性风可以抬升电离层，同时也会引起电离层的成分变化。白天由于电离层较低高度有光电离作用的补充，正暴发生，而夜间光电离作用消失同时离子结合率上升，负暴发生。另外，比较准实时神经网络法和最小二乘球谐函数法的反演结果发现，两种方法反演的总电子含量具有高度一致性。另外，准实时神经网络法能更好地反映电离层的变化，主要表现在正午和午夜，而夜间准实时神经网络法的总电子含量较低也与图 5的结果相符合。

4 结语

本文提出了一种基于人工神经网络的准实时单站电离层反演方法。在设计的神经网络中，输入层包含作为初值的前一天的硬件偏差以及穿刺点的余纬、日地坐标系下的经度以及高度角，输出层为球谐函数系数以及硬件偏差。通过不断调整权值，输出层可以利用梯度下降算法根据残差进行迭代。因此硬件偏差可以随着观测值动态调整并且总电子含量也可以准实时反演。实验分别采用本文方法和经典的最小二乘球谐函数法对位于北纬30.28°、东经109.46°的单站在2014年春分、夏至、秋分和冬至的磁静日全球导航卫星系统数据进行电离层反演，位于东经110°子午线的18个站作为参考站，通过网格法获得硬件偏差和总电子含量参考值。结果显示使用本文方法估计的硬件偏差整体上更接近参考值，反演的总电子含量相对于最小二乘球谐函数法更接近参考值。在另一个实验中，2016年10月13日电离层暴事件的总电子含量分别通过位于北纬40.52°、西经3.09°的单站使用本文方法和最小二乘球谐函数法反演。两种方法反演的总电子含量也具有高度一致性，基于神经网络的硬件偏差估计与磁静日的硬件偏差估计更接近。这些结果说明准实时神经网络方法相比最小二乘球谐函数法具有较高的精度。神经网络同样适用于准实时多站电离层反演，值得进一步推广。

致谢: 感谢IGS网站及中国陆态网提供全球导航卫星系统实测数据。