Low earth orbit enhanced SBAS ionospheric monitoring system and grid estimation method

Abstract

The application discloses a low-orbit enhanced SBAS ionospheric monitoring system and a grid estimation method, and belongs to the technical field of satellite navigation. Firstly, the global mobile characteristics and the advantages of integrated navigation of low-orbit satellites are used to compensate for the defects of insufficient performance of the SBAS ground monitoring network in the edge area. Secondly, aiming at the incomplete characteristics of the ionospheric observation of low-orbit satellites, a low-orbit enhanced SBAS grid ionospheric delay and scale factor synchronous estimation method is provided. After data preprocessing and extraction of multi-level vertical ionospheric delay, an optimization model with the minimum grid point ionospheric delay variance as the target is constructed, and the ionospheric delay and scale factor grid synchronous estimation is realized through block processing and iterative optimization of independent variables. The low-orbit enhanced SBAS ionospheric monitoring system and the grid estimation method provided by the application are helpful to guarantee the real-time performance of the SBAS ionospheric grid, improve the accuracy of the SBAS ionospheric grid estimation, and provide the scale factor to meet the demand of future low-orbit satellite navigation ionospheric precise correction.

Description

Technical Field

[0001] This invention belongs to the field of satellite navigation technology, and in particular relates to a low-orbit enhanced SBAS ionospheric monitoring system and a grid estimation method. Background Technology

[0002] Global Navigation Satellite Systems (GNSS), capable of providing continuous navigation services around the clock, free flight routes from the constraints of ground-based stations, significantly improving airspace utilization and operational efficiency in civil aviation. However, the inherent vulnerability of GNSS makes it susceptible to deception and interference, easily leading to navigation risks and safety issues. With the development of satellite navigation technology and the ever-increasing demands of navigation and positioning users, the results of single-point positioning using broadcast ephemeris data can no longer meet the accuracy and integrity requirements of high-precision users for GNSS systems. Therefore, Satellite-Based Augmentation Systems (SBAS) have emerged.

[0003] Wide Area Augmentation System (SBAS) collects GNSS observation information through a ground monitoring network, generates differential corrections and integrity parameters, and broadcasts them to users via geostationary orbit (GEO) satellites. This achieves wide-area differential correction and integrity enhancement for GNSS, further ensuring the accuracy and integrity of aviation users' navigation systems. Since the U.S. Wide Area Augmentation System (WAAS) became operational in 2003, many countries and regions around the world have successively promoted the construction of SBAS systems, aiming to provide navigation service capabilities for LPV-200 and CAT-1 aircraft to their own and surrounding areas. Currently, SBAS primarily provides single-frequency (SF) augmentation services, while dual-frequency multi-constellation (DFMC) augmentation services are under planning and testing.

[0004] SF-SBAS enhances accuracy and integrity by broadcasting ephemeris and clock corrections, ionospheric corrections, and corresponding integrity parameters to correct user ephemeris and clock errors and ionospheric delays. Compared to ephemeris and clock corrections, ionospheric corrections play a decisive role in SF-SBAS enhancement services. SBAS abstracts the diffuse three-dimensional ionosphere in space as a two-dimensional thin layer where all electrons are concentrated at a certain height, and divides this thin layer into uniformly defined ionospheric grid points (IGPs) according to certain intervals. The ground monitoring network observes the delay at the ionospheric puncture point (IPP) in real time, estimates the ionospheric delay of the IGP, and transmits the validity of the IGP and the grid ionospheric vertical delay (GIVD) to users via message types 18 and 26, thereby achieving real-time ionospheric correction at the user end.

[0005] However, due to various limitations, ground-based monitoring networks can only be deployed in limited areas over land, resulting in a severely uneven distribution of IPPs (Integrated Precipitation Points). IPPs are densely distributed within the network's coverage area, allowing for effective estimation of IGPs located within it; however, IPPs are sparsely distributed outside the network's coverage area, making it difficult to guarantee the performance of IGPs at the edge of the SBAS service area, and even leading to frequent unavailability. The rapidly developing low-Earth orbit (LEO) constellations will serve as an effective technical means to compensate for the deficiencies of ground-based monitoring networks. On the one hand, compared to medium-Earth orbit (MEO) or high-Earth orbit (HEO) GNSS satellites, LEO satellites have faster movement speeds, generating more and denser IPPs under the short-term static assumption of the ionosphere. On the other hand, due to their unique orbital altitude, LEO satellites are located in the middle of the ionosphere, enabling them to extract information from the upper ionosphere through their onboard GNSS receivers and from the lower ionosphere through their downlink dual-frequency navigation signals.

[0006] While low-Earth orbit (LEO) constellations provide extremely rich ionospheric observation information, effectively utilizing this information to build higher-precision ionospheric correction models remains a challenge. Unlike GNSS signals that traverse the entire ionosphere, LEO ionospheric observations only penetrate a portion of it. GNSS signals received by LEO satellites and LEO navigation signals received from the ground can only reflect partial characteristics of the ionosphere. One solution is to use tomography to synthesize LEO and GNSS ionospheric observation information to construct a three-dimensional model of ionospheric electron density, and then extract the vertical delay at the IGP (Integrated Geometric Point). However, tomography requires finely divided three-dimensional ionospheric voxels, consuming significant computational resources during modeling, making it unsuitable for the high real-time requirements of SF-SBAS services. Another solution is to obtain a scaling factor between the delays of the upper and lower ionospheric layers using an empirical ionospheric model, then reconstruct the incomplete LEO ionospheric observations and fuse them with GNSS ionospheric observations to establish an ionospheric grid model. However, empirical models are insufficient to accurately describe the complex ionospheric variation patterns, and low-Earth orbit ionospheric observations based on scaling factor recovery have large errors, thus affecting the accuracy of ionospheric grid models. Summary of the Invention

[0007] The purpose of this invention is to provide a low-orbit enhanced SBAS ionospheric monitoring system and grid estimation method to solve the problems of insufficient SBAS ground monitoring network coverage and poor monitoring performance in the service edge area in the above-mentioned technologies.

[0008] To achieve the above objectives, the present invention provides a low-Earth orbit enhanced SBAS ionospheric monitoring system, comprising a low-Earth orbit satellite-based monitoring platform, a ground monitoring station, and an information processing center;

[0009] The low-Earth orbit satellite-based monitoring platform relies on low-Earth orbit satellites and is equipped with an onboard GNSS receiver and a downlink navigation signal generation module. It has the ability to receive satellite-based GNSS navigation signals and broadcast independent dual-frequency navigation signals to the ground.

[0010] The ground monitoring station is equipped with multiple independent navigation receivers to receive and analyze GNSS observation information carried in GNSS navigation signals, low-orbit satellite navigation signals, and low-orbit satellite communication signals. It also has Internet access capabilities and aggregates various observation information to the information processing center for comprehensive analysis and processing.

[0011] The information processing center is equipped with high-performance computing equipment to receive and preprocess observation data from ground monitoring stations and low-orbit satellites, and to run a method for synchronous estimation of ionospheric delay and scaling factor of the low-orbit enhanced SBAS grid, generating ionospheric delay corrections and corresponding integrity parameters for SBAS.

[0012] A method for estimating the ionospheric grid of LEO-enhanced SBAS (Balanced Subsurface Assay System) utilizes three types of navigation observation data: LEO satellite and GNSS, ground station and LEO satellite, and ground station and GNSS, to estimate the ionospheric delay and scaling factor of the LEO-enhanced SBAS grid. The specific process for estimating the ionospheric delay and scaling factor of the LEO-enhanced SBAS grid is as follows:

[0013] S1. Construct geometry-free observables; construct the corresponding geometry-free observables using dual-frequency pseudorange and carrier phase observations from navigation data;

[0014] S2. Cycle slip detection: After acquiring geometrically unaffected observations, cycle slip detection is performed on the carrier phase observations to eliminate those affected by cycle slips.

[0015] S3. Filtering and smoothing: The geometrically uncarrier observations and the corresponding geometrically unpseudorange observations obtained through cycle slip detection are filtered and smoothed to improve the accuracy of tilted ionospheric delay extraction.

[0016] S4. Tilted ionospheric delay extraction: Hardware bias in the smoothing results is corrected using differential code bias products to further obtain the observed values ​​of the tilted ionospheric delay and the corresponding measurement error variance.

[0017] S5. Calculation of ionospheric puncture point location; Construct a thin double-layer ionospheric model, and calculate the puncture point location of each tilted ionospheric observation on the thin double-layer ionospheric model based on the coordinates of the ground monitoring station and GNSS / Low Orbit satellite ephemeris data;

[0018] S6. Inclined / Vertical Ionospheric Delay Conversion: Based on the assumption of ionospheric isotropy, the tilted ionospheric delay from different sources is converted into vertical ionospheric delay using a projection function.

[0019] S7. Ionospheric puncture point search; Search for ionospheric puncture points within the vicinity of the grid point to be estimated in the ionospheric region;

[0020] S8. Construct a grid synchronization estimation model for ionospheric delay and scaling factor; comprehensively consider the spatial correlation of the ionosphere, measurement error characteristics, and the relative positional relationship between grid points and puncture points, construct an optimization model with minimizing the variance of the vertical delay estimation error of grid points as the objective function and weighting factor and scaling factor as independent variables.

[0021] S9. Solve the grid synchronization estimation model for ionospheric delay and scaling factor; perform block optimization on the weighting factor and scaling factor using the characteristics of the independent variables.

[0022] Preferably, the double-layer ionospheric thin-layer model includes an ionospheric thin layer to be estimated and a hidden top thin layer. The ionospheric thin layer to be estimated is located below the orbital altitude of the low-Earth orbit satellite and is used for calculating the puncture point location and tilt / vertical delay conversion of two types of observations: ground station and low-Earth orbit satellite, and ground station and GNSS. The hidden top thin layer is located above the orbital altitude of the low-Earth orbit satellite and is used for calculating the puncture point location and tilt / vertical delay conversion of low-Earth orbit satellite and GNSS observations.

[0023] Preferably, a rectangular-to-circular search scheme is designed during the ionospheric puncture point search process. The specific process is as follows: first, a rectangular area is calculated based on the grid point position and the maximum search radius; then, the distance between the puncture point and the grid point within the rectangular area is calculated; and finally, a circular search is performed based on this.

[0024] Preferably, the specific process of constructing the ionospheric delay and scaling factor grid synchronization estimation model in S8 is as follows: Assume that the final selected ionospheric puncture points at the specified ionospheric grid points contain n1 v T n2 v B And n3 v E , where v T v B v E , representing the vertical ionospheric delays after tilt delay conversion obtained from LEO satellite and GNSS, ground station and LEO satellite, and ground station and GNSS observations, respectively, with n1, n2, and n3 representing the corresponding number; the column vectors composed of vertical ionospheric delay observations from different sources are v1, v2, and v3, respectively. The estimated value of the ionospheric treatment delay at a specified ionospheric grid point is a linear combination of all selected ionospheric puncture point observations, i.e.:

[0025]

[0026] In the formula, K i This represents the scaling factor matrix corresponding to different types of observations, with the weighting vector w selected. i The criterion is to make variance To minimize this, we can construct the following optimization problem:

[0027]

[0028] In the formula, k0 represents the scaling factor between the vertical delay of the top ionosphere at a specified IGP and the vertical delay of the overall ionosphere, and s = [1 0 0]. T constraint matrix G i The elements are composed of:

[0029]

[0030] In the formula, Δx i,j To correspond to the difference in three-dimensional coordinates between the IPP and the selected IGP in the geocentric-geofixed coordinate system, l e and l n Here are the unit vectors for the east and north directions in the northeast-north coordinate system with the selected IGP as the origin:

[0031]

[0032] The objective function in equation (2) The expression is:

[0033]

[0034] In the formula, scalar c0 and vector c i Sum matrix N i,j Indicates the spatial correlation of the ionosphere;

[0035] The expression for c0 is:

[0036]

[0037] c i The j-th element is:

[0038]

[0039] N i,j The element in row p and column q is:

[0040]

[0041] In the formula This represents the covariance of the ionospheric residuals of two spatially overlapping IPPs. The covariance of the ionospheric residuals of two IPPs at infinite distances, d decorr Indicates the distance to the ionospheric space; and d decorrBased on long-term ionospheric statistical data, the service area of ​​SBAS needs to be given in advance;

[0042] In equation (5), matrix M i,j This represents the correlation between navigation signal measurement noise and hardware correction delay residuals; for the case i = j and p = q, the element in the p-th row and q-th column is:

[0043]

[0044] For other cases, matrix M i,j The element in row p and column q is:

[0045]

[0046] In the formula, matrix K i (i = 1, 2, 3) is determined by the scaling factors k0 and k e k n The position of IPP relative to IGP determines the amplification of ionospheric observations originating from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS. For i=1, the scaling factor matrix K1 for ionospheric observations originating from LEO satellites and GNSS is:

[0047]

[0048] In the formula, diag(·) is a diagonal matrix with the corresponding vector as its diagonal elements, and all other elements in the matrix are 0; for i=2, the ionospheric observations derived from ground stations and low-Earth orbit satellites are represented by the corresponding scaling factor matrix K2:

[0049]

[0050] For i=3, i.e., ionospheric observations originating from ground stations and GNSS, the corresponding scaling factor matrix K3 is:

[0051]

[0052] Preferably, the specific process of S9 solving the ionospheric delay and scaling factor grid synchronization estimation model is as follows:

[0053] S901, Initialization;

[0054] The parameters to be estimated are initialized by selecting the reciprocal of the number of IPPs involved in the estimation as the initial value of the weight vector, and the initial value of the scaling factor as the ratio of the average top vertical delay of the parameters involved in the estimation to the average overall vertical delay.

[0055]

[0056] In the formula: Indicates a length of n i A column vector with all elements being 1, z = [k0, k e ,k n ] T This represents a column vector consisting of scaling factors, with z0 being its initial value.

[0057] S902, Fixed z t Solve sub-optimization problem 1

[0058] Given the iteration value w from the previous step 1,t w 2,t w 3,t and z t Then, first fix z. t The optimization problem degenerates into:

[0059]

[0060] In the formula, c = [c1, c2, c3] T , W -1 The expression is:

[0061]

[0062] Matrix G T The structure is as follows:

[0063]

[0064] The sub-optimization problem 1 corresponding to equation (15) is a typical Kriging interpolation problem, and its analytical solution is:

[0065] y t+1 =[W-WG(G T WG) -1 G T W]c+WG(G T WG) -1 s (18)

[0066] S903, fixed y t+1 Solve sub-optimization problem 2

[0067] Get y t+1 Then, it is fixed in equation (2). Since the constraint in equation (2) has no effect on z, the constrained optimization problem degenerates into the following unconstrained optimization problem:

[0068]

[0069] In the formula, T is a constant independent of z;

[0070] To make the objective function value To reduce z, the gradient descent method is used. t Update:

[0071] z t+1 =z t -ηL T YMY T k (20)

[0072] In the formula, Y = diag(y t+1 ), where η is a positive real number, the elements of matrix L are:

[0073]

[0074] In the formula:

[0075]

[0076]

[0077] S904. Convergence determination;

[0078] The convergence metric is constructed as follows:

[0079] Δ=||y t+1 -y t || 2 +||z t+1 -z t || 2 (twenty four)

[0080] For a given convergence tolerance threshold eps, if Δ > eps, the optimization iteration process is determined to have not converged, and the process jumps to S903 to start a new round of iteration; otherwise, the iterative optimization process is determined to have converged, and the vector z is used to determine the convergence. t+1 k0 is taken as the optimal estimate of the scaling factor at the selected IGP, and expressed as a vector y t+1 w in 1,t+1 w 2,t+1 w 3,t+1 , with vector z t+1 Construct matrix K i Estimate the vertical ionospheric delay at the selected IGP

[0081]

[0082] In the formula, matrix K i (i = 1, 2, 3) is determined by the scaling factors k0 and k e k nAnd the position of IPP relative to IGP determines the amplification of ionospheric observations from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS. i This represents a column vector consisting of different vertical ionospheric delay observations.

[0083] Therefore, the present invention employs the aforementioned low-orbit enhanced SBAS ionospheric monitoring system and grid estimation method, which has the following beneficial effects:

[0084] (1) This invention provides a low-orbit enhanced SBAS ionospheric monitoring system based on low-orbit satellites and ground monitoring stations for single-frequency SBAS air navigation services, which helps to improve the accuracy and integrity of air navigation in the edge service area of ​​SBAS;

[0085] (2) This invention provides a method for synchronous estimation of ionospheric delay and scaling factor of low-orbit enhanced SBAS grid, which can make full use of the observation information of low-orbit satellites and ground monitoring stations, and helps to solve the ionospheric modeling problem caused by incomplete ionospheric observations of low-orbit satellites.

[0086] (3) This invention simplifies complex optimization models by using the block optimization approach of independent variables, and can improve the accuracy of grid models while ensuring the real-time generation of SBAS ionospheric grids.

[0087] (4) The present invention can simultaneously estimate the overall ionospheric vertical delay grid and the ratio factor grid of the vertical delay between the top ionospheric layer and the overall ionospheric layer, and can provide a high-precision ionospheric correction model for the current GNSS constellation and the future low-Earth orbit satellite navigation constellation.

[0088] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0089] Figure 1 Diagram of the Low Orbit Enhanced SBAS Ionospheric Monitoring System;

[0090] Figure 2 A schematic diagram of navigation, communication, and ionospheric observation for the low-Earth orbit enhanced SBAS ionospheric monitoring network;

[0091] Figure 3 Flowchart of a method for synchronizing ionospheric delay and scaling factor estimation in a low-orbit enhanced SBAS grid;

[0092] Figure 4 This diagram illustrates the method for searching ionospheric puncture points. Detailed Implementation

[0093] Example

[0094] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0095] Please see Figure 1 A low-Earth orbit enhanced SBAS ionospheric monitoring system includes a low-Earth orbit satellite-based monitoring platform, a ground monitoring station, and an information processing center;

[0096] The low-Earth orbit satellite-based monitoring platform relies on low-Earth orbit satellites and is equipped with an onboard GNSS receiver and a downlink navigation signal generation module. It has the ability to receive satellite-based GNSS navigation signals and broadcast independent dual-frequency navigation signals to the ground.

[0097] The ground monitoring station is equipped with multiple independent navigation receivers to receive and analyze GNSS observation information carried in GNSS navigation signals, low-orbit satellite navigation signals, and low-orbit satellite communication signals. It also has Internet access capabilities and aggregates various observation information to the information processing center for comprehensive analysis and processing.

[0098] The information processing center deploys high-performance computing equipment to receive and preprocess observation data from ground monitoring stations and low-Earth orbit satellites. It also runs a method for synchronously estimating ionospheric delay and scaling factors using a gridded SBAS (Suspension-Based Satellite Automation System) grid, generating ionospheric delay corrections and corresponding integrity parameters for SBAS. This involves navigation signals and communication signals such as... Figure 2 As shown. The Low Earth Orbit (LEO) satellite carries a GNSS receiver capable of receiving GNSS navigation signals and is equipped with a navigation signal generation module for broadcasting dual-frequency navigation signals to the ground. The ground station (MS) can simultaneously receive both GNSS and LEO dual-frequency navigation signals. A communication link can be established between the LEO and MS to transmit the GNSS navigation information received by the LEO. The communication link between the MS and the ground information processing center is used to send MS and LEO observation information to the information processing center. Ionospheric observations are included in the dual-frequency navigation ranging signals. Using the LEO satellite orbit as a boundary, LEO-GNSS observations contain upper ionospheric information, MS-LEO observations contain lower ionospheric information, and MS-GNSS observations contain overall ionospheric information. The information processing center comprehensively utilizes LEO-GNSS, MS-LEO, and MS-GNSS observation information to establish an overall ionospheric grid model, and simultaneously estimates the scaling factor grid between the overall ionospheric delay and the upper ionospheric delay.

[0099] Please see Figure 3A method for estimating the ionospheric grid of LEO-enhanced SBAS (Balanced Low Earth Orbit Satellite System) is proposed. This method utilizes three types of navigation observation data—LEO satellite and GNSS, ground station and LEO satellite, and ground station and GNSS—to estimate the ionospheric delay and scaling factor of the LEO-enhanced SBAS grid. The specific process for estimating the ionospheric delay and scaling factor of the LEO-enhanced SBAS grid is as follows:

[0100] S1. Constructing geometry-free observables; constructing corresponding geometry-free observables using dual-frequency pseudorange and carrier phase observations from navigation data; the specific process is as follows:

[0101] The raw data arriving at the information processing center includes ephemeris data from GNSS and LEO satellites, as well as dual-frequency pseudorange and carrier phase observations received by MS and LEO navigation receivers. Without distinguishing between observation types, the pseudorange and carrier phase observations are denoted as follows: and Where k represents the epoch and q represents the frequency. Corresponding geometry-free (GF) observations are generated from pseudorange and carrier phase observations respectively:

[0102]

[0103] In the formula L GF Indicates that there are no geometric observations of the carrier phase, λ q This represents the wavelength of the carrier wave at the corresponding frequency point (q = 1, 2).

[0104] S2. Cycle slip detection: After acquiring geometrically unaffected observations, cycle slip detection is performed on the carrier phase observations to remove those affected by cycle slips; the specific process is as follows:

[0105] After obtaining the geometrically unaffected observations, cycle slip detection is required on the carrier phase observations to prevent abnormal carrier phase from affecting the subsequent filtering and smoothing process. The test statistic t for cycle slip detection... cs The structure is the difference between geometric carrier phase observations from one epoch to the next:

[0106] t cs =|L GF,k -L GF,k-1 | (27)

[0107] The corresponding detection threshold T cs for:

[0108]

[0109] In the formula, θ represents the elevation angle of the transmitter relative to the receiver, and b GF Related to the receiver sampling period T, its value is:

[0110]

[0111] If in a smooth period N CCL All geometrically carrier-free observables satisfy t cs If ≤T, it means that no cycle slip was detected, and the corresponding navigation observations will participate in the subsequent filtering and smoothing process; otherwise, the corresponding observations should be excluded from the filtering, smoothing and ionosphere extraction processes.

[0112] S3. Filtering and smoothing: Filtering and smoothing are performed on the geometrically uncarrier-free observations obtained through cycle slip detection and their corresponding geometrically unpseudorange observations to improve the accuracy of tilted ionospheric delay extraction. The specific process is as follows:

[0113] The geometrically uncarrier-free observations obtained through cycle slip detection, along with their corresponding geometrically uncarrier-free pseudorange observations, are filtered and smoothed to reduce observation noise and improve the accuracy of tilted ionospheric delay extraction. The smoothing operation is as follows:

[0114]

[0115] In the formula L CCL,k N represents the smoothed observation value at the current epoch. CCL This indicates the smoothing window length. For the three different observations—MS-LEO, MS-GNSS, and LEO-GNSS—different smoothing window lengths can be selected based on their motion characteristics to achieve a trade-off between IPP density and smoothing performance. Smoothed observation L CCL,k noise variance for:

[0116]

[0117] In the formula, σ P and σ L σ represents the standard deviation of pseudorange and carrier measurement errors, respectively. m and σ n These represent multipath error and receiver measurement thermal noise, respectively.

[0118] S4. Obtaining the tilted ionospheric delay: Using differential code bias products to correct hardware bias in the smoothing results, further obtain the observed values ​​of the tilted ionospheric delay and the corresponding measurement error variance; the specific process is as follows:

[0119] The smoothed observations include two components: ionospheric delay, receiver hardware delay, transmitter hardware delay, and small observation noise. Hardware delay is represented by differential code bias (DCB), which is corrected using a DCB product for the corresponding frequency to obtain the tilted ionospheric delay observations.

[0120]

[0121] in, f1 and f2 are the frequencies of the corresponding navigation signal points, respectively; B r and B t These are the receiver and transmitter DCBs provided by the product, respectively. Since DCB products also have uncertainties, the variance of the tilted ionospheric delay observation error... for:

[0122]

[0123] In the formula, and The variances of the DCB correction residuals at the receiver and transmitter are respectively.

[0124] S5. Ionospheric Penetration Point Location Calculation: A two-layer ionospheric thin-layer model is constructed. Based on the coordinates of the ground monitoring station and GNSS / Low Earth Orbit (LEO) satellite ephemeris data, the penetration point locations of each tilted ionospheric observation on the two-layer ionospheric thin-layer model are calculated. The two-layer ionospheric thin-layer model includes the ionospheric thin layer to be estimated and a hidden top thin layer. The ionospheric thin layer to be estimated is located below the LEO satellite orbital altitude and is used for calculating the penetration point locations and tilt / vertical delay conversion between ground station and LEO satellite, and between ground station and GNSS observations. The hidden top thin layer is located above the LEO satellite orbit and is used for calculating the penetration point locations and tilt / vertical delay conversion between LEO satellite and GNSS observations. The specific process is as follows:

[0125] The point where the line of sight formed by the receiver and transmitter of the navigation signal intersects the thin shell of the ionosphere is called the ionospheric penetration point (IPP). Therefore, the coordinates of the IPP are determined by the positions of the receiver and transmitter in the Earth-centered Earth-fixed coordinate system. The Earth's central angle ψ corresponding to the ionospheric penetration point is... IPP for:

[0126]

[0127] In the formula R e h represents the Earth's radius. s This represents the height of the ionospheric shell. For MS-LEO and MS-GNSS observations, the puncture point is located on the ionospheric shell to be modeled, and its height is h. s =350km; For LEO-GNSS observations, the puncture point is located in the hidden top ionospheric shell, and its height is determined by the following formula:

[0128] h s =(0.0027F) 107 +1.79)h o -5.52F 107 +1350 (35)

[0129] In the formula, h o F represents the orbital altitude of a low-Earth orbit satellite. 107This indicates the intensity index of solar activity.

[0130] IPP geographical latitude φ IPP The calculation formula is:

[0131] φ IPP =sin -1 (sinφ r cosψ pp +cosφ r sinψ pp cosα) (36)

[0132] In the formula, φ r α is the geographical latitude of the receiver, and α is the azimuth angle of the transmitter relative to the receiver.

[0133] IPP Geographic Longitude λ IPP The calculation needs to be divided into two cases:

[0134] When φ r >70° and tanψ IPP cosθ>tan(π / 2-φ r ) or φ r When < -70°:

[0135]

[0136] otherwise:

[0137]

[0138] In the formula λ r Indicates the geographical longitude of the recipient.

[0139] S6. Tilt / Vertical Ionospheric Delay Conversion: Based on the assumption of ionospheric isotropy, the tilted ionospheric delay from different sources is converted into the vertical ionospheric delay using a projection function; the specific process is as follows:

[0140] The tilted ionospheric delay extracted from LEO-GNSS, MS-LEO, and MS-GNSS observations is denoted as s. T s B and s E The ionospheric delay grid models the vertical ionospheric delay, and converting the tilted ionospheric delay to the vertical ionospheric delay is a necessary step. Based on the assumption of ionospheric isotropy, the conversion formula for ionospheric delay is:

[0141]

[0142] The tilt delay s obtained from MS-LEO and MS-GNSS observations B and s EThe transformation process is implemented on the thin shell of the ionosphere to be modeled, and the projection function M(θ) is of the form:

[0143]

[0144] For the tilt delay s obtained from LEO-GNSS observations T The transformation process occurs within the implicit top ionospheric thin shell, and the projection function M(θ) takes the form:

[0145]

[0146] s T s B and s E The corresponding vertical ionospheric delays after conversion are v T v B and v E .

[0147] S7. Ionospheric puncture point search; Searching for ionospheric puncture points within the vicinity of the grid point to be estimated in the ionosphere; The specific process is as follows:

[0148] For a given IGP (Integrated Gateway Point), the spatial correlation of ionospheric delay weakens with increasing distance from other IPPs (Integrated Power Points). When estimating the ionospheric delay at an IGP, only nearby IPPs need to be considered. The process of searching for nearby IPPs is as follows: Figure 4 As shown, it is necessary to calculate the maximum search radius R. max Minimum search radius R min Minimum number of IPPs N min And the target number of IPPs N obj The setup is as follows. Specifically, in the ionospheric puncture point search process, a rectangular-to-circular search scheme was designed. The process involves first calculating a rectangular area based on the grid point positions and the maximum search radius, then calculating the distance between the puncture point and the grid point within the rectangular area, and finally performing a circular search based on this. The specific process is as follows:

[0149] Given R max Next, the rectangular region defined by δφ and δλ is delineated. δφ and δλ are calculated as follows:

[0150]

[0151] In the formula, h s φ represents the height of the thin ionospheric shell to be modeled. IGP This indicates the selected geographical latitude of the IGP.

[0152] For all IPPs within the rectangular region, calculate the distance between each IPP and the IGP on the modeled ionospheric shell, and sort them in ascending order. If R min The number of IPPs within is greater than Nobj Then the fitting radius R fit Set to R min Otherwise, gradually increase R. fit Until R fit The number of IPPs equals N obj If R fit Expand to R max When the number of IPPs is greater than N min Then let R fit =R max Otherwise, mark the IGP as unavailable.

[0153] If an IGP is marked as unavailable, proceed to the next IGP for IPP selection; otherwise, execute S8, utilizing R... fit Estimate the vertical ionospheric delay at the IGP within the IPP of a circle with radius .

[0154] S8. Construct a grid synchronization estimation model for ionospheric delay and scaling factor; comprehensively considering the spatial correlation of the ionosphere, measurement error characteristics, and the relative positional relationship between grid points and puncture points, construct an optimization model with minimizing the variance of the vertical delay estimation error of grid points as the objective function and weighting factor and scaling factor as independent variables; the specific process is as follows:

[0155] Assume that the final selected IPP at the specified IGP contains n1 v T n2 v B And n3 v E , where v T v B v E Let v1, v2, and v3 represent the vertical ionospheric delays after tilt delay conversion obtained from LEO satellite and GNSS, ground station and LEO satellite, and ground station and GNSS observations, respectively. The column vectors composed of vertical ionospheric delay observations from different sources are v1, v2, and v3, respectively. The estimated ionospheric delay at the specified IGP is a linear combination of all selected IPP observations, i.e.:

[0156]

[0157] In the formula, K i This represents the scaling factor matrix corresponding to different types of observations, with the weighting vector w selected. i The criterion is to make variance To minimize this, we can construct the following optimization problem:

[0158]

[0159] In the formula, k0 represents the scaling factor between the vertical delay of the top ionosphere at a specified IGP and the vertical delay of the overall ionosphere, and s = [1 0 0]. T constraint matrix G i The elements are composed of:

[0160]

[0161] In the formula, Δx i,j To correspond to the difference in three-dimensional coordinates between the IPP and the selected IGP in the geocentric-geofixed coordinate system, l e and l n Here are the unit vectors for the east and north directions in the northeast-north coordinate system with the selected IGP as the origin:

[0162]

[0163] The objective function in equation (2) The expression is:

[0164]

[0165] In the formula, scalar c0 and vector c i Sum matrix N i,j Indicates the spatial correlation of the ionosphere;

[0166] The expression for c0 is:

[0167]

[0168] c i The j-th element is:

[0169]

[0170] N i,j The element in row p and column q is:

[0171]

[0172] In the formula This represents the covariance of the ionospheric residuals of two spatially overlapping IPPs. The covariance of the ionospheric residuals of two IPPs at infinite distances, d decorr Indicates the distance to the ionospheric space; and d decorr Based on long-term ionospheric statistical data, the service area of ​​SBAS needs to be given in advance;

[0173] In equation (5), matrix M i,jThis represents the correlation between navigation signal measurement noise and hardware correction delay residuals; for the case i = j and p = q, the element in the p-th row and q-th column is:

[0174]

[0175] For other cases, matrix M i,j The element in row p and column q is:

[0176]

[0177] In the formula, matrix K i (i = 1, 2, 3) is determined by the scaling factors k0 and k e k n The position of IPP relative to IGP determines the amplification of ionospheric observations originating from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS. For i=1, the scaling factor matrix K1 for ionospheric observations originating from LEO satellites and GNSS is:

[0178]

[0179] In the formula, diag(·) is a diagonal matrix with the corresponding vector as its diagonal elements, and all other elements in the matrix are 0; for i=2, the ionospheric observations derived from ground stations and low-Earth orbit satellites are represented by the corresponding scaling factor matrix K2:

[0180]

[0181] For i=3, i.e., ionospheric observations originating from ground stations and GNSS, the corresponding scaling factor matrix K3 is:

[0182]

[0183] S9. Solve the grid synchronization estimation model for ionospheric delay and scaling factor; utilize the characteristics of the independent variables to perform block optimization on the weighting factor and scaling factor. The specific process is as follows:

[0184] S901, Initialization;

[0185] The parameters to be estimated are initialized by selecting the reciprocal of the number of IPPs involved in the estimation as the initial value of the weight vector, and the initial value of the scaling factor as the ratio of the average top vertical delay of the parameters involved in the estimation to the average overall vertical delay.

[0186]

[0187] In the formula: Indicates a length of n i A column vector with all elements being 1, z = [k0, ke [kn] T This represents a column vector consisting of scaling factors, with z0 being its initial value.

[0188] S902, Fixed z t Solve sub-optimization problem 1

[0189] Given the iteration value w from the previous step 1,t w 2,t w 3,t and z t Then, first fix z. t The optimization problem degenerates into:

[0190]

[0191] In the formula, c = [c1, c2, c3] T , W -1 The expression is:

[0192]

[0193] Matrix G T The structure is as follows:

[0194]

[0195] The sub-optimization problem 1 corresponding to equation (15) is a typical Kriging interpolation problem, and its analytical solution is:

[0196] y t+1 =[W-WG(G T WG) -1 G T W]c+WG(G T WG) -1 s (18)

[0197] S903, fixed y t+1 Solve sub-optimization problem 2

[0198] Get y t+1 Then, it is fixed in equation (2). Since the constraint in equation (2) has no effect on z, the constrained optimization problem degenerates into the following unconstrained optimization problem:

[0199]

[0200] In the formula, T is a constant independent of z;

[0201] To make the objective function value To reduce z, the gradient descent method is used. t Update:

[0202] z t+1 =z t -ηL T YMY T k (20)

[0203] In the formula, Y = diag(y t+1 ), where η is a positive real number, the elements of matrix L are:

[0204]

[0205] In the formula:

[0206]

[0207]

[0208] S904. Convergence determination;

[0209] The convergence metric is constructed as follows:

[0210] Δ=||y t+1 -y t || 2 +||z t+1 -z t || 2 (twenty four)

[0211] For a given convergence tolerance threshold eps, if Δ > eps, the optimization iteration process is determined to have not converged, and the process jumps to S903 to start a new round of iteration; otherwise, the iterative optimization process is determined to have converged, and the vector z is used to determine the convergence. t+1 k0 is taken as the optimal estimate of the scaling factor at the selected IGP, and expressed as a vector y t+1 w in 1,t+1 w 2,t+1 w 3,t+1 , with vector z t+1 Construct matrix K i Estimate the vertical ionospheric delay at the selected IGP

[0212]

[0213] In the formula, matrix K i (i = 1, 2, 3) is determined by the scaling factors k0 and k e k n And the position of IPP relative to IGP determines the amplification of ionospheric observations from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS. i This represents a column vector consisting of different vertical ionospheric delay observations.

[0214] Therefore, this invention employs the aforementioned low-Earth orbit enhanced SBAS ionospheric monitoring system and grid estimation method. Addressing the shortcomings of insufficient coverage and poor monitoring performance in service edge areas of the SBAS ground monitoring network, it utilizes the onboard receiver and navigation enhancement module of the low-Earth orbit satellite to provide richer ionospheric observation information. This improves the estimation performance of ionospheric grid points in service edge areas while ensuring the real-time performance of the SBAS ionospheric grid. Simultaneously, it uses multi-level ionospheric observation information to obtain the scaling factor of the vertical delay between the top ionosphere and the overall ionosphere, meeting the future requirements for high-precision ionospheric correction in low-Earth orbit satellite navigation.

[0215] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for estimating the ionospheric grid of a low-orbit enhanced SBAS, characterized in that, The ionospheric delay and scaling factor of the LEO-enhanced SBAS grid are estimated using three types of navigation observation data: LEO satellite and GNSS, ground station and LEO satellite, and ground station and GNSS. The specific process for estimating the ionospheric delay and scaling factor of the LEO-enhanced SBAS grid is as follows: S1. Construct geometry-free observables; construct the corresponding geometry-free observables using dual-frequency pseudorange and carrier phase observations from navigation data; S2. Cycle slip detection: After acquiring geometrically unaffected observations, cycle slip detection is performed on the carrier phase observations to eliminate those affected by cycle slips. S3, Filtering and smoothing processing; Filtering and smoothing are performed on geometrically uncarrier-free observations obtained through cycle slip detection and their corresponding geometrically unpseudorange observations to improve the accuracy of tilted ionospheric delay extraction. S4. Obtaining the tilted ionospheric delay: Using differential code bias products to correct hardware bias in the smoothing results, the observed values ​​of the tilted ionospheric delay and the corresponding measurement error variance are obtained. S5. Calculation of ionospheric puncture point location; Construct a thin double-layer ionospheric model, and calculate the puncture point location of each tilted ionospheric observation on the thin double-layer ionospheric model based on the coordinates of the ground monitoring station and GNSS / Low Orbit satellite ephemeris data; S6. Inclined / Vertical Ionospheric Delay Conversion: Based on the assumption of ionospheric isotropy, the tilted ionospheric delay from different sources is converted into vertical ionospheric delay using a projection function. S7. Ionospheric puncture point search; Search for ionospheric puncture points within the vicinity of the grid point to be estimated in the ionospheric region; S8. Construct a grid synchronization estimation model for ionospheric delay and scaling factor; comprehensively consider the spatial correlation of the ionosphere, measurement error characteristics, and the relative positional relationship between grid points and puncture points, construct an optimization model with minimizing the variance of the vertical delay estimation error of grid points as the objective function and weighting factor and scaling factor as independent variables. S9. Solve the grid synchronization estimation model for ionospheric delay and scaling factor; perform block optimization on the weighting factor and scaling factor using the characteristics of the independent variables.

2. The method for estimating the ionospheric grid of a low-orbit enhanced SBAS according to claim 1, characterized in that: The double-layer ionospheric thin-layer model includes an ionospheric thin layer to be estimated and a hidden top thin layer. The ionospheric thin layer to be estimated is located below the orbital altitude of the low-Earth orbit satellite and is used for calculating the puncture point location and tilt / vertical delay conversion of two types of observations: ground station and low-Earth orbit satellite, and ground station and GNSS. The hidden top thin layer is located above the orbital altitude of the low-Earth orbit satellite and is used for calculating the puncture point location and tilt / vertical delay conversion of low-Earth orbit satellite and GNSS observations.

3. The method for estimating the ionospheric grid of a low-orbit enhanced SBAS according to claim 1, characterized in that: In the process of searching for ionospheric puncture points, a rectangular to circular search scheme was designed. The specific process is as follows: first, a rectangular area is calculated based on the grid point positions and the maximum search radius; then, the distance between the puncture point and the grid point within the rectangular area is calculated; and finally, a circular search is performed based on this.

4. The method for estimating the ionospheric grid of a low-orbit enhanced SBAS according to claim 1, characterized in that, The specific process of constructing the ionospheric delay and scaling factor grid synchronization estimation model in S8 is as follows: Assume that the final selected ionospheric puncture point at the specified ionospheric grid point includes... indivual , indivual as well as indivual ,in, , , These represent the vertical ionospheric delays after tilt delay conversion obtained from LEO satellite and GNSS observations, ground station and LEO satellite observations, and ground station and GNSS observations, respectively. , , These represent the corresponding numbers; the column vectors composed of vertical ionospheric delayed observations from different sources are respectively , and The estimated ionospheric treatment delay at a specified ionospheric grid point is a linear combination of observations from all selected ionospheric puncture points, i.e.: （1） In the formula, This represents the scaling factor matrix corresponding to different types of observations, and the weighting vector is selected. The criterion is to make variance To minimize this, we can construct the following optimization problem: （2） In the formula, This represents the scaling factor between the vertical delay of the top ionosphere at the specified IGP and the overall vertical delay of the ionosphere. constraint matrix The elements are composed of: （3） In the formula, This is the difference between the three-dimensional coordinates of the corresponding IPP and the selected IGP in the geocentric-geofixed coordinate system. and Here are the unit vectors for the east and north directions in the northeast-north coordinate system with the selected IGP as the origin: （4） The objective function in equation (2) The expression is: （5） In the formula, scalar ,vector sum matrix Indicates the spatial correlation of the ionosphere; The expression is: （6） The Middle The elements are: （7） The Middle Okay, number The elements of the column are: （8） In the formula This represents the covariance of the ionospheric residuals of two spatially overlapping IPPs. The covariance of the ionospheric residuals of two IPPs at infinite distances. Indicates the distance to the ionospheric space; , and Based on long-term ionospheric statistical data, the service area of ​​SBAS needs to be given in advance; In equation (5) the matrix This indicates the correlation between navigation signal measurement noise and hardware correction delay residuals; for and In the case of its first Okay, number The elements of the column are: （9） For other cases, the matrix No. Okay, number The elements of the column are: （10） In the formula matrix ( ) by scaling factor , , And the position of IPP relative to IGP determines the amplification of ionospheric observations from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS; for The scale factor matrix originates from ionospheric observations from low-Earth orbit satellites and GNSS. for: （11） In the formula, A diagonal matrix with the corresponding vector as its diagonal element, where all other elements are 0; for The scale factor matrix originates from ionospheric observations from ground stations and low-Earth orbit satellites. for: （12） for That is, the scale factor matrix corresponding to ionospheric observations from ground stations and GNSS. for: （13）。 5. The method for estimating the ionospheric grid of a low-orbit enhanced SBAS according to claim 1, characterized in that, The specific process of solving the ionospheric delay and scaling factor grid synchronization estimation model in S9 is as follows: S901, Initialization; The parameters to be estimated are initialized by selecting the reciprocal of the number of IPPs involved in the estimation as the initial value of the weight vector, and the initial value of the scaling factor as the ratio of the average top vertical delay of the parameters involved in the estimation to the average overall vertical delay. （14） In the formula: Indicates length is A column vector whose elements are all 1s. This represents a column vector composed of scaling factors. It is its initial value; S902, Fixed Solve sub-optimization problem 1; Given the iteration value from the previous step , , and Then, first fix The optimization problem (2) degenerates into: （15） In the formula, , , The expression is: （16） matrix The structure is as follows: （17） The sub-optimization problem 1 corresponding to equation (15) is a typical Kriging interpolation problem, and its analytical solution is: （18） S903, Fixed Solve sub-optimization problem 2; get Then, it is fixed in equation (2), because the constraint of equation (2) has an effect on... This is ineffective; the constrained optimization problem degenerates into the following unconstrained optimization problem: （19） In the formula To and Irrelevant constants; To make the objective function value To reduce, the gradient descent method is used. Update: （20） In the formula, , , For positive real numbers, the matrix The elements are composed of: （21） In the formula: （22） （23） S904. Convergence determination; The convergence metric is constructed as follows: （24） For a given convergence tolerance threshold value ,if If the optimization iteration process fails, it is determined that the process has not converged, and the process jumps to S903 to start a new round of iteration. Otherwise, it is determined that the iterative optimization process has converged, using vectors. In As the optimal estimate of the scaling factor at the selected IGP, and expressed as a vector In , , , with vector Constructing a matrix Estimate the vertical ionospheric delay at the selected IGP : （25） In the formula matrix ( ) by scaling factor , , And the position of the IPP relative to the IGP determines the amplification of ionospheric observations from LEO satellites and GNSS, ground stations and LEO satellites, and ground stations and GNSS. This represents a column vector consisting of different vertical ionospheric delay observations.

Images