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

Abstract

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.

Description

PPP improvement method based on GNSS tropospheric tomography technology

technical field

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

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.

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.

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.

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

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:

The PPP improvement method based on GNSS tropospheric tomography technology comprises the following steps:

S1: STD estimation of oblique path tropospheric delay based on PPP technology:

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 observation equation establishment:

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:

表示载波相位观测值，i为卫星信号频率，r和s分别代表接收机和卫星系统，

表示星地几何距离，t_r,i为接收机钟差，

表示卫星钟差，

为电离层延迟，

为对流层延迟，B_r,i和

则分别为接收机端和卫星端的相位硬件延迟，λ_i为信号频率的波长，

表示整周模糊度，

为载波相位的残差项；In the formula,

Indicates the carrier phase observation value, i is the satellite signal frequency, r and s represent the receiver and satellite system respectively,

Indicates the geometric distance between the satellite and the earth, t _r,i is the clock error of the receiver,

Indicates the satellite clock error,

is the ionospheric delay,

is the tropospheric delay, B _r,i and

are the phase hardware delays of the receiver and the satellite respectively, _λi is the wavelength of the signal frequency,

represents the ambiguity of the whole week,

is the residual term of the carrier phase;

S1.2: Zenith Tropospheric Delay ZTD Estimation:

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: 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;

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)

Therefore, the oblique path tropospheric delay STD can be expressed as follows:

为测站到卫星的方位角，ε_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,

is the azimuth angle from the station to the satellite, ε _t is the tropospheric delay residual;

S2: Refractive index inversion in the grid of the tomographic region:

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: Construction of tomographic model observation equation:

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 ∫ _s Nds (4)

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.

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 (a _ijk x _ijk ) (5)

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.

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)

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: Construction of tomographic constraint equation:

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 tropospheric tomography model construction:

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:

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 tropospheric tomography model solution:

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

Decomposed into:

B＝UΛV ^T (8)

∑＝diag(σ₁,σ₂,…,σ_r)，σ₁≥σ₂≥…≥σ_r，σ_i(i＝1,2,…,r)为矩阵A^TA的特征值的平方根，r为矩阵B的秩(r≤min(m,n))，U是由矩阵AA^T的特征向量组成的正交矩阵，V是由矩阵A^TA的特征向量组成的正交矩阵。如果矩阵B的广义逆定义为：In the formula, B∈R ^m×n , U∈R ^m×m , V∈R ^n×n ,

∑=diag(σ ₁ ,σ ₂ ,…,σ _r ), σ ₁ ≥σ ₂ ≥…≥σ _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Λ ^-1 U ^Τ (9)

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 ^-1 L=VΛ ^-1 U ^Τ L(10)

S3: Oblique-path Tropospheric Delay STD Retrieval Using GNSS Tropospheric Tomography:

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:

表示利用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,

Indicates the tropospheric total delay recovery value of the oblique path retrieved by GNSS tropospheric tomography technology;

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 tropospheric tomography technology improves PPP positioning method:

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:

表示利用GNSS层析技术反演得到的卫星s到接收机r的斜路径对流层总延迟。In the formula,

Indicates the total tropospheric delay of the oblique path from satellite s to receiver r retrieved by GNSS tomographic technology.

Therefore, taking into account the improved PPP observation equation of GNSS tomographic technology, it can be expressed as follows:

表示消除了STD项的载波相位观测值，

为经过改正后的载波相位残差，其余各项不变；In the formula,

Indicates the carrier phase observations with the STD term removed,

is the corrected carrier phase residual, and the rest remain unchanged;

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.

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;

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

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.

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: STD estimation of oblique path tropospheric delay based on PPP technology:

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 observation equation establishment:

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:

表示载波相位观测值，i为卫星信号频率，r和s分别代表接收机和卫星系统，

表示星地几何距离，t_r,i为接收机钟差，

表示卫星钟差，

为电离层延迟，

为对流层延迟，B_r,i和

则分别为接收机端和卫星端的相位硬件延迟，λ_i为信号频率的波长，

表示整周模糊度，

为载波相位的残差项；In the formula,

Indicates the carrier phase observation value, i is the satellite signal frequency, r and s represent the receiver and satellite system respectively,

Indicates the geometric distance between the satellite and the earth, t _r,i is the clock error of the receiver,

Indicates the satellite clock error,

is the ionospheric delay,

is the tropospheric delay, B _r,i and

are the phase hardware delays of the receiver and the satellite respectively, _λi is the wavelength of the signal frequency,

represents the ambiguity of the whole week,

is the residual term of the carrier phase;

S1.2: Zenith Tropospheric Delay ZTD Estimation:

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: 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;

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)

Therefore, the oblique path tropospheric delay STD can be expressed as follows:

为测站到卫星的方位角，ε_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,

is the azimuth angle from the station to the satellite, ε _t is the tropospheric delay residual;

S2: Refractive index inversion in the grid of the tomographic region:

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: Construction of tomographic model observation equation:

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 ∫ _s Nds (4)

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.

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 (a _ijk x _ijk ) (5)

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.

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)

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: Construction of tomographic constraint equation:

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 tropospheric tomography model construction:

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:

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 tropospheric tomography model solution:

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

Decomposed into:

B＝UΛV ^T (8)

∑＝diag(σ₁,σ₂,…,σ_r)，σ₁≥σ₂≥…≥σ_r，σ_i(i＝1,2,…,r)为矩阵A^TA的特征值的平方根，r为矩阵B的秩(r≤min(m,n))，U是由矩阵AA^T的特征向量组成的正交矩阵，V是由矩阵A^TA的特征向量组成的正交矩阵。如果矩阵B的广义逆定义为：In the formula, B∈R ^m×n , U∈R ^m×m , V∈R ^n×n ,

∑=diag(σ ₁ ,σ ₂ ,…,σ _r ), σ ₁ ≥σ ₂ ≥…≥σ _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Λ ^-1 U ^Τ (9)

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 ^-1 L=VΛ ^-1 U ^Τ L(10)

S3: Oblique-path Tropospheric Delay STD Retrieval Using GNSS Tropospheric Tomography:

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:

表示利用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,

Indicates the tropospheric total delay recovery value of the oblique path retrieved by GNSS tropospheric tomography technology;

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 tropospheric tomography technology improves PPP positioning method:

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:

表示利用GNSS层析技术反演得到的卫星s到接收机r的斜路径对流层总延迟。In the formula,

Indicates the total tropospheric delay of the oblique path from satellite s to receiver r retrieved by GNSS tomographic technology.

Therefore, taking into account the improved PPP observation equation of GNSS tomographic technology, it can be expressed as follows:

表示消除了STD项的载波相位观测值，

为经过改正后的载波相位残差，其余各项不变；In the formula,

Indicates the carrier phase observations with the STD term removed,

is the corrected carrier phase residual, and the rest remain unchanged;

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:

in the formula (I), the compound is shown in the specification,

representing carrier phase observations, i being the satellite signal frequency, r and s representing the receiver and satellite system respectively,

representing the geometric distance of the star to the earth, t _r,i In order to be the clock difference of the receiver,

which represents the clock error of the satellite or satellites,

in order to provide an ionospheric delay,

delay of troposphere, B _r,i And

the phase hardware delay, λ, is then the receiver side and the satellite side, respectively _i Is the wavelength of the frequency of the signal,

the number of the integer ambiguities is expressed,

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:

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,

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:

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

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 ，

∑＝diag(σ ₁ ,σ ₂ ,…,σ _r )，σ ₁ ≥σ ₂ ≥…≥σ _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:

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,

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:

in the formula (I), the compound is shown in the specification,

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:

in the formula (I), the compound is shown in the specification,

a carrier phase observation representing the removal of the STD term,

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.

Images