CN105510945A - PPP positioning method applied to satellite navigation landing outfield detection - Google Patents

Abstract

The invention provides a PPP positioning method applied to satellite navigation landing outfield detection. On the basis of mobile communication data transmission, a precise ephemeris and a precise clock difference of an IGS network station are obtained in real time, and a broadcast ephemeris and observation data are obtained by a receiver antenna; an influence on an ionized layer is eliminated based on double-frequency ionized-layer-free zero-difference observation value combination, and ambiguity and coordinate accumulation is carried out by using extended Kalman filter, and rapid convergence is carried out; and ambiguity information and a covariance matrix that are obtained each time are introduced into a Lambda algorithm to carry out resolution, thereby obtaining a correct ambiguity value and a correct receiver coordinate. According to the invention, a requirement that the cm-level positioning precision can be reached by the signal receiver can be realized. Meanwhile, because the precise track and the precise clock difference that are issued by the IGS network station are obtained in an on-line mode in real time and the high-precision positioning is realized in real time, the practical engineering requirement can be met well.

Description

PPP positioning method applied to satellite navigation landing outfield detection

Technical Field

The invention belongs to the field of satellite navigation, and relates to a high-precision positioning method.

Background

Satellite navigation has the advantages of all weather, wide coverage, low cost and the like, and has become the most important navigation means in a navigation system. In satellite navigation positioning, two positioning modes of pseudo-range point positioning and differential positioning are mainly adopted, the pseudo-range point positioning can be realized by only one receiver, the real-time performance is good, the single positioning precision is generally dozens of meters, the requirement of rough position positioning can be well met, but the precision requirement cannot meet the actual requirements of surveying and mapping, detection, calibration and the like. Obtaining a high-precision coordinate position, wherein the adopted technology is mainly a relative positioning technology, at least 2 receivers are required for relative positioning, and 1 receiver is arranged on a high-precision known point and is called a reference station; and another 1 receiver is placed on an unknown point to be mapped and is called a rover, and then the two receivers form a base line by receiving common GPS satellite signals, and perform base line calculation to calculate the position coordinates of the rover. The relative positioning technology has high positioning precision and good real-time performance, can well meet the requirement of high-precision position, but along with the increase of the length of the base line, the error correlation between the base station and the mobile station is weakened, and the positioning precision of the mobile station is reduced. When no reference point exists in a measuring area needing to be measured or a CORS system cannot cover the measuring area, long-distance joint measurement is needed to introduce a known point, the workload is very large, and the positioning accuracy is reduced. The precise single-point positioning technology of the GPS is a brand new positioning mode developed in recent years, precise ephemeris and precise clock error provided by IGS are adopted, a user receives ephemeris and observation data by using a single GPS dual-frequency receiver, and high-precision positioning of any position within thousands of square kilometers and even the global range can be realized by subsequent calculation.

Disclosure of Invention

In order to overcome the defects of the prior art, the invention provides a method for obtaining high-precision measurement in real time by using a single dual-frequency receiver.

The technical scheme adopted by the invention for solving the technical problem comprises the following steps:

(1) acquiring a precise ephemeris and a precise clock error of an IGS website in real time, and receiving a broadcast ephemeris and observation data;

(2) constructing a dual-frequency L_iPseudo-range observation value P (L) of_i) And carrier phase observation phi (L)_i) Wherein, i is 1,2,P(L_i)＝ρ+c(dt-dT)+d_orb+d_trop+d_ion/Li+d_mpath/P(Li)+(P(L_i) ρ is the geometric distance from the satellite to the receiver, (X)_s,Y_s,Z_s) Is the satellite's coordinates at the signal transmission time T, (x, y, z) are the receiver coordinates at the signal reception time T ═ T + ρ/c,c is the speed of light, dT is the satellite clock difference, dT is the receiver clock difference, d_orbIs the satellite orbital error, d_tropIs tropospheric delay, d_ion/LiIs ionospheric delay on Li, lambda_iIs the wavelength of Li, N_iIs the integer ambiguity of Li,is the initial phase of the receiver oscillator,is the initial phase of the satellite oscillator, d_mpath/P(Li)Is the multipath effect of pseudorange measurements on Li, d_mpath/φ(Li)Multipath effects of carrier phase measurements on Li, () measurement noise;

(3) calculating non-differential combined observed values of double-frequency ionosphere-free including ionosphere-free pseudo range observed value P_IFAnd ionosphere-free phase observation value phi_IF(ii) a Wherein,

$P_{IF}=\frac{f_1^{2}P\left(L_1\right)-f_2^{2}P\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+d_{orb}+M\·zpd+d_{mpath/P(L1+L2)}+\mathrm{\ϵ}\left(P\right(L1+L2\left)\right),$

the product of the tropospheric delay zpd and the mapping function M is the tropospheric delay d_trop；

(4) Construction System x_k＝(X_k,Y_k,Z_k,zpd_k,N_k) (ii) a Constructing a measurement equation y_k＝(P_IF,_k,Φ_IF,k)＝H(x_k)x_k(ii) a Build system estimation equationEstimation of equation covariance for a formation system $P_k(+)=(I-K_{k}H({\hat{x}}_k(-)\left)\right)P_k(-);$ Constructing a gain equation $K_k=P_k(-)H\left(\hat{x}\right(-\left)\right){\left(H\right({\hat{x}}_k(-)\left)P_k\right(-\left)H{\left({\hat{x}}_k(-)\right)}^T\right)}^{-1};$ Wherein, X_k,Y_k,Z_kRepresenting the receiver position at time k, zpd_kTropospheric delay at time k, N_kIs the ambiguity at time k; p_IF,kRepresenting ionospheric-free pseudorange observations at time k, phi_IF,kRepresents the ionospheric-free phase observation at time k, H (x)_k) Is a coefficient matrix of the measurement equation;representing the system estimate at time k-1,representing the system estimate at time K, K_kIs a gain matrix, y_kTo measure the observed value; i is the identity matrix, P_k(-) at time k-1Covariance of (P)_k(+) is the time kThe covariance of (a); ambiguity and coordinate accumulation are carried out by adopting extended Kalman filtering, the ambiguity obtained each time and covariance matrix thereof are substituted into Lambda algorithm, and coordinate X of the receiver is obtained by calculation_k,Y_k,Z_k。

The invention has the beneficial effects that:

(1) compared with the prior method of relying on pseudo-range single-point positioning or differential positioning, the PPP positioning method for external field detection provided by the invention meets the requirement that a single receiver can achieve cm-level positioning accuracy;

(2) compared with the past method of precise single-point positioning calculation after work, the PPP positioning method for external field detection provided by the invention obtains the precise orbit and the precise clock error published by the IGS website on line in real time, realizes real-time high-precision positioning, and better meets the actual engineering requirements.

Drawings

FIG. 1 is a schematic diagram of a precise single-point positioning calculation data acquisition process;

FIG. 2 is a diagram illustrating differences in the X direction of PPP positioning;

FIG. 3 is a diagram illustrating differences in the Y-direction of PPP positioning;

FIG. 4 is a diagram illustrating differences in the Z-direction of PPP positioning;

FIG. 5 is a diagram of pseudorange single point location X direction difference;

FIG. 6 is a diagram of pseudorange single point location Y direction difference;

fig. 7 is a diagram of differences in the pseudorange single point position Z direction.

Detailed Description

The present invention will be further described with reference to the following drawings and examples, which include, but are not limited to, the following examples.

According to the method, the precise orbit and the precise clock error issued by the IGS website are obtained on line in real time, precise single-point calculation is carried out at the receiver end, and the cm-level positioning precision is obtained in real time. The specific contents are as follows:

(1) precise track and precise clock error published by IGS website are obtained on line in real time

Referring to fig. 1, the real-time single-point GPS positioning solution needs to obtain the real-time ephemeris and precision clock error of the IGS website, a communication module is added in front of a receiver positioning solution module, the precision ephemeris and clock error of the IGS website are obtained in real time through mobile communication data transmission, the broadcast ephemeris and observation data are obtained through a receiver antenna, and then the real-time ephemeris, the precision clock error, the broadcast ephemeris and the observation data are sent to the receiver positioning solution module for solution.

(2) Constructing dual-frequency pseudorange observations and carrier phase observations

Adopting high-precision GPS satellite ephemeris and satellite clock error and dual-frequency carrier phase observed values, and adopting a non-differential model to carry out high-precision single-point positioning, wherein the equation is as follows:

P(L_i)＝ρ+c(dt-dT)+d_orb+d_trop+d_ion/Li+d_mpath/P(Li)+(P(L_i))(2)

Φ(L_i) Is the carrier phase observation (m) on Li, i ═ 1, 2; p (L)_i) Is L_iA pseudorange observation (m); ρ is the geometric distance from the satellite to the receiver, (X)_s,Y_s,Z_s) Is the coordinate of the satellite transmission time T, (x, y, z) is the receiver coordinate of the signal reception time T ═ T + ρ/c,c is the speed of light (m/s); dt is the satellite clock error(s); dT is the receiver clock difference(s); d_orbIs the satellite orbital error (m); d_tropIs tropospheric delay (m); d_ion/LiIonospheric delay (m) over Li; lambda [ alpha ]_iIs the wavelength (m) of Li; n is a radical of_iIs the integer ambiguity (week) of Li;is the initial phase of the receiver oscillator;is the initial phase of the satellite oscillator; d_mpath/P(Li)Multipath effects (m) of pseudorange measurements over Li; d_mpath/φ(Li)Multipath effects (m) that are carrier phase measurements on Li; (. cndot.) is measurement noise.

(3) Calculating non-difference combined observed value of double-frequency non-ionosphere

Delaying the troposphere by d using the precision satellite ephemeris and precision satellite clock error of the IGS_tropRepresenting the product of the pair-flow layer delay zpd and the mapping function M. The combination of dual-frequency ionosphere-free observation values is adopted to eliminate the first-order influence of the ionosphereThe mathematical model of point location will become:

$\begin{array}{c}{P}_{IF}=\frac{f_1^{2}P\left(L_1\right)-f_2^{2}P\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+d_{orb}\\ +M\·zpd+d_{mpath/P(L1+L2)}+\mathrm{\ϵ}\left(P\right(L1+L2\left)\right)\end{array}---\left(3\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}=\frac{f_1^2\mathrm{\φ}\left(L_1\right)-f_2^2\mathrm{\φ}\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+M\·zpd+\mathrm{\λ}N\\ +d_{mpath/\mathrm{\φ}(L1+L2)}+\mathrm{\ϵ}\left(\mathrm{\φ}\right(L1+L2\left)\right)\end{array}---\left(4\right)$

(4) extended Kalman Filter (EKF) based positioning

And ambiguity and coordinate accumulation are carried out by adopting Extended Kalman Filtering (EKF), and fast convergence is achieved. The construction equation is as follows:

the system equation:

x_k＝(X_k,Y_k,Z_k,zpd_k,N_k)(5)

X_k,Y_k,Z_krepresenting the receiver position at time k, zpd_kTropospheric delay at time k, N_kIs the ambiguity at time k.

The measurement equation:

y_k＝(P_IF,k,Φ_IF,k)＝H(x_k)x_k(6)

P_IF,krepresenting ionospheric-free pseudorange observations at time k, phi_IF,kRepresents the ionospheric-free phase observation at time k, H (x)_k) Is a coefficient matrix of the measurement equation.

The system estimation equation:

${\hat{x}}_k(+)={\hat{x}}_k(-)+K_k(y_k-H({\hat{x}}_k(-)\left)\right({\hat{x}}_k(-)\left)\right)---\left(7\right)$

representing the system estimate at time k-1,representing the system estimate at time K, K_kIs a gain matrix, y_kTo measure the observed value.

The system estimates the equation covariance:

$P_k(+)=(I-K_{k}H({\hat{x}}_k(-)\left)\right)P_k(-)---\left(8\right)$

i is the identity matrix, P_k(-) at time k-1Covariance of (P)_k(+) is the time kThe covariance.

Gain equation:

$K_k=P_k(-)H\left(\hat{x}\right(-\left)\right){\left(H\right({\hat{x}}_k(-)\left)P_k\right(-\left)H{\left({\hat{x}}_k(-)\right)}^T\right)}^{-1}---\left(9\right)$

the coordinates X of the receiver can be obtained by the calculation of the expressions (5) to (9)_k,Y_k,Z_k。

An embodiment of the invention comprises the following steps:

(1) precise track and precise clock error published by IGS website are obtained on line in real time

Referring to fig. 1, real-time single-point positioning calculation needs to obtain an IGS website real-time precise ephemeris and clock error, a communication module is added in front of a receiver positioning calculation module, the IGS website ultrafast ephemeris and clock error are obtained in real time through mobile communication data transmission, broadcast ephemeris and observation data are obtained through a receiver antenna, and then the observed quantities are sent to the receiver positioning calculation module for calculation.

(2) Constructing dual-frequency code and phase observations

With a dual frequency GPS receiver, the code and phase observations between the receiver and the satellite at L1 and L2 can be expressed as:

P(L_i)＝ρ+c(dt-dT)+d_orb+d_trop+d_ion/Li+d_mpath/P(Li)+(P(L_i))(10)

(3) calculating non-difference combined observed value of double-frequency non-ionosphere

The ionospheric delay error is one of the main errors of non-differential phase precision single-point positioning accuracy, and the generated distance influence can even reach hundreds of meters when the satellite elevation angle is low. The first-order influence of the ionosphere is eliminated by combining the two-frequency non-ionosphere non-difference observation values, and the observation equation is changed into:

$\begin{array}{c}{P}_{IF}=\frac{f_1^{2}P\left(L_1\right)-f_2^{2}P\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}+c(dt-dT)+d_{orb}\\ +d_{trop}+d_{mpath/P(L1+L2)}+\mathrm{\ϵ}\left(P\right(L1+L2\left)\right)\end{array}---\left(12\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}=\frac{f_1^2\mathrm{\φ}\left(L_1\right)-f_2^2\mathrm{\φ}\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}+c(dt-dT)+d_{orb}\\ +d_{trop}+\frac{{\mathrm{cf}}_1{N}_1-{\mathrm{cf}}_2{N}_2}{f_1^2-f_2^2}+d_{mpath/\mathrm{\φ}(L1+L2)}+\mathrm{\ϵ}\left(\mathrm{\φ}\right(L1+L2\left)\right)\end{array}---\left(13\right)$

delaying the troposphere by d using the precision satellite ephemeris and precision satellite clock error of the IGS_tropRepresenting the product of the streamlining delay zpd and its mapping function M. The mathematical model of single point positioning will become:

$\begin{array}{c}{P}_{IF}=\frac{f_1^{2}P\left(L_1\right)-f_2^{2}P\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+d_{orb}\\ +M\·zpd+d_{mpath/P(L1+L2)}+\mathrm{\ϵ}\left(P\right(L1+L2\left)\right)\end{array}---\left(14\right)$

$\begin{array}{c}{\mathrm{\Φ}}_{IF}=\frac{f_1^2\mathrm{\φ}\left(L_1\right)-f_2^2\mathrm{\φ}\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+M\·zpd+\mathrm{\λ}N\\ +d_{mpath/\mathrm{\φ}(L1+L2)}+\mathrm{\ϵ}\left(\mathrm{\φ}\right(L1+L2\left)\right)\end{array}---\left(15\right)$

(4) calculating ambiguities and coordinates

And ambiguity and coordinate accumulation are carried out by adopting Extended Kalman Filtering (EKF), and fast convergence is achieved. The construction equation is as follows:

x_k＝(x_k,y_k,z_k,zpd_k,N_k)(16)

y_k＝(P_IF,k,Φ_IF,k)＝H(x_k)x_k(17)

${\hat{x}}_k(+)={\hat{x}}_k(-)+K_k(y_k-H({\hat{x}}_k(-)\left)\right({\hat{x}}_k(-)\left)\right)---\left(18\right)$

$P_k(+)=(I-K_{k}H({\hat{x}}_k(-)\left)\right)P_k(-)---\left(19\right)$

$K_k=P_k(-)H\left(\hat{x}\right(-\left)\right){\left(H\right({\hat{x}}_k(-)\left)P_k\right(-\left)H{\left({\hat{x}}_k(-)\right)}^T\right)}^{-1}---\left(20\right)$

substituting the ambiguity and covariance matrix obtained each time into Lambda algorithm for fixing, and then fixing the coordinate value X of the receiver_k,Y_k,Z_kAnd degree of ambiguity N_kAnd (4) carrying in the formulas (16) to (20) to finish rapid convergence, and obtaining a cm-level positioning result.

The following describes a specific implementation of the calibration and positioning method of the present invention, taking the calibration and positioning result in a certain area as an example, and comparing the result with the pseudorange single-point positioning result. The method provided by the invention is utilized, and the specific steps of positioning in a certain actual calibration are as follows:

(1) mounting receivers at known points

The receiver is set up at a known point R (-3050576.283, 4551789.819, 3192861.287), turned on and continuously observed, and observations and positioning results, including PPP positioning and pseudorange single point positioning results, are recorded.

(2) Online real-time precision track and precision clock error parallel positioning calculation

Through mobile communication data transmission, the receiver acquires the ultrafast ephemeris and clock error of the IGS website in real time, acquires the broadcast ephemeris and observation data through the antenna of the receiver, sends the observation quantities to the receiver positioning resolving module for resolving, and acquires the position of the calibration point in real time.

(3) PPP positioning and pseudo range single point positioning

After continuous observation, the point coordinate obtained by PPP final convergence is P (-3050576.249, 4551789.876, 3192861.307), and the coordinate of the pseudo-range single-point positioning point is Q (-3050572.278, 4551794.041, 3192864.773).

(4) Alignment error comparison

Comparing errors of the PPP location and pseudorange location with the known point R (-3050576.283, 4551789.819, 3192861.287) as a true value, subtracting: the errors of PPP positioning in XYZ three directions are: (0.034, 0.057, 0.020); errors of pseudo-range single-point positioning in three directions of XYZ are respectively: (3.971,4.165,3.486). From the difference, it can be seen that: PPP positioning error is in the order of cm and pseudorange single point positioning error is around 5 m.

(5) Convergence comparison

As can be seen from fig. 2, 3 and 4, PPP can be converged to cm level after about 10min, while as can be seen from fig. 5, 6 and 7, pseudorange single-point positioning has no convergence due to large error.

Claims (1)

1. A PPP positioning method applied to satellite navigation landing external field detection is characterized by comprising the following steps:

(1) acquiring a precise ephemeris and a precise clock error of an IGS website in real time, and receiving a broadcast ephemeris and observation data;

(2) constructing a dual-frequency L_iPseudo-range observation value P (L) of_i) And carrier phase observation phi (L)_i) Wherein, i is 1,2,P(L_i)＝ρ+c(dt-dT)+d_orb+d_trop+d_ion/Li+d_mpath/P(Li)+(P(L_i) ρ is the geometric distance from the satellite to the receiver, (X)_s,Y_s,Z_s) Is the satellite's coordinates at the signal transmission time T, (x, y, z) are the receiver coordinates at the signal reception time T ═ T + ρ/c,c is the speed of light, dT is the satellite clock difference, dT is the receiver clock difference, d_orbIs the satellite orbital error, d_tropIs tropospheric delay, d_ion/LiIs ionospheric delay on Li, lambda_iIs the wavelength of Li, N_iIs the integer ambiguity of Li,is the initial phase of the receiver oscillator,is the initial phase of the satellite oscillator, d_mpath/P(Li)Is the multipath effect of pseudorange measurements on Li, d_mpath/φ(Li)Multipath effects of carrier phase measurements on Li, () measurement noise;

(3) calculating non-differential combined observed values of double-frequency ionosphere-free including ionosphere-free pseudo range observed value P_IFAnd ionosphere-free phase observation value phi_IF(ii) a Wherein,

$P_{IF}=\frac{f_1^{2}P\left(L_1\right)-f_2^{2}P\left(L_2\right)}{f_1^2-f_2^2}=\mathrm{\ρ}-cdT+d_{orb}+M\·zpd+d_{mpath/P(L1+L2)}+\mathrm{\ϵ}\left(P\right(L1+L2\left)\right),$

the product of the tropospheric delay zpd and the mapping function M is the tropospheric delay d_trop；

(4) Construction System x_k＝(X_k,Y_k,Z_k,zpd_k,N_k) (ii) a Constructing a measurement equation y_k＝(P_IF,k,Φ_IF,k)＝H(x_k)x_k(ii) a Structural system estimatorProgram for programmingEstimation of equation covariance for a formation system $P_k(+)=(I-K_{k}H({\hat{x}}_k(-)\left)\right)P_k(-);$ Constructing a gain equation $K_k=P_k(-)H\left(\hat{x}\right(-\left)\right){\left(H\right({\hat{x}}_k(-)\left)P_k\right(-\left)H{\left({\hat{x}}_k(-)\right)}^T\right)}^{-1};$ Wherein, X_k,Y_k,Z_kRepresenting the receiver position at time k, zpd_kTropospheric delay at time k, N_kIs the ambiguity at time k; p_IF,kRepresenting ionospheric-free pseudorange observations at time k, phi_IF,kRepresents the ionospheric-free phase observation at time k, H (x)_k) Is a coefficient matrix of the measurement equation;representing the system estimate at time k-1,representing the system estimate at time K, K_kIs a gain matrix, y_kTo measure the observed value; i is the identity matrix, P_k(-) at time k-1Covariance of (P)_k(+) is the time kThe covariance of (a); ambiguity and coordinate accumulation are carried out by adopting extended Kalman filtering, fast convergence is carried out, and coordinate X of the receiver is obtained through calculation_k,Y_k,Z_k。