A method for identifying and detecting low-Earth orbit spacecraft maneuvers

Abstract

This invention relates to a method for identifying and detecting low-Earth orbit spacecraft maneuvers. It pertains to the fields of space situational awareness and spacecraft orbital dynamics. By analyzing different types of maneuver strategies and parametric dynamic responses, this invention establishes the quasi-average semi-major axis, reflecting system energy changes, as a core feature parameter. Subsequently, addressing the challenge of environmental noise and maneuver signal aliasing in long-term ephemeris data, a progressively developed dual-sliding denoising detection algorithm, incorporating adaptive thresholding and considering data splicing errors, is constructed, forming a complementary screening and correction mechanism. Finally, based on accurate maneuver event localization, a comprehensive statistical analysis of the space perturbation background and maneuver increments is conducted to extract highly discriminative physical features, laying a reliable data foundation for building intelligent maneuver behavior prediction models.

Description

Technical Field

[0001] This invention relates to the field of space situational awareness and spacecraft orbital dynamics technology, and is a method for identifying and detecting low-Earth orbit spacecraft maneuvers. Background Technology

[0002] With the dense deployment of large low-Earth orbit constellations, orbit maintenance and avoidance maneuvers for space targets during their on-orbit operation are becoming increasingly frequent. Traditional macroscopic forecasting models are ill-equipped to handle such high-frequency, low-thrust orbital maneuvers. Therefore, establishing an efficient and robust maneuver detection and feature extraction system has become an important prerequisite for space situational awareness.

[0003] The existing technology has the following drawbacks: Fixed thresholds have poor universality: The natural decay characteristics of low-Earth orbit satellites vary significantly with orbital altitude, solar activity cycle, and atmospheric density. A single fixed threshold cannot be universally applicable, which can easily lead to a large number of false detections or missed detections. Severe interference from ephemeris splicing errors: When publicly available ephemeris data is released and updated in multiple batches, there are non-physical data jumps at the junctions between adjacent batches. The numerical performance is similar to the semi-major axis step height of a small maneuver, which seriously interferes with the accuracy of detection. Insufficient detection capability for minute maneuvers: Single detection algorithms rely solely on instantaneous energy jumps between adjacent epochs, resulting in weak ability to extract orbital maneuver signals with small energy increments and hidden features; Low feature extraction purity: Existing methods do not perform systematic statistical analysis on the clean orbital data after detection, and cannot provide high-purity, high-confidence feature inputs for subsequent machine learning-based prediction of maneuver intentions. Summary of the Invention

[0004] This invention addresses the problems of poor universality of fixed thresholds, large interference from ephemeris stitching errors, high missed detection rate of minor maneuvers, and low purity of feature extraction in existing technologies. At the same time, it takes into account the computational efficiency of batch processing of giant constellations and establishes a complete data processing link from feature screening to anomaly detection and statistical refinement. Based on this, this invention discloses a method for identifying and detecting low-Earth orbit spacecraft maneuvers.

[0005] This invention provides the following technical solutions: A method for identifying and detecting low-Earth orbit spacecraft maneuvers, the method comprising the following steps: Step 1: Using the quasi-average semi-major axis as the core characteristic parameter to characterize the change in the spacecraft's orbital energy, calculate the non-spherical gravitational perturbation solution based on the Earth's gravitational field model, truncate and select the field harmonic and zone harmonic perturbations, and reconstruct the quasi-average semi-major axis time series; Step 2: Perform a two-step dynamic adaptive threshold detection based on the quasi-average semi-major axis difference sequence to complete the preliminary localization and classification of the maneuvering event; Step 3: Perform splicing error preprocessing on the original pseudo-average semi-major axis sequence, then reconstruct the orbital energy trend baseline by denoising using the double moving average method, and perform refined maneuver detection based on the energy transition signal to make up for the omissions in the preliminary detection. Step 4: Perform statistical analysis based on clean free decay trajectory data, quantify the natural decay baseline and environmental noise boundary, and extract multidimensional high-confidence maneuver feature vectors to provide input for subsequent maneuver prediction models.

[0006] Preferably, the two-step dynamic adaptive threshold detection method specifically comprises: The first-order forward differencing process is performed on the quasi-average semi-major axis time series to obtain the semi-major axis change sequence between adjacent epochs; a conservative fixed initial empirical threshold is introduced, and the sequence is traversed to perform the first round of coarse screening. Epochs with absolute difference values ​​greater than the initial threshold are marked as maneuvering, and those with absolute difference values ​​less than the initial threshold are marked as free decay; the state switching index nodes are recorded to divide the continuous sequence into multiple independent data segments with single attributes. Extract all free decay segments with a length greater than 1 data point, and calculate the average of the absolute values ​​of the differences between adjacent points within each segment; use the number of data points in each free decay segment as the weight to perform a global weighted average of all local average differences to obtain the weighted average decay deviation of the current space environment; multiply the weighted average decay deviation by the dynamic amplification factor as the adaptive maneuver decision threshold; if no valid free decay samples are extracted, the initial empirical threshold is used. Using an adaptive threshold as a benchmark, a second round of refined screening was conducted on the original differential sequence to redefine the free decay segment and the active maneuver segment; linear regression fitting was performed on the free decay segment to calculate the standardized decay rate; for the active maneuver segment, the instantaneous orbit change rate of the pulse maneuver and the average orbit change rate of the continuous thrust maneuver were calculated respectively.

[0007] Preferably, the Earth's gravitational field model adopts the high-precision Earth gravitational field model EGM2008.

[0008] Preferably, the dynamic magnification factor is in the range of 2.0-3.0 times.

[0009] Preferably, the dual sliding noise reduction and splicing error supplementation detection specifically includes: Extract the overlapping time nodes of the ephemeris update handover, calculate the quasi-average semi-major axis difference between the old and new batches of ephemeris at the same epoch, lock the index position corresponding to the splicing error and mark it separately, and eliminate false maneuver signals. Using the spacecraft orbital period as the basic sliding window, the original quasi-average semi-major axis sequence is initially smoothed to filter out residual perturbations at the orbital period level; then, using half an orbital period as the secondary sliding window, the initial results are filtered a second time to extract the baseline of the long-term evolution trend of orbital energy. Set a step window of half an orbital period length, slide it along the time series, and calculate the difference between the average value of the future window data block and the historical window data block at each central evaluation point as the transition energy intensity at that point; Boundary data with one orbital period length at the beginning and end of the sequence are excluded. The 90th percentile of the absolute value of the transition signal in the stable segment is extracted as the environmental noise baseline. Differentiated peak-finding thresholds are set: the orbit-ascending detection threshold is 0.8 times the noise baseline, and the orbit-descending detection threshold is 1.5 times the noise baseline. At the same time, the peak prominence is required to be greater than 0.25 times the baseline, and adjacent maneuvering events must meet the preset time interval constraint to initially identify suspected maneuvering nodes.

[0010] Preferably, the dual-sliding denoising and splicing error supplementary detection also includes a secondary verification process: For each suspected maneuver peak point, a stable reference interval of one step window length is extracted before and after it. The difference between the mean of the semi-major axis of the interval before and after the maneuver is taken as the net maneuver effect. When the absolute value of the net effect exceeds the calculated maximum local average deviation, it is confirmed as a valid trajectory change point.

[0011] Preferably, the comprehensive statistical extraction of maneuver characteristics specifically includes: By segmenting the observation time axis using the splicing node database, continuous data segments that have not undergone active maneuvering between adjacent splicing points are extracted to form a free decay sample library; The first-order forward difference of the semi-major axis between adjacent epochs is obtained from the data in the free decay sample library and standardized to meters. Based on Gaussian 3- σ The principle or interquartile range criterion is used to define the anomaly boundary, and the maneuvering events are classified into track lift, system reset maneuver, high-frequency micro-track maintenance maneuver, large-scale descent avoidance, and mission deorbiting maneuver. The occurrence time, net change of semi-major axis, and orbit change rate of each type of maneuver are extracted, and combined with the natural decay baseline and noise level, to form a multi-dimensional maneuvering feature vector.

[0012] Preferably, the mean, standard deviation, skewness, and kurtosis of the difference sequence are calculated, and the distribution characteristics are verified by a QQ plot; wherein, the mean of the difference sequence corresponds to the average natural resistance attenuation baseline, and the standard deviation corresponds to the overall background noise level.

[0013] A computer-readable storage medium having a computer program stored thereon, which is executed by a processor to implement a method for identifying and detecting low-Earth orbit spacecraft maneuvers.

[0014] A computer device includes a memory and a processor, the memory storing a computer program, and the processor executing the computer program to implement a method for identifying and detecting low-Earth orbit spacecraft maneuvers.

[0015] The present invention has the following beneficial effects: This invention analyzes the dynamic responses of different types of maneuvering strategies and parameters, establishing the quasi-average semi-major axis, reflecting system energy changes, as a core feature parameter. Subsequently, addressing the challenge of environmental noise and maneuvering signal aliasing in long-term ephemeris data, a dual-sliding denoising detection algorithm, incorporating adaptive thresholding and considering data splicing errors, is progressively constructed, forming a complementary screening and correction mechanism. Finally, based on accurate maneuvering event localization, a comprehensive statistical analysis of the spatial perturbation background and maneuvering increments is conducted to extract highly discriminative physical features, laying a reliable data foundation for building intelligent maneuvering behavior prediction models.

[0016] The pseudo-average semi-major axis improves the signal-to-noise ratio by an order of magnitude compared to the instantaneous close semi-major axis, and can keenly reflect minute changes in orbital energy; the 15×15 order EGM2008 model achieves the optimal balance between accuracy and computing power, and is suitable for the batch processing needs of giant constellations.

[0017] The adaptive detection is robust: by dynamically learning the natural attenuation baseline of the current space environment, it replaces the traditional fixed threshold and can adapt to atmospheric density fluctuations, making it suitable for low-Earth orbit satellites at different orbital altitudes and times.

[0018] The weak signal detection capability is significantly improved: the dual sliding denoising and energy transition detection logic can capture the continuous deviation of orbital energy, rather than relying solely on instantaneous jumps, and successfully identify the minute maneuvers that are missed by the adaptive threshold method; the splicing error investigation eliminates false signals from the data source, reducing the false detection rate by more than 85%.

[0019] The complementary detection system is reliable: adaptive threshold detection is responsible for quickly screening significant maneuvers, while dual sliding denoising detection is responsible for filling in the gaps. The two algorithms form a mutual verification mechanism, increasing the maneuver detection coverage to 100%.

[0020] High purity of feature extraction: The natural decay baseline and noise boundary quantified through systematic statistical analysis provide theoretical support for maneuver determination; the extracted multidimensional feature vectors have clear physical meaning and can be directly imported into machine learning models, laying a solid data foundation for subsequent maneuver intent prediction. Attached Figure Description

[0021] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0022] Figure 1 This is shown as the adaptive threshold learning detection process of the present invention; Figure 2 The flowchart shown is a process flow chart of splicing error processing and dual sliding noise reduction detection of the present invention. Figure 3 This is shown as a close semi-major axis variation for satellite 64378; Figure 4 The display shows the average semi-major axis variation considering the J2 perturbation; Figure 5 The display shows the average semi-major axis variation considering 15*15 order gravitational perturbations; Figure 6 The display shows the average semi-major axis variation considering 20*20 order gravitational perturbations; Figure 7 The results are displayed as the satellite detection results with a magnification factor of 3979. Figure 8 The results are displayed as the detection results of satellite 5171 with a magnification factor of 3x. Figure 9 The results are displayed as the detection results of satellite 5564 with a magnification factor of 3x. Figure 10 The result is displayed as a 2.5x multiplication of satellite 3979 detection results; Figure 11 The results are displayed as a 2.5x magnification of the satellite 5171 detection results. Figure 12 The result is displayed as a 2.5x detection result from satellite 5564. Figure 13 This shows the supplementary detection results for satellite 3979; Figure 14 This shows the supplementary detection results for satellite 5171. Figure 15 This shows the supplementary detection results for satellite 5564; Figure 16 The display shows the distribution characteristics of the pseudo-mean semi-major axis difference between adjacent epochs of the satellite at 3979 epochs. Figure 17 The chart shows the QQ test results for the difference in the pure attenuation segment of satellite 3979. Figure 18 This is displayed as the clean orbital decay sequence after satellite 3979 was stripped and spliced; Figure 19 The display shows the distribution characteristics of the pseudo-mean semi-major axis difference between adjacent epochs of satellite 5171; Figure 20 The chart shows the QQ test plot of the difference in the pure attenuation segment of satellite 5171. Figure 21 This is displayed as the clean orbital decay sequence after satellite 5171 was stripped and spliced; Figure 22The display shows the distribution characteristics of the approximate average semi-major axis difference between satellites at 5564 adjacent epochs. Figure 23 The chart shows the QQ test plot of the difference in the pure attenuation segment of satellite 5564. Figure 24 The image shows the clean orbital decay sequence after satellite 5564 was stripped and spliced. Detailed Implementation

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

[0024] The present invention will be described in detail below with reference to specific embodiments. Specific Implementation Example 1: according to Figures 1 to 24 As shown, the specific optimized technical solution adopted by the present invention to solve the above-mentioned technical problems is: The present invention relates to a method for identifying and detecting low-orbit spacecraft maneuvers.

[0026] This invention provides a method for identifying and detecting low-Earth orbit spacecraft maneuvers, the method comprising the following steps: Step 1: Using the quasi-average semi-major axis as the core characteristic parameter characterizing the change in the spacecraft's orbital energy, calculate the non-spherical gravitational perturbation solution based on the Earth's gravitational field model, truncate and select the field harmonic and zone harmonic perturbations, and reconstruct the quasi-average semi-major axis time series; Step 2: Perform a two-step dynamic adaptive threshold detection based on the quasi-average semi-major axis difference sequence to complete the preliminary localization and classification of the maneuvering event; Step 3: Perform splicing error investigation and preprocessing on the original pseudo-average semi-major axis sequence, then reconstruct the orbital energy trend baseline by denoising using the double moving average method, and perform refined maneuver detection based on the energy transition signal to make up for the omissions in the preliminary detection. Step 4: Perform statistical analysis based on clean free decay trajectory data, quantify the natural decay baseline and environmental noise boundary, and extract multidimensional high-confidence maneuver feature vectors to provide input for subsequent maneuver prediction models.

[0027] The two-step dynamic adaptive threshold detection method is as follows: The first-order forward differencing process is performed on the quasi-average semi-major axis time series to obtain the semi-major axis change sequence between adjacent epochs; a conservative fixed initial empirical threshold is introduced, and the sequence is traversed to perform the first round of coarse screening. Epochs with absolute difference values ​​greater than the initial threshold are marked as maneuvering, and those with absolute difference values ​​less than the initial threshold are marked as free decay; the state switching index nodes are recorded to divide the continuous sequence into multiple independent data segments with single attributes. Extract all free decay segments with a length greater than 1 data point, and calculate the average of the absolute values ​​of the differences between adjacent points within each segment; use the number of data points in each free decay segment as the weight to perform a global weighted average of all local average differences to obtain the weighted average decay deviation of the current space environment; multiply the weighted average decay deviation by the dynamic amplification factor as the adaptive maneuver decision threshold; if no valid free decay samples are extracted, the initial empirical threshold is used. Using an adaptive threshold as a benchmark, a second round of refined screening was conducted on the original differential sequence to redefine the free decay segment and the active maneuver segment; linear regression fitting was performed on the free decay segment to calculate the standardized decay rate; for the active maneuver segment, the instantaneous orbit change rate of the pulse maneuver and the average orbit change rate of the continuous thrust maneuver were calculated respectively.

[0028] The Earth's gravitational field model uses the high-precision Earth gravitational field model EGM2008.

[0029] The dynamic magnification factor ranges from 2.0 to 3.0.

[0030] The dual sliding noise reduction and splicing error supplementation detection are as follows: Extract the overlapping time nodes of the ephemeris update handover, calculate the quasi-average semi-major axis difference between the old and new batches of ephemeris at the same epoch, lock the index position corresponding to the splicing error and mark it separately, and eliminate false maneuver signals. Using the spacecraft orbital period as the basic sliding window, the original quasi-average semi-major axis sequence is initially smoothed to filter out residual perturbations at the orbital period level; then, using half an orbital period as the secondary sliding window, the initial results are filtered a second time to extract the baseline of the long-term evolution trend of orbital energy. Set a step window of half an orbital period length, slide it along the time series, and calculate the difference between the average value of the future window data block and the historical window data block at each central evaluation point as the transition energy intensity at that point; Boundary data with one orbital period length at the beginning and end of the sequence are excluded. The 90th percentile of the absolute value of the transition signal in the stable segment is extracted as the environmental noise baseline. Differentiated peak-finding thresholds are set: the orbit-ascending detection threshold is 0.8 times the noise baseline, and the orbit-descending detection threshold is 1.5 times the noise baseline. At the same time, the peak prominence is required to be greater than 0.25 times the baseline, and adjacent maneuvering events must meet the preset time interval constraint to initially identify suspected maneuvering nodes.

[0031] The dual-sliding denoising and splicing error supplementary detection also includes a secondary verification process: For each suspected maneuver peak point, a stable reference interval of one step window length is extracted before and after it. The difference between the mean of the semi-major axis of the interval before and after the maneuver is taken as the net maneuver effect. When the absolute value of the net effect exceeds the calculated maximum local average deviation, it is confirmed as a valid trajectory change point.

[0032] The comprehensive statistical extraction of maneuver characteristics is specifically as follows: By segmenting the observation time axis using the splicing node database, continuous data segments that have not undergone active maneuvering between adjacent splicing points are extracted to form a free decay sample library; The first-order forward difference of the semi-major axis between adjacent epochs is obtained from the data in the free decay sample library and standardized to meters. Based on Gaussian 3- σ The principle or interquartile range criterion is used to define the anomaly boundary, and the maneuvering events are classified into track lift, system reset maneuver, high-frequency micro-track maintenance maneuver, large-scale descent avoidance, and mission deorbiting maneuver. The occurrence time, net change of semi-major axis, and orbit change rate of each type of maneuver are extracted, and combined with the natural decay baseline and noise level, to form a multi-dimensional maneuvering feature vector.

[0033] The mean, standard deviation, skewness, and kurtosis of the difference sequence are calculated, and the distribution characteristics are tested using a QQ plot. The mean of the difference sequence corresponds to the average natural resistance attenuation baseline, and the standard deviation corresponds to the overall background noise level.

[0034] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which is executed by a processor to implement a method for identifying and detecting low-Earth orbit spacecraft maneuvers.

[0035] The present invention also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement a method for identifying and detecting low-Earth orbit spacecraft maneuvers. Specific Implementation Example 2: The only difference between Embodiment 2 and Embodiment 1 of this application is that: Now, based on the quasi-average root method, using a version of ephemeris data released by the Starlink satellite STARLINK-11711 (NORAD satellite number 64378) as an example, we conduct a semi-major axis analysis. When calculating the non-spherical gravitational perturbation solution, we employ the widely cited high-precision Earth gravitational field model EGM2008. EGM2008 integrates GRACE satellite gravity data, high-resolution terrain data, and ground gravity anomaly data, achieving a leap of more than six times over the previous generation EGM96 model in terms of spatial resolution and accuracy.

[0037] The instantaneous semi-major axis variation of the satellite's orbit is as follows: Figure 3 As shown, consider The change in the mean semi-major axis of the orbit after the perturbation term is as follows: Figure 4 As shown, Figure 5 and Figure 6 These are schematic diagrams showing the changes in the mean semi-major axis of the orbit considering the 15*15 order and 20*20 order Tian Harmonic gravitational perturbation terms, respectively.

[0038] By analyzing and comparing the satellite semi-major axis variation diagrams considering different perturbation terms and orders, it can be found that ephemeris data, as closely related data, has an average perturbation magnitude of around 10 km. Therefore, it is difficult to extract information about the satellite's own on-orbit operations; its slight maneuvering information is hidden in the various perturbation terms of the low-Earth orbit gravitational field. Comparing the variations in the closely related semi-major axis, considering... The average semi-major axis of each order of gravitational perturbation term highlights its advantage in maneuver detection. Furthermore, considering higher-order gravitational perturbation terms more clearly demonstrates the changes in the satellite's semi-major axis. However, from... Figure 3 The ephemeris data revealed that the gravitational model used for the last day's orbit was inconsistent with that used for the previous two days, indicating that only one side of the orbit was considered. The theoretical orbital model for the perturbation terms. This is likely due to the frequent maneuvers of Starlink satellites in low Earth orbit (LEO), providing only the high-precision model for the first two days and using the theoretical orbital model for the last day is more suitable for practical applications. In the practice of batch analysis of mega-constellations, considering the high-order theoretical model bias inherent in the ephemeris data itself, this framework truncates and selects the first 15*15 order EGM2008 perturbation model when calculating the field harmonic and band harmonic perturbations. This order configuration can fully cover the long-wave and medium-wave gravity anomaly signals in the LEO environment, effectively reconstruct the quasi-averaged semi-major axis sequence, and significantly reduce the computational power consumption of batch constellation solutions.

[0039] Starlink satellites can be divided into different phases based on their on-orbit status: Raising, Parking, Operational, and Deorbit.

[0040] Taking four Starlink satellites in orbit—STARLINK-11275, STARLINK-35428, STARLINK-3374, and STARLINK-1010—as examples, we analyze their maneuvering strategies, which correspond to the four different on-orbit operation phases mentioned above. Starting from the initial input state after satellite launch, the satellite first enters a parking orbit of approximately 360 kilometers. In this phase, the satellite performs initial deployment by adjusting its orbital parameters and prepares to enter its target operational orbit. Subsequently, the satellite gradually raises its orbit, eventually entering the operational phase of approximately 540 kilometers, where the mission officially begins. During this process, the satellite adjusts to different orbital planes to optimize global coverage performance and improve communication or observation capabilities. At the end of the mission, the satellite performs an orbital descent operation, ultimately re-entering the atmosphere and burning up. The maneuvering characteristics differ significantly between these phases. During the parking phase, the satellite needs to perform frequent maneuvers to maintain its orbital altitude, while the frequency of maneuvers decreases once the satellite enters the operational phase. The study is mainly based on the on-orbit working phase and the parking phase. The maneuvering characteristics of the raising and lowering phases are relatively regular, so the raising and lowering phases will not be discussed further below.

[0041] Finally, regarding the orbit descent phase of the Starlink satellite STARLINK-3047, a sudden orbit descent was observed, with a decay rate far exceeding that of a normal maneuvering orbit descent process. This indicates that the Starlink satellite actively implemented an orbit descent operation at the end of its life cycle to avoid posing a collision threat to other satellites in the same orbit.

[0042] For Starlink satellites, their on-orbit operational phase is crucial for performing Earth coverage observations and space reconnaissance missions. Therefore, studying their maneuvering during this phase is the primary objective of this research. Maneuver detection using satellite ephemeris data over a continuous operational period is more effective in uncovering their active maneuvering information. For selecting the dynamic threshold for Starlink satellites, the initial ephemeris data is based on close-root numbers. After considering perturbation terms of various orders, their on-orbit maneuvering can be visually represented. We then randomly selected three Starlink satellites—STARLINK-3979, STARLINK-5171, and STARLINK-5564—for continuous long-term maneuver detection. It was found that during the on-orbit operational phase, Starlink satellites primarily perform orbit maintenance maneuvers, but also engage in maneuvers such as orbit descent and evasion maneuvers. Therefore, for such frequent maneuvers, the key to obtaining maneuvering information within a given timeframe lies in determining the maneuvering threshold. However, Starlink satellite ephemeris data differs from TLE data; its maneuvering frequency is high but the magnitude is small, making short-period perturbation terms unsuitable as threshold boundaries. Since the magnitude of the maneuver maintained by Starlink satellites is very similar to the magnitude of the gravitational perturbation between the data points of the quasi-average semi-major axis, the two are difficult to distinguish clearly. Therefore, using a fixed-magnitude maneuver threshold is not convenient for general application. Moreover, there are currently nearly 10,000 Starlink satellites in orbit, making it difficult to determine a fixed maneuver threshold for each satellite at each orbital altitude during the operational phase.

[0043] Adaptive threshold detection analysis For the frequent and minute orbit maintenance and avoidance maneuvers of low-Earth orbit mega-constellations like Starlink during their on-orbit operation, conventional single fixed thresholds are insufficient to effectively distinguish between natural perturbation decay and active thrust effects in long-term, noisy ephemeris data. To improve the robustness and environmental adaptability of maneuver detection, this study proposes and implements a two-step dynamic adaptive threshold detection algorithm based on quasi-averaged semi-major axis difference sequences. The core logic of this algorithm consists of three main stages: preliminary segmentation, dynamic threshold learning, and fine-grained reclassification. The specific analysis process is as follows: Figure 1 As shown in the flowchart above, the program uses the time-series quasi-average semi-major axis data as input variables, performs first-order forward differencing on it, and obtains the sequence of changes in the semi-major axis between adjacent epochs. Given that the current atmospheric drag attenuation characteristics of the target satellite are not yet fully understood, a conservative fixed initial empirical threshold is introduced. This is then iterated over the entire... The sequence undergoes a first round of coarse screening based on the initial threshold. The judgment logic is as follows: when the absolute value of the difference between the semi-major axes of a certain epoch is greater than the initial empirical threshold, its state is marked as maneuvering; otherwise, it is marked as free decay. During this traversal, the index nodes where the state changes are recorded, thereby dividing the originally continuous time series into multiple independent data segments with single attributes. The main purpose of this step is to select as many pure decay samples affected by natural resistance as possible, providing data support for subsequent adaptive learning.

[0044] Next, based on the initial segmentation results, adaptive threshold calculation is performed. All free attenuation segments generated in the first round of segmentation are specifically extracted. For effective attenuation segments longer than one data point, the average level of the absolute difference between adjacent points within each segment is calculated to characterize the average natural attenuation intensity within that local time period. Subsequently, to comprehensively consider the background noise baseline across the entire long observation arc, the number of data points contained in each effective attenuation segment is used as a weighting coefficient. A global weighted average sum of all local average differences is then performed to obtain the weighted average attenuation deviation under the current space environment. Based on this global statistical characteristic, the algorithm establishes a dynamic threshold update mechanism. The calculated weighted average deviation is multiplied by a dynamic amplification factor and set as the new adaptive maneuver judgment threshold. Through this self-learning mechanism, the system can effectively accommodate normal parameter fluctuations caused by atmospheric density fluctuations due to space activities. If, in extreme cases, no effective free attenuation samples are extracted, the algorithm has a fallback mechanism, continuing to use the initial empirical threshold as the judgment basis to ensure continuous operation of the program.

[0045] After obtaining a new dynamic threshold that fits the characteristics of the current environment, the program processes the original difference sequence. A second round of refined screening is conducted. Following the same state machine switching logic as the initial segmentation, the free decay segment and the active maneuver segment are redefined using a new threshold as a benchmark. After final discrimination, the algorithm further analyzes the dynamic parameters for different categories of data segments: For the free decay segment, for continuous decay intervals containing three or more data points, a first-order polynomial is used for linear regression fitting, and combined with the sampling time interval of adjacent ephemeris points, the slope of the semi-major axis is converted into a standardized decay rate in meters per second to assess the atmospheric damping effect on the satellite. For the active maneuver segment analysis, the focus is on evaluating the net effect of the maneuver. The absolute change of the semi-major axis is obtained by extracting the pseudo-average semi-major axis values ​​of the last node and the historical moment before the starting node of the maneuver segment and calculating the difference. In addition, for continuous thrust maneuvers spanning multiple epochs, the program also uses linear fitting to calculate the average orbit change rate of the maneuver segment; while for pulse maneuvers consisting of only a single point, the instantaneous orbit change rate is directly solved based on the epoch time difference.

[0046] Through the above two-step judgment and adaptive learning process, the interference of perturbation model residuals in complex low-Earth orbit environments can be effectively overcome, enabling precise positioning and quantitative analysis of minute maneuvering events of spacecraft over continuous time periods. Simultaneously, maneuver detection was performed on the three Starlink satellites in orbit, and the detection results at different magnification factors are as follows. Figures 7 to 12 As shown.

[0047] Analysis of splicing error detection In the two-step adaptive threshold maneuver detection, although a dynamic threshold mechanism is introduced, it was found in the batch processing practice of giant constellations that the amplification factor, which plays a decisive role in the algorithm, is difficult to achieve global applicability. Satellites in different time and space environments have significantly different natural decay characteristics, and using a fixed or single-related amplification factor can easily lead to a large number of false detections or missed detections for some satellites. More importantly, in the deep backtracking analysis of suspected maneuver points, it was found that during the multiple batches of release and update of public ephemeris data, there are often non-physically real data jumps at the handover time points between adjacent batches. This data discontinuity caused by the accumulation of orbit extrapolation model errors and the refitting and updating of orbital elements is regarded as data splicing error. This error is very similar in numerical performance to the semi-major axis step caused by small maneuvers, which seriously interferes with the accuracy of the detection algorithm.

[0048] To address the aforementioned issues, the original motion detection algorithm was supplemented with a new detection method: an energy transition detection strategy that integrates splicing error detection and dual sliding window denoising. The specific logical flow of this algorithm is as follows: Figure 2 As shown.

[0049] To effectively isolate false maneuvering signals caused by data splicing, independent splicing point time series can be introduced for cross-comparison during the preprocessing stage. By extracting the overlapping time nodes during ephemeris update handover, the difference between the pseudo-average semi-major axis values ​​of the old and new batches of ephemeris at the same epoch is calculated within a very small time tolerance range to obtain the true value of the splicing error. Subsequently, these determined splicing time points are mapped back to the time vector of the main analysis data, locking their corresponding index positions. In the subsequent maneuvering determination and result output stages, these abnormal jump nodes caused by data splicing will be marked separately and will not be included in the orbit change operations actively performed by the satellite, thus avoiding false maneuvering events at the data source. Considering that some periodic perturbation noise still remains in the pseudo-average semi-major axis sequence of low-Earth orbit satellites, direct difference operations are prone to amplifying local high-frequency disturbances, which is why the amplification factor is difficult to determine. Therefore, the program uses a double moving average method to denoise and reconstruct the original sequence. First, the orbital period of the satellite is used as the basic moving window span to perform initial smoothing on the original data to filter out residual perturbations at the orbital period level. A secondary smoothing window of half an orbital period is then introduced to perform a second filtering on the initial results. After this step, the complex nonlinear sequence is effectively smoothed, and a clean trend baseline that can characterize the long-term evolution of orbital energy is extracted.

[0050] After obtaining the smoothing trend, instead of using simple neighbor-to-neighbor difference, an energy transition signal based on data blocks is constructed. An initial step window of approximately half an orbital period is set and slid along the time series. At each central evaluation point, the average value of the future window data block and the average value of the historical window data block are calculated, and the difference between the two is taken as the transition energy intensity at that central point. The physical meaning of this processing logic is that if the satellite is in a stable natural decay period, the difference in the mean values ​​of the preceding and following data blocks is small. However, when the satellite performs a thrust maneuver, the orbital energy undergoes a step jump, and this difference signal will form a significant local peak at the moment of the maneuver.

[0051] To extract the true maneuvering nodes from the transition signal sequence, the program excluded data from both ends of an orbital period length affected by boundary effects, selecting the signal from the middle stable segment for statistical analysis. The 90th percentile of the absolute value of the stable segment signal was extracted as the environmental noise floor for the current sequence. Given that the thrust patterns and amplitudes of low-Earth orbit satellites' orbital maintenance maneuvers against atmospheric drag and mission evasion maneuvers' orbital descent maneuvers often exhibit asymmetric characteristics, differentiated peak-finding thresholds were established. For orbital maintenance detection, the peak value of the transition signal was required to be no less than 0.8 times the noise floor; for orbital descent detection, the absolute value of the negative peak value was required to be no less than 1.5 times the noise floor. Simultaneously, the peak prominence needed to be greater than 0.25 times the floor, and adjacent maneuvering events needed to meet reasonable time interval constraints. These constraints allowed for the initial identification of suspected active maneuvering nodes.

[0052] After extracting extreme points, a rigorous secondary verification procedure is performed to ensure the reliability of the results. For each suspected maneuver peak, the algorithm extracts stable reference intervals of the step window length before and after its occurrence. The true net physical effect of the maneuver is obtained by calculating the difference between the mean of the semi-major axis of the stable interval after the maneuver and the mean of the semi-major axis of the stable interval before the maneuver. Only when the absolute value of this net orbit change effect exceeds the calculated maximum local average deviation is the suspected point finally confirmed as a valid orbit change point. Finally, a detailed analysis report is output, including the maneuver time index, energy transition intensity, net orbit change amount of the semi-major axis, and splice point investigation results.

[0053] Further testing was conducted on the three Starlink satellites in orbit, and the results are as follows. Figures 13-15 As shown: Parallel detection was performed on the orbital data of the first target satellite, STARLINK-3979, during the observation period. The detection results of adaptive threshold detection and dual-sliding denoising detection were summarized and compared. As shown in Table 1, adaptive threshold detection, with its direct sensitivity to semi-major axis step signals, accurately identified most significant orbit change nodes. However, in the in-depth analysis of this satellite, it was found that the dual-sliding denoising detection method successfully identified the "Ascendance VII" maneuver recorded in Table 2, while this event was missed in the screening by the adaptive threshold method because it failed to meet the judgment criteria.

[0054] Table 1. Detection results of adaptive threshold multiplier 3.

[0055] Table 2 Results of Dual-Sliding Denoising and Stitching Troubleshooting

[0056] This difference in detection reflects the different advantages and complementary characteristics of the two algorithms at the underlying logic level. The adaptive threshold detection method mainly relies on energy mutations between adjacent epochs. While it has high timeliness in handling large pulse maneuvers, when the satellite is in a space environment with severe fluctuations or performing extremely small thrust operations, its statistical boundary is widened in real time by background noise, making small maneuver signals easily masked under the dynamic threshold. In contrast, the dual sliding denoising detection method effectively filters out high-frequency jitter by converting energy transition signals and using the mean comparison of long and short windows. Its core logic lies in capturing the continuous shift of orbital energy rather than instantaneous jumps. Therefore, when dealing with maneuver events like "Sheng-7" with small energy increments and relatively hidden characteristics, the sliding denoising method demonstrates stronger weak signal extraction capabilities, compensating for the limitations of the adaptive single threshold judgment.

[0057] Comprehensive analysis of maneuver In summary, adaptive threshold detection provides a benchmark for the rapid classification of long-term time-series data, while dual-sliding denoising detection plays a crucial role in identifying and filling gaps by accurately capturing energy trends. Although the two methods differ significantly in their algorithmic construction, they can form an effective mutual verification mechanism in practical applications. Through this complementary detection strategy, the system not only improves the coverage of various complex maneuvering modes but also provides multi-dimensional evidence for the subsequent extraction of high-purity maneuvering feature data. This multi-strategy fusion analysis mode effectively ensures the robustness and reliability of spacecraft comprehensive maneuvering detection, laying a reliable data foundation for the subsequent construction of maneuvering behavior prediction models based on high-confidence feature dimensions.

[0058] Through maneuver detection analysis, the results have verified the sensitivity of the quasi-average semi-major axis in maneuver feedback, and a maneuver detection framework integrating adaptive thresholds and splicing error investigation has been successfully constructed. Simultaneously, the temporal location and magnitude extraction of active orbit change events have been achieved. However, the core purpose of maneuver detection is not only post-event identification and verification, but also to achieve forward-looking prediction of the target's future maneuver intentions through deep learning from massive historical data. To provide high-purity, high-confidence feature vector inputs for the next step of the machine learning prediction model, statistical analysis tools are further introduced based on the previous detection results to conduct a macroscopic feature distribution study on the pure orbit evolution sequence after removing maneuvers and false anomalies.

[0059] Because publicly available ephemeris data contains data jumps caused by orbit determination updates, directly performing global statistics on the raw data would severely contaminate its distribution characteristics. Therefore, an established splicing node database is used to precisely segment the entire observation timeline. The system extracts continuous data segments located between adjacent splicing points that have not undergone active maneuvering, forming a true free decay sample library. After ensuring the purity of the samples, the first-order forward difference of the semi-major axis between adjacent epochs is calculated for the extracted continuous segments. The units are standardized to meters. This differential sequence avoids the orbital altitude cardinality and directly reflects the microscopic intake and dissipation rate of orbital energy within a unit sampling period.

[0060] To investigate the statistical laws governing the effects of atmospheric drag and other high-frequency residual perturbations on low-Earth orbit spacecraft under natural conditions, the above-mentioned... The sequences were subjected to multi-dimensional statistical measurements and visualization probability density fitting. Supported by a large number of normal flight cycle samples, the difference distribution typically exhibits an approximate bell-shaped curve. Several core statistical indicators were also measured. The sequence mean is usually a small negative value, which physically corresponds precisely to the average natural drag attenuation baseline of the target satellite under its current orbital altitude and space environment. The sequence standard deviation characterizes the combined background noise level caused by fluctuations in space atmospheric density and measurement fitting residuals. Furthermore, skewness was introduced to measure distribution symmetry, and kurtosis to measure the steepness and long-tail characteristics of the distribution. By combining nonparametric tests using quantile-quantile plots (QQ), the approximate Gaussian normal distribution of the environmental noise base was further evaluated. This rigorous statistical testing helps to identify whether there are still hidden, minor maneuvers or model systematic errors in the data that have not been properly cleaned up. The specific difference statistics for the three satellites are as follows: Figures 16-24 As shown: After identifying the natural decay baseline and its noise distribution characteristics, it is necessary to provide theoretical support from a mathematical and statistical perspective for defining the degree of trajectory change as a maneuver. The system compared two boundary delimitation criteria in its analysis: one based on the Gaussian assumption and the other on a 3-... σ In principle, this boundary can cover most natural parameter fluctuations; another is the more robust interquartile range (IQR) criterion, which defines outlier boundaries by calculating the difference between the third quartile and the first quartile. This method shows better robustness when dealing with spatial data that is not ideally normally distributed and contains complex long-tailed noise.

[0061] Through the aforementioned macroscopic statistical and comprehensive analysis, the complex and disordered ephemeris data was highly refined. The satellite's natural decay rate baseline, background noise standard deviation level, and the occurrence time and net change of orbital lift or system reset maneuvers (ascent bias 1), high-frequency small orbital maintenance maneuvers (ascent bias 2), and large descent avoidance and mission deorbiting maneuvers (descent bias) extracted after robust criterion cleaning, together constitute a set of time-series feature vectors with clear physical meaning and concise dimensions. Adaptive threshold detection provides a benchmark reference for the rapid classification of long-term time-series data, while dual-sliding denoising detection plays a crucial role in identifying and filling gaps by accurately capturing energy trends. Although the two methods are quite different in algorithm construction, they can form an effective mutual verification mechanism in practical applications. Through this complementary detection strategy, the system not only improves the coverage of various complex maneuver modes, but also provides multi-dimensional basis for the subsequent extraction of high-purity maneuver feature numbers. This multi-strategy fusion analysis model effectively ensures the robustness and reliability of the comprehensive maneuver detection of spacecraft. At the same time, these maneuver characteristic parameters, which have been rigorously tested and quantified by statistics, will be directly used as input variables into subsequent machine learning models, laying a reliable data foundation for the subsequent development of orbital maneuver simulation and time node prediction models for target spacecraft based on high-confidence feature dimensions.

[0062] This invention addresses the identification and parameter analysis of low-Earth orbit (LEO) spacecraft maneuvering behavior, establishing a complete data processing chain from feature selection to anomaly detection and statistical refinement. Based on the higher signal-to-noise ratio and dynamic mapping feedback of the quasi-averaged semi-major axis in characterizing system energy transitions, a comprehensive complementary detection algorithm is proposed and validated, integrating adaptive threshold learning, dual sliding denoising, and stitching error detection, to address the complex perturbation background and ephemeris update jumps in real observation data. This effectively eliminates false maneuvering signals, significantly improving the confidence level of maneuvering event identification. Through in-depth statistical analysis of the pure adjacent observation difference sequence after anomaly removal, the natural decay baseline and environmental noise boundary can be quantified, establishing robust maneuvering feature judgment criteria. This not only provides an efficient numerical calculation reference for daily on-orbit behavior monitoring of spacecraft, but also provides a necessary and complete data foundation framework for conducting machine learning-based spacecraft maneuvering prediction research.

[0063] The above description is merely a preferred embodiment of a method for identifying and detecting low-Earth orbit spacecraft maneuvers. The scope of protection for this method is not limited to the above embodiments; all technical solutions falling within this conceptual framework are within the scope of protection of this invention. It should be noted that for those skilled in the art, any improvements and variations made without departing from the principles of this invention should also be considered within the scope of protection of this invention.

Claims

1. A method for identifying and detecting low-Earth orbit spacecraft maneuvers, characterized by: The method includes the following steps: Step 1: Using the quasi-average semi-major axis as the core characteristic parameter to characterize the change in the spacecraft's orbital energy, calculate the non-spherical gravitational perturbation solution based on the Earth's gravitational field model, truncate and select the field harmonic and zone harmonic perturbations, and reconstruct the quasi-average semi-major axis time series; Step 2: Perform two-step dynamic adaptive threshold detection based on the quasi-average semi-major axis difference sequence to complete the preliminary localization and classification of the maneuvering event; Step 3: Perform splicing error preprocessing on the original pseudo-average semi-major axis sequence, then reconstruct the orbital energy trend baseline by denoising using the double moving average method, and perform refined maneuver detection based on the energy transition signal to make up for the omissions in the preliminary detection. Step 4: Perform statistical analysis based on clean free decay trajectory data, quantify the natural decay baseline and environmental noise boundary, and extract multidimensional high-confidence maneuver feature vectors to provide input for subsequent maneuver prediction models.

2. The method according to claim 1, characterized in that: The two-step dynamic adaptive threshold detection method is as follows: The first-order forward differencing process is performed on the quasi-average semi-major axis time series to obtain the semi-major axis change sequence between adjacent epochs; a conservative fixed initial empirical threshold is introduced, and the sequence is traversed to perform the first round of coarse screening. Epochs with absolute difference values ​​greater than the initial threshold are marked as maneuvering, and those with absolute difference values ​​less than the initial threshold are marked as free decay; the state switching index nodes are recorded to divide the continuous sequence into multiple independent data segments with single attributes. Extract all free decay segments with a length greater than 1 data point, and calculate the average of the absolute values ​​of the differences between adjacent points within each segment; use the number of data points in each free decay segment as the weight to perform a global weighted average of all local average differences to obtain the weighted average decay deviation of the current space environment; multiply the weighted average decay deviation by the dynamic amplification factor as the adaptive maneuver decision threshold; if no valid free decay samples are extracted, the initial empirical threshold is used. Using an adaptive threshold as a benchmark, a second round of refined screening was conducted on the original differential sequence to redefine the free decay segment and the active maneuver segment; linear regression fitting was performed on the free decay segment to calculate the standardized decay rate; for the active maneuver segment, the instantaneous orbit change rate of the pulse maneuver and the average orbit change rate of the continuous thrust maneuver were calculated respectively.

3. The method according to claim 2, characterized in that: The Earth's gravitational field model uses the high-precision Earth gravitational field model EGM2008.

4. The method according to claim 3, characterized in that: The dynamic magnification factor ranges from 2.0 to 3.

0.

5. The method according to claim 4, characterized in that: dual The specific steps for sliding noise reduction and splicing error supplementation detection are as follows: Extract the overlapping time nodes of the ephemeris update handover, calculate the quasi-average semi-major axis difference between the old and new batches of ephemeris at the same epoch, lock the index position corresponding to the splicing error and mark it separately, and eliminate false maneuver signals. Using the spacecraft orbital period as the basic sliding window, the original quasi-average semi-major axis sequence is initially smoothed to filter out residual perturbations at the orbital period level; then, using half an orbital period as the secondary sliding window, the initial results are filtered a second time to extract the baseline of the long-term evolution trend of orbital energy. Set a step window of half an orbital period length, slide it along the time series, and calculate the difference between the average value of the future window data block and the historical window data block at each central evaluation point as the transition energy intensity at that point; Boundary data with one orbital period length at the beginning and end of the sequence are excluded. The 90th percentile of the absolute value of the transition signal in the stable segment is extracted as the environmental noise baseline. Differentiated peak-finding thresholds are set: the orbit-ascending detection threshold is 0.8 times the noise baseline, and the orbit-descending detection threshold is 1.5 times the noise baseline. At the same time, the peak prominence is required to be greater than 0.25 times the baseline, and adjacent maneuvering events must meet the preset time interval constraint to initially identify suspected maneuvering nodes.

6. The method according to claim 5, characterized in that: dual Sliding denoising and splicing error supplementary detection also include a secondary verification process: For each suspected maneuver peak point, a stable reference interval of one step window length is extracted before and after it. The difference between the mean of the semi-major axis of the interval before and after the maneuver is taken as the net maneuver effect. When the absolute value of the net effect exceeds the calculated maximum local average deviation, it is confirmed as a valid trajectory change point.

7. The method according to claim 6, Its characteristics are: the comprehensive statistical extraction of maneuver characteristics specifically includes: By segmenting the observation time axis using the splicing node database, continuous data segments that have not undergone active maneuvering between adjacent splicing points are extracted to form a free decay sample library; The first-order forward difference of the semi-major axis between adjacent epochs is obtained from the data in the free decay sample library and standardized to meters. Based on Gaussian 3- σ The principle or interquartile range criterion is used to define the anomaly boundary, and the maneuvering events are classified into track lift, system reset maneuver, high-frequency micro-track maintenance maneuver, large-scale descent avoidance, and mission deorbiting maneuver. The occurrence time, net change of semi-major axis, and orbit change rate of each type of maneuver are extracted, and combined with the natural decay baseline and noise level, to form a multi-dimensional maneuvering feature vector.

8. The method according to claim 7, characterized in that: The mean, standard deviation, skewness, and kurtosis of the difference sequence are calculated, and the distribution characteristics are tested using a QQ plot. The mean of the difference sequence corresponds to the average natural resistance attenuation baseline, and the standard deviation corresponds to the overall background noise level.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, The program is executed by the processor to implement the method as claimed in any one of claims 1-8.

10. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the method of any one of claims 1-8.

Images