A low ionosphere electron density bottom-up inversion method

Abstract

The present application relates to a kind of low ionosphere electron density from bottom to top inversion method, comprising: in the first height range, to super phase is weighted linear fitting, obtain background phase, and extrapolate the background phase obtained by fitting to higher height;According to background phase and super phase, calculate the residual phase, according to the residual phase, calculate the total electron content of height;The total electron content of height of multiple height is arranged according to the tangent point height to form the total electron content of height matrix, obtain the total electron content profile of height;Abel weighted forward model is constructed;Based on the total electron content profile of height, optimal estimation inversion is carried out, and inversion electron density is obtained;And background and layered enhancement separation is carried out to inversion electron density.The method of the present application avoids the resolution reduction and error amplification problem caused by F area long tail in traditional Abel inversion, significantly improves the inversion accuracy, and realizes the effective separation of background and layered structure.

Description

Technical Field

[0001] This invention relates to the field of ionospheric detection and space weather monitoring technology, and in particular to a bottom-up inversion method for low ionospheric electron density based on GNSS radio occultation (RO) observation data. Background Technology

[0002] The lower ionosphere (approximately 60–130 km, including regions D (60–90 km) and E (90–130 km)) is a critical region controlling the availability of high-frequency communications, ionospheric delay and phase stability of GNSS (National Navigation Satellite System), and the timeliness of space weather monitoring. Structural changes in this layer often lead to link interruptions, positioning errors, and artifacts during disturbance periods. Due to the lack of global coverage for ground-based systems and the instability of classical inversion methods at low altitudes due to upper-layer influences, reliable low-altitude electron density profiles remain of significant value.

[0003] Existing ionospheric electron density inversion techniques primarily employ the classic Abel inversion method. This method maps tangent height information to a radial profile; however, in the 60-130 km range, geometrically weighted broadening and coupling with the F region (high ionosphere, above 150 km) lead to reduced resolution and amplified errors at lower altitudes. Specifically, the following problems exist: (1) Classical Abel inversion is highly sensitive to the F-region ionosphere and its variability. Due to the long-tail effect of the Abel function in the F-region, the E-region inversion must deal with the significant tail of the F-region weighting function, making the E-region electron density inversion very sensitive to the F-region residue and variability.

[0004] (2) Existing techniques typically rely on auxiliary data or prior models to infer the electron density in the E-region. Ionospheric inhomogeneities and auxiliary data can lead to significant biases in the inversion of the electron density in the E-region, depending on the prior information used and how it is constructed.

[0005] (3) Traditional methods require the simultaneous use of L1 and L2 dual-frequency data, but L2 data has high noise, which limits the number of available profiles. In addition, existing technologies have difficulty detecting small-scale electron density fluctuations in the E region (such as E region abrupt changes) because these structures have small amplitudes and require higher measurement accuracy, while the noise level of L2 data limits the ability to identify such fine structures.

[0006] (4) Existing quality control methods are not capable enough to handle strong local signal-to-noise ratio (SNR) variability and abnormal behavior above 90km, resulting in unreliable profiles entering the inversion process and affecting the accuracy of the final results. Summary of the Invention

[0007] To address at least some of the problems mentioned above in the prior art, the present invention provides a bottom-up inversion method for low ionospheric electron density, comprising: Within the first altitude range, a weighted linear fit is performed on the superphase to obtain the background phase, and the fitted background phase is extrapolated to a higher altitude. The remaining phase is calculated based on the background phase and superphase, and the total electron content at elevation is calculated based on the remaining phase. The total electron content at multiple elevations is arranged according to the height of the tangent point to form the total electron content at elevation. The matrix is ​​used to obtain the total electron content profile at elevation. Construct an Abel-weighted forward model; Optimal estimation inversion was performed based on the total electron content profile at elevation to obtain the inverted electron density; and Background and layered enhancement separation are performed on the inverted electron density.

[0008] Furthermore, before constructing the Abel weighted forward model, the process also includes: quality control of the total elevation electron content profile, removing unqualified total elevation electron content profiles.

[0009] Furthermore, the first altitude range is 50-80km; Background phase fitting function ; With tangent height For the independent variable to be in superphase We performed a weighted linear fit to obtain the background phase corresponding to different tangent point heights: , W represents the weighting function, and the weights... It increases with altitude.

[0010] Furthermore, the remaining phase is calculated based on the background phase and superphase, and the total elevation electron content is calculated based on the remaining phase, including: The remaining phase is defined as: , Indicates background phase, Indicates superphase; The relationship between residual phase and total electron content at elevation: , Indicates the total electron content at elevation. Indicates the remaining phase. It is the GPS L1 frequency, 1575.42MHz.

[0011] Furthermore, constructing the Abel weighted forward model includes: When the tangent point heights are approximately the same, the elevation total electron content profile and the target profile are similar. The relationship follows a discrete Abel integral; Construct a discrete Abel-type integral forward model on a uniform height grid: with a step size of Uniform height grid Above, order and ,in This represents the distance from the current location to the Earth's center. Re represents the distance from the Earth's center to the position at time t=i, and Re represents the Earth's radius; Multiply each column of the discrete Abel integral forward model by This yields the Abel weighted forward model: , Represents the weighting function. express ; , in Scene matching background provided for the FIRI-2018 empirical model; The vector form of the Abel weighted forward model is: , in ε represents the measurement noise.

[0012] Furthermore, based on the total elevation electron content profile, optimal estimation inversion is performed to obtain the inverted electron density, including: Construct the observation covariance matrix: ,in , SNR represents the lower limit of signal-to-noise ratio; The short correlation index kernel is used as the prior covariance: ; The formula for calculating the optimal estimate of electron density under zero prior mean is: , Represents electron density, hTEC indicates , , Represent the covariance matrix; When calculating, the above formula should be written as: .

[0013] Furthermore, the inverted electron density is subjected to background and layered enhancement separation, including: The inverted electron density is averaged along the height grid to obtain the initial long-wave background Ne(0); Iterative annealing of the initial long-wave background yields the long-wave background field. ; The corresponding layered reinforcement component is , This is for inverting electron density.

[0014] Furthermore, it also includes: constructing a global context field, including: Background values ​​were sampled at 100 km in each total electron content profile obtained through quality control at elevation. Using natural neighborhood interpolation combined with mild Gaussian smoothing, the data is mapped onto a 0.5°×0.5° global grid to generate a global background electron density distribution map of the low ionosphere.

[0015] Furthermore, it also includes: preprocessing the occultation observation data, calculating the tangent point position and tangent point altitude, transforming the tangent point coordinates from the geocentric inertial coordinate system to the WGS-84 geodetic coordinate system, and calculating the solar zenith angle; Occultation observation data includes superphase Signal-to-noise ratio (SNR) and cut-off point height GNSS ephemeris, LEO ephemeris, and UTC time information.

[0016] The present invention has at least the following beneficial effects: This invention employs a bottom-up background estimation and removal method. It removes the F-region background at the measurement level, and the geometric kernel for E-region electron density inversion sharpens dramatically above the tangent point height. This avoids the resolution reduction and error amplification problems caused by the long tail of the F-region in traditional Abel inversion, significantly improving inversion accuracy. Compared with the IRI model, the background electron density at 100km is on the order of 1-8×10¹ 0 (e·m⁻³) is consistent with climatological expectations.

[0017] This invention relies solely on single-frequency (L1) 50Hz superphase data, without the need for an auxiliary F-region model or prior profile, to achieve a completely self-consistent E-region electron density inversion. This avoids the bias introduced by traditional methods due to improper prior selection, and improves the universality and portability of the method.

[0018] Innovatively, the background electron density of the FIRI-2018 empirical model is used as a column weight in the physically guided embedding Abel forward matrix, rather than being imposed on the solution as a prior mean; the FIRI weighting is significantly tightened. The spectral distribution reduces the condition number, allowing weaker regularization to meet the requirements for stable solutions. The adaptive diagonal loading mechanism automatically intervenes when conditions deteriorate, ensuring numerical stability under various ionospheric conditions.

[0019] This invention, based on SNR local variability detection and negative residual ratio screening above 90km, effectively eliminates unreliable profiles affected by multipath interference, strong gradients, or thin layers, thus avoiding the easy contamination of inversion results by unreliable profiles. In 932 test profiles, this method achieved a stable result with a pass rate of 40%, without background fitting failures or missing 100km.

[0020] This invention uses an iterative annealing method to adaptively identify and separate layered enhancements (such as occasional E layers) on a scale of 2-10km from long-wave background. It does not require a preset fixed template and can retain sharp layers and stabilize the background under different local conditions, thus achieving effective separation of background and layered structure.

[0021] This invention enables the reproducible production of global background maps from raw data through standardized batch processing, clear parameter and unit definitions, and complete log recording; the 0.5° grid global map supports systematic analysis of seasonal and local time variations.

[0022] The entire method of this invention can be implemented purely in software and deployed on a general computing platform. The vectorized kernel and condition checks avoid over-computation of pathological cases, making it suitable for batch processing and quasi-business operation, and it has good practicality. Attached Figure Description

[0023] To further illustrate the above and other advantages and features of the various embodiments of the present invention, a more specific description of the embodiments of the invention will be presented with reference to the accompanying drawings. It is to be understood that these drawings depict only typical embodiments of the invention and are therefore not intended to limit its scope. In the drawings, identical or corresponding parts will be indicated by identical or similar reference numerals for clarity.

[0024] Figure 1 The flowchart of a bottom-up inversion method for low ionospheric electron density according to the present invention is shown.

[0025] Figure 2 A geometrical schematic diagram of the tangent point of the present invention is shown.

[0026] Figure 3 A superphase and background fitting plot of an embodiment of the present invention is shown.

[0027] Figure 4 A superphase and background fitting plot of an embodiment of the present invention is shown.

[0028] Figure 5 A superphase and background fitting plot of an embodiment of the present invention is shown.

[0029] Figure 6 A superphase and background fitting plot of an embodiment of the present invention is shown.

[0030] Figure 7 A cross-sectional view of hTEC according to an embodiment of the present invention is shown.

[0031] Figure 8 A cross-sectional view of hTEC according to an embodiment of the present invention is shown.

[0032] Figure 9 A cross-sectional view of hTEC according to an embodiment of the present invention is shown.

[0033] Figure 10 A cross-sectional view of hTEC according to an embodiment of the present invention is shown.

[0034] Figure 11 The diagram shows a triptych of the inversion chain of Case a, an embodiment of the present invention.

[0035] Figure 12 The inversion chain triptych of Case b, an embodiment of the present invention, is shown.

[0036] Figure 13 The inversion chain triptych of Case c, an embodiment of the present invention, is shown.

[0037] Figure 14 A 100km global background electron density distribution map of an embodiment of the present invention is shown. Detailed Implementation

[0038] It should be noted that the components in the accompanying drawings may be shown exaggerated for illustrative purposes and may not be to scale.

[0039] In this invention, the various embodiments are merely intended to illustrate the solutions of the invention and should not be construed as limiting.

[0040] In this invention, unless otherwise specified, the quantifiers “a” and “one” do not exclude scenarios involving multiple elements.

[0041] It should also be noted that, in the embodiments of the present invention, only a portion of the parts or components may be shown for clarity and simplicity. However, those skilled in the art will understand that, under the teachings of the present invention, the required parts or components can be added as needed for specific scenarios.

[0042] It should also be noted that within the scope of this invention, the terms "same", "equal", and "equal to" do not mean that the two values ​​are absolutely equal, but allow for a certain reasonable error. In other words, the terms also cover "substantially the same", "substantially equal", and "substantially equal to".

[0043] It should also be noted that in the description of this invention, the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the invention and for simplifying the description, and do not explicitly or implicitly suggest that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0044] Furthermore, the embodiments of the present invention describe the process steps in a specific order. However, this is only for the convenience of distinguishing each step, and is not a limitation on the order of each step. In different embodiments of the present invention, the order of each step can be adjusted according to the process.

[0045] Existing ionospheric electron density inversion techniques have the following main drawbacks: (1) The existing Abel inversion method has a long tail interference problem in the F region in the low ionosphere (60-130km), which leads to low resolution and amplified error in the electron density inversion in the E region.

[0046] (2) Traditional methods are highly dependent on auxiliary data and prior models, have poor self-consistency, and are difficult to adapt to the robust inversion requirements under different ionospheric conditions.

[0047] (3) Existing technologies are unable to effectively separate the background electron density of the ionosphere from the layered enhancement structure (such as the sporadic E layer), resulting in inaccurate construction of the global background field.

[0048] 4) The lack of an adaptive quality control mechanism based on the quality of observation data makes it easy for unreliable profiles to contaminate the inversion results.

[0049] (5) Existing methods are insufficient in terms of numerical stability, and are prone to non-physical solutions or noise amplification when the condition number deteriorates.

[0050] This invention provides a bottom-up inversion method for low ionospheric electron density to solve the above problems.

[0051] Figure 1 The flowchart of a bottom-up inversion method for low ionospheric electron density according to the present invention is shown.

[0052] A bottom-up inversion method for low ionospheric electron density includes the following steps: Step 1: Preprocess the occultation observation data, calculate the tangent point position and height, transform the tangent point coordinates from the geocentric inertial coordinate system (ECI coordinate system) to the WGS-84 geodetic coordinate system, and calculate geometric parameters such as the solar zenith angle.

[0053] The occultation observation data includes data acquired at a high speed of 50Hz by the COSMIC-2 satellite occultation antenna (super-phase). Signal-to-noise ratio (SNR) and cut-off point height Low-rate GNSS ephemeris and LEO ephemeris, UTC time information.

[0054] Low-rate GNSS ephemeris and LEO ephemeris refer to orbital state data with a large sampling interval (below 50Hz) and sparse sampling. They consist of two parts: GNSS ephemeris: orbital parameters of GPS / Galileo / GLONASS navigation satellites (transmitter); LEO ephemeris: orbital parameters of the COSMIC-2 low-Earth orbit satellite (receiver).

[0055] By mapping low-rate ephemeris (GNSS ephemeris and LEO ephemeris) to a 50Hz time grid through cubic interpolation, the LEO satellite position and GNSS satellite position at each 50Hz observation time are obtained, and a common spatiotemporal baseline is established.

[0056] Based on the synchronized LEO and GNSS satellite positions, the instantaneous line-of-sight unit vector is calculated to determine the tangent point position and altitude. The tangent point coordinates are then transformed from the ECI coordinate system to the WGS-84 geodetic coordinate system to calculate geometric parameters such as the solar zenith angle.

[0057] GNSS satellite (transmitter) → passes through Earth's atmosphere / ionosphere → LEO occultation satellite (COSMIC-2, receiver). Electromagnetic rays propagate outside the Earth, and the point closest to the Earth's center is the ray tangent point.

[0058] Figure 2 The radio propagation path between the GNSS satellite (point G) and the LEO satellite (point L) is shown, with the line-of-sight vector d̂ and the tangent point labeled. tangent height Earth's radius The key geometric parameters, such as the tangent height, are explained, illustrating the method for calculating the tangent height.

[0059] The line-of-sight unit vector d̂ (from GNSS to LEO) is calculated using the following formula: , in express LEO satellite position vector at time (Geocenter → LEO satellite). express GNSS satellite position vector at time (Geocenter → GNSS satellite).

[0060] , is the path scalar parameter corresponding to the tangent point in the ray parameter equation, representing the path length from the GNSS-launched satellite along the line of sight to the tangent point.

[0061] Tangent point position .

[0062] Tangent height .

[0063] The position of the tangent point in the ECI coordinate system is transformed to the Earth-fixed coordinate system using the following formula: , , Indicates the position of the tangent point in the ECI coordinate system. Indicates the position of the point of tangency in the Earth-fixed coordinate system; The WGS-84 coordinates (φ, λ, h) of the tangent point in the Earth-fixed coordinate system are obtained through a brief fixed-point iteration calculation: Radius of curvature of the Mao-You circle , Indicates the geodetic latitude of the point of tangency. Indicates the semi-major axis. Represents the square of the first eccentricity; Iterative update formula for geodetic latitude: , The magnitude of the projection of the tangent point position vector onto the equatorial plane; The geodetic height of the tangent point is calculated synchronously each time: .

[0064] Data sources can be extended from COSMIC-2 to other GNSS-RO missions (such as MetOp / GRAS, TerraSAR-X, GRACE, etc.), or multi-source data can be fused to improve coverage.

[0065] Step 2: Within the first height range, perform weighted linear fitting on the superphase to obtain the background phase, and extrapolate the fitted background phase to a higher height.

[0066] Within the altitude range of 50-80km, for superphase A weighted linear fit is performed to obtain a fitting function representing the slowly changing background in region F. The fitted background phase is then extrapolated to a higher level. Step 2 is a bottom-up background estimation.

[0067] Within the altitude range of 50-80km, for superphase A weighted linear fit is performed to obtain a fitting function representing the slowly changing background of region F. a and b here only represent the background phase. With tangent height It is an approximately linear relationship.

[0068] Superphase It consists of two parts: geometric curvature delay and plasma dispersive phase lead. Within the 50-80 km altitude window, the contribution of the F region to both parts varies with tangential altitude. Slowly changing, it can be approximated as a linearly slowly varying term. Combining the two into the background phase: , Indicates geometric curvature delay, This indicates that the plasma dispersion phase is ahead. Indicates the background phase.

[0069] With tangent height For the independent variable to be in superphase We performed a weighted linear fit to obtain the background phase corresponding to different tangent point heights: , W represents the weighting function, and the weights... It increases with altitude.

[0070] Design weight The purpose of increasing weights with altitude is to: ensure the upper part of the 50-80km altitude window is closer to the target zone, better representing the local slope of the background; mitigate the impact of lower boundary noise and outliers on the fitting results; and constrain the extrapolation error to near the upper edge of region E. Employing increasing weights with altitude emphasizes samples closer to the target zone and suppresses the leverage effect of boundaries. The reason for choosing the 50-80km altitude range is that this band is located above the lower boundary of the E / D zone, avoiding near-surface artifacts; at the same time, the far-field contribution of the F zone is manifested as a gentle baseline here.

[0071] The fitted background phase is extrapolated to a distance of over 80 km, and the residual is calculated. , This represents the background phase, and the residual is mainly dominated by the plasma delay in the E / D region.

[0072] Background phase obtained from fitting Extrapolate to higher This approach enables the early separation of the upper F-region slow projection and the stratospheric tail term in the observation domain, revealing high gradient information in the E / D region. Compared to the traditional "Abel inversion followed by filtering" approach, this "bottom-up" background removal method allows the geometric sensitivity kernel to converge faster above the tangent point, significantly compressing the long-tail influence of the F-region and improving the condition number and local resolvability of subsequent inversions.

[0073] The background fitting function can be extended from linear to quadratic or higher-order polynomials to handle the nonlinear variation of the F-region contribution under strong dynamic conditions.

[0074] Step 3: Calculate the remaining phase based on the background phase and superphase, and calculate the total electron content at elevation based on the remaining phase.

[0075] After obtaining the background phase Then, the remaining phase is defined as: .

[0076] The physical meaning of the residual phase can be understood as "the quantity dominated by the local plasma after stripping away the F-region and the slowly varying terms of the troposphere." This quantity is directly related to the electron content: for a single-frequency L1, the plasma refractive index... The first-order approximation gives a linear ratio between the phase delay and the electron content, leading to a scaling transformation that maps the phase to the electron content, resulting in the formula: , Indicates the total electron content at elevation. Indicates the remaining phase. The frequency is GPS L1 (1575.42MHz). This transformation places the observations on the TECU scale, giving the SNR-driven covariance a direct physical meaning.

[0077] Step 4: Arrange the multiple total electron contents at different elevations according to the height of the tangent point to form the total electron contents at different elevations. The matrix is ​​used to obtain the total elevation electron content profile; the quality control of the total elevation electron content profile is performed based on the signal-to-noise ratio (SNR), and unqualified total elevation electron content profiles are removed.

[0078] A series of observations of an occultation event (a GNSS signal being obscured by the Earth) show that the tangent point altitude slowly decreases from high altitude to low altitude, and each tangent point altitude corresponds to a total electron content at an elevation.

[0079] The total electron content at multiple elevations is determined by the height of the tangent point. The arrangement forms the total electron content of the elevation. The matrix is ​​used to obtain the total electron content profile at elevation.

[0080] A height-centered sliding window is used to calculate the robust median baseline and volatility index. Local range R(k) and relative change Q(k) are defined, and an absolute threshold is set. and relative threshold Simultaneously calculate the proportion of negative residuals in systems above 90km. When the threshold value exceeds the limit mentioned above, the profile is removed to prevent unreliable data from entering the inversion process.

[0081] Specifically, a fixed-size sliding window is used. R(k) is used to determine the range of values ​​within the sliding window, and Q(k) is used to determine the range of values ​​between windows after averaging the values ​​within the window. When R(k) exceeds... Or Q(k) exceeds When that happens, the section should be removed.

[0082] Quality control related thresholds , and It can adaptively adjust according to different application scenarios (such as calm / disturbed periods, day / night).

[0083] Step 5: Construct the Abel weighted forward model.

[0084] The total electron content at multiple elevations is determined by the height of the tangent point. Arrangement The matrix yields the total electron content profile at elevation. When the tangent point height is approximately [value missing], (Total electron content profile at elevation) and target profile The relationship between (the unknowns to be inverted) follows a discrete Abel integral.

[0085] Construct a discrete Abel-type integral forward model on a uniform height grid: with a step size of Uniform height grid Above, order and ,in, This represents the distance from the current location to the Earth's center. express The distance from the location to the Earth's center, Re, represents the Earth's radius. To improve accuracy, physics guidance is embedded, multiplying each column of the forward matrix (discrete Abel-type integral forward model) by the scene matching factor. This yields the Abel weighted forward model: , Indicates weight, express ; ; in The scene matching background provided for the FIRI-2018 empirical model, the FIRI weighted model focuses the sensitivity on the E-region peak of 90-120km, suppresses upper atmospheric tail leakage, and improves the matrix condition number; The vector form of the Abel weighted forward model is: , in ε represents the measurement noise.

[0086] FIRI weighting can be replaced by other empirical ionospheric models (such as IRI-2016) or data assimilation products, as long as the model can provide a reasonable background peak structure in the E region.

[0087] The Abel weighted forward model is a forward equation, given the target profile. Integrate to calculate the ray hTEC; inversion is performed given hTEC, and the solution is obtained. .

[0088] Step 6: Perform optimal estimation inversion based on the total electron content profile of elevation to obtain the inverted electron density.

[0089] First, construct the observation covariance matrix and the prior covariance, and then calculate the optimal estimated solution.

[0090] Construct the observation covariance matrix ,in , This represents the lower limit of the signal-to-noise ratio; the observation covariance is derived from the SNR and a lower limit is set to avoid extreme values.

[0091] The short correlation index kernel is used as the prior covariance: .

[0092] With zero prior mean, the formula for calculating the optimal estimated solution is: , Represents electron density, hTEC indicates , , Represent the covariance matrix; When calculating, the above formula should be written as: .

[0093] The matrix condition number cond(A) is defined as the ratio of the largest singular value to the smallest singular value of a matrix, and is used to characterize a system of linear equations. The degree of ill-conditioning is considered. A larger condition number significantly amplifies even small noises in the observation vector, causing the solution to deviate significantly from the true value. When the condition number cond(A) exceeds a threshold, an adaptive diagonal loading term λI is added to reduce noise amplification. Non-negativity constraints are applied when necessary to ensure physical feasibility.

[0094] Step 7: Separate the inverted electron density from the background and the layered enhancement.

[0095] Inversion electron density obtained It also includes a slowly changing background field and local enhancements caused by thin ionospheric layers (such as the Es layer). To adaptively separate these two components without empirical threshold constraints, a data-driven strategy based on moving average and iterative annealing is adopted.

[0096] right A moving average was performed along the height grid (window W=31) to obtain the initial long-wave background Ne(0). This window corresponds to a height half-width of several kilometers, which effectively suppresses sharp structures at the 2–10 km scale while preserving the background morphology of region E.

[0097] The initial long-wave background is subjected to a certain number of iterative annealing processes to ensure that the samples meet the quality requirements: Define residual High-value samples that meet this condition are marked as positive outliers (candidate layer enhancement), temporarily set as missing, and then filled in and smoothed by linear interpolation to obtain... .

[0098] After the missing segments are completed using linear interpolation, window smoothing is performed again to obtain the updated background estimate Ne. (i+1) (h). Iteration continues until no new outliers are detected or the maximum number of iterations is reached. The final separated long-wavelength background field is: The corresponding layered reinforcement component is .

[0099] Iterative annealing can be replaced by other adaptive filtering methods, such as wavelet decomposition and variational mode decomposition.

[0100] Optional, step 8, constructing a global context field.

[0101] Background values ​​were sampled at 100 km in each total electron content profile obtained through quality control at elevation. A global background electron density map of the low ionosphere is generated by mapping the data onto a 0.5°×0.5° global grid using natural neighborhood interpolation combined with mild Gaussian smoothing. The global grid resolution can be varied depending on the application's spatial resolution requirements and the amount of available data.

[0102] In the above steps, steps 2 to 7 constitute the core technical solution of this invention. The key difference from existing technologies lies in its removal of the F-region background during the inversion process through "bottom-up" parameter processing. Before optimal estimation, the FIRI physical model is embedded into an Abel weighted matrix instead of being used as a prior mean, which improves the condition number without injecting model bias. Furthermore, the combination of an SNR-driven adaptive error model and short-correlation priors achieves data-driven robust inversion.

[0103] Example of single-profile electron density inversion: Taking an occultation observation by the COSMIC-2 satellite on January 3, 2021 as an example, this illustrates the complete processing flow from raw data to electron density profile.

[0104] Step 1: Preprocess the occultation observation data, calculate the tangent point position and height, transform the tangent point coordinates from the geocentric inertial coordinate system (ECI coordinate system) to the WGS-84 geodetic coordinate system, and calculate geometric parameters such as the solar zenith angle.

[0105] Read the atmPhs file of this occultation and obtain the 50Hz superphase of the L1 channel. SNR and tangent height Data. Simultaneously, the corresponding POD ephemeris file is read to obtain the position and velocity information of GPS and LEO satellites.

[0106] The line-of-sight unit vector is calculated, the tangent point position and tangent point height are calculated, the tangent point coordinates in the ECI coordinate system are transformed to the Earth-Fixed coordinate system, and then transformed to the WGS-84 coordinate system.

[0107] Step 2, within the altitude range of 50-80km, perform superphase testing. A weighted linear fit is performed to obtain the background phase, and the fitted background phase is extrapolated to a higher height.

[0108] Figure 3 The superphase and background fitting plot for case N1 illustrates the superphase. The monotonic increasing trend with tangent height, and the linear background phase within the 50-80km range. The fitting effect is shown on the right, and the corresponding signal-to-noise ratio (SNR) profile is shown on the right, reflecting the validity of the measurement data. Figure 4 The superphase and background fitting plot for Case N2 demonstrates the results under relatively ideal observation conditions. and The relationship diagram, and the SNR level of a high-quality signal. Figure 5 The superphase and background fitting plot for Case N3 demonstrates the effectiveness of background fitting under moderate SNR fluctuations, illustrating the method's tolerance to stationary noise. Figure 6The superphase and background fitting plot for Case A1 shows the SNR jump at 100-110km and The abnormal situation illustrates how the quality control mechanism identifies and eliminates unreliable profiles.

[0109] Figures 3-5 The images shown are relatively successful background fitting images, among which Figure 5 The image shows obvious continuous spikes in the SNR image due to continuous interference with the received signal, but these spikes are still retained because the amplitude does not exceed the quality control threshold. Figure 6 This demonstrates that mutations in the E region cause superphase. An abnormal change occurred, so subsequent calculations will be terminated during the quality control process.

[0110] Step 3: Calculate the remaining phase based on the background phase and superphase, and calculate the total electron content at elevation based on the remaining phase.

[0111] Figure 7 The hTEC profile of Case N1 shows the height-limited electron content that rises smoothly from 60 to 120 km, with the steepest slope located between 80 and 110 km, verifying the electron accumulation characteristics of the E region. Figure 8 The hTEC profile of Case N2 shows the gradual hTEC change under weak E region conditions, serving as an amplitude consistency check for background cases. Figure 9 The hTEC profile of Case N3 shows the S-shaped curve and slope variation, demonstrating the presence of abrupt changes in the thin layer or E region. Figure 10 The hTEC profile of Case A1 shows the sharp rise and shoulder-like features produced when the E region has a strong mutation. Quality control requires the rejection of such unreliable data.

[0112] Figures 7 to 10 This section illustrates different scenarios for the results of this step. Cases N1 and N3 represent the most common and reasonable variations. N3 indicates that slight interference with the received signal caused image glitches. Case N2 represents regions with low electron density, resulting in a slower growth rate and a higher layer height than usual, only beginning to rise significantly at 100 km. This phenomenon mainly exists in regions with relatively calm ionospheric activity, such as the polar night. Case A1 represents a signal subjected to strong E-region abrupt interference; quality control requires the exclusion of cases exhibiting this phenomenon.

[0113] Step 5: Construct the Abel weighted forward model.

[0114] Step 6: Perform optimal estimation inversion based on the total electron content profile of elevation to obtain the inverted electron density.

[0115] Figures 11-13 This is a triptych of the inversion chain (Case ac). Each diagram contains three sub-diagrams: the top left is... , The image shows the SNR; the top right is the hTEC profile; the bottom shows a comparison of the inverted electron density Ne (green), long-wave background (blue), and FIRI-2018 reference (black dashed line), demonstrating the complete chain from the original observations to the final inversion results. Figures 11-13 The three cases presented are close to ideal conditions. The green curve shows the initial Ne versus the final background electron density. The comparison shows that the above smoothing measures are very effective. The main reason for the Ne jump at higher altitudes is that the Abel core increases rapidly with decreasing altitude, thus amplifying small changes in hTEC; any fine structure or noise entering this high-gain frequency band will produce slight oscillations.

[0116] Figure 14 A 100km global background electron density distribution map is shown according to an embodiment of the present invention. Batch processing and global mapping embodiment.

[0117] All COSMIC-2 occultation data from January 3, 2021 (DOY 003) were batch-processed. A total of 932 profiles were generated; after quality control, 373 (40.0%) passed, and 559 were rejected. The main reasons for rejection were: interference above 90km (290 profiles, 31.1%) and SNR anomalies (236 profiles, 25.3%).

[0118] Natural neighborhood interpolation was performed on all 100km background values ​​passing through the profile, with a grid resolution of 0.5° × 0.5°. After interpolation, a light Gaussian smoothing was applied to generate a global background electron density distribution map. Figure 14 ).

[0119] Figure 14 This shows the daytime conditions in low and mid-latitude regions in early January 2021. Distribution. From Figure 14 It can be seen that the region with a small solar zenith angle (near the equator at noon) shows continuous enhancement, while the twilight sector shows weaker enhancement. The values ​​in the low latitudes of the Southern Hemisphere (summer) are slightly higher than those in the Northern Hemisphere (winter), consistent with climatological expectations. The patches in the figure represent small-amplitude E-layer abrupt events.

[0120] Example of handling mutation events at layer E.

[0121] When the research target is E-layer mutation (Es) events, improved processing parameters are used: the SNR gating threshold is relaxed, the prior correlation length is shortened to L_corr, and a bandpass filter is applied to hTEC to isolate thin-layer scales. The focus after inversion is... Components, quantization layer thickness (half-height, full width at half maximum), and frequency of occurrence.

[0122] Taking case A1 as an example, it exhibits a highly distinctive E-layer mutation phenomenon. Through the measures described above, this type of case was preserved. For example, event A1 was detected at 105 km. Peak value reaches 8×10¹ 0 e·m⁻³, with a layer thickness of approximately 5km. This can be seen from the image of case A1 ( Figure 6 , Figure 10 It can be seen that it has typical U-shaped change characteristics, which is consistent with the Es characteristics proposed by Zeng and Sokolovskiy (2010).

[0123] The low-ionospheric electron density inversion method of the present invention has the following innovative features: Bottom-up F-region background removal technology at the measurement level: Unlike traditional methods that handle F-region interference during the inversion stage, this invention removes the slowly changing F-region background during the measurement preprocessing stage through 50-80km weighted linear fitting, making the residuals mainly dominated by E / D region signals. This "bottom-up" strategy fundamentally weakens the influence of the far-field F-region, significantly sharpens the geometric weighting kernel of the E-region, and improves the resolution and accuracy of low-altitude inversion.

[0124] The FIRI physics model embedding of an Abel weighted matrix: This invention innovatively uses the background electron density of the FIRI-2018 empirical model as the column weight of the physics-guided embedding Abel forward matrix, rather than imposing it as a prior mean on the solution. This design improves the matrix condition number (spectral tightening, reduced condition number), focuses sensitivity on the 90-120km target layer, and avoids direct dependence on model amplitude, maintaining the data-driven nature.

[0125] Joint optimization of SNR-driven adaptive observation covariance and short correlation prior: This invention establishes an observation error model based on signal-to-noise ratio, enabling low SNR data to automatically receive smaller weights; simultaneously, a 5km short correlation index kernel is used as a smoothing prior, preserving kilometer-scale structure while suppressing high-frequency noise. The combination of these two approaches achieves robust inversion with adaptive data quality.

[0126] Iterative annealing-based adaptive separation method for background-layered enhancement: This invention combines 31-point moving average with iterative outlier suppression to adaptively separate long-wavelength background from layered enhancements (such as sporadic E-layers) on a 2-10km scale without the need for a pre-set fixed template. This method operates stably under various local conditions, avoiding the sensitivity of traditional thresholding methods to layer amplitude and location.

[0127] Multi-layered quality control and feedback mechanism: This invention establishes a full-chain quality control system from preprocessing to postprocessing, including profile screening based on SNR local variability, verification of the proportion of negative residuals above 90km, and monitoring of conditional numbers after inversion. Combined with traceable log recording and explicit parameter declaration, the repeatability and auditability of the results are ensured.

[0128] While some embodiments of the present invention have been described in this application, those skilled in the art will understand that these embodiments are merely illustrative. Numerous variations, alternatives, and improvements will arise in those skilled in the art under the teachings of this invention without departing from its scope. The appended claims are intended to define the scope of the invention and thereby cover methods and structures within the scope of the claims themselves and their equivalents.

Claims

1. A bottom-up inversion method for low ionospheric electron density, characterized in that, include: Within the first altitude range, a weighted linear fit is performed on the superphase to obtain the background phase, and the fitted background phase is extrapolated to a higher altitude. The remaining phase is calculated based on the background phase and superphase, and the total electron content at elevation is calculated based on the remaining phase. The total electron content at multiple elevations is arranged according to the height of the tangent point to form the total electron content at elevation. The matrix is ​​used to obtain the total electron content profile at elevation. Construct an Abel-weighted forward model; Optimal estimation inversion was performed based on the total electron content profile at elevation to obtain the inverted electron density; and Background and layered enhancement separation are performed on the inverted electron density.

2. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Before constructing the Abel weighted forward model, the following steps are also included: quality control of the total elevation electron content profile, removing unqualified total elevation electron content profiles.

3. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, The first altitude range is 50-80km; Background phase fitting function ; With tangent height For the independent variable to be in superphase We performed a weighted linear fit to obtain the background phase corresponding to different tangent point heights: ， W represents the weight function, and the weights... As altitude increases.

4. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, The remaining phase is calculated based on the background phase and superphase, and the total elevation electron content is calculated based on the remaining phase, including: The remaining phase is defined as: , Indicates background phase, Indicates superphase; The relationship between residual phase and total electron content at elevation: ， Indicates the total electron content at elevation. Indicates the remaining phase. It is the GPS L1 frequency, 1575.42MHz.

5. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Constructing an Abel-weighted forward model includes: When the tangent point heights are approximately the same, the elevation total electron content profile and the target profile are similar. The relationship follows a discrete Abel integral; Construct a discrete Abel-type integral forward model on a uniform height grid: with a step size of Uniform height grid Above, order and ,in This represents the distance from the current location to the Earth's center. Re represents the distance from the Earth's center to the position at time t=i, and Re represents the Earth's radius; Multiply each column of the discrete Abel integral forward model by This yields the Abel weighted forward model: , Represents the weighting function. express ; ， in Scene matching background provided for the FIRI-2018 empirical model; The vector form of the Abel weighted forward model is: , in ε represents the measurement noise.

6. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Based on the elevation total electron content profile, optimal estimation inversion is performed to obtain the inverted electron density, including: Construct the observation covariance matrix: ,in , SNR represents the lower limit of signal-to-noise ratio; The short correlation index kernel is used as the prior covariance: ; The formula for calculating the optimal estimate of electron density under zero prior mean is: ， Represents electron density, hTEC indicates , , Represent the covariance matrix; When calculating, the above formula should be written as: 。 7. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Background and layered enhancement separation of the inverted electron density includes: The inverted electron density is averaged along the height grid to obtain the initial long-wave background Ne(0); Iterative annealing of the initial long-wave background yields the long-wave background field. ; The corresponding layered reinforcement component is , This is for inverting electron density.

8. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Also includes: Constructing a global context field, including: Background values ​​were sampled at 100 km in each total electron content profile obtained through quality control at elevation. Using natural neighborhood interpolation combined with mild Gaussian smoothing, the data is mapped onto a 0.5°×0.5° global grid to generate a global background electron density distribution map of the low ionosphere.

9. The method for bottom-up inversion of low ionospheric electron density according to claim 1, characterized in that, Also includes: The occultation observation data is preprocessed to calculate the tangent point position and height, and the tangent point coordinates are transformed from the geocentric inertial coordinate system to the WGS-84 geodetic coordinate system to calculate the solar zenith angle. Occultation observation data includes superphase Signal-to-noise ratio (SNR) and cut-off point height GNSS ephemeris, LEO ephemeris, and UTC time information.

Images