CN115755115A - PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology - Google Patents

PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology Download PDF

Info

Publication number
CN115755115A
CN115755115A CN202211570241.0A CN202211570241A CN115755115A CN 115755115 A CN115755115 A CN 115755115A CN 202211570241 A CN202211570241 A CN 202211570241A CN 115755115 A CN115755115 A CN 115755115A
Authority
CN
China
Prior art keywords
gnss
std
ppp
delay
tropospheric
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
CN202211570241.0A
Other languages
Chinese (zh)
Inventor
赵庆志
王鹏程
蒋朵朵
李祖锋
袁荣才
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian University of Science and Technology
Original Assignee
Xian University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian University of Science and Technology filed Critical Xian University of Science and Technology
Priority to CN202211570241.0A priority Critical patent/CN115755115A/en
Publication of CN115755115A publication Critical patent/CN115755115A/en
Withdrawn legal-status Critical Current

Links

Images

Landscapes

  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了基于GNSS对流层层析技术的PPP改善方法,涉及GNSS精密单点定位技术领域,针对目前PPP技术中观测方程包含的对流层延迟作为未知参数难以精确估计,观测方程中的未知参数较多,PPP技术的定位精度和收敛速度不够精准的问题,现提出如下方案,包括以下步骤:S1:基于PPP技术的斜路径对流层延迟STD估算;S2:层析区域网格内折射率反演;S3:利用GNSS对流层层析技术的斜路径对流层延迟STD反演;S4:GNSS对流层层析技术改善PPP定位方法。本发明引入GNSS对流层层析技术,对卫星信号路径上的STD值进行反演,并将反演的STD改正到PPP技术的观测方程中,以减少观测方程中的未知参数,进一步提升PPP技术的定位精度和收敛速度。

Figure 202211570241

The invention discloses a PPP improvement method based on GNSS tropospheric tomography technology, and relates to the technical field of GNSS precise single point positioning. It is difficult to accurately estimate the tropospheric delay contained in the observation equation in the current PPP technology as an unknown parameter, and there are many unknown parameters in the observation equation. , the positioning accuracy and convergence speed of PPP technology are not accurate enough, the following scheme is proposed, including the following steps: S1: Estimation of oblique path tropospheric delay STD based on PPP technology; S2: Refractive index inversion in the tomographic area grid; S3 : Oblique path tropospheric delay STD inversion using GNSS tropospheric tomography; S4: GNSS tropospheric tomography improves PPP positioning method. The present invention introduces GNSS tropospheric tomography technology, inverts the STD value on the satellite signal path, and corrects the inverted STD into the observation equation of the PPP technology, so as to reduce the unknown parameters in the observation equation and further improve the PPP technology. Positioning accuracy and convergence speed.

Figure 202211570241

Description

基于GNSS对流层层析技术的PPP改善方法PPP improvement method based on GNSS tropospheric tomography technology

技术领域technical field

本发明涉及GNSS精密单点定位技术领域,尤其涉及基于GNSS对流层层析技术的PPP改善方法。The invention relates to the technical field of GNSS precise single point positioning, in particular to a PPP improvement method based on GNSS tropospheric tomography technology.

背景技术Background technique

全球卫星定位导航系统最核心和基础的功能与服务是授时和定位,其中定位可分为相对定位和绝对定位。相对定位的代表性技术为实时载波相位差分定位技术,而绝对定位的代表技术是精密单点定位(Precise Point Positioning,PPP)技术。PPP技术是指利用卫星实时或事后的精密轨道和钟差产品,对观测方程中的各项误差进行精确改正,同时采取合适的参数估计策略,在仅需一台GNSS接收机的情况下,便可实现全球范围内精密绝对定位。但GNSS信号也会受到各类误差的干扰,当信号穿过50km以下的中性大气层时,会被其中的中性大气干扰产生折射,这种折射影响被称为对流层延迟。根据卫星高度角不同,对流层延迟会对定位结果产生数米甚至数十米的误差影响。因此,对流层延迟是PPP定位中最主要的误差源之一。The core and basic functions and services of the global satellite positioning and navigation system are timing and positioning, of which positioning can be divided into relative positioning and absolute positioning. A representative technology of relative positioning is a real-time carrier phase differential positioning technology, and a representative technology of absolute positioning is a precise point positioning (Precise Point Positioning, PPP) technology. PPP technology refers to the use of satellite real-time or post-event precision orbit and clock error products to accurately correct various errors in the observation equation, and at the same time adopt an appropriate parameter estimation strategy. In the case of only one GNSS receiver, it is convenient. Accurate absolute positioning can be achieved worldwide. However, GNSS signals will also be interfered by various errors. When the signal passes through the neutral atmosphere below 50km, it will be refracted by the neutral atmospheric interference. This refraction effect is called tropospheric delay. Depending on the altitude angle of the satellite, the tropospheric delay will have an error of several meters or even tens of meters in the positioning result. Therefore, tropospheric delay is one of the most important sources of error in PPP positioning.

传统PPP定位中,斜路径对流层延迟(Slant Tropospheric Delay,STD)可表示为天顶方向对流层延迟(Zenith Tropospheric Delay,ZTD)与映射函数(MappingFunctions,MF)的乘积加上大气水平梯度改正项与其映射函数乘积。ZTD可分为天顶干延迟(Zenith Hydrostatic Delay,ZHD)和天顶湿延迟(Zenith Wet Delay,ZWD)。在处理中,由于ZHD较稳定且易于模型化,故常采用经验模型进行改正,如Saastamoinen模型、Hopfiled模型等;而ZWD不确定性较强,难以建立精确模型改正,因此利用参数估计方法确定其大小;梯度改正项为水平梯度映射函数乘以水平梯度改正北方向与东方向的分量。但经验模型和未知参数估计往往存在一定的偏差,为了减少这种估计误差,本发明提出一种利用GNSS层析技术改善PPP定位精度及收敛速度的方法。In traditional PPP positioning, the Slant Tropospheric Delay (STD) can be expressed as the product of the Zenith Tropospheric Delay (ZTD) and the mapping function (MappingFunctions, MF) plus the atmospheric horizontal gradient correction term and its mapping function product. ZTD can be divided into Zenith Dry Delay (Zenith Hydrostatic Delay, ZHD) and Zenith Wet Delay (Zenith Wet Delay, ZWD). In processing, because ZHD is relatively stable and easy to model, empirical models are often used for correction, such as Saastamoinen model, Hopfield model, etc.; while ZWD has strong uncertainty, it is difficult to establish an accurate model correction, so the parameter estimation method is used to determine its size ; The gradient correction item is the horizontal gradient mapping function multiplied by the horizontal gradient to correct the north and east components. However, there is often a certain deviation between the empirical model and the unknown parameter estimation. In order to reduce this estimation error, the present invention proposes a method for improving the PPP positioning accuracy and convergence speed by using GNSS tomographic technology.

地基GNSS对流层层析技术是指利用大量卫星信号射线上的信息重构区域三维对流层延迟结构的方法。该技术首先将研究区域离散成若干个网格,并假设卫星信号在每个网格内的折射率在一定时间内为一常数且分布均匀;每条信号路径上的观测值包含了该路径上的对流层延迟信息,因此,将每条射线在不同网格内折射率与该卫星信号在该网格内截距相乘并进一步积分得到该路径上的斜路径总延迟,多个测站对多个卫星形成的所有倾斜路径观测方程即构成了GNSS对流层层析观测方程,再结合建立的层析约束方程,得到GNSS对流层层析模型,然后通过迭代算法或非迭代算法解算GNSS对流层层析模型,从而获得各网格的折射率信息。Ground-based GNSS tropospheric tomography refers to the method of reconstructing the regional three-dimensional tropospheric delay structure by using the information on a large number of satellite signal rays. This technique first discretizes the research area into several grids, and assumes that the refractive index of the satellite signal in each grid is a constant and uniformly distributed within a certain period of time; the observations on each signal path include the The tropospheric delay information, therefore, the refraction index of each ray in different grids is multiplied by the intercept of the satellite signal in the grid and further integrated to obtain the total delay of the oblique path on the path. All oblique path observation equations formed by satellites constitute the GNSS tropospheric tomographic observation equation, combined with the established tomographic constraint equations, the GNSS tropospheric tomographic model is obtained, and then the GNSS tropospheric tomographic model is solved by iterative or non-iterative algorithms , so as to obtain the refractive index information of each grid.

但是,目前PPP技术中观测方程包含的对流层延迟作为未知参数难以精确估计,观测方程中的未知参数较多,PPP技术的定位精度和收敛速度不够精准。However, the tropospheric delay contained in the observation equation in the current PPP technology is difficult to accurately estimate as an unknown parameter. There are many unknown parameters in the observation equation, and the positioning accuracy and convergence speed of the PPP technology are not accurate enough.

发明内容Contents of the invention

本发明的目的是为了解决目前PPP技术中观测方程包含的对流层延迟作为未知参数难以精确估计,观测方程中的未知参数较多,PPP技术收敛速度较慢的缺点,提出GNSS对流层层析技术改善精密单点定位精度及加快收敛速度方法。The purpose of the present invention is to solve the disadvantages that the tropospheric delay included in the observation equation in the current PPP technology is difficult to accurately estimate as an unknown parameter, there are many unknown parameters in the observation equation, and the PPP technology has a slow convergence speed. Single-point positioning accuracy and a method to speed up the convergence speed.

为了实现上述目的,本发明采用了如下技术方案:In order to achieve the above object, the present invention adopts the following technical solutions:

基于GNSS对流层层析技术的PPP改善方法,包括以下步骤:The PPP improvement method based on GNSS tropospheric tomography technology comprises the following steps:

S1:基于PPP技术的斜路径对流层延迟STD估算:S1: STD estimation of oblique path tropospheric delay based on PPP technology:

利用PPP技术估计STD时,需综合考虑观测方程的各项误差,确保得到的STD具有较可靠的初始精度,STD计算具体过程如下:When using PPP technology to estimate STD, it is necessary to comprehensively consider the errors of the observation equation to ensure that the obtained STD has a relatively reliable initial accuracy. The specific process of STD calculation is as follows:

S1.1:PPP观测方程建立:S1.1: PPP observation equation establishment:

GNSS信号在从卫星到接收机的传播中会受到各项误差的干扰,这些误差会不同程度的影响定位精度,由于GNSS伪距观测值其本身精度为米级,因此,选择改正精度更高的载波相位观测方程,在充分顾及各项误差的情况下,GNSSPPP技术中载波相位观测方程具体表达如下:The GNSS signal will be interfered by various errors in the propagation from the satellite to the receiver, and these errors will affect the positioning accuracy to varying degrees. Since the GNSS pseudo-range observation value itself has a meter-level accuracy, it is better to choose the one with higher correction accuracy. The carrier phase observation equation, in the case of taking full account of various errors, the carrier phase observation equation in GNSSPPP technology is specifically expressed as follows:

Figure BDA0003987631920000031
Figure BDA0003987631920000031

式中,

Figure BDA0003987631920000032
表示载波相位观测值,i为卫星信号频率,r和s分别代表接收机和卫星系统,
Figure BDA0003987631920000033
表示星地几何距离,tr,i为接收机钟差,
Figure BDA0003987631920000034
表示卫星钟差,
Figure BDA0003987631920000035
为电离层延迟,
Figure BDA0003987631920000036
为对流层延迟,Br,i
Figure BDA0003987631920000037
则分别为接收机端和卫星端的相位硬件延迟,λi为信号频率的波长,
Figure BDA0003987631920000038
表示整周模糊度,
Figure BDA0003987631920000039
为载波相位的残差项;In the formula,
Figure BDA0003987631920000032
Indicates the carrier phase observation value, i is the satellite signal frequency, r and s represent the receiver and satellite system respectively,
Figure BDA0003987631920000033
Indicates the geometric distance between the satellite and the earth, t r,i is the clock error of the receiver,
Figure BDA0003987631920000034
Indicates the satellite clock error,
Figure BDA0003987631920000035
is the ionospheric delay,
Figure BDA0003987631920000036
is the tropospheric delay, B r,i and
Figure BDA0003987631920000037
are the phase hardware delays of the receiver and the satellite respectively, λi is the wavelength of the signal frequency,
Figure BDA0003987631920000038
represents the ambiguity of the whole week,
Figure BDA0003987631920000039
is the residual term of the carrier phase;

S1.2:天顶对流层延迟ZTD估计:S1.2: Zenith Tropospheric Delay ZTD Estimation:

在传统PPP解算中,将基准站坐标作为常数处理,在已知其精确坐标的情况下,通过可联合地表实测气压的Saastamoinen等经验模型计算先验ZHD,将已知坐标与计算出的先验ZHD一同带入观测方程作为起算数据,随后对观测方程线性化,进一步联合多测站多卫星观测信号构建观测方程矩阵,并利用整体最小二乘或卡尔曼滤波方法进行解算,即可得到估计的天顶对流层总延迟ZTD和梯度项;In the traditional PPP solution, the coordinates of the reference station are treated as constants. When the precise coordinates are known, the prior ZHD is calculated through the empirical model such as Saastamoinen, which can be combined with the surface measured air pressure. The known coordinates and the calculated prior The experimental ZHD is brought into the observation equation as the initial calculation data, and then the observation equation is linearized, and the observation equation matrix is further combined with multi-station and multi-satellite observation signals, and the overall least squares or Kalman filter method is used to solve it, and then we can get Estimated zenith tropospheric total delay ZTD and gradient terms;

S1.3:斜路径总延迟STD计算:S1.3: Calculate the total delay STD of the oblique path:

斜路径对流层延迟可表示为天顶方向对流层延迟与映射函数的乘积加上大气水平梯度改正项与其映射函数乘积;The oblique path tropospheric delay can be expressed as the product of the tropospheric delay in the zenith direction and the mapping function plus the product of the atmospheric horizontal gradient correction term and its mapping function;

天顶对流层总延迟ZTD主要包括天顶干延迟ZHD和天顶湿延迟ZWD两部分,通过地表气象参数精确计算ZHD,随后由PPP估计出的ZTD减去ZHD可得天顶湿延迟ZWD,表示如下:The total zenith tropospheric delay ZTD mainly includes two parts: the zenith dry delay ZHD and the zenith wet delay ZWD. The ZHD is accurately calculated through the surface meteorological parameters, and then the zenith wet delay ZWD can be obtained by subtracting the ZHD from the ZTD estimated by the PPP, which is expressed as follows :

ZWD=ZTD-ZHD(2)ZWD=ZTD-ZHD(2)

因此,斜路径对流层延迟STD可表示如下:Therefore, the oblique path tropospheric delay STD can be expressed as follows:

Figure BDA0003987631920000041
Figure BDA0003987631920000041

式中,MFh表示ZHD对应的斜路径干映射函数,MFw表示ZWD对应的斜路径湿映射函数,MFg表示水平梯度映射函数,Gn和Ge分别表示水平梯度改正的北方向分量与东方向分量,

Figure BDA0003987631920000042
为测站到卫星的方位角,εt表示对流层延迟残差;In the formula, MF h represents the oblique path dry mapping function corresponding to ZHD, MF w represents the oblique path wet mapping function corresponding to ZWD, MF g represents the horizontal gradient mapping function, G n and Ge represent the north direction component and the horizontal gradient correction respectively east component,
Figure BDA0003987631920000042
is the azimuth angle from the station to the satellite, ε t is the tropospheric delay residual;

S2:层析区域网格内折射率反演:S2: Refractive index inversion in the grid of the tomographic region:

将利用PPP技术估算的斜路径对流层延迟STD用于GNSS对流层层析模型的构建中,具体步骤如下:The oblique path tropospheric delay STD estimated by PPP technology is used in the construction of the GNSS tropospheric tomography model. The specific steps are as follows:

S2.1:层析模型观测方程构建:S2.1: Construction of tomographic model observation equation:

通过将层析区域在三维方向上划分为多个离散的网格,并对穿过各网格的卫星信号射线上的折射率进行积分,可构建STD与折射率的积分方程:By dividing the tomographic area into multiple discrete grids in the three-dimensional direction, and integrating the refractive index on the satellite signal rays passing through each grid, the integral equation of STD and refractive index can be constructed:

STD=10-6·∫sNds (4)STD=10 -6s Nds (4)

式中,N为大气折射率,包括干折射率和湿折射率,s表示由卫星传播到接收机的信号的路径长度。In the formula, N is the atmospheric refractivity, including dry refractivity and wet refractivity, and s represents the path length of the signal transmitted from the satellite to the receiver.

将上式离散化,则GNSS信号斜路径上对流层总延迟STD表示卫星信号穿过每个网格的截距与该网格内大气折射率乘积之和:Discretizing the above formula, the total tropospheric delay STD on the oblique path of the GNSS signal represents the sum of the intercept of the satellite signal passing through each grid and the product of the atmospheric refractive index in the grid:

STD=∑ijk(aijk·xijk) (5)STD=∑ ijk (a ijk x ijk ) (5)

式中,xijk表示(i,j,k)网格内的待估折射率,aijk表示射线在(i,j,k)网格内的截距,STD表示由PPP技术估算的GNSS卫星信号斜路径上的对流层总延迟估计值。In the formula, x ijk represents the refraction index to be estimated in the (i, j, k) grid, a ijk represents the intercept of the ray in the (i, j, k) grid, and STD represents the GNSS satellite estimated by the PPP technique Estimated total tropospheric delay on signal slant paths.

将研究区域内所有从顶部穿出的信号射线上的STD均用式(5)表示,则构成如下层析观测方程:The STD on all the signal rays passing through the top in the study area are expressed by formula (5), then the following tomographic observation equation is formed:

y=A·x (6)y=A x (6)

式中,y为从研究区域顶部穿出的信号射线上的STD组成的列向量,A为观测方程的系数矩阵,x为未知折射率参数组成的列向量。In the formula, y is a column vector composed of STD on the signal ray passing through the top of the study area, A is the coefficient matrix of the observation equation, and x is a column vector composed of unknown refractive index parameters.

S2.2:层析约束方程构建:S2.2: Construction of tomographic constraint equation:

由于层析区域上空GNSS卫星分布不均以及站点数目不足,使得层析区域很多网格没有射线穿过,导致观测方程的系数矩阵病态,在对层析方程求解时会出现不适定问题,因此,需构建一定的约束方程建立网格在水平和垂直方向上的函数关系,可根据高斯加权函数方法或水平平滑约束的方法构建水平约束方程,依据折射率随高度增加呈指数递减的特性,建立垂直网格内的垂直函数关系。此外,可根据层析区域内无线电探空资料、数值预报再分析资料等建立层析模型的先验信息;Due to the uneven distribution of GNSS satellites over the tomography area and the insufficient number of stations, many grids in the tomography area do not have rays passing through, resulting in ill-conditioned coefficient matrix of the observation equation, and ill-posed problems will occur when solving the tomography equation. Therefore, It is necessary to construct a certain constraint equation to establish the functional relationship between the grid in the horizontal and vertical directions. The horizontal constraint equation can be constructed according to the Gaussian weighted function method or the horizontal smooth constraint method. According to the characteristic that the refractive index decreases exponentially with the increase of height, the vertical Vertical functional relationships within the grid. In addition, the prior information of the tomographic model can be established according to the radiosonde data in the tomographic area, the numerical forecast reanalysis data, etc.;

S2.3:GNSS对流层层析模型构建:S2.3: GNSS tropospheric tomography model construction:

根据上述构建的GNSS层析观测方程、水平约束方程、垂直约束方程以及先验约束方程,建立GNSS对流层层析模型:According to the GNSS tomographic observation equation, horizontal constraint equation, vertical constraint equation and prior constraint equation constructed above, a GNSS tropospheric tomographic model is established:

Figure BDA0003987631920000061
Figure BDA0003987631920000061

式中,H、V和I分别指水平、垂直和先验约束方程的系数矩阵,C表示先验约束信息或通过探空数据等方法统计得到的层析区域内的折射率值。In the formula, H, V, and I refer to the coefficient matrix of the horizontal, vertical, and prior constraint equations, respectively, and C represents the prior constraint information or the refractive index value in the tomographic region obtained by statistical methods such as sounding data.

S2.4:GNSS对流层层析模型解算:S2.4: GNSS tropospheric tomography model solution:

本发明通过奇异值分解法(SingularValueDecomposition,SVD)对GNSS对流层层析模型(7)进行解算,将层析模型系数矩阵

Figure BDA0003987631920000062
分解为:The present invention solves the GNSS troposphere tomography model (7) by singular value decomposition (SingularValueDecomposition, SVD), and the tomography model coefficient matrix
Figure BDA0003987631920000062
Decomposed into:

B=UΛVT (8)B=UΛV T (8)

式中,B∈Rm×n,U∈Rm×m,V∈Rn×n

Figure BDA0003987631920000063
∑=diag(σ12,…,σr),σ1≥σ2≥…≥σr,σi(i=1,2,…,r)为矩阵ATA的特征值的平方根,r为矩阵B的秩(r≤min(m,n)),U是由矩阵AAT的特征向量组成的正交矩阵,V是由矩阵ATA的特征向量组成的正交矩阵。如果矩阵B的广义逆定义为:In the formula, B∈R m×n , U∈R m×m , V∈R n×n ,
Figure BDA0003987631920000063
∑=diag(σ 12 ,…,σ r ), σ 1 ≥σ 2 ≥…≥σ r , σ i (i=1,2,…,r) is the square root of the eigenvalues of the matrix A T A , r is the rank of matrix B (r≤min(m,n)), U is an orthogonal matrix composed of eigenvectors of matrix A T , and V is an orthogonal matrix composed of eigenvectors of matrix A T A. If the generalized inverse of matrix B is defined as:

B-1=VΛ-1UΤ(9)B -1 = VΛ -1 U Τ (9)

那么,线性方程组Bx=L的解,即层析区域内的折射率可表示为:Then, the solution of the linear equation system Bx=L, that is, the refractive index in the tomographic region can be expressed as:

x=B-1L=VΛ-1UΤL(10)x=B -1 L=VΛ -1 U Τ L(10)

S3:利用GNSS对流层层析技术的斜路径对流层延迟STD反演:S3: Oblique-path Tropospheric Delay STD Retrieval Using GNSS Tropospheric Tomography:

GNSS信号斜路径上的对流层总延迟STD可表示为射线穿过层析网格内的折射率与信号在该网格截距的乘积之和,因此,基于GNSS对流层层析技术获取的层析区域每个网格内的折射率,根据卫星信号在每个网格内的截距即可计算卫星信号在该网格内的STD,最后通过累计方式得到卫星信号路径上的STD,具体公式如下:The total tropospheric delay STD on the oblique path of the GNSS signal can be expressed as the sum of the product of the refractive index of the ray passing through the tomographic grid and the signal intercept in the grid. Therefore, the tomographic area obtained based on GNSS tropospheric tomography The refractive index in each grid can calculate the STD of the satellite signal in the grid according to the intercept of the satellite signal in each grid, and finally obtain the STD on the satellite signal path by accumulating, the specific formula is as follows:

Figure BDA0003987631920000071
Figure BDA0003987631920000071

式中,xijk表示由层析技术反演得到的(i,j,k)网格内的折射率,aijk表示射线在(i,j,k)网格内的截距,

Figure BDA0003987631920000072
表示利用GNSS对流层层析技术反演的斜路径对流层总延迟恢复值;In the formula, x ijk represents the refractive index in the (i, j, k) grid obtained by tomographic inversion, a ijk represents the intercept of the ray in the (i, j, k) grid,
Figure BDA0003987631920000072
Indicates the tropospheric total delay recovery value of the oblique path retrieved by GNSS tropospheric tomography technology;

通过上述方法,即可计算出层析区域内不同高度角和方位角的GNSS信号斜路径上的对流层总延迟STD;Through the above method, the total tropospheric delay STD on the oblique path of the GNSS signal at different elevation angles and azimuth angles in the tomographic area can be calculated;

S4:GNSS对流层层析技术改善PPP定位方法:S4: GNSS tropospheric tomography technology improves PPP positioning method:

将GNSS对流层层析技术反演的卫星信号传播路径上的STD加入到PPP观测方程中,对传统PPP观测方程进行改进,则公式(1)可表达成如下形式:The STD on the satellite signal propagation path retrieved by GNSS tropospheric tomography technology is added to the PPP observation equation, and the traditional PPP observation equation is improved. The formula (1) can be expressed as follows:

Figure BDA0003987631920000073
Figure BDA0003987631920000073

式中,

Figure BDA0003987631920000074
表示利用GNSS层析技术反演得到的卫星s到接收机r的斜路径对流层总延迟。In the formula,
Figure BDA0003987631920000074
Indicates the total tropospheric delay of the oblique path from satellite s to receiver r retrieved by GNSS tomographic technology.

因此,顾及GNSS层析技术改善后的PPP观测方程可表达如下:Therefore, taking into account the improved PPP observation equation of GNSS tomographic technology, it can be expressed as follows:

Figure BDA0003987631920000081
Figure BDA0003987631920000081

式中,

Figure BDA0003987631920000082
表示消除了STD项的载波相位观测值,
Figure BDA0003987631920000083
为经过改正后的载波相位残差,其余各项不变;In the formula,
Figure BDA0003987631920000082
Indicates the carrier phase observations with the STD term removed,
Figure BDA0003987631920000083
is the corrected carrier phase residual, and the rest remain unchanged;

最后,进一步对上式线性化,并联合多测站多卫星载波相位观测值构建PPP技术的函数模型和误差模型,利用整体最小二乘或卡尔曼滤波方法对PPP函数模型进行解算,得到高精度的PPP定位结果,并加快PPP收敛速度。Finally, the above formula is further linearized, and the function model and error model of the PPP technology are constructed by combining the carrier phase observations of multiple stations and multiple satellites. Accurate PPP positioning results and speed up PPP convergence.

本发明中,所述的基于GNSS对流层层析技术的PPP改善方法,利用GNSS对流层层析技术精确估计PPP观测方程中所需的斜路径对流层延迟STD,首先需要利用PPP技术估计出不同测站对应多个卫星信号路径上的STD;然后,利用PPP技术获取的STD构建GNSS对流层层析观测方程,并联合层析约束方程构建GNSS对流层层析模型,通过迭代或非迭代方法对层析模型解算得到层析区域每个网格内的折射率;再次,根据GNSS卫星信号斜路径上的对流层延迟可表示为射线穿过层析区域内网格的折射率与相应截距的乘积之和,计算每条卫星信号路径上的STD;最后,将估算的STD反带入到PPP的观测方程中,减少观测方程中的未知参数,达到提高精度和收敛速度的目的;In the present invention, the PPP improvement method based on GNSS tropospheric tomography technology uses GNSS tropospheric tomography technology to accurately estimate the oblique path tropospheric delay STD required in the PPP observation equation. First, it is necessary to use PPP technology to estimate the corresponding STD on multiple satellite signal paths; then, use the STD obtained by PPP technology to construct the GNSS tropospheric tomographic observation equation, and combine the tomographic constraint equations to construct the GNSS tropospheric tomographic model, and solve the tomographic model by iterative or non-iterative methods Obtain the refractive index in each grid in the tomographic area; again, according to the tropospheric delay on the oblique path of the GNSS satellite signal, it can be expressed as the sum of the products of the refractive index of the ray passing through the grid in the tomographic area and the corresponding intercept, and calculate STD on each satellite signal path; finally, bring the estimated STD back into the PPP observation equation, reduce the unknown parameters in the observation equation, and achieve the purpose of improving accuracy and convergence speed;

本发明设计合理,针对PPP技术中观测方程包含的对流层延迟作为未知参数难以精确估计的现状,引入GNSS对流层层析技术,对卫星信号路径上的STD值进行反演,并将反演的STD改正到PPP技术的观测方程中,以减少观测方程中的未知参数,通过引入GNSS对流层层析技术计算的STD去弥补PPP技术参数估计的缺陷,进一步提升PPP技术的定位精度和收敛速度。The present invention is reasonable in design, aiming at the fact that the tropospheric delay contained in the observation equation in the PPP technology is difficult to accurately estimate as an unknown parameter, the GNSS tropospheric tomography technology is introduced to invert the STD value on the satellite signal path, and correct the inverted STD Into the observation equation of PPP technology, in order to reduce the unknown parameters in the observation equation, the STD calculated by GNSS tropospheric tomography technology is introduced to make up for the shortcomings of PPP technology parameter estimation, and further improve the positioning accuracy and convergence speed of PPP technology.

附图说明Description of drawings

图1为本发明提出的基于GNSS对流层层析技术的PPP改善方法的步骤流程图。Fig. 1 is a flow chart of the steps of the PPP improvement method based on GNSS tropospheric chromatography technology proposed by the present invention.

具体实施方式Detailed ways

下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述。显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Apparently, the described embodiments are only some of the embodiments of the present invention, but not all of them.

参照图1,本方案提供的一种实施例:基于GNSS对流层层析技术的PPP改善方法,包括以下步骤:With reference to Fig. 1, a kind of embodiment that this scheme provides: the PPP improvement method based on GNSS tropospheric tomography technology, comprises the following steps:

S1:基于PPP技术的斜路径对流层延迟STD估算:S1: STD estimation of oblique path tropospheric delay based on PPP technology:

利用PPP技术估计STD时,需综合考虑观测方程的各项误差,确保得到的STD具有较可靠的初始精度,STD计算具体过程如下:When using PPP technology to estimate STD, it is necessary to comprehensively consider the errors of the observation equation to ensure that the obtained STD has a relatively reliable initial accuracy. The specific process of STD calculation is as follows:

S1.1:PPP观测方程建立:S1.1: PPP observation equation establishment:

GNSS信号在从卫星到接收机的传播中会受到各项误差的干扰,这些误差会不同程度的影响定位精度,由于GNSS伪距观测值其本身精度为米级,因此,选择改正精度更高的载波相位观测方程,在充分顾及各项误差的情况下,GNSSPPP技术中载波相位观测方程具体表达如下:The GNSS signal will be interfered by various errors in the propagation from the satellite to the receiver, and these errors will affect the positioning accuracy to varying degrees. Since the GNSS pseudo-range observation value itself has a meter-level accuracy, it is better to choose the one with higher correction accuracy. The carrier phase observation equation, in the case of taking full account of various errors, the carrier phase observation equation in GNSSPPP technology is specifically expressed as follows:

Figure BDA0003987631920000091
Figure BDA0003987631920000091

式中,

Figure BDA0003987631920000101
表示载波相位观测值,i为卫星信号频率,r和s分别代表接收机和卫星系统,
Figure BDA0003987631920000102
表示星地几何距离,tr,i为接收机钟差,
Figure BDA0003987631920000103
表示卫星钟差,
Figure BDA0003987631920000104
为电离层延迟,
Figure BDA0003987631920000105
为对流层延迟,Br,i
Figure BDA0003987631920000106
则分别为接收机端和卫星端的相位硬件延迟,λi为信号频率的波长,
Figure BDA0003987631920000107
表示整周模糊度,
Figure BDA0003987631920000108
为载波相位的残差项;In the formula,
Figure BDA0003987631920000101
Indicates the carrier phase observation value, i is the satellite signal frequency, r and s represent the receiver and satellite system respectively,
Figure BDA0003987631920000102
Indicates the geometric distance between the satellite and the earth, t r,i is the clock error of the receiver,
Figure BDA0003987631920000103
Indicates the satellite clock error,
Figure BDA0003987631920000104
is the ionospheric delay,
Figure BDA0003987631920000105
is the tropospheric delay, B r,i and
Figure BDA0003987631920000106
are the phase hardware delays of the receiver and the satellite respectively, λi is the wavelength of the signal frequency,
Figure BDA0003987631920000107
represents the ambiguity of the whole week,
Figure BDA0003987631920000108
is the residual term of the carrier phase;

S1.2:天顶对流层延迟ZTD估计:S1.2: Zenith Tropospheric Delay ZTD Estimation:

在传统PPP解算中,将基准站坐标作为常数处理,在已知其精确坐标的情况下,通过可联合地表实测气压的Saastamoinen等经验模型计算先验ZHD,将已知坐标与计算出的先验ZHD一同带入观测方程作为起算数据,随后对观测方程线性化,进一步联合多测站多卫星观测信号构建观测方程矩阵,并利用整体最小二乘或卡尔曼滤波方法进行解算,即可得到估计的天顶对流层总延迟ZTD和梯度项;In the traditional PPP solution, the coordinates of the reference station are treated as constants. When the precise coordinates are known, the prior ZHD is calculated through the empirical model such as Saastamoinen, which can be combined with the surface measured air pressure. The known coordinates and the calculated prior The experimental ZHD is brought into the observation equation as the initial calculation data, and then the observation equation is linearized, and the observation equation matrix is further combined with multi-station and multi-satellite observation signals, and the overall least squares or Kalman filter method is used to solve it, and then we can get Estimated zenith tropospheric total delay ZTD and gradient terms;

S1.3:斜路径总延迟STD计算:S1.3: Calculate the total delay STD of the oblique path:

斜路径对流层延迟可表示为天顶方向对流层延迟与映射函数的乘积加上大气水平梯度改正项与其映射函数乘积;The oblique path tropospheric delay can be expressed as the product of the tropospheric delay in the zenith direction and the mapping function plus the product of the atmospheric horizontal gradient correction term and its mapping function;

天顶对流层总延迟ZTD主要包括天顶干延迟ZHD和天顶湿延迟ZWD两部分,通过地表气象参数精确计算ZHD,随后由PPP估计出的ZTD减去ZHD可得天顶湿延迟ZWD,表示如下:The total zenith tropospheric delay ZTD mainly includes two parts: the zenith dry delay ZHD and the zenith wet delay ZWD. The ZHD is accurately calculated through the surface meteorological parameters, and then the zenith wet delay ZWD can be obtained by subtracting the ZHD from the ZTD estimated by the PPP, which is expressed as follows :

ZWD=ZTD-ZHD(2)ZWD=ZTD-ZHD(2)

因此,斜路径对流层延迟STD可表示如下:Therefore, the oblique path tropospheric delay STD can be expressed as follows:

Figure BDA0003987631920000109
Figure BDA0003987631920000109

式中,MFh表示ZHD对应的斜路径干映射函数,MFw表示ZWD对应的斜路径湿映射函数,MFg表示水平梯度映射函数,Gn和Ge分别表示水平梯度改正的北方向分量与东方向分量,

Figure BDA0003987631920000111
为测站到卫星的方位角,εt表示对流层延迟残差;In the formula, MF h represents the oblique path dry mapping function corresponding to ZHD, MF w represents the oblique path wet mapping function corresponding to ZWD, MF g represents the horizontal gradient mapping function, G n and Ge represent the north direction component and the horizontal gradient correction respectively east component,
Figure BDA0003987631920000111
is the azimuth angle from the station to the satellite, ε t is the tropospheric delay residual;

S2:层析区域网格内折射率反演:S2: Refractive index inversion in the grid of the tomographic region:

将利用PPP技术估算的斜路径对流层延迟STD用于GNSS对流层层析模型的构建中,具体步骤如下:The oblique path tropospheric delay STD estimated by PPP technology is used in the construction of the GNSS tropospheric tomography model. The specific steps are as follows:

S2.1:层析模型观测方程构建:S2.1: Construction of tomographic model observation equation:

通过将层析区域在三维方向上划分为多个离散的网格,并对穿过各网格的卫星信号射线上的折射率进行积分,可构建STD与折射率的积分方程:By dividing the tomographic area into multiple discrete grids in the three-dimensional direction, and integrating the refractive index on the satellite signal rays passing through each grid, the integral equation of STD and refractive index can be constructed:

STD=10-6·∫sNds (4)STD=10 -6s Nds (4)

式中,N为大气折射率,包括干折射率和湿折射率,s表示由卫星传播到接收机的信号的路径长度。In the formula, N is the atmospheric refractivity, including dry refractivity and wet refractivity, and s represents the path length of the signal transmitted from the satellite to the receiver.

将上式离散化,则GNSS信号斜路径上对流层总延迟STD表示卫星信号穿过每个网格的截距与该网格内大气折射率乘积之和:Discretizing the above formula, the total tropospheric delay STD on the oblique path of the GNSS signal represents the sum of the intercept of the satellite signal passing through each grid and the product of the atmospheric refractive index in the grid:

STD=∑ijk(aijk·xijk) (5)STD=∑ ijk (a ijk x ijk ) (5)

式中,xijk表示(i,j,k)网格内的待估折射率,aijk表示射线在(i,j,k)网格内的截距,STD表示由PPP技术估算的GNSS卫星信号斜路径上的对流层总延迟估计值。In the formula, x ijk represents the refraction index to be estimated in the (i, j, k) grid, a ijk represents the intercept of the ray in the (i, j, k) grid, and STD represents the GNSS satellite estimated by the PPP technique Estimated total tropospheric delay on signal slant paths.

将研究区域内所有从顶部穿出的信号射线上的STD均用式(5)表示,则构成如下层析观测方程:The STD on all the signal rays passing through the top in the study area are expressed by formula (5), then the following tomographic observation equation is formed:

y=A·x (6)y=A x (6)

式中,y为从研究区域顶部穿出的信号射线上的STD组成的列向量,A为观测方程的系数矩阵,x为未知折射率参数组成的列向量。In the formula, y is a column vector composed of STD on the signal ray passing through the top of the study area, A is the coefficient matrix of the observation equation, and x is a column vector composed of unknown refractive index parameters.

S2.2:层析约束方程构建:S2.2: Construction of tomographic constraint equation:

由于层析区域上空GNSS卫星分布不均以及站点数目不足,使得层析区域很多网格没有射线穿过,导致观测方程的系数矩阵病态,在对层析方程求解时会出现不适定问题,因此,需构建一定的约束方程建立网格在水平和垂直方向上的函数关系,可根据高斯加权函数方法或水平平滑约束的方法构建水平约束方程,依据折射率随高度增加呈指数递减的特性,建立垂直网格内的垂直函数关系。此外,可根据层析区域内无线电探空资料、数值预报再分析资料等建立层析模型的先验信息;Due to the uneven distribution of GNSS satellites over the tomography area and the insufficient number of stations, many grids in the tomography area do not have rays passing through, resulting in ill-conditioned coefficient matrix of the observation equation, and ill-posed problems will occur when solving the tomography equation. Therefore, It is necessary to construct a certain constraint equation to establish the functional relationship between the grid in the horizontal and vertical directions. The horizontal constraint equation can be constructed according to the Gaussian weighted function method or the horizontal smooth constraint method. According to the characteristic that the refractive index decreases exponentially with the increase of height, the vertical Vertical functional relationships within the grid. In addition, the prior information of the tomographic model can be established according to the radiosonde data in the tomographic area, the numerical forecast reanalysis data, etc.;

S2.3:GNSS对流层层析模型构建:S2.3: GNSS tropospheric tomography model construction:

根据上述构建的GNSS层析观测方程、水平约束方程、垂直约束方程以及先验约束方程,建立GNSS对流层层析模型:According to the GNSS tomographic observation equation, horizontal constraint equation, vertical constraint equation and prior constraint equation constructed above, a GNSS tropospheric tomographic model is established:

Figure BDA0003987631920000121
Figure BDA0003987631920000121

式中,H、V和I分别指水平、垂直和先验约束方程的系数矩阵,C表示先验约束信息或通过探空数据等方法统计得到的层析区域内的折射率值。In the formula, H, V, and I refer to the coefficient matrix of the horizontal, vertical, and prior constraint equations, respectively, and C represents the prior constraint information or the refractive index value in the tomographic region obtained by statistical methods such as sounding data.

S2.4:GNSS对流层层析模型解算:S2.4: GNSS tropospheric tomography model solution:

本发明通过奇异值分解法(SingularValueDecomposition,SVD)对GNSS对流层层析模型(7)进行解算,将层析模型系数矩阵

Figure BDA0003987631920000131
分解为:The present invention solves the GNSS troposphere tomography model (7) by singular value decomposition (SingularValueDecomposition, SVD), and the tomography model coefficient matrix
Figure BDA0003987631920000131
Decomposed into:

B=UΛVT(8)B=UΛV T (8)

式中,B∈Rm×n,U∈Rm×m,V∈Rn×n

Figure BDA0003987631920000132
∑=diag(σ12,…,σr),σ1≥σ2≥…≥σr,σi(i=1,2,…,r)为矩阵ATA的特征值的平方根,r为矩阵B的秩(r≤min(m,n)),U是由矩阵AAT的特征向量组成的正交矩阵,V是由矩阵ATA的特征向量组成的正交矩阵。如果矩阵B的广义逆定义为:In the formula, B∈R m×n , U∈R m×m , V∈R n×n ,
Figure BDA0003987631920000132
∑=diag(σ 12 ,…,σ r ), σ 1 ≥σ 2 ≥…≥σ r , σ i (i=1,2,…,r) is the square root of the eigenvalues of the matrix A T A , r is the rank of matrix B (r≤min(m,n)), U is an orthogonal matrix composed of eigenvectors of matrix A T , and V is an orthogonal matrix composed of eigenvectors of matrix A T A. If the generalized inverse of matrix B is defined as:

B-1=VΛ-1UΤ(9)B -1 = VΛ -1 U Τ (9)

那么,线性方程组Bx=L的解,即层析区域内的折射率可表示为:Then, the solution of the linear equation system Bx=L, that is, the refractive index in the tomographic region can be expressed as:

x=B-1L=VΛ-1UΤL(10)x=B -1 L=VΛ -1 U Τ L(10)

S3:利用GNSS对流层层析技术的斜路径对流层延迟STD反演:S3: Oblique-path Tropospheric Delay STD Retrieval Using GNSS Tropospheric Tomography:

GNSS信号斜路径上的对流层总延迟STD可表示为射线穿过层析网格内的折射率与信号在该网格截距的乘积之和,因此,基于GNSS对流层层析技术获取的层析区域每个网格内的折射率,根据卫星信号在每个网格内的截距即可计算卫星信号在该网格内的STD,最后通过累计方式得到卫星信号路径上的STD,具体公式如下:The total tropospheric delay STD on the oblique path of the GNSS signal can be expressed as the sum of the product of the refractive index of the ray passing through the tomographic grid and the signal intercept in the grid. Therefore, the tomographic area obtained based on GNSS tropospheric tomography The refractive index in each grid can calculate the STD of the satellite signal in the grid according to the intercept of the satellite signal in each grid, and finally obtain the STD on the satellite signal path by accumulating, the specific formula is as follows:

Figure BDA0003987631920000133
Figure BDA0003987631920000133

式中,xijk表示由层析技术反演得到的(i,j,k)网格内的折射率,aijk表示射线在(i,j,k)网格内的截距,

Figure BDA0003987631920000134
表示利用GNSS对流层层析技术反演的斜路径对流层总延迟恢复值;In the formula, x ijk represents the refractive index in the (i, j, k) grid obtained by tomographic inversion, a ijk represents the intercept of the ray in the (i, j, k) grid,
Figure BDA0003987631920000134
Indicates the tropospheric total delay recovery value of the oblique path retrieved by GNSS tropospheric tomography technology;

通过上述方法,即可计算出层析区域内不同高度角和方位角的GNSS信号斜路径上的对流层总延迟STD;Through the above method, the total tropospheric delay STD on the oblique path of the GNSS signal at different elevation angles and azimuth angles in the tomographic area can be calculated;

S4:GNSS对流层层析技术改善PPP定位方法:S4: GNSS tropospheric tomography technology improves PPP positioning method:

将GNSS对流层层析技术反演的卫星信号传播路径上的STD加入到PPP观测方程中,对传统PPP观测方程进行改进,则公式(1)可表达成如下形式:The STD on the satellite signal propagation path retrieved by GNSS tropospheric tomography technology is added to the PPP observation equation, and the traditional PPP observation equation is improved. The formula (1) can be expressed as follows:

Figure BDA0003987631920000141
Figure BDA0003987631920000141

式中,

Figure BDA0003987631920000145
表示利用GNSS层析技术反演得到的卫星s到接收机r的斜路径对流层总延迟。In the formula,
Figure BDA0003987631920000145
Indicates the total tropospheric delay of the oblique path from satellite s to receiver r retrieved by GNSS tomographic technology.

因此,顾及GNSS层析技术改善后的PPP观测方程可表达如下:Therefore, taking into account the improved PPP observation equation of GNSS tomographic technology, it can be expressed as follows:

Figure BDA0003987631920000142
Figure BDA0003987631920000142

式中,

Figure BDA0003987631920000143
表示消除了STD项的载波相位观测值,
Figure BDA0003987631920000144
为经过改正后的载波相位残差,其余各项不变;In the formula,
Figure BDA0003987631920000143
Indicates the carrier phase observations with the STD term removed,
Figure BDA0003987631920000144
is the corrected carrier phase residual, and the rest remain unchanged;

最后,进一步对上式线性化,并联合多测站多卫星载波相位观测值构建PPP技术的函数模型和误差模型,利用整体最小二乘或卡尔曼滤波方法对PPP函数模型进行解算,得到高精度的PPP定位结果,并加快PPP收敛速度。Finally, the above formula is further linearized, and the function model and error model of the PPP technology are constructed by combining the carrier phase observations of multiple stations and multiple satellites. Accurate PPP positioning results and speed up PPP convergence.

以上所述,仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,根据本发明的技术方案及其发明构思加以等同替换或改变,都应涵盖在本发明的保护范围之内。The above is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto, any person familiar with the technical field within the technical scope disclosed in the present invention, according to the technical solution of the present invention Any equivalent replacement or change of the inventive concepts thereof shall fall within the protection scope of the present invention.

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节,也不限制该发明仅为所述的具体实施方式。显然,根据本说明书的内容,可作很多的修改和变化。本说明书选取并具体描述这些实施例,是为了更好地解释本发明的原理和实际应用,从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。The preferred embodiments of the invention disclosed above are only to help illustrate the invention. The preferred embodiments are not exhaustive in all detail, nor are the inventions limited to specific embodiments described. Obviously, many modifications and variations can be made based on the contents of this specification. This description selects and specifically describes these embodiments in order to better explain the principle and practical application of the present invention, so that those skilled in the art can well understand and utilize the present invention. The invention is to be limited only by the claims, along with their full scope and equivalents.

Claims (1)

1. The PPP improvement method based on the GNSS troposphere chromatography technology is characterized by comprising the following steps:
s1: estimating the troposphere delay STD of the inclined path based on the PPP technology:
when the STD is estimated by using the PPP technology, the errors of each item of the observation equation need to be comprehensively considered, the obtained STD is ensured to have more reliable initial precision, and the STD calculation specific process is as follows:
s1.1: PPP observation equation establishment:
the GNSS signals are interfered by various errors in the propagation from a satellite to a receiver, the errors affect the positioning accuracy to different degrees, and as the accuracy of the GNSS pseudo-range observation value is meter-level, a carrier phase observation equation with higher correction accuracy is selected, and under the condition of fully considering the various errors, the carrier phase observation equation in the GNSSPPP technology is specifically expressed as follows:
Figure FDA0003987631910000011
in the formula (I), the compound is shown in the specification,
Figure FDA0003987631910000012
representing carrier phase observations, i being the satellite signal frequency, r and s representing the receiver and satellite system respectively,
Figure FDA0003987631910000013
representing the geometric distance of the star to the earth, t r,i In order to be the clock difference of the receiver,
Figure FDA0003987631910000014
which represents the clock error of the satellite or satellites,
Figure FDA0003987631910000015
in order to provide an ionospheric delay,
Figure FDA0003987631910000016
delay of troposphere, B r,i And
Figure FDA0003987631910000017
the phase hardware delay, λ, is then the receiver side and the satellite side, respectively i Is the wavelength of the frequency of the signal,
Figure FDA0003987631910000018
the number of the integer ambiguities is expressed,
Figure FDA0003987631910000019
a residual term that is a carrier phase;
s1.2: zenith tropospheric delay ZTD estimation:
in the traditional PPP resolving, the coordinate of a reference station is used as a constant to be processed, under the condition that the accurate coordinate of the reference station is known, the prior ZHD is calculated through Saastamoinen and other experience models which can jointly measure the air pressure on the ground surface, the known coordinate and the calculated prior ZHD are brought into an observation equation to be used as calculation data, then the observation equation is linearized, an observation equation matrix is further constructed by combining multi-station and multi-satellite observation signals, and the overall least square or Kalman filtering method is utilized to carry out resolving, so that the estimated zenith convection layer total delay ZTD and the gradient term can be obtained;
s1.3: the total delay STD of the inclined path is calculated:
the diagonal path tropospheric delay can be expressed as the product of the zenith direction tropospheric delay and the mapping function plus the product of the atmospheric horizontal gradient correction term and its mapping function;
the zenith troposphere total delay ZTD mainly comprises a zenith dry delay ZHD and a zenith wet delay ZWD, the ZHD is accurately calculated through surface meteorological parameters, and the zenith wet delay ZWD can be obtained by subtracting the ZHD from the ZTD estimated by PPP and is represented as follows:
ZWD=ZTD-ZHD(2)
thus, the diagonal path tropospheric delay STD can be expressed as follows:
Figure FDA0003987631910000021
in the formula, MF h Representing diagonal path stem mapping function, MF, corresponding to ZHD w Representing inclined path wet mapping function, MF, corresponding to ZWD g Representing a horizontal gradient mapping function, G n And G e Respectively representing the north component and the east component of the horizontal gradient correction,
Figure FDA0003987631910000022
for measuring station-to-satellite azimuth angle, e t Represents the tropospheric delay residual;
s2: refractive index inversion in the tomographic region grid:
the method is characterized in that the inclined path troposphere delay STD estimated by the PPP technology is used for constructing a GNSS troposphere chromatographic model, and comprises the following specific steps:
s2.1: constructing a chromatographic model observation equation:
by dividing the tomographic region into a plurality of discrete grids in the three-dimensional direction and integrating the refractive index on the satellite signal ray passing through each grid, an integral equation of the STD with the refractive index can be constructed:
STD=10 -6 ·∫ s Nds(4)
where N is the index of refraction of the atmosphere, including the dry and wet indices of refraction, and s represents the path length of the signal propagating from the satellite to the receiver.
Discretizing the above equation, the total tropospheric delay on the GNSS signal slope path, STD, represents the sum of the product of the intercept of the satellite signal through each mesh and the index of refraction of the atmosphere within that mesh:
STD=∑ ijk (a ijk ·x ijk ) (5)
in the formula, x ijk Denotes the index of refraction to be estimated, a, within the (i, j, k) grid ijk Represents the ray intercept within the (i, j, k) grid and STD represents the estimate of the total tropospheric delay on the GNSS satellite signal diagonal path estimated by the PPP technique.
If the STD of all the signal rays emerging from the top in the investigation region is expressed by equation (5), the following tomographic equation is formed:
y=A·x(6)
where y is the column vector consisting of the STD on the signal rays passing out from the top of the investigation region, a is the coefficient matrix of the observation equation, and x is the column vector consisting of unknown refractive index parameters.
S2.2: constructing a chromatographic constraint equation:
due to uneven distribution of GNSS satellites above a chromatographic region and insufficient number of stations, a plurality of grids in the chromatographic region do not have rays to pass through, so that a coefficient matrix of an observation equation is ill-conditioned, and an ill-defined problem can occur when the chromatographic equation is solved. In addition, the prior information of the chromatographic model can be established according to radio sounding data, numerical prediction reanalysis data and the like in the chromatographic region;
s2.3: construction of a GNSS troposphere chromatography model:
according to the established GNSS chromatographic observation equation, the established horizontal constraint equation, the established vertical constraint equation and the established prior constraint equation, establishing a GNSS troposphere chromatographic model:
Figure FDA0003987631910000041
in the formula, H, V and I respectively refer to coefficient matrixes of horizontal, vertical and prior constraint equations, and C represents prior constraint information or a refractive index value in a chromatography region obtained through statistics of sounding data and other methods.
S2.4: GNSS tropospheric tomographic model solution:
the method solves the GNSS troposphere chromatographic model (7) by a Singular Value Decomposition (SVD) method, and solves the chromatographic model coefficient matrix
Figure FDA0003987631910000042
The decomposition is as follows:
B=UΛV T (8)
in the formula, B is belonged to R m×n ,U∈R m×m ,V∈R n×n
Figure FDA0003987631910000043
∑=diag(σ 12 ,…,σ r ),σ 1 ≥σ 2 ≥…≥σ r ,σ i (i =1,2, \ 8230;, r) is the matrix A T The square root of the eigenvalues of A, r the rank of matrix B (r ≦ min (m, n)), U is represented by matrix AA T Is an orthogonal matrix composed of the eigenvectors of, V is a matrix A T And the feature vectors of A form an orthogonal matrix. If the generalized inverse of matrix B is defined as:
B -1 =VΛ -1 U Τ (9)
the solution of the system of linear equations Bx = L, i.e. the refractive index in the tomographic region, can be expressed as:
x=B -1 L=VΛ -1 U Τ L(10)
s3: performing inclined path troposphere delay STD inversion by utilizing GNSS troposphere chromatography technology:
the total tropospheric delay STD on the GNSS signal slant path may be expressed as the sum of the product of the refractive index of the ray through the tomographic grid and the signal intercept at the grid. Therefore, based on the refractive index in each grid of the tomographic region obtained by the GNSS troposphere tomography, the STD of the satellite signal in each grid can be calculated according to the intercept of the satellite signal in each grid, and finally the STD on the satellite signal path is obtained in an accumulation manner, with the following specific formula:
Figure FDA0003987631910000051
in the formula, x ijk Denotes the refractive index, a, in the (i, j, k) grid obtained by the inversion of the tomography technique ijk Represents the intercept of the ray in the (i, j, k) grid,
Figure FDA0003987631910000052
representing the total delay recovery value of the troposphere of the inclined path inverted by utilizing a GNSS troposphere chromatography technology;
by the method, the total troposphere delay STD on the GNSS signal inclined path at different altitude angles and azimuth angles in the chromatographic region can be calculated;
s4: the GNSS troposphere chromatography technology improves PPP positioning method:
adding STD on a satellite signal propagation path inverted by a GNSS troposphere tomography technology into a PPP observation equation, and improving the traditional PPP observation equation, the formula (1) can be expressed as follows:
Figure FDA0003987631910000053
in the formula (I), the compound is shown in the specification,
Figure FDA0003987631910000054
the total tropospheric delay of the oblique path from the satellite s to the receiver r obtained by inversion using GNSS tomography is shown.
Therefore, the PPP observation equation after considering the improvement of GNSS tomography can be expressed as follows:
Figure FDA0003987631910000061
in the formula (I), the compound is shown in the specification,
Figure FDA0003987631910000062
a carrier phase observation representing the removal of the STD term,
Figure FDA0003987631910000063
the carrier phase residual error after being corrected is unchanged in other items;
and finally, the above formula is further linearized, a function model and an error model of the PPP technology are established by combining the multi-station multi-satellite carrier phase observation values, the PPP function model is solved by using an overall least square or Kalman filtering method, a high-precision PPP positioning result is obtained, and the PPP convergence speed is accelerated.
CN202211570241.0A 2022-12-08 2022-12-08 PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology Withdrawn CN115755115A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211570241.0A CN115755115A (en) 2022-12-08 2022-12-08 PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211570241.0A CN115755115A (en) 2022-12-08 2022-12-08 PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology

Publications (1)

Publication Number Publication Date
CN115755115A true CN115755115A (en) 2023-03-07

Family

ID=85344638

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211570241.0A Withdrawn CN115755115A (en) 2022-12-08 2022-12-08 PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology

Country Status (1)

Country Link
CN (1) CN115755115A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116340710A (en) * 2023-05-30 2023-06-27 中国科学院精密测量科学与技术创新研究院 Neutral atmospheric skew delay calculation method based on layered rapid three-dimensional ray tracing
CN117970382A (en) * 2024-03-29 2024-05-03 中国科学院国家空间科学中心 A GNSS simulation test method and system

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116340710A (en) * 2023-05-30 2023-06-27 中国科学院精密测量科学与技术创新研究院 Neutral atmospheric skew delay calculation method based on layered rapid three-dimensional ray tracing
CN116340710B (en) * 2023-05-30 2023-09-12 中国科学院精密测量科学与技术创新研究院 Neutral atmospheric skew delay calculation method based on layered rapid three-dimensional ray tracing
CN117970382A (en) * 2024-03-29 2024-05-03 中国科学院国家空间科学中心 A GNSS simulation test method and system

Similar Documents

Publication Publication Date Title
Chen et al. Voxel-optimized regional water vapor tomography and comparison with radiosonde and numerical weather model
CA3036928C (en) Localization and tracking using location, signal strength, and pseudorange data
Xia et al. GNSS troposphere tomography based on two-step reconstructions using GPS observations and COSMIC profiles
CN103760572B (en) A kind of single-frequency PPP ionosphere based on region CORS method of weighting
CN107422351A (en) A kind of GNSS decimeter grade Differential positioning methods based on virtual grid
CN107797126A (en) BDS/GPS broadcast type network RTK algorithms based on Star Network
CN105629263A (en) Troposphere atmosphere delay error correction method and correction system
CN108549095B (en) Non-differential parallel enhancement method and system for regional CORS network
Zhao et al. An improved GNSS tropospheric tomography method with the GPT2w model
CN115755115A (en) PPP (Point-to-Point protocol) improvement method based on GNSS troposphere chromatography technology
CN108196284B (en) GNSS network data processing method for fixing single-difference ambiguity between satellites
JP2010528320A (en) Reduction of distance-dependent error in real-time kinematic (RTK) positioning
CN109541647B (en) GNSS Multipath Effect Correction Method Based on Half-Celestial Spherical Grid Point Model
CN104656108A (en) Sparse reference station network zenith troposphere delay modeling method considering elevation difference
CN106646564A (en) Navigation enhancing method based on low track satellite
CN104375157A (en) Inertial navigation assisted Big Dipper single-frequency whole-cycle ambiguity calculation method under short baseline condition
CN115616643B (en) City area modeling auxiliary positioning method
CN115982564A (en) Method and device for calculating electron density of regional ionized layer and computer equipment
CN104007479A (en) Ionized layer chromatography technology and ionized layer delay correction method based on multi-scale subdivision
CN105510945A (en) PPP positioning method applied to satellite navigation landing outfield detection
CN103543454B (en) A kind of Satellite Orbit Determination system being embedded in wireless network
CN111983650A (en) GNSS-based high-precision time transfer method
Zhao et al. An optimal tropospheric tomography approach with the support of an auxiliary area
Zhao et al. A troposphere tomography method considering the weighting of input information
Zhang et al. AN IMPROVED TROPOSPHERIC TOMOGRAPHY METHOD BASED ON THE DYNAMIC NODE PARAMETERIZED ALGORITHM.

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
WW01 Invention patent application withdrawn after publication

Application publication date: 20230307

WW01 Invention patent application withdrawn after publication