A method for solving GNSS non-differential and non-combined observation equations with ambiguity weighted constraints

Abstract

The present invention proposes a method for solving GNSS non-differenced non-combined observation equations with ambiguity weighted constraints, which belongs to the field of satellite navigation technology. The method is based on multi-frequency pseudorange and phase observation data of global or regional GNSS tracking station networks, and constructs a full-rank non-differenced non-combined full-parameter estimation model by adding constraints, thereby realizing unified estimation of parameters such as receiver clock error, satellite clock error, tropospheric delay, ionospheric delay and base frequency satellite phase deviation. By assigning weights to parameters including ambiguity in the constraints, the method can flexibly select and adjust the benchmarks for all parameter estimates except tropospheric parameters, effectively solving the problem that the original non-differenced non-combined observation equations have a single benchmark selection and cannot be weighted and constrained when solving ambiguity parameters.

Description

Technical Field

[0001] The invention relates to a GNSS non-difference non-combined observation equation solving method with ambiguity weighted constraints, belonging to the technical field of satellite navigation. Background Art

[0002] The Global Navigation Satellite System (GNSS) can provide real-time, all-weather and global positioning, navigation and timing (PNT) services for land, sea and air. In recent years, PNT services have made great progress in real-time and high-precision, and high-precision real-time positioning technologies have become increasingly abundant, such as real-time differential kinematic positioning (RTK), real-time precise point positioning (Real-Time PPP and Real-Time PPP-AR), real-time kinematic positioning based on reference station network (network RTK), and PPP-RTK technology. The realization of high-precision GNSS services depends on novel algorithms, but it cannot be separated from the support of enhanced information, such as satellite orbits, clock errors, phase deviations, and atmospheric delays. Nowadays, many scientific research and production institutions use the data of GNSS observation station networks to solve the corresponding enhanced information to ensure high-performance GNSS applications.

[0003] However, when solving these enhanced information using undifferenced and non-combined observation equations, existing methods mainly rely on selecting the ambiguity corresponding to a single receiver and a single satellite as the basis for ambiguity estimation, and fail to use all ambiguities as the basis for ambiguity estimation. Summary of the invention

[0004] In view of this, the present invention provides a method for solving GNSS non-differenced non-combined observation equations with ambiguity weighted constraints. The method utilizes the observation data of GNSS global or regional tracking stations and realizes parameter solution by applying weighted constraints. The method can be used to solve the satellite clock error, satellite phase deviation and ionosphere, troposphere and other enhanced information required for precise positioning services such as PPP-AR and PPP-RTK.

[0005] In order to achieve the above object, the technical solution adopted by the present invention is:

[0006] A method for solving GNSS non-differenced non-combined observation equations with ambiguity weighting constraints, comprising the following steps:

[0007] (1) Using a single GNSS dual-frequency or multi-frequency raw observation data, construct the non-differenced and non-combined observation equation:

[0008]

[0009] Where, the superscript s represents the satellite number, the subscript r represents the receiver number, the subscript f represents the frequency band number, e represents the epoch number, P represents the pseudorange observation value in meters, and L represents the carrier phase observation value in meters. is the geometric distance from the receiver r to the satellite s, which is a known value in meters. is the model value of the tropospheric projection function from receiver r to satellite s, T r is the tropospheric delay in the zenith direction, in meters, dt r is the clock error of the receiver r, in meters, dt s is the clock error of satellite s, in meters, b r,f is the pseudorange hardware delay of receiver r corresponding to frequency f, in meters, is the pseudorange hardware delay of satellite s corresponding to frequency f, in meters, φ r,f is the phase hardware delay of the receiver r corresponding to the frequency f, in meters, is the phase hardware delay of satellite s corresponding to frequency f, in meters, is the ionospheric delay of the first frequency of the receiver r and satellite s, in meters, is the ratio of the ionospheric delay of other frequencies to the first frequency, λ f is the wavelength corresponding to the frequency f, in meters, is the integer ambiguity corresponding to frequency f between receiver r and satellite s, in meters;

[0010] (2) Add constraint type 1:

[0011]

[0012] Where, Δdt i is the receiver clock correction value of the i-th station, Δdt i is the correction value of the i-th satellite clock error, Pdt i is the weight of the receiver clock correction value of the i-th station, Pdt i is the weight of the i-th satellite clock correction value, nr is the number of stations, and ns is the number of satellites;

[0013] Add constraint type 2: Select the first and second frequencies, and add constraint b for each station and satellite respectively r,f =0 and Where f takes two selected frequencies;

[0014] Add constraint type 3: Add constraints for other frequencies besides the frequency selected in constraint type 2.

[0015] Where f is the frequency other than the frequency selected in constraint condition type 2, Pbi,f is the weight of the receiver pseudorange hardware delay of the i-th station, is the weight of the hardware delay of the ith satellite pseudorange, nr is the number of stations, and ns is the number of satellites;

[0016] Add constraint type 4: Add conditions for each frequency

[0017] Among them, Pφ i,f is the weight of the receiver phase hardware delay of the i-th station, is the weight of the hardware delay of the i-th satellite phase, nr is the number of stations, and ns is the number of satellites;

[0018] Add constraint type 5: Add conditions for each frequency at each station

[0019] Among them, ns is the number of satellites, is the weight of the fth frequency of the ith satellite;

[0020] Add constraint type 6: Add conditions for each frequency of each satellite

[0021] Among them, nr is the number of measuring stations, R i,f is the weight of the fth frequency at the ith station;

[0022] (3) Combine the observation equation and the six types of constraints to solve the parameter T r , Except T r Except for the unchanged, the specific forms of other parameter solutions are as follows:

[0023]

[0024] in, B r =b r,2 -b r,1 , α=1 / (μ2-1), β=μ2 / (μ2-1), γ f =α·(μ f +1), κ f =α·(μ f +μ2),

[0025]

[0026]

[0027]

[0028]

[0029]

[0030]

[0031]

[0032] Compared with the existing methods, the present invention has the following advantages:

[0033] (1) Epoch parameter T r , and They are independent of each other, the number of each type of constraint is clear, and the observation equation can be solved in a single epoch;

[0034] (2) Hardware delay Benchmark It can contain a combination of multiple stations and multiple satellites, and the weight of each station and satellite in the benchmark can be set by setting the Pb of each station or satellite. i,f and Enable flexible selection of benchmarks;

[0035] (3) The first type of ambiguity criterion and It can include a combination of multiple stations and multiple satellites, and can be adjusted by adjusting the weight Pφ r,f , and R i,f , to achieve flexible selection of benchmarks;

[0036] (4) The second type of ambiguity criterion It can include a combination of multiple stations and multiple satellites, and can be adjusted by adjusting the weight Pφ r,f , and R i,f , enabling flexible selection of benchmarks.

[0037] In summary, the method of the present invention utilizes a non-differenced and non-combined observation model to achieve synchronous estimation of all parameters to be estimated, and can flexibly select a benchmark for parameter estimation by adjusting the weights between parameters in each type of constraint conditions, thereby ensuring that parameters of the same type have good consistency, self-consistency and adjustable weights. DETAILED DESCRIPTION

[0038] A method for solving GNSS non-differenced non-combined observation equations with ambiguity weighting constraints, comprising the following steps:

[0039] (1) Observation data from globally or regionally distributed tracking stations of a single GNSS system are selected, including dual-frequency or higher multi-frequency pseudorange and phase observations.

[0040] (2) Use observation data to construct observation equations:

[0041]

[0042] (3) Add constraint type 1: Where Δdt i is the receiver clock correction value of the i-th station, Δdt i is the correction value of the i-th satellite clock error, Pdt i is the weight of the receiver clock correction value of the i-th station, Pdt i is the weight of the i-th satellite clock correction value, nr is the number of stations, ns is the number of satellites; the condition number is increased by 1.

[0043] (4) Add constraint type 2: Select the first and second frequencies and add constraint b to each station and satellite. r,f =0 and Where f takes two selected frequencies; the condition number is increased by (nr+ns)×2.

[0044] (5) Add constraint type 3: Add constraints for other frequencies other than the frequency selected in constraint 2. Where f is a frequency other than the frequency selected in constraint 2, Pb i,f is the weight of the receiver pseudorange hardware delay of the fth frequency at the ith station, is the pseudorange hardware delay weight of the fth frequency of the ith satellite, nr is the number of stations, and ns is the number of satellites; nf-2 constraints are added, where nf is the total number of frequencies.

[0045] (6) Add constraint type 4: Add conditions for each frequency where Pφ i,f is the weight of the receiver phase hardware delay at the fth frequency of the ith station, is the weight of the hardware delay of the phase of the fth frequency of the ith satellite, nr is the number of stations, and ns is the number of satellites; nf constraints are added, where nf is the total number of frequencies.

[0046] (7) Add constraint type 5: Add conditions for each frequency of each station. Where ns is the number of satellites, is the weight of the fth frequency of the ith satellite; the number of constraints is increased by nr×nf, where nf is the total number of frequencies.

[0047] (8) Add constraint type 6: Add conditions for each frequency of each satellite Where nr is the number of stations, R i,f is the weight of the fth frequency at the ith station; the number of constraints is increased by ns×nf, where nf is the total number of frequencies.

[0048] (9) By combining the observation equation and the six types of constraints, the column rank deficiency problem in the equation is eliminated, and the various parameters to be estimated are solved.

[0049] This method is based on the multi-frequency pseudorange and phase observation data of the global or regional GNSS tracking station network. By adding constraints, a full-rank non-differenced non-combined full parameter estimation model is constructed to achieve unified estimation of parameters such as receiver clock error, satellite clock error, tropospheric delay, ionospheric delay, and base frequency satellite phase deviation. By assigning weights to the parameters including ambiguity in the constraints, this method can flexibly select and adjust the benchmarks for all parameter estimates except tropospheric parameters, effectively solving the problem of single benchmark selection and inability to weight ambiguity parameters when solving the original non-differenced non-combined observation equation.

[0050] Here is a more specific example:

[0051] A method for solving GNSS non-differenced non-combined observation equations with ambiguity weighting constraints, comprising the following steps:

[0052] (1) Pseudorange and phase observation data from a single GNSS system global or regional tracking station network are selected, including BeiDou-2 triple-frequency, BeiDou-3 triple-frequency, GPS dual-frequency or triple-frequency, GLONASS dual-frequency, and Galileo dual-frequency, triple-frequency, or quad-frequency observation data.

[0053] (2) Using the pseudorange and phase observations of all station single systems and all frequencies in the tracking station network, the non-differenced and non-combined original observation equations are constructed. The known satellite positions and station coordinates are introduced into the observation equations, that is, the geometric distance from the satellite to the receiver is known, and the parameters to be estimated include: tropospheric delay, receiver clock error, satellite clock error, receiver pseudorange hardware delay, satellite pseudorange hardware delay, receiver phase hardware delay, satellite phase hardware delay, ionospheric delay, and integer ambiguity.

[0054] (3) The original observation equation for constructing a single system with multiple frequencies is:

[0055] The superscript s represents the satellite number, the subscript r represents the receiver number, the subscript f represents the frequency band number, the e represents the epoch number, P represents the pseudorange observation value in meters, and L represents the carrier phase observation value in meters. is the geometric distance from the receiver r to the satellite s, which is a known value in meters. is the model value of the tropospheric projection function from receiver r to satellite s, T ris the tropospheric delay in the zenith direction, in meters, dt r is the clock error of the receiver r, in meters, dt s is the clock error of satellite s, in meters, b r,f is the pseudorange hardware delay of receiver r corresponding to frequency f, in meters, is the pseudorange hardware delay of satellite s corresponding to frequency f, in meters, φ r,f is the phase hardware delay of the receiver r corresponding to the frequency f, in meters, is the phase hardware delay of satellite s corresponding to frequency f, in meters, is the ionospheric delay of the first frequency of the receiver r and satellite s, in meters, is the ratio of the ionospheric delay of other frequencies to the first frequency, λ f is the wavelength corresponding to the frequency f, in meters, is the integer ambiguity corresponding to frequency f between receiver r and satellite s, in meters.

[0056] (4) Add constraint type 1: Where Δdt i is the receiver clock correction value of the i-th station, Δdt i is the correction value of the i-th satellite clock error, Pdt i is the weight of the receiver clock correction value of the i-th station, Pdt i is the weight of the i-th satellite clock correction value, nr is the number of stations, ns is the number of satellites; the condition number is increased by 1.

[0057] (5) Add constraint type 2: Select the first and second frequencies and add constraint b for each station and satellite. r,f =0 and Where f takes two selected frequencies; the condition number is increased by (nr+ns)×2.

[0058] (6) Add constraint type 3: Add constraints for other frequencies other than the frequency selected in constraint 2. Where f is a frequency other than the frequency selected in constraint 2, Pb i,f is the weight of the receiver pseudorange hardware delay of the fth frequency at the ith station, is the pseudorange hardware delay weight of the fth frequency of the ith satellite, nr is the number of stations, and ns is the number of satellites; nf-2 constraints are added, where nf is the total number of frequencies.

[0059] (7) Add constraint type 4: Add conditions for each frequency where Pφ i,f is the weight of the receiver phase hardware delay at the fth frequency of the ith station, is the weight of the hardware delay of the phase of the fth frequency of the ith satellite, nr is the number of stations, and ns is the number of satellites; nf constraints are added, where nf is the total number of frequencies.

[0060] (8) Add constraint type 5: Add conditions for each frequency of each station. Where ns is the number of satellites, is the weight of the fth frequency of the ith satellite; the number of constraints is increased by nr×nf, where nf is the total number of frequencies.

[0061] (9) Add constraint type 6: Add conditions for each frequency of each satellite Where nr is the number of stations, R i,f is the weight of the fth frequency at the ith station; the number of constraints is increased by ns×nf, where nf is the total number of frequencies.

[0062] (10) Combining the observation equations and the six types of constraints, the column rank deficiency problem in the equations is eliminated, except that the zenith tropospheric delay is still T r In addition, the specific forms of other parameters obtained are as follows:

[0063]

[0064] in, B r =b r,2 -b r,1 , β=μ2 / (μ2-1),γ f =α·(μ f +1), κ f =α·(μ f +μ2),

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071] In summary, the present invention not only solves the model rank deficiency problem of the observation equation by adding weighted constraints to the clock error, pseudorange hardware delay and phase deviation, and performs parameter weighting on each type of constraint including ambiguity, but also realizes the real-time optimal estimation of multiple state domain parameters including satellite clock error, phase deviation, pseudorange deviation, atmospheric ionosphere and troposphere simultaneously. The present invention can flexibly select the benchmark for parameter estimation by adjusting the weights between parameters in each type of constraint (including ambiguity constraint), can select all ambiguity parameters as the benchmark for ambiguity estimation, can solve the observation equations of dual frequency and above, generate pseudorange deviation and phase deviation of the base frequency, and generate satellite clock error and receiver clock error corresponding to the ionospheric elimination combination. The present invention not only solves the problem of multiple state domain parameters such as satellite phase deviation, ionosphere and troposphere required by PPP-AR and PPP-RTK technology, but also realizes non-equal weight constraints between ambiguities.

Claims

1. A method for solving GNSS non-differenced non-combined observation equations with ambiguity weighting constraints, characterized in that: The following steps are involved: (1) Using a single GNSS dual-frequency or multi-frequency raw observation data, construct the non-differenced and non-combined observation equation: Where, the superscript s represents the satellite number, the subscript r represents the receiver number, the subscript f represents the frequency band number, e represents the epoch number, P represents the pseudorange observation value in meters, and L represents the carrier phase observation value in meters. is the geometric distance from the receiver r to the satellite s, which is a known value in meters. is the model value of the tropospheric projection function from receiver r to satellite s, T r is the tropospheric delay in the zenith direction, in meters, dt r is the clock error of the receiver r, in meters, dt s is the clock error of satellite s, in meters, b r,f is the pseudorange hardware delay of receiver r corresponding to frequency f, in meters, is the pseudorange hardware delay of satellite s corresponding to frequency f, in meters, φ r,f is the phase hardware delay of the receiver r corresponding to the frequency f, in meters, is the phase hardware delay of satellite s corresponding to frequency f, in meters, is the ionospheric delay of the first frequency of the receiver r and satellite s, in meters, is the ratio of the ionospheric delay of other frequencies to the first frequency, λ f is the wavelength corresponding to the frequency f, in meters, is the integer ambiguity corresponding to frequency f between receiver r and satellite s, in meters; (2) Add constraint type 1: Where, Δdt i is the receiver clock correction value of the i-th station, Δdt i is the correction value of the i-th satellite clock error, Pdt i is the weight of the receiver clock correction value of the i-th station, Pdt i is the weight of the i-th satellite clock correction value, nr is the number of stations, and ns is the number of satellites; Add constraint type 2: Select the first and second frequencies, and add constraint b for each station and satellite respectively r,f =0 and Where f takes two selected frequencies; Add constraint type 3: Add constraints for other frequencies besides the frequency selected in constraint type 2. Where f is the frequency other than the frequency selected in constraint condition type 2, Pb i,f is the weight of the receiver pseudorange hardware delay of the i-th station, is the weight of the hardware delay of the ith satellite pseudorange, nr is the number of stations, and ns is the number of satellites; Add constraint type 4: Add conditions for each frequency Among them, Pφ i,f is the weight of the receiver phase hardware delay of the i-th station, is the weight of the hardware delay of the i-th satellite phase, nr is the number of stations, and ns is the number of satellites; Add constraint type 5: Add conditions for each frequency at each station Among them, ns is the number of satellites, is the weight of the fth frequency of the ith satellite; Add constraint type 6: Add conditions for each frequency of each satellite Among them, nr is the number of measuring stations, R i,f is the weight of the fth frequency at the ith station; (3) Combine the observation equation and the six types of constraints to solve the parameter T r , Except T r Except for the unchanged, the specific forms of other parameter solutions are as follows: Among them, B r =b r,2 -b r,1 , α=1 / (μ2-1), β=μ2 / (μ2-1), γ f =a·(m f +1), k f =a·(m f +μ2),

Images