研究表明，在中高轨SAR成像^[1]、地质研究^[2-3]等领域，都需要对电离层电子密度进行准确反演，其中一个重要分支就是利用GPS卫星进行电离层层析成像。1988年，Austen^[4]借鉴CT成像技术提出电离层层析(CIT)，从此开辟了人们探究电离层的新方法^[5-11]。以代数重构算法为代表的电离层层析方法是研究电离层层析反演的主流方法之一，其原理是将初始值向各观测方程进行投影，投影过程中将残差分配到各个网格的电子密度中，通过不断迭代修正得到更接近于真实值的电子密度。Pryse等^[12]提出的MART(multiplicative algebraic reconstruction technique)算法则是每次迭代时将STEC实测值与STEC反演值的比值作为电子密度的缩放因子分配到有信号穿过的电离层网格内，通过乘法和乘方运算进行电子密度的修正。在后续的发展过程中, 霍星亮等^[13]通过在不同电子密度的网格中重新分配STEC实测值与STEC反演值的残差改进了代数重构方法，但依然只有被信号穿过的网格才会逐渐逼近真值。本文借鉴代数重构算法和乘法代数重构算法，将STEC实测值与STEC反演值的比值(基准尺度因子)分配到反演空间内的所有网格上。对于有信号穿过的网格，其缩放尺度即为基准尺度因子；由于电子密度在空间中具有相关性，距离较近的网格内电子密度一般差别较小，距离较远的网格可能差别较大，所以对于没有信号穿过的网格，通过计算其网格中心与观测射线的距离，结合基准尺度因子确定其对应的尺度因子，在每次迭代过程中利用尺度因子带动电子密度逐渐向真值“靠近”。

1 代数重构方法

电离层层析常用的一种方法是根据GPS卫星的观测延迟计算斜向电子总量(slant total electron content，STEC)，进行一定区域内的电子密度重构。STEC可以表示为GPS卫星到接收机之间的斜向路径电子密度积分值：

$ {\rm{STEC}} = \int_l {{N_e}(s){\rm{d}}s} $

(1)

式中，l表示观测射线，在反演中当观测仰角足够大时将观测射线假设为直线；N_e(s)表示空间某处的电子密度，ds表示GPS卫星到接收机的斜向距离元的长度。为了便于计算和进行重构，通常将电离层按照经纬度划分为均匀网格。图 1为电离层像素模型和观测几何示意图，图中d_j表示第j个网格中心到测量射线的距离。

将式(1)改写成离散形式：

$ {\rm{STE}}{{\rm{C}}_i} = \sum\limits_{j = 1}^n {{A_{ij}}{x_j} + {e_i}, i = 1, 2, ..., m} $

(2)

式中，x_j表示第j个网格的电子密度。写成矩阵形式为：

$ {\mathit{\boldsymbol{y}}_{m \times 1}} = {\mathit{\boldsymbol{A}}_{m \times n}}{\mathit{\boldsymbol{x}}_{n \times 1}} + {\mathit{\boldsymbol{e}}_{m \times 1}} $

(3)

式中，y表示m次测量得到的STEC排列成的列向量，x_n×1表示将n个三维空间网格按照特定顺序展开成的列向量，A_ij表示第i次测量在第j个网格上的截距长度。通常，由于地面接收机数量、卫星轨道以及观测仰角的限制，可利用的观测量很少，因此需要对上式病态方程进行优化求解。

Pryse等^[12]提出代数重构(ART)算法：

$ x_j^{k + 1} = x_j^k + {\lambda _k}\frac{{{y_i} - \sum\limits_{j = 1}^n {{A_{ij}}x_j^k} }}{{{{\left\| {{A_i}} \right\|}^2}}}{A_i} $

(4)

式中，i=kmodm+1，k为迭代次数，每m次迭代称为一轮迭代, A_i表示截距矩阵的第i行。1990年，Raymund等^[14]改进ART得到新算法——MART算法并运用于电离层层析中。乘法代数重构迭代式如下：

$ x_j^{k + 1} = x_j^k{\left( {{y_i}/\sum\limits_{j = 1}^n {{A_{ij}}x_j^k} } \right)^{\frac{{{\lambda _k}{A_{ij}}}}{{{\rm{max}}{A_{ij}}}}}} $

(5)

2 改进重构方法

本文提出附加尺度因子的代数重构算法(Scaling-ART)。在每次迭代完成之前计算每个网格对应的尺度因子，距离观测射线越远的网格受到尺度缩放的影响越小，通过迭代计算尺度因子、不断修正电子密度值，使其逼近真值，得到更为准确的反演结果。具体表达式如下：

$ x_j^{k + 1} = \left\{ \begin{array}{l} {s_{ij}}^kx_j^k, {A_{ij}} = 0\\ x_j^k + {\lambda _k}\frac{{{y_i} - \sum\limits_{j = 1}^n {{A_{ij}}x_j^k} }}{{{{\left\| {{A_i}} \right\|}^2}}}{A_i}, {A_{ij}} \ne 0 \end{array} \right. $

(6)

假设第k次迭代时电子密度反演值为x^k，通过对应的截距向量A_i 得到观测射线上的基本尺度因子值s_i0^k，形如MART算法中的底数部分：

$ s_{i0}^k = \frac{{{y_i}}}{{\sum\limits_{j = 1}^n {{A_{ij}}x_j^k} }} $

(7)

将此尺度作为尺度基准，未被信号穿过的网格的尺度因子通过计算网格中心与该观测射线的距离确定：

$ s_{ij}^k = f({d_{ij}}) $

(8)

式中，s_ij^k表示第k次迭代时j号网格对应第i条观测射线的尺度因子，d_ij表示第j号网格到观测射线的距离，f(d_ij)表示d_ij的单调函数。将[mind_j, maxd_j]映射到[s_i0, 1], s_i0＜1或[1, s_i0], s_i0＞1范围内，对于未被观测射线穿过的网格，尺度因子若为1表示本次迭代不发生变化，从而有效防止STEC测量误差在整个反演区域的传播。本文所用的单调函数如下：

$ f({d_{ij}}) = {\rm{exp}}(a{d_{ij}} + b) $

(9)

在迭代时，根据上述映射区间和网格中心到观测射线的距离得到一组参数a^k和b^k，之后即可计算出所有网格的尺度因子。算法流程见图 2。

3 实验验证

为了分析本文算法在电离层电子密度反演中的重构精度，进行模拟实验和真实数据验证。两个实验选取的区域为95°~120°E、27°~42°N，反演高度为100~1 000 km。进行网格划分时，经度和纬度向间隔分别为1°，高度间隔为50 km。

3.1 数值模拟实验

在模拟实验中将IRI模型值作为电子密度“真值”进行实验。选取若干接收机和GPS卫星坐标进行截距矩阵的计算，此过程中需要避免截止高度角过大，以保证测量射线能够被假设为直线。接收机和电离层顶部穿刺点分布如图 3。

STEC通过“真值”数据和截距矩阵得到。考虑到STEC测量时存在误差，生成STEC时加入一定量的随机噪声，噪声在[-10%, 10%]内均匀分布。图 4是IRI模型在100°E子午面的电子密度分布情况，电子密度峰值在300~350 km高度，随纬度增加电子密度呈现减小的趋势。

尺度因子空间分布状况如图 5，红色线段表示观测射线，观测射线穿过的电离层底部和顶部两个网格对应的高度、纬度和经度分别为6 470 km、26°N、112°E及7 370 km、30°N、118°E。不同颜色的点表示对应位置网格的尺度因子，距离观测射线越远的网格其尺度因子呈现蓝色，其值越接近于1；距离越近的尺度因子在图中呈现黄色。根据上述公式，尺度因子为1，意味着此次迭代完成时该网格的电子密度值不会发生变化。

图 6所示的是在98°E处的子午面反演得到的绝对残差图，残差定义为真值与反演结果的差值。图 6(a)表示使用本文Scaling-ART方法进行反演的残差绝对值，图 6(b)表示利用MART算法进行电子密度重构得到的残差绝对值。图中显示，27°~34°N范围内，图 6(b)比图 6(a)红色区域更大，即表明使用MART反演得到的残差比Scaling-ART的残差更大，说明Scaling-ART算法所得结果优于MART算法。根据所选取的接收机位置和GPS测量时的坐标，由于电离层穿刺点主要集中在34°~36°N附近，而穿刺点一般在电子密度最大的高层，因此在此范围内反演的效果相对较好，在残差图中表现为该区域主要呈浅色；而此图中穿刺点分布很少的低纬度地区会存在较大误差，呈现偏红的状态，最大误差接近1.4×10¹⁰ el/m³。

图 7分别是在27°N和36°N两个纬度面上的反演结果，其中图 7(a)、图 7(b)和图 7(c)、图 7(d)分别表示在27°N和36°N纬度面用Scaling-ART算法和MART算法得到的残差。可以看出，在同一纬度下，使用两种算法得到的残差最大值为1.4×10¹⁰ el/m³。使用Scaling-ART算法得到的残差红色区域面积更小，表示残差更小，反演效果优于MART算法。从图 7(a)与图 7(c)对比来看，27°N处的反演残差明显大于36°N处，与图 6结论一致，在测量信号分布较少的空间内，反演精度相对较低。

图 8表示反演区域内各经度分别用上述两种算法反演得到的电子密度与真值之间的均方根误差。误差定义如下：

$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{j = 1}^n {({x_{r, j}} - {x_{e, j}})} }}{n}} $

(10)

式中，x_r表示电子密度真值，在本实验中为IRI模型值; x_e表示使用重构方法得到的电子密度重构值(或称为估计值)。

可以看出，在反演区域的各个经度上，本文Scaling-ART算法反演得到的电子密度均方根误差相比于MART算法更小，说明改进算法的有效性和优越性。图中所示的均方根误差对于两种反演算法而言，都有随经度增加而逐渐减小的趋势，主要是实验所选取的观测射线在不同经度上分布不均匀导致。具体而言，在95°~105°E范围内，二者反演结果相差较大，而在110°E附近，反演误差差别较小。从图 3接收机和穿刺点分布来看，95°~105°E范围内穿刺点数量很少，本文算法能够通过尺度因子在整个反演区域的作用使各个网格的电子密度逐渐接近于真值；而在穿刺点分布较为均匀且密集的区域，很多网格都被信号穿过，使用原有MART算法也能达到较好的效果。在GPS获取STEC精度较高、测量射线分布稀疏的情况下，可以考虑利用本算法进行电离层电子密度的重构。

3.2 实测数据验证

实验选区与上述模拟实验相同，在118°E、34°N附近共有4台GPS接收机分别位于汉中、平凉、西安、襄阳。实验起始时间为2014-08-25 18：00，选取的时段长度为1 h，采样间隔为30 s，截止高度角为30°。根据GPS卫星坐标和地面布设的接收机经纬度计算出真实截距矩阵，STEC亦通过GPS数据处理得到。实验中在所有符合条件的测量射线中随机选取一部分和其对应的STEC进行电离层电子密度的反演，剩余测量射线和STEC用于验证反演效果。对用于测试的测量射线，计算两种算法重构之后的STEC，求出其与实测STEC的差的绝对值，err=|STEC_r-STEC_e|，之后将两种算法得到的STEC绝对误差进行比较：

$ {\rm{dif}} = {\rm{er}}{{\rm{r}}_{{\rm{MART}}}} - {\rm{er}}{{\rm{r}}_{{\rm{Scaling - ART}}}} $

(11)

图 9(a)表示两种算法得到的重构STEC与实测STEC的比较，红色星标实线表示使用MART算法的绝对误差，橙色虚线表示真值，蓝色表示本文算法重构得到的STEC。本文算法得到的STEC误差大约在0~3TECu，说明该算法能够达到较好的重构效果。图 9(b)表示式(11)所定义的dif值，dif＞0表示使用MART方法的绝对误差更大，反之使用本文算法误差更大，落在红线上的点表示两种算法得到的反演精度相当。从图 9可以看出，在用于验证的10组观测射线中，使用本文算法进行反演之后，多数情况下绝对误差小于原有的MART算法，说明本文改进算法相比原始算法具有优越性。

4 结语

本文提出在代数重构中对每个网格加入不同尺度因子的算法。数值模拟实验中，首先计算STEC“实测”值，每次迭代中结合当前电子密度重构值得到STEC重构值，将二者的比例作为测量射线穿过网格的缩放因子；再根据其他网格中心到观测射线的距离计算每个网格对应的尺度因子，通过不断的迭代修正电子密度重构值。本算法解决了原有算法中没有测量射线穿过的网格电子密度重构结果严重依赖于初始化值的问题。在实测数据验证中，使用部分实测STEC作为已知观测值，其余实测STEC用于验证，分别利用两种方法进行电子密度重构。通过比较两种方法在测试集上的STEC绝对误差，证明了本文方法的有效性以及优越性。未来将针对电离层的各向异性进一步考虑尺度因子的选取方案，以达到更好的重构效果。