A space debris short-arc data association method and device based on a factor graph

Abstract

This invention belongs to the field of aerospace and space situational awareness technology, specifically relating to a method and apparatus for associating short-arc data of space debris based on factor graphs. The method includes: acquiring and preprocessing optical angle measurement short-arc data from observation stations; obtaining initial orbital state values ​​using an angle observation initial orbit determination algorithm; constructing a cost matrix to initially screen the initial orbital state values ​​corresponding to different short-arc observation data, and generating a limited number of multiple association hypotheses through perturbation sampling expansion; constructing a nonlinear factor graph under each association hypothesis; solving for the maximum a posteriori estimate, and selecting the optimal association hypothesis based on residual and information criteria; outputting the association cluster, orbital state, and its uncertainty, and performing a consistency check. This invention enables efficient and robust multi-arc association and preliminary orbit estimation of optical angle measurement short-arc data in both ground-based and space-based observation scenarios, achieving scene adaptability, real-time performance, and high accuracy.

Description

Technical Field

[0001] This invention belongs to the field of aerospace and space situational awareness technology, specifically relating to a method and device for associating short arc data of space debris based on factor graphs. Background Technology

[0002] With the rapid development of near-Earth space activities and constellation deployments, the number of LEO (Low Earth Orbit) space debris continues to increase, posing a collision risk to spacecraft in orbit. Optical angle measurement observations are widely used for debris monitoring due to their passive nature, wide coverage, and low cost. However, due to limitations such as observation windows, target brightness, weather, and platform geometry, the actual measurements obtained are mostly short arcs (such as right ascension RA, declination DEC, or azimuth Az and elevation E1) with very short time spans (usually less than 1 minute). The short arcs have insufficient information, high noise content, and weak geometric constraints, which limits the accuracy and efficiency of traditional correlation and IOD (Initial Orbit Determination) methods in complex scenarios.

[0003] In existing technologies, methods such as Multiple Hypothesis Tracking (MHT), Joint Data Probability (JPDA), Kalman filtering, least squares, and RANSAC have been used for association and estimation. However, facing uncertainties under short-arc conditions, the explosion of association hypothesis combinations, scene differences (such as atmospheric refraction and Earth's rotation on ground, and the relative motion and attitude effects of space-based platforms), and real-time requirements, existing methods suffer from insufficient robustness and scalability. Existing technologies such as CN111457916B focus on random finite set modeling, and US7105791B1 optimizes from the perspective of illumination observation, but have not yet systematically solved the LEO short-arc data association problem within a unified factor graph framework, especially lacking general modeling and efficient optimization methods that simultaneously consider ground / space-based and mixed scenes.

[0004] Researchers have conducted extensive work and achieved several key advancements in areas such as short-arc association, angle-only observation initial trajectory determination, TLE / SGP4 propagation, and factor graph optimization. Cai et al. proposed an improved tracklet association method for short-arc optical measurements, significantly improving the accuracy and robustness of short-arc association by refining the loss function and initial value search strategy in the short-arc case. Vallado et al. systematically compiled and revised the SGP4 / SDP4 implementation specifications, demonstrating the feasibility and limitations of SGP4 with NORAD / TLE as input in engineering applications and situational awareness comparison, laying an engineering foundation for observation-catalog matching based on catalog comparison. Gooding's angle-only initial trajectory solution program based on the Lambert problem has been widely adopted and has been extended and improved in subsequent work to adapt to multi-observation point and short-arc conditions, verifying that an engineering-meaning initial trajectory solution can be obtained under angle-only observation conditions. Dellaert et al. systematically described and promoted the application of factor graphs and incremental solvers (such as iSAM2 / GTSAM) in large-scale sparse Bayesian estimation, showing that the factor graph framework has good modeling and expressive capabilities and computational efficiency in joint data association and trajectory / state estimation problems, thus providing a mature methodology and toolchain for introducing factor graphs into short-arc data association.

[0005] In summary, existing literature has yielded relatively systematic research results in areas such as short-arc observation correlation, angle-only initial orbit determination, TLE / SGP4 orbit propagation modeling, and factor graph sparsity optimization. However, under conditions of limited observational information (short arcs, angles only measured), how to uniformly incorporate heterogeneous observational data from multiple platforms (ground-based and space-based) into the factor graph framework after time synchronization and coordinate transformation to achieve highly consistent and high-confidence arc segment correlation and joint orbit estimation still presents challenges in terms of engineering implementation and algorithm coordination. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the purpose of this invention is to provide a method and device for associating short arc data of space debris based on factor graphs, which can perform efficient and robust multi-arc segment association and preliminary orbit estimation on optical angle measurement short arc data in ground-based and space-based observation scenarios, and achieve scene adaptation, real-time performance and high accuracy.

[0007] To achieve the above objectives, this invention provides a method for associating spatial fragment short arc data based on factor graphs, comprising the following steps:

[0008] S1. Acquire optical angle measurement short arc data from the station and perform preprocessing including time unification, coordinate transformation and scene correction. Output standardized short arc observation data and observation noise covariance. The station is a ground-based or space-based observation platform.

[0009] S2. For each short arc observation data, the angle observation initial orbit determination algorithm is used to obtain the initial value of the orbit state;

[0010] S3. Construct a cost matrix based on Mahalanobis distance, use the optimal matching strategy based on Hungarian algorithm to preliminarily screen the initial values ​​of the orbital state corresponding to different short arc observation data, and generate a limited number of multi-association hypotheses through perturbation sampling, with each association hypothesis corresponding to a candidate cluster.

[0011] S4. Under each correlation hypothesis, construct a nonlinear factor graph that includes orbital state variables, observation factors, dynamic factors, correlation factors, and prior factors;

[0012] S5. For nonlinear factor graphs, the Gauss-Newton method, Levenberg-Marquardt algorithm or iSAM2 algorithm are used to solve the maximum a posteriori estimate, and the optimal correlation hypothesis is selected based on the residual and information criteria.

[0013] S6 outputs the associated clusters, orbital states and their uncertainties, and performs a consistency check.

[0014] As a preferred embodiment of the present invention, in S1, the non-cooperative target observed by the station is a low-Earth orbit space target, and it is assumed that the low-Earth orbit space target maneuvers no more than once during the observation period; the observation arc is a short arc, and the observation data is only angle observation data.

[0015] As a preferred embodiment of the present invention, the preprocessing process in S1 is as follows:

[0016] S1.1 Time unification: The timestamps of all optical angle measurement short arc data are unified to Coordinated Universal Time (UTC).

[0017] S1.2 Coordinate transformation: unify the observation coordinates of the station to the geocentric inertial coordinate system, the standard form of which is the J2000.0 coordinate system;

[0018] S1.3 Scene Correction: The ground-based observation platform is a ground-based scene, and its correction process includes:

[0019] Atmospheric refraction correction uses an atmospheric refraction model to compensate for elevation angle observations;

[0020] Earth rotation and polar motion correction: By introducing the Earth rotation angular velocity parameter and polar motion correction amount, the station attitude at the time of observation is dynamically adjusted.

[0021] Station location and meteorological parameter compensation: Combining the station's precise geodetic coordinates with real-time meteorological data during observation, the system corrects the influence of meteorological conditions on atmospheric refractive index and compensates for the interference of station location errors on coordinate transformation. The precise geodetic coordinates include the station's latitude, longitude, and elevation, while the real-time meteorological data during observation includes temperature, air pressure, and humidity.

[0022] The space-based observation platform is a space-based scenario, and its calibration process includes:

[0023] Platform attitude compensation uses real-time attitude measurement data from star sensors and gyroscopes in the space-based observation platform to correct the satellite attitude at the time of observation.

[0024] Exposure delay and rolling shutter effect correction: By using pre-calibrated camera parameters, the timestamp and angle of each observation data are corrected to eliminate imaging timing errors;

[0025] Camera coordinate system error correction involves obtaining distortion parameters through camera calibration and correcting distortion in the original angle observation data.

[0026] Short-term geometric correction under J2 perturbation: A simplified model of J2 perturbation is introduced to compensate and correct the target orbit geometric position during the short-arc observation period for the short-term orbit perturbation of low Earth orbit targets caused by the Earth's equatorial uplift effect.

[0027] S1.4 After preprocessing, the output is standardized short-arc observation data and observation noise covariance that unify time and coordinate system.

[0028] As a preferred embodiment of the present invention, in S2, for each short arc observation data, when the station is a ground-based observation platform, the initial value of the orbital state in the geocentric inertial coordinate system is obtained by using the Gauss three-point method or the Gooding angle observation method; when the station is a space-based observation platform, the initial value of the orbital state is obtained by using the extended Gooding algorithm or the relative orbit approximation based on the Clohessy-Wiltshire equation.

[0029] After obtaining the initial values ​​of the orbital state, the presence of outlier observations is identified by calculating the observation residuals. When outlier observations are found, outliers are reduced in weight or removed using robust estimation methods.

[0030] As a preferred embodiment of the present invention, the specific implementation of step S3 is as follows:

[0031] S3.1 Map the initial orbital state values ​​to the observation angular domain. After mapping, the initial orbital state values ​​corresponding to different short arc observation data are the short arcs to be associated. For all short arcs to be associated, perform the following operations:

[0032] S3.1.1 Combine the short arcs to be associated in pairs to form all possible candidate association hypotheses for short arc pairs, and obtain each pair of candidate short arcs;

[0033] S3.1.2 For each pair of candidate short arcs, based on their initial orbital state values ​​and the corresponding observation noise covariance matrix, calculate the Mahalanobis distance between the short arc and the cluster center in the observation angular domain. Here, a cluster is a group of short arcs belonging to the same spatial target, and the cluster center is the statistical or logical center of the initial orbital state values ​​of this group of short arcs in the observation angular domain.

[0034] S3.1.3 Organize the Mahalanobis distances of each pair of candidate short arcs according to the dimension of "short arc number × short arc number" to obtain the cost matrix. The element values ​​of the matrix are the association costs of each pair of short arcs.

[0035] S3.2 Input the cost matrix into the Hungarian algorithm to select the initial optimal set of association hypotheses;

[0036] S3.3. Based on the screening results, and in combination with the noise distribution and geometric constraints, a limited number of multi-association hypotheses are generated using perturbation sampling. Each association hypothesis corresponds to a candidate cluster, which is a set of short arcs that are potential targets in the same space.

[0037] As a preferred embodiment of the present invention, in S4, the nonlinear factor graph specifically includes:

[0038] Orbital state variables, including orbital state variables for each candidate cluster;

[0039] The observation factors are determined by the top-center projection model when the station is a ground-based observation platform and the platform attitude-camera geometric model when the station is a space-based observation platform.

[0040] The dynamic factor adopts a two-body orbital dynamics model to propagate the orbital state from the previous moment to the next observation moment, and is used to construct a consistency constraint for the state changing over time. The process noise includes the influence of J2 perturbation.

[0041] The correlation factor imposes constraints on geometric consistency and temporal sparsity, and employs Huber or Cauchy kernel functions to enhance robustness;

[0042] Prior factors impose loose constraints on orbital states to avoid unconstrained drift.

[0043] As a preferred embodiment of the present invention, in S5, the Gauss-Newton method, the Levenberg-Marquardt algorithm, or the iSAM2 algorithm based on incremental smoothing mapping is used to perform nonlinear least squares optimization on each factor in the nonlinear factor graph to solve for the maximum a posteriori estimate. The normalized residual, the negative log-likelihood of the posterior, and the information criterion are used to score and select multiple hypotheses to obtain the optimal association hypothesis. The information criterion adopts the Akaike information criterion or the Bayesian information criterion.

[0044] As a preferred embodiment of the present invention, in S6, the correlation cluster is the candidate cluster corresponding to the optimal correlation hypothesis, the orbital state is the maximum a posteriori estimate after solving the nonlinear factor graph, and the uncertainty is represented by the covariance matrix of the orbital state.

[0045] As a preferred embodiment of the present invention, in step S6, consistency checks are performed by propagation reprojection and multi-platform cross-validation. If the check passes, the association and estimation results are considered valid; otherwise, steps S3-S5 are re-executed.

[0046] The propagation reprojection verification is as follows: based on the orbital state of the associated cluster, the theoretical orbital state at the observation time is calculated through the orbital propagation model, and then the theoretical observation value is calculated back using the observation angle domain mapping. This is compared with the actual observation value of the standardized short arc observation data, and the normalized residual is calculated. If the residual is less than the threshold corresponding to the observation noise covariance, the consistency of the orbital state of the associated cluster in the observation angle domain is verified, and the test is passed.

[0047] The multiple platforms are ground-to-ground, space-to-space, or ground-to-space hybrid observation platforms. The multi-platform cross-validation is as follows: for the multi-platform observation scenario, firstly, the spatiotemporal reference of each platform is unified, then the orbital state of different platforms for the same associated cluster is extracted, and their position and velocity deviations are calculated. If the deviation is less than the set joint threshold, the consistency of the orbital state of the associated cluster among the multiple platforms is verified, and the test is passed.

[0048] A device for associating short arc spatial fragment data based on factor graphs includes a memory, a processor, and a computer program stored in the memory and executable on the processor. The processor executes the program to implement the above-mentioned method.

[0049] The beneficial effects of this invention are:

[0050] This invention models the short-arc data association problem as a joint probabilistic optimization problem, organically integrating observational and dynamic constraints within a factor graph framework. It combines robust loss functions and optimization algorithms such as Levenberg-Marquardt and iSAM2 to effectively overcome the pain points of insufficient information, high noise content, and weak geometric constraints in short-arc scenarios. For ground-based scenarios involving atmospheric refraction and Earth rotation correction, space-based scenarios involving platform attitude and camera coordinate transformation, and mixed scenarios involving cross-platform time synchronization and coordinate fusion, it achieves seamless adaptation across multiple scenarios through customized observation factors and cross-scenario constraint factors. In short-arc noise environments, the association accuracy can reach over 95%, with orbital position errors remaining stable within the kilometer range, significantly exceeding the accuracy limit of traditional methods in complex scenarios. Furthermore, by generating association hypotheses and residuals through perturbation sampling and selecting the optimal solution using the BIC criterion, it significantly improves the ability to suppress false associations.

[0051] This invention employs an iSAM2 incremental optimization and parallel hypothesis screening strategy, combined with highly efficient computational performance that processes hundreds of arc segments in less than 10 seconds, fully meeting the real-time requirements of online Space Situational Awareness (SSA) tasks. Relying on mature open-source libraries such as GTSAM and ERFA, it integrates engineering models such as SGP4 orbit propagation and Gauss / Gooding initial orbit determination, simplifying development and maintenance processes and lowering deployment barriers. Furthermore, through propagation reprojection back-calculation and multi-platform cross-validation in the post-processing stage, it further ensures the consistency and reliability of the correlation cluster and orbit state estimation results. This not only adapts to single ground-based or space-based observation scenarios but also efficiently handles the fusion of multi-source heterogeneous data in mixed scenarios, providing a high-precision, highly adaptable, and engineering-practical solution for monitoring and determining the orbits of low Earth orbit space debris. Attached Figure Description

[0052] Figure 1 This is a flowchart illustrating the principle of this invention;

[0053] Figure 2 It is a flowchart of the top-center projection and data correction from the orbital state to the observation angle in a ground-based observation scenario;

[0054] Figure 3 This is a technical flowchart illustrating the relationship between orbital state and observation angle and attitude in a space-based observation scenario.

[0055] Figure 4 This is a technical flowchart of spatiotemporal unification and observation fusion in a ground-to-space cross-platform scenario. Detailed Implementation

[0056] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0057] Example 1: This example mainly addresses the following scenario: Suppose that a low-Earth orbit space debris or a non-cooperative target is observed by a ground-based optical telescope (ground-based observation platform) or a space-based optical sensor (space-based observation platform) within a short observation window, and multiple short arc angle data segments are observed. Due to the short observation time span and insufficient geometric constraints, traditional arc segment correlation methods are difficult to obtain stable and reliable results under these conditions.

[0058] like Figure 1 As shown, a method for associating spatial fragment short arc data based on factor graphs includes the following steps:

[0059] S1. Acquire optical angle measurement short arc data from the station and perform preprocessing including time unification, coordinate transformation and scene correction. Output standardized short arc observation data and observation noise covariance. The station is a ground-based or space-based observation platform.

[0060] S2. For each short arc observation data, the angle observation initial orbit determination algorithm is used to obtain the initial value of the orbit state;

[0061] S3. Construct a cost matrix based on Mahalanobis distance, use the optimal matching strategy based on Hungarian algorithm to preliminarily screen the initial values ​​of the orbital state corresponding to different short arc observation data, and generate a limited number of multi-association hypotheses through perturbation sampling, with each association hypothesis corresponding to a candidate cluster.

[0062] S4. Under each correlation hypothesis, construct a nonlinear factor graph that includes orbital state variables, observation factors, dynamic factors, correlation factors, and prior factors;

[0063] S5. For nonlinear factor graphs, the Gauss-Newton method, Levenberg-Marquardt algorithm or iSAM2 algorithm are used to solve the maximum a posteriori estimate, and the optimal correlation hypothesis is selected based on the residual and information criteria.

[0064] S6 outputs the associated clusters, orbital states and their uncertainties, and performs a consistency check.

[0065] In S1, the non-cooperative target observed by the station is a low-Earth orbit space target. It is assumed that the low-Earth orbit space target will not maneuver more than once during the observation period. The observation arc is a short arc, and the observation data is only angle observation data.

[0066] The preprocessing process is as follows:

[0067] S1.1 Time unification: The timestamps of all optical angle measurement short arc data are unified to Coordinated Universal Time (UTC).

[0068] S1.2 Coordinate transformation: unify the observation coordinates of the station to the geocentric inertial coordinate system (ECI), the standard form of which is the J2000.0 coordinate system;

[0069] S1.3 Scene Correction: The ground-based observation platform is a ground-based scene, and its correction process includes:

[0070] Atmospheric refraction correction uses an atmospheric refraction model (such as the standard atmospheric model) to compensate for elevation angle observations. The atmosphere bends the observed light, causing a deviation between the measured elevation angle and the true elevation angle. The model calculates the refraction angle under different altitudes and meteorological conditions to correct the original elevation angle data and eliminate this systematic error.

[0071] Earth rotation and polar motion correction: By introducing the Earth rotation angular velocity parameter and polar motion correction amount, the station attitude at the observation time is dynamically adjusted. The station has an angular velocity due to the Earth's rotation, and the Earth's polar motion will cause the attitude of the station coordinate system relative to the ECI coordinate system to change. By introducing the Earth rotation angular velocity parameter and polar motion correction amount, the station attitude at the observation time is dynamically adjusted to ensure the accuracy of coordinate system transformation.

[0072] Station location and meteorological parameter compensation: Combining the station's precise geodetic coordinates with real-time meteorological data during observation, the system corrects the influence of meteorological conditions on atmospheric refractive index and compensates for the interference of station location errors on coordinate transformation. The precise geodetic coordinates include the station's latitude, longitude, and elevation, while the real-time meteorological data during observation includes temperature, air pressure, and humidity.

[0073] The space-based observation platform is a space-based scenario, and its calibration process includes:

[0074] Platform attitude compensation: When the space-based observation platform is in orbit, there may be attitude jitter (such as orbital perturbation, attitude control system error). The real-time attitude measurement data of the star sensor and gyroscope in the space-based observation platform are used to correct the satellite attitude at the time of observation.

[0075] Exposure delay and rolling shutter effect correction: Satellite cameras have an exposure delay (the time difference from triggering exposure to imaging), and rolling shutter imaging can cause slight differences in the imaging time of different pixels in the same frame. By using pre-calibrated camera parameters (exposure delay, shutter rolling speed), the timestamp and angle of each observation data are corrected to eliminate imaging timing errors.

[0076] Camera coordinate system error correction addresses distortions in the camera's optical system (such as radial and tangential distortion). Distortion parameters are obtained through camera calibration, and the original angle observation data are corrected to ensure that the angle measurements in the camera coordinate system accurately reflect the target's line of sight.

[0077] Short-term geometric correction under J2 perturbation: A simplified J2 perturbation model is introduced to compensate for the short-term orbital perturbation of low Earth (LEO) orbit targets caused by the Earth's equatorial bulge effect, correcting the target's orbital geometric position within the short-arc observation period. Although the short-arc observation time span is small (usually less than 1 minute), J2 perturbation may still cause slight deviations in the target position calculation. By introducing a simplified J2 perturbation model to compensate for the target's orbital geometric position within the short-arc observation period, the accuracy of the initial value of subsequent initial orbit estimation is improved.

[0078] S1.4 After preprocessing, the output is standardized short-arc observation data with unified time and coordinate system, as well as observation noise covariance (a mathematical matrix describing the characteristics of random error distribution in short-arc observation data).

[0079] The J2 perturbation simplified model is a simplified calculation model for orbital perturbations in the J2 term (Earth equatorial bulge effect) in the non-spherical gravitational field of the Earth. It is used to quickly compensate for short-term orbital deviations of low Earth orbit targets caused by the Earth's oblateness.

[0080] In S2, for each short arc observation data, when the station is a ground-based observation platform, the Gauss three-point method or the Gooding angle observation method is used to obtain the initial value of the orbital state in the geocentric inertial coordinate system; when the station is a space-based observation platform, considering the influence of relative motion, the extended Gooding algorithm or the relative orbit approximation based on the Clohessy-Wiltshire equation is used to obtain the initial value of the orbital state; the space target catalog uses the TLE database published by NORAD, and the theoretical orbital state of the target at the observation time is obtained by propagation using the SGP4 model;

[0081] For example, when the station is a ground-based observation platform, the Gooding algorithm (an initial orbit determination method suitable for high-precision angle observations of short arcs) is used to obtain the initial orbital state values. This algorithm is based on spherical trigonometry and the theory of orbital conic sections. It takes the pre-processed azimuth and elevation sequences from the ground-based station (requiring multiple sets of observations with reasonable time distribution within the short arc) as input, combines them with the station's compensated and accurate geocentric coordinates, constructs the angle observation equations, and solves for the six orbital roots (semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of perihelion, and mean perihelion), ultimately converting them into initial orbital state values ​​(position vector and velocity vector) in the J2000.0 coordinate system.

[0082] When the station is a space-based observation platform, an improved initial orbit algorithm combining the Clohessy-Wiltshire (CW) relative orbit model is used to obtain the initial orbit state. This algorithm utilizes the precise orbit parameters of the space-based platform itself (injected from onboard equipment or ground-based systems) to convert pre-processed space-based angle observations into the line-of-sight vector of the target relative to the space-based platform. Then, using the CW relative orbit model (which describes the relative motion characteristics of low Earth orbit targets) and angle observations, a system of equations is constructed to solve for the target's absolute orbital root numbers, ultimately converting them into initial orbital state values ​​(position and velocity vectors) in the J2000.0 coordinate system. During this process, it is necessary to fully incorporate the short-time orbital geometry characteristics corrected by J2 perturbations to adapt to the motion patterns of low Earth orbit targets in a space-based scenario and improve the accuracy of the initial orbit solution.

[0083] After obtaining the initial values ​​of the orbital state, the presence of outlier observations is identified by calculating the observation residuals. When outlier observations are found, outliers are reduced in weight or removed using robust estimation methods.

[0084] Outlier observations are observations whose angular observation values ​​(such as azimuth / elevation angles on ground and camera field of view angles on space) significantly deviate from the normal statistical distribution or physical expectations due to equipment failure, environmental interference, model errors, etc. Each raw optical angle measurement data point is compared with the initial orbital state value, and the observation residuals are calculated. If the residual of an observation significantly exceeds the normal statistical range (e.g., more than 3 standard deviations), it is identified as an outlier observation. Robust estimation methods (such as robust loss functions and residual threshold screening) are then used to reduce its weight or remove it to avoid interfering with subsequent factor graph optimization and data association processes.

[0085] The specific implementation method of step S3 is as follows:

[0086] S3.1 Map the initial orbital state values ​​to the observation angular domain. After mapping, the initial orbital state values ​​corresponding to different short arc observation data are the short arcs to be associated.

[0087] The mapping can be performed using a scene projection model, which transforms the orbital state from the position-velocity vector of the geocentric inertial frame to the observation angular domain of the observation platform. Its core function is to convert the inertial frame state of orbital dynamics into an angular observation form that directly matches optical angle measurement, making the calculation of Mahalanobis distance more consistent with the physical characteristics of optical observation.

[0088] For all short arcs to be associated, perform the following operations:

[0089] S3.1.1 Combine the short arcs to be associated in pairs to form all possible candidate association hypotheses for short arc pairs, and obtain each pair of candidate short arcs;

[0090] S3.1.2 For each pair of candidate short arcs, based on their initial orbital state values ​​and the corresponding observation noise covariance matrix, calculate the Mahalanobis distance between the short arc and the cluster center in the observation angular domain. Here, a cluster is a group of short arcs belonging to the same spatial target, and the cluster center is the statistical or logical center of the initial orbital state values ​​of this group of short arcs in the observation angular domain (for example, it can be the mean orbit of the initial orbital state values ​​of multiple short arcs in the same cluster, or the core orbital state obtained through preliminary clustering).

[0091] S3.1.3 Organize the Mahalanobis distance of each pair of candidate short arcs according to the dimension of "short arc number × short arc number" to obtain the cost matrix. The element values ​​of the matrix are the association costs of each pair of short arcs (reflecting similarity).

[0092] S3.2 Input the cost matrix into the Hungarian algorithm to select the initial optimal set of association hypotheses;

[0093] S3.3. Based on the screening results, and in combination with the noise distribution and geometric constraints (based on the orbital dynamics of space debris and the physical constraints of the observation platform's geometric layout), a limited number of multi-association hypotheses are generated using perturbation sampling. Each association hypothesis corresponds to a candidate cluster (the total number of association hypotheses is limited to within 1000). The candidate cluster is a set of potential short arcs belonging to the same space target.

[0094] After preprocessing and initial orbital state calculations, each short arc data point generates multiple possible arc-target association hypotheses. A cost matrix is ​​constructed using Mahalanobis distance, and then the optimal matching strategy of the Hungarian algorithm is employed to filter the short arc pairs most likely belonging to the same space debris target from these hypotheses. This initial arc association screening provides a candidate set for precise association in the subsequent factor graph. Mahalanobis distance characterizes the statistical consistency between short arcs, while the Hungarian algorithm achieves globally optimal matching, automatically adjusting the threshold and re-searching when a match fails.

[0095] In S4, the factor graph consists of orbital state variable nodes and observations, dynamics, correlations, and prior factors. The observation model is adaptively selected based on the ground-based or space-based scenario. The nonlinear factor graph specifically includes:

[0096] Orbital state variables, including orbital state variables for each candidate cluster;

[0097] For observation factors, the top-center projection (station-center projection) model is used when the station is a ground-based observation platform, and the platform attitude-camera geometry model is used when the station is a space-based observation platform.

[0098] The zenith direction of the station is used as the reference core. The local zenith-north-east horizon coordinate system of the station is used as the core. The inertial orbit state of the target is projected into the observable azimuth (Az) / elevation (El) of the station. It is a special conversion model from orbit state to observation angle in the ground-based optical angle measurement scenario.

[0099] The platform attitude-camera geometry model relies on the attitude of the space-based platform (the relative orientation of the satellite body coordinate system and the inertial system) and the camera's intrinsic / extrinsic parameters (geometric imaging relationship). It is a dedicated mapping model from the orbital state to the camera's pixel angle / image plane coordinates in space-based scenarios.

[0100] The dynamics factor employs a two-body orbital dynamics model (considering only the gravity of the central celestial body) to propagate the orbital state (position and velocity) from the previous moment to the next observation moment, thereby constructing consistency constraints on the state's change over time. The process noise includes the influence of J2 perturbation (by introducing process noise to approximate the influence of unmodeled perturbations such as J2 perturbation); the consistency constraints ensure that the orbital state (position and velocity vector) of space debris at adjacent moments conforms to the classical two-body motion laws.

[0101] The correlation factor imposes constraints on geometric consistency (short arcs under the same correlation hypothesis must have consistent orbit / observation geometry) and temporal sparsity (short arcs within the same candidate cluster must have a reasonable distribution of observation time, neither too dense nor too large intervals), eliminates physically impossible short arc correlations, and uses Huber or Cauchy kernel functions (normal residuals are retained, and abnormal residuals are weighted down) to enhance robustness.

[0102] Prior factors impose loose constraints on orbital states to avoid unconstrained drift.

[0103] As a probabilistic graphical model, factor graphs can decompose complex estimation problems into variable nodes and factor nodes, and efficiently handle uncertainties and sparse structures through message passing or MAP optimization (such as Gauss-Newton, Levenberg-Marquardt, and iSAM2). Introducing factor graphs into short-arc correlations of spatial debris, and customizing observation and dynamic factors, constraints, and robust losses for different observation scenarios, can improve accuracy, efficiency, and versatility.

[0104] In S5, the Gauss-Newton method, Levenberg-Marquardt algorithm, or iSAM2 algorithm based on incremental smoothing mapping is used to perform nonlinear least squares optimization on each factor in the nonlinear factor graph to solve for the maximum a posteriori estimate. The normalized residual, negative log-likelihood of the posterior, and information criterion are used to score and select multiple hypotheses to obtain the optimal association hypothesis. The information criterion adopts the Akaike information criterion or the Bayesian information criterion.

[0105] In S6, the correlation cluster is the candidate cluster corresponding to the optimal correlation hypothesis, the orbital state is the maximum a posteriori estimate after solving the nonlinear factor graph, the error range is generally better than 1km, and the uncertainty is represented by the covariance matrix of the orbital state.

[0106] Consistency checks are performed using propagation reprojection and multi-platform cross-validation. If the checks pass, the association and estimation results are considered valid; otherwise, steps S3-S5 are repeated.

[0107] The propagation reprojection verification is as follows: based on the orbital state of the associated cluster, the theoretical orbital state at the observation time is calculated through the orbital propagation model, and then the theoretical observation value is calculated back using the observation angle domain mapping. This is compared with the actual observation value of the standardized short arc observation data, and the normalized residual is calculated. If the residual is less than the threshold corresponding to the observation noise covariance, the consistency of the orbital state of the associated cluster in the observation angle domain is verified, and the test is passed.

[0108] The multiple platforms are ground-to-ground, space-to-space, or ground-to-space hybrid observation platforms. The multi-platform cross-validation is as follows: for the multi-platform observation scenario, firstly, the spatiotemporal reference of each platform is unified, then the orbital state of different platforms for the same associated cluster is extracted, and their position and velocity deviations are calculated. If the deviation is less than the set joint threshold, the consistency of the orbital state of the associated cluster among the multiple platforms is verified, and the test is passed.

[0109] Verification Example 1: Observation and Correlation in Ground-Based Scenarios:

[0110] The implementation environment included a ground-based optical telescope system for acquiring angular observation data of the target. The observation data was acquired from the same station within the same observation window, resulting in three short arc observation segments, each lasting approximately 10 seconds, with an angular observation accuracy of approximately 1 arcsecond.

[0111] First, the raw observation data (optical angle measurement short arc data) is preprocessed, and an atmospheric refraction model is used to correct the refraction of elevation angle observations to eliminate systematic errors caused by atmospheric refraction. Second, using the angle observation data, the Gaussian initial orbit determination method (or an existing improved Gaussian method) is employed to calculate the initial orbit parameters corresponding to each short arc, obtaining preliminary orbit estimation results (initial orbit state values). Then, correlation hypotheses are generated between the short arcs. The establishment of correlation hypotheses adopts an optimal matching strategy based on the Hungarian algorithm, with its cost matrix constructed based on the projection distance of each arc segment in the top-center projection coordinate system, thereby achieving minimum distance matching between observed arc segments. Based on this, a factor graph model for the ground-based observation scenario is constructed. The observation factors are established based on the top-center projection geometry of the fixed station, and the dynamic factors use a two-body orbit dynamics model to constrain the orbit state at adjacent times. Through the above factor graph structure, observation constraints and dynamic smoothing constraints can be considered simultaneously, thereby achieving joint optimization of observed arc segments and orbit state estimation.

[0112] The factor graph was solved using the factor graph optimization library GTSAM, employing an incremental Gaussian-Newton method for nonlinear least squares optimization. After several iterations (approximately 50), the optimization converged, yielding the optimal orbital state estimate (maximum a posteriori estimate) of the target.

[0113] The final output includes arc segment association clusters and orbit estimates. Calculation results show that the association accuracy is significantly improved after processing with the method described in this embodiment, and the obtained orbit position error is better than 1 km, verifying the effectiveness and stability of this method in ground-based short-arc observation scenarios.

[0114] Figure 2 The process of center projection and data correction from the orbital state to the observation angle in the ground-based observation scenario of Verification Example 1 is presented intuitively. Figure 2With Earth as the core reference, Earth's rotation is a key influencing factor, causing the ground-based observation station to move with the Earth, thus dynamically changing the station center attitude in ECEF (Earth-Fixed Coordinate System) over time. ECI→topo(t) represents the coordinate transformation matrix at time t, used to transform the orbital state of the observed target in the J2000.0 coordinate system (ECI) to the station's local station center horizontal coordinate system (topo). The ground-based observation station (ground-based observation platform) is the main observation subject, using the top-center projection mode (AZ-EI, a dedicated projection model with the station's zenith direction as the core reference), capturing the observation direction of the observed target through the observation optical path. In the process, the original elevation angle data is first input, and elevation angle correction is achieved through atmospheric refraction angle calculation (compensating for observation deviations caused by atmospheric refraction). Then, the coordinates are unified (aligned with UTC time and ECI coordinate system reference) by combining the coordinate transformation matrix and observation equations. Finally, the core results of ground-based observation are output: azimuth AZ and elevation EI. The time reference for the entire process is calibrated by time t (the unified Coordinated Universal Time).

[0115] Verification Example 2: Observation and Correlation in a Space-Based Scenario:

[0116] The observation platform is an in-orbit satellite optical sensor system used to acquire angular observation information of the target. A total of 5 short arc observation data were collected, each arc lasting approximately 20 seconds, with an angle measurement accuracy better than 0.5 arcseconds.

[0117] First, the raw observation data is preprocessed. Since the space-based platform has its own attitude changes, the satellite attitude at the time of observation needs to be corrected using the measurement data of the star sensor and gyroscope, so as to transform the direction of the observation optical axis to the inertial coordinate system (geocentric inertial coordinate system, the same below).

[0118] Secondly, based on the angle observation data after attitude correction, an improved angle-only initial trajectory determination method (extended Gooding algorithm) is used to estimate the initial trajectory for each short arc. This method is an extension of the classic Gooding algorithm, which improves the accuracy of initial trajectory estimation when the angle observation noise is low by introducing weighted constraints and iterative optimization processes under short arc conditions.

[0119] Then, correlation hypotheses between short arcs are generated. In the space-based observation scenario, since both the platform and the target are in orbital motion, a relative velocity constraint factor can be introduced into the correlation cost matrix to comprehensively consider the relative motion direction and line-of-sight velocity information between each arc segment, thereby improving the accuracy of the correlation hypotheses.

[0120] Based on this, a factor graph model for space-based observation scenarios is constructed. Observation factors are established based on the satellite platform coordinate system and camera projection geometry, while dynamic factors employ a two-body orbital dynamics model to describe the target's orbital evolution across continuous arc segments. This factor graph structure allows for the unification of multiple short-arc observation information and dynamic constraints within a single optimization framework for joint solution.

[0121] The factor graph solution process employs a factor graph optimization method based on the incremental smooth mapping algorithm (iSAM2), which uses incremental Gaussian-Newton iterations to achieve real-time updates of the orbital state and associated variables. After several iterations, the optimal orbital state estimate and arc segment association results converge.

[0122] The results show that, under the conditions of this validation example, the accuracy of arc segment association reaches approximately 98%, and the orbit estimation accuracy is significantly improved, verifying the efficiency and robustness of the method in this embodiment under the space-based short arc observation scenario.

[0123] Figure 3 This paper fully presents the technical process from orbital state to observation angle and attitude in the space-based observation scenario of Verification Example 2. With the Earth as the core, a geocentric reference frame is provided as the global coordinate reference. This reference frame extends downwards to the inertial coordinate system, achieving a unified spatial reference. The space-based observation platform determines its position in the inertial coordinate system and orbits the Earth. The platform constructs a camera coordinate system by defining the camera attitude and uses attitude control and measurement to perform real-time attitude updates, ensuring camera pointing accuracy. The camera coordinate system forms the observation optical path along the observation direction, pointing towards the observation target (such as space debris). The observation target itself has orbital motion, and its target orbit needs to be transformed from orbital coordinates to ECI to be unified into the inertial coordinate system. Subsequently, the observation equation is constructed through the measurement modeling stage to quantify the spatial correlation between the observation optical path and the observation target, ultimately outputting the observation angle, attitude, and time label.

[0124] Verification Example 3: Observation and Correlation in Mixed Scenes:

[0125] In this validation example, the data sources include observations from ground-based optical telescopes and space-based optical sensors. Ground-based observations provide high-time-accuracy angle measurements, while space-based observations provide multi-view spatial geometric constraints. Both types of data consist of short-arc observation segments with time spans of approximately 10 seconds and 20 seconds, respectively.

[0126] To achieve unified processing of ground-based and space-based observations, cross-scene time synchronization factors and coordinate transformation factors are introduced into the factor graph framework.

[0127] The time synchronization factor is used to correct the system time deviation between different observation platforms and unify the ground-based and space-based observation times to a unified time scale (such as UTC); the coordinate transformation factor is used to describe the geometric transformation relationship between the Earth-Fixed Coordinate System (ECEF) and the satellite camera coordinate system, so as to realize the unified expression of cross-platform observation data in the inertial frame (ECI / J2000.0).

[0128] In the factor graph construction process, the ground-based observation factors adopt the top-center projection model, the space-based observation factors adopt the camera line-of-sight projection model (platform attitude-camera geometry model), and the dynamic factors use two-body motion equations to constrain adjacent state nodes. By introducing the time synchronization factor and coordinate transformation factor into the same factor graph model, orbital state, observation-related variables, and cross-scene time and coordinate errors can be estimated simultaneously within a unified optimization framework.

[0129] The optimization process employs the incremental smoothing and mapping algorithm (iSAM2) to achieve joint optimal estimation of observation information from multiple platforms through iterative solutions.

[0130] The results show that the method in this embodiment takes less than 10 seconds to process 100 observation arc segment data, has high consistency in correlation results, and the orbit estimation error is maintained within the kilometer range, verifying the real-time performance and stability of the method in this embodiment under mixed observation scenarios.

[0131] Figure 4 The technical process of spatiotemporal unification and observation fusion in the ground-based-space-based cross-platform scenario of Verification Example 3 is presented. With the Earth as the core, it provides global time reference and geocentric reference, extending downward to the inertial coordinate system as a global spatial reference.

[0132] The process is divided into two links: ground-based and space-based, which are ultimately merged.

[0133] Ground-based observation station link: Ground-based observation stations achieve coordinate unification through R_topo→ECI (transformation matrix from station-centered horizontal coordinate system to inertial coordinate system), and their observation optical path points to the observation target (such as space debris); at the same time, a time synchronization factor Δt is introduced through the time alignment link, and a unified time reference is generated after Δt correction to ensure the consistency of the time dimension.

[0134] Space-based observation platform link: The space-based observation platform is associated with the inertial coordinate system through the station center transformation platform position, and its observation optical path also points to the observation target; and coordinate unification is completed through R_cam→ECI (camera coordinate system to inertial coordinate system transformation matrix) and coordinate transformation factor, and then enters the joint observation stage.

[0135] Finally, by combining the unified time reference with the cross-platform observation fusion factor map optimization input, and through the dual unification of coordinates and time and joint observation processing, a unified spatiotemporal observation set is output.

[0136] Example 2: A spatial fragment short arc data association device based on factor graph, including a memory, a processor, and a computer program stored in the memory and run on the processor, which implements the method in Example 1 by executing the program through the processor.

Claims

1. A method for associating spatial fragment short-arc data based on factor graphs, characterized in that... Includes the following steps: S1. Acquire optical angle measurement short arc data from the station and perform preprocessing including time unification, coordinate transformation and scene correction. Output standardized short arc observation data and observation noise covariance. The station is a ground-based or space-based observation platform. S2. For each short arc observation data, the angle observation initial orbit determination algorithm is used to obtain the initial value of the orbit state; S3. Construct a cost matrix based on Mahalanobis distance, use the optimal matching strategy based on Hungarian algorithm to preliminarily screen the initial values ​​of the orbital state corresponding to different short arc observation data, and generate a limited number of multi-association hypotheses through perturbation sampling, with each association hypothesis corresponding to a candidate cluster. S4. Under each correlation hypothesis, construct a nonlinear factor graph that includes orbital state variables, observation factors, dynamic factors, correlation factors, and prior factors; S5. For nonlinear factor graphs, use the Gauss-Newton method, Levenberg-Marquardt algorithm or iSAM2 algorithm to solve for the maximum a posteriori estimate, and select the optimal correlation hypothesis based on the residual and information criteria. S6. Output the associated clusters, orbital states and their uncertainties, and perform a consistency check; The specific implementation method of step S3 is as follows: S3.1 Map the initial orbital state values ​​to the observation angular domain. After mapping, the initial orbital state values ​​corresponding to different short arc observation data are the short arcs to be associated. For all short arcs to be associated, perform the following operations: S3.1.1 Combine the short arcs to be associated in pairs to form all possible candidate association hypotheses for short arc pairs, and obtain each pair of candidate short arcs; S3.1.2 For each pair of candidate short arcs, based on their initial orbital state values ​​and the corresponding observation noise covariance matrix, calculate the Mahalanobis distance between the short arc and the cluster center in the observation angular domain. Here, a cluster is a group of short arcs belonging to the same spatial target, and the cluster center is the statistical or logical center of the initial orbital state values ​​of this group of short arcs in the observation angular domain. S3.1.3 Organize the Mahalanobis distances of each pair of candidate short arcs according to the dimension of "short arc number × short arc number" to obtain the cost matrix. The element values ​​of the matrix are the association costs of each pair of short arcs. S3.2 Input the cost matrix into the Hungarian algorithm to select the initial optimal set of association hypotheses; S3.

3. Based on the screening results, and in combination with the noise distribution and geometric constraints, a limited number of multi-association hypotheses are generated using perturbation sampling. Each association hypothesis corresponds to a candidate cluster, which is a set of short arcs that are potential targets in the same space. In S4, the nonlinear factor graph specifically includes: Orbital state variables, including orbital state variables for each candidate cluster; The observation factors are determined by the top-center projection model when the station is a ground-based observation platform and the platform attitude-camera geometry model when the station is a space-based observation platform. The dynamic factor adopts a two-body orbital dynamics model to propagate the orbital state from the previous moment to the next observation moment, and is used to construct a consistency constraint for the state changing over time. The process noise includes the influence of J2 perturbation. The correlation factor imposes constraints on geometric consistency and temporal sparsity, and employs Huber or Cauchy kernel functions to enhance robustness; Prior factors impose loose constraints on orbital states to avoid unconstrained drift.

2. The method for associating spatial fragment short arc data based on factor graphs according to claim 1, characterized in that, In S1, the non-cooperative target observed by the station is a low-Earth orbit space target, and it is assumed that the low-Earth orbit space target will not maneuver more than once during the observation period; the observation arc is a short arc, and the observation data is only angle observation data.

3. The method for associating spatial fragment short arc data based on factor graphs according to claim 1, characterized in that, In S1, the preprocessing process is as follows: S1.1 Time unification: The timestamps of all optical angle measurement short arc data are unified to Coordinated Universal Time (UTC). S1.2 Coordinate transformation: unify the observation coordinates of the station to the geocentric inertial coordinate system, the standard form of which is the J2000.0 coordinate system; S1.3 Scene Correction: The ground-based observation platform is a ground-based scene, and its correction process includes: Atmospheric refraction correction uses an atmospheric refraction model to compensate for elevation angle observations; Earth rotation and polar motion correction: By introducing the Earth rotation angular velocity parameter and polar motion correction amount, the station attitude at the time of observation is dynamically adjusted. Station location and meteorological parameter compensation: Combining the station's precise geodetic coordinates with real-time meteorological data during observation, the system corrects the influence of meteorological conditions on atmospheric refractive index and compensates for the interference of station location errors on coordinate transformation. The precise geodetic coordinates include the station's latitude, longitude, and elevation, while the real-time meteorological data during observation includes temperature, air pressure, and humidity. The space-based observation platform is a space-based scenario, and its calibration process includes: Platform attitude compensation uses real-time attitude measurement data from star sensors and gyroscopes in the space-based observation platform to correct the satellite attitude at the time of observation. Exposure delay and rolling shutter effect correction: By using pre-calibrated camera parameters, the timestamp and angle of each observation data are corrected to eliminate imaging timing errors; Camera coordinate system error correction involves obtaining distortion parameters through camera calibration and correcting distortion in the original angle observation data. Short-term geometric correction under J2 perturbation: A simplified model of J2 perturbation is introduced to compensate and correct the target orbit geometric position during the short-arc observation period for the short-term orbit perturbation of low Earth orbit targets caused by the Earth's equatorial uplift effect. S1.4 After preprocessing, the output is standardized short-arc observation data and observation noise covariance that unify time and coordinate system.

4. The method for associating spatial fragment short arc data based on factor graphs according to claim 1, characterized in that, In S2, for each short arc observation data, when the station is a ground-based observation platform, the initial value of the orbital state in the geocentric inertial coordinate system is obtained by using the Gauss three-point method or the Gooding angle observation method; when the station is a space-based observation platform, the initial value of the orbital state is obtained by using the extended Gooding algorithm or the relative orbit approximation based on the Clohessy-Wiltshire equation. After obtaining the initial values ​​of the orbital state, the presence of outlier observations is identified by calculating the observation residuals. When outlier observations are found, outliers are reduced in weight or removed using robust estimation methods.

5. The method for associating spatial fragment short arc data based on factor graphs according to claim 1, characterized in that, In S5, the Gauss-Newton method, Levenberg-Marquardt algorithm, or iSAM2 algorithm based on incremental smoothing mapping is used to perform nonlinear least squares optimization on each factor in the nonlinear factor graph to solve for the maximum a posteriori estimate. The normalized residual, negative log-likelihood of the posterior, and information criterion are used to score and select multiple hypotheses to obtain the optimal association hypothesis. The information criterion adopts the Akaike information criterion or the Bayesian information criterion.

6. The method for associating spatial fragment short arc data based on factor graphs according to claim 1, characterized in that, In S6, the correlation cluster is the candidate cluster corresponding to the optimal correlation hypothesis, the orbital state is the maximum a posteriori estimate after solving the nonlinear factor graph, and the uncertainty is represented by the covariance matrix of the orbital state.

7. The method for associating spatial fragment short arc data based on factor graphs according to claim 6, characterized in that, In step S6, consistency checks are performed through propagation reprojection and multi-platform cross-validation. If the checks pass, the association and estimation results are considered valid; otherwise, steps S3-S5 are repeated. The propagation reprojection verification is as follows: based on the orbital state of the associated cluster, the theoretical orbital state at the observation time is calculated through the orbital propagation model, and then the theoretical observation value is calculated back using the observation angle domain mapping. This is compared with the actual observation value of the standardized short arc observation data, and the normalized residual is calculated. If the residual is less than the threshold corresponding to the observation noise covariance, the consistency of the orbital state of the associated cluster in the observation angle domain is verified, and the test is passed. The multiple platforms are ground-to-ground, space-to-space, or ground-to-space hybrid observation platforms. The multi-platform cross-validation is as follows: for the multi-platform observation scenario, firstly, the spatiotemporal reference of each platform is unified, then the orbital state of different platforms for the same associated cluster is extracted, and their position and velocity deviations are calculated. If the deviation is less than the set joint threshold, the consistency of the orbital state of the associated cluster among the multiple platforms is verified, and the test is passed.

8. A device for associating spatial fragment short-arc data based on factor graphs, characterized in that: It includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method described in any one of claims 1 to 7.

Images