A method for modeling atmospheric density based on aerodynamic parameter decoupling of mega constellation satellites

Abstract

The application discloses a kind of based on giant constellation satellite aerodynamic parameter decoupling atmospheric density modeling method, by introducing subcomponent dynamic aerodynamic decoupling mechanism, the influence of satellite body attitude and solar wing independent rotation on gas dynamics interaction response is accurately separated and synthesized, the systemic inversion deviation caused by traditional whole satellite constant modeling ignoring configuration mutation is eliminated, the calculation accuracy of total aerodynamic parameter and atmospheric resistance acceleration is significantly improved, to build high-precision, full-coverage global atmospheric density numerical grid model is reversed.In this way, the grid model is calculated with the time-varying configuration of the spacecraft to be measured in a forward mechanical closed loop, making the non-conservative force prediction of the future orbit and the temporal and spatial gradient of the thermosphere environment deeply consistent, effectively suppressing the accumulation of orbit extrapolation integral error, and finally greatly improving the accuracy and timeliness of the orbit warning result, significantly reducing the false alarm rate.

Description

Technical Field

[0001] This invention relates to the field of atmospheric density modeling technology based on decoupling of aerodynamic parameters of giant constellation satellites, and in particular to an atmospheric density modeling method based on decoupling of aerodynamic parameters of giant constellation satellites, an atmospheric density modeling device based on decoupling of aerodynamic parameters of giant constellation satellites, an electronic device, and a readable storage medium. Background Technology

[0002] Thermospheric atmospheric density is a core parameter for the precise determination of low-Earth orbit (LEO) spacecraft orbits, space collision early warning, and on-orbit lifetime estimation. When using LEO satellite-borne navigation data to invert thermospheric atmospheric density, related technologies often neglect the abrupt changes in the windward cross-section caused by the movement of individual components. This leads to a significant disconnect between the calculated aerodynamic interaction response indicators and the actual stress state, and consequently introduces a significant systematic bias in the calculation of the unit mass deceleration effect caused by the drag of the thin atmosphere.

[0003] Furthermore, when using inverted atmospheric density for subsequent space target orbit control and safety assurance, the relevant technologies often rely on solving the spacecraft's orbital evolution separately using independent dynamic equations. This isolated solution mechanism disconnects the evolution of the global dynamic atmospheric grid from the real-time physical configuration changes of the spacecraft under test. When facing the strong time-varying nature and high gradient characteristics of the thin thermosphere environment in the upper atmosphere, the insufficient accuracy of non-conservative mechanics calculations often leads to a rapid accumulation of orbit extrapolation errors in the numerical integration of the dynamic equations. Ultimately, this results in significant spatiotemporal lag or a high false alarm rate in the output spacecraft future orbit warning results. Summary of the Invention

[0004] The present invention provides an atmospheric density modeling method, apparatus, electronic device, and readable storage medium based on the decoupling of aerodynamic parameters of giant constellation satellites, in order to overcome or at least partially solve the above-mentioned problems.

[0005] To solve the above-mentioned technical problems, this application is implemented as follows: In a first aspect, embodiments of this application provide an atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites, including: Receives real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles from satellite downlink; The total aerodynamic parameters characterizing the gas dynamics interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite are determined by the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. Determine the mass of the satellite and its relative velocity with respect to the atmosphere; A global atmospheric density numerical grid model is constructed using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle for the spacecraft under test; The global atmospheric density numerical grid model is invoked to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

[0006] Optionally, the step of determining the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite using the real-time spaceborne global navigation satellite system observation data, the attitude quaternion, and the solar array rotation angle includes: Using observation data from the aforementioned onboard global navigation satellite system, a pre-trained simplified dynamic model is invoked to determine the satellite's velocity sequence; The velocity sequence is smoothed to obtain a smoothed velocity sequence; Determine the velocity vector corresponding to the smoothed velocity sequence; The first angle between the plate normal vector and the velocity vector is calculated using the attitude quaternion and the velocity vector, and the effective windward area of ​​the satellite relative to the direction of incoming flow is calculated based on the first angle. The normal vector of the wing surface is calculated using the attitude quaternion and the rotation angle of the solar wing, and the effective windward area of ​​the movable solar wing relative to the direction of the incoming flow is calculated based on the normal vector of the wing surface. The effective windward area of ​​the main body and the effective windward area of ​​the solar array are respectively combined with the corresponding aerodynamic coefficients and weighted summed to obtain the total aerodynamic parameters that change dynamically with attitude.

[0007] Optionally, the atmospheric drag acceleration is determined as follows: Using the smoothed velocity sequence, the differential acceleration of the velocity at the target time point is calculated; Determine the time series variance of the rate of change of the attitude quaternion; Regularization parameters are adaptively generated based on the time series variance and a preset formula. Construct a numerical differential objective function that includes the regularization parameter, the velocity differential acceleration, and a preset differential operator; The total acceleration vector is calculated by calling the numerical differential objective function and using the smoothed velocity sequence; The higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are determined from the total acceleration vector, and the higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are deleted to obtain the non-conservative force total acceleration vector. The direction determination result is generated by using the angle between the direction of the non-conservative force total acceleration vector and the velocity vector; The acceleration amplitude of the satellite constellation in the adjacent orbital plane is obtained, and the acceleration amplitude is compared with the total non-conservative force acceleration vector to generate a spatial consistency judgment result. Based on the direction determination results and the spatial consistency determination results, the lateral thrust interference, normal thrust interference, and along-track thrust interference in the non-conservative total force acceleration vector are identified. Using a vector decomposition algorithm, the projection result of the total non-conservative force acceleration vector is generated in the opposite direction of the velocity vector. The lateral thrust interference, the normal thrust interference, and the path-direction thrust interference are then removed from the projection result to generate atmospheric drag acceleration in the opposite velocity direction.

[0008] Optionally, the step of constructing a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity includes: The discrete atmospheric density value of the satellite's current position is calculated using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Multiple target satellites are identified, and the discrete atmospheric density values ​​corresponding to the target satellites at the same time are determined as the target discrete atmospheric density values; The target discrete atmospheric density value is mapped to each target grid of a preset global latitude, longitude and altitude grid, and the discrete density set falling into the target grid is statistically analyzed. Calculate the mean and standard deviation of each discrete density set, and based on the mean and standard deviation, remove abnormal density data in the discrete density set whose deviation reaches a preset condition, and select the pure density set corresponding to the target grid. The fusion weight of each target satellite is calculated, and the pure density set is fused using a weighted average algorithm based on the fusion weight. Finally, a global atmospheric density numerical grid model labeled with data coverage and confidence interval is generated and output.

[0009] Optionally, the step of calling the global atmospheric density numerical grid model to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle includes: Acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle for the spacecraft under test; Using observation data from the target global navigation satellite system, the predicted orbital position sequence and predicted velocity sequence of the spacecraft under test within a future preset time period are calculated by extrapolation through a simplified dynamic model. The predicted three-dimensional grid coordinates corresponding to the predicted orbit position sequence are input into the global atmospheric density numerical grid model for spatial interpolation retrieval, and the spatiotemporal background atmospheric density prediction values ​​of each point on the predicted orbit position sequence are matched. The target body normal vector is calculated using the target attitude quaternion and the velocity vector direction of the predicted velocity sequence, and the predicted effective windward area of ​​the spacecraft under test in a future preset time period is calculated based on the target body normal vector. The solar array pointing normal vector is determined by using the target attitude quaternion and the target solar array rotation angle, and the predicted effective windward area of ​​the solar array of the spacecraft under test in the future preset time period is calculated based on the solar array pointing normal vector. The predicted effective windward area of ​​the satellite body and the predicted effective windward area of ​​the solar array are multiplied by the preset basic drag coefficients of the satellite body and the solar array, respectively, and then weighted and summed to calculate the predicted total dynamic aerodynamic parameters of the spacecraft under test in the future preset time period. Using the predicted atmospheric density of the spatiotemporal background, the predicted velocity sequence, and the predicted total dynamic aerodynamic parameters, the predicted atmospheric drag acceleration sequence that the spacecraft under test will experience in a future preset time period is calculated. The predicted atmospheric drag acceleration sequence is used as a non-conservative force input term for numerical integration, and the resulting orbital warning result for the spacecraft under test is output.

[0010] Optionally, the step of smoothing the velocity sequence to obtain a smoothed velocity sequence includes: Obtain the current orbital position of the satellite; The orbital region of the satellite is determined based on its current orbital position, and adjustment conditions for a preset window length are determined based on the orbital region. The adaptively adjusted target filter window length is determined based on the aforementioned adjustment conditions; The velocity sequence is smoothed by using a filter with the target filtering window length to obtain a smoothed velocity sequence.

[0011] Optionally, the adjustment condition is: When the orbital region is determined to be in a polar region, the preset window length is reduced; or, When the orbital region is determined to be in a low-latitude region, the length of the preset window is increased.

[0012] Secondly, embodiments of this application provide an atmospheric density modeling device based on the decoupling of aerodynamic parameters of giant constellation satellites, comprising: The real-time data acquisition module is used to receive real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles from satellite downlink. The total aerodynamic parameters and atmospheric drag acceleration determination module is used to determine the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite, through the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. A satellite feature determination module is used to determine the mass of the satellite and its relative velocity with respect to the atmosphere; The model building module is used to construct a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. The target data acquisition module is used to acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle of the spacecraft under test; The orbital warning result generation module is used to call the global atmospheric density numerical grid model and generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

[0013] Thirdly, embodiments of this application provide an electronic device including a processor, a memory, and a program or instructions stored in the memory and executable on the processor, wherein the program or instructions, when executed by the processor, implement the steps of the method described in the first aspect.

[0014] Fourthly, embodiments of this application provide a readable storage medium on which a program or instructions are stored, which, when executed by a processor, implement the steps of the method described in the first aspect.

[0015] Fifthly, embodiments of this application provide a chip, the chip including a processor and a communication interface, the communication interface being coupled to the processor, the processor being used to run programs or instructions to implement the method as described in the first aspect.

[0016] The embodiments of the present invention have the following advantages: This invention, by receiving multi-source downlink data from satellites in real time and introducing a dynamic aerodynamic decoupling mechanism for components, can accurately separate and synthesize the influence of satellite attitude changes and independent solar array rotation on the gas dynamics interaction response capability. This fundamentally eliminates the systematic inversion bias caused by neglecting the sudden changes in component configurations in traditional whole-satellite constant modeling, significantly improving the calculation accuracy of atmospheric drag acceleration and total aerodynamic parameters. As a result, a more realistic, accurate, and spatiotemporally complete global atmospheric density numerical grid model can be constructed. Based on this, by performing a forward mechanical closed-loop calculation with the time-varying configurations of the spacecraft under test, such as the target attitude and target solar array rotation angle, the non-conservative force prediction of the spacecraft under test in future orbital evolution is deeply consistent with the spatiotemporal gradient characteristics of the global thermosphere environment. This effectively suppresses the rapid accumulation of numerical integration errors during orbit extrapolation, ultimately significantly improving the timeliness and accuracy of spacecraft orbital early warning results and reducing the false alarm rate in terms of safety assurance. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the steps of an atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites, as provided in this embodiment of the invention. Figure 2 This is a flowchart illustrating an atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites, provided in an embodiment of the present invention. Figure 3 This is a structural block diagram of an atmospheric density modeling device based on the decoupling of aerodynamic parameters of a giant constellation satellite provided in an embodiment of the present invention; Figure 4 This is a hardware structure block diagram of an electronic device provided in an embodiment of the present invention; Figure 5 This is a schematic diagram of a computer-readable medium provided in an embodiment of the present invention. Detailed Implementation

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

[0019] The terms "first," "second," etc., used in the specification and claims of this application are used to distinguish similar objects and are not used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of this application can be implemented in orders other than those illustrated or described herein. Furthermore, the character " / " in the specification and claims, indicating at least one of the connected objects, generally indicates that the preceding and following objects are in an "or" relationship.

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

[0021] Reference Figure 1 The diagram illustrates a flowchart of an atmospheric density modeling method based on the decoupling of aerodynamic parameters of a giant constellation satellite, as provided in an embodiment of the present invention. Specifically, it may include the following steps: Step 101: Receive real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angle from the satellite downlink. Step 102: Determine the total aerodynamic parameters characterizing the gas dynamics interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite, using the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. Step 103: Determine the mass of the satellite and its relative velocity with respect to the atmosphere; Step 104: Construct a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity; Step 105: Obtain the target global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle of the spacecraft under test; Step 106: Invoke the global atmospheric density numerical grid model, and generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

[0022] In this embodiment of the invention, real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles can be received from satellite downlink to obtain real-time observation data from low-Earth orbit constellation satellites, providing basic data input for subsequent aerodynamic geometry decoupling of the satellite body and movable solar arrays, and velocity and position inversion.

[0023] Real-time satellite-borne global navigation satellite system observation data refers to the raw positioning observation information, such as carrier phase and pseudorange, received in real time by the receivers on low-Earth orbit satellites during their on-orbit operation.

[0024] Attitude quaternions refer to a combination of four mathematical parameters used to characterize the three-dimensional pointing state of a satellite in its body coordinate system relative to its inertial or orbital coordinate system, thus avoiding Euler angle deadlock.

[0025] Solar wing rotation angle refers to the rotation angle of the independently rotatable solar panel relative to the base mounting plane defined by the satellite body, used to indicate the instantaneous deflection state of the solar wing when it tracks the sun.

[0026] In this embodiment of the invention, the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite, can be determined by the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. By decoupling the dynamic projection state of different satellite components in the thin atmosphere, the overall gas dynamic cross-sectional characteristics of the entire satellite can be accurately quantified, and other irrelevant force interferences can be eliminated to invert and purify the deceleration effect generated by the friction and collision of pure atmospheric particles.

[0027] Total aerodynamic parameters refer to the product of the combined interception cross section and the aerodynamic coefficients of the entire exposed surface structure of a satellite under a specific space attitude and solar array rotation angle, which comprehensively characterize the atmospheric drag effect of the satellite's overall structure along the relative incoming flow direction (usually expressed as C). d A).

[0028] Atmospheric drag acceleration refers to the vector value of the deceleration effect per unit mass in the reverse velocity direction caused purely by atmospheric particle collisions and friction in a spacecraft in a rarefied thermosphere environment, after eliminating all interference from celestial gravity, solar radiation pressure, and orbital maintenance thrust.

[0029] In this embodiment of the invention, the mass of the satellite and its relative velocity relative to the atmosphere can be determined to obtain two state parameters necessary for calculating the thermospheric density inversion formula: the current inertial mass of the satellite and its instantaneous relative velocity with respect to the thermospheric flow. This provides dimensional support for the physical modeling of mass conservation and momentum theorem.

[0030] The mass of a satellite refers to the actual physical mass of the satellite at the current point in time in its orbit, including the total mass of its structure, payload, and remaining propellant.

[0031] The relative velocity of a satellite with respect to the atmosphere refers to the difference between the satellite's velocity vector in inertial space and the velocity vector of the thermosphere due to the Earth's rotation and atmospheric winds.

[0032] In this embodiment of the invention, a global atmospheric density numerical grid model can be constructed using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Based on aerodynamic mechanism formulas, the scattered mechanical response at a single satellite location is transformed into local atmospheric density. Furthermore, through the mapping and fusion of spatiotemporal discrete data from multiple satellites, a standard field of spatial atmospheric density numerical grid with global coverage is inverted and constructed.

[0033] The global atmospheric density numerical grid model refers to a refined numerical simulation model that divides a three-dimensional spatial grid according to preset longitude, latitude, and altitude, and labels the atmospheric density scalar value and spatiotemporal coverage characteristics of each grid node at different times.

[0034] In this embodiment of the invention, the target global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle of the spacecraft under test can be acquired in real time to collect the current instantaneous state data of the target spacecraft that needs to be used for orbit prediction and collision warning, so as to provide measured boundaries for its future configuration time-varying laws and non-conservative force prediction.

[0035] The target global navigation satellite system observation data specifically refers to the raw global navigation satellite system navigation, positioning, and ranging data of the spacecraft currently being observed and warned.

[0036] The target attitude quaternion refers to the three-dimensional attitude quaternion parameter used by the spacecraft under test to indicate its spatial orientation at the current moment.

[0037] The target solar panel deflection angle refers to the deflection angle data of the solar panels recorded at the current moment of the spacecraft under test.

[0038] In this embodiment of the invention, the global atmospheric density numerical grid model can be invoked to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle. The aforementioned refined global atmospheric density grid is used as the environmental background field and deeply integrated with the future dynamic aerodynamic evolution of the spacecraft under test components. This enables accurate forward mechanical calculation of non-conservative atmospheric drag, thereby significantly improving the accuracy of orbital extrapolation integrals and ultimately outputting reliable orbital safety warning information.

[0039] Orbital early warning results refer to the forecast and warning information, including collision probability, dangerous approach distance, or large-scale orbital fall anomalies, output after comparing the spacecraft's future precise position sequence calculated based on orbital extrapolation with dangerous space targets or dangerous collision areas in a spatiotemporal intersection.

[0040] This invention, by receiving multi-source downlink data from satellites in real time and introducing a dynamic aerodynamic decoupling mechanism for components, can accurately separate and synthesize the influence of satellite attitude changes and independent solar array rotation on the gas dynamics interaction response capability. This fundamentally eliminates the systematic inversion bias caused by neglecting the sudden changes in component configurations in traditional whole-satellite constant modeling, significantly improving the calculation accuracy of atmospheric drag acceleration and total aerodynamic parameters. As a result, a more realistic, accurate, and spatiotemporally complete global atmospheric density numerical grid model can be constructed. Based on this, by performing a forward mechanical closed-loop calculation with the time-varying configurations of the spacecraft under test, such as the target attitude and target solar array rotation angle, the non-conservative force prediction of the spacecraft under test in future orbital evolution is deeply consistent with the spatiotemporal gradient characteristics of the global thermosphere environment. This effectively suppresses the rapid accumulation of numerical integration errors during orbit extrapolation, ultimately significantly improving the timeliness and accuracy of spacecraft orbital early warning results and reducing the false alarm rate in terms of safety assurance.

[0041] For example, satellite A in a commercial mega-constellation operating in low Earth orbit is used as the data source satellite, and the space station under test in the same orbit is used as the spacecraft under test for illustration: First, the ground station system receives real-time downlink GNSS observation data from constellation satellite A, its current attitude quaternion (e.g., [0.707, 0, 0, 0.707]), and the current deflection angle of its solar array (e.g., 35 degrees).

[0042] Next, the system uses this data for decoupled calculations: on the one hand, it uses attitude and solar array rotation angle to calculate the projection of the satellite body and battery array in the thin atmosphere, and synthesizes the current accurate total aerodynamic parameter CdA of satellite A; on the other hand, it uses GNSS extrapolation velocity and eliminates forces such as Earth's gravity field, lunar and solar gravitational forces and light pressure to extract the actual atmospheric drag acceleration of satellite A caused by atmospheric particle friction.

[0043] Subsequently, the system obtains the current actual total mass of satellite A (e.g., 200 kg) and the relative atmospheric velocity corrected by the meteorological wind field.

[0044] Then, using aerodynamic formulas and combining the aforementioned acceleration, total aerodynamic parameters, mass, and velocity, the system directly calculates the atmospheric density scalar of satellite A's current three-dimensional spatial location. Through discrete density inversion of hundreds or thousands of satellites in the entire constellation within the same time period, the system maps and interpolates these density data into a three-dimensional grid divided by global latitude, longitude, and altitude, constructing a high-precision global dynamic atmospheric density numerical grid model.

[0045] When the space station under test needs to perform future collision avoidance prediction, the system acquires the target GNSS observation data, target attitude and solar array rotation angle of the space station in real time.

[0046] Finally, the system calls the previously constructed global atmospheric density numerical grid model. Based on the space station's GNSS data, it extrapolates its predicted orbit for the next three days and retrieves the spatiotemporal background atmospheric density values ​​at various points along the orbit from the grid model. Simultaneously, it calculates the space station's future dynamic total aerodynamic parameters by combining the space station's target attitude and solar array rotation angle. The density is multiplied by the space station's own aerodynamic parameters to perform forward mechanics calculations, calculating the precise drag acceleration and performing orbital dynamics integration. Ultimately, it accurately outputs the warning results of the space station's dangerous approach distance and collision trajectory with space debris for the next three days.

[0047] Based on the above embodiments, modified embodiments of the above embodiments are proposed. It should be noted that, in order to keep the description brief, only the differences from the above embodiments are described in the modified embodiments.

[0048] In an optional embodiment of the present invention, the step of determining the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite using the real-time spaceborne global navigation satellite system observation data, the attitude quaternion, and the solar array rotation angle includes: Using observation data from the aforementioned onboard global navigation satellite system, a pre-trained simplified dynamic model is invoked to determine the satellite's velocity sequence; The velocity sequence is smoothed to obtain a smoothed velocity sequence; Determine the velocity vector corresponding to the smoothed velocity sequence; The first angle between the plate normal vector and the velocity vector is calculated using the attitude quaternion and the velocity vector, and the effective windward area of ​​the satellite relative to the direction of incoming flow is calculated based on the first angle. The normal vector of the wing surface is calculated using the attitude quaternion and the rotation angle of the solar wing, and the effective windward area of ​​the movable solar wing relative to the direction of the incoming flow is calculated based on the normal vector of the wing surface. The effective windward area of ​​the main body and the effective windward area of ​​the solar array are respectively combined with the corresponding aerodynamic coefficients and weighted summed to obtain the total aerodynamic parameters that change dynamically with attitude.

[0049] In this embodiment of the invention, by decoupling the dynamic projection state of different geometric components of a satellite in a thin atmosphere, the aerodynamic geometry of the satellite body and the independently rotatable solar array is accurately separated and dynamically synthesized. This eliminates the systematic deviation caused by the neglect of component configuration abrupt changes in traditional whole-satellite constant modeling. Ultimately, it can quantify the whole-satellite gas dynamics interactive response capability that changes dynamically with the attitude and configuration of the spacecraft in real time and with high precision.

[0050] Onboard global navigation satellite system observation data refers to the raw positioning observation information, such as carrier phase and pseudorange, received in real time by the receivers on low-orbit satellites during their on-orbit operation.

[0051] A pre-trained simplified dynamics model refers to a lightweight mechanical evolution mathematical model that is pre-fitted and trained on the ground using historical orbit data. By inputting a small amount of real-time observation information, it can quickly calculate the satellite's current high-precision position and velocity state during the extrapolation stage.

[0052] A velocity sequence is a three-dimensional set of velocity values ​​arranged in chronological order to characterize the speed and direction of a satellite's motion in space at different consecutive moments.

[0053] Smoothing refers to the process of using specific mathematical filtering algorithms to denoise and filter raw data containing high-frequency measurement noise or coarse noise in order to extract a smooth curve that reflects the true low-frequency physical evolution trend.

[0054] The smoothed velocity sequence refers to a set of high signal-to-noise ratio velocity values ​​that, after smoothing and noise filtering, can more realistically and continuously reflect the dynamic characteristics of satellites.

[0055] A velocity vector is a spatial physical quantity that has both magnitude and direction. It is used to precisely indicate the instantaneous direction and speed of a satellite's motion relative to an inertial or orbital coordinate system at a specific moment.

[0056] Attitude quaternions refer to a combination of four mathematical parameters used to characterize the three-dimensional pointing state of the satellite's body coordinate system relative to the inertial or orbital coordinate system, which can avoid Euler angle deadlock problems.

[0057] The flat plate normal vector refers to the unit geometric vector perpendicular to the windward side of the satellite body (such as the main surface of a flat plate configuration satellite), used to characterize the spatial geometric orientation of the main force-bearing surface of the satellite body.

[0058] The first included angle refers to the spatial angle between the normal vector of the satellite's flat plate and the velocity vector of the current satellite motion, which is used to measure the degree to which the satellite is tilted towards the wind.

[0059] The direction of inflow refers to the direction in which atmospheric particles flow into the satellite relative to its high-speed surface. Physically, it is usually completely opposite to the direction of the satellite's velocity vector relative to the atmosphere.

[0060] The effective windward area of ​​a satellite relative to the direction of the incoming flow refers to the geometric cross-sectional area of ​​the satellite's structure on a two-dimensional projection plane perpendicular to the direction of the incoming flow. Its size changes dynamically with the deflection of the satellite's attitude.

[0061] Solar wing rotation angle refers to the rotation angle of the independently rotatable solar panel relative to the base mounting plane defined by the satellite body, used to indicate the instantaneous deflection state of the solar wing when it tracks the sun.

[0062] The wing surface normal vector refers to the unit geometric vector perpendicular to the main light-receiving / wind-facing plane of the movable solar panel, used to characterize the real-time spatial geometric orientation of the solar array surface.

[0063] The effective windward area of ​​a movable solar array relative to the direction of incoming flow refers to the geometric cross-sectional area of ​​the solar panel structure capable of independent deflection on a two-dimensional projection plane perpendicular to the direction of incoming flow. Its size is modulated by both the overall satellite attitude and the rotation angle of the solar array itself.

[0064] The aerodynamic coefficient refers to a dimensionless characteristic constant (such as the drag coefficient C) that characterizes the drag effect generated when the material and geometry of a spacecraft's external surface undergo aerodynamic collisions and friction with particles in a rarefied thermosphere. d ).

[0065] Weighted summation is a mathematical method that involves multiplying each independent physical quantity by its corresponding weighting coefficient (such as its respective basic resistance coefficient) and then summing the results.

[0066] Total aerodynamic parameters refer to the product of the overall interception section and the aerodynamic coefficient of the satellite's entire exposed surface structure under a specific space attitude and solar array rotation angle, which comprehensively characterize the atmospheric drag effect of the satellite's overall structure along the relative incoming flow direction.

[0067] For example, in a specific implementation, the ground data processing system first acquires the raw GNSS observation information collected by the satellite's onboard receiver in real time using downlink telemetry data from the satellite, and directly inputs it into a simplified dynamic model that has been trained on the ground using historical orbit determination data to quickly calculate the satellite's spatial velocity sequence over a continuous time period. To eliminate high-frequency discrete noise introduced by abrupt changes in satellite position or receiver measurements, the system uses a filter with an adaptive window length to smooth the velocity sequence, purifying it into a high signal-to-noise ratio velocity sequence, and extracts the instantaneous three-dimensional velocity vector corresponding to each sampling point. Subsequently, the system extracts the attitude quaternion and solar array rotation angle data for the same time period, and uses the attitude quaternion and velocity vector... By performing spatial geometric projection analysis, the first angle between the normal vector of the satellite body plate and the velocity vector is calculated, and then the effective windward area of ​​the satellite body structure relative to the incoming flow direction is calculated. At the same time, the system combines the attitude quaternion with the recorded angle data of the solar array independently tracking the sun's deflection to deduce the current normal vector of the wing surface, and calculates the effective windward area of ​​the movable solar array under the combined effect of the current attitude and rotation angle based on the normal vector. Finally, the system multiplies and weighted sums the calculated effective windward area of ​​the satellite body and the effective windward area of ​​the solar array with the preset basic drag coefficients of the satellite body and the solar array on the ground, respectively, to synthesize the overall aerodynamic parameters of the entire satellite in real time, which can dynamically evolve with the satellite's flight attitude and component configuration.

[0068] This invention, through the introduction of a component dynamic aerodynamic decoupling mechanism, enables precise separation, calculation, and dynamic weighted synthesis of the effects of changes in the satellite's flight attitude and the independent servo rotation of the solar array on the aerodynamic cross-section. This fundamentally overcomes the systematic inversion bias caused by setting the overall aerodynamic area or drag coefficient of the satellite as a fixed constant and completely ignoring the abrupt changes in the servo configuration of movable components in existing technologies. The calculated total aerodynamic parameters perfectly match the high dynamic characteristics of the attitude-configuration dual time-varying of modern flat-panel giant constellation satellites, significantly improving the quantitative accuracy of the satellite's overall aerodynamic interaction response capability in a thin thermosphere environment. This lays a solid and reliable configuration feature data foundation for subsequent higher-precision, longer-time-series extraction of space atmospheric drag acceleration and dynamic thermosphere atmospheric density inversion.

[0069] In an optional embodiment of the present invention, the atmospheric drag acceleration is determined as follows: Using the smoothed velocity sequence, the differential acceleration of the velocity at the target time point is calculated; Determine the time series variance of the rate of change of the attitude quaternion; Regularization parameters are adaptively generated based on the time series variance and a preset formula. Construct a numerical differential objective function that includes the regularization parameter, the velocity differential acceleration, and a preset differential operator; The total acceleration vector is calculated by calling the numerical differential objective function and using the smoothed velocity sequence; The higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are determined from the total acceleration vector, and the higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are deleted to obtain the non-conservative force total acceleration vector. The direction determination result is generated by using the angle between the direction of the non-conservative force total acceleration vector and the velocity vector; The acceleration amplitude of the satellite constellation in the adjacent orbital plane is obtained, and the acceleration amplitude is compared with the total non-conservative force acceleration vector to generate a spatial consistency judgment result. Based on the direction determination results and the spatial consistency determination results, the lateral thrust interference, normal thrust interference, and along-track thrust interference in the non-conservative total force acceleration vector are identified. Using a vector decomposition algorithm, the projection result of the total non-conservative force acceleration vector is generated in the opposite direction of the velocity vector. The lateral thrust interference, the normal thrust interference, and the path-direction thrust interference are then removed from the projection result to generate atmospheric drag acceleration in the opposite velocity direction.

[0070] This invention accurately extracts and purifies the true drag acceleration generated by atmospheric friction in the rarefied thermosphere from the satellite's motion state. By introducing a numerical differential regularization mechanism based on attitude stability adaptive adjustment, a robust inversion of the total acceleration vector with high precision is achieved. Furthermore, by combining orbital dynamics force model subtraction, multi-criteria thrust interference spatial collaborative identification, and vector decomposition algorithms, lateral, normal, and path-direction thrust interference caused by orbit maintenance or collision avoidance is effectively removed. This fundamentally eliminates the dynamic coupling between thrust interference and atmospheric drag in non-conservative forces, ultimately obtaining an atmospheric drag acceleration vector in the reverse velocity direction with extremely high purity.

[0071] Velocity differential acceleration refers to the observed acceleration value obtained by performing primary time differential calculations on a smoothed velocity sequence, which includes high-frequency differential amplification noise and is used to preliminarily characterize the rate of change of satellite velocity over time.

[0072] The rate of change of attitude quaternions refers to how quickly the satellite attitude quaternions change over time within a continuous time window, and is used to quantify the high-dynamic fluctuations or maneuvers of satellite attitude in space.

[0073] Time series variance refers to the dispersion of all sample values ​​of the rate of change of attitude quaternions relative to their statistical mean within a set sliding time window. Physically, it is used to quantitatively assess the stability of a satellite's current attitude.

[0074] The preset formula refers to the algorithmic calculation formula that is pre-set to establish a positive correlation between regularization parameters and attitude stability indicators (such as time series variance).

[0075] The regularization parameter refers to the coefficient used in the numerical differential objective function to coordinate the weights between the data fitting residual term and the smoothing constraint term of the solution (such as the Tikhonov regularization parameter λ). Its magnitude determines the degree of penalty for high-frequency fluctuations in acceleration.

[0076] Preset differential operators refer to mathematical difference matrices (such as first-order or second-order difference operators L) that are pre-defined in the algorithm and are used to constrain and calculate the roughness or fluctuation of the total acceleration sequence to be determined in the regularization term.

[0077] The numerical differential objective function refers to a mathematically optimized functional that includes a penalty for velocity observation fitting error and an acceleration smoothing constraint, used to suppress the malignant amplification of high-frequency noise during the numerical differential process.

[0078] The total acceleration vector refers to the resultant acceleration vector caused by the combined action of all conservative forces (various types of gravity) and non-conservative forces (drag, thrust, light pressure, etc.) when a satellite moves in space.

[0079] Higher-order gravitational acceleration refers to the higher-order harmonic acceleration of Earth's gravity caused by the Earth's non-spherical shape and non-uniform density distribution (such as the gravitational acceleration calculated based on the EGM2008 gravity field model).

[0080] Three-body gravitational acceleration refers to the point mass perturbation gravitational acceleration exerted on low-orbit satellites by other large celestial bodies in the solar system, such as the Sun and the Moon.

[0081] Solar radiation pressure acceleration refers to the weak mechanical acceleration caused by the collision of solar radiation photons with the exposed surface of a satellite and the resulting momentum exchange.

[0082] The non-conservative force total acceleration vector refers to the resultant acceleration vector composed purely of surface forces or internal forces such as atmospheric drag and satellite thruster thrust, after deducting all conservative gravitational perturbations such as higher-order gravitational fields and three-body gravity, as well as solar radiation pressure acceleration, from the total acceleration of the satellite.

[0083] The direction determination result refers to the geometric direction determination conclusion generated by using the three-dimensional spatial angle relationship between the total acceleration vector of non-conservative forces and the instantaneous velocity vector of the satellite to determine whether the current force orientation deviates from the pure aerodynamic drag direction.

[0084] A neighboring orbital plane satellite constellation refers to a collection of other background reference satellites in the same giant constellation that operate in the same or adjacent orbital planes and are within a predetermined spatiotemporal proximity.

[0085] Acceleration amplitude refers to the scalar magnitude (i.e., the magnitude of the vector) of the non-conservative force acceleration experienced by a constellation of satellites in a nearby orbital plane at the same or similar time, and is used as a benchmark for the local space environment.

[0086] The spatial consistency judgment result refers to the comparison between the non-conservative total acceleration vector of the current position of a single satellite and the overall acceleration amplitude characteristics of the satellite constellation in the adjacent orbital plane, to determine whether it conforms to the global spatial gradual change law of the density of the thin atmosphere, thereby identifying local anomalous sudden changes (such as autonomous maneuvering to reduce torque) and generating a collaborative discrimination conclusion.

[0087] Lateral thrust interference refers to the deceleration or acceleration interference caused by the lateral thrust generated by the ion thruster perpendicular to the orbital plane when a satellite performs orbital offset, orbital plane adjustment, or attitude control.

[0088] Normal thrust interference refers to the geocentric / backward thrust interference generated by the thrusters when a satellite makes minor adjustments to its orbital altitude or phase. This thrust is perpendicular to the velocity direction and lies within the orbital plane.

[0089] Thrust interference along the trajectory refers to the thrust generated by the thrusters along the direction (or opposite direction) of the satellite's velocity when the satellite frequently performs orbit maintenance (such as overcoming drag to raise the orbit).

[0090] Vector decomposition algorithm refers to a spatial geometric analytical calculation method that projects a complex three-dimensional acceleration vector into multiple independent components according to the selected orthogonal coordinate system basis vectors (such as along the running velocity and its normal and lateral directions).

[0091] The projection result refers to the scalar or vector component of the total non-conservative force acceleration vector, which is precisely mapped onto the opposite direction of the satellite's instantaneous velocity vector (i.e., the direction of theoretical pure atmospheric drag) through vector decomposition.

[0092] Atmospheric drag acceleration in the opposite direction of velocity refers to the actual atmospheric drag acceleration generated purely by the friction of rarefied atmospheric particles that hinders satellite motion, after completely stripping away and eliminating interference from active thrust along the track direction and other directional forces from the projection results in the opposite direction of velocity.

[0093] For example, in the specific implementation, the ground computing system first uses the smoothed velocity sequence output from the previous stage to perform time difference calculation to obtain the velocity differential acceleration value at the target time point. Simultaneously, the system extracts the attitude quaternions for the same period and calculates their rate of change. It then uses a sliding window to statistically analyze the time series variance to quantitatively assess the satellite's attitude stability in real time. Next, the system substitutes this variance into a preset positive correlation mapping formula to adaptively generate a regularization parameter (i.e., automatically reducing the regularization constraint to capture real mechanical abrupt changes during frequent attitude maneuvers and increasing the constraint to filter noise during stable attitudes). This parameter is then combined with the velocity differential acceleration and a preset first- or second-order difference matrix to construct a Tikhonov regularized numerical differential objective function. The system then calls an optimization solver to calculate the high signal-to-noise ratio total acceleration vector using the smoothed velocity. Subsequently, the system accurately calculates the current time by calling a standard high-order gravity field model (such as EGM2008), the lunar and solar positions calculated from ephemeris data, and the radiation pressure model. The system first extracts the higher-order gravitational field acceleration, three-body gravitational acceleration, and solar radiation pressure acceleration from the previous position, and then directly subtracts these from the total acceleration vector to purify the non-conservative force total acceleration vector. Next, the system initiates a dual collaborative interference suppression process: on the one hand, it uses the non-conservative force total acceleration and velocity vector to calculate the spatial angle to generate a direction judgment result; on the other hand, it acquires the acceleration amplitude of satellites in adjacent orbital planes in the same constellation in real time for spatiotemporal comparison to generate a spatial consistency judgment result. Then, the system integrates these two judgment results to intelligently identify the lateral, normal, and track-direction thrust interference in the non-conservative force total acceleration caused by frequent orbit maintenance or collision avoidance by commercial satellites. Finally, the system uses a vector decomposition algorithm to spatially project the non-conservative force total acceleration vector in the opposite direction of the satellite's instantaneous velocity, and completely removes and strips the identified thrust interference components from the projected components, accurately generating and outputting a high-purity, thrust-free atmospheric drag acceleration in the opposite velocity direction.

[0094] This invention, through the construction of a joint framework of adaptive numerical differentiation for attitude stability and multi-criteria spatial collaborative thrust suppression, effectively solves the industry problem of existing technologies being unable to cope with high-frequency and complex dynamic interference of giant commercial satellite constellations. The beneficial effects of this invention are twofold: First, by introducing a regularization parameter that adaptively adjusts with attitude variance, it can prevent acceleration waveform distortion during large-angle bias maneuvers and strongly suppress amplified noise from numerical differentiation during stable operation, significantly improving the inversion quality of total acceleration. Second, by introducing dual criteria of direction and spatial consistency, it breaks through the technical blind spot of traditional technologies that rely solely on amplitude to distinguish between local atmospheric density lateral gradient changes and active thrust interference. It can accurately identify and completely eliminate lateral, normal, and extremely difficult-to-identify along-track thrust interference introduced by frequent orbit maintenance and collision avoidance, fundamentally eliminating the severe contamination of atmospheric drag by thrust interference in non-conservative forces, ensuring the extremely high purity and mechanical authenticity of the final extracted atmospheric drag acceleration in the reverse velocity direction.

[0095] In an optional embodiment of the present invention, the step of constructing a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity includes: The discrete atmospheric density value of the satellite's current position is calculated using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Multiple target satellites are identified, and the discrete atmospheric density values ​​corresponding to the target satellites at the same time are determined as the target discrete atmospheric density values; The target discrete atmospheric density value is mapped to each target grid of a preset global latitude, longitude and altitude grid, and the discrete density set falling into the target grid is statistically analyzed. Calculate the mean and standard deviation of each discrete density set, and based on the mean and standard deviation, remove abnormal density data in the discrete density set whose deviation reaches a preset condition, and select the pure density set corresponding to the target grid. The fusion weight of each target satellite is calculated, and the pure density set is fused using a weighted average algorithm based on the fusion weight. Finally, a global atmospheric density numerical grid model labeled with data coverage and confidence interval is generated and output.

[0096] This invention transforms multi-source, fragmented, and noise-infused discrete atmospheric density observations from single satellites into a refined and highly reliable global standardized space environment model. By introducing a spatiotemporal grid mapping mechanism, a statistically based anomaly density removal strategy, and a weighted averaging algorithm based on multi-satellite fusion weights, it achieves the refinement and collaborative fusion of heterogeneous observation data within the space grid. This results in the construction and output of a global atmospheric density numerical grid model that accurately reflects the dynamic evolution of the thermospheric atmosphere and is labeled with complete quality evaluation indicators (data coverage and confidence intervals).

[0097] Atmospheric drag acceleration refers to the vector value of the deceleration effect per unit mass in the reverse velocity direction caused purely by atmospheric particle collisions and friction in a spacecraft in a rarefied thermosphere environment, after eliminating all interference from celestial gravity, solar radiation pressure, and orbital maintenance thrust.

[0098] Total aerodynamic parameters refer to the product of the overall interception section and the aerodynamic coefficient of the satellite's entire exposed surface structure under a specific space attitude and solar array rotation angle, which comprehensively characterize the atmospheric drag effect of the satellite's overall structure along the relative incoming flow direction.

[0099] Satellite mass refers to the actual physical mass of the satellite at the current point in time in its orbit, including the total mass of its main structure, payload, and remaining propellant.

[0100] Relative velocity refers to the difference between the velocity vector of a satellite in inertial space and the velocity vector of the thermosphere due to the Earth's rotation and atmospheric winds.

[0101] Discrete atmospheric density values ​​refer to the scalar values ​​of local thermosphere atmospheric density calculated independently using aerodynamic inversion formulas based on the mechanical response parameters of a single satellite at a specific instantaneous location and time point.

[0102] Multiple target satellites refer to a collection of satellites in different orbits that operate in a giant constellation, are within the same observation window, and can provide effective atmospheric density inversion data.

[0103] The target discrete atmospheric density value refers specifically to the combination of discrete atmospheric density observation values ​​produced by multiple selected target satellites at the same specific physical moment.

[0104] The global latitude, longitude, and altitude grid refers to a set of three-dimensional discrete geometric grids that divide the space of the Earth's thermosphere according to preset longitude resolution, latitude resolution, and altitude layer thickness.

[0105] The target grid refers to the specific grid cell in the three-dimensional spatial grid into which at least one target discrete atmospheric density value falls within its spatial boundary at the current time node.

[0106] Discrete density set refers to the mathematical set of all discrete atmospheric density values ​​of targets that fall exactly within a specific target grid at all physical locations within the same time window.

[0107] The mean refers to the arithmetic mean of all density scalar data in a discrete density set within a target grid, used to characterize the central tendency of atmospheric density within that grid.

[0108] The standard deviation refers to the average square root of the deviation of each density data in the discrete density set from its arithmetic mean. Physically, it is used to quantify the dispersion and noise level of multi-source observation data within the grid.

[0109] Abnormal density data refers to gross error density values ​​in discrete density sets that deviate from the normal statistical distribution range and are severely distorted due to residual effects of single-star non-conservative force removal, sensor malfunctions, or maneuvering disturbances.

[0110] Preset conditions refer to pre-defined statistical criteria used to identify and eliminate outliers (such as the three standard deviation boundary in the Laida criterion, or the Lameda criterion).

[0111] A pure density set refers to a set of high signal-to-noise ratio observation data that remains after all anomalous density data have been completely removed from a discrete density set according to preset conditions. This set can truly reflect the atmospheric density characteristics of the network pattern.

[0112] Fusion weight refers to the relative weight coefficients that are pre-set or dynamically calculated based on the basic orbit determination accuracy, aerodynamic decoupling accuracy, or sensor aging status of each target satellite, and are used to measure the reliability of the satellite data.

[0113] The weighted average algorithm refers to a mathematical calculation method that multiplies each density data in the pure density set by the fusion weight of its respective satellite, then sums them up and divides by the sum of the total weight coefficients to obtain a comprehensive average value that takes into account the credibility of all parties.

[0114] Data coverage refers to the proportion of the total number of grids in a global atmospheric density numerical grid model where actual, valid, and clean data has been incorporated and fused within a specific time period. It is used to characterize the spatial completeness of the model.

[0115] The confidence interval refers to the numerical range that includes the true atmospheric density value at a preset probability level, calculated based on the statistical distribution characteristics of multi-source observations within the grid. It is used to characterize the statistical reliability and uncertainty of the grid density inversion result.

[0116] The global atmospheric density numerical grid model refers to a refined numerical simulation model that divides a three-dimensional spatial grid according to preset longitude, latitude, and altitude, and labels the atmospheric density scalar value and spatiotemporal coverage characteristics of each grid node at different times.

[0117] For example, in the specific implementation, the ground-based processing system first uses the atmospheric drag acceleration in the inverse velocity direction, the total dynamic aerodynamic parameters of the entire satellite, the current mass of the satellite, and the relative atmospheric velocity calculated by the previous stage, and substitutes them into the hydrodynamic inversion equation to calculate the discrete atmospheric density value of a single satellite at its current orbital position; then, the system selects multiple target satellites in the entire mega-constellation that are at the same sampling moment, and aggregates their discrete atmospheric density values ​​at this time into a unified target discrete atmospheric density value sequence; next, the system establishes a preset global latitude, longitude, and altitude three-dimensional geographic information grid, and accurately maps the target discrete atmospheric density values ​​to the corresponding target grids according to the real-time three-dimensional orbital coordinates of each satellite, and counts and classifies the discrete density sets falling within each target grid; then... The system calculates the arithmetic mean and standard deviation of the discrete density set within each grid, and uses a preset three-standard-deviation criterion to automatically identify and remove gross errors and anomalous density data with excessive deviations, thus filtering and purifying the pure density set corresponding to each target grid. Then, based on the geometric accuracy of the orbit determination of each satellite and the confidence level of attitude decoupling, the system calculates the fusion weight of each target satellite in its respective grid in real time. Finally, based on this fusion weight, the system uses a weighted average algorithm to perform multi-source assimilation and fusion processing on the pure density set within each grid, calculates the final representative density value of each grid point, and statistically analyzes the data fullness and variance characteristics of the entire grid, generating and outputting a refined global atmospheric density numerical grid model labeled with data coverage and confidence intervals.

[0118] This invention, through its pioneering global spatiotemporal grid mapping-multi-satellite statistical denoising-weighted fusion assimilation architecture, effectively solves the problems of data sparsity, high noise, and lack of quality quantification indicators caused by the reliance on single-satellite inversion in existing technologies. The beneficial effects of this invention are as follows: First, leveraging the advantage of simultaneous observation of multiple satellites in a giant constellation, it achieves a leap from single-satellite linear observation to global planar grid monitoring through spatial grid mapping, significantly improving the spatiotemporal coverage of the space atmosphere; second, it introduces an anomaly data removal mechanism based on mean and standard deviation, which adaptively eliminates gross errors introduced by residual maneuvering control of individual commercial satellites and occasional sensor malfunctions, ensuring the purity of the grid background data; third, it employs a weighted average algorithm based on fusion weights, fully considering the differences in the performance of individual satellites from different sources, making the inverted grid density field more physically realistic; fourth, the final model innovatively labels the data coverage and confidence intervals, providing crucial quantitative basis for the mechanical input uncertainty of subsequent orbital warnings and numerical extrapolation integration.

[0119] In an optional embodiment of the present invention, the step of calling the global atmospheric density numerical grid model to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle includes: Acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle for the spacecraft under test; Using observation data from the target global navigation satellite system, the predicted orbital position sequence and predicted velocity sequence of the spacecraft under test within a future preset time period are calculated by extrapolation through a simplified dynamic model. The predicted three-dimensional grid coordinates corresponding to the predicted orbit position sequence are input into the global atmospheric density numerical grid model for spatial interpolation retrieval, and the spatiotemporal background atmospheric density prediction values ​​of each point on the predicted orbit position sequence are matched. The target body normal vector is calculated using the target attitude quaternion and the velocity vector direction of the predicted velocity sequence, and the predicted effective windward area of ​​the spacecraft under test in a future preset time period is calculated based on the target body normal vector. The solar array pointing normal vector is determined by using the target attitude quaternion and the target solar array rotation angle, and the predicted effective windward area of ​​the solar array of the spacecraft under test in the future preset time period is calculated based on the solar array pointing normal vector. The predicted effective windward area of ​​the satellite body and the predicted effective windward area of ​​the solar array are multiplied by the preset basic drag coefficients of the satellite body and the solar array, respectively, and then weighted and summed to calculate the predicted total dynamic aerodynamic parameters of the spacecraft under test in the future preset time period. Using the predicted atmospheric density of the spatiotemporal background, the predicted velocity sequence, and the predicted total dynamic aerodynamic parameters, the predicted atmospheric drag acceleration sequence that the spacecraft under test will experience in a future preset time period is calculated. The predicted atmospheric drag acceleration sequence is used as a non-conservative force input term for numerical integration, and the resulting orbital warning result for the spacecraft under test is output.

[0120] This invention utilizes a refined global atmospheric density numerical grid model constructed through a previous inversion as a high-precision environmental background field. This model is then deeply integrated with the time-varying configuration and three-dimensional motion state of the spacecraft under test in the future, and a closed-loop mechanical calculation is performed. By accurately predicting the future dynamic effective windward area of ​​the spacecraft's components and the matching spatiotemporal background atmospheric density, a precise forward physical extrapolation of the non-conservative atmospheric drag acceleration within a preset future time period is achieved. This effectively suppresses the accumulation of errors during numerical integration orbit extrapolation, ultimately outputting a highly reliable spacecraft space orbit safety early warning result with a low false alarm rate.

[0121] The spacecraft under test refers to a specific target satellite or space station that is currently being observed and predicted by ground-based space surveillance systems, or that urgently needs to undergo space collision avoidance and orbital drop assessment.

[0122] The target global navigation satellite system observation data specifically refers to the raw global navigation satellite system navigation, positioning, and ranging data of the spacecraft being observed and warned of, which is currently descending at its current orbital point.

[0123] The target attitude quaternion refers to the combination of four mathematical parameters used by the spacecraft under test at the current moment to indicate the three-dimensional orientation of its body coordinate system relative to the inertial coordinate system.

[0124] The target solar array rotation angle refers to the instantaneous rotation angle of the independently deflectable solar array panel of the spacecraft under test relative to its defined mounting surface, as recorded at the current moment.

[0125] A simplified dynamics model refers to a pre-established lightweight mechanical evolution mathematical model that, after inputting initial navigation observation data, can quickly extrapolate the future orbital state of a spacecraft during the extrapolation phase using analytical or semi-analytical algorithms.

[0126] Predicted orbital position sequence refers to the set of three-dimensional spatial position coordinates of the spacecraft under test at different consecutive times within a preset time period, calculated by extrapolating and integrating the dynamic model over time.

[0127] Predicted velocity sequence refers to the set of three-dimensional velocity values ​​of a spacecraft under test at different consecutive moments within a predetermined time period, calculated by extrapolation, reflecting its future speed and orientation.

[0128] Predicting 3D grid coordinates refers to converting and mapping the predicted orbital position of the spacecraft under test into the standard geographical longitude, latitude, and spatial altitude 3D spatial grid node coordinates used in the preceding global atmospheric density numerical grid model.

[0129] The global atmospheric density numerical grid model refers to a refined numerical simulation model that divides a three-dimensional spatial grid according to preset longitude, latitude, and altitude, and labels the atmospheric density scalar value and spatiotemporal coverage characteristics of each grid node at different times.

[0130] Spatial interpolation retrieval refers to the process of performing multidimensional mathematical interpolation between adjacent discrete grid nodes of a global atmospheric density numerical grid model based on the predicted three-dimensional grid coordinates of the spacecraft, thereby matching and extracting the continuous spatiotemporal background atmospheric density at that location.

[0131] The predicted atmospheric density of the spatiotemporal background refers to the estimated scalar value of the local thermosphere atmospheric density at a specific time and location point of the spacecraft under test in the future predicted orbit, which is retrieved and output by the spatial interpolation of the grid model.

[0132] The velocity vector direction of the predicted velocity sequence refers to the geometric direction of the three-dimensional velocity vector of the spacecraft under test in inertial space or orbital coordinate system at a future moment, representing the instantaneous convection direction during the spacecraft's future flight.

[0133] The target body normal vector refers to the unit vector perpendicular to the geometric principal surface of the spacecraft under test, representing the instantaneous orientation of the spacecraft's structure in space.

[0134] The predicted effective windward area of ​​the spacecraft body refers to the extrapolated predicted two-dimensional projected plane geometric cross-sectional area of ​​the spacecraft body structure at a specific future moment, perpendicular to the direction of the relative incoming flow (opposite to the velocity vector).

[0135] The solar array pointing normal vector is a unit geometric vector perpendicular to the main wind- and light-receiving surface of the movable solar array of the spacecraft under test, used to characterize the spatial geometric orientation of the solar array at future moments.

[0136] The predicted effective windward area of ​​the solar array refers to the two-dimensional projected plane geometric cross-sectional area of ​​the solar panel structure perpendicular to the relative incoming flow direction at a specific future moment, after comprehensively considering the future attitude of the spacecraft under test and the angle of independent follow-up deflection of the solar array.

[0137] The basic drag coefficient of a satellite body refers to the dimensionless aerodynamic drag characteristic constant that is pre-set based on the physical properties of the external materials and fixed geometry of the spacecraft under test.

[0138] The basic drag coefficient of a solar panel refers to a dimensionless aerodynamic drag characteristic constant that is pre-set based on the scattering characteristics of the solar panel material surface and its flat shape.

[0139] Weighted summation is a mathematical calculation method that involves multiplying the predicted effective windward area of ​​the solar panel and the predicted effective windward area of ​​the solar array by their respective basic drag coefficients, and then summing the products.

[0140] Predicted dynamic total aerodynamic parameters refer to the estimated combination of the comprehensive interception section and the aerodynamic coefficients of the entire exposed surface shape of the spacecraft under test at a specific time and in a specific attitude configuration in the future, which comprehensively characterizes the atmospheric drag effect of the spacecraft along the relative incoming flow direction. (This is often expressed as predicted CdA).

[0141] Predicted atmospheric drag acceleration sequence refers to the set of atmospheric drag deceleration vectors that the spacecraft under test will experience at each moment in a predetermined time period in the future, calculated using the forward aerodynamic formula, combined with the predicted background atmospheric density, predicted flight speed, and predicted dynamic total aerodynamic parameters.

[0142] Non-conservative force inputs refer to the perturbational external force terms in the spacecraft orbital dynamics equations, which are composed of external surface forces such as atmospheric drag and solar radiation pressure, in addition to various conservative celestial body gravitational forces.

[0143] Numerical integration refers to the numerical calculation process of using a computer to perform time-step accumulation and iteration on the total differential dynamics equations of a spacecraft that include all conservative and non-conservative force inputs, thereby solving its time-series state.

[0144] The orbital early warning result refers to the prediction and warning information, including collision probability, safe approach distance, or orbital fall anomaly, output after comparing the spacecraft's future precise position sequence obtained by numerical integration orbital extrapolation with dangerous space targets or dangerous areas in a spatiotemporal manner.

[0145] For example, in a specific implementation, the space safety monitoring system first receives the target GNSS observation data, target attitude quaternions, and target solar array rotation information of the spacecraft under test in real time via the space-to-ground telemetry and control link. Then, the system inputs the target GNSS observation data as initial boundary conditions into a simplified dynamic model, performs orbital extrapolation to calculate the predicted orbital position sequence and predicted velocity sequence of the spacecraft under test within a preset future time period. Next, the system converts the predicted orbital position sequence into standard three-dimensional grid coordinates and directly inputs it into the global atmospheric density numerical grid model inverted from the previous stage. Multidimensional spatial interpolation is performed between spatiotemporally adjacent grid nodes to accurately match the predicted spatiotemporal background atmospheric density values ​​at various points along the spacecraft's future orbit. During this period, the system uses the downlink target attitude quaternions in conjunction with the velocity vector direction of the predicted velocity sequence to calculate the future target body normal vector, and thereby resolves the predicted effective windward area of ​​the spacecraft under test within the future time period. Simultaneously, the system combines the target attitude quaternions... The system determines the future solar wing pointing normal vector based on the current solar wing rotation angle, and then analyzes the predicted effective windward area of ​​the movable solar wing in the future. Next, the system multiplies the predicted effective windward area of ​​the main body and the predicted effective windward area of ​​the solar wing with preset basic drag coefficients for the main body and solar wing, respectively, and performs a weighted summation to synthesize online the predicted dynamic total aerodynamic parameters that reflect the highly dynamic changes of the spacecraft's configuration. Subsequently, the system calls the aerodynamic mechanism equations, multiplies the previously retrieved predicted values ​​of the spatiotemporal background atmospheric density, the extrapolated predicted velocity sequence, and the synthesized predicted dynamic total aerodynamic parameters to accurately deduce the predicted atmospheric drag acceleration sequence that the spacecraft under test will experience in the future. Finally, the system uses this predicted atmospheric drag acceleration sequence as a key non-conservative force input, assembles it together with a standard conservative gravitational field model into a high-order orbit integrator for numerical integration, completes the final refined orbit extrapolation prediction, and performs a rendezvous and comparison with the ephemeris of space debris to output a collision or descent trajectory warning result for the spacecraft under test.

[0146] This invention, through its pioneering architecture that integrates environmental background grids with the future dynamic configuration of spacecraft for forward physical closed-loop mechanical calculations, effectively solves the industry pain points of existing technologies in orbit extrapolation, such as the fragmentation of non-conservative force physical elements and the exponential vicious accumulation of extrapolation errors over time. The beneficial effects of this invention are as follows: First, it breaks through the blind spot of relying on static constant configurations or semi-empirical environmental models in traditional extrapolation. It introduces a high-precision measured inversion grid field through spatial interpolation retrieval and innovatively introduces the prediction and weighted dynamic synthesis of the effective windward area of ​​the future spacecraft body and solar array components. This makes the calculated predicted dynamic total aerodynamic parameters and predicted atmospheric drag acceleration sequences have extremely high mechanical realism and time-varying tracking capability. Second, it enables the spacecraft under test to accurately sense the spatiotemporal gradient characteristics of the global thermosphere environment and the mechanical fluctuations caused by its own configuration changes at every moment on the future extrapolated orbit. This effectively suppresses the non-conservative force accumulation error in the long-term orbit numerical integration extrapolation process, and ultimately significantly improves the timeliness and three-dimensional spatiotemporal accuracy of the orbital early warning results for low-Earth orbit spacecraft, greatly reducing the false alarm rate and missed alarm rate at the space security level.

[0147] In an optional embodiment of the present invention, the step of smoothing the velocity sequence to obtain a smoothed velocity sequence includes: Obtain the current orbital position of the satellite; The orbital region of the satellite is determined based on its current orbital position, and adjustment conditions for a preset window length are determined based on the orbital region. The adaptively adjusted target filter window length is determined based on the aforementioned adjustment conditions; The velocity sequence is smoothed using a filter with the target filtering window length to obtain a smoothed velocity sequence. The adjustment conditions are as follows: When the orbital region is determined to be in a polar region, the preset window length is reduced; or, When the orbital region is determined to be in a low-latitude region, the length of the preset window is increased.

[0148] This invention, in its embodiment, introduces a dynamic control mechanism that adaptively adjusts the filtering intensity based on the space orbit region when performing denoising and smoothing on satellite velocity sequences. Because the density of the Earth's thermosphere in polar regions (high latitudes) experiences drastic and high-frequency local spatial abrupt changes (high gradient characteristics) due to the influence of solar storms, geomagnetic storms, and auroral particle deposition, while its evolution is relatively gradual in low latitudes, this embodiment dynamically shortens the polar region filtering window to retain the true details of high-frequency atmospheric dynamic fluctuations by identifying the satellite's spatial geographic region in real time. Simultaneously, it lengthens the low-latitude region filtering window to powerfully filter out measurement noise. This ensures the quality of global velocity smoothing while fundamentally preventing the malicious smoothing and masking of key spatial dynamic features in polar regions by traditional fixed-window filters.

[0149] The satellite's current orbital position refers to the satellite's three-dimensional spatial geometric coordinates relative to the Earth's central coordinate system (such as WGS-84) or inertial coordinate system at the current specific sampling time.

[0150] The orbital region refers to the regional division of the thermosphere above the Earth based on geographical latitude and longitude or spatial altitude. In this embodiment, it specifically refers to the polar regions and low-latitude regions divided according to latitude.

[0151] The preset window length refers to the initial number of regular time steps or velocity sampling points that are pre-set in the filter algorithm when performing sliding smoothing calculations.

[0152] Adjustment conditions refer to the logical criteria pre-written in the control program, which are used to adaptively trigger mathematical control rules for stretching or shortening the filter window based on the latitude characteristics of the satellite's current orbital region.

[0153] Adaptive adjustment refers to the process by which an algorithm can automatically calculate and change the internal control parameters of a filter in real time based on the input satellite orbit position without human intervention.

[0154] The target filtering window length refers to the actual number of sampling points that are ultimately generated after the preset window length undergoes adaptive calculation based on adjustment conditions and is directly delivered to the filter to perform current point smoothing calculation.

[0155] A filter is a mathematical algorithm structure used to process discrete time series data (such as a moving average filter, a Gaussian filter, or a Kalman filter), which suppresses high-frequency noise by performing statistical calculations on the data within a window.

[0156] A velocity sequence is a three-dimensional set of velocity values ​​arranged in chronological order to characterize the speed and direction of a satellite's motion in space at different consecutive moments.

[0157] Smoothing filtering refers to the mathematical calculation process of using filters to denoise raw data containing high-frequency measurement noise or abrupt coarse noise in order to extract a smooth curve that reflects the true low-frequency physical evolution trend.

[0158] The smoothed velocity sequence refers to the set of high signal-to-noise ratio velocity values ​​that have been processed by a filter with a target filtering window length after the original velocity sequence has been denoised, thus eliminating measurement noise while preserving local characteristics.

[0159] The polar regions refer to the space orbits above the Earth's high latitude poles (usually referring to the regions between 60 and 90 degrees of geomagnetic or geographical latitude), which are most intensely affected by solar and geomagnetic activity.

[0160] Low latitude regions refer to the space orbital regions above the Earth's equator and its sides at mid-to-low latitudes (usually referring to the geographical latitude between 0 and 30 degrees). The thermosphere in this region is relatively stable with small gradient changes.

[0161] For example, in a specific implementation, after the ground data center or the onboard autonomous processing system obtains the calculated raw velocity sequence, it first extracts the three-dimensional coordinates of the satellite's current orbital position at the same sampling time and converts them into geographic latitude. Then, the system determines in real time which specific orbital region the satellite is in based on the current latitude, and calls the preset window length adjustment conditions accordingly: when the system determines that the satellite has entered the high-latitude polar region, it automatically triggers the control command to reduce the window size, adaptively shortening the preset window length of the filter to a smaller target filtering window length, thereby improving the filter's sensitivity to changes in spatial position; when the system determines that the satellite is operating in a low-latitude region, it automatically triggers the control command to increase the window size, stretching the preset window length to a wider target filtering window length, thereby maximizing the filter's noise filtering capability. Next, the system adjusts the sliding step size of the core filter, using the filter with the instantaneous target filtering window length to perform sliding smoothing filtering on the raw velocity sequence in that region. Finally, the system reassembles the filtered data from each region according to the time sequence, and outputs a high-precision smoothed velocity sequence that eliminates high-frequency noise from the onboard receiver and perfectly preserves the details of local abrupt changes in the polar region.

[0162] This invention, through a pioneering spatial latitude-driven adaptive modulation technique for filter windows, effectively addresses the industry pain point of traditional fixed-length smoothing filtering mechanisms in thermospheric atmospheric inversion, where high-frequency denoising and abrupt signal preservation are mutually exclusive. The beneficial effects of this invention are threefold: First, it breaks the limitations of traditional technologies that use a single fixed filter window globally, by introducing orbital region-triggered adjustment conditions, enabling the filter to possess spatial awareness capabilities. Second, in polar regions, actively reducing the window length effectively avoids the malicious problem of over-smoothing causing the real, high-gradient atmospheric density rise signals caused by solar storms and geomagnetic storms to be mistakenly interpreted as measurement noise, ensuring the authenticity of polar region physical feature inversion. Third, in low-latitude regions, increasing the window length significantly suppresses thermal noise and multipath interference from conventional GNSS receivers on spacecraft, greatly improving the overall signal-to-noise ratio of the smoothed velocity sequence across the entire orbit, thus providing a highly spatially adaptive dynamic data foundation for subsequent high-precision, distortion-free numerical grid construction of thermospheric atmospheric density globally.

[0163] To enable those skilled in the art to better understand the embodiments of the present invention, an example is used below to illustrate the embodiments of the present invention.

[0164] refer to Figure 2 , Figure 2 This is a flowchart illustrating an atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites, provided in an embodiment of the present invention. Step 1: Data Stream Access Real-time acquisition of onboard GNSS observation data and attitude quaternions by receiving satellite downlink data through ground station network. And the solar panel rotation angle θ. A distributed message queue is used to process the data stream.

[0165] Step 2: Simplify dynamic orbit determination and smoothing Unlike the post-processing approach of traditional precise orbit determination (POD), this invention employs a simplified dynamic model to estimate satellite position and velocity. The calculated velocity sequence is then smoothed using a Savitzky-Golay filter. A key improvement lies in the dynamic adjustment of the filter window length based on the satellite's orbital position: the window is shortened in polar regions (where atmospheric changes are drastic) and extended in low-latitude regions, balancing real-time performance with smoothness.

[0166] Step 3: Refined aerodynamic modeling of individual components.

[0167] Establish a three-dimensional geometric model of the satellite body and movable solar array components; combine real-time attitude quaternions and solar array pointing angles to calculate the effective windward area of ​​the satellite body and solar arrays relative to the incoming flow direction. and .

[0168] 3.1) Body projection, based on attitude quaternions Calculate the plate normal vector, and combine it with the velocity vector. Calculate the windward area of ​​the body. .

[0169] 3.2) Solar wing projection, based on attitude quaternions The normal vector of the solar array is calculated using the solar array pointing angle α. The solar array has a large area but a thin thickness, and its projected area is sensitive to the angle of attack and altitude. If α is ignored, when the solar array changes from "tailwind" to "crosswind," the inverted density will produce a spurious increase, resulting in the synthesized total aerodynamic parameters that dynamically change with attitude.

[0170] Among them, C d,body The fundamental aerodynamic drag coefficient of a satellite is a dimensionless characteristic constant that characterizes the drag effect generated when the surface materials and geometry of a satellite exchange momentum with atmospheric particles in the thin thermosphere. d,wing The basic aerodynamic drag coefficient of a solar array is a dimensionless characteristic constant that characterizes the drag effect generated when the surface material and flat structure of a movable solar array exchange momentum with atmospheric particles in the thin thermosphere.

[0171] Step 4: Adaptive acceleration calculation based on attitude stochastic characteristics 4.1) Velocity differential acceleration, based on the satellite velocity obtained in step 2. The differential acceleration of velocity is calculated using the 5-point difference formula. :

[0172] Δt refers to the discrete time step of the velocity sequence sampling, that is, the time interval between two adjacent velocity observation data points.

[0173] This method can effectively suppress high-frequency noise to a certain extent while maintaining time resolution.

[0174] 4.2) Acceleration calculation and denoising based on the random characteristics of attitude: To further suppress differential amplification noise, the Tikhonov regularization algorithm is introduced to smooth the velocity sequence and obtain the total acceleration. .

[0175] The objective function of the Tikhonov regularized numerical differentiation algorithm Represented as:

[0176] in, The value of the differential acceleration is the velocity. For the acceleration to be determined, The total acceleration vector of the satellite at time t refers to the total acceleration vector under non-conservative forces, which is the resultant acceleration vector obtained by subtracting all conservative forces, such as higher-order gravitational fields, three-body gravity, and solar radiation pressure, from the total acceleration vector. For differential operators; the regularization parameter λ is inversely proportional to the attitude stability index, that is, the more unstable the satellite attitude, the weaker the regularization constraint, and vice versa.

[0177] Time series for calculating attitude quaternion change rate Regularization parameters and There is a positive correlation:

[0178] λ0 refers to the preset basic regularization parameter, which is the benchmark weight coefficient used in the numerical differential functional to constrain the smoothness of the acceleration curve when the satellite is in a completely stationary state (the variance of the attitude change rate is close to zero).

[0179] k refers to the adaptive regularization adjustment coefficient, which is a pre-set proportional gain constant used to quantify the influence of the time series variance of the attitude quaternion change rate on the regularization parameter mapping.

[0180] When a satellite is detected to be performing a "bias maneuver," the variance increases, and the regularization parameter is automatically reduced. To ensure the acceleration waveform is not distorted; increase during steady-state conditions. To reduce noise. Based on variance The feature adaptively adjusts the regularization parameter of the numerical derivative to extract the total acceleration vector. .

[0181] 4.3) Force model subtraction: The total satellite acceleration is the sum of conservative and non-conservative forces. Higher-order gravitational field accelerations are subtracted. (such as the EGM2008 gravity field model), three-body gravitational acceleration and solar radiation pressure acceleration Then, the drag acceleration is obtained. :

[0182] Step 5: Suppression of space-dependent thrust disturbances and vector decoupling Satellite A detects a sudden change in acceleration at a certain moment.

[0183] Using vector projection, Projected to The opposite direction is atmospheric drag acceleration. The residual component or anomalous amplitude component is used as the thrust estimate. And remove them.

[0184] α refers to the three-dimensional spatial angle between the total non-conservative force acceleration vector and the satellite's instantaneous velocity vector.

[0185] α1 refers to the preset first directional deviation angle threshold (e.g., 160° or 170°), which is used to determine whether the total acceleration vector of non-conservative forces deviates significantly from the theoretical reverse velocity direction (180° for pure atmospheric drag direction). When α>α1 and there are no sudden changes in the amplitude of neighboring satellites, it indicates that the included angle is close to 180°. At this time, the system determines that it is not affected by lateral or normal thrust.

[0186] α2 refers to the preset second directional deviation angle threshold (e.g., 90° or 120°), which is used to further screen and determine whether there is high-angle tilted lateral thrust interference or positive along-track thrust interference; When α < α2, it indicates that the non-conservative force vector has deviated towards the velocity direction or deviated significantly from the reverse velocity direction. Based on this, the system identifies and marks the corresponding thrust disturbance component.

[0187] Step 6: Atmospheric density inversion According to the atmospheric drag formula:

[0188] Where ρ is the atmospheric density. The aerodynamic coefficient, For windward area, For satellite quality, This is the satellite's velocity relative to the atmosphere. Because... Given that the atmospheric density ρ and V at the current location can be solved by inverse problem, we can obtain the answer. 2 rel It refers to the square of the magnitude of the satellite's relative velocity vector with respect to the atmosphere, that is, the square of the satellite's instantaneous relative velocity with respect to the atmosphere.

[0189] Step 7: Global Grid Fusion Leveraging the global distribution of giant constellation satellites, the retrieved discrete density values ​​are mapped onto a global latitude-longitude-height grid (e.g., the global grid is mapped by latitude and longitude). ,high Divide the grid at intervals ).exist At any given time, the satellite inversion density values ​​falling within each grid are statistically analyzed. A weighted average algorithm (with weights depending on thrust stripping residuals, orbit determination variance, and satellite attitude stability) is used to generate a real-time three-dimensional global atmospheric density field.

[0190] Average density within the grid The calculation is as follows:

[0191] N refers to the total number of valid observation data samples in the discrete density set that falls within the current target grid (i.e., the number of data in the pure density set, or the total number of target satellites after removing outliers).

[0192] ρi refers to the i-th local discrete atmospheric density value inverted by the i-th target satellite in the current target grid within the pure density set.

[0193] Wi refers to the fusion weight corresponding to the i-th target satellite, which is used to measure the relative reliability of the satellite's inverted density data in the weighted assimilation process.

[0194] At time t, the grid Internally received The density values ​​of each satellite. Calculate the mean of this data set. and standard deviation σ. If the atmospheric density inversion value ρ of a certain satellite s satisfy

[0195] If the data is affected by attitude maneuvers or multipath effects, it is determined that the data is discarded. The final output is the mean. This serves as the real-time density value for the grid. This method leverages the redundancy of multi-star data to reduce the random errors in single-star inversion. times ( (Number of satellites within the grid).

[0196] The final output is a global atmospheric density numerical grid model, with data coverage and confidence intervals labeled.

[0197] It should be noted that, for the sake of simplicity, the method embodiments are all described as a series of actions. However, those skilled in the art should understand that the embodiments of the present invention are not limited to the described order of actions, because according to the embodiments of the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to the embodiments of the present invention.

[0198] Reference Figure 3 The diagram illustrates a structural block diagram of an atmospheric density modeling device based on the decoupling of aerodynamic parameters of a giant constellation satellite, provided in an embodiment of the present invention. Specifically, it may include the following modules: The real-time data acquisition module 301 is used to receive real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles from satellite downlink. The total aerodynamic parameters and atmospheric drag acceleration determination module 302 is used to determine the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite, through the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. Satellite feature determination module 303 is used to determine the mass of the satellite and the relative velocity of the satellite relative to the atmosphere; Model building module 304 is used to build a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. The target data acquisition module 305 is used to acquire the target global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle of the spacecraft under test; The orbital warning result generation module 306 is used to call the global atmospheric density numerical grid model and generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

[0199] As the device embodiment is basically similar to the method embodiment, the description is relatively simple, and relevant parts can be found in the description of the method embodiment.

[0200] In addition, embodiments of the present invention also provide an electronic device, such as... Figure 4 As shown, it includes a processor 401, a communication interface 402, a memory 403, and a communication bus 404, wherein the processor 401, the communication interface 402, and the memory 403 communicate with each other through the communication bus 404. Memory 403 is used to store computer programs; When the processor 401 executes the program stored in the memory 403, it implements the atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites as described in any of the above embodiments: The communication bus mentioned above can be a Peripheral Component Interconnect (PCI) bus or an Extended Industry Standard Architecture (EISA) bus, etc. This communication bus can be divided into address bus, data bus, control bus, etc. For ease of illustration, only one thick line is used to represent it in the diagram, but this does not mean that there is only one bus or one type of bus.

[0201] The communication interface is used for communication between the aforementioned terminal and other devices.

[0202] The memory may include random access memory (RAM) or non-volatile memory, such as at least one disk storage device. Optionally, the memory may also be at least one storage device located remotely from the aforementioned processor.

[0203] The processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0204] like Figure 5 As shown, in another embodiment of the present invention, a computer-readable storage medium 501 is also provided, which stores instructions that, when run on a computer, cause the computer to execute the atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites described in the above embodiment.

[0205] This application embodiment also provides a chip, which includes a processor and a communication interface. The communication interface and the processor are coupled. The processor is used to run programs or instructions to implement the various processes of the above-described embodiment of the atmospheric density modeling method based on the decoupling of aerodynamic parameters of giant constellation satellites, and can achieve the same technical effect. To avoid repetition, it will not be described again here.

[0206] It should be understood that the chip mentioned in the embodiments of this application may also be referred to as a system-on-a-chip, system chip, chip system, or system-on-a-chip, etc.

[0207] It should be noted that, in this document, the terms include, encompass, or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that includes a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "including a…" does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element. Furthermore, it should be noted that the scope of the methods and apparatuses in the embodiments of this application is not limited to performing functions in the order shown or discussed, but may also include performing functions substantially simultaneously or in the reverse order, depending on the functions involved. For example, the described methods may be performed in a different order than described, and various steps may be added, omitted, or combined. Additionally, features described with reference to certain examples may be combined in other examples.

[0208] Through the above description of the embodiments, those skilled in the art can clearly understand that the methods of the above embodiments can be implemented by means of software plus necessary general-purpose hardware platforms. Of course, they can also be implemented by hardware, but in many cases the former is a better implementation method. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) and includes several instructions to cause a terminal (which may be a mobile phone, computer, server, air conditioner, or network device, etc.) to execute the methods described in the various embodiments of this application.

[0209] The embodiments of this application have been described above with reference to the accompanying drawings. However, this application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of this application.

Claims

1. An atmospheric density modeling method based on decoupling of aerodynamic parameters of giant constellation satellites, characterized in that, include: Receives real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles from satellite downlink; The total aerodynamic parameters characterizing the gas dynamics interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite are determined by the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. Determine the mass of the satellite and its relative velocity with respect to the atmosphere; A global atmospheric density numerical grid model is constructed using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle for the spacecraft under test; The global atmospheric density numerical grid model is invoked to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

2. The method according to claim 1, characterized in that, The step of determining the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite using the real-time observation data of the spaceborne global navigation satellite system, the attitude quaternion, and the solar array rotation angle includes: Using observation data from the aforementioned onboard global navigation satellite system, a pre-trained simplified dynamic model is invoked to determine the satellite's velocity sequence; The velocity sequence is smoothed to obtain a smoothed velocity sequence; Determine the velocity vector corresponding to the smoothed velocity sequence; The first angle between the plate normal vector and the velocity vector is calculated using the attitude quaternion and the velocity vector, and the effective windward area of ​​the satellite relative to the direction of incoming flow is calculated based on the first angle. The normal vector of the wing surface is calculated using the attitude quaternion and the rotation angle of the solar wing, and the effective windward area of ​​the movable solar wing relative to the direction of the incoming flow is calculated based on the normal vector of the wing surface. The effective windward area of ​​the main body and the effective windward area of ​​the solar array are respectively combined with the corresponding aerodynamic coefficients and weighted summed to obtain the total aerodynamic parameters that change dynamically with attitude.

3. The method according to claim 2, characterized in that, The atmospheric drag acceleration is determined as follows: Using the smoothed velocity sequence, the differential acceleration of the velocity at the target time point is calculated; Determine the time series variance of the rate of change of the attitude quaternion; Regularization parameters are adaptively generated based on the time series variance and a preset formula. Construct a numerical differential objective function that includes the regularization parameter, the velocity differential acceleration, and a preset differential operator; The total acceleration vector is calculated by calling the numerical differential objective function and using the smoothed velocity sequence; The higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are determined from the total acceleration vector, and the higher-order gravitational field acceleration, the three-body gravitational acceleration, and the solar radiation pressure acceleration are deleted to obtain the non-conservative force total acceleration vector. The direction determination result is generated by using the angle between the direction of the non-conservative force total acceleration vector and the velocity vector; The acceleration amplitude of the satellite constellation in the adjacent orbital plane is obtained, and the acceleration amplitude is compared with the total non-conservative force acceleration vector to generate a spatial consistency judgment result. Based on the direction determination results and the spatial consistency determination results, the lateral thrust interference, normal thrust interference, and along-track thrust interference in the non-conservative total force acceleration vector are identified. Using a vector decomposition algorithm, the projection result of the total non-conservative force acceleration vector is generated in the opposite direction of the velocity vector. The lateral thrust interference, the normal thrust interference, and the path-direction thrust interference are then removed from the projection result to generate atmospheric drag acceleration in the opposite velocity direction.

4. The method according to claim 3, characterized in that, The steps of constructing a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity include: The discrete atmospheric density value of the satellite's current position is calculated using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. Multiple target satellites are identified, and the discrete atmospheric density values ​​corresponding to the target satellites at the same time are determined as the target discrete atmospheric density values; The target discrete atmospheric density value is mapped to each target grid of a preset global latitude, longitude and altitude grid, and the discrete density set falling into the target grid is statistically analyzed. Calculate the mean and standard deviation of each discrete density set, and based on the mean and standard deviation, remove abnormal density data in the discrete density set whose deviation reaches a preset condition, and select the pure density set corresponding to the target grid. The fusion weight of each target satellite is calculated, and the pure density set is fused using a weighted average algorithm based on the fusion weight. Finally, a global atmospheric density numerical grid model labeled with data coverage and confidence interval is generated and output.

5. The method according to claim 4, characterized in that, The step of calling the global atmospheric density numerical grid model to generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle includes: Acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle for the spacecraft under test; Using observation data from the target global navigation satellite system, the predicted orbital position sequence and predicted velocity sequence of the spacecraft under test within a future preset time period are calculated by extrapolation through a simplified dynamic model. The predicted three-dimensional grid coordinates corresponding to the predicted orbit position sequence are input into the global atmospheric density numerical grid model for spatial interpolation retrieval, and the spatiotemporal background atmospheric density prediction values ​​of each point on the predicted orbit position sequence are matched. The target body normal vector is calculated using the target attitude quaternion and the velocity vector direction of the predicted velocity sequence, and the predicted effective windward area of ​​the spacecraft under test in a future preset time period is calculated based on the target body normal vector. The solar array pointing normal vector is determined by using the target attitude quaternion and the target solar array rotation angle, and the predicted effective windward area of ​​the solar array of the spacecraft under test in the future preset time period is calculated based on the solar array pointing normal vector. The predicted effective windward area of ​​the satellite body and the predicted effective windward area of ​​the solar array are multiplied by the preset basic drag coefficients of the satellite body and the solar array, respectively, and then weighted and summed to calculate the predicted total dynamic aerodynamic parameters of the spacecraft under test in the future preset time period. Using the predicted atmospheric density of the spatiotemporal background, the predicted velocity sequence, and the predicted total dynamic aerodynamic parameters, the predicted atmospheric drag acceleration sequence that the spacecraft under test will experience in a future preset time period is calculated. The predicted atmospheric drag acceleration sequence is used as a non-conservative force input term for numerical integration, and the resulting orbital warning result for the spacecraft under test is output.

6. The method according to claim 2, characterized in that, The step of smoothing the velocity sequence to obtain a smoothed velocity sequence includes: Obtain the current orbital position of the satellite; The orbital region of the satellite is determined based on its current orbital position, and adjustment conditions for a preset window length are determined based on the orbital region. The adaptively adjusted target filter window length is determined based on the aforementioned adjustment conditions; The velocity sequence is smoothed by using a filter with the target filtering window length to obtain a smoothed velocity sequence.

7. The method according to claim 6, characterized in that, The adjustment conditions are as follows: When the orbital region is determined to be in a polar region, the preset window length is reduced; or, When the orbital region is determined to be in a low-latitude region, the length of the preset window is increased.

8. An atmospheric density modeling device based on decoupling of aerodynamic parameters of giant constellation satellites, characterized in that, include: The real-time data acquisition module is used to receive real-time satellite-borne global navigation satellite system observation data, attitude quaternions, and solar array rotation angles from satellite downlink. The total aerodynamic parameters and atmospheric drag acceleration determination module is used to determine the total aerodynamic parameters characterizing the gas dynamic interaction response capability of the satellite, and the atmospheric drag acceleration characterizing the deceleration effect per unit mass of the satellite, through the real-time observation data of the satellite-borne global navigation satellite system, the attitude quaternion, and the solar array rotation angle. A satellite feature determination module is used to determine the mass of the satellite and its relative velocity with respect to the atmosphere; The model building module is used to construct a global atmospheric density numerical grid model using the atmospheric drag acceleration, the total aerodynamic parameters, the satellite mass, and the relative velocity. The target data acquisition module is used to acquire the target's global navigation satellite system observation data, target attitude quaternion, and target solar array rotation angle of the spacecraft under test; The orbital warning result generation module is used to call the global atmospheric density numerical grid model and generate orbital warning results for the spacecraft under test based on the target global navigation satellite system observation data, the target attitude quaternion, and the target solar array rotation angle.

9. An electronic device, characterized in that, It includes a processor, a memory, and a program or instructions stored in the memory and executable on the processor, wherein the program or instructions, when executed by the processor, implement the method as described in any one of claims 1-7.

10. A readable storage medium, characterized in that, The readable storage medium stores a program or instructions that, when executed by a processor, implement the method as described in any one of claims 1-7.

Images