CN113093224B - 一种边缘增强的电离层层析方法 - Google Patents

Abstract

本发明提供一种边缘增强的电离层层析方法，包括如下方法步骤：步骤1)，利用GPS码观测值与载波相位观测值联合解算测量的电离层TEC；步骤2)，建立虚拟参考站，通过三角形线性插值计算虚拟TEC值；步骤3)，将待反演区域按照经度、纬度、高度方向均匀划分成一个个网格，求出全部观测信号在每个网格中的截距后，按照网格和观测次数的顺序排列组成观测矩阵A；步骤4)，采用IDW插值方法作为水平约束，使没有信号射线穿过的网格参与层析迭代；步骤5)，构造边缘函数，通过ART迭代算法，求解电离层电子密度分布。本发明充分考虑实际反演中边缘地区观测信号稀缺的情况，对ART算法的修正式加权以限制对边缘地区的修正，防止过度拟合。

Description

一种边缘增强的电离层层析方法

技术领域

本发明涉及航空航天技术领域，特别涉及一种边缘增强的电离层层析方法。

背景技术

电离层是指地面以上约50千米至1000千米的空间大气，该范围内的地球大气分子由于受到太阳紫外线辐射与地球上层大气原子、分子、宇宙射线中高速粒子相互作用而发生了电离，生成大量自由电子和离子。

电离层是日地系统一个的重要组成部分，也是人类生存的近地空间环境中的重要组成部分，对人类生产与活动具有十分重要的影响。一方面，正常情况下电离层的反射有助于远距离无线电通讯的实现和发展；另一方面，电离层的剧烈变化(如电离层扰动、电离层中的不规则体)会对短波通讯、航空航天和卫星导航等领域产生严重的影响，可能致使航天器受损、卫星和地面通讯设备的失灵以及高压输电网损害等。常用的电离层观测手段如数字斜测仪、电离层垂测仪、大功率散射雷达等，一般只能获得观测站所在位置上电离层相关物理量随高度和时间变化的信息，不能对电离层的顶部进行探测，并且造价昂贵、维护代价较大。

随着卫星技术的发展，将卫星探测和计算机层析技术(ComputerizedTomography，CT)相结合的电离层层析技术(Computerized Ionospheric Tomography，CIT)可实现三维甚至四维电离层电子密度的重构。基于 GNSS的电离层层析成像技术是一种全天候、大范围的电离层探测技术，随着全球范围内GNSS永久监测网的布设和各区域、局域GNSS综合应用网的建成，对电离层进行高时空分辨率的动态监测成为现实。

当前常用的电离层层析算法主要分为函数基和网格基两类。在层析问题中，由于接收机数量、观测角度等因素的限制，函数基算法用少量参数去拟合大区域的分布将很难求解，通过施加约束得到的结果也会过于平滑，掩盖了电离层小尺度扰动。网格基算法中，待反演区域按照经纬高的方向划分成一个个网格，利用经验模型给每个网格赋初始值，再通过迭代使其接近真实值，其中迭代算法如ART(Algebra Reconstruction Technique，代数重构算法)等尤其适用于这种投影数据不完整的情况。

由于前述原因，在实际层析过程中有很多网格没有任何观测信号穿过，使得反演成为一个严重的病态问题，尤其是在待反演区域的边缘地带，相比于中心地带观测值较少，传统的ART算法已经不能很好的适用，未知数太多计算效率低下，边缘地区常出现过度的拟合，难以实现高精度高效率的电离层电子密度反演。

因此，为了克服现有技术中的问题，需要一种边缘增强的电离层层析方法。

发明内容

本发明的一个目的在于提供一种边缘增强的电离层层析方法，所述方法包括如下方法步骤：

步骤1)，利用GPS码观测值与载波相位观测值联合解算测量的电离层TEC；

步骤2)，建立虚拟参考站，利用步骤1)得到的测量的电离层TEC 值，通过三角形线性插值计算虚拟TEC值；

步骤3)，将待反演区域按照经度、纬度、高度方向均匀划分成一个个网格，

计算观测信号射线，求出每条射线在其穿过的网格中的截距，求出全部观测信号在每个网格中的截距后，按照网格和观测次数的顺序排列组成观测矩阵A；

步骤4)，采用IDW插值方法作为水平约束，使没有信号射线穿过的网格参与层析迭代；

步骤5)，构造边缘函数，通过ART迭代算法，求解电离层电子密度分布。

优选地，利用载波相位L₄组合的观测值平滑伪距的方法，解算步骤1) 中测量的电离层TEC：

设第i个历元码观测P₄组合的预报值为：

(P₄)_prd(i)＝(P₄)_sm(i-1)+[L₄(i-1)]，

则第i个历元码观测P₄组合的平滑值为：

(P₄)_sm(i)＝ω_i(P₄)_obs(i)+(1-ω_i)[(P₄)_sm(i-1)+L₄(i)-L₄(i-1)]，

当i＝1时：(P₄)_sm＝(P₄)_obs，ω_i＝1.0，

设码观测P₄组合的平滑值为P₄’，则

其中，ω_i是相应时刻对应的权重因子，DCB为码观测中仪器偏差延迟，prd下标表示预报值，sm下标表示平滑值，obs下标表示观测值，f₁和f₂分别为两种载波的频率。

优选地，所述码观测P₄组合为P₂-P₁，P₁和P₂是两种频率下的伪距观测值；

所述载波相位L₄组合为L₁-L₂，L₁和L₂是两种频率下的载波相位观测值。

优选地，通过如下方法计算虚拟TEC值：

给定三个基准站S0，S1和S2的位置坐标(x₀,y₀)，(x₁,y₁)，(x₂,y₂)，以及每个基准站对应的测量的电离层TEC值V0，V1和V2；

插值求解三个基准站形成的三角形内一点坐标S对应的虚拟TEC值 V。

优选地，步骤3)中观测信号射线通过如下方程表述：

其中，

反演时间段内的卫星坐标为(X_S,Y_S,Z_S)，地面观测站的坐标 (X_D,Y_D,Z_D)，K是射线方程的截距。

优选地，步骤3)中每条射线在其穿过的网格中的截距通过如下方式表述：

其中，X₁，X₂，Y₁，Y₂，Z₁，Z₂是相邻交点的坐标。

优选地，步骤4)中采用如下方法进行IDW插值：

设有n个点，其中第i个点的平面坐标为(x_i,y_i)，垂直高度为z_i，则IDW插值函数为：

其中，

是点(x,y)到点(x_i,y_i)的水平距离，p为幂参数。

优选地，步骤5)中边缘函数表述为：

其中，x_c，y_c是中心网格的位置，x_i，y_i是当前网格的位置，a是幂参数。

优选地，步骤5)中ART迭代算法表述为：

其中，a_i表示观测矩阵A的第i行，a_i ^T表示观测矩阵A的第i行的转置，q_i是TEC向量中的第i个值，λ为每一步迭代的松弛因子，f_e为边缘函数，x^k是电离层电子密度向量。

本发明提供的一种边缘增强的代数重构算法，针对边缘地区过度拟合的情况，提出一个边缘函数减小对边缘地区的修正，采用反距离加权算法加以约束，克服通过网格的信号路径数量不足的问题，以提高电离层层析反演精度，同时采用并行计算方法，提升计算效率。同时使用经验正交函数(Empirical Orthogonal Function，EOF)，提高垂直方向的分辨率，减少未知量来正则化病态问题。

本发明充分考虑电离层电子密度分布特性，使电离层层析反演算法本身更符合空间天气规律；使用SVD技术和IDW插值算法，减少了未知数的个数，削弱了病态性，提高了水平和垂直方向的分辨率。

本发明充分考虑实际反演中边缘地区观测信号稀缺的情况，对ART 算法的修正式加权以限制对边缘地区的修正，防止过度拟合；使用并行计算方法可大幅度提升计算效率，节省运算时间。

应当理解，前述大体的描述和后续详尽的描述均为示例性说明和解释，并不应当用作对本发明所要求保护内容的限制。

附图说明

参考随附的附图，本发明更多的目的、功能和优点将通过本发明实施方式的如下描述得以阐明，其中：

图1示意性示出了电离层层析的示意图。

图2示出了本发明边缘增强的电离层层析方法的流程图。

图3示出了本发明虚拟参考站的示意图。

图4示出了本发明插值计算虚拟TEC值的示意图。

图5示出了现有技术迭代算法观测数分布示意图。

图6示出了本发明边缘增强的电离层层析方法的ART迭代算法流程图。

具体实施方式

通过参考示范性实施例，本发明的目的和功能以及用于实现这些目的和功能的方法将得以阐明。然而，本发明并不受限于以下所公开的示范性实施例；可以通过不同形式来对其加以实现。说明书的实质仅仅是帮助相关领域技术人员综合理解本发明的具体细节。

在下文中，将参考附图描述本发明的实施例。在附图中，相同的附图标记代表相同或类似的部件，或者相同或类似的步骤。

在电离层层析技术(CIT)中，需要解算方程：TEC向量＝观测矩阵A ×电子密度X，通过解算方程求解电子密度X。

TEC定义为沿卫星传播路径S对电子密度Ne的积分

也即单位底面积沿信号传播路径贯穿整个电离层的一个柱体中所含的电子数，常用单位是(1TECU＝10¹⁶el/m²)。

如图1所示电离层层析的示意图，卫星1的卫星信号穿过电离层，地面观测站2接收卫星信号，从卫星信号获取数据解算电离层的中电子含量(TEC)构成TEC向量，构造观测矩阵A，求解上述方程电子密度 X，完成电离层层析。

为了解决现有技术中电离层层析过程中存在技术问题，如图2所示本发明边缘增强的电离层层析方法的流程图，根据本发明的实施例，一种边缘增强的电离层层析方法，包括如下方法步骤：

步骤S101，获取卫星数据，解算TEC值。

根据本发明的实施例，利用GPS码观测值(伪距观测值)与载波相位观测值联合解算测量的电离层TEC。

卫星信号通过电离层时，涉及到两种速度，一种称为群速度V_g，即不同频率的一组GPS信号作为一个整体在电离层中的传播的速度；一种称为相速度V_p，即单一频率的电磁波相位传播的速度。其相应的速度对应的折射率分别为群折射率n_g和相折射率n_p，在忽略高阶项后，其大小可简化如下式：

其中，f是载波信号的频率。当GPS卫星发射的码信号与载波相位信号穿过电离层时会产生距离延迟和相位超前。

距离延迟公式如下：

相位超前公式如下：

在电离层层析中，通常利用双频接收机的观测数据，根据两个不同频率的观测值获取高精度电离层TEC值。

伪距P1和P2的码观测方程表示为：

P₁＝ρ-cdt_k+cdt_s+I₁+T+S_p1+R_p1+M+ε，

P₂＝ρ-cdt_k+cdt_s+I₂+T+S_p2+R_p2+M+ε，

式中，P₁和P₂是两种频率下的伪距观测值；ρ是GPS卫星和观测站的几何距离；cdt_k，cdt_s分别是观测站接收机和卫星的钟差引起的测距误差；I₁和I₂是电离层延迟量；T是对流层延迟量；S_P1和S_P2是卫星的硬件偏差；R_P1和R_P2是接收机的硬件偏差；M是卫星信号传播到测站所受的多径效应的影响，ε是观测噪声。

对同一个接收机而言，二者作差得到：

P₂-P₁＝I₂-I₁+DCB，

其中，DCB称为码观测中仪器偏差延迟，其值为(R_p2-R_p1)+(S_p2-S_p1)， (S_p2-S_p1)称为卫星差分码偏差，(R_p2-R_p1)称为接收机差分码偏差。，P₂-P₁称为码观测P₄组合。

在GPS的观测过程中，多路径效应的影响主要通过采用扼流圈天线以及选择合适的高度角进行抑制。

通过伪距观测解算电离层TEC的过程表示为：

其中，f₁和f₂分别为两种载波的频率，P₂-P₁称为码观测P₄组合，DCB 称为码观测中仪器偏差延迟。

载波L₁和L₂的相位观测方程表示为：

L₁＝ρ-cdt_k+cdt^s-I₁+T-λ₁N₁+S_I1+R_I1+M+ε，

L₂＝ρ-cdt_k+cdt^s-L₂+T-λ₂N₂+S_I2+R_I2+M+ε，

式中，L₁和L₂是两种频率下的载波相位观测值；ρ是GPS卫星和测站的几何距离；cdt_k，cdt_s分别是测站接收机和卫星的钟差引起的测距误差；I₁和I₂是两种频率下的电离层延迟量；T是对流层延迟量；λ₁和λ₂是两种载波的波长；N₁和N₂是两种载波的初始整周模糊度；S₁₁和S₁₂是卫星的硬件偏差；Rl1和Rl2是接收机的硬件偏差；M是卫星信号传播到测站所受的多径效应的影响，ε是观测噪声。

对同一个接收机而言，二者作差得到：

L₂-L₁＝I₁-L₂-λ₂N₂+λ₁N₁+IFB，

其中，L₁和L₂是两种频率下的载波相位观测值，L₁-L₂称为载波相位 L₄组合，IFB称为相位观测中仪器延迟偏差，其值为(R_l2-R₁₁)+(S_l2 -S_l1)。

通过载波相位观测解算电离层TEC的过程表示为：

依据码观测值形成的P₄组合可以得到绝对TEC值，但是精度相对较低，根据载波相位L₄组合求解的TEC精度很高，但其反映的是TEC的相对变化。

本发明利用载波相位L₄组合的观测值平滑伪距的方法，解算测量的电离层TEC：

设第i个历元码观测P₄组合的预报值为：

(P₄)_prd(i)＝(P₄)_sm(i-1)+[L₄(i-1)]，

则第i个历元码观测P₄组合的平滑值为：

(P₄)_sm(i)＝ω_i(P₄)_obs(i)+(1-ω_i)[(P₄)_sm(i-1)+L₄(i)-L₄(i-1)]，

当i＝1时：(P₄)_sm＝(P₄)_obs，ω_i＝1.0，

设码观测P₄组合的平滑值为P₄'，则

其中，ω_i是相应时刻对应的权重因子，DCB为码观测中仪器偏差延迟，prd下标表示预报值，sm下标表示平滑值，obs下标表示观测值，f₁和f₂分别为两种载波的频率。

步骤S102、计算虚拟TEC值。

根据本发明的实施例，在步骤S101中，由于观测站数量和角度的限制，导致实际的观测数量有限，解算出来的实际测量的电离层TEC受限。层析所需的观测数严重不足的问题，需要丰富观测值。

根据本发明的实施例，建立虚拟参考站，利用步骤S101得到的测量的电离层TEC值，通过三角形线性插值计算虚拟TEC值。

如图3所示本发明虚拟参考站的示意图。虚拟参考站(Virtual ReferenceStation，VRS)技术是建立全球导航卫星系统连续运行参考站网 (ContinuouslyOperating Reference Stations，CORS)中最常用的技术。其核心思想是利用主基准站和周围基准站的观测信息在用户站附近的粗略位置生成虚拟的载波相位观测值，按照单基准站实时动态技术(Real Time Kinematic，RTK)计算用户精确位置。该方法的优点是精度较高且分布均匀、空间服务范围广、实时性强，在电离层层析中，虚拟参考站的设计是为了改善空间中信号射线的分布。对于每个VRS，从附近的观测站中选取多个插值站，生成虚拟观测，减少无观测信号穿过的网格数量，提高重构后电离层电子密度分布的精度。

图4示出了本发明插值计算虚拟TEC值的示意图，通过如下方法计算虚拟TEC值：

给定三个基准站S0，S1和S2的位置坐标(x₀,y₀)，(x₁,y₁)，(x₂,y₂)，以及每个基准站对应的测量的电离层TEC值V0，V1和V2；

插值求解三个基准站形成的三角形内一点坐标S对应的虚拟TEC值V。图4中，虚线线段平行于S2-S1，只需求解出S1’对应的数值V1’， S2’对应的数值V2’，子线段S1’-S占S1’-S2的比例就可以用简单的线性插值求解得到虚拟TEC值V。

设S1’-S占S1-S2的比例为t₁，S1’-S0占S0-S1的比例为t₂：

又已知各点的坐标，则有：

t₂(x₁-x₀)＝x-x₀+t₁(x₂-x₁)

t₂(y₁-y₀)＝y-y₀+t₁(y₂-y₁)，

反解即可得到t₁和t₂，进而有：

V₁'＝(1-t₂)×V₀+t₂×V₁

V₂'＝(1-t₂)×V₀+t₂×V₂，

令

步骤S103、构造观测矩阵A。

反演某个区域的电离层电子密度X，要把这个区域按照经纬高三个方向，化成一个网格，由于测站数量和角度的限制，导致实际的观测数量有限，有很多网格没有任何信号穿过，本发明通过步骤S102利用插值的方法计算得到虚拟的TEC值，有了虚拟的TEC值，该网格才能在后续的迭代计算中得到更新。

将待反演区域按照经纬高方向划分为一个个网格，由精密星历内插得到所需历元上的卫星坐标，结合已知的地面站坐标求出每条观测射线在格网中的截距长度，建立观测矩阵。

具体地，根据本发明的实施例，将待反演区域按照经度、纬度、高度方向均匀划分成一个个网格。

计算观测信号射线，求出每条射线在其穿过的网格中的截距，求出全部观测信号在每个网格中的截距后，按照网格和观测次数的顺序排列组成观测矩阵A。

实施例中，将待反演区域按照经度、纬度、高度方向均匀划分成一个个网格，并假定每个网格中的电子密度均匀分布且在反演时间段内不发生变化。

观测信号射线通过如下方程表述：

其中，

根据IGS发布的精密星历通过插值得到反演时间段内的卫星坐标为 (X_s,Y_s,Z_s)，地面观测站的坐标(X_D,Y_D,Z_D)，K是射线方程的截距。

构造观测矩阵的重点就是要求出每条射线在其穿过的网格中的截距大小，确定了射线方程以后，再确定相应的高度、经度和纬度面的方程，即可求得交点坐标。

其中，高度面是指以地球质心为圆心，以地球半径R和大地高H之和为半径的球面，方程为：

X²+Y²+Z²＝(R+H)²，

经度面是指过WGS84坐标系Z轴线的平面，精度J是该平面与零度子午线的夹角，方程为：

tanJ·X-Y＝0，

纬度面是指母线与旋转轴Z轴夹角为纬度W，顶点为地球质心的圆锥面，方程为：

Z²＝tan² W·(X²+Y²)，

得到每条射线在其穿过的网格中的截距的如下表述：

其中，X₁，X₂，Y₁，Y₂，Z₁，Z₂是相邻交点的坐标。

求出全部观测信号在每个网格中的截距后，再按照网格和观测次数的顺序排列组成矩阵，没有被信号穿过的网格截距为0。

在上文中阐述到电离层层析需要解算方程：TEC向量＝观测矩阵A×电子密度X，求解电子密度X。在电离层层析中，测量的TEC的数量远远小于未知的网格电子密度数，TEC和电子密度之间的矩阵关系中，通过几个EOF的线性组合可以显著地将未知数的数量减少为组合的几个未知系数。

本发明的一些优选的实施例中，可以从一个特定的大数据集中(如利用IRI模型生成大量的剖面)提取几个主要的正交函数，由它们的线性组合近似为电离层的垂直剖面，利用SVD方法构造EOF，降低观测矩阵和电子密度矩阵的维数。

在一些实施例中，从特定的数据矩阵(X_m×n)中提取EOF如下所示：

奇异值包含在对角阵Σ中，EOF则在U矩阵中，通过奇异值的大小评估每个EOF的贡献比来选择优势EOF。生成的电子密度剖面可以表示为一个矩阵X_m×n，其中m为垂直方向的网格数，n为由IRI模型生成的垂直剖面数。

因为电子密度剖面随当地时间变化，因此要计算每小时的EOF。将奇异值分解应用到数据矩阵后，只需要从U矩阵中选取了几个最大的EOF即可，未知数的数量将大大减少。

电离层层析技术(CIT)中，斜路径的总电子含量(STEC)可近似为路径长度与GPS信号路径上每个网格的电子密度相乘的和。

在一些实施例中，以斜路径的总电子含量(STEC)为例，利用SVD 方法构造EOF，降低观测矩阵和电子密度矩阵的维数。

观察到的第i个STEC以矩阵形式可以表示为:

k是总网格的个数，a_ij是第i条射线穿过第j个网格的截距，x_j是第 j个网格的电子密度，e_i是测量和估计误差。在该矩阵表达式中使用EOF，可以转换成：

A_ik·E_kc·E_kc ^T·X_k1＝B_i1＝STEC_i1，

H_icY_ci＝B_i1＝STEC_i1，

其中E是来自式

中U矩阵的EOF，c 是EOF的数量，H矩阵为变换后的观测矩阵，Y为变换后的电子密度矩阵。通过上述变换，信号路径的总数不变，但观测矩阵和电子密度矩阵的维数降低。

本发明采用奇异值分解法(Singular Value Decomposition，SVD)构造一组经验正交函数(Empirical Orthogonal Function，EOF)，将电离层的垂直剖面近似为标准正交函数的线性组合，减少未知数的数量。

步骤S104、采用IDW插值方法作为水平约束，使没有信号射线穿过的网格参与层析迭代。

在电离层层析问题中，由于前述的接收机数量、分布和观测角度等等原因的限制，使得划分的网格中有许多没有被观测信号穿过，没有观测值就没有截距，在后续的迭代中将无法参与，得不到任何改善，始终保持着经验模型给出的初值，误差较大。

考虑到实际中电离层电子密度应当是平滑变化的，某个网格的电子密度和周围网格的电子密度有着密切的关系，本发明反距离加权插值法 (Inverse Distance Weight，IDW)让没有信号穿过的网格也能得到改善，参与后续的迭代计算。

具体地，采用IDW方法作为水平方向上的约束，基于距输出网格的距离来控制已知网格对内插值的影响，改善无观测信号网格反演迭代中始终保持初值的情况。

IDW插值使用一组采样点的线性权重组合来确定网格元值(像元值)，权重是一种反距离函数，进行插值处理的表面应当是具有局部因变量的表面。此方法假定所映射的变量因受到与其采样位置间的距离的影响而减小，主要依赖于反距离函数的幂值。幂参数可基于距输出点的距离来控制已知点对内插值的影响，通常是一个正实数，默认值为2。

设有n个点，其中第i个点的平面坐标为(x_i,y_i)，垂直高度为z_i，采用如下方法进行IDW插值：

设有n个点，其中第i个点的平面坐标为(x_i,y_i)，垂直高度为z_i，则 IDW插值函数为：

其中，

是点(x,y)到点(x_i,y_i)的水平距离 p为幂参数。

步骤S105、构造边缘函数，通过ART迭代算法，求解电离层电子密度分布。

现有的各种常用电离层层析算法，无论是迭代重构算法、非迭代重构算法，还是卡尔曼滤波算法、神经网络算法以及模式参数拟合算法，在评价反演结果时考虑的都是反演地区的中心部分。即使存在着接收机数量、观测角度有限等问题，中心部分的网格被观测信号穿过的概率也比较大，因此观测值相对丰富，反演精度相对较高。

现有技术中ART算法的迭代式如下：

其中，a_i ^T表示观测矩阵A的第i行的转置，y_i是TEC向量中的第i 个值，λ为每一步迭代的松弛因子，x^k是电离层电子密度向量。从几何观点来看，每一步迭代相当于将电离层电子密度向量x^k向第i个方程所代表的超平面进行投影。每一步迭代对应于电离层中的一次TEC测量，即针对一个方程进行，修正的依据是利用第k次迭代计算的电离层电子密度求出的TEC与实际测量的TEC之差，对电离层电子密度分布图像做相应的修正。

以2014年12月3日14时，20°N-25°N，106°E-120°E，50-1000KM 范围内，基于香港CORS站的观测得到的半小时内各网格观测数的分布情况，如图5所示现有技术迭代算法观测数分布示意图(网格分辨率为 0.5°×1.0°×50KM)。

在边缘部分，只有依靠虚拟参考站以及约束得到的少量观测值情况下，会出现过度修正导致反演得到的电子密度远小于实际值，在峰值密度高度附近该现象尤为明显。

针对边缘拟合易过度的情况，依据网格的空间位置构造边缘函数f_e调整对不同位置网格的修正程度。

根据本发明的实施例，边缘函数表述为：

其中，x_c，y_c是中心网格的位置，x_i，y_i是当前网格的位置，a是幂参数，依据实际反演的情况取1到2之间。

据此，减小对边缘网格的修正，更新后的ART迭代算法表述为：

其中，a_i 表示观测矩阵A的第i行，a_i ^T表示观测矩阵A的第i行的转置，q_i是TEC向量中的第i个值，λ为每一步迭代的松弛因子，f_e为边缘函数，x^k是电离层电子密度向量。如图6所示本发明边缘增强的电离层层析方法的ART迭代算法流程图。

本发明提供的一种边缘增强的代数重构算法，针对边缘地区过度拟合的情况，提出一个边缘函数减小对边缘地区的修正，采用反距离加权算法加以约束，克服通过网格的信号路径数量不足的问题，以提高电离层层析反演精度，同时采用并行计算方法，提升计算效率。同时使用经验正交函数(Empirical Orthogonal Function，EOF)，提高垂直方向的分辨率，减少未知量来正则化病态问题。

本发明充分考虑电离层电子密度分布特性，使电离层层析反演算法本身更符合空间天气规律；使用SVD技术和IDW插值算法，减少了未知数的个数，削弱了病态性，提高了水平和垂直方向的分辨率。

本发明充分考虑实际反演中边缘地区观测信号稀缺的情况，对ART 算法的修正式加权以限制对边缘地区的修正，防止过度拟合；使用并行计算方法可大幅度提升计算效率，节省运算时间。

结合这里披露的本发明的说明和实践，本发明的其他实施例对于本领域技术人员都是易于想到和理解的。说明和实施例仅被认为是示例性的，本发明的真正范围和主旨均由权利要求所限定。

Claims (7)

1.一种边缘增强的电离层层析方法，其特征在于，所述方法包括如下方法步骤：

步骤1)，利用GPS码观测值与载波相位观测值联合解算测量的电离层TEC；

步骤2)，建立虚拟参考站，利用步骤1)得到的测量的电离层TEC值，通过三角形线性插值计算虚拟TEC值；

步骤3)，将待反演区域按照经度、纬度、高度方向均匀划分成一个个网格，

计算观测信号射线，求出每条射线在其穿过的网格中的截距，求出全部观测信号在每个网格中的截距后，按照网格和观测次数的顺序排列组成观测矩阵A；

步骤4)，采用IDW插值方法作为水平约束，使没有信号射线穿过的网格参与层析迭代；

步骤5)，构造边缘函数，通过ART迭代算法，求解电离层电子密度分布；

所述步骤5)中边缘函数表述为：

其中，x_c，y_c是中心网格的位置，x_i，y_i是当前网格的位置，a是幂参数；

所述步骤5)中ART迭代算法表述为：

其中，a_i表示观测矩阵A的第i行，

表示观测矩阵A的第i行的转置，q_i是TEC向量中的第i个值，λ为每一步迭代的松弛因子，f_e为边缘函数，x^k是电离层电子密度向量。

2.根据权利要求1所述的方法，其特征在于，利用载波相位L₄组合的观测值平滑伪距的方法，解算步骤1)中测量的电离层TEC：

设第i个历元码观测P₄组合的预报值为：

(P₄)_prd(i)＝(P₄)_sm(i-1)+[L₄(i-1)]，

则第i个历元码观测P₄组合的平滑值为：

(P₄)_sm(i)＝ω_i(P₄)_obs(i)+(1-ω_i)[(P₄)_sm(i-1)+L₄(i)-L₄(i-1)]，

当i＝1时：(P₄)_sm＝(P₄)_obs，ω_i＝1.0，

设码观测P₄组合的平滑值为P₄'，则

其中，ω_i是相应时刻对应的权重因子，DCB为码观测中仪器偏差延迟，prd下标表示预报值，sm下标表示平滑值，obs下标表示观测值，f1和f2分别为两种载波的频率。

3.根据权利要求2所述的方法，其特征在于，所述码观测P₄组合为P₂-P₁，P₁和P₂是两种频率下的伪距观测值；

所述载波相位L₄组合为L₁-L₂，L₁和L₂是两种频率下的载波相位观测值。

4.根据权利要求1所述的方法，其特征在于，通过如下方法计算虚拟TEC值：

给定三个基准站S0，S1和S2的位置坐标(x₀,y₀)，(x₁,y₁)，(x₂,y₂)，以及每个基准站对应的测量的电离层TEC值V0，V1和V2；

插值求解三个基准站形成的三角形内一点坐标S对应的虚拟TEC值V。

5.根据权利要求1所述的方法，其特征在于，步骤3)中观测信号射线通过如下方程表述：

其中，

反演时间段内的卫星坐标为(X_s,Y_s,Z_s)，地面观测站的坐标(X_D,Y_D,Z_D)，K是射线方程的截距。

6.根据权利要求1所述的方法，其特征在于，步骤3)中每条射线在其穿过的网格中的截距通过如下方式表述：

其中，X₁，X₂，Y₁，Y₂，Z₁，Z₂是相邻交点的坐标。

7.根据权利要求1所述的方法，其特征在于，步骤4)中采用如下方法进行IDW插值：

设有n个点，其中第i个点的平面坐标为(x_i,y_i)，垂直高度为z_i，则IDW插值函数为：

其中，

是点(x,y)到点(x_i,y_i)的水平距离，p为幂参数。

Images