CN102147475B - 利用gps信号同时确定三维几何位置和重力位的方法及其装置 - Google Patents

Abstract

本发明涉及一种同时确定三维几何位置和重力位的方法及其装置，尤其是涉及一种利用GPS信号同时确定三维几何位置和重力位的方法及其装置。本发明在已有的GPS定位方法的基础上，加进利用GPS信号确定重力位(包括海拔高程)的方法，后者是一种新发明的方法，其技术方案是：建立GPS信号与重力位之间的观测方程，利用多次观测数据，基于最小二乘原理求解重力位。因此，本发明具有如下优点：不需要繁重的水准测量以及重力测量，也不需要事先确定大地水准面再利用GPS定位方法确定海拔高，从而节省了大量的人力、财力，而且在也不受山区及跨海洋测量限制。

Description

利用GPS信号同时确定三维几何位置和重力位的方法及其装置

技术领域

本发明涉及一种同时确定三维几何位置和重力位的方法及其装置，尤其是涉及一种利用GPS信号同时确定三维几何位置和重力位的方法及其装置。

背景技术

利用GPS信号确定三维几何位置已是非常成熟的技术。为了确定重力位或海拔高，按传统方法，需要繁重的水准测量以及重力测量，或者需要事先确定大地水准面再利用GPS定位方法确定海拔高，这不仅需要大量的人力、财力，而且在山区会受到很大限制。由于GPS信号包含三维几何位置信息，从而可以确定三维几何位置，这是目前非常成熟的方法。又由于GPS信号包含电磁波信号(简称信号)频移信息，因而也可利用GPS信号确定重力位(包括海拔高程)。这是本发明方法的基础。

发明内容

本发明主要是解决现有技术所存在的需要繁重的水准测量以及重力测量，或者需要事先确定大地水准面再利用GPS定位方法确定海拔高，这不仅需要大量的人力、财力，而且在山区会受到很大限制等的技术问题；提供了一种不需要繁重的水准测量以及重力测量，也不需要事先确定大地水准面再利用GPS定位方法确定海拔高，从而节省了大量的人力、财力，而且在也不受山区及跨海洋测量限制的利用GPS信号同时确定三维几何位置和重力位的方法及其装置。

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

一种利用GPS信号同时确定三维几何位置和重力位的方法，其特征在于，包括以下步骤：

步骤1，利用GPS接收机采集观测至少4颗卫星的、持续至少20分钟的相位以及信号频移观测值，观测值的采样率为1秒，每个观测值对应有观测时刻；

步骤2，由三维几何位置确定模块利用步骤1中的相位观测值以及精密星历确定地面测站三维几何位置；

步骤3，由卫星引力位确定模块利用EGM2008重力场模型以及精密星历确定GPS卫星在对应于观测值时刻的引力位；

步骤4，由测站重力位确定模块将步骤2中确定的三维几何位置以及步骤3中卫星处引力位信息引入GPS信号相位-频移观测方程之中，将测站的重力位以及时钟误差作为待定参数引入GPS信号相位-频移观测方程之中，同时加进一阶多普勒频移、电离层频移、对流层频移等模型改正，建立GPS信号相位-频移观测方程，并基于最小二乘原理获取测站的重力位；

步骤5，由海拔高程确定模块基于求解的测站重力位，以大地水准面上的重力位为基准，利用Bruns公式确定测站的海拔高程；

步骤6，由信息输出模块输出三维几何位置和绝对重力位及海拔高程。

在上述的利用GPS信号同时确定三维几何位置和重力位的方法，所述的步骤4中，

GPS信号的相位-频移观测方程

${\mathrm{\Φ}}_A=\frac{f}{c}{\mathrm{\ρ}}_A=\frac{f_0}{c}{\mathrm{\ρ}}_A+\frac{\mathrm{\Δf}}{c}{\mathrm{\ρ}}_A$

其中，Φ_A是测站A处的GPS信号的相位观测值，f₀和f分别是信号发射和接收频率，ρ_A是星站间距离，c是真空中光速。

在上述的利用GPS信号同时确定三维几何位置和重力位的方法，所述的步骤4中，Δf＝f-f₀是信号频移量，由如下方程确定

Δf＝Δf_d+Δf_i+Δf_t+Δf_c+ΔΔf_g                    (1)

其中，Δf_d、Δf_i和Δf_t分别是经典多普勒频移、电离层频移和对流层频移，均由模型值给出；时钟误差Δf_c是待定参数；Δf_g是重力频移，由下式给出

$\mathrm{\Δ}{f}_g=f\frac{W_A-W_S^i}{c^2}$

其中测站处的重力位W_A是待定参数，W_S ⁱ是沿轨飞行卫星在t_i时刻的重力位，可表示成

$W_S^i=W[x_S\left(t_i\right),y_S\left(t_i\right),z_S\left(t_i\right)]$

其中，x_S(t_i)，y_S(t_i)，z_S(t_i)是卫星在t_i时刻的在地心惯性系中的坐标。

在上述的利用GPS信号同时确定三维几何位置和重力位的方法，包括依次连接的GPS接收机、三维几何位置确定模块、卫星引力位确定模块、测站重力位确定模块、海拔高程确定模块以及信息输出模块。

因此，本发明具有如下优点：1.在现有的利用GPS技术确定三维几何位置的基础上，同时实现利用GPS信号确定重力位(包括海拔高程)，实现GPS两位测量***2.不需要繁重的水准测量以及重力测量，也不需要事先确定大地水准面再利用GPS定位方法确定海拔高，从而节省了大量的人力、财力，而且在也不受山区及跨海洋测量限制。

附图说明

附图1是本发明的一种原理图；

附图2是本发明的一种工作流程图；

具体实施方式

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

实施例：

利用GPS信号同时确定三维几何位置和重力位的方法，包括以下步骤：

步骤1，利用GPS接收机采集观测至少4颗卫星的、持续至少20分钟的相位以及信号频移观测值，观测值的采样率为1秒，每个观测值对应有观测时刻；

步骤2，由三维几何位置确定模块利用步骤1中的相位观测值以及精密星历确定地面测站三维几何位置；

步骤3，由卫星引力位确定模块利用EGM2008重力场模型以及精密星历确定GPS卫星在对应于观测值时刻的引力位；

步骤4，由测站重力位确定模块将步骤2中确定的三维几何位置以及步骤3中卫星处引力位信息引入GPS信号相位-频移观测方程之中，将测站的重力位以及时钟误差作为待定参数引入GPS信号相位-频移观测方程之中，同时加进一阶多普勒频移、电离层频移、对流层频移等模型改正，建立GPS信号相位-频移观测方程，并基于最小二乘原理获取测站的重力位；

本步骤的一阶多普勒频移、电离层频移、对流层频移等模型改正，具体操作步骤为：

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

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

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

首先介绍上述步骤4.1至4.3所需用到的模型和相关公式：

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

$\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{\ρ}}_{\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{\ρ}}_{\mathrm{trop}1}}{t_2-t_1}=-\frac{{\mathrm{\Δ\ρ}}_{\mathrm{trop}2}-{\mathrm{\Δ\ρ}}_{\mathrm{trop}1}}{(t_2-t_1)c}f---\left(8\right)$

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

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

$\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。

本发明的具体相关计算步骤及原理如下：

利用GPS信号确定三维几何位置的模型省略，可在一般性教科书中找到。确定重力位(包括海拔高程)的模型如下：

I)首先建立GPS信号的相位观测方程

${\mathrm{\Φ}}_A=\frac{f}{c}{\mathrm{\ρ}}_A=\frac{f_0}{c}{\mathrm{\ρ}}_A+\frac{\mathrm{\Δf}}{c}{\mathrm{\ρ}}_A---\left(1\right)$

其中，Φ_A是测站A处的GPS信号的相位观测值，f₀和f分别是信号发射和接收频率，ρ_A是星站间距离，c是真空中光速，Δf＝f-f₀是信号频移量，由如下方程确定

Δf＝Δf_d+Δf_i+Δf_t+Δf_c+Δf_g                  (2)

其中，Δf_d、Δf_i和Δf_t分别是经典多普勒频移、电离层频移和对流层频移，均由模型值给出；时钟误差Δf_c是待定参数；Δf_g是重力频移，由下式给出

$\mathrm{\Δ}{f}_g=f\frac{W_A-W_S^i}{c^2}---\left(3\right)$

其中测站处的重力位W_A是待定参数，W_S ⁱ是沿轨飞行卫星在t_i时刻的重力位，可表示成

$W_S^i=W[x_S\left(t_i\right),y_S\left(t_i\right),z_S\left(t_i\right)]---\left(4\right)$

其中，x_S(t_i)，y_S(t_i)，z_S(t_i)是卫星在t_i时刻的在地心惯性系中的坐标(由精密星历给出)。由方程(1)-(4)，可写出如下观测方程

${\mathrm{\Φ}}_A=\frac{f_0}{c}{\mathrm{\ρ}}_A+\frac{{\mathrm{\Δf}}_d+{\mathrm{\Δf}}_i+{\mathrm{\Δf}}_t}{c}{\mathrm{\ρ}}_A+\frac{1}{c}{\mathrm{\ρ}}_A{\mathrm{\Δf}}_c+\frac{W_A-W[x_S\left(t_i\right),y_S\left(t_i\right),z_S\left(t_i\right)]}{c^3}{\mathrm{\ρ}}_{A}f---\left(5\right)$

II)由于卫星轨道(精密星历)已知，因而可根据最新发布的重力场模型EGM2008确定卫星在t_i时刻的重力位：

W[x_S(t_i)，y_S(t_i)，z_S(t_i)]＝V[x_S(t_i)，y_S(t_i)，z_S(t_i)]+Q[x_S(t_i)，y_S(t_i)，z_S(t_i)]      (6)

其中Q[x_S(t_i)，y_S(t_i)，z_S(t_i)]是离心力位，由下式给出

$Q[x_S\left(t_i\right),y_S\left(t_i\right),z_S\left(t_i\right)]=\frac{{\mathrm{\ω}}^2}{2}[x_S^2\left(t_i\right)+y_S^2\left(t_i\right)]---\left(7\right)$

其中ω＝7.292115×10^-5rad/s是地球自转角速率，V[x_S(t_i)，y_S(t_i)，z_S(t_i)]是地球引力位，可根据后面陈述的第III)步骤确定。

III)基于EGM2008模型，可在地固参考系中的球面坐标(r，θ，λ)下确定地球引力位：

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

其中GM是地心常数(已知)，a＝6378136.572m是地球赤道平均半径，P_nm(cosθ)是正规化缔合Legendre函数，C_nm和S_nm是正规化球谐系数，由EGM2008模型给出。

IV)将地心惯性系下的坐标(x_S，y_S，z_S)转化为地固坐标系下的坐标(r，θ，λ)即可由方程(8)确定在任意时刻t_i在轨卫星所在处的引力位(这里省略时刻标记t_i)，进而根据方程(6)和(7)可确定在任意时刻t_i在轨卫星所在处的重力位。坐标转换方程如下：

x_S＝rcosλsinθ

y_S＝rsinλsinθ

z_S＝rcosθ                             (9)

其逆变换为

$r=\sqrt{{x_S}^2+{y_S}^2+{z_S}^2}$

$\mathrm{\θ}={\cos }^{-1}\left(\frac{z_S}{r}\right)$

$\mathrm{\λ}={\tan }^{-1}\left(\frac{y_S}{x_S}\right)---\left(10\right)$

V)对于每个观测值均可建立一个观测方程(5)，于是得到大量的多余观测方程。基于最小二乘原理求解，即可得到地面上测站A处的重力位W_A。

VI)利用大地水准面上的重力位常数W₀＝62636855.580，即可采用如下的Bruns公式

$H_A=\frac{W_0-W_A}{\stackrel{\‾}{\mathrm{\γ}}}---\left(11\right)$

确定测站A点的海拔高程H_A，其中γ是参考椭球面上的平均正常重力。

根据上述模型，可利用GPS信号确定地面测站的重力位及海拔高程。再结合已有的确定三维几何位置的方法，即可实现GPS两位测量，即利用GPS信号同时确定三维几何位置和重力位(包括海拔高程)。

在本发明中，所采用的装置，包括依次连接的GPS接收机1、三维几何位置确定模块2、卫星引力位确定模块3、测站重力位确定模块4、海拔高程确定模块5以及信息输出模块6。

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

Claims (3)

1.一种利用GPS信号同时确定三维几何位置和重力位的方法，其特征在于，包括以下步骤：

步骤1，利用GPS接收机采集观测至少4颗卫星的、持续至少20分钟的相位以及信号频移观测值，观测值的采样率为1秒，每个观测值对应有观测时刻；

步骤2，由三维几何位置确定模块利用步骤1中的相位观测值以及精密星历确定地面测站三维几何位置；

步骤3，由卫星引力位确定模块利用EGM2008重力场模型以及精密星历确定GPS卫星在对应于观测时刻的引力位；

步骤4，由测站重力位确定模块将步骤2中确定的三维几何位置以及步骤3中卫星在对应于观测时刻的引力位信息引入GPS信号相位-频移观测方程之中，将测站的重力位以及时钟误差作为待定参数引入GPS信号相位-频移观测方程之中，同时加进一阶多普勒频移、电离层频移、对流层频移模型改正，建立GPS信号相位-频移观测方程，并基于最小二乘原理获取测站的重力位；

步骤5，由海拔高程确定模块基于求解的测站重力位，以大地水准面上的重力位为基准，利用Bruns公式确定测站的海拔高程；

步骤6，由信息输出模块输出测站所在点的三维几何位置、重力位及海拔高程。

2.根据权利要求1所述的利用GPS信号同时确定三维几何位置和重力位的方法，其特征在于，所述的步骤4中，

GPS信号的相位-频移观测方程

${\mathrm{\Φ}}_A=\frac{f}{c}{\mathrm{\ρ}}_A=\frac{f_0}{c}{\mathrm{\ρ}}_A+\frac{\mathrm{\Δf}}{c}{\mathrm{\ρ}}_A$

其中，Φ_A是测站A处的GPS信号的相位观测值，f₀和f分别是信号发射和接收频率，ρ_A是星站间距离，c是真空中光速，Δf＝f-f₀是信号频移量；

所述的步骤4中，Δf＝f-f₀是信号频移量，f₀和f分别是信号发射和接收频率，由如下方程确定

Δf＝Δf_d+Δf_i+Δf_t+Δf_c+Δf_g

其中，Δf_d、Δf_i和Δf_t分别是经典多普勒频移、电离层频移和对流层频移；时钟误差Δf_c是待定参数；Δf_g是重力频移，由下式给出

$\mathrm{\Δ}{f}_g=f\frac{W_A-W_S^i}{c^2}$

其中c是真空中光速，而测站处的重力位W_A是待定参数，

是沿轨飞行卫星在t_i时刻的重力位，可表示成

$W_S^i=W[x_S\left(t_i\right),y_S\left(t_i\right),z_S\left(t_i\right)]$

其中，x_S(t_i)，y_S(t_i)，z_S(t_i)是卫星在t_i时刻的在地心惯性系中的坐标，其中，“t_i时刻”的含义是：在第i个观测时刻点，其中，i＝1，2，…，N。

3.一种利用权利要求1所述的方法的装置，其特征在于，包括依次连接的GPS接收机(1)、三维几何位置确定模块(2)、卫星引力位确定模块(3)、测站重力位确定模块(4)、海拔高程确定模块(5)以及信息输出模块(6)。

Images