Space debris removal task planning method based on dynamic clustering and reachable graph network

Abstract

The invention discloses a space debris removal task planning method based on a dynamic clustering and reachable graph network, which comprises the following steps: acquiring initial positions of all candidate debris, and calculating six orbit elements of all candidate debris on each time step of a task time window to form a dynamic ephemeris database; clustering all candidate fragments corresponding to each time step by adopting a DBSCAN clustering method, and taking each obtained cluster as a candidate path; eliminating a candidate path of which the spacecraft speed meets the relative speed constraint of all candidate fragments in all candidate paths to obtain a feasible path set; taking candidate paths in the feasible path set as nodes, and generating all possible node pairs; carrying out feasibility verification on the node pairs, and forming a reachable graph network by adopting the node pairs passing the feasibility verification; and according to the reachable graph network, performing path optimization by adopting a genetic algorithm, and generating a globally optimal candidate fragment removal task path for removing the most candidate fragments.

Description

Technical Field

[0001] This invention relates to task path planning technology, specifically to a spatial fragmentation removal task planning method based on dynamic clustering and reachability graph networks. Background Technology

[0002] In recent years, with the increasing frequency of human space activities, the number of defunct spacecraft, rocket debris, and space debris generated by collisions left in near-Earth orbit has increased dramatically. This high-speed space debris has become a serious threat to the safe operation of spacecraft in orbit and may trigger the "Kessler effect," where debris from collisions further exacerbates the risk of subsequent collisions, creating a vicious cycle that ultimately renders certain orbital regions unusable. Therefore, conducting Active Debris Removal (ADR) missions is of paramount strategic importance for ensuring the safety of space assets and maintaining the sustainable development of the space environment.

[0003] For space debris removal missions, the mainstream orbital design schemes can be divided into two categories: one is pulse propulsion, which uses chemical thrusters to generate instantaneous high thrust to achieve rapid orbital maneuvering; the other is continuous low-thrust propulsion, which uses electric propulsion and other high specific impulse systems to achieve efficient orbital transfer through long-term operation. In terms of mission mode, to improve economic efficiency, a "multi-target" mode is usually adopted, in which a single or a small number of servicing spacecraft sequentially removes multiple debris targets.

[0004] However, planning such multi-target, multi-spacecraft, long-duration debris removal missions is an extremely complex global optimization problem, a typical NP-hard problem, with its core difficulty lying in: The vast combinatorial search space: Mission planning requires not only determining which targets to eliminate from a huge set of candidate debris (ranging from hundreds to thousands) (target selection), but also determining the optimal order in which to eliminate these targets (sequence planning), and allocating a suitable subset of targets to each mission spacecraft (mission assignment). As the number of targets increases, the possible combinatorial paths grow exponentially, making exhaustive search impossible.

[0005] Highly Dynamic and Coupled Transfer Costs: Unlike the classic static Traveling Salesman Problem (TSP), the "distance" between any two objectives in a space mission (i.e., the fuel consumption required for orbital transfer, typically quantified as a velocity increment ΔV) changes dynamically over time. Transfer costs depend not only on the departure and arrival orbits but also strongly on the departure time and transfer duration. This time dependence tightly couples path planning with time planning, significantly increasing the difficulty of solving the problem.

[0006] High-precision orbital dynamics models are required because low-Earth orbit missions are significantly affected by various perturbations, such as Earth's non-spherical gravity (mainly J2 perturbation) and atmospheric drag. Over mission cycles lasting months or even years, these perturbations can cause long-term effects such as natural precession of the orbital plane. Traditional methods, such as solving the Lambert problem based on two-body models, produce significant errors when calculating long-duration transfers, failing to meet the requirements for accurate planning. While high-precision numerical integration models are accurate, the computational burden is enormous; in scenarios where global optimization algorithms require tens of thousands or even millions of iterations, the computational time becomes unacceptable.

[0007] Currently, the methods for solving this type of problem have the following limitations: Heuristic and Greedy Algorithms: Greedy strategies designed based on dynamic laws such as "minimum difference in right ascension of ascending nodes" can quickly provide feasible solutions, but they are prone to getting trapped in local optima and cannot guarantee the global quality of the solution.

[0008] Traditional evolutionary algorithms (such as genetic algorithms and ant colony algorithms), while possessing global search capabilities, are highly dependent on the efficiency and accuracy of distance evaluation when dealing with such problems. Using a high-precision model results in excessive evaluation costs, leading to slow convergence; while using a coarse approximation model (such as ignoring perturbations or simplifying maneuver patterns) can cause significant deviations between the planning results and the actual optimal solution due to model mismatch.

[0009] Mixed Integer Nonlinear Programming (MINLP): Although this method can accurately model the problem, its solution efficiency drops sharply with the increase of variable dimensionality (especially integer variables), making it difficult to apply to scenarios with large-scale target sets.

[0010] In summary, existing technologies generally suffer from a trade-off between efficiency and optimality when dealing with global optimization planning problems for large-scale space debris removal tasks. Summary of the Invention

[0011] To address the aforementioned shortcomings in existing technologies, the spatial debris removal task planning method based on dynamic clustering and reachability graph networks provided by this invention solves the technical problem of balancing efficiency and accuracy in existing spatial debris removal task planning methods.

[0012] To achieve the above-mentioned objectives, the technical solution adopted by this invention is as follows: A spatial debris removal task planning method based on dynamic clustering and reachability graph networks is provided, which includes the following steps: S1. Obtain the initial positions of all candidate fragments and calculate the orbital six-element number of all candidate fragments at each time step of the mission time window to form a dynamic ephemeris database. S2. Based on the dynamic ephemeris database, the DBSCAN clustering method is used to cluster all candidate fragments corresponding to each time step, and each cluster is used as a candidate path. S3. Eliminate all candidate paths where none of the spacecraft velocities satisfy the relative velocity constraints of all its candidate debris, and obtain a set of feasible paths; S4. Using the candidate paths in the feasible path set as nodes, generate all possible node pairs. There are no overlapping candidate fragments in the node pairs, and the time of the starting node is earlier than that of the ending node. S5. Perform feasibility verification on node pairs, and use the node pairs that pass the feasibility verification to form a reachable graph network; S6. Based on the reachability graph network, a genetic algorithm is used to optimize the path and generate a candidate fragment removal task path that removes the most candidate fragments and is globally optimal.

[0013] The beneficial effects of this invention are as follows: This solution can automatically and systematically identify all spatially highly clustered debris clusters throughout the entire mission cycle and form a detailed list of candidate path points, thereby transforming the original, unstructured spatiotemporal point cloud data into a structured set of candidate mission targets with clear physical meaning.

[0014] This approach avoids the pitfalls of traditional methods getting stuck in local optima by constructing a reachability graph network, enabling the search for the globally optimal solution across the entire task space. By clustering fragments, constructing the reachability graph network, and finally combining it with a genetic algorithm, it can efficiently handle complex scenarios involving hundreds or even thousands of fragment targets and task durations spanning months or even years. While ensuring the physical feasibility of the planning scheme and satisfying multiple engineering constraints, it can quickly find the globally optimal or near-optimal task sequence to maximize task benefits (e.g., clearing the maximum number of fragments or minimizing task costs).

[0015] Furthermore, methods for verifying the feasibility of node pairs include: S51. Use the Lambert problem solver to calculate the transfer trajectory between each node pair. When the solver generates the velocity vectors of the two nodes of the node pair, proceed to step S52; otherwise, delete the node pair. S52. Calculate the semi-major axis a of the transfer trajectory based on the position and velocity vector of the starting node in the node pair. t Then calculate the average semi-major axis a of all candidate fragments. a ; S53. Determine whether the transfer paths between node pairs satisfy the economic constraints. If yes, proceed to step S54; otherwise, delete the node pair. The economic constraints are: in, Preset capacity; S54. Determine whether the absolute value of the difference between the velocity vector of a node in a node pair and its corresponding center velocity is less than or equal to its velocity margin. If yes, proceed to step S55; otherwise, delete the corresponding node pair. S55. Check if the perigee distance of the transfer orbit of the node pair is greater than the sum of the Earth's radius and the minimum safe altitude. If so, keep the node pair; otherwise, delete the node pair.

[0016] Furthermore, the space debris removal task planning method also includes the opportunity objective of searching for node pairs: A high-precision numerical integrator is used for trajectory propagation along the transfer trajectory of the node pair, and the transfer time is recorded. At each point in time during the propagation process, calculate the distance between the spacecraft's position and the positions of all candidate debris that do not belong to a node pair in the dynamic ephemeris database; Determine if the distance to the candidate fragment is less than a preset proximity threshold. If so, record it as an opportunity target; otherwise, do not record it.

[0017] Furthermore, the space debris removal mission planning method also includes calculating the pulse maneuvers required by the spacecraft at the starting node of the node pair: in, The pulse maneuver required by the spacecraft at node i, the starting point of the node pair; Let i be the velocity vector of the starting node i in the node pair; Let be the center velocity of node i, the starting node of the node pair; It is the Euclidean norm; Record the attribute information of the valid edges between node pairs: the index of the target node in the node pair, and the transition time t. f Pulse maneuver The velocity vectors of the starting node i and the target node j in the node pair, as well as the number and ID list of the corresponding opportunity targets.

[0018] The beneficial effects of the above technical solution are as follows: by determining the node pairs in this way, a complex, multi-constraint trajectory design problem can be systematically and automatically transformed into a discrete network model that includes all feasible task segments and their costs and benefits. This model serves as the data foundation for subsequent global optimal path search, greatly reducing the complexity of the final optimization problem.

[0019] When verifying the feasibility of node pairs, this plan comprehensively considers multiple constraints such as time, location, speed, fuel (indirectly through ΔV), and orbital altitude. The planning scheme is complete and engineering feasible.

[0020] Furthermore, the initial population generation method of the genetic algorithm in step S6 includes: S61. Based on the reachable graph network, select a node with a non-zero out-degree as the starting point of the task path, and select a node from the valid neighbors of the starting point of the task path that is not repeated with the node in the task path as the next node. S62. Select a node from the valid neighbors of the next node that does not overlap with the node in the task path as its next node, and repeat the current operation until the length of the task path is equal to the preset length. S63. Repeat steps S61 and S62 until all possible task paths are generated as the initial population, with each task path representing one individual.

[0021] Furthermore, the mutation operation of the genetic algorithm is as follows: randomly select a cutoff point in an individual, retain its first half, take the cutoff point as the next node, and repeat step S62 until the path length is equal to the preset length. The expression for the fitness function in a genetic algorithm is: in, The fitness function; , and All are preset weights; The total non-repeating candidate fragment revenue; =K is the candidate path length reward, and K is the preset length; The time span between the starting node and the target node in an individual; and These are the k-th and (k-1)-th nodes on the path corresponding to the individual, respectively. For nodes The set of candidate fragments covered; for and Opportunity targets on the edge between; and These are the symbols for performing a union operation on the fragment set of all nodes on the path and the chance target set of all edges, respectively. This is an extraction operation used to obtain the set of opportunity target fragments recorded on an edge; It is the cardinality of the set.

[0022] The beneficial effects of the above technical solution are as follows: Through the design of the fitness function of the above multi-objective, the genetic algorithm can be guided to find the comprehensive optimal solution that not only clears more fragments, but also has a long task sequence, full planning, and better meets the actual needs of engineering, effectively avoiding the "short-sighted" behavior of simple greedy algorithms.

[0023] Furthermore, step S3 further includes: S31. Based on the total number M of candidate fragments in the candidate path, define M velocity spheres in the three-dimensional velocity space. Each velocity sphere has the velocity of a candidate fragment as its center and the velocity threshold as its radius. S32. A minimization-maximization optimization problem to determine whether the intersection of M spheres is non-empty: in, =[v x ,v y ,v z [ ] represents the three-dimensional velocity vector of the spacecraft to be optimized; To minimize the parameter operator; Let R be the velocity of the m-th candidate fragment in the cluster; R is the radius of the smallest enclosing sphere of the intersection of the M velocity spheres. It is the Euclidean norm; S33. Take the arithmetic mean of the velocities of all candidate fragments within the cluster as... The initial value is used to solve the minimization-maximization optimization problem using a constrained nonlinear programming solver, yielding the optimal value. and radius R; S34. Determine whether the radius R is greater than or equal to the velocity threshold. If so, then there is no spacecraft whose velocity satisfies the relative velocity constraint of all its candidate fragments, and it is removed. Otherwise, it is retained. S35. Use all the retained candidate paths to form a set of feasible paths.

[0024] The beneficial effects of the above technical solution are as follows: it determines the "optimal center velocity" and "velocity margin" for each feasible path, providing key quantitative input for global path planning; it transforms the complex constraint determination problem into an efficient numerical optimization problem, improving screening efficiency and automation; it achieves precise screening from "geometrically possible" to "physically feasible", significantly improving the effectiveness of subsequent planning; and it enhances the robustness and reliability of the entire planning framework.

[0025] Furthermore, step S1 further includes: Obtain the Cartesian state data of all candidate fragments at the initial moment as the initial position, and generate the time sampling point sequence of each candidate fragment according to the task's time window and discrete time step. Based on the time sampling point sequence, using The perturbation model calculates the orbital six-roots number of each candidate fragment at each time step sampling point, and uses the orbital six-roots numbers of all candidate fragments at all sampling points to form a dynamic ephemeris database.

[0026] Furthermore, the methods for obtaining the six base numbers of the orbital path include: exist In long-term perturbation theory, the orbital energy, orbital shape, and orbital inclination, which are the six orbital roots, do not undergo first-order changes during long-term evolution. Let the semi-major axis a at any time step t be... t eccentricity e t Track inclination angle i t All are equal to their initial values; The methods for calculating the right ascension of the ascending node, the argument of perigee, and the mean perigee angle in the six roots of the orbit at any time step t include: Calculation by The long-term average rate of change of the right ascension Ω of the ascending node and the argument ω of the perigee caused by the perturbation: in, and These are the long-term average rates of change of the right ascension Ω of the ascending node and the argument ω of the perigee, respectively; for Term coefficient; p is the radius of the Earth's equator; p is the semi-major diameter. Based on the average angular velocity n of the orbit, Using the term coefficient and the Earth's equatorial radius, calculate the average angular velocity: in, The average angular velocity; This is the square root operator. Calculate the right ascension of the ascending node at any time step t. Perigeal argument Peace Angle : , , ; right , and Perform a modulo operation to obtain the angle-normalized parameters.

[0027] The beneficial effects of the above technical solution are as follows: This solution uses the above method to solve the orbital six roots, avoiding complex and time-consuming numerical integration. The long-term evolution results of the orbital six roots can be obtained quickly through simple algebraic operations, which greatly improves the computational efficiency and makes it very suitable for batch forecasting of a large number of targets over a long period of time. Attached Figure Description

[0028] Figure 1This is a flowchart of a spatial debris removal task planning method based on dynamic clustering and reachability graph networks.

[0029] Figure 2 This is a schematic diagram of the clustering of candidate path points found when 1T=63 minutes in the example.

[0030] Figure 3 This is a schematic diagram of the velocity feasible region analysis within a path point. Detailed Implementation

[0031] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0032] refer to Figure 1 , Figure 1 A flowchart of a spatial debris removal task planning method based on dynamic clustering and reachability graph networks is shown; Figure 1 As shown, the method S includes steps S1 to S6.

[0033] In step S1, the initial positions of all candidate fragments are obtained, and the orbital six-root numbers of all candidate fragments at each time step of the mission time window are calculated to form a dynamic ephemeris database. In implementation, step S1 of this solution further includes: Obtain the Cartesian state data of all candidate fragments at the initial moment as the initial position, and generate the time sampling point sequence of each candidate fragment according to the task's time window and discrete time step. Based on the time sampling point sequence, using The perturbation model calculates the orbital six-roots number of each candidate fragment at each time step sampling point, and uses the orbital six-roots numbers of all candidate fragments at all sampling points to form a dynamic ephemeris database.

[0034] The methods for obtaining the six base numbers of the orbital path include: exist In long-term perturbation theory, the orbital energy, orbital shape, and orbital inclination, which are the six orbital roots, do not undergo first-order changes during long-term evolution. Let the semi-major axis a at any time step t be... t eccentricity e t Track inclination angle i t All are equal to their initial values; The methods for calculating the right ascension of the ascending node, the argument of perigee, and the mean perigee angle in the six roots of the orbit at any time step t include: Calculation by The long-term average rate of change of the right ascension Ω of the ascending node and the argument ω of the perigee caused by the perturbation: in, and These are the long-term average rates of change of the right ascension Ω of the ascending node and the argument ω of the perigee, respectively; for Term coefficient; p is the radius of the Earth's equator; p is the semi-major diameter. Based on the average angular velocity n of the orbit, Using the term coefficient and the Earth's equatorial radius, calculate the average angular velocity: in, The average angular velocity; This is the square root operator. Calculate the right ascension of the ascending node at any time step t. Perigeal argument Peace Angle : , , ; right , and Perform a modulo operation to obtain the angle-normalized parameters.

[0035] In step S2, based on the dynamic ephemeris database, the DBSCAN clustering method is used to cluster all candidate fragments corresponding to each time step, and each cluster is used as a candidate path.

[0036] In one embodiment of the present invention, the detailed implementation process of DBSCAN clustering is as follows: Data preparation: Extract the 3D position coordinates P of all N fragments at the current time step t from the dynamic ephemeris database. i (t)=[x i ,y i ,z i ], forming an N x3 The data matrix.

[0037] DBSCAN Algorithm Parameter Mapping: Directly mapping the physical constraints of the task to the core parameters of the DBSCAN algorithm: Neighborhood Radius (Epsilon, ε): This parameter is set to the maximum approach distance ε required by the task. distFor example, 30 kilometers. In the algorithm, it defines a spherical space with radius ε centered on a certain fragment. Minimum number of points (MinPts): This parameter is set to the minimum number of fragments required by the task to form a valuable clearing cluster, for example, 2. It defines how many other points (including itself) a point must contain within its ε-neighborhood to be considered a "core point".

[0038] Core point identification and cluster growth: The DBSCAN algorithm traverses all fragment points and makes the following judgment based on the above parameters: For any fragment P, calculate the number of fragments in its ε-neighborhood. If the number of fragments in the neighborhood is greater than or equal to MinPts, then fragment P is marked as a core point. Starting from any unvisited core point, a new cluster is created. Then, through density reachability, all points reachable from this core point (i.e., within the ε-neighborhood of a core point) (including other core points and non-core "boundary points") are assigned to this cluster. This process is iterated until the cluster can no longer grow.

[0039] Candidate pathpoint generation: After the clustering process is completed at time step t, all fragment sets belonging to the same cluster are identified as potential purge opportunities. For each identified cluster k, the system calculates its geometric center (i.e., the arithmetic mean of all fragment locations within the cluster) as the representative location of the cluster, denoted as P. k (t). Simultaneously, record the time t when the cluster occurs, and the ID list D of all fragments contained in the cluster. k These three sets of information {t, P} k (t),D k Together, they form a candidate path point.

[0040] By repeating the above process at all time sampling points, all spatially highly clustered debris clusters throughout the entire mission cycle can be automatically and systematically identified, forming a detailed list of candidate waypoints. This list forms the basis for subsequent velocity feasibility analysis and mission planning. This step transforms the raw, unstructured spatiotemporal point cloud data into a structured set of candidate mission targets with clear physical meaning.

[0041] The advantages of using DBSCAN are: Strong shape adaptability: It can identify debris clusters of arbitrary shapes, which is very consistent with the distribution characteristics of linear, arc-shaped or other irregular shapes that may be formed when debris is scattered along the orbit after a satellite disintegration event. This is more advantageous than algorithms such as K-Means, which can only identify spherical clusters.

[0042] No need to preset the number of clusters: DBSCAN can automatically discover all clustered regions that meet the density conditions based on the distribution characteristics of the data itself, without the need for manual pre-setting of how many clusters to search, which makes the search process more automated and objective.

[0043] We obtained 345 fragments, set the task time window to 24 hours, and the discrete time step to one minute. We then used steps S1-S2 of this scheme to cluster the 345 fragments. The clustering results can be referenced. Figure 2 ;exist Figure 2 Each five-pointed star represents a cluster formed at 63 minutes.

[0044] In step S3, candidate paths that do not contain a spacecraft whose velocity satisfies the relative velocity constraints of all its candidate debris are eliminated, resulting in a set of feasible paths. In implementation, the preferred step S3 of this scheme further includes: S31. Based on the total number M of candidate fragments in the candidate path, define M velocity spheres in the three-dimensional velocity space. Each velocity sphere has the velocity of a candidate fragment as its center and a velocity threshold as its radius. The preferred velocity threshold in this scheme is 150 m / s.

[0045] S32. A minimization-maximization optimization problem to determine whether the intersection of M spheres is non-empty: in, =[v x ,v y ,v z [ ] represents the three-dimensional velocity vector of the spacecraft to be optimized; To minimize the parameter operator; Let R be the velocity of the m-th candidate fragment in the cluster; R is the radius of the smallest enclosing sphere of the intersection of the M velocity spheres. It is the Euclidean norm; The minimization-maximization optimization problem designed in this scheme can efficiently determine whether the intersection of M spheres is non-empty and find an "optimal" point in the intersection as the recommended velocity for the spacecraft; the goal of the optimization problem is to find an optimal spacecraft velocity. This minimizes the maximum distance from the velocity to the velocity vector of all fragments within the cluster.

[0046] S33. Take the arithmetic mean of the velocities of all candidate fragments within the cluster as... The initial value is used to solve the minimization-maximization optimization problem using a constrained nonlinear programming solver (preferably the fmincon function in MATLAB) to obtain the optimal value. and radius R; S34. Determine whether the radius R is greater than or equal to the velocity threshold. If so, then there is no spacecraft whose velocity satisfies the relative velocity constraint of all its candidate fragments, and it is removed. Otherwise, it is retained. S35. Use all the retained candidate paths to form a set of feasible paths.

[0047] For each feasible candidate path that passes the evaluation, the system will record and store the following key information: 1) time_min: the time of occurrence (i.e., the time corresponding to the starting node among the feasible candidate paths); 2) position_km: spatial location (cluster centroid); 3) covered_debris_ids: the list of covered debris IDs; 4) feasible_velocity_ms: the optimal spacecraft velocity obtained from the optimization solution. 5) velocity_margin_ms: velocity margin, defined as velocity threshold minus radius R. This value represents the spacecraft's available velocity adjustment space at this point without violating any constraints, and is an important parameter for subsequent orbit transfer planning.

[0048] Figure 3 This invention demonstrates a key calculation result from an embodiment of the present invention, specifically the calculation result of step 3. Assuming information on 345 fragments is obtained, the task time window is set to 24 hours, and the discrete time step is one minute, after executing steps S1 to S3 of this scheme, the following result is formed: Figure 3 The results shown depict the three-dimensional velocity space of a feasible pathpoint identified by the system at t=197.0 minutes. This cluster comprises fragments 88, 139, 166, 261, and 298.

[0049] In this three-dimensional coordinate system, the three coordinate axes Vx, Vy, and Vz represent the three components of the velocity vector, with units of m / s. Each element in the diagram has a clear physical meaning: The five small dots of different colors in the diagram represent the actual velocity vectors of the five fragments within the cluster at that particular moment. The coordinates V of these points are... m =[vx m ,vy m ,vz m [It is precisely extracted from the dynamic ephemeris database generated in step S1.]

[0050] The translucent spheres of the same color surrounding each small colored dot represent the permissible velocity range of that fragment. The center of each sphere is the velocity vector V of the corresponding fragment. m Its radius is the maximum relative velocity constraint V preset by the task. t(In this embodiment, it is 150 m / s). Any velocity point V located inside the sphere. sc All satisfy ||V sc -V m ||≤150m / s conditions.

[0051] The dark area formed by the overlapping of all the translucent spheres in the diagram represents the core objective of this scheme—the globally permissible velocity domain of this pathpoint. Any velocity point located within this overlapping region simultaneously satisfies the relative velocity constraints with all five fragments within the cluster. The non-emptiness of this region provides intuitive geometric proof that this pathpoint is considered "velocity feasible."

[0052] The yellow pentagram located at the center of the intersection region of all spheres represents the optimal spacecraft recommended velocity V calculated by this invention through a minimax optimization problem. sc This point is the point with the smallest maximum distance to the velocity vectors of all debris, meaning it is the geometric center of the globally permissible velocity domain. Choosing this point as the spacecraft's target velocity provides the maximum velocity margin for subsequent orbital maneuvers and control.

[0053] To facilitate understanding of the minimization-maximization optimization problem in this scheme, the components of the formula are explained below: Represents two three-dimensional vectors and The Euclidean distance or norm between them is expressed as the relative velocity (rate) between the spacecraft and the m-th debris, which is the basic unit for evaluating the cost of a maneuver.

[0054] For the maximum value operator, it is applicable to a given... Calculate the relative velocity of the spacecraft with the velocities of all M fragments within the cluster, and return the maximum value. This expression represents the velocity of the spacecraft if it chooses a specific velocity... This will present the most severe speed matching challenge. In other words, it represents the "worst-case speed cost" or "maximum relative speed" for performing a cleanup operation on the entire cluster. Regardless of... Regardless of the choice, the spacecraft must be able to overcome at least this maximum speed difference.

[0055] This is the Argument of the Minimum operator. It doesn't seek the minimum value of the "maximum relative velocity" itself, but rather the input parameter that minimizes this "maximum relative velocity," i.e., the optimal spacecraft velocity vector. This operator represents a global search process. The system searches in all possible three-dimensional velocity spaces R. 3 The search is underway to find an "ideal" spacecraft speed. The definition of this "ideal" is: when the spacecraft is at this speed, the worst-case scenario it faces (i.e., the maximum speed difference with a fragment in the cluster) is minimized.

[0056] In step S4, candidate paths from the feasible path set are used as nodes to generate all possible node pairs. These node pairs do not contain overlapping candidate fragments, and the starting node's time is earlier than the ending node's time. Specifically: Check if the fragment cluster sets represented by the two nodes i and j in the node pair have an intersection. If there is an intersection, the two nodes are considered to be mutually exclusive in a task sequence and no edge is established. The time t for checking the endpoint j j Is it later than time t of starting point i? i Only when the transition time t f =t j -t i The transition is only meaningful when the value is greater than 0; extract the position of the starting point i, r1=P. i And the position of the endpoint j, r2=P j .

[0057] In step S5, the feasibility of node pairs is verified, and the node pairs that pass the feasibility verification are used to construct a reachable graph network; the method for verifying the feasibility of node pairs includes: S51. Use the Lambert problem solver to calculate the transfer trajectory between each node pair. When the solver generates the velocity vectors of the two nodes of the node pair, proceed to step S52; otherwise, delete the node pair. S52. Calculate the semi-major axis a of the transfer trajectory based on the position and velocity vector of the starting node in the node pair. t Then calculate the average semi-major axis a of all candidate fragments. a ; S53. Determine whether the transfer paths between node pairs satisfy the economic constraints. If yes, proceed to step S54; otherwise, delete the node pair. The economic constraints are: in, Preset capacity; S54. Determine whether the absolute value of the difference between the velocity vector of a node in a node pair and its corresponding center velocity is less than or equal to its velocity margin. If yes, proceed to step S55; otherwise, delete the corresponding node pair. S55. Check if the perigee distance of the transfer orbit of the node pair is greater than the sum of the Earth's radius and the minimum safe altitude. If so, keep the node pair; otherwise, delete the node pair.

[0058] Space debris removal task planning methods also include opportunity objectives for searching node pairs: A high-precision numerical integrator is used for trajectory propagation along the transfer trajectory of the node pair, and the transfer time is recorded. At each point in time during the propagation process, calculate the distance between the spacecraft's position and the positions of all candidate debris that do not belong to a node pair in the dynamic ephemeris database; Determine whether the distance to the candidate fragment is less than a preset proximity threshold (preferably 30 kilometers). If so, record it as an opportunity target; otherwise, do not record it.

[0059] Space debris removal mission planning methods also include calculating the pulse maneuvers required by the spacecraft at the starting node of the node pair: in, The pulse maneuver required by the spacecraft at node i, the starting point of the node pair; Let i be the velocity vector of the starting node i in the node pair; Let be the center velocity of node i, the starting node of the node pair; It is the Euclidean norm; Record the attribute information of the valid edges between node pairs: the index of the target node in the node pair, and the transition time t. f Pulse maneuver The velocity vectors of the starting node i and the target node j in the node pair, as well as the number and ID list of the corresponding opportunity targets.

[0060] In step S6, based on the reachability graph network, a genetic algorithm is used to optimize the path and generate a candidate fragment removal task path that removes the most candidate fragments and is globally optimal.

[0061] In one embodiment of the present invention, the initial population generation method of the genetic algorithm in step S6 includes: S61. Based on the reachable graph network, select a node with a non-zero out-degree as the starting point of the task path, and select a node from the valid neighbors of the starting point of the task path that is not repeated with the node in the task path as the next node. S62. Select a node from the valid neighbors of the next node that does not overlap with the node in the task path as its next node, and repeat the current operation until the length of the task path is equal to the preset length. S63. Repeat steps S61 and S62 until all possible task paths are generated as the initial population, with each task path representing one individual.

[0062] The mutation operation of the genetic algorithm is as follows: randomly select a cutoff point in an individual, retain its first half, take the cutoff point as the next node, and repeat step S62 until the path length is equal to the preset length. The expression for the fitness function in a genetic algorithm is: in, The fitness function; , and All are preset weights; The total non-repeating candidate fragment revenue; =K is the candidate path length reward, and K is the preset length; The time span between the starting node and the target node in an individual; and These are the k-th and (k-1)-th nodes on the path corresponding to the individual, respectively. For nodes The set of candidate fragments covered; for and Opportunity targets on the edge between; and These are the symbols for performing a union operation on the fragment set of all nodes on the path and the chance target set of all edges, respectively. This is an extraction operation used to obtain the set of opportunity target fragments recorded on an edge; It is the cardinality of the set.

[0063] In summary, this solution transforms the complex and continuous trajectory design problem into a structured graph theory optimization problem, which can efficiently and systematically provide globally optimal planning solutions for multi-objective and multi-constraint debris removal tasks, significantly improving the efficiency of task planning and task benefits.

Claims

1. A spatial fragmentation removal task planning method based on dynamic clustering and reachability graph networks, characterized in that, Including the following steps: S1. Obtain the initial positions of all candidate fragments and calculate the orbital six-element number of all candidate fragments at each time step of the mission time window to form a dynamic ephemeris database. S2. Based on the dynamic ephemeris database, the DBSCAN clustering method is used to cluster all candidate fragments corresponding to each time step, and each cluster is used as a candidate path. S3. Eliminate all candidate paths where none of the spacecraft velocities satisfy the relative velocity constraints of all its candidate debris, and obtain a set of feasible paths; S4. Using the candidate paths in the feasible path set as nodes, generate all possible node pairs. There are no overlapping candidate fragments in the node pairs, and the time of the starting node is earlier than that of the ending node. S5. Perform feasibility verification on node pairs, and use the node pairs that pass the feasibility verification to form a reachable graph network; S6. Based on the reachability graph network, a genetic algorithm is used to optimize the path and generate a candidate fragment removal task path that removes the most candidate fragments and is globally optimal.

2. The space debris removal task planning method according to claim 1, characterized in that, Methods for verifying the feasibility of node pairs include: S51. Use the Lambert problem solver to calculate the transfer trajectory between each node pair. When the solver generates the velocity vectors of the two nodes of the node pair, proceed to step S52; otherwise, delete the node pair. S52. Calculate the semi-major axis a of the transfer trajectory based on the position and velocity vector of the starting node in the node pair. t Then calculate the average semi-major axis a of all candidate fragments. a ; S53. Determine whether the transition paths between node pairs satisfy economic constraints. If so, proceed to step [step number missing]. S54, otherwise delete the node pair; the economic constraint is: in, Preset capacity; S54. Determine whether the absolute value of the difference between the velocity vector of a node in a node pair and its corresponding center velocity is less than or equal to its velocity margin. If yes, proceed to step S55; otherwise, delete the corresponding node pair. S55. Check if the perigee distance of the transfer orbit of the node pair is greater than the sum of the Earth's radius and the minimum safe altitude. If so, keep the node pair; otherwise, delete the node pair.

3. The space debris removal task planning method according to claim 2, characterized in that, It also includes the opportunity objective of searching for node pairs: A high-precision numerical integrator is used for trajectory propagation along the transfer trajectory of the node pair, and the transfer time is recorded. At each point in time during the propagation process, calculate the distance between the spacecraft's position and the positions of all candidate debris that do not belong to a node pair in the dynamic ephemeris database; Determine if the distance to the candidate fragment is less than a preset proximity threshold. If so, record it as an opportunity target; otherwise, do not record it.

4. The space debris removal task planning method according to claim 3, characterized in that, It also includes calculating the pulse maneuvers required by the spacecraft at the starting node of the node pair: in, The pulse maneuver required by the spacecraft at node i, the starting point of the node pair; Let i be the velocity vector of the starting node i in the node pair; Let be the center velocity of node i, the starting node of the node pair; It is the Euclidean norm; Record the attribute information of the valid edges between node pairs: the index of the target node in the node pair, and the transition time t. f Pulse maneuver The velocity vectors of the starting node i and the target node j in the node pair, as well as the number and ID list of the corresponding opportunity targets.

5. The space debris removal task planning method according to claim 3, characterized in that, The initial population generation method of the genetic algorithm in step S6 includes: S61. Based on the reachable graph network, select a node with a non-zero out-degree as the starting point of the task path, and select a node from the valid neighbors of the starting point of the task path that is not repeated with the node in the task path as the next node. S62. Select a node from the valid neighbors of the next node that does not overlap with the node in the task path as its next node, and repeat the current operation until the length of the task path is equal to the preset length. S63. Repeat steps S61 and S62 until all possible task paths are generated as the initial population, with each task path representing one individual.

6. The space debris removal task planning method according to claim 5, characterized in that, The mutation operation of the genetic algorithm is as follows: randomly select a cutoff point in an individual, retain its first half, take the cutoff point as the next node, and repeat step S62 until the path length is equal to the preset length. The expression for the fitness function in a genetic algorithm is: in, The fitness function; , and All are preset weights; The total non-repeating candidate fragment revenue; =K is the candidate path length reward, and K is the preset length; The time span between the starting node and the target node in an individual; and These are the k-th and (k-1)-th nodes on the path corresponding to the individual, respectively. For nodes The set of candidate fragments covered; for and Opportunity targets on the edge between; and These are the symbols for performing a union operation on the fragment set of all nodes on the path and the chance target set of all edges, respectively. This is an extraction operation used to obtain the set of opportunity target fragments recorded on an edge; It is the cardinality of the set.

7. The space debris removal task planning method according to claim 3, characterized in that, Step S3 further includes: S31. Based on the total number M of candidate fragments in the candidate path, define M velocity spheres in the three-dimensional velocity space. Each velocity sphere has the velocity of a candidate fragment as its center and the velocity threshold as its radius. S32. A minimization-maximization optimization problem to determine whether the intersection of M spheres is non-empty: in, =[v x ,v y ,v z [ ] represents the three-dimensional velocity vector of the spacecraft to be optimized; To minimize the parameter operator; Let R be the velocity of the m-th candidate fragment in the cluster; R is the radius of the smallest enclosing sphere of the intersection of the M velocity spheres. It is the Euclidean norm; S33. Take the arithmetic mean of the velocities of all candidate fragments within the cluster as... The initial value is used to solve the minimization-maximization optimization problem using a constrained nonlinear programming solver, yielding the optimal value. and radius R; S34. Determine whether the radius R is greater than or equal to the velocity threshold. If so, then there is no spacecraft whose velocity satisfies the relative velocity constraint of all its candidate fragments, and it is removed. Otherwise, it is retained. S35. Use all the retained candidate paths to form a feasible path set.

8. The space debris removal task planning method according to claim 1, characterized in that, Step S1 further includes: Obtain the Cartesian state data of all candidate fragments at the initial moment as the initial position, and generate the time sampling point sequence of each candidate fragment according to the task's time window and discrete time step. Based on the time sampling point sequence, using The perturbation model calculates the orbital six-roots number of each candidate fragment at each time step sampling point, and uses the orbital six-roots numbers of all candidate fragments at all sampling points to form a dynamic ephemeris database.

9. The space debris removal task planning method according to claim 8, characterized in that, Methods for obtaining the six base numbers of a orbital include: exist In long-term perturbation theory, the orbital energy, orbital shape, and orbital inclination, which are the six orbital roots, do not undergo first-order changes during long-term evolution. Let the semi-major axis a at any time step t be... t eccentricity e t Track inclination angle i t All are equal to their initial values; The methods for calculating the right ascension of the ascending node, the argument of perigee, and the mean perigee angle in the six roots of the orbit at any time step t include: Calculation by The long-term average rate of change of the right ascension Ω of the ascending node and the argument ω of the perigee caused by the perturbation: in, and These are the long-term average rates of change of the right ascension Ω of the ascending node and the argument ω of the perigee, respectively; for Term coefficient; p is the radius of the Earth's equator; p is the semi-major diameter. Based on the average angular velocity n of the orbit, Using the term coefficient and the Earth's equatorial radius, calculate the average angular velocity: in, The average angular velocity; This is the square root operator. Calculate the right ascension of the ascending node at any time step t. Perigeal argument Peace Angle : ， ， ； right , and Perform a modulo operation to obtain the angle-normalized parameters.

Images