CN109917424A - 多因子约束下的nwp反演对流层延迟的残差改正方法 - Google Patents

Abstract

本发明公开了一种多因子约束下的NWP反演对流层延迟的残差改正方法，首先获取作业区域里所有连续运行参考站的连续一年的ZTD和该区域NWP提供的连续一年的再分析资料，利用积分法反演该区域上述连续运行参考站的ZTD，然后基于最小绝对残差法的多项式拟合法拟合NWP反演ZTD的残差，获取当前时刻的NWP预报资料和观测站的初始坐标，利用上述残差模型计算当前时刻改正后的NWP_ZTD的残差，利用当前时刻的NWP预报资料反演当前时刻的NWP_ZTD的初值，将上述NWP_ZTD的初值与改正后的NWP_ZTD的残差相加，即得到当前时刻精确的NWP_ZTD。本方法首次根据温度、湿度、位置等多因素与对流层残差的时间序列间的关系来建模，改正NWP反演的对流层延迟的残差提高对流层延迟的估计精度至亚厘米级。

Description

多因子约束下的NWP反演对流层延迟的残差改正方法

技术领域

本发明涉及一种对流层延迟的残差改正方法，尤其是多因子约束下的NWP反演对流层延迟的残差改正方法，利用多种气象因素时间序列的变化规律建模，改正数值天气预报模型(NWP)反演的对流层延迟的残差，从而提高精密单点定位(PPP)和长距离基线实时精密动态定位(RTK)的收敛速度和定位精度，属于全球卫星导航与定位技术领域。

背景技术

对流层延迟是全球卫星导航***(Global Navigation Satellite System,简称GNSS)导航和定位的主要误差源之一，其对电磁波信号产生的影响是非色散的折射。由于对流层折射的影响，在天顶方向可使电磁波的传播路径差s达到2.3m，当卫星高度角为10°时可达20m，严重制约着GNSS精密单点定位和中长距离基线差分定位的模糊度收敛速度和定位精度。

在GNSS数据处理的现有技术中，总结起来有三种方法可用来消弱对流层延迟：1.通过经验模型估计天顶对流层延迟(ZTD)，但经验模型估计的ZTD的精度差，最大误差可达20cm；2.在数学模型中将ZTD作为参数与其他未知参数(如站点坐标和模糊度)共同求解。然而，对流层延迟与位置强相关，方程需要较强的卫星几何约束和较长时间的观测才能得到稳定解；3.利用数值天气预报模型的再分析或预报资料反演ZTD。

现如今NWP预报资料反演ZTD的残差大概在±6cm上下浮动，平均残差和均方根值(RMS)分别在1cm和3cm左右，不能直接用于GNSS精密定位中。在GNSS精密定位解算中，通常将此ZTD作为初值，将ZTD的残差作为未知数，与模糊度，位置等参数共同求解。此时，NWP反演ZTD的精度将直接影响模糊度的收敛速度。要实现实时PPP或长距离RTK，需要提高NWP反演ZTD的精度，即改正NWP反演的ZTD的残差，针对这一问题，目前还未有切实可行的方法。

发明内容

本发明的目的在于，为了克服现有技术中存在的不足，提供了一种多因子约束下的NWP反演对流层延迟的残差改正方法，用于解决现有的方法得到的对流层延迟精度差的技术问题。

为了解决背景技术所存在的问题，本发明的多因子约束下的NWP反演对流层延迟的残差改正方法，包括以下步骤如下：

第一步，获取作业区域里所有连续运行参考站的连续一年的ZTD，记做IGS_ZTD，并剔除粗差；该数据可在各省市连续运行参考站网中心网站获得。

第二步，获取该区域NWP提供的连续一年的再分析资料，包括气压、温度、相对湿度及位势高度，该数据可由欧洲中尺度天气预报中心(ECMWF)再分析资料中ERA-Interim(Jan1979-present)产品的分层气象数据获得。

第三步，利用积分法反演该区域上述连续运行参考站的ZTD，记做NWP_ZTD。采用分层积分法反演ZTD的计算公式为：

上式中，ZTD_grid为格网点在测站所在高度上的格网点的ZTD值，H_IGS为参考站的高程，H_top为NWP资料顶层高度，N是大气折射率(无量纲)，N_i表示第i个积分区域范围内的大气折射率，△H_i是在第i个积分区域范围内的高度。

其中大气折射率N的近似表达公式如下：

上式中，k₁＝77.604K/mbar，k₂＝64.79K/mbar，k₃＝377600.0K²/mbar，P表示压强，e表示水汽压，单位均是mbar，N_try表示干性气体的折射率，N_wet表示湿性气体的折射率。

第四步，在相对湿度、温度和位置的约束下，基于最小绝对残差法的多项式拟合法拟合NWP反演ZTD的残差：

f(RH,T₀)＝p₀₀+p₁₀*RH+p₀₁*T₀+p₂₀*RH²+p₁₁*RH*T₀+p₀₂*T₀ ²+p₃₀*RH³

+p₂₁*RH²*T₀+p₁₂*RH*T₀ ²+p₀₃*T₀ ³+p₃₁*RH³*T₀

+p₂₂*RH²*T₀ ²+p₁₃*RH*T₀ ³+p₀₄*T₀ ⁴

式中：RH代表相对湿度，单位为％；T₀代表温度，单位为℃；P₀₀、P₁₀、P₀₁、P₂₀、P₁₁、P₀₂、P₃₀、P₂₁、P₁₂、P₀₃、P₃₁、P₂₂、P₁₃、P₀₄是多项式的系数；f(RH,T₀)代表拟合的ZTD残差。

第五步，获取当前时刻的NWP预报资料和观测站的初始坐标，利用上述残差模型计算当前时刻改正后的NWP_ZTD的残差。

第六步，利用当前时刻的NWP预报资料，采用积分法反演当前时刻的NWP_ZTD的初值。

第七步，将上述NWP_ZTD的初值与改正后的NWP_ZTD的残差相加，即可得到当前时刻精确的NWP_ZTD。

第八步，将此精确的NWP_ZTD直接用于GNSS实时精密定位中精密单点定位和长距离基线实时精密动态定位。

与现有技术相比，本发明的有益效果是：

第一，首次根据温度、湿度、位置等多因素与对流层残差的时间序列间的关系来建模，改正NWP反演的对流层延迟的残差提高对流层延迟的估计精度至亚厘米级。

第二，将本方法应用在PPP和RTK中，可显著缩短模糊度收敛时间，提高定位精度。

第三，由于NWP模型可提供覆盖全球的气象数据，所以本方法可估计全球NWP模型覆盖下的任意一点的ZTD，适用范围很广。

附图说明

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

图2是SCIP站采用本发明方法前后对流层延迟对比图。

图3是PIN1站采用本发明方法前后对流层延迟对比图。

图4是SCIP站采用本发明方法的RTK解算情况。

图5是PIN1站采用本发明方法的RTK解算情况。

说明：图2中(a)和(b)为SCIP站2016年采用本发明方法改正前、后的ZTD残差对比；图3中(a)和(b)为PIN1站2016年采用本发明方法改正前、后的ZTD残差对比；图4是以SCIP站为流动站、基线长为148km的基线；图5是以PIN1站为流动站、基线长为216km的基线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为证明基于多因子约束下的NWP反演对流层延迟的残差改正方法对提高NWP反演对流层延迟精度的有效性，下面选取IGS MGEX网SCIP站(32.91°N,118.49°W)与PIN1站(33.61°N,116.46°W)2015年和2016年两整年的对流层数据以及SCIP站、PIN1站、ROCK站(34.24°N,118.68°W)三个站点2016年1月13日数据组成的两条长基线进行测试，一条基线长度为148km(在SCIP站与ROCK站之间，其中ROCK站为基站，SCIP站为流动站)，另一条基线长度为216km(在PIN1站与ROCK站之间，其中ROCK站为基站，PIN1站为流动站)。同时，选取三个站点相应区域相应时间的欧洲中尺度天气预报中心(ECMWF)再分析资料中ERA-Interim(Jan 1979-present)产品的分层气象数据，其平面分辨率为0.125°×0.125°，垂直分辨率为37层(顶层高度约为47km)，时间分辨率为6小时。

实施例1

第一步，在IGS MGEX网中获取SCIP站和ROCK站2015年连续一年的ZTD，记做IGS_ZTD，并剔除粗差。

第二步，在欧洲中尺度天气预报中心(ECMWF)再分析资料中ERA-Interim(Jan1979-present)产品中，获取平面分辨率为0.125°×0.125°、垂直分辨率为37层、时间分辨率为6小时的SCIP站和ROCK站区域2015年连续一年的分层气象数据，气象数据包括气压、温度、相对湿度及位势高度。

第三步，利用积分法反演SCIP站区域气象数据得到SCIP站2015年的ZTD，记做NWP_ZTD。采用分层积分法反演ZTD的计算公式为：

式中ZTD_grid为格网点在测站所在高度上的格网点的ZTD值，H-IGS为参考站的高程，H-top为NWP资料顶层高度，N是大气折射率(无量纲)，N_i表示第i个积分区域范围内的大气折射率，△H_i是在第i个积分区域范围内的高度。

式中，大气折射率N的近似表达公式如下：

式中，k₁＝77.604K/mbar，k₂＝64.79K/mbar，k₃＝377600.0K²/mbar，P表示压强，e表示水汽压，单位均是mbar，N_try表示干性气体的折射率，N_wet表示湿性气体的折射率。

第四步，在相对湿度、温度和位置的约束下，基于最小绝对残差法的多项式拟合法拟合NWP反演的2015年ZTD的残差，得到残差模型如下：

f(RH,T₀)＝p₀₀+p₁₀*RH+p₀₁*T₀+p₂₀*RH²+p₁₁*RH*T₀+p₀₂*T₀ ²+p₃₀*RH³

+p₂₁*RH²*T₀+p₁₂*RH*T₀ ²+p₀₃*T₀ ³+p₃₁*RH³*T₀

+p₂₂*RH²*T₀ ²+p₁₃*RH*T₀ ³+p₀₄*T₀ ⁴

式中：RH代表相对湿度，单位为％；T₀代表温度，单位为℃；P₀₀、P₁₀、P₀₁、P₂₀、P₁₁、P₀₂、P₃₀、P₂₁、P₁₂、P₀₃、P₃₁、P₂₂、P₁₃、P₀₄为多项式的系数；f(RH,T₀)代表拟合的ZTD残差。

第五步，获取SCIP站2016年的NWP预报资料，利用上述残差模型计算当前时刻改正后的NWP_ZTD的残差。

第六步，利用SCIP站2016年的NWP预报资料，采用积分法反演当前时刻的NWP_ZTD的初值。

第七步，将上述NWP_ZTD的初值与改正后的NWP_ZTD的残差相加，即可得到当前时刻精确的NWP_ZTD。SCIP站2016年的ZTD采用上述方法改正前、后残差对比如图2所示，图(a)为改正前每天的ZTD，主要集中在20mm处上下浮动，图(b)为改正后每天的ZTD，主要集中在0mm处上下浮动，从图中可看出，改正后的ZTD精度明显提高。

在IGS MGEX网中下载SCIP站和ROCK站2016年1月13日对应的RTK解算所需要的卫星观测数据文件、精密星历文件、广播星历文件、精密钟差文件，利用上述精确估计的NWP_ZTD对SCIP站与ROCK站组成的基线分别进行标准RTK解算，标准NWP约束下的RTK解算以及改正NWP约束下的RTK解算如图4所示，在三种方法下，模糊度解算时间依次缩短，解算成功率依次提高，明显体现了本发明方法的优越性。

实施例2

第一步，在IGS MGEX网中获取PIN1站和ROCK站2015年连续一年的ZTD，记做IGS_ZTD，并剔除粗差。

第二步，在欧洲中尺度天气预报中心(ECMWF)再分析资料中ERA-Interim(Jan1979-present)产品中，获取平面分辨率为0.125°×0.125°、垂直分辨率为37层、时间分辨率为6小时的PIN1站和ROCK站区域2015年连续一年的分层气象数据，气象数据包括气压、温度、相对湿度及位势高度。

第三步，利用积分法反演PIN1站区域气象数据得到SCIP站2015年的ZTD，记做NWP_ZTD。采用分层积分法反演ZTD的计算公式为：

上式中：ZTD_grid为格网点在测站所在高度上的格网点的ZTD值，H-IGS为参考站的高程，H-top为NWP资料顶层高度，N是大气折射率(无量纲)，N_i表示第i个积分区域范围内的大气折射率，△H_i是在第i个积分区域范围内的高度。

其中，大气折射率N的近似表达公式如下：

式中，k₁＝77.604K/mbar，k₂＝64.79K/mbar，k₃＝377600.0K²/mbar，P表示压强，e表示水汽压，单位均是mbar，N_try表示干性气体的折射率，N_wet表示湿性气体的折射率。

第四步，在相对湿度、温度和位置的约束下，基于最小绝对残差法的多项式拟合法拟合NWP反演的2015年ZTD的残差，得到残差模型如下：

f(RH,T₀)＝p₀₀+p₁₀*RH+p₀₁*T₀+p₂₀*RH²+p₁₁*RH*T₀+p₀₂*T₀ ²+p₃₀*RH³

+p₂₁*RH²*T₀+p₁₂*RH*T₀ ²+p₀₃*T₀ ³+p₃₁*RH³*T₀

+p₂₂*RH²*T₀ ²+p₁₃*RH*T₀ ³+p₀₄*T₀ ⁴

上式中：RH代表相对湿度，单位为％；T₀代表温度，单位为℃；P₀₀、P₁₀、P₀₁、P₂₀、P₁₁、P₀₂、P₃₀、P₂₁、P₁₂、P₀₃、P₃₁、P₂₂、P₁₃、P₀₄为多项式的系数；f(RH,T₀)代表拟合的ZTD残差。

第五步，获取PIN1站2016年的NWP预报资料，利用上述残差模型计算当前时刻改正后的NWP_ZTD的残差。

第六步，利用PIN1站2016年的NWP预报资料，采用积分法反演当前时刻的NWP_ZTD的初值。

第七步，将上述NWP_ZTD的初值与改正后的NWP_ZTD的残差相加，即可得到当前时刻精确的NWP_ZTD。PIN1站2016年的ZTD采用上述方法改正前、后残差对比如图2所示，图(a)为改正前每天的ZTD，主要集中在20mm处上下浮动，图(b)为改正后每天的ZTD，主要集中在0mm处上下浮动，从图中可看出，改正后的ZTD精度明显提高。

在IGS MGEX网中下载PIN1站和ROCK站2016年1月13日对应的RTK解算所需要的卫星观测数据文件、精密星历文件、广播星历文件、精密钟差文件，利用上述精确估计的NWP_ZTD对PIN1站与ROCK站组成的基线分别进行标准RTK解算，标准NWP约束下的RTK解算以及改正NWP约束下的RTK解算如图5所示，在三种方法下，模糊度解算时间依次缩短，解算成功率依次提高，明显体现了本发明方法的优越性。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (1)

1.多因子约束下的NWP反演对流层延迟的残差改正方法，其特征在于，包含以下具体步骤：

第一步，获取作业区域里所有连续运行参考站的连续一年的ZTD，记做IGS_ZTD，并剔除粗差，该数据可在各省市连续运行参考站网中心网站获得；

第二步，获取该区域NWP提供的连续一年的再分析资料，包括气压、温度、相对湿度及位势高度，该数据可由欧洲中尺度天气预报中心(ECMWF)再分析资料中ERA-Interim产品的分层气象数据获得；

第三步，利用积分法反演该区域上述连续运行参考站的ZTD，记做NWP_ZTD，采用分层积分法反演ZTD的计算公式为：

上式中，ZTD_grid为格网点在测站所在高度上的格网点的ZTD值，H_IGS为参考站的高程，H_top为NWP资料顶层高度，N是大气折射率(无量纲)，N_i表示第i个积分区域范围内的大气折射率，△H_i是在第i个积分区域范围内的高度；

其中大气折射率N的近似表达公式如下：

上式中，k₁＝77.604K/mbar，k₂＝64.79K/mbar，k₃＝377600.0K²/mbar，P表示压强，e表示水汽压，单位均是mbar，N_try表示干性气体的折射率，N_wet表示湿性气体的折射率；

第四步，在相对湿度、温度和位置的约束下，基于最小绝对残差法的多项式拟合法拟合NWP反演ZTD的残差：

f(RH,T₀)＝p₀₀+p₁₀*RH+p₀₁*T₀+p₂₀*RH²+p₁₁*RH*T₀+p₀₂*T₀ ²+p₃₀*RH³

+p₂₁*RH²*T₀+p₁₂*RH*T₀ ²+p₀₃*T₀ ³+p₃₁*RH³*T₀

+p₂₂*RH²*T₀ ²+p₁₃*RH*T₀ ³+p₀₄*T₀ ⁴

式中：RH代表相对湿度，单位为％；T₀代表温度，单位为℃；P₀₀、P₁₀、P₀₁、P₂₀、P₁₁、P₀₂、P₃₀、P₂₁、P₁₂、P₀₃、P₃₁、P₂₂、P₁₃、P₀₄是多项式的系数；f(RH,T₀)代表拟合的ZTD残差；

第五步，获取当前时刻的NWP预报资料和观测站的初始坐标，利用上述残差模型计算当前时刻改正后的NWP_ZTD的残差；

第六步，利用当前时刻的NWP预报资料，采用积分法反演当前时刻的NWP_ZTD的初值；

第七步，将上述NWP_ZTD的初值与改正后的NWP_ZTD的残差相加，即可得到当前时刻精确的NWP_ZTD。