A thermal spherical density perturbation inversion method based on orbital drag response parameter difference

Abstract

This invention belongs to the field of atmospheric density inversion technology, specifically relating to a thermospheric density perturbation inversion method based on orbital drag response parameter difference. The steps include: for a low Earth orbit target, acquiring the target's orbital data product; constructing candidate target pairs, calculating orbital similarity scores, and obtaining a set of valid target pairs; propagating the targets in the valid target pairs to a preset unified epoch to obtain orbital geometric information; calculating the logarithmic drag difference component to obtain the baseline correction differential perturbation amount and acquiring robust dispersion, retaining target pairs with robust dispersion less than a preset dispersion threshold; inverting the relative thermospheric density perturbation amount at the midpoint of the target pair; generating a thermospheric density perturbation field at a preset region and preset time resolution, and outputting it according to a preset format. This invention relies on publicly available TLE data and achieves thermospheric density perturbation inversion without the need for dedicated equipment or single-satellite calibration.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric density inversion technology, specifically relating to a thermal density perturbation inversion method based on the difference of orbital drag response parameters. Background Technology

[0002] Thermospheric density (thermosphere density) in the low Earth orbit (LEO) altitude range of 200–700 km is a core parameter for space weather monitoring, LEO spacecraft orbit forecasting, and space debris collision early warning. Geomagnetic disturbances and changes in solar radiation can cause short-term, drastic fluctuations in thermospheric density. Abnormal atmospheric density will significantly exacerbate changes in atmospheric drag, directly causing inaccurate satellite orbit decay calculations and increasing the risk of debris collisions. Therefore, accurate and large-scale real-time inversion of thermospheric density disturbances is an essential technical requirement in the fields of aerospace operations and space weather.

[0003] Current methods for obtaining thermosphere density mainly fall into the following three categories: Empirical model methods, such as NRLMSISE-00, JB2008, and DTM2013, are based on historical observations. These methods are highly efficient, but they exhibit significant lag in response to geomagnetic events, resulting in large errors and making them unsuitable for real-time, high-precision applications. Precise orbit determination (POD): Based on the inversion density of orbital perturbations measured by an on-orbit satellite GNSS receiver, a typical example is the ESA Swarm mission. This method offers high accuracy but requires on-orbit equipment and has limitations in spatial coverage and temporal resolution. Single-target TLE drag analysis: Extracting orbital decay rate and inverting density from the two-line orbital elements (TLE) time series of a single target. (TLE...) These are empirical drag parameters used in the SGP4 fitting process, which incorporate multiple factors including actual atmospheric density, satellite attitude, surface-to-mass ratio, orbit fitting arc, maneuvering, and orbit determination error; therefore, based on a single target... When directly inverting absolute atmospheric density, the inversion results are greatly affected by the systematic bias of the target's own characteristics.

[0004] Current monitoring of thermospheric density disturbances faces the following challenges: First, insufficient spatial coverage, making it difficult to achieve continuous monitoring across a wide range and multiple altitudes using only a small number of satellites equipped with dedicated payloads; second, limited real-time performance, with traditional empirical models exhibiting lag in response to short-term events such as geomagnetic disturbances and changes in solar radiation; and third, the inversion of single-target orbital drag is easily affected by target attitude, surface-to-mass ratio, drag coefficient, orbital fitting error, and maneuvering behavior, making it difficult to distinguish between individual target differences and actual atmospheric density disturbances.

[0005] The shortcomings of the existing three types of technologies have limited the implementation of low-cost, wide-coverage, and high-time-efficiency thermospheric density perturbation detection. There is an urgent need for a technical solution that does not rely on satellite-specific hardware payloads, can eliminate the inherent parameter errors of the satellite itself, and can achieve full-domain density perturbation inversion based on massive amounts of publicly available orbital data. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the purpose of this invention is to provide a thermosphere density perturbation inversion method based on orbital drag response parameter difference. Relying on publicly available TLE data, the method achieves thermosphere density perturbation inversion by pairing and differencing targets in the same shell and orbital plane, static solar robustness correction and confidence weighted fusion, without the need for special equipment and single-star calibration.

[0007] To achieve the above objectives, this invention provides a method for inverting thermal layer density perturbation based on the difference in orbital drag response parameters, comprising the following steps: S1. For targets in low Earth orbit, acquire the target's orbital data products, analyze them to obtain the target identifier, epoch, orbital state parameters and drag response parameters, and perform quality control to obtain a valid set of low Earth orbit targets. S2. Determine the orbital shell, orbital plane, and in-plane phase angle of the effective low-Earth orbit target based on the orbital state parameters. Within the same orbital shell and the same orbital plane, construct candidate target pairs based on the in-plane phase difference and orbital altitude difference, calculate the orbital similarity score, and obtain the set of effective target pairs. S3. Propagate the target of the effective target alignment to the preset unified epoch to obtain orbital geometric information including target position, altitude, velocity, relative atmospheric velocity and target alignment point position; S4. Calculate the logarithmic drag difference component based on the drag response parameters of the two targets in the target pair, determine the target pair baseline using robust statistics within the static daily reference period, obtain the baseline correction differential disturbance, obtain the robust dispersion using the logarithmic drag difference component, and retain target pairs with robust dispersion less than the preset dispersion threshold. S5. Based on the target midpoint position, altitude, and unified epoch, obtain the background thermosphere density using a preset thermosphere atmospheric model. Based on the baseline correction differential perturbation, orbital geometry information, and background thermosphere density, establish the mapping relationship between the baseline correction differential perturbation and the relative perturbation of thermosphere density, and invert to obtain the relative perturbation of thermosphere density at the target midpoint position. S6. Perform spatiotemporal fusion of the relative perturbations of thermal layer density for multiple target pairs to generate thermal layer density perturbation fields under preset regions and preset time resolutions, and output them in a preset format.

[0008] In a preferred embodiment of the present invention, in S1, the track data product is the target TLE data, and the drag response parameter is the TLE empirical drag parameter. ; The orbital state parameters include orbital elements, position and velocity parameters, and orbital geometry parameters. The orbital elements include the semi-major axis, orbital inclination, eccentricity, right ascension of the ascending node, argument of perigee, mean perigee, and mean motion. The position and velocity parameters include the target's position and velocity at the corresponding epoch. The orbital geometry parameters include orbital altitude, orbital plane markings, and in-plane phase angles.

[0009] As a preferred embodiment of the present invention, in S1, the quality control includes: removing targets with abnormal trajectories, suspected maneuvering targets, and targets with missing data, as well as removing... Invalid or Less than preset Threshold records; Orbital anomaly targets include targets whose orbital altitude is outside the preset altitude range, whose eccentricity is greater than the preset eccentricity threshold, or whose orbital elements exceed the preset physical reasonable range; Missing data targets include missing target identifiers, epochs, orbital elements, or... The goal is to obtain at least one of the following data; Invalid includes It can be a null value, a non-numeric value, or equal to zero; The criteria for identifying a suspected maneuvering target are as follows: a target is identified as a suspected maneuvering target if at least one of the following conditions is met: The target is that the change in the semi-major axis between two consecutive TLE updates is greater than a preset semi-major axis change threshold; The change in average motion between adjacent epochs is greater than the preset threshold for average motion change. The change in adjacent epochs is greater than the preset value. Change threshold It is a positive constant; At least one of the following parameters—orbital altitude, orbital inclination, eccentricity, right ascension of the ascending node, or mean motion—experiences an anomalous jump between adjacent epochs; The TLE data time interval exceeds the preset time interval threshold.

[0010] As a preferred embodiment of the present invention, in S2, constructing candidate target pairs specifically involves: selecting targets with orbital altitudes ranging from 200km to 700km from the effective low-Earth orbit target set; determining the orbital shell, orbital plane, and in-plane phase angle of each target based on orbital state parameters; and pairing them up in pairs according to preset pairing similarity constraints within the same orbital shell and the same orbital plane to form candidate target pairs. The preset pairing similarity constraints include: the phase difference between the two targets in the same target pair in the orbital plane is less than a preset phase difference threshold, and the orbital height difference is less than a preset height difference threshold; and at least one of the following constraints is satisfied: the orbital inclination difference between the two targets in the same target pair is less than a preset inclination difference threshold, and the relative velocity difference is less than a preset velocity difference threshold.

[0011] As a preferred embodiment of the present invention, in S2, for the target pair The orbital similarity score is calculated as follows: ; In the formula, This represents the orbital similarity score between target i and target j; This represents the difference in orbital altitude between target i and target j; This represents the difference in orbital inclination angle between target i and target j; This represents the difference in right ascension between the ascending nodes of target i and target j; This represents the epoch time difference between the orbital data products of target i and target j; They are respectively , , , The corresponding normalized scaling parameter; They are respectively , , , The corresponding weighting coefficients; Calculate orbital similarity scores for each candidate target pair, sort them from highest to lowest orbital similarity scores, and retain the top N candidate target pairs as the set of effective low-Earth orbit target pairs.

[0012] As a preferred embodiment of the present invention, in S3, the SGP4 orbital propagation model is used to propagate each target to a preset unified epoch. Before propagation, the drag response parameters of the target are time-smoothed. The time-smoothing process adopts a fixed-length sliding window mid-range filter, and the sliding window length is 3 hours to 24 hours.

[0013] As a preferred embodiment of the present invention, in S4, for the target pair The logarithmic resistance difference component is calculated according to the following formula: ; In the formula, This represents the logarithmic difference in drag between target i and target j at time t; This indicates that target i at time t ; This indicates that target j at time t ; The baseline correction differential perturbation is calculated as follows: ; In the formula, This represents the baseline correction differential perturbation between target i and target j at time t; Indicates the target relative to the baseline, using the static day reference time period. the median; based on To obtain robust dispersion: ; In the formula, MAD represents the robust dispersion of the logarithmic resistance difference component between target i and target j over a static daily reference period; MAD represents the median absolute deviation. Indicates the reference time period for quiet days; when If the dispersion is less than the preset dispersion threshold, the corresponding target pair is retained; otherwise, the corresponding target pair is discarded.

[0014] As a preferred embodiment of the present invention, in S5, the preset thermospheric model is an empirical atmospheric model, a physical atmospheric model, or a combination of both, and the target pair is obtained based on the preset thermospheric model. The midpoint of the background thermal layer density at time t ; Will Mapped to the relative perturbation of thermal layer density: ; In the formula, Indicates target pair The relative perturbation of the thermal layer density at the midpoint of the position at time t; Indicates target pair The differential perturbation mapping coefficients.

[0015] As a preferred embodiment of the present invention The robust linear regression relationship between the logarithmic density perturbation and the baseline correction difference perturbation, as determined by the empirical model, is as follows: ; In the formula, Indicates target pair The logarithmic density perturbation of the empirical model at time t; For the target The bias term; For the target The residual term at time t represents the difference between the logarithmic density perturbation of the empirical model and the predicted value of the linear regression. Calculate using the following formula: ; In the formula, median represents the median; Calculate the reference time period for quiet days The median absolute deviation is used as the residual robust dispersion. When the residual robust dispersion is less than a preset residual threshold, the corresponding differential perturbation mapping coefficient is retained. .

[0016] As a preferred embodiment of the present invention, after obtaining the relative perturbation of the thermal layer density, the density enhancement factor, the model-constrained density estimate, and the normalized density perturbation index are calculated based on the relative perturbation of the thermal layer density. The density enhancement factor is calculated as follows: ; In the formula, Indicates target pair Density enhancement factor at time t; The model constraint density estimate is calculated as follows: ; In the formula, Indicates target pair The estimated value of the model constraint density at the midpoint of the model at time t; The normalized density perturbation index is calculated as follows: ; In the formula, Indicates target pair The normalized density perturbation exponent at time t; This indicates the preset lower limit of dispersion.

[0017] As a preferred embodiment of the present invention, in S6, for the relative perturbation of thermal layer density of multiple target pairs falling into the same preset spatiotemporal grid g, a robust weighted fusion method is used for spatiotemporal fusion, and the fused relative perturbation of thermal layer density is... Represented as: ; In the formula, wmedian represents the weighted median; Indicates target pair The fusion weight at time t.

[0018] As a preferred embodiment of the present invention, the confidence index of the spatiotemporal grid g is calculated. : ; In the formula, This represents the number of target pairs participating in the fusion within the spatiotemporal grid g; Indicates the target's quantitative scale parameter; denoted as the median absolute deviation of the relative perturbation of the thermal layer density within the spatiotemporal grid g; is the average orbital similarity score of the target pairs within the spatiotemporal grid g; Will Normalized to the range of 0 to 1, the result is compared with a preset confidence threshold. When the normalized value is greater than the confidence threshold, the fusion result of the spatiotemporal grid g is retained.

[0019] The beneficial effects of this invention are: This invention employs a core architecture of target pairing within the same orbital shell and on the same orbital surface, logarithmic drag difference component differentiation, and robust correction based on a static solar reference. Through pairing and differentiation, it effectively eliminates systematic deviations caused by inherent characteristics such as the target surface mass ratio and drag coefficient, eliminating the need for individual calibration of the ballistic parameters of each satellite. The data source is publicly available TLE orbital data of low-Earth orbit targets, without relying on dedicated on-orbit equipment such as accelerometers and GNSS on satellites. Relying on a massive number of low-Earth orbit constellation targets, it can achieve thermosphere density perturbation inversion with full coverage at altitudes from 200km to 700km, resulting in low data acquisition costs and outstanding spatial coverage capabilities.

[0020] This invention employs a robust statistical and confidence-weighted fusion mechanism, achieving robust handling of noise and outliers based on median absolute deviation. It utilizes a spatiotemporal grid confidence index to quantitatively evaluate and screen the results, effectively improving the reliability of the inversion results. This mechanism is highly compatible, incorporating not only TLE empirical drag parameters but also various drag-related observational data such as semi-major axis attenuation rate and non-conservative acceleration. Combined with multiple data fusion algorithms, it further enhances inversion accuracy. The inversion products can generate diverse products such as atmospheric three-dimensional density reconstruction, orbital lifetime correction, and space weather anomaly warnings. It can quickly capture short-term anomalies in thermosphere density caused by geomagnetic disturbances, demonstrating excellent practicality and scalability. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the principle of the method of the present invention; Figure 2 This is a schematic diagram of the time series of the group resistance disturbance surrogate index in Example 1; Figure 3 This is a schematic diagram of the time series of the grid median perturbation rate in Example 1; Figure 4 This is a time box plot showing the median perturbation ratio distribution of the grid during different periods of the quiet day, the perturbation period, and the recovery period in Example 1; Figure 5 This is a spatial representation of the grid median perturbation ratio in Example 1. Detailed Implementation

[0022] The embodiments of the present invention will be further described below with reference to the accompanying drawings: Example 1: As Figure 1 As shown, a method for inverting thermal spherical density perturbation based on the difference in orbital drag response parameters includes the following steps: S1. For targets in low Earth orbit, acquire the target's orbital data products, analyze them to obtain the target identifier, epoch, orbital state parameters and drag response parameters, and perform quality control to obtain a valid set of low Earth orbit targets. S2. Determine the orbital shell, orbital plane, and in-plane phase angle of the effective low-Earth orbit target based on the orbital state parameters. Within the same orbital shell and the same orbital plane, construct candidate target pairs based on the in-plane phase difference and orbital altitude difference, calculate the orbital similarity score, and obtain the set of effective target pairs. S3. Propagate the target of the effective target alignment to the preset unified epoch to obtain orbital geometric information including target position, altitude, velocity, relative atmospheric velocity and target alignment point position; S4. Calculate the logarithmic drag difference component based on the drag response parameters of the two targets in the target pair, determine the target pair baseline using robust statistics within the static daily reference period, obtain the baseline correction differential disturbance, obtain the robust dispersion using the logarithmic drag difference component, and retain target pairs with robust dispersion less than the preset dispersion threshold. S5. Based on the target midpoint position, altitude, and unified epoch, obtain the background thermosphere density using a preset thermosphere atmospheric model. Based on the baseline correction differential perturbation, orbital geometry information, and background thermosphere density, establish the mapping relationship between the baseline correction differential perturbation and the relative perturbation of thermosphere density, and invert to obtain the relative perturbation of thermosphere density at the target midpoint position. S6. Perform spatiotemporal fusion of the relative perturbations of thermal layer density for multiple target pairs to generate thermal layer density perturbation fields under preset regions and preset time resolutions, and output them in a preset format.

[0023] In S1, the orbital data product is TLE data for the target. The drag response parameter is a parameter characterizing the orbital energy decay, orbital decay change, or orbital fitting drag correction caused by atmospheric drag on a low Earth orbit target. For TLE data, the drag response parameter is the TLE empirical drag parameter. ; The orbital state parameters include orbital elements, position and velocity parameters, and orbital geometry parameters, which are obtained directly or calculated. The orbital elements include the semi-major axis, orbital inclination, eccentricity, right ascension of the ascending node, argument of perigee, mean perigee, and mean motion. The position and velocity parameters include the target's position and velocity at the corresponding epoch. The orbital geometry parameters include orbital altitude, orbital plane markings, and in-plane phase angles.

[0024] Quality control includes: removing targets with abnormal trajectories, suspected maneuvering targets, targets with missing data, and removing... Invalid or Less than preset Threshold records; Orbital anomaly targets include targets whose orbital altitude is outside the preset altitude range, whose eccentricity is greater than the preset eccentricity threshold, or whose orbital elements exceed the preset physical reasonable range; Missing data targets include missing target identifiers, epochs, orbital elements, or... The goal is to obtain at least one of the following data; Invalid includes It can be a null value, a non-numeric value, or equal to zero; The criteria for identifying a suspected maneuvering target are as follows: a target is identified as a suspected maneuvering target if at least one of the following conditions is met: The target is that the change in the semi-major axis between two consecutive TLE updates is greater than a preset semi-major axis change threshold; The change in average motion between adjacent epochs is greater than the preset threshold for average motion change. The change in adjacent epochs is greater than the preset value. Change threshold It is a positive constant to prevent singularities (arguments are 0) in logarithmic calculations. At least one of the following parameters—orbital altitude, orbital inclination, eccentricity, right ascension of the ascending node, or mean motion—experiences an anomalous jump between adjacent epochs; the TLE data time interval exceeds a preset time interval threshold.

[0025] Among them, between two adjacent observation epochs, the values ​​of any orbital parameter such as orbital altitude, orbital inclination, eccentricity, right ascension of the ascending node, and mean motion deviate from the natural evolution law of the orbit, showing a large amplitude and abrupt jump, which is different from the gradual change brought about by atmospheric drag and celestial perturbation, and is called an anomalous jump.

[0026] In S2, the construction of candidate target pairs specifically involves: selecting targets with orbital altitudes ranging from 200km to 700km from the effective set of low-Earth orbit targets; determining the orbital shell, orbital plane, and in-plane phase angle of each target based on orbital state parameters; pairing targets in pairs within the same orbital shell and orbital plane according to preset pairing similarity constraints to form candidate target pairs; the orbital shell is a spatial layer divided according to orbital altitude, and the orbital altitudes of all low-Earth orbit targets within the same shell are within the same preset range; The preset pairing similarity constraints include: the phase difference between the two targets in the same target pair in the orbital plane is less than a preset phase difference threshold, and the orbital height difference is less than a preset height difference threshold; and at least one of the following constraints is satisfied: the orbital inclination difference between the two targets in the same target pair is less than a preset inclination difference threshold, and the relative velocity difference is less than a preset velocity difference threshold.

[0027] For target pair The orbital similarity score is calculated as follows: ; In the formula, This represents the orbital similarity score between target i and target j; This represents the difference in orbital altitude between target i and target j; This represents the difference in orbital inclination angle between target i and target j; This represents the difference in right ascension between the ascending nodes of target i and target j; This represents the epoch time difference between the orbital data products of target i and target j; They are respectively , , , The corresponding normalized scaling parameter; They are respectively , , , The corresponding weighting coefficients, and ; Calculate orbital similarity scores for each candidate target pair, sort them from highest to lowest orbital similarity scores, and retain the top N candidate target pairs as the set of effective low-Earth orbit target pairs.

[0028] In S3, the SGP4 orbital propagation model is used to propagate each target to a preset unified epoch. Before propagation, the drag response parameters of the target are smoothed over time. The time smoothing process uses a fixed-length sliding window with a median filter, and the sliding window length is 3 to 24 hours.

[0029] A sliding window mean filter is used for time smoothing to remove short-term noise, outliers, and random fluctuations in the drag response parameters, and to suppress interference from random errors and instantaneous anomalies in the TLE data. At the same time, the true trend of parameter changes is preserved, which improves the stability of subsequent calculations of logarithmic drag difference components and differential perturbation and the reliability of inversion results.

[0030] In S4, for the target pair The logarithmic resistance difference component is calculated according to the following formula: ; In the formula, This represents the logarithmic difference in drag between target i and target j at time t; This indicates that target i at time t ; This indicates that target j at time t ; The baseline correction differential perturbation is calculated as follows: ; In the formula, This represents the baseline correction differential perturbation between target i and target j at time t; Indicates the target relative to the baseline, using the static day reference time period. The median, i.e. median means taking the median. Indicates the reference time period for quiet days; based on To obtain robust dispersion: ; In the formula, The MAD (Median Absolute Deviation) represents the robust dispersion of the logarithmic resistance difference between target i and target j over a static daily reference period; the MAD represents the median absolute deviation for any variable x. ; when If the dispersion is less than the preset dispersion threshold, the corresponding target pair is retained; otherwise, the corresponding target pair is discarded.

[0031] In S5, the preset thermospheric atmospheric model is an empirical atmospheric model, a physical atmospheric model, or a combination of both. The target pair is obtained based on the preset thermospheric atmospheric model. The midpoint of the background thermal layer density at time t ; Will Mapped to the relative perturbation of thermal layer density: ; In the formula, Indicates target pair The relative perturbation of the thermal layer density at the midpoint of the position at time t; Indicates target pair The differential perturbation mapping coefficients.

[0032] The robust linear regression relationship (robust linear regression model) between the logarithmic density perturbation and the baseline correction difference perturbation was determined through an empirical model. ; In the formula, Indicates target pair The logarithmic density perturbation of the empirical model at time t; For the target The bias term; For the target The residual term at time t represents the difference between the logarithmic density perturbation of the empirical model and the predicted value of the linear regression. Calculate using the following formula: ; Calculate the reference time period for quiet days The median absolute deviation is used as the residual robust dispersion. When the residual robust dispersion is less than a preset residual threshold, the corresponding differential perturbation mapping coefficient is retained. .

[0033] After obtaining the relative perturbation of the thermal layer density, the density enhancement factor, the model-constrained density estimate, and the normalized density perturbation index are calculated based on the relative perturbation of the thermal layer density. The density enhancement factor is calculated as follows: ; In the formula, Indicates target pair Density enhancement factor at time t; The model constraint density estimate is calculated as follows: ; In the formula, Indicates target pair The estimated value of the model constraint density at the midpoint of the model at time t; The normalized density perturbation index is calculated as follows: ; In the formula, Indicates target pair The normalized density perturbation exponent at time t; This indicates a preset lower limit for the dispersion, used to prevent the denominator of the formula from taking the value of zero, thus ensuring computational stability.

[0034] The density enhancement factor is used to intuitively reflect the thermal density perturbation ratio caused by perturbation events. The model-constrained density estimate is used to correct the bias of the empirical model to obtain absolute density data that is closer to real observations. The normalized density perturbation index is used to standardize the perturbation intensity under different scenarios for cross-scenario comparative analysis.

[0035] In S6, for the relative perturbations of thermal layer density of multiple target pairs falling into the same preset spatiotemporal grid g, a robust weighted fusion method is used for spatiotemporal fusion. The fused relative perturbations of thermal layer density are... Represented as: ; In the formula, wmedian represents the weighted median; Indicates target pair The fusion weight at time t.

[0036] Determine according to the following formula: ; Confidence index for calculating spatiotemporal grid g : ; In the formula, This represents the number of target pairs participating in the fusion within the spatiotemporal grid g; Indicates the target's quantitative scale parameter; denoted as the median absolute deviation of the relative perturbation of the thermal layer density within the spatiotemporal grid g; is the average orbital similarity score of the target pairs within the spatiotemporal grid g; It is a pre-defined positive integer used to define the transition benchmark from "insufficient" to "sufficient" target pair quantity, controlling the growth curve of the sample quantity term, when the number of target pairs within the spatiotemporal grid is much greater than... When the sample size term approaches 1, it no longer provides a significant gain to the confidence level.

[0037] Will Normalized to the range of 0 to 1, the result is compared with a preset confidence threshold. If the normalized value is greater than the confidence threshold, the fusion result of the spatiotemporal grid g is retained; otherwise, the fusion result of the spatiotemporal grid is discarded.

[0038] The preset output format may include: spatiotemporal grid center time, spatiotemporal grid center location, height layer, relative perturbation of thermal layer density after fusion, density enhancement factor after fusion, normalized density perturbation index after fusion, model constraint density estimate after fusion, number of target pairs participating in fusion, and confidence index of spatiotemporal grid.

[0039] In this embodiment, various thresholds can be determined by combining the characteristics of the low-Earth orbit target, historical observation data of the static reference period, experimental calibration results, and engineering experience: thresholds for orbital parameters and data quality (orbital height, eccentricity, semi-major axis variation, average motion variation, data time interval, etc.). Thresholds for phase difference, altitude difference, inclination difference, and velocity difference are calibrated based on the normal operating range of low Earth orbit targets, the natural evolution of orbits, and maneuvering characteristics. Statistical thresholds such as logarithmic drag parameter variation threshold, dispersion threshold, residual threshold, and confidence threshold are determined through robust statistical results, including the statistical distribution of massive historical data during undisturbed periods and the median absolute deviation. Orbit similarity scoring-related scale parameters are optimally set based on the results of paired screening experiments. All thresholds can be flexibly fine-tuned according to the actual monitoring area, orbit target type, and application scenario. Other preset parameters can be determined similarly using the above methods.

[0040] Traditional single-target TLE drag analysis methods are easily affected by individual factors such as target attitude, surface-to-mass ratio, drag coefficient, and orbit fitting noise. When estimating thermosphere density directly from single-target drag response parameters, it is easy to incorporate individual target system biases. However, the method in this embodiment uses similar targets in the same orbital shell and orbital plane, combined with logarithmic drag difference component differentiation and static solar baseline correction, which can effectively eliminate or weaken relatively stable individual biases between targets, significantly improve the robustness of thermosphere density disturbance identification, and is more suitable for the full-domain monitoring and disturbance inversion of large-scale low-Earth orbit target groups.

[0041] To verify the effectiveness of the method in this embodiment, a segment of low-Earth orbit target data including the calm period, the disturbance period, and the recovery period was selected for processing. First, the baseline correction differential disturbance of multiple target pairs was processed. Grouped by time and spatial grid, within the same grid, weighted by target-to-orbit similarity score or spatiotemporal grid confidence index, By performing robust weighted statistics, a group resistance disturbance proxy index is constructed, forming a... Figure 2 The time series shown indicates that during the disturbance period, the index significantly increases relative to the level of the quiescent baseline period, suggesting that the method based on the difference of the target-resistance response parameters can effectively capture thermosphere density disturbance events.

[0042] For multiple target pairs falling into the same preset spatiotemporal grid g, the relative perturbation of the thermal layer density is first determined. The density enhancement factor of each target pair was calculated. Then construct the fusion weights A weighted median algorithm is used to robustly fuse the relative perturbations of thermal layer density for all target pairs within the grid, resulting in a fused grid value. Finally, the exponentiation of the fusion result yields the grid median perturbation factor. Finally obtained Time series, such as Figure 3 As shown. During the disturbance period The rise relative to the unperturbed baseline indicates that the method in this embodiment can output model constraint scaling results related to thermal layer density perturbation.

[0043] The median perturbation multiple of the grid was statistically analyzed during different periods of the quiet day, the disturbance period, and the recovery period, and the results are as follows: Figure 4 The box-line distribution is shown. This result reflects the overall distributional differences in the median perturbation multiple of the grid at different time periods, and can be used to assess the density perturbation level during the perturbation period and the recovery period relative to a quiet day.

[0044] At a certain unified epoch, the spatial representation of the median perturbation ratio of each spatiotemporal grid is obtained as follows: Figure 5The figure shows the spatial distribution of the gridded density perturbation. It also shows the median perturbation ratio at the midpoint of the target pair. It can be used to create a thermal density disturbance field to display products.

[0045] Example 2: The difference between this example and Example 1 is that in S1, the orbital data product is precise ephemeris data, and the drag response parameters are obtained based on the target's position and velocity at multiple epochs through orbital dynamics fitting, orbital decay rate estimation, or drag acceleration correction parameter estimation. In S3, the target of the effective target pair is interpolated to a preset unified epoch, that is, the position and velocity at the preset unified epoch are obtained by time interpolation. In subsequent calculations, the estimated drag response parameters are used instead. Participate in the calculation.

[0046] Example 3: A thermal spherical density perturbation inversion system based on orbital drag response parameter difference, comprising the following modules connected in sequence: The data acquisition and preprocessing module is used to acquire orbital data products of low Earth orbit targets, parse them to obtain target identifiers, epochs, orbital state parameters and drag response parameters, and perform quality control to obtain a valid set of low Earth orbit targets. The target pair construction module determines the orbital shell, orbital plane, and in-plane phase angle of effective low-Earth orbit targets based on orbital state parameters. Within the same orbital shell and orbital plane, candidate target pairs are constructed based on the in-plane phase difference and orbital altitude difference. The orbital similarity score is calculated to obtain a set of effective target pairs. The orbit propagation and geometry reconstruction module is used to propagate the target of the effective target pair to a preset unified epoch to obtain orbital geometry information including the target position, altitude, velocity, relative atmospheric velocity, and the position of the target center point; The differential perturbation calculation module calculates the logarithmic drag difference component based on the drag response parameters of the two targets in the target pair, determines the target pair baseline using robust statistics within the static daily reference period, obtains the baseline correction differential perturbation, obtains the robust dispersion using the logarithmic drag difference component, and retains target pairs with robust dispersion less than the preset dispersion threshold. The density perturbation inversion module obtains the background thermosphere density using a preset thermosphere atmospheric model based on the target midpoint position, altitude, and unified epoch. Based on the baseline correction differential perturbation, orbital geometry information, and background thermosphere density, it establishes a mapping relationship between the baseline correction differential perturbation and the relative perturbation of thermosphere density, and inverts to obtain the relative perturbation of thermosphere density at the target midpoint position. The spatiotemporal fusion module is used to perform spatiotemporal fusion of the relative perturbations of thermal density of multiple target pairs, generate thermal density perturbation fields under preset regions and preset time resolutions, and output them in a preset format.

[0047] Example 4: A thermal layer density perturbation inversion device based on orbital drag response parameter difference, comprising: One or more processors; Memory, used to store one or more computer programs; When one or more programs are executed by one or more processors, the one or more processors perform the method in Embodiment 1 or Embodiment 2.

[0048] Example 5: A computer-readable storage medium having executable instructions stored thereon, which, when executed by a processor, cause the processor to perform the method in Example 1 or Example 2.

[0049] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can make equivalent substitutions or modifications based on the technical solution and concept of the present invention within the scope of the technology disclosed in the present invention, and such modifications should also be considered to fall within the scope of protection of the present invention.

Claims

1. A method for inverting thermal layer density perturbation based on the difference in orbital drag response parameters, characterized in that... Includes the following steps: S1. For targets in low Earth orbit, acquire the target's orbital data products, analyze them to obtain the target identifier, epoch, orbital state parameters and drag response parameters, and perform quality control to obtain a valid set of low Earth orbit targets. TLE data targeting track data products, with drag response parameters based on TLE empirical drag parameters. ; The orbital state parameters include orbital elements, position and velocity parameters, and orbital geometry parameters. The orbital elements include the semi-major axis, orbital inclination, eccentricity, right ascension of the ascending node, argument of perigee, mean perigee, and mean motion. The position and velocity parameters include the target's position and velocity at the corresponding epoch. The orbital geometry parameters include orbital altitude, orbital plane markings, and in-plane phase angles. S2. Determine the orbital shell, orbital plane, and in-plane phase angle of the effective low-Earth orbit target based on the orbital state parameters. Within the same orbital shell and the same orbital plane, construct candidate target pairs based on the in-plane phase difference and orbital altitude difference, calculate the orbital similarity score, and obtain the set of effective target pairs. S3. Propagate the target of the effective target alignment to the preset unified epoch to obtain orbital geometric information including target position, altitude, velocity, relative atmospheric velocity and target alignment point position; S4. Calculate the logarithmic drag difference component based on the drag response parameters of the two targets in the target pair, determine the target pair baseline using robust statistics within the static daily reference period, obtain the baseline correction differential disturbance, obtain the robust dispersion using the logarithmic drag difference component, and retain target pairs with robust dispersion less than the preset dispersion threshold. For target pair The logarithmic resistance difference component is calculated according to the following formula: ； In the formula, This represents the logarithmic difference in drag between target i and target j at time t; This indicates that target i at time t ; This indicates that target j at time t ; It is a positive constant; The baseline correction differential perturbation is calculated as follows: ； In the formula, This represents the baseline correction differential perturbation between target i and target j at time t; Indicates the target relative to the baseline, using the static day reference time period. the median; based on To obtain robust dispersion: ； In the formula, MAD represents the robust dispersion of the logarithmic resistance difference component between target i and target j over a static daily reference period; MAD represents the median absolute deviation. Indicates the reference time period for quiet days; when If the dispersion is less than the preset threshold, the corresponding target pair is retained; otherwise, the corresponding target pair is discarded. S5. Based on the target midpoint position, altitude, and unified epoch, obtain the background thermosphere density using a preset thermosphere atmospheric model. Based on the baseline correction differential perturbation, orbital geometry information, and background thermosphere density, establish the mapping relationship between the baseline correction differential perturbation and the relative perturbation of thermosphere density, and invert to obtain the relative perturbation of thermosphere density at the target midpoint position. S6. Perform spatiotemporal fusion of the relative perturbations of thermal layer density for multiple target pairs to generate thermal layer density perturbation fields under preset regions and preset time resolutions, and output them in a preset format.

2. The thermal layer density perturbation inversion method based on the difference in orbital drag response parameters according to claim 1, characterized in that, In S1, quality control includes: removing targets with abnormal trajectories, suspected maneuvering targets, and targets with missing data, as well as removing... Invalid or Less than preset Threshold records; Orbital anomaly targets include targets whose orbital altitude is outside the preset altitude range, whose eccentricity is greater than the preset eccentricity threshold, or whose orbital elements exceed the preset physical reasonable range; Missing data targets include missing target identifiers, epochs, orbital elements, or... The goal is to obtain at least one of the following data; Invalid includes It can be a null value, a non-numeric value, or equal to zero; The criteria for identifying a suspected maneuvering target are as follows: a target is identified as a suspected maneuvering target if at least one of the following conditions is met: The target is that the change in the semi-major axis between two consecutive TLE updates is greater than a preset semi-major axis change threshold; The change in average motion between adjacent epochs is greater than the preset threshold for average motion change. The change in adjacent epochs is greater than the preset value. Threshold for change; At least one of the following parameters—orbital altitude, orbital inclination, eccentricity, right ascension of the ascending node, or mean motion—experiences an anomalous jump between adjacent epochs; The TLE data time interval exceeds the preset time interval threshold.

3. The thermal layer density perturbation inversion method based on the difference in orbital drag response parameters according to claim 1, characterized in that, In S2, constructing candidate target pairs specifically involves: selecting targets with orbital altitudes ranging from 200km to 700km from the effective low-Earth orbit target set; and determining the orbital shell, orbital plane, and in-plane phase angle of each target based on orbital state parameters. Within the same orbital shell and the same orbital plane, candidate target pairs are formed by pairing them up according to preset pairing similarity constraints. The preset pairing similarity constraints include: the phase difference between the two targets in the same target pair within the orbital plane is less than a preset phase difference threshold, and the orbital height difference is less than a preset height difference threshold; And satisfy at least one of the following constraints: the difference in inclination angle between the two target orbits of the same target pair is less than a preset inclination angle difference threshold, and the difference in relative velocity is less than a preset velocity difference threshold.

4. The thermal layer density perturbation inversion method based on the difference in orbital drag response parameters according to claim 1, characterized in that, In S2, for the target pair The orbital similarity score is calculated as follows: ； In the formula, This represents the orbital similarity score between target i and target j; This represents the difference in orbital altitude between target i and target j; This represents the difference in orbital inclination angle between target i and target j; This represents the difference in right ascension between the ascending nodes of target i and target j; This represents the epoch time difference between the orbital data products of target i and target j; , , and They are respectively , , , The corresponding normalized scaling parameter; , , and They are respectively , , , The corresponding weighting coefficients; Calculate orbital similarity scores for each candidate target pair, sort them from highest to lowest orbital similarity scores, and retain the top N candidate target pairs as the set of effective low-Earth orbit target pairs.

5. The thermal layer density perturbation inversion method based on the difference in orbital drag response parameters according to claim 1, characterized in that, In S3, the SGP4 orbital propagation model is used to propagate each target to a preset unified epoch. Before propagation, the drag response parameters of the target are smoothed over time. The time smoothing process uses a fixed-length sliding window with a mean value filter, and the sliding window length is 3 to 24 hours.

6. The thermal layer density perturbation inversion method based on the difference of orbital drag response parameters according to claim 1, characterized in that, In S5, the preset thermospheric atmospheric model is an empirical atmospheric model, a physical atmospheric model, or a combination of both. The target pair is obtained based on the preset thermospheric atmospheric model. The midpoint of the background thermal layer density at time t ; Will Mapped to the relative perturbation of thermal layer density: ； In the formula, Indicates target pair The relative perturbation of the thermal layer density at the midpoint of the position at time t; Indicates target pair The differential perturbation mapping coefficients.

7. The thermal layer density perturbation inversion method based on the difference in orbital drag response parameters according to claim 6, characterized in that, The robust linear regression relationship between the logarithmic density perturbation and the baseline correction difference perturbation, as determined by the empirical model, is as follows: ； In the formula, Indicates target pair The logarithmic density perturbation of the empirical model at time t; For the target The bias term; For the target The residual term at time t represents the difference between the logarithmic density perturbation of the empirical model and the predicted value of the linear regression. Calculate using the following formula: ； In the formula, median represents the median; Calculate the reference time period for quiet days The median absolute deviation is used as the residual robust dispersion. When the residual robust dispersion is less than a preset residual threshold, the corresponding differential perturbation mapping coefficient is retained. .

8. The thermal layer density perturbation inversion method based on the difference of orbital drag response parameters according to claim 1, characterized in that, After obtaining the relative perturbation of the thermal layer density, the density enhancement factor, the model-constrained density estimate, and the normalized density perturbation index are calculated based on the relative perturbation of the thermal layer density. The density enhancement factor is calculated as follows: ； In the formula, Indicates target pair Density enhancement factor at time t; The model constraint density estimate is calculated as follows: ； In the formula, Indicates target pair The estimated value of the model constraint density at the midpoint of the model at time t; The normalized density perturbation index is calculated as follows: ； In the formula, Indicates target pair The normalized density perturbation exponent at time t; This indicates the preset lower limit of dispersion.

9. The thermal layer density perturbation inversion method based on the difference of orbital drag response parameters according to claim 1, characterized in that, In S6, for the relative perturbation of thermal layer density of multiple target pairs falling into the same preset spatiotemporal grid g, a robust weighted fusion method is used for spatiotemporal fusion, and the fused relative perturbation of thermal layer density is... Represented as: ； In the formula, wmedian represents the weighted median; Indicates target pair The fusion weight at time t.

10. The thermal layer density perturbation inversion method based on the difference of orbital drag response parameters according to claim 9, characterized in that, Confidence index for calculating spatiotemporal grid g : ； In the formula, This represents the number of target pairs participating in the fusion within the spatiotemporal grid g; Indicates the target's quantitative scale parameter; denoted as the median absolute deviation of the relative perturbation of the thermal layer density within the spatiotemporal grid g; is the average orbital similarity score of the target pairs within the spatiotemporal grid g; Will Normalized to the range of 0 to 1, the result is compared with a preset confidence threshold. When the normalized value is greater than the confidence threshold, the fusion result of the spatiotemporal grid g is retained.

Images