CN102147247B - 通过提取gps信号重力频移确定海拔高程的方法及其装置 - Google Patents

Abstract

本发明涉及一种确定海拔高程的方法及其装置，尤其是涉及一种通过提取GPS信号重力频移确定海拔高程的方法及其装置。基于IGS提供的GPS观测数据(obs文件和sp3文件)，可提取重力频移，并基于重力频移信号进一步确定单站重力位及海拔高。此外，可以通过重力场模型计算重力频移的理论值，验证所提取的重力频移信号的有效性。因此，本发明具有如下优点：1.实施相距较远的两点之间的水准测量不需要通过若干个(或多个)中间测站不断传递，节省了人力物力资源，并且施测周期短，测量精准；2.能在山区或近距离高程差非常大的两点之间实施海拔高程差测量；3.能够实现全球统一的高程基准。

Description

通过提取GPS信号重力频移确定海拔高程的方法及其装置

技术领域

本发明涉及一种确定海拔高程的方法及其装置，尤其是涉及一种通过提取GPS信号重力频移确定海拔高程的方法及其装置。

背景技术

精确确定地面点的海拔高程(简称海拔高)是地球科学领域的一项重要任务。海拔高是沿重力铅垂线至大地水准面的距离。按通常的定义，大地水准面是最接近于平均海水面的重力等位面(Listing，1872，Grafarend，199obs)。精确提供地面点的海拔高对研究与地形变化相应的全球气候环流问题、各类工程建设、地球空间信息服务等具有重要的科学价值和应用价值。

利用传统的方法确定海拔高需要采用重力测量联合水准测量的方法来实现，但存在如下缺陷：一是实施相距较远的两点之间的水准测量需要通过若干个(或多个)中间测站不断传递，不仅极大地消耗人力物力资源，施测周期长，而且测量误差随着传递路线的增长而增加；二是在山区(或近距离高程差非常大的两点之间)难于实施水准测量；三是难以实现全球统一高程基准，因为平均海水面并非等位面，两个分离的大陆难以实现(高精度的)水准联接。

为了克服上述方法的缺陷，Bjerhamma(1985)基于广义相对论提出了相对论大地测量的概念[1]，通过移动精密时钟来测定地面两点之间的重力位差(简称时钟移动法)，进而将重力位差转换为海拔高。1993年，申文斌等[2]提出了与Bjerhammar方法等价的利用电磁波信号的重力频移确定地面上两点之间的重力位差的设想(简称频移法)；Brumberg & Groten(2obsobs2)[3]也研究了类似的方法。时钟法和频移法都是广义相对论的推论。由于电磁波信号具有能量，它从一点传播到另一点时会产生能量得失：从具有较高重力位的地方向具有较低重力位的地方传播时，能量增加，因而频率增大；反之，则频率减小(这里采用物理学中的定义，因而重力位取负值，与物理大地测量中的定义正好反一个正负号)。这种信号频率的移动与重力位差之间具有确定的关系，即重力位频移方程，简称重力频移。重力频移方程的有效性已被很多物理学实验所证实。频移法比时钟移动法有更高的可行性和技术优势，因为可直接利用GPS信号进行重力频移观测。

基于重力频移确定重力位差(或高程差)的精度取决于测定频移的精度，或者说取决于守时固有频率的精度。近年来，时间与频率科学发展迅速，国际上已制造出守时稳定度为10^-16量级的时钟(Heavner et al.，2005)[5]，这相当于利用频移法单次测定海拔高程的精度可以达到米级。如果施以多次观测，可望达到更高的精度。国内外有好几个研究小组一直在致力于“光频标”的研制(Reichert et al.，2000；Diddams et al.，2001；Ma et al.，2004)，并已取得重要突破，他们相互比对了各自的“光频标”，其稳定度已达到10^-18量级。科学家预言，在不远的将来，即可实现稳定度为10^-18量级的“光钟”，这为利用频移法测定厘米级重力位差或海拔高提供了基础。

发明内容

本发明主要是解决现有技术所存在的实施相距较远的两点之间的水准测量需要通过若干个(或多个)中间测站不断传递，不仅极大地消耗人力物力资源，施测周期长，而且测量误差随着传递路线的增长而增加等的技术问题；提供了一种实施相距较远的两点之间的水准测量不需要通过若干个(或多个)中间测站不断传递，节省了人力物力资源，并且施测周期短，测量精准的一种通过提取GPS信号重力频移确定海拔高程的方法及其装置。

本发明还有一目的是解决现有技术所存在的在山区或近距离高程差非常大的两点之间难于实施水准测量等的技术问题；提供了一种能在山区或近距离高程差非常大的两点之间实施海拔高程差测量的一种通过提取GPS信号重力频移确定海拔高程的方法及其装置。

本发明再有一目的是解决现有技术所存在的难以实现全球统一高程基准，因为平均海水面并非等位面，两个分离的大陆难以实现高精度的水准联接等的技术问题；提供了一种能够实现全球统一的高程基准的一种通过提取GPS信号重力频移确定海拔高程的方法及其装置。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，包括以下步骤：

步骤1，由参数配置模块根据用户需求配置各分析参数，包括设置低阶重力位模型，经典多普勒频移模型，对流层模型，牛顿内插阶次以及微分时间间隔；

步骤2，由GPS参数读取模块从GPS文件中读入GPS观测文件，包括obs文件和精密星历sp3文件，并获取总频移量；

步骤3，由卫星定位模块根据精密星历通过牛顿内插确定卫星位置，牛顿内插选定步骤1中的内插阶次，并基于低阶重力场模型求卫星引力位；

步骤4，由旋转改正模块实施地球旋转改正；

步骤5，由地面站重力位确定模块根据步骤2中的总频移量确定地面站重力位；

步骤6，由地面站海拔高确定模块以大地水准面上的重力位为基准，基于布隆斯公式，用地面站重力位确定地面站海拔高程。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤5中，确定地面站重力位的具体步骤如下：

步骤5.1，由经典多普勒频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出经典多普勒频移值；

步骤5.2，由对流层频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出对流层频移值；

步骤5.3，由电离层频移确定单元通过obs文件中L1和L2载波的频移值，利用双频消除法计算电离层频移；

步骤5.4，由重力频移确定单元根据步骤2中的总频移量去掉步骤5.1、步骤5.2以及步骤5.3中的各频移值，获取到重力频移值；

步骤5.5，由地面站重力位确定单元根据步骤5.4中获得的重力频移以及步骤3所述的通过低阶重力场模型求得的卫星引力位，基于重力频移方程，确定地面站重力位。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤2中，总频移量Δf可表示为

Δf＝Δf_d+Δf_V+Δf_h+Δf_i+Δf_t+Δf_c+Δf_r

其中，Δf_d和Δf_k分别是经典多普勒频移和横向经典多普勒频移，Δf_V是引力频移，Δf和Δf_t分别是电离层和对流层影响，Δf_c是原子钟误差，Δf_r是随机误差，引力频移与横向经典多普勒频移之和是重力频移，记为Δf_g。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤5.1中，经典多普勒频移是由于卫星与地面观测站有相对速度导致地面站接收到的卫星信号的频率的变化，具体公式为：

$\mathrm{\Δf}=f-f_s=-\frac{f_s}{c}\frac{\mathrm{ds}}{\mathrm{dt}}$

其中c是真空中光速，f_s和f分别为发射频率和接收频率，单位时间内测站与卫星间的距离变化为ds/dt。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤3中，牛顿法插值通过构造各阶差商产生多项式系数，以待求点和插值点几何位置关系为自变量形成多项式，具体公式如下：

差商公式：

$f[x_0,x_1...x_n]=\frac{f[x_0,x_1...x_{n-2},x_n]-f[x_0,x_1...x_{n-2},x_{n-1}]}{x_n-x_{n-1}}$

牛顿插值公式：

f(x)＝f(x₀)+f[x₀，x₁](x-x₀)

      +f[x₀，x₁，x₂](x-x₀)(x-x₁)+...

      +f[x₀，x₁...x_n](x-x₀)(x-x₁)...(x-x_n-1)

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤4中，地球旋转改正即坐标变换，由如下方程给出

X_t＝Xcosa+Ysina

Y_t＝-Xsina+Ycosa

式中，a＝2πt/(3600×24)，t为卫星信号从卫星传送到地面站所花费的时间，(X，Y)是卫星在惯性系中的位置，(X_t，Y_t)是卫星在地固系中的位置。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤5.2中，对流层频移是GPS信号在传播过程中的对流层延迟影响，可以用对流层模型差分获得，可选模型有(1)SAASTMOINEN；(2)MODIFIELDHOPIFIELD；(3)SIMPLIFIED HOPIFIELD。计算对流层相位距离延迟的简化的Hopifield模型公式如下：

${\mathrm{\Δ\ρ}}_{\mathrm{trop}}={\mathrm{\ΔS}}_d+{\mathrm{\ΔS}}_w=\frac{K_d}{\sin {(E^2+6.25)}^{\frac{1}{2}}}+\frac{K_w}{\sin {(E^2+2.25)}^{\frac{1}{2}}}$

其中

$K_d=155.2\×{10}^{-7}\frac{P_s}{T_s}(h_d-h_s),$ $K_w=155.2\×{10}^{-7}\·\frac{4810}{{T_S}^2}{e}_s(h_w-h_s)$

式中，Δρ_trop为相位延迟，E为高度角，h_d＝401386+148.72(T-273.16)，h_w＝11000m，h_s是地面站高程，以米为单位，P_s是气压，以毫巴为单位；

其中

e_s＝RH·exp(-37.2465+0.21366T-0.0002568908T²)，T_s为干温，T是观测时刻的温度，采用绝对温标，单位K，RH为相对湿度。在没有测定测站气象元素的情况下，采用给定的标准气压参数计算，其中P_s取1013.25mBar，T_s取291.15K，RH取50。

根据以上模型可以求出时刻t₁的对流层相位距离延迟Δρ_trop(t₁)和时刻t₂的对流层延迟Δρ_trop(t₂)，最后能得到对流层频移的改正公式：

${\mathrm{\Δf}}_{\mathrm{trop}}=-\frac{f}{c}\frac{{\∂\mathrm{\Δ\ρ}}_{\mathrm{trop}}}{\∂t}\≈-\frac{f}{c}\frac{{\mathrm{\Δ\ρ}}_{\mathrm{trop}2}-{\mathrm{\Δ\ρ}}_{\mathrm{trop}1}}{t_2-t_1}=-\frac{{\mathrm{\Δ\ρ}}_{\mathrm{trop}2}-{\mathrm{\Δ\ρ}}_{\mathrm{trop}1}}{\left(t_2-t_1\right)c}f.$

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤5.3中，电离层频移的影响可用双频消除法消去，这是由于GPS信号的电离层频移与频率本身成反比，而其他所有频移与频率本身成正比，如下式：

$f_j^{\mathrm{obs}}=\mathrm{\α}{\tilde{f}}_j+\frac{\mathrm{\χ}}{{\tilde{f}}_j},$ j＝1，2

$\left\{\begin{array}{c}\mathrm{\α}=\frac{k(k\stackrel{\‾}{f_1^{\mathrm{obs}}}-\stackrel{\‾}{f_2^{\mathrm{obs}}})}{{\tilde{f}}_1(k^2-1)}\\ \mathrm{\χ}=\frac{\stackrel{\‾}{f_1^{\mathrm{obs}}}-k\stackrel{\‾}{f_2^{\mathrm{obs}}}}{1-k^2}{\tilde{f}}_1\end{array}\right.$

其中，χ电离层频移项，α为其它频移项之和，k＝f₁/f₂，信号的真实发射频率分别用

和

表示。通过L1和L2载波上的D1和D2观测值列出方程组，可以解出α和χ；根据α可提取重力频移信号。

在上述的通过提取GPS信号重力频移确定海拔高程的方法，所述的步骤5.5中，所述的重力位频移方程，由下式给出：

$W-W_S=c^2\frac{\mathrm{\Δf}}{f_0}$

式中，f₀是固有发射频率，Δf是重力频移，卫星重力位W_S可由低阶重力场模型求取。低阶重力场模型W_S是输入到程序中的已知量，它是低阶引力位V和离心力位Q之和：W_S＝V+Q。

卫星低阶引力位，即取8阶的获取方法如下：

$V(r_t,{\mathrm{\θ}}_t,{\mathrm{\λ}}_t)=\frac{\mathrm{GM}}{r_t}[1+\underset{n=1}{\overset{8}{\mathrm{\Σ}}}{\left(\frac{a}{r_t}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{nm}}\cos m{\mathrm{\λ}}_t+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin m{\mathrm{\λ}}_t){\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\cos {\mathrm{\θ}}_t\right)]$

其中，GM是地心常数(已知)，a＝6378136.572m是地球赤道平均半径，P_nm(cosθ_t)是正规化缔合Legendre函数，C_nm和S_nm是正规化球谐系数，是已知值，该值由EGM2008模型给出，(r_t，θ_t，λ_t)是卫星在t时刻的球面坐标，其中r_t，θ_t，λ_t分别是卫星在t时刻的站星距离，余纬和经度。

离心力位的获取方法如下：

$Q=\frac{{\mathrm{\ω}}^2}{2}{\left(r_t\sin {\mathrm{\θ}}_t\right)}^2$

其中，ω＝7.292115×10^-5rad/s是地球自转角速率。

一种利用通过提取GPS信号重力频移确定海拔高程的方法的装置，其特征在于，包括依次连接的参数配置模块、GPS参数读取模块、卫星定位模块、旋转改正模块、地面站重力位确定模块以及地面站海拔高确定模块，所述的地面站重力位确定模块包括依次连接的经典多普勒频移确定单元、对流层频移确定单元、电离层频移确定单元、重力频移确定单元以及地面站重力位确定单元。

因此，本发明具有如下优点：1.实施相距较远的两点之间的水准测量不需要通过若干个(或多个)中间测站不断传递，节省了人力物力资源，并且施测周期短，测量精准；2.能在山区或近距离高程差非常大的两点之间实施海拔高程差测量；3.能够实现全球统一的高程基准。

附图说明

图1是本发明的工作流程图；

图2是本发明的重力频移模型值；

图3为本发明的重力频移计算值。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

一种通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，包括以下步骤：

步骤1，由参数配置模块根据用户需求配置各分析参数，包括设置低阶重力位模型，经典多普勒频移模型，对流层模型，牛顿内插阶次以及微分时间间隔；

步骤2，由GPS参数读取模块从GPS文件中读入GPS观测文件，包括obs文件和精密星历sp3文件，并获取总频移量；

步骤3，由卫星定位模块根据精密星历通过牛顿内插确定卫星位置，牛顿内插选定步骤1中的内插阶次，并基于低阶重力场模型求卫星引力位；

步骤4，由旋转改正模块实施地球旋转改正；

步骤5，由地面站重力位确定模块根据步骤2中的总频移量确定地面站重力位；

步骤5.1，由经典多普勒频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出经典多普勒频移值；

步骤5.2，由对流层频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出对流层频移值；

步骤5.3，由电离层频移确定单元通过obs文件中L1和L2载波的频移值，利用双频消除法计算电离层频移；

步骤5.4，由重力频移确定单元根据步骤2中的总频移量去掉步骤5.1、步骤5.2以及步骤5.3中的各频移值，获取到重力频移值；

步骤5.5，由地面站重力位确定单元根据步骤5.4中获得的重力频移以及步骤3所述的通过低阶重力场模型求得的卫星引力位，基于重力频移方程，确定地面站重力位。

步骤6，由地面站海拔高确定模块以大地水准面上的重力位为基准，基于布隆斯公式，用地面站重力位确定地面站海拔高程。

下面介绍一下本发明所涉及到的理论基础：

GPS信号频移由多部分组成。考虑引起信号频移的主要原因，总频移观测量Δf可表示为

Δf＝Δf_d+Δf_V+Δf_h+Δf_i+Δf_t+Δf_c+Δf_r         (1)

其中，Δf_d和Δf_k分别是经典多普勒频移和横向经典多普勒频移，Δf_V是引力频移，Δf_i和Δf_t分别是电离层和对流层影响，Δf_c是原子钟误差，Δf_r是随机误差。引力频移与横向经典多普勒频移之和是重力频移，记为Δf_g。重力频移是目标信息，而其它频移是要尽量消除或减弱的对象。

经典多普勒频移是由于卫星与地面观测站有相对速度导致地面站接收到的卫星信号的频率的变化，具体公式为：

$\mathrm{\Δf}=f-f_s=-\frac{f_s}{c}\frac{\mathrm{ds}}{\mathrm{dt}}---\left(2\right)$

其中c是真空中光速，其中f_s和f分别为发射频率和接收频率，单位时间内观测站与卫星间的距离变化为ds/dt。

地面站的坐标由obs文件的头文件中提供，卫星坐标由sp3精密星历用牛顿内插法(取适当阶次)得到。牛顿法插值通过构造各阶差商产生多项式系数，以待求点和插值点几何位置关系为自变量形成多项式，计算速度快，效率高。具体公式如下：

差商公式：

$f[x_0,x_1...x_n]=\frac{f[x_0,x_1...x_{n-2},x_n]-f[x_0,x_1...x_{n-2},x_{n-1}]}{x_n-x_{n-1}}---\left(3\right)$

牛顿插值公式：

f(x)＝f(x₀)+f[x₀，x₁](x-x₀)

     +f[x₀，x₁，x₂](x-x₀)(x-x₁)+...

      +f[x₀，x₁...x_n](x-x₀)(x-x₁)...(x-x_n-1)       (4)

卫星的位置是在惯性坐标系中给出的；地面站坐标则通常在地固坐标系中表示。因此，需要做地球自转改正，实际上是坐标变换。坐标变换由如下方程给出

X_t＝Xcosa+Ysina    (5)

Y_t＝-Xsina+Ycosa

式中，a＝2πt/(3600×24)，t为卫星信号从卫星传送到地面站所花费的时间，(X，Y)是卫星在惯性系中的位置，(X_t，Y_t)是卫星在地固系中的位置。

对流层频移是GPS信号在传播过程中的对流层延迟的影响，可以用对流层模型差分获得，可选模型有(1)SAASTMOINEN；(2)MODIFIELDHOPIFIELD；(3)SIMPLIFIED HOPIFIELD。计算对流层相位距离延迟的简化的Hopifield模型公式如下：

${\mathrm{\Δ\ρ}}_{\mathrm{trop}}={\mathrm{\ΔS}}_d+{\mathrm{\ΔS}}_w=\frac{K_d}{\sin {(E^2+6.25)}^{\frac{1}{2}}}+\frac{K_w}{\sin {(E^2+2.25)}^{\frac{1}{2}}}---\left(6\right)$

其中

$K_d=155.2\×{10}^{-7}\frac{P_s}{T_s}(h_d-h_s),K_w=155.2\×{10}^{-7}\·\frac{4810}{{T_S}^2}{e}_s(h_w-h_s)---\left(7\right)$

式中，Δρ_trop为相位延迟，E为高度角，h_d＝401386+148.72(T-273.16)，h_w＝11000m，h_s是地面站高程，以米为单位，P_s是气压，以毫巴为单位；

其中

e_s＝RH·exp(-37.2465+0.21366T-0.0002568908T²)，T_s为干温，T是观测时刻的温度，采用绝对温标，单位K，RH为相对湿度。在没有测定测站气象元素的情况下，采用给定的标准气压参数计算，其中P_s取1013.25mBar，T_s取291.15K，RH取50。

根据以上模型可以求出时刻t₁的对流层相位距离延迟Δρ_trop(t₁)和时刻t₂的对流层延迟Δρ_trop(t₂)，最后能得到对流层频移的改正公式：

${\mathrm{\Δf}}_{\mathrm{trop}}=-\frac{f}{c}\frac{{\∂\mathrm{\Δ\ρ}}_{\mathrm{trop}}}{\∂t}\≈-\frac{f}{c}\frac{{\mathrm{\Δ\ρ}}_{\mathrm{trop}2}-{\mathrm{\Δ\ρ}}_{\mathrm{trop}1}}{t_2-t_1}=-\frac{{\mathrm{\Δ\ρ}}_{\mathrm{trop}2}-{\mathrm{\Δ\ρ}}_{\mathrm{trop}1}}{\left(t_2-t_1\right)c}f---\left(8\right)$

电离层频移的影响可用双频消除法消去，这是由于GPS信号的电离层频移与频率本身成反比，而其他所有频移与频率本身成正比，如下式：

$f_j^{\mathrm{obs}}=\mathrm{\α}{\tilde{f}}_j+\frac{\mathrm{\χ}}{{\tilde{f}}_j},$ j＝1，2                            (9)

$\left\{\begin{array}{c}\mathrm{\α}=\frac{k(k\stackrel{\‾}{f_1^{\mathrm{obs}}}-\stackrel{\‾}{f_2^{\mathrm{obs}}})}{{\tilde{f}}_1(k^2-1)}\\ \mathrm{\χ}=\frac{\stackrel{\‾}{f_1^{\mathrm{obs}}}-k\stackrel{\‾}{f_2^{\mathrm{obs}}}}{1-k^2}{\tilde{f}}_1\end{array}\right.---\left(10\right)$

其中，χ电离层频移项，α为其它频移项之和，k＝f₁/f₂，信号的真实发射频率分别用

和

表示。通过L1和L2载波上的D1和D2观测值列出方程组，可以解出α和χ；根据α可提取重力频移信号。

信号频率的移动与重力位差之间具有确定的关系，即重力位频移方程，由下式给出：

$W-W_S=c^2\frac{\mathrm{\Δf}}{f_0}---\left(11\right)$

式中，f₀是固有发射频率，Δf是重力频移，卫星重力位W_S可由低阶重力场模型求取。低阶重力场模型W_S是输入到程序中的已知量，它是低阶引力位V和离心力位Q之和：W_S＝V+Q。

卫星低阶引力位，即取8阶的获取方法如下：

$V(r_t,{\mathrm{\θ}}_t,{\mathrm{\λ}}_t)=\frac{\mathrm{GM}}{r_t}[1+\underset{n=1}{\overset{8}{\mathrm{\Σ}}}{\left(\frac{a}{r_t}\right)}^n\underset{m=0}{\overset{n}{\mathrm{\Σ}}}({\stackrel{\‾}{C}}_{\mathrm{nm}}\cos m{\mathrm{\λ}}_t+{\stackrel{\‾}{S}}_{\mathrm{nm}}\sin m{\mathrm{\λ}}_t){\stackrel{\‾}{P}}_{\mathrm{nm}}\left(\cos {\mathrm{\θ}}_t\right)]$

其中，GM是地心常数(已知)，a＝6378136.572m是地球赤道平均半径，P_nm(cosθ_t)是正规化缔合Legendre函数，C_nm和S_nm是正规化球谐系数，是已知值，该值由EGM2008模型给出，(r_t，θ_t，λ_t)是卫星在t时刻的球面坐标，其中r_t，θ_t，λ_t分别是卫星在t时刻的站星距离，余纬和经度。

离心力位的获取方法如下：

$Q=\frac{{\mathrm{\ω}}^2}{2}{\left(r_t\sin {\mathrm{\θ}}_t\right)}^2$

其中，ω＝7.292115×10^-5rad/s是地球自转角速率。

在本发明中，可以对整个计算过程进行多次获取，并对每次计算结果求平均值，并进行精度评估；将计算结果以″.gda″文件形式进行保存，该文件可直接被显示模块调用显示。

在本发明中，显示模块实现了以下几大项功能：

1.列表显示分析结果

日期视图中列表显示出含有数据的日期；卫星视图中列表显示出日期视图中选定日期中含有数据的卫星；历元视图中列表显示出选定日期选定卫星中含有数据的历元及每个历元的历元号，时间，高度角，重力频移，经典多普勒频移，对流层频移，电离层频移，总频移。显示模块获取用户选中的日期以及该日期内存在观测数据的卫星，从分析结果中获取对应的所有历元的数据，列表显示出来。

2.根据选中历元数据绘制图表

可绘制的图形项为重力频移，如图3，经典多普勒频移，对流层频移，电离层频移，总频移。

3.查看单站重力位/海拔高及其精度分析结果

显示模块可从分析结果中直接读取所有历元的精度评估结果；也可采用用户选择的历元所涵盖的数据进行精度评估；还可采用用户指定的高度截止角下限作为限制条件从所有历元中所筛选的数据进行精度评估。

在本发明中，所采用的装置，包括依次连接的参数配置模块1、GPS参数读取模块2、卫星定位模块3、旋转改正模块4、地面站重力位确定模块5以及地面站海拔高确定模块6，所述的地面站重力位确定模块5包括依次连接的经典多普勒频移确定单元7、对流层频移确定单元8、电离层频移确定单元9、重力频移确定单元10以及地面站重力位确定单元11。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (9)

1.一种通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，包括以下步骤： 

步骤1，由参数配置模块根据用户需求配置各分析参数，包括设置低阶重力位模型，经典多普勒频移模型，对流层模型，牛顿内插阶次以及微分时间间隔； 

步骤2，由GPS参数读取模块从GPS文件中读入GPS观测文件，包括obs文件和精密星历sp3文件，并获取总频移量； 

步骤3，由卫星定位模块根据精密星历通过牛顿内插确定卫星位置，牛顿内插选定步骤1中的内插阶次，并基于低阶重力场模型求卫星引力位； 

步骤4，由旋转改正模块实施地球旋转改正； 

步骤5，由地面站重力位确定模块根据步骤2中的总频移量确定地面站重力位；其中，步骤5中的“确定地面站重力位”的具体步骤由以下步骤完成： 

步骤5.1，由经典多普勒频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出经典多普勒频移值； 

步骤5.2，由对流层频移确定单元通过obs文件中给出的地面站位置和sp3精密星历计算出对流层频移值； 

步骤5.3，由电离层频移确定单元通过obs文件中L1和L2载波的频移值，利用双频消除法计算电离层频移； 

步骤5.4，由重力频移确定单元根据步骤2中的总频移量去掉步骤5.1、步骤5.2以及步骤5.3中的各频移值，获取到重力频移值； 

步骤5.5，由地面站重力位确定单元根据步骤5.4中获得的重力频移以及步骤3所述的通过低阶重力场模型求得的卫星引力位，基于重力频移方程，确定地面站重力位； 

步骤6，由地面站海拔高确定模块以大地水准面上的重力位为基准，基于布隆斯公式，用地面站重力位确定地面站海拔高程。 

2.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤2中，总频移量Δf可表示为 

Δf＝Δf_d+Δf_V十Δf_h+Δf_i+Δf_t+Δf_c+Δf_r； 

其中，Δf_d和Δf_h分别是经典多普勒频移和横向经典多普勒频移，Δf_V是引力频移，Δf_i和Δf_t分别是电离层和对流层影响，Δf_c是原子钟误差，Δf_r是随机误差，引力频移与横向经典多普勒频移之和是重力频移，记为Δf_g。 

3.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤5.1中，经典多普勒频移是由于卫星与地面观测站有相对速度导致地面站接收到的卫星信号的频率的变化，具体公式为： 

其中c是真空中光速，f_s和f分别为发射频率和接收频率，单位时间内测站与卫星间的距离变化为ds／dt。 

4.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤3中，牛顿法插值通过构造各阶差商产生 多项式系数，以待求点和插值点几何位置关系为自变量形成多项式，具体公式如下： 

差商公式： 

牛顿插值公式： 

f(x)＝f(x₀)+f[x₀，x₁](x-x₀) 

+f[x₀，x₁，x₂](x-x₀)(x-x₁)+...； 

+f[x₀，x₁...x_n](x-x₀)(x-x₁)...(x-x_n-1) 

由此可完成步骤3。 

5.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤4中，地球旋转改正即坐标变换，由如下方程给出 

X_t＝Xcosa+Ysina 

Y_t＝-Xsina+Ycosa 

式中，a＝2πt／(3600×24)，t为卫星信号从卫星传送到地面站所花费的时间，(X，Y)是卫星在惯性系中的位置，(X_t，Y_t)是卫星在地固系中的位置。 

6.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤5.2中，对流层频移是GPS信号在传播过程中的对流层延迟影响，采用对流层模型差分获得，可选模型有SAASTMOINEN；MODIFIELD HOPIFIELD；SIMPLIFIED HOPIFIELD；计算对流层相位距离延迟的简化的Hopifield模型公式如下： 

其中 

式中，Δρ_trop为相位延迟，E为高度角，h_d＝401386+148.72(T- 273.16)，h_w＝11000m，h_s是地面站高程，以米为单位，P_s是气压，以毫巴为单位； 

其中 

e_s＝RH·exp(- 37.2465+0.21366T- 0.0002568908T²)， 

T_s为干温，T是观测时刻的温度，采用绝对温标，单位K，RH为相对湿度，在没有测定测站气象元素的情况下，采用给定的标准气压参数计算，其中P_s取1013.25mBar，T_s取291.15K，RH取50； 

根据以上模型求出时刻t₁的对流层相位距离延迟Δρ_trop(t₁)和时刻t₂的对流层延迟Δρ_trop(t₂)，最后能得到对流层频移的改正公式： 

其中f_s是卫星发射信号的频率，c是真空中光速。 

7.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤5.3中，电离层频移的影响可用双频消除法消去，这是由于GPS信号的电离层频移与频率本身成反比，而其他所有频移与频率本身成正比，如下式： 

其中，χ电离层频移项，α为其它频移项之和，

信号的真实发射频率分别用

和

表示；通过L1和L2载波上的D1和D2观测值列出方程组，可以解出α和χ；根据α可提取重力频移信号。 

8.根据权利要求1所述的通过提取GPS信号重力频移确定海拔高程的方法，其特征在于，所述的步骤5.5中，所述的重力位频移方程，由下式给出： 

式中，f₀是固有发射频率，Δf是重力频移，c是真空中光速；卫星重力位W_S可由低阶重力场模型求取；低阶重力场模型W_S是输入到程序中的已知量，它是低阶引力位V和离心力位Q之和：W_S＝V+Q； 

卫星低阶引力位，即取8阶的获取方法如下： 

其中，GM是地心常数，a＝6378136.572m是地球赤道平均半径， 

是正规化缔合Legendre函数，和

是正规化球谐系数，是已知值，该值由EGM2008模型给出，(r_t，θ_t，λ_t)是卫星在t时刻的球面坐标， 其中r_t，θ_t，λ_t分别是卫星在t时刻的站星距离，余纬和经度； 

离心力位的获取方法如下： 

其中，ω＝7.292115×10^-5rad／s是地球自转角速率。 

9.一种利用权利要求1所述的方法的装置，其特征在于，包括依次连接的参数配置模块、GPS参数读取模块、卫星定位模块、旋转改正模块、地面站重力位确定模块以及地面站海拔高确定模块，所述的地面站重力位确定模块包括依次连接的经典多普勒频移确定单元、对流层频移确定单元、电离层频移确定单元、重力频移确定单元以及地面站重力位确定单元。 

Images