CN103454695B - 一种gps电离层tec层析方法 - Google Patents

Abstract

本发明提出了一种GPS电离层TEC层析方法，所述方法通过获取GPS观测数据，对照GPS精密星历，利用拉格朗日插值方法内插卫星任意时刻坐标，进而探测并修复GPS周跳，计算GPS相位平滑伪距、TEC值和层析投影矩阵，最后进行电离层TEC层析。本发明方法基于GPS观测数据，给出的层析投影矩阵求解方法，方便、快捷、高效地实现了电离层TEC层析计算。

Description

一种GPS电离层TEC层析方法

技术领域

本发明属于大气观测技术领域，具体指的是一种GPS电离层TEC层析方法。

背景技术

电离层是地球高层大气被电离了的部分，包括从离地面约50公里开始一直伸展到约1000公里高度的地球高层大气空域。电离层是一种微弱的电离气体，其中存在相当多的自由电子和离子，能改变电磁波的传播速度，使电磁波发生折射、反射和散射，产生极化面的旋转并受到不同程度的吸收，从而对短波通讯、卫星通讯和导航定位等产生严重影响。

电离层电子浓度总含量（TotalElectronContent，TEC）又称电离层电子浓度柱含量、积分含量、总电子含量等，是单位面积内电子浓度沿卫星信号传播路径的积分。TEC是一个表征电离层物理特征的重要参量，对于电离层物理理论、电磁波传播理论及应用等研究都具有十分重要的意义。

CT是计算机体层摄影术（ComputedTomography，CT）的简称，在CT的发展过程中还曾被称为：计算机辅助断层摄影、计算机横断面体层摄影、计算机横断面轴向体层摄影等。CT技术可以在不破坏物质结构的前提下，根据物体周边所获取的物理观测量的投影数据，运用一定的数学变换，通过计算机处理，重建物体某一层面上的二维图像，并依据一系列二维图像构建物体某些部分的三维图像，该技术广泛应用于医学、海洋学、声学遥感等诸多研究领域。

GPS卫星信号在传输过程中受到电离层影响产生延迟和弯曲，TEC与GPS信号电离层延迟密切相关。因此通过一定的数学模型，利用GPS卫星信号电离层延迟信息反演电离层TEC，并进一步建立一定观测区域内电离层TEC三维图像，即是GPS电离层TEC层析技术。电离层层析技术是一种全天候、大范围的电离层探测技术，对于监测和研究电离层不同尺度不均匀性、电离层环境、电离层时空变化、电离层异常变化、建立电离层预报模型以及电离层全范围的监测均具有重要的科学意义和应用价值。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足，提出一种GPS电离层TEC层析方法。所述方法针对传统GPS电离层TEC层析方法计算繁琐、不利于编程的问题，给出了层析投影矩阵的新的求解方法，实现了基于GPS观测数据方便、快捷、高效地进行电离层TEC层析计算。

为了解决上述技术问题，本发明所采用的技术方案是：

一种GPS电离层TEC层析方法，具体步骤如下：

步骤1，获取GPS观测数据；

步骤2，对照GPS精密星历，利用拉格朗日插值方法内插卫星任意时刻坐标；

步骤3，探测并修复GPS周跳；

步骤4，计算GPS相位平滑伪距；

步骤5，计算电离层电子浓度总含量TEC值；

步骤6，由卫星和测站坐标，计算层析投影矩阵A，A的表达式为：

$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \·\·\·& a_{1n}\\ a_{21}& a_{22}& \·\·\·& a_{2n}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ a_{m1}& a_{m2}& \·\·\·& a_{\mathrm{mn}}\end{array}\right]$

其中，a_ij是第i条GPS信号沿倾斜路径穿过第j个层析网格的射线长度，i={1,2,...,m},m为测站卫星观测个数；j={1,2,...,n}，n为总的网格数；

步骤7，根据TEC值和层析投影矩阵，进行电离层TEC层析，TEC层析方程为：

TEC=A·X+ξ

其中，TEC=[TEC₁,TEC₂,...,TEC_i,...,TEC_m]^T，X=[ρ₁,ρ₂,...,ρ_i,...,ρ_m]^T，ξ=[ξ₁,ξ₂,...ξ_i,...,ξ_m]^T；TEC_i为第i个卫星对应的TEC值，ξ_i为TEC_i误差，ρ_i是第i个网格内的电子密度,[·]^T为矩阵转置运算。

本发明的有益效果是：本发明提出了一种GPS电离层TEC层析方法，所述方法通过获取GPS观测数据，对照GPS精密星历，利用拉格朗日插值方法内插卫星任意时刻坐标，进而探测并修复GPS周跳，计算GPS相位平滑伪距、TEC值和层析投影矩阵，最后进行电离层TEC层析。本发明方法基于GPS观测数据，给出的层析投影矩阵求解方法，方便、快捷、高效地实现了电离层TEC层析计算。

附图说明

图1是GPS电离层TEC层析流程图。

图2是测站与卫星位置示意图。

图3是投影矩阵计算流程图。

具体实施方式

下面结合附图，对本发明提出的一种GPS电离层TEC层析方法进行详细说明：

如图1所示，一种GPS电离层TEC层析方法包含下面7个步骤：

1、下载GPS数据文件，主要包括：观测数据文件（O文件），并将它们转化为Rinex格式；从IGS网站下载对应于观测时段的SP3精密星历文件。

2、利用拉格朗日插值方法内插卫星任意时刻坐标。

3、GPS周跳的探测与修复。

4、GPS相位平滑伪距。

5、计算电离层电子浓度总含量TEC值。

6、计算层析投影矩阵；

6.1确定层析网格的范围

设定层析网格的范围为：经度范围(L₀,L₁)，纬度范围(B₀,B₁)，经度与纬度网格的间隔范围为G°；高度范围(H₀,H₁)，高度的间隔为gh;

设测站S的坐标为：(X_s,Y_s,Z_s)，t时刻GPS卫星j的坐标为：所求测站与GPS卫星连线同层析网格的交点坐标为：

直线方程：

$\frac{X_s^j\left(t\right)-X_s}{X_t^j-X_s}=\frac{Y_s^j\left(t\right)-Y_s}{Y_t^j-Y_s}=\frac{Z_s^j\left(t\right)-Z_s}{Z_t^j-Z_s}=k---\left(1\right)$

这里，令 $\mathrm{\Δ}{X}_s^j=X_t^j-X_s,\mathrm{\Δ}{Y}_s^j=Y_t^j-Y,\mathrm{\Δ}{Z}_s^j=Z_t^j-Z.$

空间极坐标形式：

$\left\{\begin{array}{c}{R}_s^j\left(t\right)=\sqrt{X_s^j{\left(t\right)}^2+Y_s^j{\left(t\right)}^2+Z_s^j{\left(t\right)}^2}\\ L_s^j\left(t\right)=\arctan \frac{Y_s^j\left(t\right)}{X_s^j\left(t\right)}\\ B_s^j\left(t\right)=\arcsin \frac{Z_s^j\left(t\right)}{R_s^j\left(t\right)}\end{array}\right.---\left(2\right)$

6.2求测站与GPS卫星连线与不同高度面的交点

设同测站S与GPS卫星j连线相交的高度面的高度为h，则式（2）中的取值为：

$R_s^j\left(h\right)=h+r---\left(3\right)$

其中，r为WGS-84椭球的长半径。

令：

$\left\{\begin{array}{c}{X}_s^j\left(t\right)=\mathrm{\Δ}{X}_s^j\·k+X_s\\ Y_s^j\left(t\right)=\mathrm{\Δ}{Y}_s^j\·k+Y_s\\ Z_s^j\left(t\right)=\mathrm{\Δ}{Z}_s^j\·k+Z_s\end{array}\right.---\left(4\right)$

将（3）、（4）代入（2）式，变换为一元二次方程的形式：

a·k²+b·k+c=0（5）

这里， $a=\mathrm{\Δ}{X}_s^{j2}+\mathrm{\Δ}{Y}_s^{j2}+\mathrm{\Δ}{Z}_s^{j2},$ $b=2(\mathrm{\Δ}{X}_s^j\·X_s+\mathrm{\Δ}{Y}_s^j\·Y_s+\mathrm{\Δ}{Z}_s^j\·Z_s),$ $c=(X_s^2+Y_s^2+Z_s^2)-R_s^j{\left(h\right)}^2.$

解（5）式，得到方程的两个解k₁、k₂，对应于k₁、k₂的交点坐标为：卫星j与交点A、C之间的距离为d_k1、d_k2。

如图2所示，图2中，外圈大圆表示卫星j的运行轨道，S表示测站位置，内圈虚线圆表示高度为h的层析网格面，直线AB表示测站S与卫星j连线同高度面h的交点，O-XYZ表示空间直角坐标系，（L,B,R）表示交点A的极坐标（地理坐标）。

若交点处于测站与卫星之间，则要满足d_k1<d_k2,由此确定k的唯一值，其对应的交点坐标为：转化为地理坐标

6.3确定同测站与卫星连线相交的经、纬度网格线

将测站S的坐标(X_s,Y_s,Z_s)、卫星j的坐标(X^j,Y^j,Z^j)分别转化为相对应的地理坐标，记为：(L_s,B_s,H_s)，(L^j,B^j,H^j)。

将L_s、L^R、L₀、L₁表示为如式（6）所示的矩阵形式：

[L_s,L^R,L₀,L₁]（6）

将式（6）中矩阵元素按从小到大进行排序，得到：

[L₁,L₂,L₃,L₄]（7）

设测站S与卫星j连线与经度网格相交的范围为(L_min,L_max)。根据式（7），可以取L_min=L₂，L_max=L₃。同理可以求出纬度范围(B_min,B_max)。

根据网格的间隔范围G°，求出(L_min,L_max)及(B_min,B_max)范围内存在的经、纬度网格及其度数。

6.4求测站与GPS卫星连线与经度网格的交点

设测站S与卫星j连线与经度网格（度数为L）的交点坐标为由式（2）得：

$\tan L=\frac{Y_s^j\left(l\right)}{X_s^j\left(l\right)}=\mathrm{mk}---\left(8\right)$

将式（8）代入式（1），得：

$X_s^j\left(l\right)=\frac{\mathrm{\&Delta;}{Y}_s^j\&CenterDot;X_s-\mathrm{\&Delta;}{X}_s^j\&CenterDot;Y_s}{\mathrm{\&Delta;}{Y}_s^j-\mathrm{mk}\&CenterDot;\mathrm{\&Delta;}{X}_s^j}---\left(9\right)$

$Y_s^j\left(l\right)=\mathrm{mk}\&CenterDot;X_s^j\left(l\right)---\left(10\right)$

$Z_s^j\left(l\right)=\frac{X_s^j\left(l\right)-X_s}{\mathrm{\&Delta;}{X}_s^j}---\left(11\right)$

将交点坐标转化为该点相对应地理坐标存储于计算机中。

6.5求测站与GPS卫星连线与纬度网格的交点

此时，纬度网格为已知值，令其值为：B。测站S与卫星j连线与度数为B的纬度网格的交点坐标为由式（2）得：

$R_s^j\left(b\right)=\frac{Z_s^j\left(b\right)}{\sin B}---\left(12\right)$

将式（12）代入式（1）变换为如式（5）所示的一元二次方程形式，按6.2中所述方法解方程、确定参数k的唯一值并求出测站S与卫星j连线与纬度网格的交点坐标为将交点转化为相对应的地理坐标

6.6剔除不属于投影范围内的交点

根据式（13）条件及各个交点的地理坐标(L,B,H)剔除不属于投影范围内的交点。

$\left\{\begin{array}{c}{L}_0\&le;L\&le;L_1\\ B_0\&le;B\&le;B_1\\ H_0\&le;H\&le;H_1\end{array}\right.---\left(13\right)$

6.7排序

由各个交点的空间直角坐标及测站S的坐标(X_s,Y_s,Z_s)，计算各个交点同测站S之间的距离，并根据距离值由小到大进行排列。

6.8计算网格内GPS信号的路径长度

由各个交点的地理坐标(L,B,H)，根据式（14）计算交点所在的网格位置。

$\left\{\begin{array}{c}{L}_i=\mathrm{fix}[(L-L_0)/G]\\ B_i=\mathrm{fix}[(B-B_0)/G]\\ H_i=\mathrm{fix}[(H-H_0)/\mathrm{gh}]\end{array}\right.---\left(14\right)$

式中，L_i为网格位置的经度编号，B_i为网格位置的纬度编号，H_i为网格位置的高度编号；函数fix()表示朝零方向取整。

根据6.7的排序结果，求网格间的距离：

$d_i=\sqrt{{(X(i+1)-X\left(i\right))}^2+{(Y(i+1)-Y\left(i\right))}^2+{(Z(i+1)-Z\left(i\right))}^2}---\left(15\right)$

A_j（L_i,B_i,H_i）=d_i（16）

这里(X(i+1),Y(i+1),Z(i+1))，(X(i),Y(i),Z(i))表示经过6.7排序后相邻两点的坐标。

6.9组成投影矩阵

将矩阵A_j中元素按行进行排列成A′_j，并将同一历元中卫星相对应的行矩阵排序组成投影矩阵。

$A=\left[\begin{array}{c}{A}_1^{\&prime;}\\ A_2^{\&prime;}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ A_j^{\&prime;}\end{array}\right]---\left(17\right)$

步骤6.1-6.9的实施流程参见图3。

7、根据TEC值和层析投影矩阵，进行电离层TEC层析。

根据第5步计算的TEC值及第6步计算的层析投影矩阵，组成式（18）：

TEC=A·X+ξ（18）

其中，TEC=[TEC₁,TEC₂,...,TEC_i,...,TEC_m]^T，X=[ρ₁,ρ₂,...,ρ_i,...,ρ_n]^T，ξ=[ξ₁,ξ₂,...ξ_i,...,ξ_m]^T；TEC_i为第i个卫星对应的TEC值，ξ_i为TEC_i误差，ρ_i是第i个网格内的电子密度,[·]^T为矩阵转置运算；m为测站卫星观测个数；n为总的网格数。

解式（18），即可得到TEC层析结果。

Claims (1)

1.一种GPS电离层TEC层析方法，其特征在于，具体步骤如下：

步骤1，获取GPS观测数据；

步骤2，对照GPS精密星历，利用拉格朗日插值方法内插卫星任意时刻坐标；

步骤3，探测并修复GPS周跳；

步骤4，计算GPS相位平滑伪距；

步骤5，计算电离层电子浓度总含量TEC值；

步骤6，由卫星和测站坐标，计算层析投影矩阵A，A的表达式为：

其中，a_mn是第m条GPS信号沿倾斜路径穿过第n个层析网格的射线长度，m为测站卫星观测个数；n为总的网格数；层析投影矩阵A的计算过程步骤如下：

步骤6.1，确定层析网格的范围；

设定层析网格的范围为：经度范围(L₀,L₁)，纬度范围(B₀,B₁)，经度与纬度网格的间隔范围为G°；高度范围(H₀,H₁)，高度的间隔为gh；

设测站S的空间直角坐标为：(X_s,Y_s,Z_s)，t时刻GPS卫星j的空间直角坐标为：测站S与GPS卫星j连线与层析网格交点的空间直角坐标为：(X_s,Y_s,Z_s)、和三点满足直线方程：

设定：

将用空间极坐标表示，这里

步骤6.2，求测站与GPS卫星连线与不同高度面的交点；

设测站S与GPS卫星j连线相交的高度面的高度为h，其交点坐标为则的取值为：r为WGS-84椭球的长半径；

根据直线方程设

将用一元二次方程a·k²+b·k+c＝0进行表示；

方程中，

解方程a·k²+b·k+c＝0，得到k的两个解k₁、k₂；对应于k₁、k₂交点的空间直角坐标为：卫星j与交点k₁、k₂之间的距离为d_k1、d_k2 ；

如果交点位于测站S和卫星j之间，需要满足d_k1<d_k2,由此确定k的唯一值，其对应交点的空间直角坐标为：

将转化为相应的地理坐标

步骤6.3，确定同测站与卫星连线相交的经、纬度网格线；

将测站S的空间直角坐标(X_s,Y_s,Z_s)、卫星j的空间直角坐标(X^j,Y^j,Z^j)分别转化为相对应的地理坐标(L_s,B_s,H_s)和(L^j,B^j,H^j)。

将L_s、L^R、L₀、L₁表示为矩阵形式：[L_s,L^R,L₀,L₁]，L^R表示测站S与卫星j连线与不同高度面的交点地理坐标的经度值；并将该矩阵元素按从小到大进行排序，得到：[L′₁,L′₂,L′₃,L′₄]；

设测站S与卫星j连线与经度网格相交的范围为(L_min,L_max)；取L_min＝L′₂，L_max＝L′₃；同理求出纬度范围(B_min,B_max)；

根据网格的间隔范围G°，求出(L_min,L_max)及(B_min,B_max)范围内存在的经、纬度网格及其度数；

步骤6.4，求测站与GPS卫星连线与经度网格的交点；

设测站S与卫星j连线与经度网格交点的空间直角坐标为所述网格度数为W，则和W满足条件：

根据和通过可以计算得到的值，并将其转化为相对应的地理坐标

步骤6.5，求测站与GPS卫星连线与纬度网格的交点；

此时，纬度网格值为B，测站S与卫星j连线与度数为B的纬度网格交点的空间直角坐标为对应的极坐标为和B满足关系式：

按步骤6.2中所述方法确定参数k的唯一值并求出测站S与卫星j连线与纬度网格的交点的空间直角坐标将其为对应的地理坐标

步骤6.6，剔除不属于投影范围内的交点；

依据投影范围和各个交点的地理坐标(L,B,H)，剔除不属于投影范围内的交点；

步骤6.7，排序；

由各个交点的空间直角坐标及测站S的空间直角坐标(X_s,Y_s,Z_s)，计算各个交点同测站S之间的距离，并根据距离值由小到大进行排列；

步骤6.8，计算网格内GPS信号的路径长度；

根据各个交点的地理坐标(L,B,H)，计算交点所在的网格位置；其计算公式为：

式中，L_i为网格位置的经度编号，B_i为网格位置的纬度编号，H_i为网格位置的高度编号；函数fix()表示朝零方向取整；

根据步骤6.7的排序结果，求网格间的距离d_i：

将距离d_i存入A_j的相应位置A_j(L_i,B_i,H_i)＝d_i,A_j表示卫星j的网格距离矩阵，i表示第i个层析网格；A_j(L_i,B_i,H_i)表示卫星j的第i个网格距离矩阵元素；

这里(X(i+1),Y(i+1),Z(i+1))，(X(i),Y(i),Z(i))表示经过步骤6.7排序后相邻两点的空间直角坐标；

步骤6.9，组成投影矩阵；

将矩阵A_j中元素按行进行排列成A′_j，并将同一历元中卫星相对应的行矩阵排序组成投影矩阵A，

步骤7，根据TEC值和层析投影矩阵，进行电离层TEC层析，TEC层析方程为：TEC＝A·X+ξ；

其中，TEC＝[TEC₁,TEC₂,...,TEC_i,...,TEC_m]^T，X＝[ρ₁,ρ₂,...,ρ_i,...,ρ_m]^T，ξ＝[ξ₁,ξ₂,...ξ_i,...,ξ_m]^T；TEC_i为第i个卫星对应的TEC值，ξ_i为TEC_i误差，ρ_i是第i个网格内的电子密度,[i]^T为矩阵转置运算。