Low earth orbit constellation multi-satellite cooperative decommissioning and deorbiting method

Abstract

This invention relates to the field of spacecraft orbit control and discloses a method for coordinated deorbiting of multiple satellites in a low-Earth orbit constellation. The method includes: the satellite to be deorbited constructing a deorbit crossing envelope and sending deorbit intention data frames to physically overlapping associated satellites; the associated satellites generating a double-layer exclusion corridor parameter containing an inner rigid safety boundary and an outer elastic service boundary; the satellite to be deorbiting generating a global deorbit cost function with a dynamic penalty term; when the outer elastic boundary is a hard constraint with no feasible solution, controlled contraction is performed only on the outer boundary, and the contraction amount is converted into a numerical penalty, while maintaining the inner rigid boundary and solving for a cooperative compromise solution; the compromise solution is converted into an action sequence for execution, and in case of an anomaly, the command is revoked and a minimum orbit decay process is initiated. This invention decouples collision avoidance conditions in layers, expands the solution space through constraint relaxation, and balances the physical collision avoidance baseline with on-orbit service continuity.

Description

Technical Field

[0001] This invention relates to the field of spacecraft orbit control, specifically a method for coordinated deorbiting of multiple satellites in a low-Earth orbit constellation. Background Technology

[0002] As the scale of low-Earth orbit satellite constellations expands, the decommissioning and deorbiting of satellites nearing the end of their lifespan has become a routine operation for maintaining orbital resources. During the satellite deorbiting process, the satellite to be deorbited crosses multiple orbital planes, resulting in multiple encounters with normally orbiting satellites within the constellation, thus posing a collision risk.

[0003] Currently, conventional deorbit avoidance strategies primarily involve setting a uniform safety distance envelope for satellites and using it as a hard constraint for trajectory optimization. In rendezvous scenarios with densely packed space nodes within a constellation, the satellite to be deorbited simultaneously faces space constraints from multiple in-orbit satellites. A single hard safety boundary constraint creates mutual compression within the optimization model, leading to solution space shrinkage or even a situation where no feasible solution exists. To obtain a feasible trajectory, existing methods typically indiscriminately reduce the safety boundaries of all associated satellites. This approach not only violates the physical collision avoidance baseline but also interferes with the routine mission execution of surrounding in-orbit satellites, causing damage to the overall operational continuity of the constellation.

[0004] Furthermore, the aging of satellite platform components during the decommissioning phase poses a risk of sudden malfunctions during deorbit maneuvers. Current deorbit control procedures focus on early trajectory planning, lacking command interception and intervention procedures during the execution phase. When anomalies occur during satellite deorbiting, particularly attitude control failures of varying degrees, existing methods lack corresponding graded physical degradation and fallback measures. This can lead to decommissioned satellites with impaired maneuverability remaining in orbit, evolving into space debris and posing a threat to other spacecraft at the same orbital altitude. Summary of the Invention

[0005] The technical problem solved by this invention is that during the decommissioning and deorbiting process of existing low-Earth orbit constellation satellites, due to the high density of space nodes in the constellation, there is a risk of rendezvous and collision between the satellite to be deorbited and the satellites in orbit. At the same time, existing deorbiting avoidance strategies usually adopt a single hard safety boundary constraint. In complex scenarios with multiple satellites densely intersecting, the optimization process often becomes infeasible due to the limited solution space. It is difficult to balance physical collision avoidance safety with the continuity of constellation operation services, and it lacks adaptive coordination and multi-level backup actions to deal with sudden anomalies.

[0006] To address the above problems, the present invention provides the following technical solution:

[0007] This invention provides a method for the coordinated decommissioning and deorbiting of multiple satellites in a low-Earth orbit constellation, comprising the following steps:

[0008] The satellite to be deorbited generates initial deorbit candidate schemes, and the deorbit crossing envelope is constructed by extrapolating orbital dynamics.

[0009] Select the associated satellites that physically overlap with the deorbit crossing envelope, and send the deorbit intention data frame in a directional manner;

[0010] The associated satellites assess rendezvous risks and generate a two-layer four-dimensional spatiotemporal exclusion corridor parameter that includes an inner rigid safety boundary and an outer elastic service boundary.

[0011] The satellite to be deorbited receives the double-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data, and generates a constellation-level global deorbit cost function with a dynamic constraint relaxation penalty term.

[0012] When the outer elastic service boundary has no feasible solution as a hard constraint, controlled contraction is performed only on the outer elastic service boundary. The contraction amount is converted into a numerical penalty and injected into the constellation-level global deorbit cost function. The inner rigid safety boundary is kept from relaxation and a cooperative compromise solution is obtained.

[0013] The cooperative compromise solution is converted into an action sequence and sent to the associated satellite for execution. If an anomaly is confirmed, the action sequence that has not yet entered the irreversible execution window is canceled, and the process of ensuring orbit decay is initiated.

[0014] Furthermore, the specific steps for generating initial deorbit candidate schemes for the satellite to be deorbited and constructing the deorbit crossing envelope through orbital dynamics extrapolation include:

[0015] The satellite to be deorbited performs orbital nonlinear dynamics integral calculations based on the initial deorbit candidate schemes to obtain the predicted position vector during the deorbiting process;

[0016] The orbital dynamics state transition matrix is ​​calculated by extrapolation using orbital dynamics, and the time evolution sequence of the position error covariance matrix corresponding to the predicted position vector is extracted.

[0017] By combining the time evolution sequence of the predicted position vector and the position error covariance matrix, a spatial error boundary covering the descending trajectory is constructed in three-dimensional space, and the derailment crossing envelope is constructed using the three-dimensional derailment crossing envelope calculation formula.

[0018] Furthermore, the specific steps for filtering out associated satellites that physically overlap with the deorbiting crossing envelope and directionally sending deorbiting intention data frames include:

[0019] Extract the constellation's global approximate ephemeris data stored on-board and estimate the space occupancy range of the remaining active satellite nodes within the constellation;

[0020] Nodes whose spatial occupancy range and deorbit crossing envelope physically overlap are identified, and these physically overlapping nodes are confirmed as associated satellites.

[0021] The node physical identifier, the timestamp sequence corresponding to the system's global time variable, the time evolution sequence of the predicted position vector of the satellite to be deorbited, and the time evolution sequence of the position error covariance matrix are combined and packaged to generate and send a deorbit intention data frame.

[0022] Furthermore, the specific steps for assessing the rendezvous risk of the associated satellites include:

[0023] The satellite analyzes the received deorbit intention data frames and performs a time-parameterized relative trajectory search within a preset rendezvous time window to calculate the closest approach time.

[0024] Extract the relative position and relative velocity near the closest approach time, and superimpose the position error covariance of the satellite to be deorbited with its own position error covariance to generate a joint position error covariance matrix;

[0025] The joint position error covariance matrix is ​​symmetric and positive definite by using the joint position error covariance matrix symmetry formula and the joint position error covariance matrix positive definite formula. The rendezvous risk is then quantitatively assessed based on the processing results.

[0026] Furthermore, the specific steps for generating the parameters of the two-layer four-dimensional spatiotemporal repulsion corridor, which includes an inner rigid safety boundary and an outer elastic service boundary, include:

[0027] Based on the rendezvous risk and the equivalent radius of satellite physical collision avoidance obtained from the assessment, the inner rigid safety boundary is established using the inner rigid safety boundary shape matrix calculation formula;

[0028] Combining task interference avoidance margin, an outer elastic service boundary that can be controlled to shrink within a preset range is constructed outside the inner rigid safety boundary using the outer elastic service boundary shape matrix calculation formula;

[0029] The inner rigid safety boundary and the outer elastic service boundary are combined and encapsulated with their own position state sequences to generate a two-layer four-dimensional spatiotemporal repulsion corridor parameter.

[0030] Furthermore, the specific steps for the satellite to be deorbited to receive the dual-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data, and generate a constellation-level global deorbit cost function with a dynamic constraint relaxation penalty term, include:

[0031] The parameters of the two-layer four-dimensional spatiotemporal exclusion corridor for satellites to be deorbited are analyzed, and heterogeneous cost weights are assigned to each satellite node in combination with node cost data.

[0032] By combining the various physical consumption quantities and cost weights and normalizing them, the weighted cost value of each node is calculated.

[0033] The weighted maneuver cost of the satellite to be deorbited, the weighted secondary generation cost of the associated satellite set, and the dynamic constraint relaxation penalty term are integrated, and a constellation-level global deorbit cost function is generated using the constellation-level global deorbit cost function calculation formula.

[0034] Furthermore, before performing controlled shrinkage only on the outer elastic service boundary when no feasible solution is found as a hard constraint, the method further includes:

[0035] The outer elastic service boundary is transformed into a spatial avoidance matrix constraint for optimization. If no feasible solution is detected continuously, the Lagrange multiplier values ​​of each constraint in the last iteration are extracted.

[0036] By combining the values ​​of Lagrange multipliers, the amount of constraint violation, and the constraint scale normalization factor, the stagnation score of each associated satellite is calculated using the stagnation score calculation formula.

[0037] The associated satellite with the highest blocking score is identified as the core blocking node, and controlled contraction is prepared to be performed on the core blocking node.

[0038] Furthermore, the specific steps of performing controlled contraction only on the outer elastic service boundary, converting the contraction amount into a numerical penalty injected into the constellation-level global deorbit cost function, maintaining the inner rigid safety boundary from relaxation, and obtaining a cooperative compromise solution include:

[0039] Controlled contraction is performed on the outer elastic service boundary for the core blocking node. The outer boundary insufficient space is calculated using the formula for calculating the outer boundary insufficient space, and the contraction amount is obtained.

[0040] The inner boundary retention constraint formula is used to ensure that the inner rigid safety boundary cannot be relaxed, while the shrinkage amount is converted into a numerical penalty by increasing the dynamic penalty benchmark coefficient.

[0041] Numerical penalties are injected into the constellation-level global deorbit cost function and the optimization process is restarted until the convergence condition is met to obtain a cooperative compromise solution.

[0042] Furthermore, the specific steps of converting the cooperative compromise solution into an action sequence and sending it to the associated satellite for execution include:

[0043] The cooperative compromise solution is analyzed and then converted into discrete control instructions applicable to different hardware configurations.

[0044] For associated satellites using active orbit control, discrete control commands are converted into ignition pulse commands, and corresponding action timing sequences are generated and sent to the associated satellites for execution.

[0045] For correlated satellites using differential drag control, discrete control commands are converted into aerodynamic drag control commands, and corresponding action timings are generated and sent to the correlated satellites for execution.

[0046] Furthermore, the specific steps for canceling the action sequence that has not yet entered the irreversible execution window when an anomaly is confirmed, and switching to the minimum track attenuation process, include:

[0047] When the system monitoring log generates an anomaly flag, the action sequence that has not yet entered the irreversible execution window is forcibly intercepted and revoked.

[0048] If the attitude of the satellite to be deorbited is controllable, the orientation of the drag-increasing components is adjusted to increase the equivalent windward area along the velocity direction, forming a continuous increase in aerodynamic drag to reduce the orbital altitude and enter the minimum orbital attenuation process.

[0049] If the attitude of the satellite to be deorbited completely fails, the passive deorbiting device will be triggered to enter the ballistic decay trajectory, downgrade to a one-way passive avoidance mode, and switch to the minimum orbit decay process.

[0050] This invention provides a method for the coordinated decommissioning and deorbiting of multiple satellites in a low-Earth orbit constellation. It offers the following advantages:

[0051] 1. This invention generates a two-layer, four-dimensional spatiotemporal exclusion corridor parameter, comprising an inner rigid safety boundary and an outer elastic service boundary, thereby decoupling collision avoidance conditions in a layered manner. The inner rigid boundary ensures the physical collision avoidance baseline, while the outer elastic boundary takes into account the service continuity of on-orbit satellites. This layered structure overcomes the limitations of a single rigid safety boundary in complex rendezvous scenarios, ensuring that the underlying physical collision avoidance constraints are not exceeded during optimization calculations, thus improving the safety of multi-satellite coordinated deorbiting processes.

[0052] 2. When multi-satellite intersection leads to an infeasible solution due to the outer elastic service boundary as a hard constraint, this invention utilizes Lagrange multipliers to calculate the stagnation score to locate the core stagnation node. Controlled contraction is then performed only on the outer boundary of this node, and the contraction amount is converted into a numerical penalty injected into the global deorbit cost function. This step transforms the unsolvable hard-constraint optimization problem into a soft-constraint optimization problem with a dynamic penalty term. By relaxing constraints, the solution space is expanded, leading to a feasible deorbit compromise solution, while controlling the overall constellation's service performance loss.

[0053] 3. This invention introduces a multi-level anomaly degradation and backup process during the deorbiting phase, intercepting and canceling abnormal actions by monitoring the irreversible execution window. Based on this, and according to the controllable attitude of the satellite to be deorbited, it selectively initiates a drag-increasing process to adjust the orientation of drag-increasing components or triggers a ballistic decay process using a passive deorbiting device. This tiered physical backup mechanism allows the satellite to escape its operational orbit even in the event of a sudden anomaly or failure, relying on the natural decay mechanism of its orbit, thus reducing the risk of long-term space debris accumulation. Attached Figure Description

[0054] Figure 1 This is a diagram of the low-Earth orbit constellation multi-satellite collaborative system architecture of the present invention;

[0055] Figure 2 This is a flowchart illustrating the overall workflow of the low-Earth orbit constellation multi-satellite coordinated deorbiting method of the present invention.

[0056] Figure 3 This is a diagram of the three-dimensional spatial boundary envelope and traversal trajectory according to an embodiment of the present invention;

[0057] Figure 4 This is a contour distribution of the cost function and an optimization trajectory diagram according to an embodiment of the present invention. Detailed Implementation

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

[0059] See attached document Figure 1 The operating environment for implementing the low-Earth orbit constellation multi-satellite coordinated deorbiting method of the present invention includes the satellite to be deorbited, associated satellites, and inter-satellite communication network.

[0060] The satellite to be deorbited is configured to trigger a decommissioning and orbit reduction procedure and generate a deorbiting intention data frame; the associated satellite is configured to receive the deorbiting intention data frame, assess the rendezvous risk with the satellite to be deorbited, and provide feedback on the two-layer four-dimensional spatiotemporal exclusion corridor parameters; the inter-satellite communication network is configured to transmit the deorbiting intention data frame, the two-layer four-dimensional spatiotemporal exclusion corridor parameters, and action timing data between the satellite to be deorbited and the associated satellite to establish an inter-satellite data connection.

[0061] The satellites to be deorbited and their associated satellites are all on-orbit nodes in a low-Earth orbit constellation. To support the multi-satellite coordinated deorbit planning, the on-orbit nodes participating in the coordination are equipped with onboard computing units, inter-satellite communication transceiver units, attitude control actuators, and orbit control actuators.

[0062] The onboard computing unit is used to perform orbital dynamics extrapolation, rendezvous risk assessment, associated node screening, and finite-dimensional constraint optimization solutions; the inter-satellite communication transceiver unit is used to send and receive deorbit intention data frames, exclusion corridor parameters, and confirmed action sequences; the attitude control actuator is used to adjust the satellite's attitude to form a predetermined thrust direction or aerodynamic drag configuration; and the orbit control actuator is used to perform active orbital maneuvers according to the action sequence.

[0063] For small satellite nodes with limited computing power, the collaborative optimization solution can be completed by the satellite to be deorbited or the designated master satellite node; if the satellite-to-ground link is available, the ground-aided computing system can also provide the initial solution, parameter injection, or offline verification results. The ground-aided computing system is not a necessary component for the on-orbit collaborative execution of this method. Other associated satellites perform rendezvous risk assessment, constraint parameter feedback, and action timing confirmation.

[0064] For microsatellite nodes with differential drag control capabilities, the orbit control actuator and attitude control actuator can be partially merged in terms of control functions, generating an increase in aerodynamic drag by changing the satellite's frontal area.

[0065] The system has a preset failure condition baseline. When the active orbit control actuator of the satellite to be deorbited fails, the inter-satellite communication transceiver unit is interrupted, or the coordinated action sequence cannot be confirmed in a closed loop, the system stops the multi-satellite coordinated planning process and switches to the minimum orbit decay process. Under the minimum orbit decay process, the system degrades to a one-way passive avoidance mode.

[0066] If the attitude control actuators of the satellite to be deorbited remain controllable, the satellite will increase its frontal area through attitude adjustments, forming a differential drag configuration or a drag-enhanced deorbit configuration. If the attitude control actuators also fail, the satellite will extrapolate its passive decay trajectory according to a ballistic trajectory. When the satellite to be deorbited is unable to continue transmitting updated data, the associated satellites will independently extrapolate the passive decay trajectory of the satellite to be deorbited based on the most recently received deorbit intention data frame, constellation-shared ephemeris, ground-injected space situational awareness data, or local observation data, and calculate a local avoidance scheme accordingly.

[0067] In the system's logical role definition, the role of a satellite node is dynamically determined based on its operational status and mission instructions. A satellite node that enters a planned decommissioning state and initiates a deorbiting procedure is defined as a satellite awaiting deorbiting. Within the spatial range traversed by the satellite awaiting deorbiting, on-orbit satellites whose predicted trajectories overlap with its own are defined as associated satellites. When an associated satellite reaches its design life and receives a decommissioning instruction, the logical role of that node switches to that of a satellite awaiting deorbiting.

[0068] The system establishes a unified spatiotemporal reference. The spatial reference adopts the Earth inertial coordinate system, and the position and velocity parameters of the satellite to be deorbited and its associated satellites are all mapped to the Earth inertial coordinate system for calculation. The time reference adopts the constellation's globally synchronized timestamp, and the internal clocks of each satellite node within the constellation are periodically aligned through GNSS timing, time synchronization frames from the inter-satellite communication network, or time parameters injected from the ground.

[0069] Each node writes a unified timestamp into the deorbit intention data frame and repulsion corridor parameters, and limits the time synchronization error to within a preset threshold. The preset threshold is determined based on the constellation's relative velocity, the rendezvous time window length, and the precision of the action timing control. When a short-term interruption occurs in the inter-satellite communication network, each node uses the time reference of the most recent synchronization and the local clock drift model to maintain trajectory prediction consistency for a limited time.

[0070] See attached document Figure 2 This invention provides a method for the coordinated decommissioning and deorbiting of multiple satellites in a low-Earth orbit constellation, comprising the following steps.

[0071] S1, Initial Deorbit Intent Generation and Envelope Calculation. The satellite to be deorbited generates initial deorbit candidate schemes, and based on these initial deorbit candidate schemes, an deorbit crossing envelope is constructed through orbital dynamics extrapolation.

[0072] S2, Targeted Screening and Data Distribution. The satellite to be deorbited compares its ephemeris data with that of surrounding nodes, filters out associated satellites that physically overlap with the deorbiting envelope, and then sends a packaged deorbiting intent data frame to the selected associated satellites.

[0073] S3, Risk Assessment and Corridor Generation. After receiving the deorbiting intention data frame, the associated satellite assesses the rendezvous risk with the satellite to be deorbited and generates a two-layer four-dimensional spatiotemporal exclusion corridor parameter that includes an inner rigid safety boundary and an outer elastic service boundary.

[0074] S4, Global Cost Construction. The satellite to be deorbited receives two-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data, which are then integrated to generate a constellation-level global deorbit cost function with dynamic constraint relaxation penalty terms.

[0075] S5, Flexible Degradation and Replanning. When the satellite to be deorbited fails to find a continuous feasible solution by treating the outer elastic service boundary as a hard constraint, controlled contraction is performed only on the outer elastic service boundary. The contraction amount is then converted into a numerical penalty and injected into the constellation-level global deorbit cost function. Under the premise that the inner rigid safety boundary cannot be relaxed, a cooperative compromise solution is obtained.

[0076] S6, Scheme Confirmation and Guaranteed Execution. The satellite to be deorbited will convert the cooperative compromise solution into the corresponding action sequence and send it to the associated satellites so that the corresponding orbit control actuators can perform physical maneuvers; in the event of communication abnormalities or the inability to confirm the action sequence in a closed loop, the action sequence that has not yet entered the irreversible execution window will be terminated or canceled, and the guaranteed orbit decay process will be initiated.

[0077] The technical details and mathematical models of the above steps will be explained below.

[0078] After the satellite to be deorbited triggers the decommissioning procedure, the system executes step S1, which is the initial deorbiting intention generation and envelope calculation step. This step is used to form the predicted trajectory of the satellite to be deorbited during the descent process and its spatial uncertainty boundary, and to provide a data foundation for subsequent screening of associated satellites and risk assessment. This step includes the following sub-steps.

[0079] S101, extrapolation calculation of nonlinear dynamics of predicted trajectory.

[0080] In this embodiment, the onboard computing unit of the satellite to be deorbited extracts the initial space state vector and remaining propellant mass data currently in orbit, and generates initial deorbit candidate schemes accordingly. The initial space state vector includes at least the position and velocity parameters of the satellite to be deorbited in the Earth's inertial coordinate system; the remaining propellant mass data is used to constrain subsequent deorbit maneuver capabilities.

[0081] In the low Earth orbit (LEO) environment, factors such as thin atmospheric drag, Earth's non-spherical gravitational perturbation, and thrust execution deviations cause the actual trajectory of a satellite undergoing deorbiting maneuvers to deviate from an ideal two-body orbit. The onboard computing unit calculates continuous thrust control input commands based on pre-stored control strategies and performs orbital nonlinear dynamic integration calculations in the Earth inertial coordinate system based on these commands, obtaining a three-dimensional spatial coordinate sequence distributed over time during the deorbiting process. This three-dimensional spatial coordinate sequence constitutes the predicted position vector of the satellite to be deorbited.

[0082] For modeling the perturbation forces of the space environment and selecting the parameters of the numerical integrator, those skilled in the art can perform adaptive calculations based on the actual constellation's orbital altitude, attitude and orbit control accuracy, and onboard computing resources. For example, a standard fourth-order Runge-Kutta numerical integrator can be used for step-by-step solutions. The numerical integral extrapolation of orbital dynamics is a well-known technique in this field and will not be elaborated here.

[0083] S102, Initialization and time propagation of the orbital state error covariance matrix.

[0084] To quantify the spatial position deviation of the predicted position vector of the satellite to be deorbited due to environmental perturbation uncertainties, navigation and positioning errors, and thrust execution errors, the onboard computing unit sets an initial state error covariance matrix based on the current navigation and positioning measurement accuracy.

[0085] Specifically, the initial state error covariance matrix can be set as a six-dimensional state covariance matrix, whose diagonal elements are mainly based on the three-axis position variance and three-axis velocity variance output by the GNSS receiver on the satellite or the ground orbit determination results. At the same time, the attitude and acceleration measurement errors introduced by the star sensor and the inertial measurement unit can be mapped into position and velocity uncertainties through the orbital dynamics model and superimposed on the initial state error covariance matrix.

[0086] The onboard computing unit utilizes the orbital dynamics state transition matrix to propagate and solve the initial state error covariance matrix along the time axis, obtaining a six-dimensional state error covariance matrix that evolves over time. The formula for calculating its covariance time propagation is as follows:

[0087] ;

[0088] In the formula, express The six-dimensional state error covariance matrix at time t; Indicates the initial time; Represents the initial state error covariance matrix; Indicates from arrive The orbital dynamics state transition matrix; This represents the transpose of the state transition matrix; This represents the system process noise covariance matrix caused by unmodeled environmental perturbations and control execution noise during this time period.

[0089] When constructing the three-dimensional deorbiting pass envelope, the onboard computing unit uses the position error covariance matrix extraction formula, from... Extracting the position error covariance matrix The formula for extracting the position error covariance matrix is:

[0090] ;

[0091] In the formula, This represents the position error covariance matrix of the satellite to be deorbited; This represents the selection matrix for extracting the three-dimensional position components from the six-dimensional state vector. In this embodiment, Specifically The block matrix is ​​used to project the six-dimensional state space onto the three-dimensional position space; Represents a third-order identity matrix; Represents a third-order zero matrix; This represents the six-dimensional state error covariance matrix of the satellite to be deorbited; This represents the transpose of the selection matrix used to extract the three-dimensional position components from the six-dimensional state vector. This represents the matrix transpose operator. This represents the matrix product term used to project the three-dimensional position error covariance matrix from the six-dimensional state error covariance matrix.

[0092] S103, Safety Margin Addition and Off-Track Crossing Envelope Construction.

[0093] As state errors accumulate and spread with increasing orbit extrapolation time, a statistically significant spatial error boundary forms around the predicted trajectory of the satellite to be deorbited. The onboard computing unit utilizes the predicted position vector of the satellite to be deorbited, the position error covariance matrix obtained from propagation calculations, and a pre-set confidence interval to construct a spatial error boundary in three-dimensional physical space that envelops the deorbiting trajectory, thus forming the deorbit crossing envelope.

[0094] The envelope model is generated using the three-dimensional off-track crossing envelope calculation formula, which is as follows:

[0095] ;

[0096] In the formula, This represents the three-dimensional deorbiting envelope of the satellite to be deorbited; Represents an arbitrary position vector in three-dimensional physical space; Symbols indicating membership of elements in a set; Represents the three-dimensional set of real numbers; Indicates a condition separator; This represents the predicted position vector of the satellite to be deorbited; Represents a system-wide time variable; This represents the matrix transpose operator. This represents the position error covariance matrix of the satellite to be deorbited; This represents the inverse matrix of the position error covariance matrix of the satellite to be deorbited; Indicates the threshold of the envelope expansion factor; Represents the position deviation vector; The transpose of the position deviation vector; Represents the squared term of the distance from Maharanobis; The symbol represents the less than or equal to relation.

[0097] This formula is used to limit the system's global time variable. The condition that the squared term of the Mahalanobis distance is not greater than the threshold of the envelope expansion factor is met. The set of three-dimensional spatial points is used to obtain the uncertainty boundary of the space occupancy of the satellite to be deorbited at that moment.

[0098] Envelope expansion factor threshold The physical meaning of is the scale parameter that controls the collision avoidance alarm limit in three-dimensional space. Since the inverse of the covariance matrix cancels out the physical dimension of the position deviation in the calculation, the squared term of the Mahalanobis distance is a dimensionless statistical distance.

[0099] Since the three-dimensional position error corresponds to a three-degree-of-freedom chi-square distribution under the Gaussian approximation condition, the envelope expansion factor threshold... (Also a dimensionless pure number) can be determined according to the chi-square quantile corresponding to the preset confidence level. When using a spatial coverage confidence interval of approximately 95% to 99.9%, The acceptable range is greater than or equal to 7.815 and less than or equal to 16.27; when further expansion of the warning envelope is required in engineering applications, It can also be increased to 25. When it is desired to cover approximately 99% of the three-dimensional position error, The preferred value is 11.345. By using the above value method, a balance can be achieved between risk coverage capability and onboard computing load for the off-orbit crossing envelope.

[0100] After completing the three-dimensional deorbit crossing envelope calculation for the satellite to be deorbited, the system proceeds to step S2, namely the directional screening and data distribution step. This step is used to screen out related satellites from the constellation nodes that may be affected by the deorbiting process of the satellite to be deorbited, and send them the data required for subsequent risk assessment, avoiding excessive communication resource consumption caused by broadcasting across the entire network. Step S2 includes the following sub-steps.

[0101] S201, preliminary screening of spatial overlap of associated satellites based on ephemeris comparison.

[0102] The onboard computing unit of the satellite to be deorbited extracts the constellation's global approximate ephemeris data stored onboard. As a preferred method, the onboard computing unit first reads the orbital plane inclination and right ascension parameters of the ascending node of each node, and eliminates nodes that are in completely non-intersecting orbits with the satellite to be deorbited or that have no possibility of rendezvous within a preset time window, in order to reduce the amount of computation required for subsequent spatial overlap determination.

[0103] For the remaining nodes, the onboard computing unit estimates the space occupancy range of the remaining active satellite nodes within the constellation based on ephemeris data. This space occupancy range can be dynamically determined based on the corresponding node's orbit prediction error, maneuver response time, relative velocity, and mission safety margin. As a simplified implementation, the space occupancy range can be quantified as a physical sphere with the node center as its center and a preset screening radius as its radius. The preset screening radius is only used for initial screening of associated satellites and is not used as the final collision avoidance distance. Its specific value can be configured according to constellation density, ephemeris update cycle, and onboard computing resources.

[0104] If, within the expected deorbiting time window, the spatially occupied region of an active satellite node physically overlaps with the three-dimensional deorbiting crossing envelope of the satellite to be deorbited (i.e., their mathematical intersection is not empty), then the satellite to be deorbited includes this active satellite node in the associated satellite set. Through this method, the system identifies associated satellites whose deorbiting crossing process overlaps with that of the satellite to be deorbited.

[0105] S202, Standardized encapsulation and targeted transmission of off-track intention data frames.

[0106] The onboard computing unit of the satellite to be deorbited combines and packages the node physical identifier code, the timestamp sequence corresponding to the system's global time variables, the time evolution sequence of the predicted position vector of the satellite to be deorbited, and the time evolution sequence of the position error covariance matrix of the satellite to be deorbited, according to the preset communication protocol format within the constellation, to form a deorbit intention data frame.

[0107] The satellite to be deorbited controls its own inter-satellite communication transceiver unit to transmit the deorbiting intention data frame to each associated satellite in the associated satellite set via physical radio channels or optical inter-satellite links. After receiving the deorbiting intention data frame, each associated satellite performs subsequent rendezvous risk assessment and exclusion corridor parameter generation based on the timestamp sequence, predicted position vector, and position error covariance matrix contained therein.

[0108] After the deorbiting satellite completes the screening of associated satellites and sends the deorbiting intention data frame, the system executes step S3, namely the risk assessment and corridor generation step. This step is used to enable the associated satellite to calculate the rendezvous risk based on the received deorbiting intention data frame and to transform its own safety constraints into two-layer four-dimensional spatiotemporal exclusion corridor parameters that can be used by the satellite to be deorbited for subsequent optimization. This step includes the following sub-steps.

[0109] S301, Time Delay Compensation and Time Parameterization Nearest Approach Analysis.

[0110] In this embodiment, the associated satellite receives deorbit intention data frames through its own inter-satellite communication transceiver unit. The inter-satellite communication transceiver unit can be a Ka-band phased array antenna, a laser communication terminal, or other inter-satellite communication payloads suitable for low-Earth orbit constellations.

[0111] Due to the inter-satellite distance between the satellite to be deorbited and its associated satellites, the deorbiting intention data frame experiences a propagation delay during transmission. To mitigate the impact of time misalignment on the rendezvous risk assessment results, the onboard computing unit of the associated satellite extracts the timestamp sequence from the deorbiting intention data frame and compares it with its own system's local clock to obtain the communication delay difference. The communication delay difference can be determined by dividing the current physical distance between the two satellites by the speed of light in a vacuum, or it can be estimated by combining the inter-satellite link ranging results with the time synchronization frame.

[0112] Correlated satellites compensate for the time alignment of their own ephemeris data based on the communication delay difference, and calculate their own position status under the synchronization time to reduce position assessment errors caused by inconsistent time references under high-speed relative motion conditions.

[0113] After obtaining synchronization, the associated satellites perform a time-parameterized relative trajectory search within a preset rendezvous time window, calculating the closest approach time, closest approach distance, and relative state covariance between the satellite to be deorbited and the associated satellites. The closest approach time is calculated using the formula:

[0114] ;

[0115] In the formula, Indicates the closest point in time; The operator that represents the independent variable that minimizes the objective function; Represents a system-wide time variable; Indicates the start time of the risk assessment time window; Indicates the length of the risk assessment time window; Indicates the time window for risk assessment; This represents the predicted position vector of the satellite to be deorbited; Indicates the first The predicted position vectors of the associated satellites; Indicates the associated satellite sequence number index; Indicates the satellite to be deorbited and the first Predicted relative position vectors between the associated satellites; The Euclidean norm operator is used to represent the Euclidean norm. Indicates the satellite to be deorbited and the first The predicted relative distances between the associated satellites.

[0116] Correlated satellites calculate rendezvous risk based on their relative positions, relative velocities, and joint position error covariance matrix near the closest approach time. For approximate Keplerian orbit scenarios, correlated satellites can assist in preliminary geometric screening using the minimum orbital intersection distance algorithm; however, in this embodiment, the final rendezvous risk assessment is based on the closest approach analysis results after time synchronization.

[0117] S302, Positive definiteness of state covariance matrix and numerical stabilization processing.

[0118] During the rendezvous risk assessment process, the associated satellites need to superimpose the position error covariance of the satellite to be deorbited with their own position error covariance to generate a joint position error covariance matrix. Considering that the joint position error covariance matrix may lose its positive definiteness due to floating-point rounding errors, nonlinear truncation errors, and numerical propagation errors during multi-orbit sub-long-period extrapolation, the onboard computing unit of the associated satellite performs regularization processing on the joint position error covariance matrix to avoid numerical divergence in the subsequent matrix inversion process.

[0119] As a preferred approach, this regularization is achieved by symmetrizing the joint position error covariance matrix and applying eigenvalue lower bound constraints. Let the joint position error covariance matrix before positive definiteness be... The associated satellites are first symmetricized using the joint position error covariance matrix symmetry formula, which is as follows:

[0120] ;

[0121] In the formula, This represents the joint position error covariance matrix after symmetry processing; This represents the joint position error covariance matrix before positive definiteness. This represents the transpose of the joint position error covariance matrix before positive definite transformation; This represents the matrix transpose operator. Represents the divisor in the symmetric averaging operation; This represents the symmetry treatment term used to average the joint position error covariance matrix and its transpose before positive definite transformation.

[0122] The associated satellites are then used to construct a positive definite matrix using the joint position error covariance matrix positive definiteness formula. The joint position error covariance matrix positive definiteness formula is as follows:

[0123] ;

[0124] In the formula, This represents the joint position error covariance matrix after positive definite processing; This represents the joint position error covariance matrix after symmetry processing; Represents the regularization coefficient; Represents a third-order identity matrix; This represents a regularized diagonal matrix used to improve the minimum eigenvalue of the joint position error covariance matrix.

[0125] Regularization coefficient The regularization coefficient is determined using the following formula:

[0126] ;

[0127] In the formula, Represents the regularization coefficient; This represents the function for selecting the maximum value. This represents the lower bound of the regularization coefficient; Indicates the lower limit of the preset feature value; Represents the minimum eigenvalue function; This represents the joint position error covariance matrix after symmetry processing; This represents the smallest eigenvalue of the joint position error covariance matrix after symmetry processing; This represents the difference between the preset lower limit of eigenvalues ​​and the minimum eigenvalue of the joint position error covariance matrix after symmetry processing.

[0128] After the above processing The minimum eigenvalue is not less than This ensures that the subsequent boundary matrix construction and matrix inversion operations have a stable numerical basis.

[0129] S303, Mapping and Construction of Inner Layer Rigid Safety Boundary Parameters.

[0130] By combining the joint position error covariance matrix after positive definite processing and the physical envelope size of the associated satellites, the onboard computing unit of the associated satellites establishes an inner space no-entry boundary for basic collision avoidance. This boundary is used to define a rigid safety zone that satellites to be deorbited are prohibited from entering.

[0131] Geometrically, the inner rigid safety boundary can be considered as a safety boundary obtained by approximating the probabilistic error covariance ellipsoid with the satellite entity's circumscribed envelope using the Minkowski summation. The boundary parameters are determined using the inner rigid safety boundary shape matrix calculation formula, which is:

[0132] ;

[0133] In the formula, Represents the shape matrix of the inner rigid safety boundary; Indicates the inner layer's safety expansion coefficient; This represents the joint position error covariance matrix after positive definite processing; This represents the joint position error covariance matrix term after being amplified according to the inner layer safety expansion coefficient; Indicates the satellite's physical collision avoidance equivalent radius; The term represents the square of the satellite's physical collision avoidance equivalent radius; Represents a third-order identity matrix; This represents the isotropic physical envelope matrix term formed by the equivalent radius of the satellite's physical collision avoidance. This represents the matrix addition operator.

[0134] Inner layer safety expansion coefficient The scale used to control the confidence interval of the basic collision avoidance probability. When the joint position error approximately follows a three-dimensional Gaussian distribution, It can be set according to the three-degree-of-freedom chi-square quantiles. As a preferred option, The value range is set to be greater than or equal to 7.815 and less than or equal to 11.345, respectively, to correspond to approximately 95% to 99% of the three-dimensional spatial coverage confidence interval. In the low-risk screening stage or when the onboard computing power is limited, a lower engineering screening value can be used, but the confidence threshold that meets the collision avoidance safety requirements should be used in the final action timing confirmation stage.

[0135] Satellite physical collision avoidance equivalent radius The value is determined by the maximum circumscribed circle radius after the satellite's three-dimensional structure is unfolded, combined with the attitude control pointing error margin, and is generally between 5m and 25m.

[0136] S304, Outer elastic service boundary parameter mapping and double-nested corridor parameter generation.

[0137] In addition to ensuring physical collision avoidance safety, when performing Earth observation, communication, or other on-orbit missions, the associated satellites also need to reduce the risks of optical obstruction, radio beam overlap, or mission service interruption caused by close flybys of satellites awaiting deorbiting. Therefore, the onboard computing unit of the associated satellites further constructs an outer elastic service boundary outside the inner rigid safety boundary.

[0138] The parameters are generated using the outer elastic service boundary shape matrix calculation formula, which is as follows:

[0139] ;

[0140] In the formula, Represents the shape matrix of the outer elastic service boundary; Indicates the scaling factor for the task protection zone; Represents the shape matrix of the inner rigid safety boundary; This represents the shape matrix of the inner rigid safety boundary after being magnified according to the task protection zone scaling factor; Indicates the margin for avoiding task interference; The term representing the squared margin of task interference avoidance; Represents a third-order identity matrix; This represents an isotropic task service protection matrix term formed by the task interference avoidance margin.

[0141] Task Protection Zone Scaling Factor This is used to extend the service protection area outward from the collision avoidance baseline, and its value ranges from 1.2 to 3.0 depending on the current task sensitivity. Task interference avoidance margin This represents the physical spacing constant required to avoid optical obstruction or radio beam overlap. Its value can be set according to the characteristics of the onboard payload, link beamwidth, observation field of view, and constellation density. Preferably, this represents the mission interference avoidance margin. The range can be set from 5km to 30km; different mission interference avoidance margins can be set for communication payloads, remote sensing payloads, or highly sensitive scientific payloads.

[0142] As can be seen from the above matrix construction relationship, the outer elastic service boundary covers the inner rigid safety boundary in space and serves as a service protection constraint that can be controlled to shrink in subsequent optimization solutions.

[0143] After completing the above calculations, the onboard computing unit of the associated satellite combines and encapsulates the inner rigid safety boundary shape matrix, the outer elastic service boundary shape matrix, and its own position state sequence to generate a two-layer four-dimensional spatiotemporal exclusion corridor parameter containing the inner rigid safety boundary and the outer elastic service boundary, and feeds it back to the satellite to be deorbited through the inter-satellite communication transceiver unit.

[0144] In the optimization solution process, the satellite to be deorbited is described using the inner rigid safety boundary constraint formula. The inner rigid safety boundary corresponding to each associated satellite, and the constraint formula for the inner rigid safety boundary are as follows:

[0145] ;

[0146] In the formula, Indicates the global time variable of the satellite to be deorbited. The predicted position vector at the corresponding time; Indicates the first The correlation of satellites in the system's global time variable The predicted location vector at the corresponding time (which uses the same parameter as the predicted location vector in the aforementioned risk assessment stage). Indicates the associated satellite sequence number index; Indicates the satellite to be deorbited and the first The relative position vector between the associated satellites; Represents the transpose of this relative position vector; Indicates the first The shape matrix of the inner rigid safety boundary corresponding to each associated satellite; Indicates the first The inverse matrix of the inner rigid safety boundary shape matrix corresponding to each associated satellite; Indicates the first The normalized squared distance term of the inner rigid safety boundary corresponding to each associated satellite; This indicates the threshold for boundary determination.

[0147] Among them, the inner rigid safety boundary serves as a physical collision avoidance constraint that cannot be relaxed, while the outer elastic service boundary serves as a task service constraint that can be controlled and contracted within a preset range.

[0148] After the satellite to be deorbited receives the two-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data fed back by the associated satellites, the system executes step S4, the global cost construction step. This step is used to map the maneuver cost of the satellite to be deorbited, the secondary mission cost of the associated satellites, and the relaxation penalty of the outer elastic service boundary to the same optimization objective, so as to form the evaluation benchmark for subsequent collaborative planning. Step S4 includes the following sub-steps.

[0149] S401, Quantification of heterogeneous state value of nodes and allocation of weight coefficients.

[0150] The inter-satellite communication transceiver unit of the satellite to be deorbited receives the two-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data. The node cost data includes at least the remaining design lifetime percentage of the associated satellites, the current mission priority, data communication load, maneuverability status, and acceptable service interruption margin.

[0151] Because satellites awaiting deorbiting and associated satellites currently performing on-orbit missions differ in terms of system retention value, mission priority, and maneuver resources, the onboard computing unit of the satellite awaiting deorbiting assigns heterogeneous cost weights to the satellite and associated satellites based on the remaining design life percentage of each node, the current data communication traffic, mission priority, and maneuverability status.

[0152] To facilitate onboard optimization, the onboard computing unit normalizes each cost weight to a dimensionless range of 0 to 1. In this embodiment, the relinquishment weight of the satellite to be deorbited is mapped to a low weight range of 0.01 to 0.1, enabling it to assume primary responsibility for trajectory adjustment during collaborative avoidance; the protection weight of each associated satellite is mapped to a high weight range of 0.5 to 1.0 based on mission priority, in order to reduce the probability of its normal on-orbit service being interrupted.

[0153] S402, Calculation of the cost of secondary tasks in the system.

[0154] When assessing the global impact of joint evasion maneuvers, the onboard computing unit of the satellite to be deorbited calculates the resource consumption and mission disturbance of each node. The costs of secondary missions include, but are not limited to, the propellant mass cost consumed by the orbital control actuators in performing maneuvers, the cost of the ground coverage gap caused by deviation from the rated operating orbit, the cost of communication link interruption, and the cost of attitude recovery.

[0155] Because the aforementioned costs correspond to different physical dimensions, the onboard computing unit normalizes each physical consumption quantity, converting it into a dimensionless relative cost score between 0 and 1. Normalization methods can include maximum value normalization, nominal task threshold normalization, or segmented normalization based on task priority; specific normalization rules can be configured according to the constellation's mission type and onboard computing resources. The onboard computing unit multiplies the normalized relative cost score by the corresponding cost weight to obtain the weighted cost value of each node.

[0156] S403 integrates the predefined dynamic constraint relaxation penalty term with the constellation-level global deorbit cost function.

[0157] To address the solution stagnation caused by hard constraints on the outer elastic service boundary in congested space scenarios, the onboard computing unit of the satellite to be deorbited pre-includes a dynamic constraint relaxation penalty term in the global cost function. This penalty term is used in the subsequent elastic degradation and replanning phase to transform the controlled shrinkage of the outer elastic service boundary into a numerical cost, which participates in the optimization solution together with the maneuver cost and mission disturbance cost.

[0158] The onboard computing unit integrates the weighted maneuver cost of the satellite to be deorbited, the weighted secondary generation cost of the associated satellite set, and the dynamic constraint relaxation penalty term, and generates a constellation-level global deorbit cost function using the constellation-level global deorbit cost function calculation formula. The constellation-level global deorbit cost function calculation formula is as follows:

[0159] ;

[0160] In the formula, This represents the constellation-level global deorbit cost function. This represents the weighting coefficient for the maneuver of the satellite to be deorbited; This indicates the value of the maneuvering of the satellite to be deorbited; This represents the thrust control variables for the satellite to be deorbited. This represents the weighted maneuver value of the satellite to be deorbited. This represents the summation operator for a set of associated satellites. Indicates the total number of associated satellites; Indicates the first The protection weight coefficient for each associated satellite; Indicates the first The secondary mission disturbance cost of each associated satellite (used to quantify the resource consumption or normal service interruption loss caused by the associated satellite due to evasion maneuvers). Indicates the first The weighted secondary mission disturbance value of each associated satellite; This represents the weighted sum of secondary mission perturbation costs for the associated satellite set; Indicates the dynamic penalty baseline coefficient; Indicates the first The mission sensitivity penalty weight for each associated satellite; Indicates the first The constant slack variables used by each associated satellite within the replanning time window; Indicates the first The squared terms of the constant slack variables of each associated satellite; Indicates the first A mission sensitivity-weighted slack penalty term for each associated satellite; This represents the penalty term for slack in dynamic constraints.

[0161] Dynamic penalty benchmark coefficient This is used to control the penalty intensity introduced into the objective function when the outer flexible service boundary is contracted in a controlled manner. During the basic planning phase, the system relaxes the constant variables corresponding to each associated satellite. The outer flexible service boundary is used as a hard constraint for solving; when the optimizer continuously outputs a state code indicating no feasible solution, the system enters the elastic degradation and replanning phase. Set to greater than initial penalty value .

[0162] No. The slack variable for each associated satellite represents the geometric distance by which the outer flexible service boundary of that satellite is allowed to shrink. When using a constant-value form within the reprogramming time window, the slack variable for the first satellite... The constant slack variables corresponding to the associated satellites are denoted as . The formula for dynamic slack variable constraints is:

[0163] ;

[0164] In the formula, This represents the lower limit of the slack variable; Indicates the first The constant slack variable used by each associated satellite within the replanning time window is used to characterize the geometric distance by which the outer flexible service boundary of the associated satellite is allowed to shrink within the current replanning time window. Indicates the first The maximum allowable shrinkage margin for each associated satellite within the replanning time window.

[0165] By constructing the cost function as described above, satellites to be deorbited can simultaneously evaluate their own maneuvering costs, related satellite mission disturbances, and outer service boundary relaxation costs within the same optimization objective, providing a numerical solution basis for subsequent flexible degradation and replanning.

[0166] After constructing the constellation-level global deorbit cost function with dynamic constraint relaxation penalty terms, the system executes step S5, the elastic degradation replanning step. This step fully defines the logical closed loop from baseline state assessment to anomaly degradation handling. Specifically, in this step, the system first uses the outer elastic service boundary as a hard constraint for baseline solution. If no feasible solution is found, a replanning mechanism is triggered to controllably soften the outer service protection constraints while maintaining the inner rigid safety boundary from being breached. This step includes the following sub-steps.

[0167] S501, Setting the Inequality Hard Constraint Solution Based on the Outer Elastic Service Boundary. In this embodiment, the onboard computing unit of the satellite to be deorbited transforms the outer elastic service boundary in the received two-layer four-dimensional spatiotemporal exclusion corridor parameters into a space avoidance matrix constraint condition. This transformation process involves constructing a normalized squared distance term of the satellite to be deorbited relative to the outer elastic service boundary and linearizing it when necessary, making it directly callable by the constraint optimization solver. In the initial planning stage, to verify whether an ideal conflict-free path exists in the current rendezvous scenario, the onboard computing unit sets the outer elastic service boundary as a conventional inequality hard constraint. This hard constraint setting serves as a prerequisite for triggering elastic degradation. The system first attempts to minimize the constellation-level global deorbit cost function while simultaneously satisfying the inner rigid safety boundary constraint and the outer elastic service boundary constraint. To obtain the thrust ignition timing or attitude drag control timing, the onboard computing unit calls the constraint optimization solver. The constraint optimization solver can employ sequential quadratic programming, sequential convex optimization, interior point method, or other nonlinear programming solution methods suitable for real-time operation by the onboard processor. When onboard computing resources are limited, the system can reduce the computational load by shortening the prediction time window, reducing the dimensionality of control variables, using piecewise constant thrust control variables, calling a pre-set initial solution library, or having the master control node solve the problem centrally.

[0168] S502, continuous infeasible solution state detection and core bottleneck node targeted localization. In scenarios with dense intersections of space nodes, the outer elastic service boundaries corresponding to different associated satellites may overlap in time and space, causing the optimizer to fail to obtain feasible trajectories that satisfy all hard constraints within the set iteration range. When the onboard computing unit detects that the optimizer is in a continuous state without feasible solutions, it will detect the infeasible solution state without feasible solutions and target localization of the core bottleneck nodes. In each of the subsequent restarts, a "no feasible solution" status code is output, or the maximum number of iterations is reached. If the internal mechanism fails to ensure that all outer elastic service boundary constraints meet the preset clearance threshold, an elastic degradation and replanning mechanism is triggered. and Based on the onboard processor's computing power, rendezvous time window length, and action sequence preparation time, the onboard computing unit, as a preferred method, extracts the Lagrange multiplier values ​​of each constraint in the final iteration to locate the main constraints causing the solution stall. This is then combined with the constraint violation amount and constraint scale normalization factor to calculate the stall score. The sluggishness score of each associated satellite is determined using the sluggishness score calculation formula, which is as follows:

[0169] ;

[0170] In the formula, Indicates the first Obstruction score of each associated satellite; Indicates the associated satellite sequence number index; Indicates the first Lagrange multipliers for the constraints corresponding to each associated satellite; Indicates the first The absolute value of the Lagrange multipliers corresponding to the constraints of each associated satellite; Indicates the first The amount of violation of constraints or insufficient airspace corresponding to each associated satellite; Indicates the first The product of the resistance strength of each associated satellite's corresponding constraint conditions; Indicates the first The scale normalization factor for the constraints corresponding to each associated satellite; This represents a small positive number used to prevent the denominator of the score calculation formula from being zero; Indicates the first The system will adjust the normalized denominator terms for the constraints corresponding to each associated satellite. The most correlated satellite is identified as the core blocking node. If multiple correlated satellites have the same or similar blocking scores, the system can further sort them according to mission priority, nearest proximity distance, or duration of constraint violation.

[0171] S503, Calculation of controlled boundary contraction under baseline collision avoidance constraints. When the satellite to be deorbited encounters a continuous infeasibility problem in the optimization solution, it requests controlled contraction of the outer elastic service boundary corresponding to the core blocking node. Controlled contraction only applies to the outer elastic service boundary; the inner rigid safety boundary remains a non-relaxable constraint. For the [missing information]... For each associated satellite, the system uses a relative position vector calculation formula to determine the relative position vector of the satellite to be deorbited relative to the associated satellite. The relative position vector calculation formula is as follows:

[0172] ;

[0173] In the formula, Indicates the first The relative position vectors of the associated satellites; Indicates the global time variable of the satellite to be deorbited. Position vector at the corresponding time; Indicates the first The correlation of satellites in the system's global time variable Position vector at the corresponding time; Indicates the satellite to be deorbited and the first The relative position vectors between the associated satellites. The corresponding relative distances are calculated using the relative distance calculation formula, which is:

[0174] ;

[0175] In the formula, Indicates the first The relative distances between the associated satellites; Indicates the first The relative position vectors of the associated satellites; The Euclidean norm operator is used to represent the Euclidean norm. Indicates the first The Euclidean norm of the relative position vectors corresponding to each associated satellite. At that time, the system uses the relative azimuth unit vector calculation formula to determine the current relative azimuth unit vector. The relative azimuth unit vector calculation formula is as follows:

[0176] ;

[0177] In the formula, Indicates the first The relative azimuth unit vector corresponding to each associated satellite; Indicates the first The relative position vectors of the associated satellites; Indicates the first The relative distances between the associated satellites; Indicates the first The unit vector is obtained by normalizing the relative position vectors corresponding to each associated satellite. The system uses the formula for calculating the inner layer equivalent boundary radius to determine the inner layer equivalent boundary radius along the relative azimuth. The formula for calculating the inner layer equivalent boundary radius is:

[0178] ;

[0179] In the formula, Indicates the first The equivalent inner boundary radius corresponding to each associated satellite; Indicates the first The relative azimuth unit vector corresponding to each associated satellite; Indicates the first The transpose of the relative azimuth unit vectors corresponding to each associated satellite; Indicates the first The shape matrix of the inner rigid safety boundary corresponding to each associated satellite; Indicates the first The inverse matrix of the inner rigid safety boundary shape matrix corresponding to each associated satellite; Indicates the first The inverse quadratic term of the inner rigid safety boundary corresponding to each associated satellite along the relative azimuth; It represents the exponent of the inverse quadratic term when the negative one-half power is applied; This indicates the conversion of the inverse quadratic term of the inner rigid safety boundary along the relative orientation into a calculation term with the dimension of radius (length). The system uses the formula for calculating the outer equivalent boundary radius to determine the outer equivalent boundary radius along the relative orientation. The formula for calculating the outer equivalent boundary radius is:

[0180] ;

[0181] In the formula, Indicates the first The equivalent outer boundary radius corresponding to each associated satellite; Indicates the first The relative azimuth unit vector corresponding to each associated satellite; Indicates the first The transpose of the relative azimuth unit vectors corresponding to each associated satellite; Indicates the first The outer elastic service boundary shape matrix corresponding to each associated satellite; Indicates the first The inverse matrix of the outer elastic service boundary shape matrix corresponding to each associated satellite; Indicates the first The inverse quadratic term of the outer elastic service boundary corresponding to each associated satellite along the relative azimuth; It represents the exponent of the inverse quadratic term when the negative one-half power is applied; This indicates the conversion of the inverse quadratic term of the outer elastic service boundary along the relative orientation into a calculation term with the dimension of length (radius). The system uses the formula for calculating the signed clearance distance of the outer boundary to determine the signed clearance distance of the outer boundary. The formula for calculating the signed clearance distance of the outer boundary is:

[0182] ;

[0183] In the formula, Indicates the first The signed clearance distance of the outer boundary corresponding to each associated satellite; Indicates the first The relative distances between the associated satellites; Indicates the first The equivalent outer boundary radius corresponding to each associated satellite; Indicates the first The difference between the relative distances of the associated satellites and the radius of the outer equivalent boundary. The system uses the inner boundary signed clearance distance calculation formula to determine the inner boundary signed clearance distance. The inner boundary signed clearance distance calculation formula is as follows:

[0184] ;

[0185] In the formula, Indicates the first The signed clearance distance of the inner boundary corresponding to each associated satellite; Indicates the first The relative distances between the associated satellites; Indicates the first The equivalent inner boundary radius corresponding to each associated satellite; Indicates the first The difference between the relative distances of the associated satellites and the radius of the inner equivalent boundary. During the elastic degradation phase, the controlled contraction of the outer elastic service boundary is described using the outer boundary relaxation constraint formula. With constant relaxation variables within the replanning time window, the outer boundary relaxation constraint formula is:

[0186] ;

[0187] In the formula, Indicates the first The signed clearance distance of the outer boundary corresponding to each associated satellite; Indicates the first The constant slack variables used by each associated satellite within the replanning time window; Represents a global time variable in the system. At the corresponding time, for the first The net distance between the outer boundaries of the associated satellites after applying constant relaxation variables; This represents the lower limit of the clearance distance after the outer boundary relaxes. During the elastic degradation phase, the inner boundary preservation constraint formula is simultaneously used to ensure that the inner rigid safety boundary cannot relax. The inner boundary preservation constraint formula is:

[0188] ;

[0189] In the formula, Indicates the first The signed clearance distance of the inner boundary corresponding to each associated satellite; This represents the lower limit of the signed clearance distance at the inner boundary. The system calculates the corresponding clearance shortfall based on the signed clearance distance at the outer boundary, and determines it using the outer boundary clearance shortfall calculation formula. The formula for calculating insufficient clearance at the outer boundary is:

[0190] ;

[0191] In the formula, Indicates the first The correlation of satellites in the system's global time variable The outer boundary clearance is insufficient at the corresponding moment; This represents the function for selecting the maximum value. This indicates the lower limit of insufficient clearance at the outer boundary; Indicates the first The outer boundary of each associated satellite is the negative of the signed clearance distance. When When it is negative, This represents the minimum relaxation distance required to re-satisfy the outer boundary relaxation constraints; when When it is a non-negative value, Pick If constant slack variables are used within a replanning time window, the system uses the constant slack calculation formula to calculate the first... The constant slack variable for each associated satellite, the formula for calculating the constant slack is:

[0192] ;

[0193] in, The maximum allowable shrinkage margin of the replanning time window is determined using the formula:

[0194] ;

[0195] In the formula, Indicates the first The constant slack variables used by each associated satellite within the replanning time window; This represents the function for selecting the minimum value. This represents the function for selecting the maximum value. Represents a system-wide time variable; Indicates the start time of the risk assessment time window; Indicates the length of the risk assessment time window; Indicates the time window for risk assessment; Indicates the first The correlation of satellites in the system's global time variable The outer boundary clearance is insufficient at the corresponding moment; Indicates the first The maximum outer boundary clearance of a single associated satellite within the risk assessment time window; Indicates the first The maximum allowable shrinkage margin for each associated satellite within the replanning time window; Indicates the first The correlation of satellites in the system's global time variable The maximum allowable contraction margin at the corresponding time point; Indicates the first The minimum value among the maximum permissible shrinkage margins of each associated satellite at each moment within the risk assessment time window. If time-by-time relaxation is used, the system calculates the dynamic relaxation variable using the time-by-time relaxation calculation formula, and then... The correlation of satellites in the system's global time variable The time-by-time dynamic slack variable corresponding to the given time is denoted as The formula for calculating the relaxation amount at each time step is:

[0196] ;

[0197] In the formula, Indicates the first The correlation of satellites in the system's global time variable The time-by-time dynamic slack variables at the corresponding time point are distinct from the constant slack variables that take constant values ​​within the replanning time window. ; This represents the function for selecting the minimum value. Indicates the first The correlation of satellites in the system's global time variable The outer boundary clearance is insufficient at the corresponding moment; Indicates the first The correlation of satellites in the system's global time variable The maximum allowable shrinkage margin at the corresponding time point. In this embodiment, to reduce the onboard computational load, a constant slack variable with a constant value within the replanning time window is preferably used. Instead of using time-varying dynamic slack variables. . No. The correlation of satellites in the system's global time variable The maximum permissible shrinkage margin at a given time is determined using the maximum permissible shrinkage margin calculation formula, which is:

[0198] ;

[0199] In the formula, Indicates the first The correlation of satellites in the system's global time variable The maximum allowable contraction margin at the corresponding time point; Indicates the first The equivalent outer boundary radius corresponding to each associated satellite; Indicates the first The equivalent inner boundary radius corresponding to each associated satellite; Indicates the first The difference between the outer equivalent boundary radius and the inner equivalent boundary radius corresponding to each associated satellite. This setting ensures that even if the outer elastic service boundary undergoes controlled contraction, it will not cross the inner rigid safety boundary, thus maintaining the physical collision avoidance baseline.

[0200] S504, numerical penalty dynamic injection and cooperative compromise solution optimization convergence. In obtaining the... After slackling the constant variables of the associated satellites, the onboard computing unit of the satellite to be deorbited updates the dynamic penalty baseline coefficients used to control the penalty intensity. As a preferred approach, the onboard computing unit will dynamically penalize the baseline coefficient. From the initial penalty value Begin, according to the preset ratio The iterative increment is represented by the following increment logic:

[0201] ;

[0202] In the formula, Indicates the first The dynamic penalty benchmark coefficient corresponding to the next elastic solution iteration; Indicates the first The dynamic penalty benchmark coefficient corresponding to the next elastic solution iteration; This indicates the number of internal iterations in the elastic solution; This indicates the increment rate of the dynamic penalty benchmark coefficient; This represents the multiplication operator; This represents the dynamic penalty baseline coefficient updated according to an incremental ratio. As a preferred option, A scalar value between 1.2 and 2.0 can be used. The onboard computation unit transforms the contraction corresponding to the constant relaxation variable into a numerical penalty and injects it into the constellation-level global deorbit cost function with a dynamic constraint relaxation penalty term. Subsequently, the onboard computation unit restarts the constraint optimization solution process based on the updated soft constraints. With the controlled relaxation of the outer elastic service boundary constraints and the update of the dynamic constraint relaxation penalty term, the optimizer determines the cooperative compromise solution that reduces the global comprehensive cost from the feasible solution set that satisfies the inner rigid safety boundary constraints, action execution constraints, and termination conditions. If the dynamic penalty baseline coefficient... If the optimizer still has no feasible solution output after the penalty is increased to the set upper limit threshold, the system determines that there is no executable cooperative evasion path under the current physical space and control resource conditions, and switches to the backup trajectory decay process.

[0203] After the multi-satellite collaborative planning is completed, the system executes step S6, namely the scheme confirmation and backup execution step. In this step, the satellite to be deorbited converts the collaborative compromise solution into corresponding action timing sequences and sends them to the associated satellites so that the corresponding orbit control actuators can perform physical maneuvers. If communication is abnormal, the action timing sequences cannot be confirmed in a closed loop, or the actuators of the satellite to be deorbited fail, the system terminates the collaborative execution process, cancels the action timing sequences that have not yet entered the irreversible execution window, and then switches to the backup orbit decay process. This step includes the following sub-steps.

[0204] S601, Cooperative compromise solution transformation and node action timing confirmation.

[0205] The onboard computing unit of the satellite to be deorbited will analyze the obtained cooperative compromise solution into discrete control commands and generate the corresponding action timing. For satellite nodes using active orbit control, this analysis process can be based on the equivalent impulse principle and pulse width modulation logic to discretize the continuous thrust curve output by the optimizer into ignition pulse commands that include thruster start-up time, start-up time, and duty cycle parameters.

[0206] For satellite nodes employing differential drag control, this analytical process can be translated into aerodynamic drag control commands that include attitude pointing angle, frontal area configuration, and attitude holding duration.

[0207] The satellite awaiting deorbiting transmits its action sequence to the associated satellites participating in the avoidance maneuver through its own inter-satellite communication transceiver unit, based on the cooperative compromise solution. Under normal operating conditions, each satellite node activates its underlying hardware according to the confirmed action sequence. Specifically, the attitude control actuator adjusts the satellite's attitude according to the action sequence to establish a predetermined thrust direction; the orbit control actuator performs active orbital maneuvers according to the thrust direction. For differential drag control nodes, the attitude control actuator adjusts the satellite's frontal area to perform aerodynamic drag modulation.

[0208] S602, Abnormal operating condition communication blockage judgment and termination of coordinated ignition sequence.

[0209] The system configures a periodic health status heartbeat packet detection mechanism for on-orbit nodes to set a disability condition baseline. In this embodiment, the heartbeat packet timeout threshold is the consecutive loss of acknowledgment data frames for three calibration communication cycles.

[0210] If, during the execution of the predetermined action sequence, the inter-satellite communication transceiver unit of the satellite to be deorbited experiences a communication interruption due to electromagnetic interference or other reasons, triggering the aforementioned timeout threshold, or if the internal sensors detect a failure in the orbit control actuator such as fuel leakage, insufficient thrust, valve abnormality, or attitude instability, the system monitoring log will generate an abnormality flag.

[0211] Upon detecting the aforementioned disabling conditions, and confirming that the deorbited satellite and its associated satellites can no longer maintain synchronization in a spatiotemporal state, the system terminates the multi-satellite collaborative planning mechanism. For ignition sequences that have not yet entered the irreversible execution window, the system cuts off or cancels the corresponding action sequence; for action sequences that have already entered the irreversible execution window, each node executes according to pre-agreed safety shutdown rules, which include completing the minimum safety pulse, entering attitude maintenance mode, broadcasting anomaly status indicators, and switching to the local avoidance calculation process.

[0212] S603, passive aerodynamic drag mode switching and baseline attenuation process execution.

[0213] When the active orbit control actuator of the satellite to be deorbited fails, inter-satellite communication is interrupted, or the timing of coordinated actions cannot be confirmed in a closed loop, the system terminates the multi-satellite coordinated planning mechanism and enters the minimum orbit decay process. Under the minimum orbit decay process, the system degrades to a one-way passive avoidance mode.

[0214] For microsatellite nodes with differential drag control capabilities, the orbit control actuator and attitude control actuator can be functionally merged, generating incremental aerodynamic drag by changing the frontal area. The control principle no longer relies on the working propellant to generate active thrust, but instead adjusts the satellite's aerodynamic shape to change the orbital decay rate.

[0215] In this embodiment, the guaranteed orbit decay process determines the execution sequence based on the residual capability of the actuators. If the satellite to be deorbited still possesses attitude control capabilities, it adjusts the orientation of the solar panels or drag-increasing components via the attitude control flywheel to increase the equivalent windward area along the velocity direction, creating a continuous increase in aerodynamic drag to reduce orbital altitude. If the attitude control actuator fails, but the satellite to be deorbited is equipped with passive deorbiting devices such as drag-increasing sails or deployable drag plates with independent or pre-set triggering capabilities, the system triggers these passive deorbiting devices when the release conditions are met.

[0216] If the current orbital altitude and space environment conditions prevent attitude drag enhancement or passive deorbiting devices from meeting the preset deorbiting deadline, the system will prioritize perigee reduction maneuvers before active control capabilities completely fail. If active control capabilities are insufficient to complete perigee reduction maneuvers, the satellite to be deorbited will be marked as a passively decaying target requiring long-term unilateral avoidance by associated satellites. During this period, associated satellites will independently calculate local avoidance schemes based on the passive decay trajectory of the satellite to be deorbited under atmospheric drag, ensuring the operational safety of the constellation's physical network architecture through unilateral concession maneuvers.

[0217] It should be noted that the low-Earth orbit constellation multi-satellite cooperative deorbiting method of this invention does not presuppose the existence of feasible cooperative solutions in all space rendezvous scenarios. When the inner rigid safety boundaries inevitably overlap in the physical space and time dimensions, or when the satellite to be deorbited completely loses its attitude control, orbit control, and communication capabilities, the system determines that the current cooperative plan is not executable and marks the satellite to be deorbited as a passive risk target. The associated satellites and the ground space situational awareness system then perform continuous tracking and unilateral avoidance.

[0218] To aid in understanding the physical transformation logic and calculation process of this invention, a specific application example is given below, using simulated orbital operation data in a space environment, along with numerical substitution calculations.

[0219] At the current moment corresponding to the system's global time variable, it is assumed that the diagonal elements of the position error covariance matrix of the satellite to be deorbited in the Earth's inertial coordinate system are all 0.01, with units of 0.01. The system-set envelope expansion factor threshold. A value of 9 is used to construct the off-track crossing envelope during the rapid screening phase of the project; when entering the final action timing confirmation phase, the system can... The value was increased to 11.345, corresponding to a 99% three-dimensional confidence level, to enhance the conservatism of the final collision avoidance confirmation.

[0220] An envelope model is generated using the three-dimensional off-track crossing envelope calculation formula. Substituting the numerical values ​​into the three-dimensional off-track crossing envelope calculation formula yields:

[0221] ;

[0222] After identifying the first associated satellite that causes physical overlap, the associated satellite performs positive definite processing on the calculated joint position error covariance matrix to obtain a positive definite joint position error covariance matrix with all diagonal elements being 0.02.

[0223] In this embodiment, to illustrate the boundary generation process in the final action timing confirmation stage, an inner layer safety expansion coefficient is set. The value is 7.815, and the equivalent radius for satellite physical collision avoidance is 0.015 km. The boundary parameters are determined using the inner rigid safety boundary shape matrix calculation formula. Substituting the values ​​into the formula yields:

[0224] ;

[0225] Based on the sensitivity of the current on-orbit communication mission, the mission protection zone scaling factor is set to 2.0, and the mission interference avoidance margin is set to 5.0 km. Parameters are generated using the outer elastic service boundary shape matrix calculation formula. Substituting the values ​​into the formula yields:

[0226] ;

[0227] When the optimizer sets the outer elastic service boundary as an inequality hard constraint and outputs a no-solution status code, it triggers the elastic degradation replanning mechanism. The optimizer extracts the constraint conflict requirement of the first associated satellite as 3.5 km. Combining the boundary geometry characteristics, the system obtains the maximum allowable shrinkage margin of the first associated satellite as 4.635 km.

[0228] The specific value is determined using the controlled boundary shrinkage calculation formula. Substituting the value into the controlled boundary shrinkage calculation formula yields:

[0229] ;

[0230] In the formula, This represents the constant slack variable used by the first associated satellite within the replanning time window; This represents the function for selecting the minimum value. This represents the numerical value indicating the insufficient outer boundary clearance corresponding to the first associated satellite; This represents the maximum allowable shrinkage margin for the first associated satellite within the replanning time window; This represents the smaller of the value of the insufficient outer boundary clearance corresponding to the first associated satellite and the value of the maximum allowable shrinkage margin of the first associated satellite within the replanning time window.

[0231] Therefore, the system determines the contraction of the outer elastic service boundary of the first associated satellite to be 3.5 km without intruding on the inner rigid safety boundary. This dynamic relaxation variable is then introduced into the numerical penalty module for global optimization, and the system finally obtains a cooperative compromise solution that balances avoidance safety and control costs.

[0232] To verify the engineering application performance of the above method, a multi-satellite collaborative architecture for a low Earth orbit constellation was constructed and tested in a simulation environment. A constellation operation network containing multiple nodes was configured in a low Earth orbit environment at an altitude of 700 km. One satellite node facing thrust performance degradation was extracted as a satellite to be deorbited, and three on-orbit satellites with spatial overlap during the deorbiting phase were selected as associated satellites based on space ephemeris data.

[0233] Two control architectures were set up during the testing process for control variable analysis. The control group adopted a hard constraint solution method with a fixed distance threshold, requiring a fixed envelope spacing of 5.0 km between on-orbit nodes; the experimental group ran the low-Earth orbit constellation multi-satellite coordinated deorbiting method of this embodiment of the invention, triggering elastic degradation replanning when the constraint optimizer outputs a state code of no feasible solution, and injecting a dynamic constraint relaxation penalty term for the outer elastic service boundary into the objective function.

[0234] The computation process runs on a semi-physical simulation platform equipped with a spaceborne real-time operating system control architecture. The system performs orbital nonlinear dynamics integral and constraint optimization solutions, extracts control convergence states, continuous thrust input sequences, and spatial position distribution data of each intersection node during the optimization iteration process, and completes the mapping of data in the three-dimensional physical coordinate system and the two-dimensional parameter solution space.

[0235] Figure 3 The diagram illustrates the geometric relative relationships when a satellite awaiting deorbiting crosses the collision avoidance zone of a associated satellite, including the inner rigid safety boundary, the outer elastic service boundary, the position of the associated satellite, the relative trajectory of the decommissioned satellite, the closest approach point, the closest approach distance, and the outer boundary crossing point. Figure 3 Using the associated satellite as the relative coordinate origin, it reflects the state of the satellite to be deorbited, where the relative crossing trajectory enters the outer elastic service boundary but does not enter the inner rigid safety boundary.

[0236] Figure 4 The diagram illustrates the cost contour lines, the outer soft-constraint feasible boundary, the maximum allowable shrinkage margin, the stagnation region in the hard-constraint solution, the optimization iteration path, and the final cooperative compromise solution. Figure 4 It reflects the morphological distribution of the constellation-level global deorbit cost function with dynamic constraint relaxation penalty term in the space composed of thrust control variables and dynamic relaxation variables, and embodies the process of the system transitioning from a hard-constraint infeasible state to a soft-constraint feasible state.

[0237] Based on simulation results and Figure 3 , Figure 4 The visualization results show that, in a high-density rendezvous scenario, the decommissioning trajectory of the satellite to be deorbited crosses the outer elastic service boundary but does not touch the inner rigid safety boundary. If the hard constraint calculation criterion with a fixed distance threshold in the control group is used, the optimization solver may terminate the calculation due to the degradation of the feasible solution domain caused by the spatial geometric overlap leading to hard constraint conflicts, thus causing the orbit avoidance control process to stall.

[0238] By analyzing the evolution trend of the objective function with the control dimension during the parameter optimization process, the system detects the continuous optimization of the optimizer. If no feasible solution is found after each restart, or after the maximum number of iterations... If the internal service boundary clearance threshold cannot be met, elastic degradation and replanning will be performed. Figure 4 In the solution trajectory, this process manifests as the optimization path shifting from a region of low relaxation variables to a region of non-zero relaxation variables. This displacement indicates that the system actively introduces non-zero dynamic relaxation variables, transforming the originally unsatisfactory outer service protection constraints into numerical penalty terms in the constellation-level global deorbit cost function, while maintaining the inner rigid safety boundary from being breached.

[0239] After several iterations of optimization, the computation point enters a feasible region that satisfies the inner rigid safety boundary constraints and action execution constraints. Convergence is achieved at a position on the cost contour line where the cost is low and the termination condition is met. This allows for the calculation of control parameters that meet the requirements of the cooperative compromise solution while ensuring the physical collision avoidance safety baseline. Simulation and semi-physical simulation results show that this computational architecture expands the feasible region of the optimization solution under the condition of controlled contraction of the outer boundary, reduces the impact of cooperative avoidance actions on the continuity of associated satellite missions, and reduces repeated invalid solutions caused by hard constraint conflicts in crowded orbital environments.

[0240] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for coordinated deorbiting of multiple satellites in a low-Earth orbit constellation, characterized in that, include: The satellite to be deorbited generates initial deorbit candidate schemes, and the deorbit crossing envelope is constructed by extrapolating orbital dynamics. Select the associated satellites that physically overlap with the deorbit crossing envelope, and send the deorbit intention data frame in a directional manner; The associated satellites assess rendezvous risks and generate a two-layer four-dimensional spatiotemporal exclusion corridor parameter that includes an inner rigid safety boundary and an outer elastic service boundary. The satellite to be deorbited receives the double-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data, and generates a constellation-level global deorbit cost function with a dynamic constraint relaxation penalty term. When the outer elastic service boundary has no feasible solution as a hard constraint, controlled contraction is performed only on the outer elastic service boundary. The contraction amount is converted into a numerical penalty and injected into the constellation-level global deorbit cost function. The inner rigid safety boundary is kept from relaxation and a cooperative compromise solution is obtained. The cooperative compromise solution is converted into an action sequence and sent to the associated satellite for execution. In the event of communication failure or inability to confirm the action sequence in a closed loop, the action sequence that has not yet entered the irreversible execution window is canceled, and the process of ensuring orbit decay is initiated.

2. The method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for generating initial deorbit candidate schemes for the satellite to be deorbited, and constructing the deorbit crossing envelope through orbital dynamics extrapolation, are as follows: The satellite to be deorbited performs orbital nonlinear dynamics integration based on the initial deorbit candidate scheme to obtain the predicted position vector during the deorbiting process; Using the orbital dynamics state transition matrix for covariance time propagation, the time evolution sequence of the position error covariance matrix corresponding to the predicted position vector is extracted; By combining the time evolution sequence of the predicted position vector and the position error covariance matrix, a spatial error boundary covering the descending trajectory is constructed in three-dimensional space, and the derailment crossing envelope is constructed using the three-dimensional derailment crossing envelope calculation formula.

3. The method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for selecting associated satellites that physically overlap with the deorbit crossing envelope and directionally sending deorbit intention data frames are as follows: Extract the constellation's global approximate ephemeris data stored on-board and estimate the space occupancy range of the remaining active satellite nodes within the constellation; Nodes that physically overlap with the space-occupied interval and the off-orbit crossing envelope are selected, and the nodes that physically overlap are identified as the associated satellites; The node physical identifier, the timestamp sequence corresponding to the system global time variable, the time evolution sequence of the predicted position vector of the satellite to be deorbited, and the time evolution sequence of the position error covariance matrix are combined and packaged to generate and send the deorbit intention data frame in a targeted manner.

4. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for assessing the rendezvous risk of associated satellites are as follows: The deorbiting intention data frame received by the associated satellite is parsed and a time-parameterized relative trajectory search is performed within a preset rendezvous time window to calculate the closest approach time; Extract the relative position and relative velocity near the closest approach time, and superimpose the position error covariance of the satellite to be deorbited with its own position error covariance to generate a joint position error covariance matrix; The joint position error covariance matrix is ​​symmetrically processed using the joint position error covariance matrix symmetry formula; then, the eigenvalue lower bound constraint is regularized using the joint position error covariance matrix positive definite formula, and the rendezvous risk is quantitatively evaluated based on the regularization result.

5. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for generating the parameters of the two-layer four-dimensional spatiotemporal repulsion corridor, which includes an inner rigid safety boundary and an outer elastic service boundary, are as follows: Based on the rendezvous risk and the equivalent radius of satellite physical collision avoidance obtained from the assessment, the insurmountable inner rigid safety boundary is established using the inner rigid safety boundary shape matrix calculation formula. Combining task interference avoidance margin, the outer boundary parameters are generated using the outer elastic service boundary shape matrix calculation formula; based on the outer boundary parameters, the outer elastic service boundary, which can be controlled to shrink within a preset range, is constructed outside the inner rigid safety boundary; The inner rigid safety boundary and the outer elastic service boundary are combined and encapsulated with their own position state sequences to generate the double-layer four-dimensional spatiotemporal repulsion corridor parameters.

6. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for the satellite to be deorbited to receive the dual-layer four-dimensional spatiotemporal exclusion corridor parameters and node cost data, and generate a constellation-level global deorbit cost function with a dynamic constraint relaxation penalty term, are as follows: The satellite to be deorbited analyzes the parameters of the two-layer four-dimensional spatiotemporal exclusion corridor and assigns heterogeneous cost weights to each satellite node in combination with the node cost data. The physical consumption quantities are normalized, and the normalized relative cost score is multiplied by the cost weight to calculate the weighted cost value of each node. The weighted maneuver cost of the satellite to be deorbited, the weighted secondary generation cost of the associated satellite set, and the dynamic constraint relaxation penalty term are integrated, and the constellation-level global deorbit cost function is generated using the constellation-level global deorbit cost function calculation formula.

7. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, Before performing controlled shrinkage only on the outer elastic service boundary, when the outer elastic service boundary has no feasible solution as a hard constraint, the process further includes: The outer elastic service boundary is transformed into a spatial avoidance matrix constraint for optimization. If the state of no feasible solution is continuously detected, the Lagrange multiplier values ​​of each constraint in the last iteration are extracted. Combining the Lagrange multiplier values, constraint violation amounts, and constraint scale normalization factors, the stagnation score of each associated satellite is calculated using the stagnation score calculation formula. The associated satellite with the highest obstruction score is identified as the core obstruction node, and the controlled contraction is prepared to be executed against the core obstruction node.

8. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 7, characterized in that, The specific steps for performing controlled contraction only on the outer elastic service boundary, converting the contraction amount into a numerical penalty injected into the constellation-level global deorbit cost function, maintaining the inner rigid safety boundary from relaxation, and obtaining the cooperative compromise solution are as follows: Controlled contraction is performed on the outer elastic service boundary for the core blocking node. The outer boundary clearance deficiency is calculated using the formula for calculating the outer boundary clearance deficiency, and a constant relaxation variable is obtained as the contraction amount. The inner boundary retention constraint formula is used to ensure that the inner rigid safety boundary cannot be relaxed, while the shrinkage amount is converted into the numerical penalty by incrementally increasing the dynamic penalty benchmark coefficient. The numerical penalty is injected into the constellation-level global deorbit cost function and the optimization process is restarted until the convergence condition is met to obtain the cooperative compromise solution.

9. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for converting the cooperative compromise solution into an action sequence and sending it to the associated satellite for execution are as follows: The cooperative compromise solution is analyzed and converted into discrete control instructions applicable to different hardware configurations. For associated satellites using active orbit control, the discrete control commands are converted into ignition pulse commands, and the corresponding action timing is generated and sent to the associated satellite for execution. For associated satellites employing differential drag control, the discrete control commands are converted into aerodynamic drag control commands, and the corresponding action timing is generated and sent to the associated satellites for execution.

10. A method for coordinated deorbiting of multiple low-Earth orbit constellations according to claim 1, characterized in that, The specific steps for canceling the action sequence that has not yet entered the irreversible execution window and switching to the backup track attenuation process in the event of communication abnormality or inability to confirm the action sequence in a closed loop are as follows: When the system monitoring log generates an anomaly flag, the multi-star collaborative planning mechanism is terminated, and the action sequence that has not yet entered the irreversible execution window is cut off or canceled. If the attitude of the satellite to be deorbited is controllable, the drag-increasing component is adjusted to increase the equivalent windward area along the velocity direction, forming a continuous aerodynamic drag increment to reduce the orbital altitude and enter the minimum orbital attenuation process. If the attitude control actuator of the satellite to be deorbited fails, and the satellite to be deorbited is equipped with a passive deorbiting device with independent triggering capability or preset triggering capability, then the passive deorbiting device is triggered; if the passive deorbiting device cannot be triggered, then the trajectory is extrapolated according to the ballistic passive attenuation trajectory, and downgraded to a one-way passive avoidance mode, and the process of ensuring the minimum orbital attenuation is initiated.

Images