A Fast Estimation Method for Real-Time Beidou Phase Fractional Bias Based on Parallel Computing

Abstract

The present invention discloses a method for rapidly estimating the fractional phase bias of Beidou in real time based on parallel computing, comprising the following steps: Step 1, data stream time synchronization and validity verification; Step 2, according to the synchronized real-time observations obtained in Step 1 and updating the required external files, multi-station parallel computing of the MW average value and the ionosphere-free combined floating ambiguity; Step 3, according to the multi-station MW average value obtained in Step 2, as the input observation value, rapid estimation of the time-varying wide-lane fractional phase bias; Step 4, according to the multi-station ionosphere-free combined floating ambiguity obtained in Step 2, as the input observation value, rapid estimation of the time-varying narrow-lane fractional phase bias; Step 5, linearly converting the fractional phase bias of each frequency. The present invention realizes the rapid estimation of the time-varying phase fractional bias of Beidou in multi-station network solutions, which not only improves the efficiency of real-time Beidou fractional phase bias estimation, but also improves the accuracy and reliability of the estimated Beidou fractional phase bias.

Description

Technical Field

[0001] The present invention relates to the technical field of satellite navigation and positioning, and in particular to a method for rapidly estimating the fractional phase bias of Beidou in real time based on parallel computing. Background Art

[0002] The Beidou Navigation Satellite System (BDS) is a global navigation satellite system (GNSS) independently built and operated by China, and is an important national space infrastructure. Currently, the annual output value of the industry based on Beidou location services has reached as high as 400 billion yuan, and shows a stable and rapid growth trend. Reliable spatial position information is required for unmanned driving, smart city construction, natural resource surveys, etc.

[0003] Precise Point Positioning (PPP) technology [1] can obtain precise three-dimensional coordinates in the geocentric reference frame globally. However, since its ambiguity parameters absorb the pseudorange and phase hardware delays, the accuracy and reliability of the PPP floating-point solution that loses the integer-week characteristic need to be further improved.

[0004] To achieve reliable fixing of PPP ambiguities, scholars at home and abroad have successively proposed fractional phase bias models [2] , integer phase clock models [3] , clock bias decoupling models [4] etc., and compared the above three methods in terms of the type of server correction and the degree of freedom of the client model. Since the construction period of the BDS system is relatively short and three types of hybrid orbit satellites are adopted, including geostationary earth orbit satellites (GEO), inclined geosynchronous orbit satellites (IGSO), and medium earth orbit satellites (MEO), there are certain differences between the satellite observation values, satellite attitude models, satellite space environments, etc. and the rest of the GNSS systems. Ignoring the differences between the two reduces the BDS PPP positioning performance.

[0005] At the same time, in order to meet the timeliness and long-term stability of real-time position services, the efficiency and numerical stability of BDS time-varying fractional phase bias estimation are particularly important. In recent years, many scholars have adopted methods such as introducing external matrix libraries, non-rigorous network solution models, and reducing the fractional phase bias update time to improve the efficiency of fractional phase bias estimation. However, the above technical means rely on hardware platforms, are difficult to achieve cross-platform operation, reduce the accuracy and reliability of fractional phase bias, are not conducive to reliable fixing of BDS PPP ambiguities, and the rigorous network solution model has problems such as high-dimensional matrix inversion and time-consuming variance-covariance.

[0006] Therefore, a reliable estimation of the real-time BDS phase fractional bias to achieve a reliable BDS PPP fixed solution has important practical value and research significance for promoting the Beidou system in China to provide high-quality navigation, positioning, and timing services. Summary of the Invention

[0007] The technical problem to be solved by the present invention is to provide a method for rapidly estimating the real-time Beidou phase fractional bias based on parallel computing, which can achieve rapid and robust estimation of the real-time BDS phase fractional bias.

[0008] To solve the above technical problem, the present invention provides a method for rapidly estimating the real-time Beidou phase fractional bias based on parallel computing, including the following steps:

[0009] Step 1: Data stream time synchronization and validity verification;

[0010] Step 2: According to the synchronized real-time observations obtained in Step 1 and update the external required files, perform parallel computing of the MW average value and the ionosphere-free combined floating ambiguity for multiple stations;

[0011] Step 3: Using the MW average value of multiple stations obtained in Step 2 as the input observation value, perform rapid estimation of the time-varying wide-lane phase fractional bias;

[0012] Step 4: Using the ionosphere-free combined floating ambiguity of multiple stations obtained in Step 2 as the input observation value, perform rapid estimation of the time-varying narrow-lane phase fractional bias;

[0013] Step 5: Linearly transform the phase fractional biases of each frequency.

[0014] Preferably, in Step 1, the data stream time synchronization and validity verification specifically include the following steps:

[0015] Step 11: Obtain the local time, set the phase fractional bias update rate and the data stream waiting time, and use them as the time reference for obtaining the synchronized real-time data stream epoch by epoch;

[0016] Step 12: Use the broadcast ephemeris and real-time orbit / clock correction data to restore the real-time satellite orbit and clock data, and monitor whether there are any missing and abnormal data;

[0017] Step 13: Detect and update the differential code bias (DCB) files of multiple systems, the predicted Earth orientation parameter (EOP) files, the antenna files, and the ocean tide files in BLQ format, as the external input data for subsequent data preprocessing;

[0018] Step 14: For the problem of missing observations at some stations, especially the phenomenon of data interruption in consecutive epochs, set a waiting time, and the data of this station will not participate in the subsequent phase fractional bias calculation during the waiting time.

[0019] Preferably, in step 2, the multi-station parallel calculation of floating-point ambiguities specifically includes the following steps:

[0020] Step 21, delete the satellites with missing dual-frequency pseudorange and phase observation values, and form the following dual-frequency linear combination observation values ​​one by one:

[0021]

[0022] In the formula, f 1 and f 2 Represent the first and second frequencies of BDS, P 1 / P 2 and L 1 / L 2 They are the pseudorange observation values ​​and phase observation values ​​of the first and second frequencies of BDS, PI is the difference between the pseudorange observation values ​​of the first and second frequencies, LI is the difference between the phase observation values ​​of the first and second frequencies, and the dual-frequency carrier phase L 1 / L 2 and pseudorange observation value P 1 / P 2 The constructed MW combined observations, PC is the dual-frequency ionospheric-free pseudorange combined observations, and LC is the dual-frequency ionospheric-free phase combined observations;

[0023] Step 22, taking into account the orbital altitudes and PI value thresholds of the three orbital types of satellites GEO / IGSO / MEO to detect pseudorange gross errors, using the MW combination and the LI combination to detect cycle slips, considering the nearly geostationary characteristics of the GEO satellite, the ionospheric change rate threshold of the LI linear combination to detect cycle slips is set to half of that of MEO / IGSO, and the observation data exceeding the threshold is removed, and the number of arcs and the switching time of the new and old arcs of each continuous observation arc are recorded;

[0024] Step 23: Correct various error models, including tropospheric delay, antenna phase center correction, phase winding, tidal effect, and satellite multipath spatial error. The GEO satellite attitude adopts the zero bias mode, the IGSO / MEO satellite adopts the dynamic bias mode, and the moving sliding window is used to calculate the epoch-by-epoch MW average m mw and its error σ mw :

[0025]

[0026] In the formula, k represents the window size within the continuous arc segment, and the maximum moving window size is set to 120;

[0027] Step 24: Set the prior variances of the observations of GEO / IGSO / MEO satellites according to formula (3). Since the environments of the receiver and satellites in different orbits are different, the empirical variances of the pseudorange and phase observations of the three types of satellites are set to be 3 / 0.3 / 0.3 m and 0.03 / 0.003 / 0.003 m respectively, where e represents the satellite elevation angle; 0 ;

[0028]

[0029] Step 25: Use the PC and LC after correcting various model errors in Step 23, adopt the ionosphere-free combined PPP model in formula (4), use the IGG-III weight function model, and judge and reduce the weights of abnormal observations according to the standardized residuals, so as to obtain a reliable ionosphere-free combined float ambiguity A if ;

[0030]

[0031] In the formula, dt r represents the receiver clock error of station r, Z r and m represent the zenith tropospheric wet delay and wet delay projection function of the station respectively, ε PC and ε LC represent the unmodeled errors of the PC and LC observations respectively;

[0032] Step 26: When there are observation data of multiple stations, according to the steps of Step 21 - Step 25, adopt OpenMP parallel processing to obtain the synchronized MW average value of multiple stations and the ionosphere-free combined float ambiguity.

[0033] Preferably, in Step 3, the fast estimation of the time-varying wide-lane phase fractional bias specifically includes the following steps:

[0034] Step 31: Use the MW average value of multiple stations obtained in Step 2 as the input observation value, then the wide-lane float ambiguity is expressed as:

[0035]

[0036] In the formula, represents the MW average value of multiple stations, λ w represents the wide-lane wavelength, B w,* represents the receiver wide-lane phase fractional bias, represents the satellite wide-lane phase fractional bias, N w is the corresponding integer wide-lane ambiguity;

[0037] Assume that a total of q receivers and p satellites are observed. Taking the wide-lane integer ambiguity and the wide-lane phase fractional biases of satellites and receivers as estimation parameters, the following function model for solving the wide-lane phase fractional bias of the network is constructed:

[0038]

[0039] In the formula, I represents the identity matrix, and the subscript represents its dimension. e p represents a column vector with all elements being 1, and the subscript represents its number of rows. represents the Kalman filter design matrix for solving the wide-lane phase fractional bias of the network;

[0040] Step 32: The ambiguity of the continuous arc segment is a constant, and its process noise should be 0. Therefore, the ambiguity parameter of the continuous arc segment is set to a constant model. Due to the time-varying characteristics of the fractional bias parameters of receivers and satellites, the wide-lane phase fractional bias of receivers and the wide-lane fractional bias of satellite receivers are respectively set to random walks, and their spectral densities are 1.0e -4 and 1.0e -6 m 2 / s. Using the diagonal property of the state transition matrix, scalar calculations for Kalman filter time update are performed, that is, based on the parameter estimation of the previous epoch, the state prediction of the current epoch is carried out to achieve fast numerical calculations;

[0041] Step 33: Arrange in ascending order according to the mean square error σ mw of the MW smoothed value, and form an undirected connected graph with receivers, satellites, and wide-lane floating-point ambiguities. Introduce the Kruskal minimum spanning tree algorithm to generate a robust independent wide-lane ambiguity reference;

[0042] Step 34: Take the wide-lane ambiguity reference value as a virtual observation value, and set its virtual variance to 1.0e -9 , and set the satellite phase fractional bias to 0, that is, is added as a constraint to the wide-lane phase fractional bias network solution model, and the wide-lane phase fractional bias of each satellite relative to this reference constraint is estimated using the Kalman filter;

[0043] Step 35: Regard the MW smoothed value of each satellite at each station and the selected wide-lane ambiguity reference value as a single observation value, and implement Kalman filter scalar calculations according to formula (7), avoiding the time-consuming problem of high-dimensional matrix inversion in the standard filter; and use OpenMP parallel computing to achieve fast update of the filter variance-covariance matrix, further accelerating the calculation efficiency of the network solution phase fractional bias;

[0044]

[0045] In the formula, is the variance of the observations of the \(i\)-th satellite at the \(k\)-th epoch. and represent the predicted value and the filtered value of the Kalman filter state respectively. and are the variance matrices corresponding to the states respectively, \(I\) represents the identity matrix, and \(k\) k,i represents the Kalman filter gain matrix, and \(l\) k,i represents the observation of the \(i\)-th satellite at the \(k\)-th epoch, and \(a\) k,i is in formula (6) The row vector corresponding to the observation of the \(i\)-th satellite at this station in the matrix;

[0046] Step 36: Based on the wide-lane ambiguity estimation deviation of 0.15 cycles, the standard deviation of 0.15 cycles, and the probability of 1000 that the ambiguity is fixed as an integer, if the wide-lane ambiguity obtained by filtering a certain satellite according to Step 35 meets the above requirements, then fix it as an integer. Traverse all the estimated wide-lane ambiguities until no new wide-lane ambiguity can be fixed, and output the satellite wide-lane phase fractional deviation at this time The fixed wide-lane ambiguity obtained in this step is the premise of Step 4.

[0047] Preferably, in Step 4, the rapid estimation of the time-varying narrow-lane phase fractional deviation specifically includes the following steps:

[0048] Step 41: Taking the wide-lane integer ambiguity fixed in Step 3 as a known value, correct the ionosphere-free float ambiguity obtained in Step 2 for each satellite, as shown in formula (8):

[0049]

[0050] In the formula, is the floating-point value of the corrected narrow-lane ambiguity, \(\lambda\) w and \(\lambda\) n are the BDS wide-lane and narrow-lane wavelengths. Taking as the input observation, taking the narrow-lane integer ambiguity and the satellite and receiver narrow-lane phase fractional deviations as the estimation parameters, the design matrix of the network solution narrow-lane phase fractional deviation function model is the same as formula (6);

[0051] Step 42: Refer to the steps in Step 3 to make the narrow-lane phase fractional deviation estimation consistent with the wide-lane phase fractional deviation estimation, which is convenient for program extension and use, and output the satellite narrow-lane phase fractional deviation \(B\) 1 s .

[0052] Preferably, in Step 5, the linear conversion of the fractional deviations of each frequency phase specifically includes the following steps:

[0053] Using the wide-lane and narrow-lane phase fractional biases obtained in steps 3 and 4, according to the coefficients of the corresponding frequencies of the MW combination and the LC combination in formula (1), turn them into the fractional biases of each frequency phase according to formulas (9) and (10) and mark the pseudorange observation value channels, so that users can use a flexible function model to realize wide-area single-site ambiguity fixing;

[0054]

[0055]

[0056] In the formula, and are the fractional biases of the first frequency and the second frequency after conversion, that is, the real-time Beidou phase fractional bias, and DCB is the satellite pseudorange differential code delay on the corresponding pseudorange observation value channel.

[0057] The beneficial effects of the present invention are as follows: The present invention adopts a strict data quality control strategy, and based on parallel computing and single-observation iterative update algorithms, it realizes the rapid estimation of time-varying phase fractional biases in multi-site network solutions, which not only improves the efficiency of real-time Beidou phase fractional bias estimation, but also improves the accuracy and reliability of the estimated phase fractional biases. Description of the Drawings

[0058] Figure 1 is a schematic diagram of the method flow of the present invention.

[0059] Figure 2 is a schematic diagram of the process of per-site parallel data preprocessing and floating-point ambiguity calculation of the present invention.

[0060] Figure 3 is a schematic diagram of the process of rapid estimation of time-varying wide-lane and narrow-lane fractional biases of the present invention.

[0061] Figure 4 is a schematic diagram of the selection of an independent ambiguity reference of the present invention.

[0062] Figure 5 is the pseudocode of the present invention for realizing rapid update of filtering variance-covariance using OpenMP. Detailed Embodiments

[0063] As Figure 1 shown, a method for rapid estimation of real-time Beidou phase fractional bias based on parallel computing includes the following steps:

[0064] (1) Data stream time synchronization and validity verification;

[0065] This step specifically includes:

[0066] (1-1) Obtain the local time, set the phase decimal deviation update rate and the data stream waiting time, and use them as the time reference for obtaining the synchronous real-time data stream epoch by epoch;

[0067] (1-2) Use the broadcast ephemeris and real-time orbit / clock correction data to restore the real-time satellite orbit and clock data, and monitor whether there are any missing or abnormal data;

[0068] (1-3) Detect and update the differential code bias (DCB) files of multiple systems, predict the Earth Orientation Parameters (EOP) files, antenna files, and ocean tide files in BLQ format, and use them as external input data for subsequent data preprocessing;

[0069] (1-4) For the problem of missing observations at some stations, especially the phenomenon of data interruption in multiple consecutive epochs, set the waiting time, and the data of this station will not participate in the subsequent phase decimal deviation calculation during the waiting time. After completing this step, synchronous real-time observations can be obtained and the externally required files can be updated, preparing reliable data support for the parallel calculation of each station in step 2.

[0070] (2) Parallel calculation of float ambiguities for multiple stations;

[0071] The process of parallel data preprocessing and float ambiguity calculation for each station is as Figure 2 shown, and this step specifically includes:

[0072] (2-1) Delete the satellites with missing dual-frequency pseudorange and phase observations, and form the following dual-frequency linear combination observations for each satellite:

[0073]

[0074] In the formula, f 1 and f 2 represent the frequencies of the first and second frequencies of BDS respectively, P 1 / P 2 and L 1 / L 2 are the pseudorange observations and phase observations of the first and second frequencies of BDS respectively, PI is the difference between the pseudorange observations of the first and second frequencies, LI is the difference between the phase observations of the first and second frequencies, the MW (Melbourne-Wu¨bbena) observation is the combined observation constructed by the dual-frequency carrier phase L 1 / L 2 and the pseudorange observation P 1 / P 2 The combined observation PC is the dual-frequency ionosphere-free pseudorange combined observation, and LC is the dual-frequency ionosphere-free phase combined observation.

[0075] (2-2) Detect pseudorange gross errors by taking into account the orbital altitudes of GEO / IGSO / MEO satellites and the PI value threshold (60 m), detect cycle slips using the MW combination and the LI combination. Considering the near-geostationary characteristics of GEO satellites, the ionospheric change rate threshold for cycle slip detection using the LI linear combination is set to half of that for MEO / IGSO. Reject the observation data exceeding the threshold, and record the number of arc segments and the switching time between new and old arc segments for each continuous observation arc segment;

[0076] (2-3) Correct various error models, including tropospheric delay, antenna phase center correction, phase winding, tidal effect, satellite multipath spatial error (this error correction is only for BDS-2 IGSO / MEO satellites), etc. Among them, the GEO satellite attitude adopts the zero-bias mode, the IGSO / MEO satellites adopt the dynamic-bias mode, and a moving sliding window is used to calculate the epoch-by-epoch MW average m mw and its mean square error σ mw :

[0077]

[0078] In the formula, k represents the window size within the continuous arc segment (the maximum moving window size is set to 120).

[0079] (2-4) Set the prior variances of the observations of GEO / IGSO / MEO satellites according to formula (3). Since the environments of the receiver and satellites in different orbits are different, the empirical variances of the pseudorange and phase observations of the three types of satellites are respectively set as σ 0 to be 3 / 0.3 / 0.3 m and 0.03 / 0.003 / 0.003 m, where e represents the satellite elevation angle.

[0080]

[0081] (2-5) Using the PC and LC after correcting various model errors in step (2-3), adopt the ionosphere-free combination PPP model as in formula (4), use the IGG-III weight function model, and judge and reduce the weights of abnormal observations according to the standardized residuals, so as to obtain a reliable ionosphere-free combination floating-point ambiguity A if .

[0082]

[0083] In the formula, dt r represents the receiver clock error of station r, Z r and m respectively represent the wet delay of the troposphere at the zenith of the station and the wet delay projection function, ε PC and ε LC respectively represent the unmodeled errors of the PC and LC observations.

[0084] (2 - 6) When there is observation data from multiple stations, following the steps in (2 - 1) to (2 - 5) and using OpenMP parallel processing, the synchronized MW average values of multiple stations and the ionosphere - free combined float ambiguities A can be obtained. if The MW average values and ionosphere - free combined float ambiguities obtained in this step are respectively used as the input observed values for the network solution models in subsequent steps (III) and (IV).

[0085] (III) Fast estimation of time - varying wide - lane phase fractional biases;

[0086] The process of fast estimation of time - varying wide - lane fractional biases is as Figure 3 shown, and this step specifically includes:

[0087] (3 - 1) Taking the MW average values of multiple stations calculated in step (II) as input observed values, the wide - lane float ambiguity can be expressed as:

[0088]

[0089] In the formula, represents the MW average values of multiple stations, λ w represents the wide - lane wavelength, B w,* represents the receiver wide - lane phase fractional bias, represents the satellite wide - lane phase fractional bias, N w is the corresponding integer - cycle wide - lane ambiguity.

[0090] Assuming that q receivers and p satellites are observed in total, taking the wide - lane integer ambiguity and satellite and receiver wide - lane phase fractional biases as estimation parameters, the following network solution wide - lane phase fractional bias function model is constructed:

[0091]

[0092] In the formula, I represents the identity matrix, the subscript represents its dimension, e p represents a column vector with all elements being 1, the subscript represents its number of rows, represents the Kalman filter design matrix of the network solution wide - lane phase fractional bias.

[0093] (3 - 2) The ambiguity of a continuous arc segment is a constant, and its process noise should be 0. Set the ambiguity parameter of the continuous arc segment as a constant model. Due to the time - varying characteristics of the fractional bias parameters of the receiver and satellite, set the receiver wide - lane phase fractional bias and satellite receiver wide - lane fractional bias as random walks respectively, and their spectral densities are 1.0e -4 and 1.0e -6 m 2 / s, Using the diagonal property of the state transition matrix, scalar calculations for Kalman filter time updates are performed, that is, based on the parameter estimation of the previous epoch, the state prediction of the current epoch is carried out to achieve fast numerical calculations.

[0094] (3-3) Arrange in ascending order according to the mean square error σ of the MW smoothing value, and form an undirected connected graph with receivers, satellites, and wide-lane floating-point ambiguities. Introduce the Kruskal minimum spanning tree algorithm to generate a robust independent wide-lane ambiguity reference. mw For the schematic diagram of the selection of the independent ambiguity reference, the squares and circles represent stations and satellites respectively, and the thick black lines indicate that the ambiguity is selected as the independent ambiguity reference. Figure 4 For the schematic diagram of the selection of the independent ambiguity reference, the squares and circles represent stations and satellites respectively, and the thick black lines indicate that the ambiguity is selected as the independent ambiguity reference;

[0095] (3-4) Take the wide-lane ambiguity reference value as a virtual observation value, and set its virtual variance to 1.0e -9 , and set the satellite phase fractional bias to 0, that is Add it as a constraint to the wide-lane phase fractional bias network solution model, and use Kalman filter to estimate the wide-lane phase fractional bias of each satellite relative to this reference constraint.

[0096] (3-5) Regard the MW smoothing value of each satellite at each station and the selected wide-lane ambiguity reference value as a single observation value, and implement Kalman filter scalar calculations according to formula (7), avoiding the time-consuming problem of inverting high-dimensional matrices in standard filtering; and use OpenMP parallel computing to achieve fast update of the filter variance-covariance matrix, and its pseudocode is as Figure 5 shown to further accelerate the calculation efficiency of the network solution phase fractional bias;

[0097]

[0098] In the formula, is the variance of the observation noise of the i-th satellite at the k-th epoch, and respectively represent the Kalman filter state prediction value and the filter value, and are the variance matrices corresponding to the states respectively, I represents the identity matrix, k k,i represents the Kalman filter gain matrix, l k,i represents the observation value of the i-th satellite at the k-th epoch, a k,i is the row vector corresponding to the observation value of the i-th satellite at this station in the matrix in formula (6).

[0099] (3 - 6) Based on the bias (0.15 cycles), standard deviation (0.15 cycles) of the wide - lane ambiguity estimation value, and the probability (1000) of the ambiguity being fixed as an integer, if the wide - lane ambiguity obtained by filtering a certain satellite according to (3 - 5) meets the above requirements, it is fixed as an integer. All estimated wide - lane ambiguities are traversed until no new wide - lane ambiguity can be fixed, and the fractional phase deviation of the satellite wide - lane at this time is output. The fixed wide - lane ambiguity obtained in this step is the premise of step (iv).

[0100] (iv) Fast estimation of the time - varying fractional phase deviation of the narrow - lane;

[0101] Similar to the previous step, the process of fast estimation of the time - varying fractional deviation of the narrow - lane is Figure 3 relatively consistent, and the main difference lies in the input data. This step specifically includes:

[0102] (4 - 1) Taking the wide - lane integer ambiguity fixed in step (iii) as a known value, the ionosphere - free floating ambiguity obtained in step (ii) is corrected for each satellite, as shown in Equation (8):

[0103]

[0104] In the formula, is the floating - point value of the corrected narrow - lane ambiguity, λ w and λ n are the BDS wide - lane and narrow - lane wavelengths. Taking as the input observation value, with the narrow - lane integer ambiguity and the fractional phase deviations of the satellite and receiver narrow - lanes as the estimation parameters, the design matrix of the network solution fractional phase deviation function model of the narrow - lane is the same as Equation (6).

[0105] (4 - 2) - (4 - 6) Refer to the steps in (iii) to make the estimation of the fractional phase deviation of the narrow - lane consistent with that of the wide - lane, which is convenient for program extension and use, and output the fractional phase deviation of the satellite narrow - lane at this time.

[0106] (v) Linearly transform the fractional phase deviations of each frequency;

[0107] Using the fractional phase deviations of the wide - lane and narrow - lane obtained in step (iii) and step (iv), according to the coefficients of the corresponding frequencies of the MW combination and LC combination in (1), they are transformed to the fractional phase deviations of each frequency according to Equations (9) and (10) and the pseudorange observation value channels are marked, so that users can use a flexible function model to achieve wide - area single - station ambiguity fixing.

[0108]

[0109]

[0110] In the formula, and are respectively the first and second frequency phase decimal deviations after conversion, that is, the real-time Beidou phase decimal deviation, and DCB is the satellite pseudorange differential code deviation on the corresponding pseudorange observation value channel.

[0111] In summary, the present invention adopts a strict data quality control strategy, and based on parallel computing and a single observation value iterative update algorithm, realizes the rapid estimation of time-varying phase decimal deviations in multi-station network solutions, which not only improves the efficiency of real-time Beidou phase decimal deviation estimation, but also improves the accuracy and reliability of the estimated phase decimal deviation.

Claims

1. A fast estimation method for real-time Beidou phase fractional bias based on parallel computing, characterized in that, it includes the following steps: Step 1, data stream time synchronization and validity verification; Step 2, based on the synchronous real-time observations obtained in Step 1 and update the required external files, and perform parallel computing of MW mean values and ionosphere-free combined floating ambiguities at multiple stations; Step 3, taking the MW mean values of multiple stations obtained in Step 2 as input observations, and perform fast estimation of time-varying wide-lane phase fractional bias; specifically including the following steps: Step 31, taking the MW mean values of multiple stations obtained in Step 2 as input observations, then the wide-lane floating ambiguity is expressed as: In the formula, represents the multi-station MW average value, and λ w represents the wide-lane wavelength, B w represents the fractional deviation of the receiver's wide-lane phase, represents the fractional deviation of the satellite's wide-lane phase, and N w is the corresponding wide-lane integer ambiguity; Assume that a total of q receivers and p satellites are observed, and taking the wide-lane integer ambiguity and the wide-lane phase fractional biases of satellites and receivers as estimation parameters, construct the following network solution wide-lane phase fractional bias function model: where \(I\) represents the identity matrix, and the subscript represents its dimension, \(e\) p represents a column vector with all elements being 1, and the subscript represents its number of rows, represents the Kalman filter design matrix for the fractional deviation of the network solution wide-lane phase; \(B\) w,1 ... \(B\) w,q is the fractional deviation of the wide-lane phase of the first to the \(q\)th receivers, and is the fractional deviation of the wide-lane phase of the first to the \(p\)th satellites; Step 32: The ambiguity of the continuous arc segment is a constant, and its process noise should be 0. Therefore, the ambiguity parameter of the continuous arc segment is set to a constant model. Due to the time-varying characteristics of the fractional bias parameters of the receiver and the satellite, the fractional bias of the receiver wide-lane phase and the fractional bias of the satellite receiver wide-lane are respectively set to random walks, and their spectral densities are 1.0e -4 and 1.0e -6 m 2 / s. Using the diagonal property of the state transition matrix, scalar calculations for Kalman filter time updates are performed, that is, based on the parameter estimation of the previous epoch, the state prediction of the current epoch is carried out to achieve fast numerical calculations; Step 33: Ascendingly sort according to the mean square error σ of the MW smoothing value mw to form an undirected connected graph with the receiver, satellite, and wide-lane floating-point ambiguity, and introduce the Kruskal minimum spanning tree algorithm to generate a robust independent wide-lane ambiguity reference; Step 34: Use the wide-lane ambiguity reference value as a virtual observation, and set its virtual variance to 1.0e -9 , set the satellite phase fractional deviation to 0, that is Add it as a constraint to the wide-lane phase fractional deviation network solution model, and use Kalman filtering to estimate the wide-lane phase fractional deviation of each satellite relative to this reference constraint; Step 35, regarding the MW smoothed value of each satellite at each station and the selected wide-lane ambiguity reference value as a single observation value, and implementing Kalman filter scalar calculation according to formula (7), avoiding the time-consuming problem of high-dimensional matrix inversion in the standard filter; and using OpenMP parallel computing to achieve fast update of the filter variance-covariance matrix, further accelerating the calculation efficiency of the network solution phase fractional bias; Wherein, is the variance of the observations of the i-th satellite at the k-th epoch, and respectively represent the predicted value and the filtered value of the Kalman filter state, and are the variance matrices corresponding to the states respectively, I represents the identity matrix, and k k,i represents the Kalman filter gain matrix, and l k,i represents the observation of the i-th satellite at the k-th epoch, and a k,i is the row vector corresponding to the observation of the i-th satellite at this station in the matrix in formula (6); Step 36: Based on the wide-lane ambiguity estimation value deviation of 0.15 cycles, the standard deviation of 0.15 cycles, and the probability of 1000 that the ambiguity is fixed as an integer, if the wide-lane ambiguity obtained by filtering a certain satellite according to Step 35 meets the above requirements, then fix it as an integer. Traverse all the estimated wide-lane ambiguities until no new wide-lane ambiguity can be fixed, and output the fractional deviation of the satellite wide-lane phase at this time Step 4, taking the ionosphere-free combined floating ambiguities of multiple stations obtained in Step 2 as input observations, and perform fast estimation of time-varying narrow-lane phase fractional bias; Step 5, linearly transform the fractional biases of each frequency phase.

2. The fast estimation method for real-time Beidou phase fractional bias based on parallel computing according to claim 1, characterized in that, in Step 1, the data stream time synchronization and validity verification specifically includes the following steps: Step 11, obtain the local time, set the phase fractional bias update rate and the data stream waiting time, and use this as the time reference for obtaining synchronous real-time data streams epoch by epoch; Step 12, use the broadcast ephemeris and real-time orbit / clock correction data to restore the real-time satellite orbit and clock data, and monitor whether there are missing and abnormal data; Step 13, detect and update multi-system differential code bias DCB files, predicted Earth orientation parameter EOP files, antenna files, and ocean tide files in BLQ format as external input data for subsequent data preprocessing; Step 14, aiming at the problem of missing observations at some stations, where there are continuous epochs with interrupted observation data, set the waiting time, and the data of this station will not participate in the subsequent phase fractional bias calculation during the waiting time.

3. The fast estimation method for real-time Beidou phase fractional bias based on parallel computing according to claim 1, characterized in that, in Step 2, the parallel computing of floating ambiguities at multiple stations specifically includes the following steps: Step 21, delete the satellites with missing dual-frequency pseudorange and phase observations, and form the following dual-frequency linear combination observations for each satellite: where f 1 and f 2 represent the first frequency and the second frequency of BDS respectively, P 1 / P 2 and L 1 / L 2 are the pseudorange observation values and phase observation values of the first frequency and the second frequency of BDS respectively, PI is the difference between the pseudorange observation values of the first frequency and the second frequency, LI is the difference between the phase observation values of the first frequency and the second frequency. The MW combined observation value is constructed from the dual-frequency carrier phase L 1 / L 2 and the pseudorange observation value P 1 / P 2 . PC is the dual-frequency ionosphere-free pseudorange combined observation value, and LC is the dual-frequency ionosphere-free phase combined observation value; Step 22: Detect pseudo-range gross errors by considering the orbital altitudes and PI value thresholds of GEO / IGSO / MEO satellites, detect cycle slips using the MW combination and the LI combination. Considering the near-geostationary characteristics of GEO satellites, the ionospheric change rate threshold for detecting cycle slips using the LI linear combination is set to half of that of MEO / IGSO. Reject the observation data exceeding the threshold, and record the number of arc segments in each continuous observation arc segment and the switching time between the old and new arc segments; Step 23, correct various error models, including tropospheric delay, antenna phase center correction, phase wrapping, tidal effect, satellite multipath spatial error, where the GEO satellite attitude adopts a zero-bias mode, the IGSO / MEO satellite adopts a dynamic bias mode, and a moving sliding window is used to calculate the MW average value m for each epoch mw and its variance σ mw : where k represents the window size within the continuous arc segment, and the maximum moving window size is set to 120; Step 24. Set the prior variances of the observations of GEO / IGSO / MEO satellites according to formula (3). Since the environments of the receiver and satellites in different orbits are different, the empirical variance values of the pseudorange and phase observations of the three types of satellites are set to be 3 / 0.3 / 0.3 m and 0.03 / 0.003 / 0.003 m respectively, where e represents the satellite elevation angle; 0 ​ Step 25: Using the PC and LC after correcting various model errors in Step 23, adopting the ionosphere-free combined PPP model as in Formula (4), and using the IGG-III weight function model, judge and reduce the weights of abnormal observations according to the standardized residuals, so as to obtain a reliable ionosphere-free combined float ambiguity A if ; where dt r represents the receiver clock error of station r, Z r and m represent the wet tropospheric delay and the wet delay projection function at the zenith of the station respectively, ε PC and ε LC represent the unmodeled errors of the PC and LC observations respectively; Step 26: When there is observation data from multiple stations, follow the steps of Step 21 - Step 25 and use OpenMP for parallel processing to obtain the synchronized MW average values of multiple stations and the ionosphere-free combined floating ambiguities.

4. The real-time Beidou phase fractional bias rapid estimation method based on parallel computing as claimed in claim 1, characterized in that, in Step 4, the rapid estimation of the time-varying narrow-lane phase fractional bias specifically includes the following steps: Step 41: Taking the fixed wide-lane integer ambiguity in Step 3 as a known value, correct the ionosphere-free floating ambiguity obtained in Step 2 for each satellite, as shown in formula (8): In the formula, is the floating-point value of the ambiguity of the corrected narrow lane, λ w and λ n are the wide-lane and narrow-lane wavelengths of BDS. Taking as the input observation value, taking the narrow-lane integer ambiguity and the fractional deviations of the satellite and receiver narrow-lane phases as the estimation parameters, the design matrix of the network solution narrow-lane phase fractional deviation function model is the same as formula (6); Step 42: Refer to the steps in Step 3 to make the narrow-lane phase fractional deviation estimation consistent with the wide-lane phase fractional deviation estimation, which is convenient for program extension and use, and output the satellite narrow-lane phase fractional deviation at this time 5. The real-time Beidou phase fractional bias rapid estimation method based on parallel computing as claimed in claim 1, characterized in that, in Step 5, the linear conversion of the fractional biases of each frequency phase specifically includes the following steps: Using the wide-lane and narrow-lane fractional biases of the phase obtained in Step 3 and Step 4, according to the coefficients of the corresponding frequencies of the MW combination and the LC combination in formula (1), convert them to the fractional biases of each frequency phase according to formula (9) and formula (10) and mark the pseudo-range observation value channels, so that users can use a flexible function model to achieve wide-area single-station ambiguity fixing; Wherein, and are respectively the decimal deviation of the first frequency and the second frequency phase after conversion, that is, the real-time BDS phase decimal deviation, and DCB is the satellite pseudorange differential code delay on the corresponding pseudorange observation value channel.

Images