GPS电离层层析成像技术以其能够反演电子密度三维结构日益受到广泛重视.多年来，国内外学者已在GPS电离层层析反演方法与技术方面开展了诸多研究，取得了许多重要成果(Austen et al.，1988; Bust and Mitchell，2008; Raymund et al.，1990; Pryse et al.，1993; Xu and Zou，2003，Jin et al.，2006; Wen et al.，2007a; Yao et al.，2014a;Yuan et al.，2007; Wen et al.，2007b; 徐继生等，2005; 施闯等，2010).其中，以代数重构为代表的行迭代类重构算法是一类重要的层析反演方法，包括：加法代数重构(ART)(Austen et al.，1988)、乘法代数重构(MART)(Raymund et al.，1990)、同时迭代重构(SIRT)(Pryse et al.，1993)等.该类算法通过将GPS射线在观测方程组成的超平面内进行投影迭代，逐步缩小观测值与投影重构值之间的差距，进而估算出最终的电子密度结果.在层析投影迭代重构反演过程中，层析算法通常都假定待反演区域的电离层电子密度按照一定的空间格网间隔离散分布，每个像素格网内的电子密度均匀分布，在此基础上，以GPS TEC观测射线和层析系统像素格网交叉形成的截距长度为权重比，对TEC层析投影重建值和GPS TEC实测值之间的误差进行分配(Austen et al.，1988; Bust and Mitchell，2008; Raymund et al.，1990; Pryse et al.，1993).然而，需要指出的是，投影重建的TEC贡献来源于GPS射线截距与格网像素内电子密度估算值之间的乘积，GPS射线截距在迭代反演过程中保持不变，电子密度误差是TEC实测值与其投影重构估算值差异的主要来源，也就是说，在GPS层析系统迭代反演误差中，电子密度误差是决定因素，GPS射线截距对误差起放大作用.因此，传统ART、SIRT、MART等GPS层析迭代反演算法中仅以射线截距作为误差分配准则的唯一要素是不合理的.

为提高层析反演电离层电子密度的精度和计算效率，国内外一些学者考虑在GPS层析迭代模型中引入电子密度参量.Mitchell等(1997)和Pryse等(1998)给出了一种加权迭代重构算法，在经典ART算法的基础上，引入了当前电离层电子密度相对于GPS观测视线方向的最大电子密度的比值作为迭代重构时的权重，使得电离层电子密度的修正值成比例于电离层初值模型提供的电子密度；Wen等(2007b)提出一种修正的迭代重构算法，在GPS每条观测射线的每轮迭代修正中引入了电子密度参数作为迭代参量加速收敛.实际上，上述改进方法的数学原理类似，本质上都是确保反演的电离层电子密度垂直剖面结构相对于背景电离层模型保持稳定，弥补水平方向观测射线的不足，提高GPS电离层层析系统反演计算效率与精度(Das and Shukla，2011)，但在电离层扰动异常情况下有可能畸变电子密度反演结构和演变形态.

此外，现有GPS电离层层析迭代算法中的松弛因子通常都取经验的固定常数实现控制噪声对层析反演结果的影响.然而，已有研究表明：这种迭代反演算法解决不适定问题并不总是得到好的结果，主要原因是迭代算法会产生极限环，迭代过程在极限点间摆动，进而导致反演解算结果的不稳定，甚至出现反演质量恶化(Ahn et al.，2006; Byrne and Graham-Eagle，1992; Kak et al.，2002).换句话说，该类方法将松弛因子取值为固定经验常数，忽视了松弛因子对电离层电子密度反演精度和反演结果的平滑程度加以调节并使两者之间取得平衡的作用.传统取固定经验常数作为松弛因子的策略方法，将导致噪声误差在电子密度像素格网系统内传播的不平衡进而影响GPS电离层层析反演的电子密度精度.需要注意的是，在利用卫星信号开展对地观测时受到了仪器噪声和传播噪声两类误差的影响(Armstrong，1998; Asmar et al.，2005; Rawer，1962).仪器噪声包括了卫星和接收机硬件设备引起的白噪声及与设备信噪比变化密切相关的相位信号波动；传播噪声则是卫星信号穿越大气介质时由于大气密度变化导致的折射率波动进而引起的随机相位信号波动(Armstrong，1998; Asmar et al.，2005).因此，对于GPS电离层层析观测系统而言，GPS卫星射线穿越不同方位与不同高度的电离层层析像素格网时受不同像素格网内电子密度变化引起的传播噪声影响是不同的.针对迭代重构层析算法中松弛因子的不同选择对反演结果的影响，Qu等(2007)、Jiang和Wang(2003)利用矩阵的奇异值定理研究了松弛因子的选取对迭代重构算法收敛性的影响；Elfving等(2012，2010)则基于分析迭代重构算法的“半收敛”现象，提出选取松弛因子的方法与策略；Yao等(2014b)提出利用Elfving等(2012)提供的松弛因子选择策略，进行电离层电子密度联合迭代重构自适应反演研究.但上述研究是从纯数学角度研究与讨论松弛因子的选取问题，没有从电离层实际物理变化状态中考虑松弛因子的选取策略.

针对上述问题，本文一方面，提出以GPS观测射线穿过层析像素格网系统形成的射线截距与对应的像素格网内电子密度乘积为组合自变量(代表该像素格网内电子密度对观测射线方向的TEC贡献)，寻求合适的GPS电离层层析迭代表达式，重新在不同像素格网内更合理分配TEC实测值与反演值之间的差距；另一方面，考虑到GPS信号穿越电离层介质时的传播噪声与电子密度变化相关，提出基于电离层电子密度变化重构松弛因子的研究思路，削弱噪声对反演电子密度的影响，最终达到提高电离层电子密度反演精度.

代数重构算法(ART)是一种迭代的级数展开法，由20世纪初的数学家Kaczmarz(1937)在求解稀疏相容的线性方程组时提出，后被Gordon等(1970)将其引入图像重建领域.Austen等(1988)首次利用该算法开展电离层电子密度层析三维反演研究.该算法通常在迭代之前利用经验电离层模型给电离层区域内的每个像素赋予一个初值，然后采用迭代的方式逐步改善待重构电子密度估值.每一步迭代对应于一条斜距GPS TEC测量(如第i条)，即针对一个方程进行，完成一次全部的斜距GPS TEC测量迭代称之为一轮迭代，修正的依据是利用第k次迭代计算的电离层电子密度求出的斜距TEC与实际测量的斜距GPS TEC之差，对电离层电子密度分布做相应的修正，使迭代结果逐渐收敛.在第k次迭代中，该方法以如下方式修正电离层电子密度：

(1)

其中

式中，x表示GPS电离层层析系统像素格网代表的电子密度；y表示GPS观测射线对应的电离层TEC；i表示一条GPS观测射线，共有m条射线数(i=1，…，m)；k是层析反演迭代次数；λ为层析迭代反演松弛因子，通常取为固定经验常数且0＜λ＜2；M 表示层析模型的迭代表达式；j表示GPS层析系统像素格网，总数量为n(j=1，…，d，…，n)，d表示其中的第d个像素格网；a_jⁱ(a_dⁱ)是GPS第i条观测射线穿越第j(d)个电子密度像素格网时所形成的截距，射线没有穿越的像素格网截距大小则为0.从几何观点来看，每一步迭代相当于将所有电子密度 x 集合组成的向量向第i个方程所代表的超平面进行投影，投影程度由λ控制，当λ=1时为完全投影.ART算法不保证反演结果为正，因此利用其进行电离层层析重建时，在迭代过程中必须加上电离层电子密度各元素非负的约束条件.图 1给出了一条GPS观测射线穿越电离层层析系统像素格网时的示意图.

图 1

Fig. 1

图 1 GPS射线穿越电离层层析像素格网示意图 Fig. 1 Illustration of the GPS ray passing through an example reconstruction grid of ionospheric tomography

根据式(1)可以看出，电离层电子密度误差是ART层析算法重构TEC值(反演的电子密度积分值)与GPS TEC实测值之间存在差异的主因，GPS射线穿过层析像素格网系统形成的射线截距主要起放大作用.考虑上述问题，本文提出引入GPS射线截距与电子密度的乘积为组合自变量(该组合变量是对应的层析像素格网中的电子密度在整条GPS观测射线方向TEC的组成部分，物理意义明确)，构造新的 GPS电离层层析迭代模型，重新在不同电子密度像素格网内合理分配GPS TEC实测值与TEC反演值之间的误差.参照式(1)，提出新的电离层层析迭代表达式如下：

(2)

如前所述，松弛因子的主要作用是平衡与调节电离层层析反演的电子密度精度与结果的平滑程度.松弛因子取值较大时，反演的电子密度结果将较为平滑，掩盖了电离层局部变化特征；松弛因子取值较小时，又可能受噪声干扰导致电子密度反演精度降低.除仪器设备引起的白噪声之外，GPS电离层层析反演噪声还包括GPS信号传播时由于介质电子密度变化引起的传播噪声.据此，本文设计了一组与电子密度变化相关的松弛因子λ的表达式：

(3)

综上所述，新的GPS电离层层析反演模型如下：

(4)

式中，λ_d^i，k与M_dⁱ的表达式分别见式(2)与式(3).

本文实验研究采用了2011年11月27日中国陆态网络131个地面基准站的高精度GPS观测数据，地面观测站分布大体均匀，观测数据以国际GNSS服务组织(IGS)标准采样(30S)获取，卫星截止高度角为15°.分别利用各个观测站的双频载波相位平滑码数据，并采用“两步法”消除站星硬件延迟影响(Li et al.，2012)，计算出对应各个观测站与卫星观测射线方向上的斜电离层TEC值.

电离层电子密度层析反演区域所选定的经度范围为75°E—135°E，纬度范围为10°N—55°N，高度范围为90~990 km；同时在对选定的电离层反演区域进行三维空间离散化时在经度和纬度方向上分别取为5°与2.5°的空间间隔，在电离层垂直高度方向上取30 km的空间间隔.

此外，IRI2007电离层模型提供的电子密度作为层析反演时的电子密度初值；为验证本文提出的电离层层析反演新算法的可靠性，分别利用北京(40.3°N，116.2°E)和三亚(18.3°N，109.6°E)的电离层测高仪获得的两站上空电离层电子密度剖面信息与新算法反演的电离层电子密度结果进行分析对比.

这里仅以部分结果为代表展示，其他时段结果与此类似.图 2分别给出了北京站和三亚站上空在2011年11月27日北京时间10 : 00、12 : 00、14 : 00与16 : 00的电子密度沿高度方向上的剖面变化结果.同时，采用电离层测高仪观测获取的电子密度剖面信息作为参考“真值”和不同的电离层层析反演算法结果进行比较.

图 2

Fig. 2

图 2 利用地基GPS观测数据分别采用代数重构算法(ART)和本文提出的新算法反演北京(40.3°N，116.2°E)上空电离层电子密度与电离层测高仪观测所得的电子密度剖面之间的比较结果 Fig. 2 Comparison of the ionospheric electron density profiles over Beijing station(40.3°N，116.2°E)derived from the ART and new proposed algorithms using the ground-based GPS observations with the profiles obtained by the ionosondes

根据图 2可以看出，在北京站上空，本文提出的GPS电离层层析反演新方法给出的峰值电子密度及其底部电离层电子密度都与电离层测高仪观测获得的结果符合一致，尤其在12 : 00、 14 : 00与16 : 00，新方法反演的峰值电子密度结果与测高仪提供的峰值电子密度结果基本相吻合，相对于传统ART算法反演的峰值电子密度结果，精度有了较好的提高.

然而，无论是新方法还是传统ART算法，反演的电离层电子密度峰值高度以上的顶部电子密度结果与电离层测高仪获得的电子密度变化趋势一致，但数值结果存在明显差异.需要说明的是，电离层测高仪能够通过直接观测获取峰值电子密度及其高度以下的底部电子密度准确结果，电子密度峰值高度以上的顶部电子密度却是根据一定的算法获得(Liu et al.，2014; Hu et al.，2014)，因此，本文比较不同方法获得的顶部电子密度时，仅进行定性的趋势变化比较而不开展定量的结果差异比较.

从图 3可以看出，在低纬度三亚观测站上空，由于电离层变化活动复杂，传统ART层析反演算法给出的电离层电子密度相对于电离层测高仪观测得到的电子密度“真值”差异较大，无论是在底部电子密度剖面，还是在顶部电子密度剖面，都和测高仪观测获得的电离层电子密度演变形态不符合；然而，新的电离层层析反演方法给出的电子密度结果，从总体上看与电离层测高仪观测获得的电子密度变化趋势保持一致，此外，在12 : 00与16 : 00，无论是相对于电离层测高仪获得的峰值电子密度还是底部电子密度结果，新方法反演的电子密度均与其相接近.这也从一个侧面说明了本文提出的新的电离层层析反演方法的合理性及其相对于传统ART算法的优越性.

图 3

Fig. 3

图 3 利用地基GPS观测数据分别采用代数重构算法(ART)和本文提出的新算法反演三亚(18.3°N，109.6°E)上空电离层电子密度与电离层测高仪观测所得的电子密度剖面之间的比较结果 Fig. 3 Comparison of the ionospheric electron density profiles over Sanya station(18.3°N，109.6°E)derived from the ART and new proposed algorithms using the ground-based GPS observations with the profiles obtained by the ionosondes

但必须说明的是，在部分时段(如10 : 00与14 : 00)，新方法反演的三亚测站上空的峰值电子密度结果相对于电离层测高仪观测获得的峰值电子密度仍存在一定的差异，这也说明对于在电离层变化活动较为复杂的低纬区域，新方法存在一定的局限性，这也是本文工作后续仍需研究与改进的地方.另外，在本文实验研究中电离层层析反演的低纬地区，GPS观测数据相对稀疏且处于层析反演的边缘区域，这可能也是导致该区域的电离层电子密度反演精度降低的原因之一.

考虑本文研究中采用的GPS数据在低纬区域分布稀疏并且处于反演区域边缘，本小节给出的电离层电子密度随高度与纬度剖面的演变图中，纬度范围限定在20°N—45°N，高度范围限定在90~990 km，并且以110°E上空电离层电子密度的剖面反演结果为例进行讨论.

图 4给出了分别采用ART算法(a—d)与新算法(e—h)反演得到的电离层电子密度随高度与纬度在不同时间内的剖面演变图.可以看出，从北京时间10 : 00到16 : 00，新算法反演的电离层电子密度随纬度变化的峰结构向中纬度方向的扩展程度比ART算法反演的电子密度峰结构变化较大，这与图 2中反映出的北京上空利用新算法反演的电子密度较ART反演的电子密度大并且更接近于电离层测高仪观测得到的电子密度一致(尤其在12 : 00、 14 : 00与16 : 00的结果较为明显)；此外，图 4中也可以看出，在北京时间12 : 00与14 : 00，ART算法反演的低纬区域在900 km左右和100 km左右的电子密度值大于新算法反演的在该区域电子密度值(即图 4b中12 : 00时刻和图 4c中14 : 00时刻，在上述两处代表电子密度量级的颜色相对偏蓝色而不是紫色)，表明在上述区域中ART反演的电子密度值大于新算法反演的电子密度值，这与图 3中反映的三亚上空电子密度结果变化一致，即ART算法在该处反演的电子密度存在较大异常值，和实际的电离层形态变化不符合.综上所述，在图 4中，新算法更好地反映了电离层电子密度随高度与纬度在不同时间的演变形态.

图 4

Fig. 4

图 4 ART算法(a—d)和新算法(e—h)反演的电离层电子密度随高度与纬度变化剖面图比较(反演时间段从左至右分别是2011年11月27日10 : 00、12 : 00、14 : 00与16 : 00) Fig. 4 Ionospheric imaging of the electron density distribution reconstructed from the ART algorithm(a—d)and the new algorithm(e—h)in a height versus latitude plane(Inversion results from left panel to right panel are respectively presented at 10 : 00，12 : 00，14 : 00 and 16 : 00 Beijing time on November 27，2011)

进一步以电离层测高仪观测提供的电子密度为“真值”，定量统计不同时刻GPS电离层层析方法反演的峰值电子密度误差的绝对值、底部电离层电子密度平均误差的绝对值以及在对应时间段内的RMS值，表 1与表 2分别给出电离层电子密度在北京与三亚上空于地方时10 : 00、12 : 00、14 : 00 与16 : 00时刻相关的重构反演统计结果.

表 1(Table 1)

表 1 ART与新方法反演重构的北京站(40.3°N，116.2°E)上空不同时刻的电子密度相对于电离层测高仪观测“真值” 的误差统计分析，电子密度单位为1012el·m-3 Table 1 Error statistics analysis of the ionospheric electrondensity derived from the ART and new algorithm in comparisonwith those from ionosonde measurements over Beijing Station(40.3°N，116.2°E). The unit of ionospheric electron density is 1012el·m-3 北京时间 峰值电子密度误差绝对值 底部电子密度平均误差绝对值 ART 新方法 ART 新方法 10:00 0.128 0.083 0.079 0.029 12:00 0.040 0.004 0.045 0.031 14:00 0.180 0.009 0.095 0.028 16:00 0.126 0.024 0.080 0.052 RMS 0.129 0.044 0.085 0.048

表 1 ART与新方法反演重构的北京站(40.3°N，116.2°E)上空不同时刻的电子密度相对于电离层测高仪观测“真值” 的误差统计分析，电子密度单位为1012el·m-3 Table 1 Error statistics analysis of the ionospheric electrondensity derived from the ART and new algorithm in comparisonwith those from ionosonde measurements over Beijing Station(40.3°N，116.2°E). The unit of ionospheric electron density is 1012el·m-3

从表 1可以看出，在北京站上空，相对于电离层测高仪观测获得的电子密度“真值”，无论是峰值电子密度误差的绝对值、底部电子密度平均误差的绝对值还是在对应时间段内的RMS值，新方法的反演误差结果都小于ART算法的反演误差结果；从表 2结果可以看出，在三亚站上空，除在14 : 00时与16 : 00时新方法反演的峰值电子密度误差略大于ART算法反演的峰值电子密度误差，其他结果都显示新方法的反演误差及RMS值小于ART算法的反演误差及RMS值，这与前述的图 3与图 4结果相符合，从而进一步说明了新算法相对于传统ART算法的优越性.此外，也可以从结果上看出，由于低纬地区电离层复杂多变且观测数据稀疏，新方法在中纬度地区反演电离层电子密度的精度优于在低纬度地区的反演精度.

表 2(Table 2)

表 2 ART与新方法反演重构的三亚站(18.3°N，109.6°E)上空不同时刻的电子密度相对于电离层测高仪观测“真值” 的误差统计分析，电子密度单位为1012el·m-3 Table 2 Error statistics analysis of the ionospheric electrondensity derived from the ART and new algorithm incomparison with those from ionosonde measurements over Sanya Station(18.3°N，109.6°E). The unit of ionospheric electron density is 1012el·m-3 北京时间 峰值电子密度误差绝对值 底部电子密度平均误差绝对值 ART 新方法 ART 新方法 10:00 0.640 0.571 0.544 0.208 12:00 0.181 0.020 0.217 0.114 14:00 0.048 0.294 0.313 0.145 16:00 0.189 0.220 0.448 0.219 RMS 0.347 0.340 0.599 0.272

表 2 ART与新方法反演重构的三亚站(18.3°N，109.6°E)上空不同时刻的电子密度相对于电离层测高仪观测“真值” 的误差统计分析，电子密度单位为1012el·m-3 Table 2 Error statistics analysis of the ionospheric electrondensity derived from the ART and new algorithm incomparison with those from ionosonde measurements over Sanya Station(18.3°N，109.6°E). The unit of ionospheric electron density is 1012el·m-3

需要补充说明的是，本文提出的电离层层析反演改进新算法是利用部分中国陆态网络数据进行的探索性研究，尽管在一定程度上提高了电离层电子密度反演结果精度，但根据前述讨论可以看出，在电离层变化活动复杂地区(如低纬度区域)，其改进精度仍然有限，更优的GNSS电离层层析反演算法有待利用更长时间的观测数据进行深入探讨和研究.

本文提出一种新的电离层电子密度三维层析反演方法，相对于传统的ART层析反演算法，新方法能够更合理地在GPS电离层层析像素格网内调整与分配TEC实测值与层析反演值之间的差异，而且通过设计了一组新的松弛因子控制与削弱了噪声对电子密度反演结果的影响，从而有效提高电离层电子密度三维反演的精度.此外，新方法在低纬度地区反演的电离层电子密度精度相对于在中纬度地区降低，这一方面与低纬度地区电离层复杂多变有关，另一方面也与处于层析反演边缘区域的低纬地区的GPS观测数据较少相关，这是后续研究有待进一步深入的地方.

感谢中国大陆构造环境监测网络提供GPS观测数据，感谢中国科学院地质与地球物理研究所刘立波研究员提供的电离层测高仪资料.