A real-time ionospheric robust extraction method based on PPP-B2b and BIE weighted estimation in complex environment

Abstract

The present application relates to a kind of real-time ionosphere robust extraction method under complex environment based on PPP-B2b and BIE weighted estimation, comprising: collecting broadcast ephemeris and SSR information, the orbit, clock correction number in broadcast ephemeris and SSR information are matched with data identification, for estimating and broadcast stable satellite end uncalibrated phase delay product, for constructing dual-frequency non-combination PPP observation model, obtain the variance-covariance matrix of float ambiguity vector in state vector;Using the variance-covariance matrix of float ambiguity vector searches candidate vector in the whole space, calculates the posterior probability weight of each candidate solution in candidate vector to observation likelihood, for weighted fusion estimation to ionosphere parameter, obtains robust ionosphere slant delay;Robust ionosphere slant delay is mapped into the vertical total electron content at puncture point, generates the high-precision VTEC grid product that entire service area is covered and time-space continuous.

Description

Technical Field

[0001] This invention relates to the fields of high-precision data processing for satellite navigation and space environment detection, and in particular to a robust real-time ionospheric extraction method for complex environments based on PPP-B2b and BIE weighted estimation. Background Technology

[0002] Real-time, high-precision Total Electron Content (TEC) information plays an irreplaceable and crucial role in improving the accuracy of single-frequency Precise Point Positioning (PPP) of Global Navigation Satellite System (GNSS) from meters to decimeters and even centimeters, as well as in applications such as real-time space weather monitoring. Currently, the technical approach for generating such ionospheric products from GNSS reference station network data is relatively clear: First, an uncombined (UC) PPP model is constructed using dual-frequency or multi-frequency observations from each reference station; then, the geometric distance between the satellite and receiver is obtained with centimeter- to millimeter-level accuracy through fixed carrier phase integer ambiguity resolution (PPP-AR); finally, this high-precision geometric distance is used to separate the ionospheric slant TEC (STEC) along the signal propagation path from the observation equations, and further, regional or global Vertical TEC (VTEC) grid products are generated through spatial modeling.

[0003] In this technological chain, rapid and reliable ambiguity fixation is a prerequisite and core bottleneck for obtaining high-precision ionospheric information. Currently, PPP-AR generally adopts Integer Least Squares (ILS) estimation theory and relies on algorithms such as LAMBDA to achieve efficient search and fixation. This method is relatively effective in open-sky environments with good observation conditions. However, in complex observation environments such as urban canyons and dense forests, severe multipath effects and frequent signal obstruction will significantly degrade the carrier phase observation quality, leading to a sharp increase in the uncertainty of the floating-point ambiguity solution in the non-combined PPP model, and its probability distribution tends to diffuse. At this time, the "either / or" fixation strategy followed by traditional ILS faces a fundamental challenge: if ambiguity reliability is forcibly fixed through statistical testing, it is very easy to cause incorrect fixation, which will lead to undetectable systematic gross errors in the inverted ionospheric delay, causing the product to fail completely; if a less accurate floating-point solution is forced to be used due to the failure of the test, the extracted ionospheric delay will contain huge noise and cannot meet the requirements of high-precision applications.

[0004] To alleviate the difficulties of ambiguity fixation in complex environments, strategies such as Partial Ambiguity Resolution (PAR) have been proposed. These strategies improve the success rate of fixing subset ambiguities by sacrificing observations from some low-quality satellites. However, this method comes at the cost of reduced spatial geometric strength and ionospheric puncture point sampling density, which lowers the accuracy and spatial resolution of the final product. It is a compromise rather than a fundamental solution.

[0005] In recent years, the Best Integer Equivariant (BIE) estimation theory has been proven to be superior to traditional fixed integer and floating-point solutions in terms of mean square error when the uncertainty of the floating-point ambiguity solution is large. This is achieved by probability-weighting all possible integer candidate solutions. Currently, BIE research and applications primarily focus on directly improving the convergence speed and accuracy of the positioning results themselves, with its output considered a better solution to the positioning problem. Meanwhile, the PPP-B2b service provided by China's BeiDou-3 global satellite navigation system broadcasts real-time precise orbit and clock corrections via geostationary orbit (GEO) satellites, providing a robust data foundation for users to achieve real-time precise point positioning without relying on the internet.

[0006] However, existing technologies have yet to provide a solution for how to deeply integrate the real-time, precise observation conditions provided by PPP-B2b services with the statistically optimal characteristics of the BIE algorithm in handling ambiguity and uncertainty, and how to leverage the advantages of this integration to pre-emptively guide the highly robust and continuous estimation process of the ionosphere, a key space environment parameter. Specifically, in complex environments, designing effective algorithmic mechanisms to utilize the ambiguity probability distribution information output by BIE to ensure that the ionospheric parameter extraction process does not experience a "performance cliff" due to ambiguity uncertainty, and achieving robustness throughout the entire process from data acquisition to product generation, remains a critical technical challenge that urgently needs to be overcome. Summary of the Invention

[0007] To address the problems existing in the prior art, the purpose of this invention is to provide a robust real-time ionospheric extraction method for complex environments based on PPP-B2b and BIE weighted estimation. The ambiguity integer candidate solutions and their posterior probability weights output by the BIE algorithm are innovatively injected into the state estimation process of ionospheric parameters in the non-combined PPP model as a new "probabilistic virtual observation constraint". This directly transforms the statistically optimal characteristics of handling ambiguity uncertainty into an improvement in the accuracy and robustness of ionospheric parameter estimation.

[0008] To achieve the above objectives, the present invention provides the following solution: A robust real-time ionospheric extraction method for complex environments based on PPP-B2b and BIE weighted estimation includes: Collect broadcast ephemeris data from GNSS satellites and SSR information from PPP-B2b signals in BeiDou GEO satellites, and perform data identification matching on the orbit and clock error correction values ​​in the broadcast ephemeris and SSR information to recover the precise orbit and clock error; Based on the precise orbit and clock error estimation, a stable uncalibrated phase delay product from the satellite end is broadcast to construct a dual-frequency non-combined PPP observation model and obtain the variance-covariance matrix of the floating-point ambiguity vector in the state vector. Candidate vectors are searched throughout the space using the variance-covariance matrix of the floating-point ambiguity vector. The posterior probability weight of each candidate solution in the candidate vector with respect to the observation likelihood is calculated and used to perform weighted fusion estimation of ionospheric parameters to obtain robust ionospheric slant delay. The robust ionospheric slant delay is mapped to the vertical total electron content at the puncture point to generate a high-precision VTEC grid product that covers the entire service area and is spatiotemporally continuous.

[0009] Optionally, restoring the precision orbit and clock error includes: Calculate the orbital correction in the star-fixed coordinate system at the current moment, and define the unit vector of the star-fixed coordinate system to convert the orbital correction to the geocentric geofixed coordinate system, correct the satellite orbit of the broadcast ephemeris, and obtain the precise orbit; Calculate the satellite clock error correction for the current moment, and use the satellite clock error correction to correct the satellite clock error of the broadcast ephemeris to obtain the precise clock error.

[0010] Optionally, estimating and broadcasting stable uncalibrated phase delay products at the satellite end based on the precise orbit and clock bias includes: The server performs static PPP model calculations on the observation data of the precise orbit and clock error and the regional reference station network to obtain the floating-point ambiguity of each satellite in each station; Spatial domain adjustment is performed on the floating-point ambiguities of all regional reference stations under the same satellite to estimate and broadcast the uncalibrated phase delay product at the satellite end.

[0011] Optionally, constructing the dual-frequency non-combined PPP observation model includes: ; ; in, This represents the pseudorange observation of satellite s at frequency j by receiver r. This represents the carrier phase observation of satellite s at frequency j by receiver r. The distance between satellite s and receiver r represents the geometric distance, and c represents the speed of light in a vacuum. Represents the receiver clock bias, which includes receiver code hardware delay. The precise clock bias representing the recovered satellite s, The ionospheric coefficient representing frequency j, The ionospheric slant delay represents the signal path delay between satellite s and receiver r at the first frequency. The tropospheric delay represents the signal path delay between satellite s and receiver r. Noise and multipath error representing pseudorange observations, The carrier wavelength representing frequency j, This represents the ambiguity parameter obtained by receiver r after recovering the integer characteristics of satellite s at frequency j. This represents the phase hardware delay of receiver r at frequency j. This represents the phase hardware delay of satellite s at frequency j, i.e., the UPD product estimated and broadcast by the server. This represents the noise and multipath error of the carrier phase observations.

[0012] Optionally, searching for candidate vectors throughout the space using the variance-covariance matrix of the floating-point ambiguity vector includes: Using the floating-point ambiguity vector as the expectation and the variance-covariance matrix of the floating-point ambiguity vector as the metric, candidate vectors are searched throughout the space. During the search for the candidate vectors, the search range is limited by preset search conditions: ; in, Let n be the floating-point ambiguity vector, where n is the number of ambiguity parameters. floating-point ambiguity vector The variance-covariance matrix, Let be the i-th integer candidate vector, where each component is an integer. The threshold is determined based on the search space volume. Given a set containing K candidate integer sets, Let z be an n-dimensional integer space, and let z be the set of all possible integer solutions.

[0013] Optionally, obtaining the robust ionospheric slack delay includes: Calculate the posterior probability weight of each candidate solution in the candidate vector with respect to the observation likelihood: ; in, Representative candidate integer vector The posterior probability weight, i.e., the relative probability that the candidate solution is correct. The square of the Mahalanobis distance is represented by K, and the total number of candidate integer pairs is represented by K. The ionospheric parameters are estimated by weighted fusion estimation using the weighted virtual observation defusion method or the multi-hypothesis Kalman filter method, and the robust ionospheric slant delay is obtained.

[0014] Optionally, the weighted virtual observation defusion method is used to perform weighted fusion estimation of ionospheric parameters using the posterior probability weights, including: Each candidate solution in the candidate vector is substituted as an integer constraint into the original phase observation model to eliminate ambiguity parameters, generate a virtual phase observation sequence, and perform a fast solution on the virtual phase observation sequence to obtain the ionospheric oblique delay estimate. Probability-weighted fusion estimation is performed based on the ionospheric slant delay estimate and the posterior probability weight: ; In the formula, For the final ionospheric slant delay obtained by BIE-weighted estimation of satellite s, For candidate solutions The posterior probability weights, To use candidate solutions The ionospheric oblique delay estimate obtained by virtual observation solution, where K is the total number of candidate solutions.

[0015] Optionally, the weighted fusion estimation of ionospheric parameters using the posterior probability weights via multi-hypothesis Kalman filtering includes: The candidate solution with the highest posterior probability for the number of targets is selected to initialize the parallel filtering hypothesis for the number of targets. The ambiguity state of the parallel filtering hypothesis is fixed, and the initial probability of the parallel filtering hypothesis is set as the posterior probability weight. In each filtering epoch, each parallel filtering hypothesis is independently predicted and corrected for state, used to calculate the innovation observation residuals and innovation covariance for each parallel filtering hypothesis, and dynamically update the probabilities of each parallel filtering hypothesis to reflect the degree of support of the current observations for each hypothesis. ; in, Representative hypothesis Prior probabilities before measurement update Representative hypothesis The posterior probability after measurement update Representative hypothesis The information vector is the difference between the observed value and the predicted value; Representative hypothesis The new information covariance matrix, Represents the determinant of a matrix; The probabilities of each parallel filtering hypothesis, dynamically updated, are used to perform a weighted sum of all hypothetical states for the ionospheric parameters: ; in, This represents the final ionospheric slant delay estimate for the current epoch satellite s. Assumption The ionospheric slant delay state value of satellite s.

[0016] Optionally, generating the high-precision VTEC grid product includes: The robust ionospheric slant delay is converted into the total slant electron content, and the relationship between the total slant electron content and the robust ionospheric slant delay is obtained. Using a monolayer ionospheric model, the relationship between the oblique total electron content and the robust ionospheric oblique delay is mapped to the vertical total electron content at the puncture point: ; Wherein, VTEC represents the vertical total electron content. The relationship between total oblique electron content and robust ionospheric oblique delay. Represents the robust ionospheric slack delay between satellite s and receiver r. This represents the zenith distance of the signal path at the height of a single layer of the ionosphere. The projection function is represented by the satellite elevation angle, azimuth angle, and single-layer height. The carrier frequency is the first frequency. The high-precision VTEC grid product is generated by aggregating the total vertical electron content at the puncture points corresponding to all reference stations within the region via the server.

[0017] Optionally, the method further includes: The high-precision VTEC grid product is broadcast to target users via satellite PPP-B2b link or mobile Internet for single-frequency PPP ionospheric delay correction and space weather monitoring.

[0018] The beneficial effects of this invention are as follows: (1) Achieved “performance continuity” in ionosphere extraction: By introducing the probability weighting mechanism of BIE, this invention completely avoids the performance precipitate that occurs when the ambiguity fixation fails in traditional methods. In complex environments, even if a deterministic ambiguity fixation solution cannot be obtained, the system can still output an ionosphere estimate that is statistically optimal and has significantly better accuracy than a pure floating-point solution, ensuring the continuity and availability of product output.

[0019] (2) Significantly improves robustness of extraction under complex environments: Traditional ILS fixation is susceptible to fluctuations in observation quality, leading to incorrect fixation and gross errors in ionospheric solutions. The multi-hypothesis weighted or filtering method of this invention effectively suppresses catastrophic bias caused by a single erroneous integer solution by comprehensively evaluating all reasonable ambiguity possibilities, thus greatly enhancing the robustness of the ionospheric extraction process to observational anomalies and model errors.

[0020] (3) Fully leverage the synergistic advantages of PPP-B2b and BIE: This invention does not simply use PPP-B2b and BIE side by side, but creatively uses the intermediate information of "ambiguity probability distribution" output by BIE as the key input for optimizing ionospheric parameter estimation, forming an innovative technical path of precise observation conditions, probability-optimal ambiguity processing and highly robust ionospheric extraction, and realizing the synergistic effect of "1+1>2".

[0021] (4) Ensures the accuracy and reliability of the final product: Since the BIE estimation is optimal in the sense of mean square error, and this method ensures that the ionospheric parameters can make full use of this optimality, the final real-time ionospheric grid product has higher accuracy and reliability in complex environments, and can more effectively serve high-requirement applications such as single-frequency PPP accuracy improvement and real-time space weather monitoring. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a flowchart of a robust real-time ionosphere extraction method in complex environments based on PPP-B2b and BIE weighted estimation, according to an embodiment of the present invention. Figure 2 This is a system architecture and data flow diagram according to an embodiment of the present invention; Figure 3 This is a schematic diagram illustrating the robust ionospheric extraction principle based on BIE probability weighting in an embodiment of the present invention. Figure 4 This is a flowchart of the weighted fusion estimation process for ionospheric parameters according to an embodiment of the present invention. Detailed Implementation

[0024] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0025] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0026] like Figure 1 , 2 As shown, this embodiment discloses a robust real-time ionospheric extraction method for complex environments based on PPP-B2b and BIE weighted estimation. The method includes: collecting broadcast ephemeris data from GNSS satellites and SSR information from PPP-B2b signals in BeiDou GEO satellites; performing data identification matching on the orbit and clock correction values ​​in the broadcast ephemeris and SSR information to recover precise orbits and clock errors; estimating and broadcasting stable uncalibrated phase delay products from the satellite end based on precise orbits and clock errors to construct a dual-frequency non-combined PPP observation model; obtaining the variance-covariance matrix of the floating-point ambiguity vector in the state vector; using the variance-covariance matrix of the floating-point ambiguity vector to search for candidate vectors throughout the space; calculating the posterior probability weight of each candidate solution in the candidate vector for the observation likelihood; using this weighted fusion estimation of ionospheric parameters to obtain robust ionospheric slant delay; mapping the robust ionospheric slant delay to the vertical total electron content at the puncture point; and generating a high-precision VTEC grid product that covers the entire service area and is spatiotemporally continuous.

[0027] Specifically, this invention proposes a robust real-time ionosphere extraction method for complex environments based on PPP-B2b and BIE weighted estimation, comprising: Step 1: Real-time data reception and IOD matching.

[0028] Step 2: PPP-B2b real-time precision track and clock error recovery.

[0029] Step 3: Server-side UPD estimation and client-side non-combined PPP floating-point solution.

[0030] Step 4: Generate ambiguity probability candidate set based on BIE.

[0031] Step 5: Robust estimation of ionospheric parameters based on probability weighting.

[0032] Step 6: Conversion of ionospheric oblique delay to vertical total electron content.

[0033] Step 7: Generation and broadcasting of wide-area real-time ionospheric products.

[0034] Its specific contents include: Step 1: The user terminal (or server reference station) receives multi-frequency (at least dual-frequency) pseudorange observations broadcast by BeiDou and other GNSS satellites. Carrier phase observations In addition to broadcasting ephemeris; at the same time, it receives PPP-B2b signals broadcast via BeiDou GEO satellites, analyzes the real-time state space representation (SSR) information contained therein, and then performs data identification (IOD) matching on the orbit and clock error corrections in the broadcast ephemeris and SSR information to ensure the consistency of the reference benchmark and lay the foundation for subsequent precision product recovery.

[0035] Step 2: Based on the successfully matched SSR corrections and broadcast ephemeris, restore the satellite orbit and clock error with centimeter-level accuracy in real time.

[0036] (1) Precision track restoration: The formula for calculating the orbital correction at the current time t in the star-fixed coordinate system is as follows: ; In the formula, Representing the current moment, Represents a reference time; Represents the orbital correction vector of the satellite in the star-fixed coordinate system at the current time t; Represents the radial, tangential, and normal corrections of the satellite's orbit at the current time t, where ; Represents the reference time given in the SSR information. The orbital correction vector; Represents the reference time given in the SSR information. The vector of the rate of change of the orbital correction.

[0037] Construct the transformation matrix from the star-fixed coordinate system to the geocentric Earth-fixed coordinate system. Using the satellite positions and velocities calculated from the broadcast ephemeris, define the unit vector for the star-fixed coordinate system, as shown in the following formula: ; ; ; In the formula, The satellite position represents the value calculated from the broadcast ephemeris. The satellite velocity represented by the broadcast ephemeris calculation; represents the unit vector in the star-fixed coordinate system, where .

[0038] The formula for converting the orbital corrections to the Earth-centered Earth-fixed coordinate system is as follows: ; In the formula, Represents the orbital correction vector of the satellite in the geocentric-ground-fixed coordinate system at the current time t; represents the unit vector in the star-fixed coordinate system, where ; This represents the orbital correction vector of the satellite in the star-fixed coordinate system at the current time t.

[0039] Correct the satellite's orbit Corrected to satellite orbits calculated from broadcast ephemeris The formula for obtaining real-time precise orbits is as follows: ; In the formula, Represents the precise orbital position of the satellite (ECEF coordinate system) after correction using PPP-B2b service at the current time t. This represents the satellite's orbital position calculated from the broadcast ephemeris.

[0040] (2) Precision clock error recovery: The satellite clock error correction at the current time is calculated using the SSR clock error correction polynomial model, and the formula is as follows: ; In the formula, The clock error correction value represents the current time. Representing the current moment, Reference time representing the clock error correction; The polynomial coefficients representing the clock error correction, where .

[0041] Correct the satellite's clock bias The satellite clock bias is corrected to the broadcast ephemeris calculation to obtain the real-time precise clock bias, and the formula is as follows: ; In the formula, The precise clock difference of the satellite at the current time t, This represents the satellite clock bias calculated from the broadcast ephemeris at the current time t. Represents the speed of light in vacuum.

[0042] Step 3: This step consists of two parallel and collaborative processes: the server-side and the user-side.

[0043] (1) Server-side UPD estimation: The server uses observation data from the regional reference station network (coordinates known) and the precise product recovered in step two to calculate the floating-point ambiguity of each station and each satellite using a static PPP model. By comparing all reference stations of the same satellite Perform spatial domain adjustment to estimate and broadcast stable satellite-end uncalibrated phase delay (UPD) products. The server-side observation equations are the same as those for the user side (see below), but since the coordinates of the reference station are precisely known, the coordinates can be used as fixed values, thereby more accurately estimating the UPD at the satellite end.

[0044] (2) User-side non-combined PPP floating-point solution: At the receiver, using the recovered precise orbit and clock error, a dual-frequency uncombined PPP observation equation is constructed. For satellite s and frequency j, the linearized observation equation is: ; ; In the formula, This represents the pseudorange observation of satellite s at frequency j by receiver r; This represents the carrier phase observation of satellite s at frequency j by receiver r; The distance between satellite s and receiver r represents the geometric distance; c represents the speed of light in a vacuum. Represents the receiver clock bias, which includes receiver code hardware delay; Represents the precise clock error of satellite s recovered from step two; The ionospheric coefficient representing frequency j, ,in The carrier frequency is the first frequency. The ionospheric slant delay at the first frequency represents the signal path delay between satellite s and receiver r. The tropospheric delay represents the signal path between satellite s and receiver r; This represents noise and multipath error in pseudorange observations; The carrier wavelength representing frequency j; This represents the floating-point ambiguity of receiver r for satellite s at frequency j, including the integer part and hardware delay; This represents the phase hardware delay of receiver r at frequency j; This represents the phase hardware delay of satellite s at frequency j, i.e., the UPD product estimated and broadcast by the server. This represents the noise and multipath error of the carrier phase observations.

[0045] UPD products broadcast by the application server Correcting the phase observations restores the integer properties of the ambiguity parameters. The corrected phase observation equation is: ; In the formula, This is the ambiguity parameter recovered from the integer characteristics (still a floating-point number, but with the integer part being an integer). For simplicity, the receiver phase delay is typically... The new clock bias parameters are absorbed into the receiver clock bias parameters. .

[0046] Finally, the state vector Includes receiver location and receiver clock bias. Tropospheric delay, ionospheric slant delay of each satellite and floating-point ambiguity vector (Where m is the number of satellites). The state estimate is obtained by solving the problem using a Kalman filter. and its variance-covariance matrix, in particular, the floating-point ambiguity vector. The variance-covariance matrix is ​​denoted as .

[0047] Step 4: This step abandons the traditional single-integer fixed strategy and instead adopts the Optimal Integer Equivariant (BIE) estimation framework to generate a weighted set of integer candidate solutions to quantify the uncertainty of ambiguity.

[0048] (1) Searching for candidate integer sets: using floating-point solutions For the expectation, For measurement, in integer space Search for the K highest probability integer candidate vectors So that they satisfy: ; In the formula, This is the floating-point ambiguity vector obtained in step three (with dimension n, where n is the number of ambiguity parameters). floating-point ambiguity vector The variance-covariance matrix; Let be the i-th integer candidate vector (with dimension n), where each component is an integer; The threshold is determined based on the search space volume; Let K be a set containing K candidate integer sets.

[0049] (2) Calculate the posterior probability weights: for each candidate solution Calculate its posterior probability weights based on observation likelihood. : ; In the formula, Representative candidate integer vector The posterior probability weight is a scalar between 0 and 1, representing the relative probability that the candidate solution is correct. The square of the Mahalanobis distance is defined as follows: K represents the total number of candidate integer groups.

[0050] This probability weight set Instead of being used to select a fixed "winning" solution, it is used as a complete probability distribution describing the uncertainty of ambiguity and is directly input into subsequent ionospheric parameter estimation.

[0051] Step 5: Using the probability weights generated in Step 4, perform weighted fusion estimation of ionospheric parameters to achieve a continuous performance transition, such as... Figure 3 As shown. Two preferred implementation methods are provided: (1) Method 1: Weighted virtual observation solution fusion method: Construct virtual observations for each candidate solution. This is treated as a definite integer constraint, substituted into the original phase observation equation, and the ambiguity parameter is eliminated to generate a set of corresponding "virtual" phase observation sequences. Specifically, for each candidate solution... Subtract a fixed ambiguity from the phase observations: ; In the formula, It is a candidate solution The integer value (in weeks) corresponding to the frequency j of satellite s. Thus, the virtual observation value... The corresponding observation equation no longer contains the ambiguity parameter.

[0052] Independent calculations are performed, and a simplified Kalman filter (or least squares) is re-run based on each set of virtual observations to quickly obtain a set of corresponding ionospheric oblique delay estimates. At this point, the state vector no longer contains ambiguity parameters, but only receiver position, clock error, tropospheric delay, and ionospheric delay.

[0053] Probabilistic weighted fusion yields a final robust ionosphere extraction result that is a probabilistically weighted average of the results corresponding to each candidate solution: ; In the formula, This is the final ionospheric slant delay obtained for satellite s based on BIE weighted estimation; It is a candidate solution The posterior probability weights; Is it using candidate solutions? The ionospheric oblique delay estimate obtained from virtual observation calculations; K is the total number of candidate solutions.

[0054] (2) Method 2: Multiple Hypothesis Kalman Filtering: like Figure 4 As shown, multiple hypotheses are initialized, and the one with the highest posterior probability is selected. candidate solutions Initialize M parallel filtering hypotheses This parallel filtering assumption This is derived from the weighted candidate solution set generated by the BIE algorithm in step S4. Each hypothesis... Having an independent state vector And fix its ambiguity state as The initial probabilities of each hypothesis are: (After normalization, the probability weights of the first M candidates are considered and normalized again.)

[0055] Parallel time and measurement updates are implemented. At each filtered epoch, all hypotheses independently perform state prediction (time update) and state correction based on the original observations (measurement update). For each hypothesis... Calculate its new information (observation residuals). and new information covariance .

[0056] Dynamic probability updates dynamically adjust the probabilities of each hypothesis based on the magnitude of the innovation, reflecting the degree of support the current observations provide for each hypothesis. The update formula is: ; In the formula, Representative hypothesis Prior probabilities before measurement update; Representative hypothesis The posterior probability after measurement update; Representative hypothesis The information vector (the difference between the observed value and the predicted value); Representative hypothesis The new information covariance matrix; Represents the determinant of the matrix. After updating... Normalization is performed so that .

[0057] Weighted state output, final state estimate for each epoch (especially ionospheric state) The sum of all assumed states is: ; Specifically, the ionospheric parameters include: ; In the formula, The final ionospheric slant delay estimate represents the current epoch satellite s; Assumption The posterior probability at the current epoch; Assumption The ionospheric slant delay state value of satellite s.

[0058] Step Six: Extract the robust ionospheric slack delay data in step five, in units of length. The ionospheric monolayer model (SLM) is used to map the vertical total electron content (VTEC) at the puncture point. This conversion involves the following two steps: (1) Delay the slope (meters) converted to total oblique electron content (STEC): According to the physical definition of ionospheric delay, its relationship with STEC is as follows: ; In the formula, The first frequency is the carrier frequency (unit: Hz). The unit of STEC is TECU (1 TECU = 10). 16 (electronics / square meter).

[0059] (2) Using SLM to map STEC to vertical total electron content (VTEC): ; In the formula, VTEC represents the total vertical electron content (unit: TECU). This represents the robust ionospheric skew delay between satellite s and receiver r obtained in step S5; This represents the zenith distance of the signal path at the height of a single layer of the ionosphere; The projection function is represented by the satellite elevation angle, azimuth angle, and single-layer height. , The average radius of the Earth The elevation angle of the satellite relative to the receiver.

[0060] Step 7: The server aggregates real-time VTEC data for puncture points generated by all reference stations within the region. Using spatial interpolation algorithms such as spherical harmonic function models or inverse distance weighting, it generates a spatiotemporally continuous, real-time, high-precision VTEC grid product covering the entire service area. This product can be broadcast to a massive number of users via satellite PPP-B2b links or mobile internet for applications such as single-frequency PPP ionospheric delay correction and space weather monitoring.

[0061] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A robust real-time ionospheric extraction method for complex environments based on PPP-B2b and BIE weighted estimation, characterized in that, include: Collect broadcast ephemeris data from GNSS satellites and SSR information from PPP-B2b signals in BeiDou GEO satellites, and perform data identification matching on the orbit and clock error correction values ​​in the broadcast ephemeris and SSR information to recover the precise orbit and clock error; Based on the precise orbit and clock error estimation, a stable uncalibrated phase delay product from the satellite end is broadcast to construct a dual-frequency non-combined PPP observation model and obtain the variance-covariance matrix of the floating-point ambiguity vector in the state vector. Candidate vectors are searched throughout the space using the variance-covariance matrix of the floating-point ambiguity vector. The posterior probability weight of each candidate solution in the candidate vector with respect to the observation likelihood is calculated and used to perform weighted fusion estimation of ionospheric parameters to obtain robust ionospheric slant delay. The robust ionospheric slant delay is mapped to the vertical total electron content at the puncture point to generate a high-precision VTEC grid product that covers the entire service area and is spatiotemporally continuous.

2. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, Restoring the precision orbit and clock error includes: Calculate the orbital correction in the star-fixed coordinate system at the current moment, and define the unit vector of the star-fixed coordinate system to convert the orbital correction to the geocentric geofixed coordinate system, correct the satellite orbit of the broadcast ephemeris, and obtain the precise orbit; Calculate the satellite clock error correction for the current moment, and use the satellite clock error correction to correct the satellite clock error of the broadcast ephemeris to obtain the precise clock error.

3. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, Based on the aforementioned precise orbit and clock bias estimation, stable uncalibrated phase delay products for satellite delivery are included: The server performs static PPP model calculations on the observation data of the precise orbit and clock error and the regional reference station network to obtain the floating-point ambiguity of each satellite in each station; Spatial domain adjustment is performed on the floating-point ambiguities of all regional reference stations under the same satellite to estimate and broadcast the uncalibrated phase delay product at the satellite end.

4. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, The construction of the dual-frequency non-combined PPP observation model includes: ； ； in, This represents the pseudorange observation of satellite s at frequency j by receiver r. This represents the carrier phase observation of satellite s at frequency j by receiver r. The distance between satellite s and receiver r represents the geometric distance, and c represents the speed of light in a vacuum. Represents the receiver clock bias, which includes receiver code hardware delay. The precise clock bias representing the recovered satellite s, The ionospheric coefficient representing frequency j, The ionospheric slant delay represents the signal path delay between satellite s and receiver r at the first frequency. The tropospheric delay represents the signal path delay between satellite s and receiver r. Noise and multipath error representing pseudorange observations, The carrier wavelength representing frequency j, This represents the ambiguity parameter obtained by receiver r after recovering the integer characteristics of satellite s at frequency j. This represents the phase hardware delay of receiver r at frequency j. This represents the phase hardware delay of satellite s at frequency j, i.e., the UPD product estimated and broadcast by the server. This represents the noise and multipath error of the carrier phase observations.

5. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, Searching for candidate vectors throughout the space using the variance-covariance matrix of the floating-point ambiguity vector includes: Using the floating-point ambiguity vector as the expectation and the variance-covariance matrix of the floating-point ambiguity vector as the metric, candidate vectors are searched throughout the space. During the search for the candidate vectors, the search range is limited by preset search conditions: ； in, Let n be the floating-point ambiguity vector, where n is the number of ambiguity parameters. floating-point ambiguity vector The variance-covariance matrix, Let be the i-th integer candidate vector, where each component is an integer. The threshold is determined based on the search space volume. Given a set containing K candidate integer sets, Let z be an n-dimensional integer space, and let z be the set of all possible integer solutions.

6. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, Obtaining the robust ionospheric slack delay includes: Calculate the posterior probability weight of each candidate solution in the candidate vector with respect to the observation likelihood: ； in, Representative candidate integer vector The posterior probability weight, i.e., the relative probability that the candidate solution is correct. The square of the Mahalanobis distance is represented by K, and the total number of candidate integer pairs is represented by K. The ionospheric parameters are estimated by weighted fusion estimation using the weighted virtual observation defusion method or the multi-hypothesis Kalman filter method, and the robust ionospheric slant delay is obtained.

7. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 6, characterized in that, The weighted virtual observation defusion method utilizes the posterior probability weights to perform weighted fusion estimation of ionospheric parameters, including: Each candidate solution in the candidate vector is substituted as an integer constraint into the original phase observation model to eliminate ambiguity parameters, generate a virtual phase observation sequence, and perform a fast solution on the virtual phase observation sequence to obtain the ionospheric oblique delay estimate. Probability-weighted fusion estimation is performed based on the ionospheric slant delay estimate and the posterior probability weight: ； In the formula, For the final ionospheric slant delay obtained by BIE-weighted estimation of satellite s, For candidate solutions The posterior probability weights, To use candidate solutions The ionospheric oblique delay estimate obtained by virtual observation solution, where K is the total number of candidate solutions.

8. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 6, characterized in that, The multi-hypothesis Kalman filtering method uses the posterior probability weights to perform weighted fusion estimation of ionospheric parameters, including: The candidate solution with the highest posterior probability for the number of targets is selected to initialize the parallel filtering hypothesis for the number of targets. The ambiguity state of the parallel filtering hypothesis is fixed, and the initial probability of the parallel filtering hypothesis is set as the posterior probability weight. In each filtering epoch, each parallel filtering hypothesis is independently predicted and corrected for state, used to calculate the innovation observation residuals and innovation covariance for each parallel filtering hypothesis, and dynamically update the probabilities of each parallel filtering hypothesis to reflect the degree of support of the current observations for each hypothesis. ； in, Representative hypothesis Prior probabilities before measurement update Representative hypothesis The posterior probability after measurement update Representative hypothesis The information vector is the difference between the observed value and the predicted value; Representative hypothesis The new information covariance matrix, Represents the determinant of a matrix; The probabilities of each parallel filtering hypothesis, dynamically updated, are used to perform a weighted sum of all hypothetical states for the ionospheric parameters: ； in, This represents the final ionospheric slant delay estimate for the current epoch satellite s. Assumption The ionospheric slant delay state value of satellite s.

9. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, Generating the high-precision VTEC grid product includes: The robust ionospheric slant delay is converted into the total slant electron content, and the relationship between the total slant electron content and the robust ionospheric slant delay is obtained. Using a monolayer ionospheric model, the relationship between the oblique total electron content and the robust ionospheric oblique delay is mapped to the vertical total electron content at the puncture point: ； Wherein, VTEC represents the vertical total electron content. The relationship between total oblique electron content and robust ionospheric oblique delay. Represents the robust ionospheric slack delay between satellite s and receiver r. This represents the zenith distance of the signal path at the height of a single layer of the ionosphere. The projection function is represented by the satellite elevation angle, azimuth angle, and single-layer height. The carrier frequency is the first frequency. The high-precision VTEC grid product is generated by aggregating the total vertical electron content at the puncture points corresponding to all reference stations within the region via the server.

10. The robust real-time ionosphere extraction method under complex environments based on PPP-B2b and BIE weighted estimation as described in claim 1, characterized in that, The method also includes: The high-precision VTEC grid product is broadcast to target users via satellite PPP-B2b link or mobile Internet for single-frequency PPP ionospheric delay correction and space weather monitoring.

Images