Fast and precise positioning method and system

Abstract

The present application provides a fast and precise positioning method and system. The method includes: acquiring observation data of navigation satellites and LEO augmentation satellites at a current epoch; respectively acquiring navigation telegrams of the navigation satellites and the LEO augmentation satellites, and obtaining precise orbit and clock bias; correcting errors received in the positioning process according to the acquired navigation telegrams; normalizing by taking a type of satellite navigation system as reference to obtain unified linear observation equations, and calculating observation values of positioning and velocity measurement parameters; calculating estimated values of positioning and velocity measurement parameters at the current epoch through a state equation according to the calculated observation values of positioning and velocity measurement parameters and estimated values of positioning and velocity measurement parameters at the previous epoch; generating and saving positioning and velocity measurement results at the current epoch according to the estimated values of positioning and velocity measurement parameters.

Claims

1. A fast and precise positioning method, characterized in that the method comprises following steps implemented by a user receiver: a step 1 of acquiring observation data of navigation satellites and Low Earth Orbit (LEO) augmentation satellites at a current epoch and preprocessing the observation data; a step 2 of respectively acquiring navigation telegrams of the navigation satellites and the LEO augmentation satellites, and obtaining precise orbit and clock bias of the navigation satellites and precise orbit and clock bias of the LEO augmentation satellites according to the acquired navigation telegrams of the LEO augmentation satellites; a step 3 of correcting errors received in the positioning process according to the acquired navigation telegrams by model correction and parameter estimation; a step 4 of normalizing by taking a type of satellite navigation system as reference to obtain unified linear observation equations, and calculating observation values of positioning and velocity measurement parameters; a step 5 of obtaining estimated values of positioning and velocity measurement parameters at the current epoch through a state equation according to the calculated observation values of positioning and velocity measurement parameters and estimated values of positioning and velocity measurement parameters at the previous epoch; a step 6 of generating and saving positioning and velocity measurement results at the current epoch according to the estimated values of positioning and velocity measurement parameters at the current epoch, the positioning and velocity measurement results defining a navigation map dataset; and a step 7 of displaying in real-time positioning and velocity of a moving carrier according to the navigation map dataset on a screen.

2. The positioning method as claimed in claim 1, wherein the unified linear observation equations comprise positioning observation equations, and if receiver clock bias cδ{tilde over (t)}.sub.a.sup.G corresponding to global positioning system GPS is taken as reference, then the positioning observation equations of the satellite navigation systems other than the GPS are:
ρ.sub.LC,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sub.S−d.sub.ρ.sub.LC.sub.,a.sub.G)
ϕ.sub.i,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sub.S−d.sub.ρ.sub.LC.sub.,a.sub.G)+(d.sub.ρ.sub.LC.sup.S,s−d.sub.ρ.sub.LC.sub.,a.sub.S+d.sub.ϕ.sub.LC,a.sup.S−d.sub.ϕ.sub.LC.sup.S,s)+λ.sub.LCN.sub.LC,a.sup.S,s wherein m.sub.a.sup.S,s is a mapping function corresponding to a tropospheric delay parameter in a zenith direction of a monitoring station in a GNSS system, T.sub.a is the tropospheric delay parameter in the zenith direction of the monitoring station, ρ.sub.LC,a.sup.S,s and ϕ.sub.LC,a.sup.S,s are respectively a pseudo-range and phase observation values of ionosphere-free combination, d.sub.ρ.sub.LC.sub.,a.sup.s and d.sub.ϕ.sub.LC.sub.,a.sup.s are respectively hardware delays of the pseudo-range and phase observation values of ionosphere-free combination at the receiver side, d.sub.ρ.sub.LC.sub.,a.sup.s and d.sub.ϕ.sub.LC.sub.,a.sup.s are respectively hardware delays of the pseudo-range and phase observation values of ionosphere-free combination at the satellite side; ${\lambda }_{\mathrm{LC}}=\frac{c}{f_1+f_2}$ is a wavelength of the observation values of ionosphere-free combination, c is a velocity of light in vacuum, f.sub.1 is a carrier frequency at frequency point 1, f.sub.2 is a carrier frequency at frequency point 2, N.sub.LC,a.sup.S,s is a corresponding integer ambiguity parameter, and in the equations, d.sub.ρ.sub.LC.sub.,a.sup.S−d.sub.ρ.sub.LC.sub.,a.sup.G is a difference between the hardware delays of the pseudo-range at the receiver side of the GPS and any of the satellite navigation systems other than the GPS, that is, a code bias; R.sub.a,0.sup.S,s is a distance between station and satellite calculated according to the initial coordinates of the station and satellite, l, m, and n are linearization coefficients, and respectively are $\frac{x^s-x_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $\frac{y^s-y_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $\frac{z^s-z_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $x^s,$ $y^s$ and z.sup.s are coordinates of the satellite, x.sub.a, y.sub.a and z.sub.a are the initial coordinates of the monitoring station, and Δx.sub.a, Δy.sub.a and Δz.sub.a are respectively correction values thereof.

3. The positioning method as claimed in claim 2, wherein the unified linear observation equations comprise a velocity measurement observation equation, which is: ${\stackrel{.}{\varphi }}_{i,a}^{S,s}=l\left({\stackrel{.}{x}}_a-{\stackrel{.}{x}}^S\right)+m\left({\stackrel{.}{y}}_a-{\stackrel{.}{y}}^S\right)+n\left({\stackrel{.}{z}}_a-{\stackrel{.}{z}}^S\right)+m_a^{S,s}{\stackrel{.}{T}}_a+c\delta {\stackrel{\stackrel{.}{~}}{t}}_a$ ${\stackrel{.}{\varphi }}_{i,a}^{S,s}=\frac{{\varphi }_{i,a}^{S,s}\left(k+1\right)-{\varphi }_{i,a}^{S,s}\left(k-1\right)}{2\Delta t}$ wherein {dot over (ϕ)}.sub.i,a.sup.S,s denotes the rate of phase change between the monitoring station and the satellite in the unit of cycle/s, Δt denotes the sampling interval, and {dot over (x)}.sup.s, {dot over (y)}.sup.s and ż.sup.s are rates of the satellite, {dot over (x)}.sub.a, {dot over (y)}.sub.a and ż.sub.a are the rates of the monitoring station, δ{tilde over ({dot over (t)})}.sub.a denotes a receiver clock velocity, {dot over (T)}.sub.a denotes a rate of change of the troposphere.

4. The positioning method as claimed in claim 2, wherein the navigation satellites comprise at least one of the US Global Positioning System GPS, China Beidou, EU Galileo, and Russian GLONASS satellite navigation systems.

5. A fast and precise positioning system, characterized in that the system is arranged in a user receiver and comprises: a satellite observation data receiving and processing apparatus configured for acquiring observation data of navigation satellites and Low Earth Orbit (LEO) augmentation satellites at each epoch and preprocessing the observation data; a satellite navigation telegram receiving and processing apparatus configured for respectively acquiring navigation telegrams of the navigation satellites and the LEO augmentation satellites at each epoch, and obtaining precise orbit and clock bias of the navigation satellites and precise orbit and clock bias of the LEO augmentation satellites according to the acquired navigation telegrams of the LEO augmentation satellites; a positioning error correcting apparatus configured for correcting errors received in the positioning process according to the acquired navigation telegrams by model correction and parameter estimation; a positioning and velocity measurement parameter observation value calculating apparatus configured for normalizing by taking a type of satellite navigation system as reference to obtain unified linear observation equations, and calculating observation values of positioning and velocity measurement parameters; a positioning and velocity measurement parameter estimated value calculating apparatus configured for obtaining estimated values of positioning and velocity measurement parameters at a current epoch through a state equation according to the calculated observation values of positioning and velocity measurement parameters and saved estimated values of positioning and velocity measurement parameters at the previous epoch; a positioning and velocity measurement result saving apparatus configured for generating and saving positioning and velocity measurement results at the current epoch according to the estimated values of positioning and velocity measurement parameters at the current epoch, the positioning and velocity measurement results defining a navigation map dataset; and a screen configured for displaying in real-time positioning and velocity of a moving carrier according to the navigation map dataset.

6. The positioning system as claimed in claim 5, wherein the satellite navigation telegram receiving and processing apparatus comprises a navigation satellite navigation telegram receiving and processing unit and a LEO augmentation satellite navigation telegram receiving and processing unit.

7. The positioning system as claimed in claim 5, wherein the positioning error correcting apparatus comprises a navigation satellite error correcting unit and a LEO augmentation satellite error correcting unit.

8. The positioning system as claimed in claim 5, wherein the unified linear observation equations comprise positioning observation equations, and if receiver clock bias cδ{tilde over (t)}.sub.a.sup.G corresponding to global positioning system GPS is taken as reference, then the positioning observation equations of the satellite navigation systems other than the GPS are:
ρ.sub.LC,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sub.S−d.sub.ρ.sub.LC.sub.,a.sub.G)
ϕ.sub.i,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sub.S−d.sub.ρ.sub.LC.sub.,a.sub.G)+(d.sub.ρ.sub.LC.sup.S,s−d.sub.ρ.sub.LC.sub.,a.sub.S+d.sub.ϕ.sub.LC,a.sup.S−d.sub.ϕ.sub.LC.sup.S,s)+λ.sub.LCN.sub.LC,a.sup.S,s wherein m.sub.a.sup.S,s is a mapping function corresponding to a tropospheric delay parameter in a zenith direction of a monitoring station in a GNSS system, T.sub.a is the tropospheric delay parameter in the zenith direction of the monitoring station, ρ.sub.LC,a.sup.S,s and ϕ.sub.LC,a.sup.S,s are respectively a pseudo-range and phase observation values of ionosphere-free combination, d.sub.ρ.sub.LC.sub.,a.sup.s and d.sub.ϕ.sub.LC.sub.,a.sup.s are respectively hardware delays of the pseudo-range and phase observation values of ionosphere-free combination at the receiver side, d.sub.ρ.sub.LC.sup.S,s and d.sub.ϕ.sub.LC.sup.S,s are respectively hardware delays of the pseudo-range and phase observation values of ionosphere-free combination at the satellite side; ${\lambda }_{\mathrm{LC}}=\frac{c}{f_1+f_2}$ is a wavelength of the observation values of ionosphere-free combination, c is a velocity of light in vacuum, f.sub.1 is a carrier frequency at frequency point 1, f.sub.2 is a carrier frequency at frequency point 2, N.sub.LC,a.sup.S,s is a corresponding integer ambiguity parameter, and in the equations, d.sub.ρ.sub.LC.sub.,a.sup.S−d.sub.ρ.sub.LC.sub.,a.sup.G is a difference between the hardware delays of the pseudo-range at the receiver side of the GPS and any of the satellite navigation systems other than the GPS, that is, a code bias; R.sub.a,0.sup.S,s is a distance between station and satellite calculated according to the initial coordinates of the station and satellite, l, m, and n are linearization coefficients, and are respectively $\frac{x^s-x_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $\frac{y^s-y_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $\frac{z^s-z_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$ $x^s,$ $y^s$ and z.sup.s are coordinates of the satellite, x.sub.a, y.sub.a and z.sub.a are the initial coordinates of the monitoring station, and Δx.sub.a, Δy.sub.a and Δz.sub.a are respectively correction values thereof.

9. The positioning system as claimed in claim 8, wherein the unified linear observation equations comprise a velocity measurement observation equation, which is: ${\stackrel{.}{\varphi }}_{i,a}^{S,s}=l\left({\stackrel{.}{x}}_a-{\stackrel{.}{x}}^S\right)+m\left({\stackrel{.}{y}}_a-{\stackrel{.}{y}}^S\right)+n\left({\stackrel{.}{z}}_a-{\stackrel{.}{z}}^S\right)+m_a^{S,s}{\stackrel{.}{T}}_a+c\delta {\stackrel{\stackrel{.}{~}}{t}}_a$ ${\stackrel{.}{\varphi }}_{i,a}^{S,s}=\frac{{\varphi }_{i,a}^{S,s}\left(k+1\right)-{\varphi }_{i,a}^{S,s}\left(k-1\right)}{2\Delta t}$ wherein {dot over (ϕ)}.sub.i,a.sup.S,s denotes a rate of phase change between the monitoring station and the satellite in a unit of cycle/s, Δt denotes a sampling interval, and {dot over (x)}.sup.s, {dot over (y)}.sup.s and ż.sup.s are the rates of the satellite, {dot over (x)}.sub.a, {dot over (y)}.sub.a and ż.sub.a are the rates of the monitoring station, δ{tilde over ({dot over (t)})}.sub.a denotes a receiver clock velocity, {dot over (T)}.sub.a denotes a rate of change of a troposphere.

10. A non-transitory computer-readable storage medium storing at least one executable instruction, wherein the executable instruction is configured to cause a processor to perform a fast and precise positioning method, and the method comprises following steps implemented by a user receiver: a step 1 of acquiring observation data of navigation satellites and Low Earth Orbit (LEO) augmentation satellites at a current epoch and preprocessing the observation data; a step 2 of respectively acquiring navigation telegrams of the navigation satellites and the LEO augmentation satellites, and obtaining a precise orbit and a clock bias of the navigation satellites and a precise orbit and a clock bias of the LEO augmentation satellites according to the acquired navigation telegrams of the LEO augmentation satellites; a step 3 of correcting errors received in the positioning process according to the acquired navigation telegrams by model correction and parameter estimation; a step 4 of normalizing by taking a type of satellite navigation system as reference to obtain unified linear observation equations, and calculating observation values of positioning and velocity measurement parameters; a step 5 of obtaining estimated values of positioning and velocity measurement parameters at the current epoch through a state equation according to the calculated observation values of positioning and velocity measurement parameters and estimated values of positioning and velocity measurement parameters at the previous epoch; a step 6 of generating and saving positioning and velocity measurement results at the current epoch according to the estimated values of positioning and velocity measurement parameters at the current epoch, the positioning and velocity measurement results defining a navigation map dataset; and a step 7 of displaying in real-time positioning and velocity of a moving carrier according to the navigation map dataset on a screen.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The drawings are only used to illustrate example embodiments, and are not considered as limitation to the present application. And throughout the drawings, the same reference symbols are used to denote the same components.

(2) In the drawings:

(3) FIG. 1 is a flowchart of a fast and precise positioning method according to an embodiment of the present application;

(4) FIG. 2 is a schematic structural diagram of a rapid and precise positioning system according to an embodiment of the present application;

(5) FIG. 3 is a flow chart of a fast and precise positioning method according to an embodiment of the present application.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

(6) Hereinafter, exemplary embodiments of the present disclosure will be described in more detail with reference to the accompanying drawings. Although the drawings show exemplary embodiments of the present disclosure, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein.

(7) I. Unified Linear Observation Equations of Medium-, High-, Low-Orbit Augmentation Satellites Obtained by Normalizing by Taking a Type of Satellite Navigation System as Reference

(8) To implement the positioning method provided by the present application, it is first necessary to construct and linearize unified observation equations for the medium-, high-, low-orbit augmentation satellites, and the receiver obtains the observation values of positioning and velocity measurement parameters according to the constructed linear observation equations. Wherein multi-frequency information sources of the medium-, high-, low-orbit constellation include multi-frequency information sources of at least one of all existing satellite navigation systems and the LEO augmentation satellite navigation system. The navigation satellites and the LEO augmentation satellites have the same positioning methods, and observation values of the two can be put together for adjustment solution. The mathematical model of the observation equations themself is a nonlinear equation, so it is necessary to perform Taylor expansion of the equation, and a linear equation can be obtained after discarding the second-order terms. The observation values of the navigation satellites and the LEO augmentation satellites can be expressed as a linear equation system of the positions and the receiver clock biases. Using differential observation values, the observation equations related to the monitoring station velocity term and the rate of change of the receiver clock biases can be obtained. By combining these two types of observation equations, the optimal estimation of the three parameters of PVT can be obtained.

(9) The basic observation values of the navigation satellites acquired by the receiver from navigation telegrams include two type of pseudo-ranges ρ and carrier phases ϕ at multiple frequency points. The observation values of the pseudo-range and phase from the satellite s to the monitoring station a at the frequency point i can be expressed as:
ρ.sub.i,a.sup.s=R.sub.a.sup.s+m.sub.a.sup.sT.sub.a+cδt.sub.a−cδt.sup.s+γ.sub.iI.sub.a.sup.s+d.sub.ρ.sub.t.sub.,a−d.sub.ρ.sub.t.sup.s
ϕ.sub.i,a.sup.s=R.sub.a.sup.s+m.sub.a.sup.sT.sub.a+cδt.sub.a−cδt.sup.s+γ.sub.iI.sub.a.sup.s+d.sub.ϕ.sub.t.sub.,a−d.sub.ϕ.sub.t.sup.s+λ.sub.iN.sub.i,a.sup.s  (1)
wherein R.sub.a.sup.s is the geometric distance between the satellite and the monitoring station, T.sub.a is the tropospheric delay parameter in the zenith direction of the monitoring station, the mapping function corresponding to T.sub.a is m.sub.a.sup.s, c is the velocity of light in vacuum, δt.sup.s and δt.sub.a are respectively the satellite clock bias and the receiver clock bias,

(10) ${\gamma }_i=\frac{f_1^2}{f_2^2},$
in which f.sub.i is the carrier frequency at the frequency point i, the wavelength corresponding to f.sub.i is

(11) ${\lambda }_i=\frac{c}{f_i},$
I.sub.a.sup.s is the oblique ionospheric delay, d.sub.ρ.sub.i.sub.,a and d.sub.ϕ.sub.i.sub.,a are respectively the hardware delays of the pseudo-range and the carrier phase at the receiver side, d.sub.ρ.sub.i.sup.s and d.sub.ϕ.sub.i.sup.s are the hardware delays of the pseudo-range and the carrier phase at the satellite side, N.sub.i,a.sup.s is the integer ambiguity parameter. In the above equation, error corrections such as antenna phase deviation and variation, phase Wind-Up, and relativistic effect of satellite clock bias, and errors such as multipath and observation value noises are ignored.

(12) In the equation (1), the zenith tropospheric delay parameter T.sub.a and the receiver clock bias δt.sub.a are only related to the monitoring station, the satellite clock bias δt.sup.s is only related to the satellite, the oblique ionospheric delay parameter I.sub.a.sup.s is related to the monitoring station and the satellite, and the hardware delay parameters of the pseudo-range and phase at the satellite side or the receiver side are mainly related to the monitoring station, the satellite, types of observation values, tracking frequency and the like, respectively.

(13) In navigation satellite data processing, different types of combinations of phase and pseudo-range observation values are often constructed as needed, wherein because the influence of the first-order ionosphere is eliminated in the ionosphere-free combination, it is widely used to construct observation equations for high-precision data processing. The observation equation can be expressed as:

(14) ${\rho }_{\mathrm{LC},a}^s=\frac{f_1^2}{f_1^2-f_2^2}{\rho }_{1,a}^s-\frac{f_2^2}{f_1^2-f_2^2}{\rho }_{2,a}^s=R_a^s+m_a^s{T}_a+c\delta t_a-c\delta t^s+d_{{\rho }_{\mathrm{LC}},a}-d_{{\rho }_{\mathrm{LC}}}^s$${\varphi }_{\mathrm{LC},a}^s=\frac{f_1^2}{f_1^2-f_2^2}{\varphi }_{1,a}^s-\frac{f_2^2}{f_1^2-f_2^2}{\varphi }_{2,a}^s=R_a^s+m_a^s{T}_a+c\delta t_a-c\delta t^s+d_{{\varphi }_{\mathrm{LC}},a}-d_{{\varphi }_{\mathrm{LC}}}^s+{\lambda }_{\mathrm{LC}}{N}_{\mathrm{LC},a}^s.$
wherein ρ.sub.LC,a.sup.s and ϕ.sub.LC,a.sup.s are respectively pseudo-range and phase observation values of the ionosphere-free combination, d.sub.ρ.sub.LC.sub.,a and d.sub.ϕ.sub.LC.sub.,a are respectively the hardware delays of the pseudo-range and phase observation values of the ionosphere-free combination at the receiver side, the values of which are:

(15) 0$\begin{array}{cc}{d}_{{\rho }_{\mathrm{LC}},a}=\frac{f_1^2}{f_1^2-f_2^2}{d}_{{\rho }_1,a}-\frac{f_2^2}{f_1^2-f_2^2}{d}_{{\rho }_2,a}& \left(3\right)\end{array}$$d_{{\varphi }_{\mathrm{LC}},a}=\frac{f_1^2}{f_1^2-f_2^2}{d}_{{\varphi }_1,a}-\frac{f_2^2}{f_1^2-f_2^2}{d}_{{\varphi }_2,a}$

(16) The hardware delays d.sub.ρ.sub.LC.sup.s and d.sub.ϕ.sub.LC.sup.s of the pseudo-range and phase observation values of the ionosphere-free combination at the satellite side are respectively:

(17) $\begin{array}{cc}{d}_{{\rho }_{\mathrm{LC}}}^s=\frac{f_1^2}{f_1^2-f_2^2}{d}_{{\rho }_1}^s-\frac{f_2^2}{f_1^2-f_2^2}{d}_{{\rho }_2}^s& \left(4\right)\end{array}$$d_{{\varphi }_{\mathrm{LC}}}^s=\frac{f_1^2}{f_1^2-f_2^2}{d}_{{\varphi }_1}^s-\frac{f_2^2}{f_1^2-f_2^2}{d}_{{\varphi }_2}^s$
wherein

(18) ${\lambda }_{\mathrm{LC}}=\frac{c}{f_1+f_2}$
is the wavelength of the observation values of ionosphere-free combination, N.sub.LC,a.sup.s the corresponding integer ambiguity parameter, the value of which is:

(19) $\begin{array}{cc}{N}_{\mathrm{LC},a}^s=N_{1,a}^s+\frac{f_2}{f_1-f_2}\left(N_{1,a}^s-N_{2,a}^s\right)& \left(5\right)\end{array}$

(20) Taking into account the correlation of each parameter in the equation (2) with the monitoring stations, the satellites and the signal frequencies, etc., for multi-system observations, the equation (2) can be extended to:
ρL.sub.C,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδt.sub.a−cδt.sup.S,s+d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sup.S,s
ϕL.sub.C,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδt.sub.a−cδt.sup.S,s+d.sub.ϕ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sup.S,s+λ.sub.LCN.sub.LC,a.sup.S,s
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδt.sub.a−cδt.sup.S,s+d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sup.S,s
ϕ.sub.i,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδt.sub.a−cδt.sup.S,s+d.sub.ϕ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sup.S,s+λ.sub.LCN.sub.LC,a.sup.S,s  (6)
wherein S denotes the GNSS system. For GPS, Galileo, QZSS and Beidou navigation satellite systems and the like that use code division multiple access technology, the carrier frequencies of different satellites thereof are the same, so the hardware delays of the pseudo-range and carrier phase observation values at the receiver side are the same for all single-system satellites. However, because the GLONASS system uses frequency division multiple access technology, its corresponding hardware delays of the pseudo-range and phase at the receiver side are also related to the satellite (frequency), and different GLONASS satellites (frequency) correspond to different hardware delays at the receiver side.

(21) Since in the PVT model, the clock biases of the navigation satellites are the same and they will absorb the hardware delay d.sub.ρ.sub.LC.sup.s of the pseudo-range at the satellite side during estimation, and the hardware delay of the pseudo-range d.sub.ρ.sub.LC.sub.,a at the receiver side will be absorbed by the receiver clock bias, the observation equations at this time are:
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.S
ϕ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.S+(d.sub.ρ.sub.LC.sup.S,s−dρ.sub.LC,a.sup.s+d.sub.ϕ.sub.LC,a.sup.s−dϕ.sub.LC.sup.S,s)+λ.sub.LCN.sub.LC,a.sup.S,s
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.S  (7)
wherein cδ{tilde over (t)}.sub.a=cδ{tilde over (t)}.sub.a+d.sub.ρ.sub.LC,a. When the multi-mode navigation system is processed in combination, generally only one receiver clock bias parameter is estimated, but the foregoing indicates that the receiver clock bias parameter will absorb the hardware delay of the observation value of the pseudo-range at the receiver side, and this delay parameter is related to the signal frequency and the navigation system, which thus results in that different systems correspond to different receiver clock biases cδ{tilde over (t)}.sub.a. If the receiver clock bias cδ{tilde over (t)}.sub.a.sup.G corresponding to the GPS system is taken as reference, the observation equations of other systems can be rewritten as:
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)
ϕ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G++(d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)+(d.sub.ρ.sub.LC.sup.S,s−d.sub.ρ.sub.LC.sub.,a.sup.s+d.sub.ϕ.sub.LC.sub.,a.sup.s−d.sub.ϕ.sub.LC.sub.,a.sup.S,s)+λ.sub.LCN.sub.LC,a.sup.S,s
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)  (8)
wherein d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G is the difference between the hardware delays of the pseudo-range at the receiver side of the GPS and any of the satellite navigation systems other than the GPS, that is, the code bias. If considering the difference of the time references between the different navigation systems, it is necessary to introduce an additional constant bias parameter. This constant bias parameter and the Differential Code Bias (DCB) parameter will constitute the Inter-System Bias (ISB) parameter. For navigation systems using code division multiple access, all satellites thereof correspond to the same ISB parameters, and because the GLONASS system uses frequency division multiple access technology, different satellites (frequencies) thereof correspond to different ISB parameters, and at this time the ISB parameters are actually the combination of different inter-system code biases, time reference differences and Inter-Frequency Biases (IFB) of different satellites in the GLONASS system. As a satellite navigation system, the LEO augmentation satellite constellation has the same positioning mathematical model as the related GNSS system. The LEO augmentation satellite navigation system can be regarded as a new navigation system, and only needs to estimate additional ISB parameters.
The GNSS observation equation itself is a nonlinear equation, and related parameter estimation methods are generally applicable to linear systems, so it needs to be Taylor expanded. The GNSS observation equation is expanded at the approximate coordinates of the monitoring station according to the Taylor's formula, and its second-order terms are discarded, so that linear expressions about position and time are obtained as follows:
ρ.sub.LC,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)
ϕ.sub.LC,a.sup.S,s=R.sub.a,0.sup.S,s−lΔx.sub.a−mΔy.sub.a−nΔz.sub.a+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G++(d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)+(d.sub.ρ.sub.LC.sup.S,s−d.sub.ρ.sub.LC.sub.,a.sup.s+d.sub.ϕ.sub.LC.sub.,a.sup.s−d.sub.ϕ.sub.LC.sup.S,s)+λ.sub.LCN.sub.LC,a.sup.S,s.sup.−
ρ.sub.LC,a.sup.S,s=R.sub.a.sup.S,s+m.sub.a.sup.S,sT.sub.a+cδ{tilde over (t)}.sub.a.sup.G+(d.sub.ρ.sub.LC.sub.,a.sup.s−d.sub.ρ.sub.LC.sub.,a.sup.G)   (9)
wherein R.sub.a,0.sup.S,s is the distance between station and satellite calculated according to the initial coordinates of the station and satellite, l, m and n are linearization coefficients, and are respectively

(22) $\frac{x^s-x_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$$\frac{y^s-y_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$$\frac{z^s-z_a}{\sqrt{{\left(x^s-x_a\right)}^2+{\left(y^s-y_a\right)}^2+{\left(z^s-z_a\right)}^2}},$
and x.sup.s, y.sup.s and z.sup.s are the coordinates of the satellite, x.sub.a, y.sub.a and z.sub.a are the initial coordinates of the monitoring station, Δx.sub.a, Δy.sub.a and Δz.sub.a respectively are respectively correction values thereof.

(23) In the equation (9), only the timing and positioning functions are completed, and the velocity measurement observation equation is:

(24) $\begin{array}{cc}{\stackrel{.}{\varphi }}_{i,a}^{S,s}=l\left({\stackrel{.}{x}}_a-{\stackrel{.}{x}}^S\right)+m\left({\stackrel{.}{x}}_a-{\stackrel{.}{x}}^S\right)+n\left({\stackrel{.}{x}}_a-{\stackrel{.}{x}}^S\right)+m_a^{S,s}{\stackrel{.}{T}}_a+c\delta {\stackrel{\stackrel{.}{~}}{t}}_a& \left(10\right)\end{array}$${\stackrel{.}{\varphi }}_a^{S,s}=\frac{{\stackrel{.}{\varphi }}_a^{S,s}\left(k+1\right)-{\stackrel{.}{\varphi }}_a^{S,s}\left(k-1\right)}{2\Delta t}$
wherein {dot over (ϕ)}.sub.i,s.sup.S,s; denotes the rate of phase change between the monitoring station and the satellite in the unit of cycle/s, Δt denotes the sampling interval, and is {dot over (x)}.sup.s, {dot over (y)}.sup.s and ż.sup.s are the rates of the satellite, {dot over (x)}.sub.a, {dot over (y)}.sub.a and ż.sub.a are the rates of the monitoring station, δ{tilde over ({dot over (t)})}.sub.a denotes the receiver clock velocity, {dot over (T)}.sub.a denotes the rate of change of the troposphere.
II. Constructing Positioning and Velocity Measurement Parameter State Equation with Root-Mean-Square Filtering Algorithm

(25) After establishing position and time observation equations and velocity observation equations, the root-mean-square filtering algorithm is used to perform state estimation on the positioning and velocity measurement parameters. Due to the addition of the LEO augmentation satellite observation values, the rapid convergence of the PPP can be realized and parameters information with higher precision can be obtained.

(26) The main steps of root-mean-square information filtering will be given below, and its state equation is:
x.sub.k=Φ(t.sub.k,t.sub.k-1)x.sub.k-1+Γ(t.sub.k,t.sub.k-1)u.sub.k-1
wherein x.sub.k-1 has priori value x.sub.k-1 and priori variance P.sub.k-1, and is the position, velocity or clock bias parameter to be estimated. The square root of the priori variance (Cholesky decomposition) is calculated to construct a virtual observation equation:
b.sub.k-1=R.sub.k-1x.sub.k-1
wherein x.sub.k-1=x.sub.k-1+η.sub.k-1, E[η.sub.k-1]=0, E[η.sub.k-1η.sub.k-1.sup.T]=P.sub.k-1=R.sub.k-1.sup.−1R.sub.k-1.sup.−T. And in the original equation:
ū.sub.k-1=u.sub.k-1+α.sub.k-1

(27) E[α.sub.k-1]=0, E[α.sub.k-1α.sub.k-1.sup.T]=Q, thereby constructing the virtual observation equation of state noise:
b.sub.u.sub.k−1=R.sub.uū.sub.k-1=R.sub.uu.sub.k-1+α.sub.k-1.

(28) The flitering observation equation is:
y.sub.k-1=H.sub.k-1x+ε.sub.k-1
wherein E[ε]=0, E(εε.sup.T)=I.

(29) According to the minimum variance criterion, the observation update function of the root-mean-square information filtering algorithm can be constructed:
Ĵ.sub.k-1=∥R.sub.k-1x.sub.k-1−b.sub.k-1∥.sup.2+(H.sub.k-1x.sub.k-1−y.sub.k-1).sup.2+∥R.sub.uu.sub.k-1−b.sub.u.sub.k−1∥.sup.2,

(30) If it is written in matrix form, then:

(31) ${\hat{J}}_{k-1}={.Math.\left[\begin{array}{c}{\stackrel{\_}{R}}_{k-1}\\ H_{k-1}\end{array}\right]x_{k-1}-\left[\begin{array}{c}{\stackrel{\_}{b}}_{k-1}\\ y_{k-1}\end{array}\right].Math.}^2+{.Math.R_u{u}_{k-1}-{\stackrel{\_}{b}}_{u_{k-1}}.Math.}^2,$

(32) By orthogonally changing the above equation, the following can be obtained:

(33) ${\hat{J}}_{k-1}={.Math.\left[\begin{array}{c}{\hat{R}}_{k-1}\\ 0\end{array}\right]x_{k-1}-\left[\begin{array}{c}{\hat{b}}_{k-1}\\ e_{k-1}\end{array}\right].Math.}^2+{.Math.R_u{u}_{k-1}-{\stackrel{\_}{b}}_{u_{k-1}}.Math.}^2.$

(34) It is also possible to construct the state update function of the root-mean-square information filtering algorithm according to the minimum variance criterion:
J.sub.k=(e.sub.k-1).sup.2∥{circumflex over (R)}.sub.k-1Φ.sup.−1(t.sub.k,t.sub.k-1)(x.sub.k−Γ(t.sub.k,t.sub.k-1)u.sub.k-1)−{circumflex over (b)}.sub.k-1∥.sup.2+∥R.sub.uu.sub.k-1−b.sub.u.sub.k−1∥.sup.2

(35) If it is written in matrix form, then:

(36) ${\stackrel{\_}{J}}_k={\left(e_{k-1}\right)}^2+.Math.\left[\begin{array}{cc}{R}_u& 0\\ -{\tilde{R}}_k\Gamma \left(t_k,t_{k-1}\right)& {\tilde{R}}_k\end{array}\right]\left[\begin{array}{c}{u}_{k-1}\\ x_k\end{array}\right]-\left[\begin{array}{c}{\stackrel{\_}{b}}_{u_{k-1}}\\ {\hat{b}}_{k-1}\end{array}\right].Math.$
wherein {tilde over (R)}.sub.k={circumflex over (R)}.sub.k-1Φ.sup.−1(t.sub.k, t.sub.k-1), and by orthogonally transforming, the following can be obtained:

(37) $\begin{array}{cc}{\stackrel{\_}{J}}_k={\left(e_{k-1}\right)}^2+.Math.\left[\begin{array}{cc}{\stackrel{\_}{R}}_{u_k}& {\stackrel{\_}{R}}_{{\mathrm{ux}}_k}\\ 0& {\stackrel{\_}{R}}_k\end{array}\right]\left[\begin{array}{c}{u}_{k-1}\\ x_k\end{array}\right]-\left[\begin{array}{c}{\tilde{b}}_{u_{k-1}}\\ {\stackrel{\_}{b}}_{k-1}\end{array}\right].Math..& \left(11\right)\end{array}$

(38) When using medium and high orbit information sources to solve positioning and velocity measurement parameters, due to the limitations of the satellite constellation, the accuracy of the solution and the convergence time often cannot meet the requirements of fast and high-precision positioning. The use of medium-, high-, low-orbit multi-frequency information source fusion positioning can augment the geometric structure of the visible satellites, achieve rapid convergence, and thereby improve the accuracy of positioning solution. FIG. 1 shows a flowchart of a positioning method according to an embodiment of the present application. The method may be executed by a positioning system which may be implemented in at least one of software and hardware. Wherein the medium-, high-, low-orbit constellation multi-frequency information sources include the multi-frequency information sources of at least one of all existing satellite navigation systems and the LEO augmentation satellite navigation system. As shown in FIG. 1, the positioning method according to an embodiment of the present application includes a step S110 to a step S160.

(39) In the step S110, observation data of navigation satellites and LEO augmentation satellites are acquired at a current epoch and preprocessed. The process is as follows: acquiring multi-system multi-band observation values and LEO augmentation satellite observation values through receiver tracking and observations, and preprocessing the data.

(40) Wherein the navigation satellites include at least one of the US GPS, China Beidou, EU Galileo, and Russian GLONASS satellite navigation systems.

(41) In the step S120, navigation telegrams of the navigation satellites and the LEO augmentation satellites are acquired, and precise orbits and clock biases of the navigation satellites and of the LEO augmentation satellites are obtained at the same time according to the acquired navigation telegrams of the LEO augmentation satellites. The process is: acquiring the navigation telegrams of the navigation satellites and the LEO augmentation satellites, and using the number of orbits and the clock bias coefficients provided by the navigation telegrams to interpolate to obtain the satellite position and the satellite clock bias at the current time point. Wherein because the LEO augmentation satellites have different characteristics from the navigation satellites, the navigation telegrams of the LEO augmentation satellites are different from the navigation telegrams of the navigation satellites. For example, the navigation telegrams of the LEO augmentation satellites have more types of parameters. Therefore, the calculation of the orbits and clock biases of the LEO augmentation satellites is also different from the calculation of the orbits and clock biases of the navigation satellites. For example, compared with the calculation of the orbits of the navigation satellites, in the calculation of the orbits of the LEO augmentation satellites, more perturbation factors need to be considered. Because the accuracy of the orbits and satellite clock biases of the LEO augmentation satellites and the navigation satellites given by the broadcast ephemeris generally cannot meet the high-precision positioning requirements. In an embodiment, in order to obtain real-time orbits and real-time clock biases with high precision, State Space Representation (SSR) correction information can be received in real time through the network.

(42) In the step S130, errors received in the positioning process are corrected according to the acquired navigation telegrams.

(43) Errors that can be corrected by the error model are corrected in the step 130, and then errors that cannot be corrected by the error model are corrected through the calculation of the positioning and velocity measurement parameter observation values in the step S140 and the estimation of the positioning and velocity measurement parameters in the step S150. For the LEO augmentation satellites and the navigation satellites, some errors are different, and the errors need to be corrected correspondingly according to different satellite navigation systems.

(44) In the positioning process, the positioning result is often affected by multiple terms of errors, and weakening each error is the basis for obtaining high-precision positioning results. According to the correlation, these errors can be divided into errors related to the monitoring stations, errors related to the satellites, and errors related to satellite signal propagation. The commonly used methods to weaken the positioning errors include model correction and parameter estimation. For some error terms, the physical characteristics of which has been understood, their effects can be accurately eliminated by using correction formulas, such as relativistic effects, earth rotation effects, etc.; for some error terms that can be fitted with a model, their effects can be eliminated by using model values obtained by the fitting model, such as solid earth tide correction, tropospheric correction, etc.; and for some other error terms with unknown physical characteristics and poor model fitting, parameter estimation methods can be used to eliminate their impact on positioning, such as receiver clock biases, etc.:

(45) In the step S140, observation values of positioning and velocity measurement parameters are calculated according to unified linear observation equations obtained by normalizing by taking a type of satellite navigation system as reference. The process is as follows: according to the obtained observation data and navigation telegrams, the position of the receiver through the above formula (9) is calculated, and the clock bias of the receiver can also be calculated at the same time; the velocity of the receiver can also be calculated through the above formula (10).

(46) In the step S150, according to the calculated observation values of positioning and velocity measurement parameters and estimated values of positioning and velocity measurement parameters at the previous epoch, estimating positioning and velocity measurement parameters at the present epoch through a state equation, to obtain estimated values of positioning and velocity measurement parameters at the current epoch. The process is as follows: according to the calculated observation values of positioning and velocity measurement parameters and estimated values of positioning and velocity measurement parameters at the previous epoch, calculating estimated values of positioning and velocity measurement parameters at the present epoch through the above formula (11), and saving the calculated estimated values of positioning and velocity measurement parameters.

(47) In the step S160, according to the estimated values of positioning and velocity measurement parameters at the current epoch, positioning and velocity measurement results at the current epoch are generated and saved, and the processing returns to the step S110.

(48) FIG. 2 shows a positioning system according to an embodiment of the present application. Wherein the medium-, high-, low-orbit constellation multi-frequency information sources include the multi-frequency information sources of at least one of all existing satellite navigation systems and the LEO augmentation satellite navigation system.

(49) As shown in FIG. 2, the positioning system according to an embodiment of the present application includes: a satellite observation data receiving and processing apparatus 11, a satellite navigation telegram receiving and processing apparatus 12, a positioning error correcting apparatus 13, a positioning and velocity measurement parameter observation value calculating apparatus 14, a positioning and velocity measurement parameter estimated value calculating apparatus 15 and a positioning and velocity measurement result saving apparatus 16.

(50) Wherein the satellite observation data receiving and processing apparatus 11 is configured for acquiring observation data of navigation satellites and LEO augmentation satellites at each epoch and preprocessing the data.

(51) The satellite navigation telegram receiving and processing apparatus 12 is configured for acquiring navigation telegrams of the navigation satellites and the LEO augmentation satellites at each epoch, and at the same time obtaining precise orbit and clock bias of the navigation satellites and of the LEO augmentation satellites according to the acquired navigation telegrams of the LEO augmentation satellites. In an embodiment, the satellite navigation telegram receiving and processing apparatus 12 includes a navigation satellite navigation telegram receiving and processing unit and a LEO augmentation satellite navigation telegram receiving and processing unit.

(52) The positioning error correcting apparatus 13 is configured for correcting errors received in the positioning process according to the acquired navigation telegrams. In an embodiment, the positioning error correcting apparatus 13 includes a navigation satellite error correcting unit and a LEO augmentation satellite error correcting unit.

(53) The positioning and velocity measurement parameter observation value calculating apparatus 14 is configured for normalizing by taking a type of satellite navigation system as reference to obtain unified linear observation equations, and calculating observation values of positioning and velocity measurement parameters according the unified linear observation equations.

(54) The positioning and velocity measurement parameter estimated value calculating apparatus 15 is configured for estimating positioning and velocity measurement parameters at the present epoch according to the calculated observation values of positioning and velocity measurement parameters and the saved estimated values of positioning and velocity measurement parameters at the previous epoch, to obtain estimated values of positioning and velocity measurement parameters at the current epoch through a state equation.

(55) The positioning and velocity measurement result saving apparatus 16 configured for generating and saving positioning and velocity measurement results at the current epoch according to the estimated values of positioning and velocity measurement parameters at the current epoch.

(56) FIG. 3 is a schematic diagram of the working principle of a fast and precise positioning method provided by an embodiment of the present application. Wherein the implementation of the positioning method may include a navigation satellite constellation, a LEO constellations, a ground operation and control system, and a user receiver having a screen capable of displaying data. Wherein: the navigation satellite constellation includes at least one of the US GPS, China Beidou, EU Galileo and Russian GLONASS satellite navigation systems, and is configured for broadcasting navigation satellite signals. The LEO constellation includes multiple LEO satellites distributed on multiple orbital planes, and the multiple LEO satellites broadcast navigation direct signals based on high-precision time-frequency references through specific frequency bands to provide stable coverage of the world or specific service areas, and is configured for broadcasting navigation direct signals and navigation augmentation information. The ground operation and control system performs business calculation and processing, and controls and manages the satellites and the constellations. The user receiver receives navigation direct signals broadcasted by the navigation satellites and the LEO satellites, as well as the navigation augmentation information broadcasted by the LEO satellites, and performs precise positioning, velocity measurement and timing based on the navigation direct signals from the navigation satellites and the LEO satellites and the navigation augmentation information.