虚拟射线辅助的三维电离层层析方法

范玉荣， 符杰林， 安 涛， 王俊义， 林基明

(1.桂林电子科技大学 认知无线电与信息处理教育部重点实验室，广西 桂林 541004； 2.中国科学院 上海天文台，上海 200030; 3.桂林电子科技大学 广西高校卫星导航与位置感知重点实验室,广西 桂林 541004)

电离层内大量存在的自由电子会影响无线电信号在电离层中的传播，因而电离层电子密度分布信息不仅对电离层物理具有重要意义，而且对星地通信和卫星导航也具有重要意义[1]。传统的二维电离层建模假定电离层区域所有的自由电子都集中在一定高度的球壳薄层上，虽然也能达到较好的精度，但无法反映电离层的空间分布情况[2]。

GNSS电离层层析成像根据反演区域内大量信号传播路径上的斜向总电子含量(slant total electronic density, 简称STEC)来反演特定区域的电子密度分布，能够有效反演得到三维电离层信息。相比传统电离层探测方法，GNSS电离层层析具有探测持续时间长、探测范围广及成本低等优势，近年来在电离层探测中被广泛应用[2-3]。然而，由于GNSS卫星观测视角和地面观测站数量及分布的限制，基于GNSS的电离层层析通常存在因数据不足而引起的不适定问题[2-10]。为解决不适定问题带来的影响，国内外学者提出了各种解决方法。文献[4]采用同时迭代重建技术(SIRT)进行电离层电子密度的反演，减少了噪声对解算结果的影响；闻德保等[5]通过代数重构算法(ART)对松弛因子进行了改进，提出了IART算法；姚宜斌等[6]针对联合迭代重构算法迭代收敛慢且易受噪声影响的问题，利用上一轮迭代的电子密度反演结果，自适应地调整松弛因子和加权参数，提出了ASIRT算法；文献[7]提出了一种TV-MART算法，用总变差(TV)最小化结合MART算法对电子密度进行反演，减少了由噪声引起的不稳定性；文献[8]用自搜索松弛因子对IART算法进行改进，减少了由于松弛因子不合适而引起的误差。以上方法在一定程度上解决了反演过程中的不适定问题，提高了反演精度，但仍存在由于观测站分布不均和观测角度有限带来的很多格网依旧缺乏射线穿过而得不到有效修正的问题，尤其是靠近边缘部分的格网，电离层层析成像的病态性依然突出[9]。

文献[10]提出一种约束自适应联合迭代反演方法，通过对观测矩阵附加平滑约束，对无射线穿过的格网电子密度进行了一定的修正，但当无射线穿过的格网较多时，其精度提升并不明显。文献[11]将边缘射线用于电离层层析中，提高了射线格网覆盖率，显著提高了反演精度，但该反演精度取决于处于层析范围内的边缘射线的总电子含量观测值的精度，而这个精度较难控制。

鉴于此，通过在地面缺乏观测站的区域设立虚拟观测站，建立一系列虚拟观测射线，再建立二维薄壳电离层模型，并根据投影函数计算出这些虚拟观测射线的STEC值。这些虚拟观测射线能够有效扩大射线的覆盖范围，使得原本无射线穿过的格网的电子密度也能得到有效修正，从而得到更加准确的电离层电子密度分布。

1 系统模型及问题分析

电离层层析成像使用GPS射线的电离层STEC反演电离层电子密度分布。STEC可从双频GPS接收机的伪距和载波相位观测值中提取[12]：

(1)

其中：P4,sm为载波相位平滑伪距之后的观测值；c为光速；DCBi、DCBj分别为卫星和接收机的差分码偏差(differential code biases, 简称DCB)。

获得的STEC观测值可表示为电子密度沿信号路径上的积分，

(2)

其中：Ne为电子密度；l为信号路径；r为t时刻经度、纬度和高度所组成的位置向量。将反演区域按照经度、纬度、高度方向上划分为三维格网，则每条路径上的STEC观测值可表示为

(3)

其中：m为穿过电离层的射线总数；n为格网总数；aij为第i条射线在第j个格网内的截距；xj为第j个格网的电子密度；ei为观测噪声与其他观测误差的和。电离层层析模型的简化示意图如图1所示。

图1 电离层层析模型简化示意图

式(3)的矩阵形式可表示为

y=Ax+e，

(4)

其中：y为GNSS信号射线传播路径上电离层STEC构成的m维列向量；A为GNSS信号射线穿过格网时的截距构成的m×n维投影矩阵；x为所有格网像素中心电子密度构成的n维列向量；e为噪声和观测误差构成的列向量。

由于观测卫星数有限，式(4)中的投影矩阵A通常是一个稀疏矩阵，运用迭代重构算法对式(4)进行解算。由于MART迭代速度较快，且能解决反演值为负的问题，在电离层层析中被广泛使用，其反演迭代公式为

(5)

在实际观测中，受地面观测站分布不均和卫星观测角度的限制，部分格网无任何射线穿过，因而这些格网的电子密度依赖于初始经验值，得不到有效修正。通过增加虚拟射线的方式增加射线穿过格网的覆盖率，通过事先建立的薄壳电离层模型结合投影函数获得虚拟射线的STEC观测值。

2 利用虚拟射线改进的反演方法

薄壳电离层模型由于其结构简单且拥有不错的精度而在电离层建模方法中被广泛应用[13]。本研究采用低阶球谐函数模型进行电离层建模。电离层球谐函数模型[14]为

(6)

根据式(6)用足量的STEC观测值建立方程组，利用最小二乘法可解算出球谐函数未知球谐系数。通过二维薄壳模型和投影函数可得到区域内任意位置任意仰角的STEC值，且观测仰角越高，其精度越高[15]。

在地面观测站分布较少的区域设置虚拟观测站，将其周围仰角较高的一系列射线作为虚拟射线，利用二维薄壳模型和投影函数计算虚拟射线的STEC值，与实际观测射线一起依据式(4)、(5)进行反演重构。设置虚拟射线的具体操作流程如下：

1)将反演区域底部按经纬度划分为2°×1.5°的格网，将上下左右的格网内均无观测站的格网的中心位置设为虚拟观测站位置。

2)将反演区域顶部按经纬度划分为0.5°×0.5°的格网，每个格网中心点与地面虚拟观测站连线，从中选择仰角不低于60°的射线作为虚拟射线。为使虚拟射线通过尽可能多的格网，将虚拟观测站位置高度设置为-50 km。

3)计算虚拟射线在450 km高度面的穿刺点位置及天顶角，利用之前得到的薄壳电离层模型和投影函数计算得到虚拟射线的STEC值。

4)将虚拟射线与实测射线按式(5)构建矩阵，用MART进行解算。

3 数据处理及分析

采用欧洲地区2018年9月20日的IGS观测站的GPS观测数据进行实验。反演区域设为纬度35°N～55°N，经度5°W～25°E，高度范围为100～1 000 km。在纬度、经度和高度上的空间分辨率分别取1°、2°、20 km。IGS观测站及虚拟观测站的分布如图2所示。其中，深色·表示实际观测站，浅色· 表示虚拟观测站，深色×表示测高仪观测站。

图2 观测站分布图

3.1 DCB及薄壳电离层模型精度评估

DCB是提取STEC观测值中的主要误差来源，精确估计DCB可提高STEC精度，同时也可得到更加精确的薄壳电离层模型。单独用24 h的GPS观测数据解算得到当天的卫星和接收机DCB信息，在之后的薄壳电离层模型解算和电离层层析中都利用这些DCB信息以得到更加精确的STEC观测值。为验证解算得到的DCB的精度，将解算得到的DCB与中科院(CAS)发布的DCB数据进行对比，得到如图3所示的差值折线图。

图3 解算得到的DCB与CAS的DCB差值

由于CAS给出的接收机DCB数据有限，图3(b)只包含CAS给出了数据的部分观测站的DCB差值。从图3可看出，解算得到的GPS卫星DCB和接收机DCB与CAS发布的DCB差值均在1 ns内，其中卫星DCB的最大误差为0.537 ns，观测站接收机DCB的最大误差为0.945 ns。这表明，解算得到的DCB具有较高的精度，从而保证了通过提取得到的STEC的可靠性。

用30 min内GPS观测数据进行电离层薄壳模型解算，将解算得到的薄壳模型与CODE发布的全球电离层图GIM进行对比。解算得到的每个时刻的电离层薄壳模型与GIM的差值均值和均方根如表1所示。从表1可看出，解算得到的薄壳模型与GIM的误差基本在1.3 TECU内。这表明，得到的薄壳模型精度较高，同时也保证了添加的虚拟射线的STEC精度。

表1 薄壳模型与GIM的差值均值与均方根

3.2 模拟数据实验

由于对反演区域整体真实的电离层状态未知，直接用实测数据进行电离层电子密度反演无法全面评判算法的有效性和稳定性，因而用数值模拟实验验证算法的有效性和稳定性。根据观测站位置和卫星位置计算射线在每个格网中的截距并构建投影矩阵A。其中，观测站坐标和卫星坐标均为2019年9月20日UT11:45～12:15观测间隔内的真实值。利用IRI2016模型获得2019年9月20日UT12:00时刻待反演区域内各格网中心点处的电离层电子密度xIRI。将A与xIRI相乘，得到各条射线上的模拟STEC。为了更接近真实情况，在模拟STEC数据上增加了最大为模拟STEC均值10%的随机噪声。同时，将初始电离层电子密度设置为0.8xIRI。用xIRI作为真实值进行评估。

图4为添加了虚拟射线后的射线分布图，其中浅色线为实际观测射线，深色线为添加的虚拟观测射线。从图4可看出，虚拟射线填补了由于实际观测站和卫星位置分布限制导致的格网覆盖空白，使得射线格网覆盖率大大提升。经统计，添加虚拟射线后格网覆盖率达到了90%以上，比添加前的格网覆盖率提高了约30%。

图4 射线分布图

图5为2种方法迭代30次的均方根误差曲线。其中：3DCIT表示传统电离层层析方法；3DCIT-VR表示添加了虚拟射线的电离层层析方法。从图5可看出，添加虚拟射线后，由于格网覆盖率的优势，其一步迭代后的误差更小，且误差下降速度更快。最终收敛时，3DCIT-VR、3DCIT的均方根误差分别为0.063 9×1011、0.120 6×1011el/m3。

图5 2种方法迭代误差

图6为2种方法反演得到的电离层电子密度误差经度切片。从图6(a)可看出，未添加虚拟射线时，由于大量格网无射线穿过而未得到任何修正，因而其误差较大，尤其在地面测站分布较少的区域呈现出明显的不连续性。从图6(b)可看出，3DCIT-VR的误差明显更小，且更加均匀，这表明添加虚拟射线能够提高电离层层析精度。

图6 2种方法反演得到的电离层电子密度误差经度切片

3.3 实测数据实验

为进一步验证本方法的优越性，用30 min内实际观测数据进行电离层层析反演。先用30 min内观测数据建立区域薄壳电离层模型，再算得虚拟射线的STEC值，最后将虚拟观测数据与实测数据一同进行迭代反演，得到电离层电子密度分布。将反演结果与未添加虚拟射线时的反演结果进行比较，并使用PQ052和AT138两个电离层测高仪的实测数据对层析结果进行验证。

图7为用2种方法在UT6:00、UT12:00和UT18:00三个时刻反演得到的电子密度与相应时刻PQ052、AT138两个电离层观测站的观测数据对比曲线。由图2可知，PQ052位于观测站较为丰富的区域，AT138周围则缺乏实际观测站。从图7(a)可看出，2种方法得到的反演结果较接近。这表明添加虚拟射线并未降低实际射线丰富区域的反演精度，且有时对反演精度还有一些改善。从图7(b)可看出，3DCIT-VR反演得到的电离层电子密度剖面更接近电离层实际观测值。这也表明了添加虚拟射线可在一定程度上提高电离层反演精度，从而进一步证明了本方法的优越性。从图7还可看出，反演得到的电离层电子密度剖面峰值高度与实测值峰值高度存在一些差异，这是因为观测射线仰角偏高导致垂直分辨率不高，未来可进一步引入掩星观测数据进行反演以改善这个问题。

4 结束语

通过在缺乏地面观测站的位置设置虚拟观测站，以添加虚拟射线的方式增加可用射线数量，进而提高格网覆盖率。本方法在保证射线丰富区域反演精度的同时，使无射线穿过的格网电离层电子密度也能得到有效修正。通过模拟数据实验及实测数据实验验证了本方法的有效性，并利用2个测高仪观测站(PQ052、AT138)分别对射线丰富区域和射线欠丰富区域的电离层反演效果做了验证。实验结果表明，本方法效果显著，解决了缺乏射线穿过的格网电离层电子密度得不到有效反演的问题。

图7 2种方法反演得到的电子密度与PQ052和AT138两个电离层观测站观测数据对比