CN103323888B - Gnss大气探测数据中对流层延迟误差的消除方法 - Google Patents

Abstract

本发明公开了一种GNSS大气探测数据中对流层延迟误差的消除方法，其特征在于，融合GPT模型、UNB3m模型和经典对流层模型，所述经典对流层模型为Hopfield或Saastamoinen模型，包括如下步骤：以GPT模型计算测站的温度与气压；再以UNB3m模型计算测站的水气压；最后通过使用GMF映射函数（Global Mapping Function）计算站星视线方向的对流层延迟改正。本发明方法适用于缺少实测气象数据的情况，能够较精确的给出对流层延迟的先验值，可应用于GNSS定位、GNSS气象学等领域。

Description

GNSS大气探测数据中对流层延迟误差的消除方法

技术领域

本发明属于GNSS大气探测技术，GNSS技术为大气探测提供了一个强有力的工具，其具有无校准、全天候、精度高、垂直分辨率高、全球均匀覆盖等特点；特别涉及精确确定GNSS大气探测的对流层延迟误差的消除方法。

背景技术

利用GNSS技术获取地球大气温、压、湿等相关参数。在GNSS观测值中，电磁信号穿过大气层产生的延迟是主要误差源之一，简称大气误差。在GNSS各种误差中,包括卫星钟差、接收机钟差、大气延迟等误差；其中GNSS大气误差主要包括对流层延迟和电离层延迟两部分。其中，电离层延迟与电磁信号的传播频率相关，采用双频无电离层组合观测值的方法可以消除95％的电离层延迟误差。对流层延迟比电离层延迟情况更加复杂，观测值组合法不适用于对流层延迟改正。因此，在GNSS非差以及长基线的差分数据处理中，必须采用合适的模型改正对流层延迟，或者将其视为未知数做参数估计。

为有效消除对流层延迟对GNSS高精度数据处理的影响，国内外学者很早就开始针对对流层延迟改正展开研究，提出如Hopfield和Saastamoinen等经验模型。近年来，根据长期观测的气象数据，一些新的对流层延迟经验模型相继被提出，如欧洲的EGNOS模型和北美的UNB3m模型等。1998年，IGS中心对流层研究工作组(IGSTWG)成立，将对流层延迟纳入研究内容之中。目前，IGS分析中心提供全球300多个跟踪站的对流层延迟值，内容包括天顶对流层延迟、N方向水平梯度延迟、E方向水平梯度延迟以及各延迟分量的中误差等。

Hopfield和Saastamoinen模型需要测站位置的气象参数，但是一般测站不具备实测气象参数的条件，如果使用标准气象参数替代会降低模型的精度；EGNOS和UNB3m模型的适用范围小，全球范围内的模型精度不均匀；参数估计对流层延迟的方法受其他待估参数影响，精度与时效性难以同时保证，如CODE中心提供的高精度对流层延迟数据有五天的滞后时间。

发明内容

针对现有的快速准确得到GNSS对流层延迟技术存在的不足，本发明提出了一种融合了GPT模型、UNB3m模型和经典对流层模型Hopfield或Saastamoinen模型的新的全球对流层延迟模型，有效提高模型方法确定对流层延迟的精度且适用于全球范围。

本发明采用以下技术方案：GNSS大气探测数据中对流层延迟误差的消除方法，包括如下步骤：

步骤一、以GPT模型计算测站的温度与气压

GPT模型的输入参数为：对GNSS大气探测数据测量的测站经度、测站纬度、大地高和年积日，GPT模型的输出值为：温度、气压和高程异常，所述高程异常为大地高与正常高的差；

首先以9阶9次球谐函数求出测站位置温度与气压的年平均值a₀、年变化幅度A和高程异常N：

其中：P_nm表示勒让德多项式；A_nm和B_nm表示球谐函数的系数，通过专家天气图ECMWF中心三年的分析数据和重力模型EGM96重力模型(静态重力模型)数据以最小二乘的方法解算得到；和λ表示测站的纬度和经度；

然后在年平均值和年变化幅度的基础上，以余弦函数的形式求出具体某一年积日的温度与气压值。

$\left[\begin{array}{c}{P}_0\\ T_0\end{array}\right]=a_0+A\·\cos \left(\frac{\mathrm{doy}-28}{365.25}2\mathrm{\π}\right)---\left(2\right)$

其中，P₀和T₀表示平均海平面的温度与气压值，doy表示年积日。最后根据测站高程进行梯度改正：

P＝P₀[1-0.000226(H_大地-N)]^5.225    (3)

T＝T₀-0.0065(H_大地-N)(4)

其中：P和T表示测站位置的温度与气压值，H_大地-N表示测站高于海平面的高度，即测站的正常高；

步骤二、以UNB3m模型(全球对流层天顶延迟模型)计算测站的水气压

首先查表1得到测站附近纬度的相对湿度年平均值和年变化幅度，线性内插到测站所在纬度位置，然后以余弦函数计算某一年积日的相对湿度大小，余弦函数同公式(2)：

表1UNB3m模型相对湿度的均值和振幅经验值

最后通过相对湿度计算水气压，具体计算公式参考IERS2003规范：

$e=\frac{\mathrm{RH}}{100}\·e_s\·f_w---\left(5\right)$

e_s＝0.01exp(1.2378847×10^-5T²-1.9121316×10^-2T    (6)

+33.93711047-6.3431645×10³T^-1)

f_w＝1.00062+3.14×10^-6P+5.6×10^-7(T-273.15)²    (7)

其中：e表示水气压，RH表示相对湿度，e_s表示饱和水气压，f_w表示缩放因子；

步骤三、根据步骤一和二得到的经验气象参数，以经典对流层模型Hopfield或Saastamoinen分别计算天顶对流层延迟的干分量ZHD和湿分量ZWD改正，以Saastamoinen模型计算：

ZHD＝0.002277P    (8)

$\mathrm{ZWD}=0.002277(\frac{1255}{T}+0.05)e---\left(9\right)$

以Hopfield模型计算：

ZHD＝77.64P/T    (10)

$\mathrm{ZWD}=-12.96\frac{e}{T}+371800\frac{e}{T^2}---\left(11\right)$

步骤四、使用GMF映射函数(Global Mapping Function全球映射函数模型)计算站星视线方向的对流层延迟改正：

$m{\left(\mathrm{\ϵ}\right)}_h=\frac{1+\frac{a_h}{1+\frac{b_h}{1+c_h}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_h}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_h}{\sin \left(\mathrm{\ϵ}\right)+c_h}}}+\frac{1}{\sin \left(\mathrm{\ϵ}\right)}-\frac{1+\frac{a_{\mathrm{ht}}}{1+\frac{b_{\mathrm{ht}}}{1+c_{\mathrm{ht}}}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+c_{\mathrm{ht}}}}}---\left(12\right)$

$m{\left(\mathrm{\ϵ}\right)}_w=\frac{1+\frac{a_w}{1+\frac{b_w}{1+c_w}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_w}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_w}{\sin \left(\mathrm{\ϵ}\right)+c_w}}}---\left(13\right)$

ΔT＝m(ε)_h·ZHD+m(ε)_w·ZWD    (14)

其中：ε表示卫星高度角，a,b,c表示模型中的球谐函数系数，H表示大地高，以米为单位，下标h对应干分量系数，下标ht对应干分量高程修正系数，下标w对应湿分量系数，ΔT表示站星视线方向的对流层延迟总量。

本发明的有益效果：(1)GPT/UNB3m融合经典模型确定全球对流层延迟的方法有效提高了对流层延迟模型的精度，扩大了适用范围，过程中不再使用标准气象参数，充分体现了对流层延迟的时空变化特性。(2)GPT/UNB3m融合经典模型相比于北美UNB3m模型和欧洲EGNOS模型，与CODE中心计算的对流层延迟的变化曲线更吻合，特别体现在南半球中高纬度地区，提升了对流层延迟模型的普适性。本发明融合了GPT模型、UNB3m模型和经典对流层模型(Hopfield或Saastamoinen模型)，其中GPT/UNB3m模型是用来提供气象参数用于Hopfield模型，得到了很好的结果。本发明方法适用于缺少实测气象数据的情况，能够较精确的给出对流层延迟的先验值，可应用于GNSS定位、GNSS气象学等领域。

附图说明

图1测试所用IGS连续运行跟踪站分布图(深色点表示测站)；

图2为计算测站位置的天顶对流层延迟的流程；

图3BJFS(北京)站不同对流层延迟模型精度比较；

图4SANT(智利)站不同对流层延迟模型精度比较；

图5LHAZ(拉萨)站不同对流层延迟模型精度比较；

图6OHI3(南极)站不同对流层延迟模型精度比较。

具体实施方式

下面结合附图对本发明做更进一步地解释。

本发明提出一种新的确定全球对流层延迟的方法，包括如下步骤：

(1)以GPT模型计算测站的温度与气压：

GPT模型的输入参数为：测站经度、测站纬度、大地高和年积日，GPT模型的输出值为：温度、气压和高程异常，所高程异常为大地高与正常高的差；

首先以9阶9次球谐函数求出测站位置温度与气压的年平均值a₀、年变化幅度A和高程异常N：

其中：P_nm表示勒让德多项式；A_nm和B_nm表示球谐函数的系数，通过ECMWF中心三年的分析数据和EGM96重力模型数据以最小二乘的方法解算得到；和λ表示测站的纬度和经度；

然后在年平均值和年变化幅度的基础上，以余弦函数的形式求出具体某一年积日的温度与气压值。

$\left[\begin{array}{c}{P}_0\\ T_0\end{array}\right]=a_0+A\·\cos \left(\frac{\mathrm{doy}-28}{365.25}2\mathrm{\π}\right)---\left(2\right)$

其中，P₀和T₀表示平均海平面的温度与气压值，doy表示年积日。最后根据测站高程进行梯度改正：

P＝P₀[1-0.000226(H_大地-N)]^5.225    (3)

T＝T₀-0.0065(H_大地-N)    (4)

其中：P和T表示测站位置的温度与气压值，H_大地-N表示测站高于海平面的高度，即测站的正常高；

(2)以UNB3m模型计算测站的水气压：

首先查表1得到测站附近纬度的相对湿度年平均值和年变化幅度，线性内插到测站所在纬度位置，然后以余弦函数计算某一年积日的相对湿度大小，余弦函数同公式(2)：

表1UNB3m模型相对湿度的均值和振幅经验值

最后通过相对湿度计算水气压，具体计算公式参考IERS2003：

$e=\frac{\mathrm{RH}}{100}\·e_s\·f_w---\left(5\right)$

e_s＝0.01exp(1.2378847×10^-5T²-1.9121316×10^-2T    (6)

+33.93711047-6.3431645×10³T^-1)

f_w＝1.00062+3.14×10^-6P+5.6×10^-7(T-273.15)²    (7)

其中：e表示水蒸气局部气压，简称水气压，RH表示相对湿度，e_s表示饱和水气压，f_w表示缩放因子；

(3)根据前两步得到的经验气象参数，以经典对流层计算模型Hopfield或Saastamoinen(简称H或S模型)分别计算天顶对流层延迟的干分量和湿分量改正：

以Saastamoinen模型计算：

ZHD＝0.002277P    (8)

$\mathrm{ZWD}=0.002277(\frac{1255}{T}+0.05)e---\left(9\right)$

以Hopfield模型计算：

ZHD＝77.64P/T(10)

$\mathrm{ZWD}=-12.96\frac{e}{T}+371800\frac{e}{T^2}---\left(11\right)$

(4)使用GMF映射函数(Global Mapping Function)计算站星视线方向的对流层延迟改正STD：

$m{\left(\mathrm{\ϵ}\right)}_h=\frac{1+\frac{a_h}{1+\frac{b_h}{1+c_h}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_h}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_h}{\sin \left(\mathrm{\ϵ}\right)+c_h}}}+\frac{1}{\sin \left(\mathrm{\ϵ}\right)}-\frac{1+\frac{a_{\mathrm{ht}}}{1+\frac{b_{\mathrm{ht}}}{1+c_{\mathrm{ht}}}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+c_{\mathrm{ht}}}}}---\left(12\right)$

$m{\left(\mathrm{\ϵ}\right)}_w=\frac{1+\frac{a_w}{1+\frac{b_w}{1+c_w}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_w}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_w}{\sin \left(\mathrm{\ϵ}\right)+c_w}}}---\left(13\right)$

ΔT＝m(ε)_h·ZHD+m(ε)_w·ZWD    (14)

其中：ε表示卫星高度角，a,b,c表示模型中的球谐函数系数，下标h对应干分量系数，下标ht对应干分量高程修正系数，下标w对应湿分量系数，ΔT表示站星视线方向的对流层延迟总量。

为了真实有效地反映各类对流层延迟模型的精度，选取分布全球的一共244个IGS连续运行跟踪站进行测试，测站分布如图1所示。总共设计了四套获取天顶对流层延迟(ZTD)的方案。

方案1：从网上下载IGS的CODE分析中心最终对流层延迟分析产品，数据以SINEX文件格式存储，提供2小时间隔的测站天顶对流层延迟信息；

方案2：以图2中的流程计算测站位置的天顶对流层延迟，其中经典对流层模型选择Saastamoinen模型；

方案3：以图2中的流程计算测站位置的天顶对流层延迟，其中经典对流层模型选择Hopfield模型；

方案4：以UNB3m模型计算测站位置的天顶对流层延迟。

欧洲定轨中心CODE是权威的GNSS数据处理中心，负责各种基于GPS/GLONASS组合的数据解算。天顶对流层延迟数据是它的产品之一，以CODwwwwd.TRO.Z命名并存放于FTP服务器免费提供下载。CODE使用参数估计的方式计算得到对流层天顶延迟，相比于模型计算值更为精确，因此将CODE每天12组对流层延迟值(以2小时间隔)的平均值作为每日对流层延迟的参考真值，分别对方案2、3、4计算得到的每日对流层延迟值的误差做统计分析。分析数据的时间跨度为1408周到1668周，总共5年。具体结果如表2所示。

表2不同对流层延迟模型误差的均方根比较及误差分布

其中RMS代表中误差。

由表2看出，方案2和3整体略好于方案4。三种方案的误差均方根非常接近，误差分布也差不多，80％以上的历元对流层模型误差小于5cm。因此，这三种方案都可以作为经验模型计算天顶对流层延迟值。但是具体到某些站点时，三种方案在精度和适用范围上有差别。取部分不同纬度、不同大地高、具有代表性的站点，对不同方案得到的天顶对流层延迟值作分析，具体如图3～6和表3所示。图3～6中曲线中断部分表示CODE提供的数据缺失，相应历元的模型计算值也没有画出。

表3不同对流层延迟算法精度统计(单位：cm)

通过比较分析，实验结果如下：

(1)如图3和表3所示，在BJFS(北京)站，明显GPT/UNB3m融合经典模型要优于UNB3m模型。UNB3m的对流层延迟值振幅略小，不能精确的表达出北京地区对流层延迟值的年度变化。

(2)如图4和表3所示，在SANT(智利)站，对流层延迟的模型值都偏高，但UNB3m的误差更大，相比而言GPT/UNB3m融合经典模型能较好的反映对流层延迟信息。

(3)如图5和表3所示，在LHAZ(拉萨)站，由于海拔很高的缘故，对流层延迟值误差普遍较大，几种方案之间无明显差别。

(4)如图6和表3所示，在OHI3(南极)站，对流层延迟无明显周期性变化，GPT/UNB3m融合经典模型能较好的体现出这一特点，而UNB3m模型不能。这也验证了GPT/UNB3m融合经典模型能较好的反映南半球中高纬度地区温度气压真实值这一特点。

本发明的实施例不构成对本发明范围的限定，而是以审定的权利要求及其等同物确定。

Claims (2)

1.一种GNSS大气探测数据中对流层延迟误差的消除方法，其特征在于，融合GPT模型、UNB3m模型和经典对流层模型，所述经典对流层模型为Hopfield或Saastamoinen模型，包括如下步骤：

步骤一、以GPT模型计算测站的温度与气压，

GPT模型的输入参数为：测站经度、测站纬度、大地高和年积日，GPT模型的输出值为：温度、气压和高程异常，所述高程异常为大地高与正常高的差；

首先以9阶9次球谐函数求出测站位置温度与气压的年平均值a₀、年变化幅度A和高程异常N：

其中：P_nm表示勒让德多项式；A_nm和B_nm表示球谐函数的系数，通过ECMWF中心三年的分析数据和EGM96重力模型数据以最小二乘的方法解算得到；和λ表示测站的纬度和经度；m和n代表0至9的整数；

然后在年平均值和年变化幅度的基础上，以余弦函数的形式求出具体某一年积日的温度与气压值：

$\left[\begin{array}{c}{P}_0\\ T_0\end{array}\right]=a_0+A\·\cos \left(\frac{\mathrm{doy}-28}{365.25}2\mathrm{\π}\right)---\left(2\right)$

其中：P₀和T₀表示平均海平面的温度与气压值，doy表示年积日；

最后根据测站高程进行梯度改正：

P＝P₀[1-0.000226(H_大地-N)]^5.225    (3)

T＝T₀-0.0065(H_大地-N)    (4)

其中：P和T表示测站位置的温度与气压值，H_大地-N表示测站高于海平面的高度，即测站的正常高；

步骤二、以UNB3m模型计算测站的水气压

首先查表1得到测站附近纬度的相对湿度年平均值和年变化幅度，线性内插得到测站所在纬度位置，然后以余弦函数计算某一年积日的相对湿度大小，余弦函数同公式(2)：

表1 UNB3m模型相对湿度的均值和振幅经验值

最后通过相对湿度计算水气压，具体计算公式参考IERS2003：

$e=\frac{\mathrm{RH}}{100}\·e_s\·f_w---\left(5\right)$

e_s＝0.01exp(1.2378847×10^-5T²-1.9121316×10^-2T

                                            (6)

              +33.93711047-6.3431645×10³T^-1

f_w＝1.00062+3.14×10^-6P+5.6×10^-7(T-273.15)²    (7)

其中：e表示水气压，RH表示相对湿度，e_s表示饱和水气压，f_w表示缩放因子；

步骤三、根据步骤一和二得到的经验气象参数，以经典对流层模型Hopfield或Saastamoinen分别计算天顶对流层延迟的干分量ZHD和湿分量ZWD改正，以Saastamoinen模型计算：

ZHD＝0.002277P    (8)

$\mathrm{ZWD}=0.002277(\frac{1255}{T}+0.05)e---\left(9\right)$

以Hopfield模型计算：

ZHD＝77.64P/T    (10)

$\mathrm{ZWD}=-12.96\frac{e}{T}+371800\frac{e}{T^2}---\left(11\right)$

步骤四、使用GMF映射函数(Global Mapping Function)计算站星视线方向的对流层延迟改正：

$m{\left(\mathrm{\ϵ}\right)}_h=\frac{1+\frac{a_h}{1+\frac{b_h}{1+c_h}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_h}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_h}{\sin \left(\mathrm{\ϵ}\right)+c_h}}}+\frac{1}{\sin \left(\mathrm{\ϵ}\right)}-\frac{1+\frac{a_{\mathrm{ht}}}{1+\frac{b_{\mathrm{ht}}}{1+c_{\mathrm{ht}}}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_{\mathrm{ht}}}{\sin \left(\mathrm{\ϵ}\right)+c_{\mathrm{ht}}}}}---\left(12\right)$

$m{\left(\mathrm{\ϵ}\right)}_w=\frac{1+\frac{a_w}{1+\frac{b_w}{1+c_w}}}{\sin \left(\mathrm{\ϵ}\right)+\frac{a_w}{\sin \left(\mathrm{\ϵ}\right)+\frac{b_w}{\sin \left(\mathrm{\ϵ}\right)+c_w}}}---\left(13\right)$

ΔT＝m(ε)_h·ZHD+m(ε)_w·ZWD    (14)

其中：ε表示卫星高度角a,b,c表示模型中的球谐函数系数，H表示大地高，以米为单位，下标h对应干分量系数，下标ht对应干分量高程修正系数，下标w对应湿分量系数，ΔT表示站星视线方向的对流层延迟总量。

2.根据权利要求1所述的GNSS大气探测数据中对流层延迟误差的消除方法，其特征在于，利用GPT模型和UNB3m模型来提供测站气象参数，融合经典对流层模型即Hopfield或Saastamoinen模型进行大气延迟误差改正。