CN115561793A - Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation - Google Patents
Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation Download PDFInfo
- Publication number
- CN115561793A CN115561793A CN202211076810.6A CN202211076810A CN115561793A CN 115561793 A CN115561793 A CN 115561793A CN 202211076810 A CN202211076810 A CN 202211076810A CN 115561793 A CN115561793 A CN 115561793A
- Authority
- CN
- China
- Prior art keywords
- phase
- deviation
- time
- ambiguity
- value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 27
- 230000001360 synchronised effect Effects 0.000 claims abstract description 10
- 239000011159 matrix material Substances 0.000 claims description 30
- 238000001914 filtration Methods 0.000 claims description 26
- 238000007667 floating Methods 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 claims description 17
- 150000001875 compounds Chemical class 0.000 claims description 12
- 238000012937 correction Methods 0.000 claims description 8
- 230000002159 abnormal effect Effects 0.000 claims description 6
- 238000013461 design Methods 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 238000007781 pre-processing Methods 0.000 claims description 5
- 239000005433 ionosphere Substances 0.000 claims description 4
- 238000005295 random walk Methods 0.000 claims description 4
- 230000001174 ascending effect Effects 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 238000012544 monitoring process Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 230000003595 spectral effect Effects 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 239000005436 troposphere Substances 0.000 claims description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- 238000004804 winding Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000003908 quality control method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
- G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
- G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
- Radio Relay Systems (AREA)
Abstract
The invention discloses a real-time Beidou phase decimal deviation rapid estimation method based on parallel computation, which comprises the following steps of: step 1, data stream time synchronization and validity check; step 2, updating external required files according to the synchronous real-time observed values obtained in the step 1, and calculating the MW average value and the ionosphere-free combined floating-point ambiguity in parallel by multiple observation stations; step 3, taking the MW average value of the multiple measuring stations obtained in the step 2 as an input observation value, and quickly estimating the phase decimal deviation of the time-varying wide lane; step 4, taking the ionosphere-free combined floating-point ambiguity of the multiple measuring stations obtained in the step 2 as an input observation value, and performing rapid estimation on the time-varying narrow lane phase decimal deviation; and 5, linearly converting the decimal deviation of each frequency phase. The method and the device realize the rapid estimation of the Beidou time-varying phase decimal deviation solved by the multi-station network, improve the efficiency of the real-time Beidou phase decimal deviation estimation, and improve the precision and the reliability of the estimated Beidou phase decimal deviation.
Description
Technical Field
The invention relates to the technical field of satellite navigation and positioning, in particular to a real-time Beidou phase decimal deviation fast estimation method based on parallel computation.
Background
The Beidou Satellite Navigation System (BDS) is a Global Navigation Satellite System (GNSS) which is independently constructed and operated in China and is set as an important national space foundation. At present, the annual production value of the industry based on Beidou position service reaches four billion yuan, a stable high-speed growth situation is presented, and unmanned driving, smart city construction, natural resource investigation and the like all need reliable spatial position information.
Precision Point-to-Point Positioning (PPP) technology [1] Accurate three-dimensional coordinates under a geocentric reference frame can be acquired in a global range, however, as the ambiguity parameter absorbs pseudo range and phase hardware delay, the precision and reliability of PPP floating solution losing the characteristics of the whole cycle need to be further improved.
In order to realize reliable fixation of PPP ambiguity, scholars at home and abroad successively provide a phase decimal deviation model [2] Integer phase Zhong Moxing [3] Clock difference decoupling model [4] And the three methods are compared in the aspects of server correction number type, client model freedom degree and the like. The BDS system has a short construction period, and three types of satellites with mixed orbits are adopted, including geostationary orbit satellites (GEO), inclined geosynchronous orbit satellites (IGSO) and medium earth orbit satellites (MEO), so that the satellite observation value quality, the satellite attitude model, the satellite space environment and the like have certain differences with other GNSS systems, and the BDS PPP positioning performance is reduced in spite of the difference between the two.
Meanwhile, in order to meet the timeliness and long-term stability of the real-time location service, the efficiency and numerical stability of BDS time-varying phase decimal deviation estimation are important. In recent years, many scholars adopt methods such as introducing an external matrix library, a non-strict network solution model, reducing the updating time of phase fractional deviation and the like to improve the estimation efficiency of the phase fractional deviation. However, the above technical means depend on a hardware platform, cross-platform operation is difficult to achieve, accuracy and reliability of phase fractional deviation are reduced, reliable and fixed BDS PPP ambiguity is not facilitated, and a strict network solution model has the problems of time consumption of high-dimensional matrix inversion and variance-covariance and the like.
Therefore, the method can reliably estimate the real-time BDS phase decimal deviation to realize reliable BDS PPP fixed solution, and has important practical value and research significance for promoting China to provide high-quality navigation, positioning and time service for Beidou.
Disclosure of Invention
The invention aims to solve the technical problem of providing a real-time Beidou phase decimal deviation fast estimation method based on parallel computation, which can realize fast and steady estimation of real-time BDS phase decimal deviation.
In order to solve the technical problem, the invention provides a real-time Beidou phase decimal deviation fast estimation method based on parallel computation, which comprises the following steps of:
and 5, linearly converting the decimal deviation of each frequency phase.
Preferably, in step 1, the data stream time synchronization and validity check specifically includes the following steps:
step 11, acquiring local time, setting phase decimal deviation update rate and data stream waiting time, and taking the phase decimal deviation update rate and the data stream waiting time as time references for acquiring synchronous real-time data streams from epoch to epoch;
step 12, recovering the orbit and clock error data of the real-time satellite by using the broadcast ephemeris and the real-time orbit/clock error correction number, and monitoring whether the data is missing or abnormal;
step 13, detecting and updating a multi-system differential code deviation DCB file, a forecast earth orientation parameter EOP file, an antenna file and a sea tide file in a BLQ format, and taking the files as external input data of subsequent data preprocessing;
and step 14, setting waiting time aiming at the problem of observation value loss of part of observation stations, particularly the phenomenon of observation value data interruption of a plurality of continuous epochs, wherein the observation station data does not participate in the subsequent phase decimal deviation calculation in the waiting time.
Preferably, in step 2, the parallel calculation of the floating ambiguity by the multiple measurement stations specifically includes the following steps:
step 21, deleting the dual-frequency pseudo range and phase observation value missing satellites, and forming the following dual-frequency linear combination observation values one by one:
in the formula (f) 1 And f 2 Representing the first and second frequencies, P, of the BDS, respectively 1 /P 2 And L 1 /L 2 The first frequency pseudo-range observed value and the second frequency pseudo-range observed value and the phase observed value of the BDS are respectively, PI is the difference value of the pseudo-range observed values of the first frequency and the second frequency, LI is the difference value of the phase observed values of the first frequency and the second frequency, and the phase L of the dual-frequency carrier wave is 1 /L 2 And pseudorange observations P 1 /P 2 Constructing a MW combined observation value, wherein PC is a dual-frequency ionosphere-free pseudo-range combined observation value, and LC is a dual-frequency ionosphere-free phase combined observation value;
step 22, regarding the orbit height and PI value threshold of the GEO/IGSO/MEO three orbit type satellites to carry out pseudo-range gross error detection, detecting the cycle slip by using MW combination and LI combination, considering the nearly geostationary property of the GEO satellites, setting the ionosphere change rate threshold of the LI linear combination detection cycle slip to be half of the MEO/IGSO, removing the observation data exceeding the threshold, and recording the number of arc sections of each continuous observation arc section and the switching time of new and old arc sections;
step 23, correcting various error models, including tropospheric delay, antenna phase center correction, phase winding, tidal effect, and satellite multipath space error, wherein the GEO satellite attitude adopts a zero bias mode, the IGSO/MEO satellite adopts a dynamic bias mode,and calculating the average value m of the epoch-by-epoch MW by adopting a moving sliding window mw And wherein the error σ mw ：
Wherein k represents the window size within a continuous arc segment, and the maximum moving window size is set to 120;
step 24, setting prior variances of observed values of the GEO/IGSO/MEO three orbit type satellites according to a formula (3), and setting variance empirical values sigma of pseudo ranges and phase observed values of the three types of satellites respectively due to different environments of the receiver and the different orbit satellites 0 3/0.3/0.3m and 0.03/0.003/0.003m, wherein e represents the satellite height angle;
step 25, using the PC and LC corrected by the errors of the models in the step 23, adopting the ionosphere-free combined PPP model in the formula (4), using the IGG-III weight function model, judging and reducing the weight of the abnormal observed value according to the standardized residual error, thereby obtaining the reliable ionosphere-free combined floating ambiguity A if ；
In the formula (dt) r Indicating the receiver clock error, Z, of the station r r And m represents the wet delay and wet delay projection functions of the zenith troposphere of the survey station, epsilon PC And ε LC Respectively representing the unmodeled errors of the PC and LC observed values;
and step 26, when there are observation data of a plurality of observation stations, according to the steps from step 21 to step 25, adopting OpenMP parallel processing to obtain the synchronized multi-station MW average value and the ionosphere-free combined floating ambiguity.
Preferably, in step 3, the fast estimation of the time-varying wide lane phase fractional deviation specifically includes the following steps:
step 31, taking the average value of the MW of the multiple observation stations obtained in step 2 as an input observation value, and expressing the wide lane floating ambiguity as follows:
in the formula (I), the compound is shown in the specification,
denotes the mean value of MW, lambda, of the multiple stations w Indicating a wide-lane wavelength, B w,* Indicating the receiver wide-lane phase fractional deviation,
indicating satellite wide-lane phase fractional deviation, N w Is the corresponding whole-circumference wide lane ambiguity;
assuming that q receivers and p satellites are observed, and taking the whole-cycle ambiguity of the wide lane and the phase decimal deviation of the satellites and the wide lane of the receivers as estimation parameters, constructing a network solution wide lane phase decimal deviation function model as follows:
in which I denotes the identity matrix, the subscript denotes its dimension, e p A column vector with elements of 1, the subscripts indicating the number of rows,
a Kalman filtering design matrix representing network solution wide lane phase decimal deviation;
step 32, the ambiguity of the continuous arc section is a constant, the process noise should be 0, therefore the ambiguity parameter of the continuous arc section is set as a constant model, the small deviation of the wide lane phase of the receiver and the small deviation of the wide lane of the satellite receiver are respectively set as random walk due to the time-varying characteristic of the small deviation parameter of the receiver and the satellite, and the spectral density is respectively 1.0e -4 And 1.0e -6 m 2 The method comprises the following steps that (1) the scalar calculation of Kalman filtering time updating is carried out by utilizing the diagonal characteristic of a state transition matrix, namely, the state prediction of the current epoch is carried out on the basis of the parameter estimation of the last epoch, so that the rapid numerical calculation is realized;
step 33, mean error σ of smoothed value according to MW mw Performing ascending order arrangement, forming an undirected connected graph by using the receiver, the satellite and the wide lane floating ambiguity, and introducing a Kruskal minimum spanning tree algorithm to generate a steady independent wide lane ambiguity standard;
step 34, taking the wide lane ambiguity reference value as a virtual observation value, and setting the virtual variance of the wide lane ambiguity reference value as 1.0e -9 The satellite phase fractional deviation is set to 0, i.e.
Adding the constraint to a wide lane phase decimal deviation net solution model, and estimating wide lane phase decimal deviation of each satellite relative to the reference constraint by using Kalman filtering;
step 35, regarding the MW smooth value of each satellite of each station and the selected wide lane ambiguity reference value as a single observation value, and realizing Kalman filtering scalar calculation according to a formula (7), thereby avoiding the time-consuming problem of high-dimensional matrix inversion of standard filtering; the filtering variance-covariance matrix is rapidly updated by adopting OpenMP parallel computation, and the computation efficiency of network solution phase decimal deviation is further accelerated;
in the formula (I), the compound is shown in the specification,
is the variance of the ith satellite observation for the kth epoch,
and
respectively representThe Kalman filter state prediction value and the filtered value,
and
are variance matrixes corresponding to states, I represents an identity matrix, k k,i Representing a Kalman filter gain matrix,/ k,i An ith satellite observation, a, representing the kth epoch k,i Is as in formula (6)
A row vector corresponding to the observation value of the ith satellite of the observation station is arrayed;
step 36, based on the widelane ambiguity estimation value deviation of 0.15 cycle, the standard deviation of 0.15 cycle and the probability 1000 of ambiguity fixed as an integer, if the widelane ambiguity obtained by filtering a certain satellite according to the step 35 meets the requirements, fixing the widelane ambiguity as an integer, traversing all the estimated widelane ambiguities until no new widelane ambiguity can be fixed, and outputting the phase decimal deviation of the satellite widelane at the moment
The fixed widelane ambiguity obtained in this step is a prerequisite for step 4.
Preferably, in the step 4, the fast estimation of the time-varying narrow lane phase fractional deviation specifically includes the following steps:
step 41, taking the whole-cycle ambiguity of the wide lane fixed in step 3 as a known value, and correcting the non-ionosphere floating ambiguity obtained in step 2 one by one satellite, as shown in formula (8):
in the formula (I), the compound is shown in the specification,
i.e. the corrected narrow lane ambiguity floating point value, lambda w And λ n For BDS wide and narrow lane wavelengths, to
For inputting an observed value, the whole cycle ambiguity of the narrow lane and the phase decimal deviation of the satellite and the narrow lane of the receiver are taken as estimation parameters, and the design matrix of the constructed network solution narrow lane phase decimal deviation function model is the same as the formula (6);
step 42, referring to the step in the step 3, keeping the phase decimal deviation estimation of the narrow lane consistent with the phase decimal deviation estimation of the wide lane, facilitating the program expansion and use, and outputting the phase decimal deviation B of the narrow lane of the satellite at the moment 1 s 。
Preferably, in step 5, the step of linearly converting the phase fractional deviation of each frequency specifically includes the following steps:
turning the phase decimal deviations of the wide lane and the narrow lane obtained in the steps 3 and 4 to each frequency phase decimal deviation and marking a pseudo-range observation value channel according to a formula (9) and a formula (10) by using the phase decimal deviations of the wide lane and the narrow lane obtained in the steps 4 and according to coefficients of corresponding frequencies of MW combinations and LC combinations in the formula (1), so that a user can realize wide-area single-station ambiguity fixation by adopting a flexible function model;
in the formula (I), the compound is shown in the specification,
and
the converted first frequency phase decimal deviation and the converted second frequency phase decimal deviation are respectively the real-time Beidou phase decimal deviation, and the DCB is the satellite pseudo-range differential code delay on the corresponding pseudo-range observation value channel.
The invention has the beneficial effects that: the method adopts a strict data quality control strategy, and realizes the rapid estimation of the time-varying phase decimal deviation of the multi-station network solution based on parallel calculation and a single observation value iteration updating algorithm, thereby not only improving the efficiency of the real-time Beidou phase decimal deviation estimation, but also improving the precision and the reliability of the estimated phase decimal deviation.
Drawings
FIG. 1 is a schematic flow chart of the method of the present invention.
FIG. 2 is a schematic flow chart of the station-by-station parallel data preprocessing and the floating ambiguity calculation according to the present invention.
Fig. 3 is a schematic flow chart of the time-varying wide lane and narrow lane fractional deviation fast estimation method of the present invention.
FIG. 4 is a diagram of independent ambiguity reference selection according to the present invention.
FIG. 5 is a pseudo code for filter variance-covariance fast update using OpenMP according to the present invention.
Detailed Description
As shown in fig. 1, a real-time Beidou phase fractional deviation fast estimation method based on parallel computation includes the following steps:
data stream time synchronization and validity check;
the method specifically comprises the following steps:
(1-1) acquiring local time, setting a phase decimal deviation update rate and data stream waiting time, and taking the phase decimal deviation update rate and the data stream waiting time as a time reference for acquiring synchronous real-time data streams from epoch to epoch;
(1-2) recovering real-time satellite orbit and clock error data by using a broadcast ephemeris and a real-time orbit/clock error correction number, and monitoring whether the data is missing or abnormal;
(1-3) detecting and updating a multi-system Differential Code Bias (DCB) file, an Earth Orientation Parameter (EOP) file, an antenna file and a tide file in a BLQ format, wherein the DCB file is used as external input data for subsequent data preprocessing;
and (1-4) aiming at the problem of observation value missing of part of the observation stations, particularly the phenomenon of observation value data interruption occurring in a plurality of continuous epochs, setting waiting time, wherein the observation station data does not participate in the subsequent phase decimal deviation calculation in the waiting time. After the step is completed, the synchronous real-time observation value can be obtained, the external required file can be updated, and reliable data support is prepared for the station-by-station parallel calculation in the step 2.
(II) parallel computing floating ambiguity by multiple measuring stations;
the flow of the station-by-station parallel data preprocessing and the floating ambiguity calculation is shown in fig. 2, and the steps specifically include:
(2-1) deleting the dual-frequency pseudo range and phase observation value missing satellites, and forming the following dual-frequency linear combination observation values one by one satellite:
in the formula (f) 1 And f 2 Representing the first and second frequencies, P, of the BDS, respectively 1 /P 2 And L 1 /L 2 Respectively being a first frequency pseudorange observation value and a second frequency pseudorange observation value and a phase observation value of the BDS, PI being a difference value of the pseudorange observation values of the first frequency and the second frequency, LI being a difference value of the phase observation values of the first frequency and the second frequency, MW (Melbourne-Wu) observation value being a dual-frequency carrier phase L 1 /L 2 And pseudorange observations P 1 /P 2 And C is a dual-frequency ionosphere-free pseudo-range combined observed value, and LC is a dual-frequency ionosphere-free phase combined observed value.
(2-2) performing pseudo-range coarse difference detection by considering the orbit height and PI value threshold (60 m) of the GEO/IGSO/MEO three orbit type satellites, detecting the cycle slip by using MW combination and LI combination, considering the nearly static ground characteristic of the GEO satellites, setting the ionospheric change rate threshold of the LI linear combination detection cycle slip to be half of that of the MEO/IGSO, removing the observation data exceeding the threshold, and recording the number of arc sections of each continuous observation arc section and the switching time of new and old arc sections;
(2-3) correcting various error models including tropospheric delay, antenna phase center correction, phase wrapping, tidal effects, satellite multipath spatial error (this error correction is only for BDS-2IGSO/MEO satellites), etc., where the GEO satellite attitudeAdopting a zero bias mode, adopting a dynamic bias mode by an IGSO/MEO satellite, and calculating an epoch-by-epoch MW average value m by adopting a mobile sliding window mw And wherein the error σ mw ：
Where k represents the window size within a continuous arc segment (the maximum moving window size is set to 120).
(2-4) setting prior variances of observed values of the GEO/IGSO/MEO three orbit type satellites according to a formula (3), and respectively setting variance empirical values sigma of pseudo ranges and phase observed values of the three types of satellites due to different environments of a receiver and different orbit satellites 0 3/0.3/0.3m and 0.03/0.003/0.003m, where e represents the satellite height angle.
(2-5) correcting the PC and LC after various model errors in the step (2-3), adopting an ionosphere-free combined PPP model in the formula (4), judging and reducing the weight of an abnormal observation value according to a standardized residual error by using an IGG-III weight function model, and thus obtaining the reliable ionosphere-free combined floating-point ambiguity A if 。
In the formula (dt) r Indicating the receiver clock error, Z, of the station r r And m represents the wet delay and wet delay projection functions of the zenith troposphere of the survey station, epsilon PC And ε LC Indicating PC and LC observations, respectively, without modeling error.
(2-6) when observation data of a plurality of observation stations exist, openMP parallel processing is adopted according to the steps (2-1) - (2-5), and the synchronized multi-observation-station MW average value and the ionosphere-free combined floating ambiguity A can be obtained if Mean MW and ionosphere free combined float ambiguities obtained from this stepAnd (5) respectively using the degrees as input observed values of the network solution model in the subsequent step (three) and the step (four).
(III) quickly estimating the time-varying wide lane phase decimal deviation;
the flow of time-varying wide lane decimal deviation fast estimation is shown in fig. 3, and the steps specifically include:
(3-1) taking the average value of the MW of the multiple stations calculated in the step (two) as an input observation value, the wide-lane floating ambiguity can be expressed as:
in the formula (I), the compound is shown in the specification,
denotes the mean value of MW, lambda, of the multiple stations w Indicating a wide-lane wavelength, B w,* Indicating the receiver wide-lane phase fractional deviation,
indicating satellite wide-lane phase fractional deviation, N w Is the corresponding whole-circumference wide-lane ambiguity.
Assuming that q receivers and p satellites are observed, and taking the whole-cycle ambiguity of the wide lane and the phase decimal deviation of the satellites and the wide lane of the receivers as estimation parameters, constructing a network solution wide lane phase decimal deviation function model as follows:
in which I denotes the identity matrix, the subscript denotes its dimension, e p A column vector with elements of 1, the subscripts indicating the number of rows,
and a Kalman filtering design matrix representing the network solution wide lane phase decimal deviation.
(3-2) the ambiguity of the continuous arc segment is a constant, the process noise should be 0, and the continuous arc segment will be connectedThe ambiguity parameter of the continuous arc section is set as a constant model, and the small number deviation of the wide lane phase of the receiver and the small number deviation of the wide lane of the satellite receiver are respectively set as random walk due to the time-varying characteristic of the small number deviation parameter of the receiver and the satellite, and the spectral density of the random walk is 1.0e -4 And 1.0e -6 m 2 And/s, performing scalar calculation of Kalman filtering time updating by using the diagonal characteristic of the state transition matrix, namely performing state prediction of the current epoch on the basis of parameter estimation of the last epoch so as to realize rapid numerical calculation.
(3-3) mean error σ of smoothed values in terms of MW mw Performing ascending order arrangement, forming a non-directional connected graph by using floating ambiguity of a receiver, a satellite and a wide lane, introducing a Kruskal minimum spanning tree algorithm to generate a stable independent wide lane ambiguity standard, wherein FIG. 4 is an independent ambiguity standard selection schematic diagram, a square frame and a circle respectively represent a station and the satellite, and a black bold line represents that the ambiguity is selected as the independent ambiguity standard;
(3-4) As a virtual observation value, a value of a widelane ambiguity reference is set to have a virtual variance of 1.0e -9 The satellite phase fractional deviation is set to 0, i.e.
Adding the constraint to a wide lane phase decimal deviation net solution model, and estimating wide lane phase decimal deviation of each satellite relative to the reference constraint by using Kalman filtering;
(3-5) taking the MW smooth value of each satellite of each station and the selected widelane ambiguity reference value as a single observation value, and realizing Kalman filtering scalar calculation according to a formula (7), thereby avoiding the time-consuming problem of high-dimensional matrix inversion of standard filtering; the filtering variance-covariance matrix is quickly updated by adopting OpenMP parallel computation, and pseudo codes of the filtering variance-covariance matrix are shown in FIG. 5, so that the computation efficiency of network solution phase decimal deviation is further accelerated;
in the formula (I), the compound is shown in the specification,
for the variance of the ith satellite observation noise for the kth epoch,
and
respectively representing a Kalman filtering state prediction value and a filtering value,
and
are respectively a variance matrix corresponding to the states, I represents an identity matrix, k k,i Representing a Kalman filter gain matrix,/ k,i An ith satellite observation, a, representing the kth epoch k,i Is in the formula (6)
And (5) a row vector corresponding to the ith satellite observation value of the observation station is formed.
(3-6) based on the widelane ambiguity estimation value deviation (0.15 week), the standard deviation (0.15 week) and the probability (1000) that the ambiguity is fixed as an integer, if the widelane ambiguity obtained by filtering a certain satellite according to (3-5) meets the requirements, fixing the widelane ambiguity as an integer, traversing all estimated widelane ambiguities until no new widelane ambiguity can be fixed, and outputting the phase decimal deviation of the satellite at the moment
The fixed widelane ambiguity obtained in this step is a precondition for step (four).
(IV) quickly estimating the time-varying narrow lane phase decimal deviation;
similar to the previous step, the flow of time-varying narrow lane decimal deviation fast estimation is more consistent with that of fig. 3, and the main difference is that data is input, and the step specifically includes:
(4-1) correcting the ionosphere-free floating ambiguity obtained in the step (two) one by taking the wide-lane integer ambiguity fixed in the step (three) as a known value as shown in a formula (8):
in the formula (I), the compound is shown in the specification,
i.e. the corrected narrow lane ambiguity floating point value, lambda w And λ n BDS wide and narrow lane wavelengths. To be provided with
And (3) taking the narrow-lane whole-cycle ambiguity and the satellite and receiver narrow-lane phase fractional deviation as estimation parameters for inputting an observed value, wherein the design matrix of the constructed network solution narrow-lane phase fractional deviation function model is the same as the formula (6).
(4-2) - (4-6) referring to the step (III), so that the narrow lane phase decimal deviation estimation and the wide lane phase decimal deviation estimation are kept consistent, the program expansion and use are convenient, and the satellite narrow lane phase decimal deviation at the moment is output
(V) linearly converting each frequency phase decimal deviation;
turning the phase decimal deviations of the wide lane and the narrow lane obtained in the third step and the fourth step to each frequency phase decimal deviation and marking a pseudo-range observation value channel according to the formulas (9) and (10) by using the phase decimal deviations of the wide lane and the narrow lane obtained in the third step and the fourth step and according to the coefficient of the corresponding frequency of the MW combination and the LC combination in the step (1), so that a user can realize wide-area single-station ambiguity fixation by adopting a flexible function model;
in the formula (I), the compound is shown in the specification,
and
the converted first frequency phase decimal deviation and the converted second frequency phase decimal deviation are respectively the real-time Beidou phase decimal deviation, and the DCB is the satellite pseudo-range differential code deviation on the corresponding pseudo-range observation value channel.
In conclusion, the method adopts a strict data quality control strategy, and based on parallel computation and a single observation value iteration updating algorithm, the time-varying phase decimal deviation fast estimation of the multi-station network solution is realized, so that the efficiency of the real-time Beidou phase decimal deviation estimation is improved, and the precision and the reliability of the estimated phase decimal deviation are improved.
Claims (6)
1. A real-time Beidou phase decimal deviation fast estimation method based on parallel computing is characterized by comprising the following steps:
step 1, data stream time synchronization and validity check;
step 2, updating external required files according to the synchronous real-time observed values obtained in the step 1, and calculating the MW average value and the ionosphere-free combined floating-point ambiguity in parallel by multiple observation stations;
step 3, taking the MW average value of the multiple measuring stations obtained in the step 2 as an input observation value, and quickly estimating the phase decimal deviation of the time-varying wide lane;
step 4, taking the ionosphere-free combined floating-point ambiguity of the multiple measuring stations obtained in the step 2 as an input observation value, and performing rapid estimation on the time-varying narrow lane phase decimal deviation;
and 5, linearly converting the decimal deviation of each frequency phase.
2. The parallel computation-based real-time Beidou phase decimal deviation fast estimation method according to claim 1, characterized in that in the step 1, the data stream time synchronization and validity check specifically comprises the following steps:
step 11, acquiring local time, setting phase decimal deviation update rate and data stream waiting time, and taking the phase decimal deviation update rate and the data stream waiting time as time references for acquiring synchronous real-time data streams from epoch to epoch;
step 12, recovering the orbit and clock error data of the real-time satellite by using the broadcast ephemeris and the real-time orbit/clock error correction number, and monitoring whether the data is missing or abnormal;
step 13, detecting and updating a multi-system differential code deviation DCB file, a forecast earth orientation parameter EOP file, an antenna file and a sea tide file in a BLQ format, and taking the files as external input data of subsequent data preprocessing;
and step 14, aiming at the problem of observation value loss of part of observation stations, continuously generating an observation value data interruption phenomenon in a plurality of epochs, and setting waiting time, wherein the observation station data does not participate in subsequent phase decimal deviation calculation in the waiting time.
3. The parallel computation-based real-time Beidou phase decimal deviation fast estimation method according to claim 1, wherein in the step 2, the parallel computation of the floating ambiguity by the multi-observation station specifically comprises the following steps:
step 21, deleting the double-frequency pseudo range and phase observation value missing satellites, and forming the following double-frequency linear combination observation values one by one:
in the formula, f 1 And f 2 Representing the first and second frequencies, P, of the BDS, respectively 1 /P 2 And L 1 /L 2 First and second frequency pseudorange observations and phase observations of the BDS, respectively, PI being a difference between the pseudorange observations at the first and second frequencies, LI being a difference between the phase observations at the first and second frequencies, from the dual-frequency carrier phase L 1 /L 2 And pseudorange observations P 1 /P 2 The MW combined observed value of the structure is PC, the double-frequency non-ionized layer pseudo-range combined observed value and LCAn frequent ionosphere-free phase combination observed value;
step 22, regarding the orbit height and PI value threshold of the GEO/IGSO/MEO three orbit type satellites to carry out pseudo-range gross error detection, detecting the cycle slip by using MW combination and LI combination, considering the nearly geostationary property of the GEO satellites, setting the ionosphere change rate threshold of the LI linear combination detection cycle slip to be half of the MEO/IGSO, removing the observation data exceeding the threshold, and recording the number of arc sections of each continuous observation arc section and the switching time of new and old arc sections;
step 23, correcting various error models including tropospheric delay, antenna phase center correction, phase winding, tidal effect and satellite multipath space error, wherein the GEO satellite attitude adopts a zero-bias mode, the IGSO/MEO satellite adopts a dynamic-bias mode, and a moving sliding window is adopted to calculate an epoch-by-epoch MW average value m mw And wherein the error σ mw ：
Wherein k represents the window size within a continuous arc segment, and the maximum moving window size is set to 120;
step 24, setting prior variances of observed values of the GEO/IGSO/MEO three orbit type satellites according to a formula (3), and setting variance empirical values sigma of pseudo ranges and phase observed values of the three types of satellites respectively due to different environments of the receiver and the different orbit satellites 0 3/0.3/0.3m and 0.03/0.003/0.003m, wherein e represents the satellite height angle;
step 25, using the PC and LC corrected by the errors of the models in the step 23, adopting the ionosphere-free combined PPP model in the formula (4), using the IGG-III weight function model, judging and reducing the weight of the abnormal observed value according to the standardized residual error, thereby obtaining the reliable ionosphere-free combined floating ambiguity A if ；
In the formula (dt) r Indicating the receiver clock error, Z, of the station r r And m represents the wet delay and wet delay projection functions of the zenith troposphere of the survey station, epsilon PC And ε LC Respectively representing the unmodeled errors of the PC and LC observed values;
and step 26, when there are observation data of a plurality of observation stations, according to the steps from step 21 to step 25, adopting OpenMP parallel processing to obtain the synchronized multi-station MW average value and the ionosphere-free combined floating ambiguity.
4. The parallel computation-based real-time Beidou phase decimal deviation fast estimation method according to claim 1, characterized in that in the step 3, the time-varying wide lane phase decimal deviation fast estimation specifically comprises the following steps:
step 31, taking the average value of the MW of the multiple observation stations obtained in step 2 as an input observation value, and expressing the wide lane floating ambiguity as follows:
in the formula (I), the compound is shown in the specification,
denotes the mean value of MW, lambda, of the multiple stations w Indicating a wide-lane wavelength, B w,* Indicating the receiver wide-lane phase fractional deviation,
indicating satellite wide-lane phase fractional deviation, N w Is the corresponding whole-circumference wide lane ambiguity;
assuming that q receivers and p satellites are observed together, and taking the whole-cycle ambiguity of the wide lane and the phase decimal deviation of the wide lane of the satellites and the receivers as estimation parameters, constructing a wide lane phase decimal deviation function model as follows:
in which I denotes the identity matrix, the subscript denotes its dimension, e p A column vector with elements of 1, the subscripts indicating the number of rows,
a Kalman filtering design matrix representing network solution wide lane phase decimal deviation;
step 32, the ambiguity of the continuous arc section is a constant, the process noise should be 0, therefore the ambiguity parameter of the continuous arc section is set as a constant model, the small deviation of the wide lane phase of the receiver and the small deviation of the wide lane of the satellite receiver are respectively set as random walk due to the time-varying characteristic of the small deviation parameter of the receiver and the satellite, and the spectral density is respectively 1.0e -4 And 1.0e -6 m 2 The method comprises the following steps that (1) the scalar calculation of Kalman filtering time updating is carried out by utilizing the diagonal characteristic of a state transition matrix, namely, the state prediction of the current epoch is carried out on the basis of the parameter estimation of the last epoch, so that the rapid numerical calculation is realized;
step 33, mean error σ of smoothed value according to MW mw Performing ascending order arrangement, forming an undirected connected graph by using the floating ambiguity of the receiver, the satellite and the wide lane, and introducing a Kruskal minimum spanning tree algorithm to generate a steady independent wide lane ambiguity standard;
step 34, taking the wide lane ambiguity reference value as a virtual observation value, and setting the virtual variance of the wide lane ambiguity reference value as 1.0e -9 The satellite phase fractional deviation is set to 0, i.e.
Adding the constraint to a wide lane phase decimal deviation net solution model, and estimating wide lane phase decimal deviation of each satellite relative to the reference constraint by using Kalman filtering;
step 35, regarding the MW smooth value of each satellite of each station and the selected wide lane ambiguity reference value as a single observation value, and realizing Kalman filtering scalar calculation according to a formula (7), thereby avoiding the time-consuming problem of high-dimensional matrix inversion of standard filtering; the filtering variance-covariance matrix is rapidly updated by adopting OpenMP parallel computation, and the computation efficiency of network solution phase decimal deviation is further accelerated;
in the formula (I), the compound is shown in the specification,
is the variance of the ith satellite observation for the kth epoch,
and
respectively representing a Kalman filtering state prediction value and a filtering value,
and
are respectively a variance matrix corresponding to the states, I represents an identity matrix, k k,i Representing a Kalman filter gain matrix,/ k,i An ith satellite observation, a, representing the kth epoch k,i Is as in formula (6)
A row vector corresponding to the observation value of the ith satellite of the observation station is formed in a matrix;
step 36, based on the widelane ambiguity estimation value deviation of 0.15 week, the standard deviation of 0.15 week and the probability 1000 of ambiguity fixed as an integer, if the widelane ambiguity obtained by filtering a certain satellite according to step 35 meets the requirements, fixing the widelane ambiguity as an integerTraversing all estimated widelane ambiguities until no new widelane ambiguities can be fixed, and outputting the phase decimal deviation of the satellite widelane at the moment
5. The parallel computation-based real-time Beidou phase decimal deviation fast estimation method according to claim 1, characterized in that in the step 4, the fast estimation of the time-varying narrow lane phase decimal deviation specifically comprises the following steps:
step 41, taking the whole-cycle ambiguity of the wide lane fixed in step 3 as a known value, and correcting the ionosphere-free floating ambiguity obtained in step 2 one by one satellite, as shown in formula (8):
in the formula (I), the compound is shown in the specification,
i.e. the corrected narrow lane ambiguity floating point value, lambda w And λ n For BDS wide and narrow lane wavelengths, to
For inputting an observed value, the whole cycle ambiguity of the narrow lane and the phase decimal deviation of the satellite and the narrow lane of the receiver are taken as estimation parameters, and the design matrix of the constructed network solution narrow lane phase decimal deviation function model is the same as the formula (6);
step 42, referring to the step in the step 3, keeping the phase decimal deviation estimation of the narrow lane consistent with the phase decimal deviation estimation of the wide lane, facilitating the program expansion and use, and outputting the phase decimal deviation B of the narrow lane of the satellite at the moment 1 s 。
6. The real-time Beidou phase decimal deviation fast estimation method based on parallel computing as claimed in claim 1, characterized in that in step 5, the step of linearly converting the phase decimal deviation of each frequency specifically comprises the following steps:
turning the phase decimal deviations of the wide lane and the narrow lane obtained in the steps 3 and 4 to each frequency phase decimal deviation and marking a pseudo-range observation value channel according to a formula (9) and a formula (10) by using the phase decimal deviations of the wide lane and the narrow lane obtained in the steps 4 and according to coefficients of corresponding frequencies of MW combinations and LC combinations in the formula (1), so that a user can realize wide-area single-station ambiguity fixation by adopting a flexible function model;
in the formula (I), the compound is shown in the specification,
and
the converted first frequency phase decimal deviation and the converted second frequency phase decimal deviation are respectively the real-time Beidou phase decimal deviation, and the DCB is the satellite pseudo-range differential code delay on the corresponding pseudo-range observation value channel.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211076810.6A CN115561793A (en) | 2022-09-05 | 2022-09-05 | Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211076810.6A CN115561793A (en) | 2022-09-05 | 2022-09-05 | Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115561793A true CN115561793A (en) | 2023-01-03 |
Family
ID=84738556
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211076810.6A Pending CN115561793A (en) | 2022-09-05 | 2022-09-05 | Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115561793A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20210223353A1 (en) * | 2018-09-13 | 2021-07-22 | L3Harris Technologies, Inc. | Opportunistic adjustable rate cross-ambiguity function geolocation |
-
2022
- 2022-09-05 CN CN202211076810.6A patent/CN115561793A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20210223353A1 (en) * | 2018-09-13 | 2021-07-22 | L3Harris Technologies, Inc. | Opportunistic adjustable rate cross-ambiguity function geolocation |
| US11988765B2 (en) * | 2018-09-13 | 2024-05-21 | L3Harris Technologies, Inc. | Opportunistic adjustable rate cross-ambiguity function geolocation |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109581452B (en) | GNSS reference station carrier phase integer ambiguity resolution method | |
| CN108549095B (en) | Non-differential parallel enhancement method and system for regional CORS network | |
| AU2008260578B2 (en) | Distance dependant error mitigation in real-time kinematic (RTK) positioning | |
| CN108196284B (en) | GNSS network data processing method for fixing single-difference ambiguity between satellites | |
| CN104714244A (en) | Multi-system dynamic PPP resolving method based on robust self-adaption Kalman smoothing | |
| EP2663878A1 (en) | Navigation system and method for resolving integer ambiguities using double difference ambiguity constraints | |
| CN111965673A (en) | Time frequency transfer method of single-frequency precise single-point positioning algorithm based on multiple GNSS | |
| CN111638535B (en) | Hybrid ambiguity fixing method for GNSS real-time precise point positioning | |
| EP3430437A1 (en) | Navigation satellite wide-lane bias determination system and method | |
| CN112394376B (en) | Large-scale GNSS network observation data non-differential whole network parallel processing method | |
| CN111290005A (en) | Differential positioning method and device for carrier phase, electronic equipment and storage medium | |
| Wang et al. | Evaluating the impact of CNES real-time ionospheric products on multi-GNSS single-frequency positioning using the IGS real-time service | |
| Liu et al. | Comparison of convergence time and positioning accuracy among BDS, GPS and BDS/GPS precise point positioning with ambiguity resolution | |
| CN113358017B (en) | Multi-station cooperative processing GNSS high-precision deformation monitoring method | |
| CN112859120A (en) | Continuous GNSS carrier phase time and frequency transfer method | |
| CN115933356B (en) | High-precision time synchronization system and method for virtual atomic clock | |
| CN115373005A (en) | High-precision product conversion method between satellite navigation signals | |
| CN114779301B (en) | Satellite navigation real-time precise single-point positioning method based on broadcast ephemeris | |
| Kuang et al. | Real-time GPS satellite orbit and clock estimation based on OpenMP | |
| CN115902968A (en) | PPP terminal positioning method based on Beidou third GEO broadcast enhancement information | |
| CN115561793A (en) | Real-time Beidou phase decimal deviation rapid estimation method based on parallel computation | |
| CN112230254B (en) | Correction method and device for GPS carrier phase multipath error | |
| Wang et al. | Comparison of three widely used multi‐GNSS real‐time single‐frequency precise point positioning models using the International GNSS Service real‐time service | |
| CN114460615B (en) | Beidou three-new frequency point positioning method and system with additional virtual observation value | |
| Kuang et al. | Galileo real-time orbit determination with multi-frequency raw observations |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |