A method and system for detecting orbital anomalies based on anomaly probability

Abstract

This invention discloses a method and system for detecting orbital anomalies based on anomaly probability. The method includes: acquiring satellite orbit data; reconstructing the satellite orbit data using a wavelet analysis method combined with forecasts to obtain the current orbital inclination and semi-major axis information after wavelet transform; calculating sub-criteria and comprehensive criteria and combining them with a limiting threshold to perform anomaly detection, obtaining the detection result; and introducing anomaly probability to map the comprehensive criteria to obtain the probability interval of the detection result. The system includes: a data acquisition module, a wavelet analysis module combined with forecasts, a detection module, and a probability calculation module. By using this invention, it is possible to quickly and accurately determine whether an orbital anomaly has occurred. This invention, as a method and system for detecting orbital anomalies based on anomaly probability, can be widely applied in the field of satellite orbital maneuver detection.

Description

Technical Field

[0001] This invention relates to the field of satellite orbital maneuver detection, and in particular to a method and system for detecting orbital anomalies based on anomaly probability. Background Technology

[0002] In today's world, where space security is increasingly valued, the ever-growing number of spacecraft in orbit, especially non-cooperative ones, places higher demands on space situational awareness capabilities. Maintaining real-time tracking of space targets is crucial for improving Earth orbit situational awareness (SSA) capabilities, and timely detection of satellite maneuvers is key to achieving real-time tracking. However, detecting satellite maneuvers, especially in real-time, is an extremely challenging task. Unknown maneuvers by non-cooperative targets could be caused by routine orbital maintenance, orbital deviations due to collisions with space debris, or even orbital interception aimed at military strikes. Rapidly identifying maneuvers and determining the intent of the target is essential to allowing sufficient time for decision-making and response, thus protecting the safety of spacecraft.

[0003] The main challenge in rapid detection of orbital maneuvers is that, under limited thrust, a satellite's orbital state changes very little in a short period, especially for satellites in medium and high orbits. Given the limited accuracy of the data, these changes are often overwhelmed by noise. To address this issue, a certain amount of time is usually required to accumulate data to obtain state changes that are not overwhelmed by noise, but this does not meet the real-time requirements. Summary of the Invention

[0004] To address the aforementioned technical problems, the purpose of this invention is to provide a method and system for detecting track anomalies based on anomaly probability, which can quickly and accurately determine whether an anomaly has occurred in the track.

[0005] The first technical solution adopted in this invention is: a method for detecting orbital anomalies based on anomaly probability, comprising the following steps:

[0006] Acquire satellite orbit data;

[0007] Based on the wavelet analysis method combined with forecasts, satellite orbit data is spliced ​​and reconstructed to obtain the current orbit inclination and semi-major axis information after wavelet transformation.

[0008] Based on the current orbital inclination and semi-major axis information after wavelet transform, sub-criteria and comprehensive criteria are calculated and combined with the limit threshold to perform anomaly detection and obtain the detection results.

[0009] By introducing anomaly probability to map the comprehensive criterion, the probability range of the detection result is obtained.

[0010] Furthermore, it also includes:

[0011] Calculate the change in orbital inclination based on satellite orbital data;

[0012] When the change in track inclination exceeds a preset threshold, it is determined to be track anomaly.

[0013] This optimized procedure allows for a quick determination of whether there has been any significant abnormal change in the track's inclination angle.

[0014] Furthermore, the step of stitching and reconstructing satellite orbit data based on the wavelet analysis method combined with forecasts to obtain the current orbit inclination and semi-major axis information after wavelet transform specifically includes:

[0015] Based on satellite orbit data, read current and historical data;

[0016] Based on the data at the current moment, a forecasting model is used to make a forecast, and the predicted orbit data is obtained.

[0017] The current data, historical data, and predicted orbit data are stitched together to obtain mixed data;

[0018] The mixed data is decomposed using the wavelet transform method to obtain the decomposed data.

[0019] After removing the high-frequency components from the decomposed data, the current orbital inclination and semi-major axis information are reconstructed using wavelet transform.

[0020] Through this optimized step, the wavelet analysis signal processing method combined with forecast can effectively reduce the inundation effect of noise on the abnormal signal.

[0021] Furthermore, the sub-criteria include the change in the semi-major axis and the change in the orbital inclination.

[0022] Furthermore, the calculation formula for the comprehensive criterion is as follows:

[0023]

[0024] In the above formula, |Δa| k |Δi| represents the change in the semi-major axis. k This represents the change in orbital inclination, with δ1 and δ2 representing the corresponding threshold values.

[0025] Furthermore, the step of mapping the anomaly probability to the comprehensive criterion to obtain the probability interval of the detection result specifically includes:

[0026] A quadratic function is used to construct the functional relationship between the probability of anomalies and the comprehensive criterion;

[0027] Calculate the mean and variance of the comprehensive criterion;

[0028] The function is segmented based on the mean and variance of the comprehensive criteria to obtain the probability interval of the detection result.

[0029] The second technical solution adopted in this invention is: a track anomaly detection system based on anomaly probability, comprising:

[0030] The data acquisition module is used to acquire satellite orbit data;

[0031] By combining the wavelet analysis module of the forecast, the satellite orbit data is stitched together and reconstructed based on the wavelet analysis method combined with the forecast, and the orbit inclination and semi-major axis information at the current moment after wavelet transformation are obtained.

[0032] The detection module is used to calculate the sub-criteria and comprehensive criteria and combine them with the limit threshold to perform anomaly detection and obtain the detection results;

[0033] The probability calculation module is used to map the anomaly probability to the comprehensive criterion to obtain the probability range of the detection result.

[0034] The beneficial effects of the method and system of this invention are as follows: By combining the wavelet analysis signal processing method of forecast, this invention can effectively reduce the influence of noise on the submersion of abnormal signals. At the same time, by adopting a comprehensive criterion method with multiple factors, the probability of missed detection and false detection can be reduced. In addition, the concept of abnormal probability is introduced, which transforms the previous judgment of maneuver from a simple "0-1" whether judgment into a probability judgment on the "0-1" interval, reducing the possible impact of false alarms. Attached Figure Description

[0035] Figure 1 This is a flowchart of the steps of a track anomaly detection method based on anomaly probability according to the present invention;

[0036] Figure 2 This is a schematic diagram of the wavelet transform decomposition process;

[0037] Figure 3 This is a schematic diagram of the reconstruction result obtained from wavelet transform decomposition;

[0038] Figure 4 This is a schematic diagram of the wavelet analysis method combined with forecasting in this invention;

[0039] Figure 5 It is the curve relationship between the probability of anomalies and the comprehensive statistics;

[0040] Figure 6 This is a structural block diagram of a track anomaly detection system based on anomaly probability according to the present invention. Detailed Implementation

[0041] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. The step numbers in the following embodiments are only for ease of explanation and do not limit the order of the steps. The execution order of each step in the embodiments can be adapted according to the understanding of those skilled in the art.

[0042] like Figure 1 As shown, this invention provides a method for detecting orbital anomalies based on anomaly probability, the method comprising the following steps:

[0043] S1. Obtain satellite orbit data;

[0044] S2.1 Calculate the change in orbital inclination based on satellite orbital data;

[0045] Satellite orbit data preprocessing;

[0046] To detect orbital anomalies in real time, a segment of orbital data from a normally operating satellite is needed to generate threshold values. The specific process is as follows:

[0047] 1) Select a continuous time period of track data as the criterion data segment. Generally, the total number of data points should be greater than 100.

[0048] 2) Remove abnormal data. This means ensuring that the data used comes entirely from stably operating Kepler orbit satellites. This step usually needs to be done manually.

[0049] 3) Assume that after screening, m sets of data are ultimately obtained that can serve as the standard segment for judging track anomalies. Let the three sub-thresholds δ1 and δ2 be the maximum values ​​of Δa and Δi among all data points within the standard segment, respectively:

[0050]

[0051] Where Δa is the change in the semi-major axis, Δi is the change in the orbital inclination, and j = 1, 2, ..., m.

[0052] In satellite orbit determination, only the satellite's velocity and position information at a certain moment is needed to determine the satellite's orbital elements. Let's denote the position vector as r and the velocity vector as v.

[0053] First, calculate the orbital angular momentum h:

[0054] h = r × v

[0055] Calculate the right ascension vector N of the ascending node:

[0056] N = K × h

[0057] Where K = [0, 0, 1] T It is the unit vector of the Z-axis in the ECI coordinate system.

[0058] Calculate the eccentricity vector e and the magnitude of eccentricity e:

[0059]

[0060] e = e

[0061] Where μ = 398600 * 10 9 m 3 / s 2 It is the Earth's gravitational constant.

[0062] Calculate the semi-pipe diameter p:

[0063]

[0064] Calculate the semi-major axis a:

[0065]

[0066] Calculate the orbital inclination angle i:

[0067]

[0068] Calculate the right ascension of the ascending node:

[0069]

[0070] Where N represents the cross product of h and the z-axis unit vector (N = [0; 0; 1] × h), N x The magnitude of the vector projected onto the x-axis.

[0071] Calculate the argument of perigee:

[0072]

[0073] Calculate the true anterior angle:

[0074]

[0075] In particular, if the eccentricity is close to 0 (e < 10) -14 If the perigee argument is:

[0076] w=0

[0077] The true nearest angle is:

[0078]

[0079] This gives us the orbital elements of the satellite at the current moment.

[0080] S2.2 When the change in track inclination exceeds a preset threshold, it is determined to be track anomaly.

[0081] Preliminary determination of track inclination;

[0082] By using three consecutive position vectors, it is possible to quickly determine whether there has been a significant change in the tilt angle of the satellite orbit.

[0083] h1 = r1 × r2

[0084] h2=r2×r3

[0085]

[0086] Compare sin(Δi) with a preset tilt angle change threshold. If the value is greater than the threshold, the orbit can be considered to have shifted. The threshold is calculated using a preset initial value and historical data.

[0087] S3. Based on the wavelet analysis method combined with the forecast, the satellite orbit data is spliced ​​and reconstructed to obtain the current orbit inclination and semi-major axis information after wavelet transformation;

[0088] Wavelet analysis is a post-processing analysis method that requires a certain amount of data accumulation to perform effective analysis. Furthermore, the reliability of denoising results is often higher in the middle of the data stream. In inter-satellite maneuver detection, some high-precision, high-speed sensors can accumulate sufficient data in a short time, but this is difficult to achieve with other data sources. During satellite tracking, if a satellite maneuvers and data updates are not timely enough, effective tracking may fail, leading to target loss. By predicting the orbit and combining this prediction information with wavelet analysis, and applying it to real-time scenarios, relatively ideal results can be achieved through simulation analysis.

[0089] Wavelet transform is a time-scale (time-frequency) analysis method for signals. It features multi-resolution analysis and can more effectively identify abrupt changes in signals, extracting dynamic signals from noise.

[0090] The precise definition of the wavelet function is as follows:

[0091] Let ψ(t) be a square-integrable function. If its Fourier transform ψ(w) satisfies the following condition:

[0092]

[0093] Then ψ(t) is called the wavelet mother function. After scaling and translating it, we obtain the wavelet basis functions ψ. a,b (t), any belonging to L 2 The continuous wavelet transform of the signal s(t) of (R) is:

[0094]

[0095] In the formula: a is the scaling factor, b is the translation factor, and "<>" is the inner product symbol. To facilitate computer computation, the continuous wavelet transform needs to be discretized. This discretization mainly targets the continuous scaling parameter a and the continuous translation parameter b. The most commonly used is a binary dynamic sampling network: a0 = 2, b0 = 1. The binary wavelet transform is obtained as follows:

[0096]

[0097] Among them, ψ j,k (t)=2 -j / 2 ψ(2 -j tk); m,n∈N

[0098] By employing a multi-level wavelet transform at a selected scale, the signal is decomposed into a low-frequency component and several high-frequency components. For the signal decomposition of satellite orbital elements, the high-frequency components primarily represent noise information, while the low-frequency components mainly contain abrupt changes in motion. By appropriately processing the high-frequency components, noise signals can be suppressed, and abrupt changes can be extracted. The wavelet transform decomposition process is as follows: Figure 2 As shown, CA represents the low-frequency component and CD represents the high-frequency component. The satellite orbit is decomposed into three levels of wavelet decomposition to obtain signals at different levels. The reconstruction results are shown below. Figure 3 As shown:

[0099] S3.1. Based on satellite orbit data, read the current time data and historical data;

[0100] S3.2. Based on the current time data, use the forecasting model to make a forecast and obtain the predicted orbit data;

[0101] S3.3. Combine the current data, historical data, and predicted orbit data to obtain mixed data;

[0102] Because wavelet analysis suffers from a "boundary problem," abrupt changes near the endpoints are difficult to identify effectively. Therefore, wavelet transform requires data accumulation from future timeframes to better extract the abrupt changes at the current moment, which severely impacts its application in real-time maneuver signal extraction.

[0103] To address this, a wavelet analysis method combining forecasting is proposed. Specifically, orbit forecasting is performed at the current time point, and wavelet transform decomposition at the current point is achieved by splicing a segment of historical data with a segment of forecast data. This processing procedure refers to... Figure 4 .

[0104] S3.4. Decompose the mixed data based on the wavelet transform method to obtain the decomposed data;

[0105] S3.5 Remove the high-frequency components from the decomposed data and reconstruct the current orbital inclination and semi-major axis information after wavelet transform.

[0106] S4. Calculate the sub-criteria and comprehensive criteria based on the current orbital inclination and semi-major axis information after wavelet transform, and combine them with the limit threshold to perform anomaly detection and obtain the detection results;

[0107] The sub-criteria include the change in semi-major axis and the change in orbital inclination.

[0108] When a satellite undergoes anomalies, parameters such as orbital elements often change accordingly. Among these variables, Δa and Δi are relatively stable when there are no anomalies, but undergo significant abrupt changes when there are anomalies. Therefore, they are selected as the sub-criteria for anomaly identification.

[0109] For a certain time k, the sub-criterion is calculated as follows:

[0110] 1) Change in semi-major axis |Δa|:

[0111] |Δa|=|a k -a k-1 |

[0112] 2) Calculate the change in track inclination |Δi|:

[0113] |Δi|=|i k -i k-1 |

[0114] The calculation formula for the comprehensive criterion is as follows:

[0115]

[0116] In the above formula, |Δa| k |Δi| represents the change in the semi-major axis. k This represents the change in orbital inclination, with δ1 and δ2 representing the corresponding threshold values.

[0117] The comprehensive criterion is compared with the critical value 1. If it is greater than the critical value 1, it is considered that an anomaly has occurred at that moment; if it is less than the critical value 1, it is considered that no anomaly has occurred at that moment.

[0118] S5. Introduce the anomaly probability to map the comprehensive criterion and obtain the probability interval of the detection result.

[0119] S5.1. Construct the functional relationship between the probability of movement and the comprehensive criterion using a quadratic function;

[0120] S5.2 Calculate the mean and variance of the comprehensive criterion;

[0121] S5.3. Based on the mean and variance of the comprehensive criterion, the function is segmented to obtain the probability interval of the detection result.

[0122] Considering that the orbital inclination increment, semi-major axis increment, and average velocity dispersion are relatively stable when the target does not undergo maneuvering, and assuming that the observation error follows a Gaussian distribution with a mean of 0, we propose to use a Gaussian distribution to describe the distribution characteristics of these variables, and the 3σ criterion can be used to determine the threshold for orbital maneuvering.

[0123] By statistically analyzing the mean and variance of the comprehensive criterion quantity under the condition of no maneuver, if the comprehensive criterion quantity falls near the mean, the probability of no target maneuver is highest. Conversely, if it deviates significantly from the mean, it indicates the presence of some kind of maneuver or intervention, causing it to deviate from its original stable state, which means the probability of the target maneuver increases. Based on this approach, a mapping expression for the maneuver probability value with respect to the comprehensive criterion quantity was designed.

[0124]

[0125] In fact, the change in maneuver probability is the same whether the statistic X moves away from the mean μ in the positive or negative direction. Furthermore, quadratic functions have the characteristic that the derivative gradually increases and the amplitude changes slowly at first and then rapidly as it moves away from the axis of symmetry. Therefore, a quadratic function is proposed to describe the functional relationship between the maneuver probability value and the comprehensive criterion quantity. Further, based on the mean μ and variance σ of the comprehensive criterion quantity... 2 Statistically, functions with X > μ are divided into four segments: 1) [μ, μ+σ); 2) [μ+σ, μ+2σ); 3) [μ+2σ, μ+3σ); 4) [μ+3σ, +∞). Based on the symmetry of quadratic functions, functions with X ∈ [0, μ) can also be determined accordingly. The specific relational expression for the function f1(X) is as follows:

[0126]

[0127] In the region [μ, μ+σ), the statistic X falls near the mean, indicating a low probability of target maneuvering. However, even when X = μ, the probability of maneuvering cannot be completely considered zero; a lower probability value P should be set. (1) When X = μ + σ, a lower probability of maneuver P can generally be set. (2) And P (2) >P (1) Therefore, the piecewise functional expression for [μ, μ+σ) can be solved by the following system of equations:

[0128]

[0129] The following matrix equation is obtained by rearranging:

[0130]

[0131] It can be seen that when P( 1 ) and P( 2 When P( 1 ) and P( 2 The value of P( ) can be determined based on the specific tracking environment and further engineering experience. Here, we will temporarily set P( ) as . 1 ) and P( 2 Setting the values ​​to 0.1 and 0.2, we can solve for...

[0132]

[0133] Secondly, in the region [μ+σ, μ+2σ), the statistic is slightly farther from the mean, and should have a slightly higher maneuver probability value than in the region [μ, μ+σ). Furthermore, as the statistic X gradually increases, the rate of increase in the maneuver probability also gradually increases. That is, when X = μ+2σ, the maneuver probability value P... (3) P should be satisfied (3) >P (2) +(P (2) -P (1) Similarly, the relational expression for the region [μ+σ, μ+2σ) can be obtained as follows:

[0134]

[0135] The following matrix equation can be obtained by rearranging:

[0136]

[0137] If P (2) and P (3) Setting the values ​​to 0.2 and 0.5 respectively, we can solve for...

[0138]

[0139] Similarly, when the statistic X falls within the region [μ+2σ, μ+3σ), it means the target is highly likely to maneuver and should have a high probability of maneuvering. Specifically, when X = 3σ, it means the target is extremely likely to maneuver, and should have a particularly high probability of maneuvering P. (4) Generally, P (4) We can approximate it as 1. Then the function expression for this region can be calculated as follows:

[0140]

[0141] The following matrix equation can be obtained by rearranging:

[0142]

[0143] When P( 4When ) = 1, we can solve for ,

[0144]

[0145] When the statistic X falls into the region [μ+3σ,+∞), it is considered that the target has maneuvered, and the probability of maneuvering is 1.

[0146] Historical TLE data for satellite 41028U in 2021 were selected. The mean and variance of the semi-major axis change Δa and the orbital inclination change Δi were statistically analyzed. Among these, μ... Δa =64.70056, μ Δi =0.0018,

[0147] The mean μ of the composite statistic X is obtained from the two distributions in the following manner. final and variance

[0148]

[0149]

[0150] Where ω1=ω2=0.5, μ is calculated. final =0.0896, Substitute it into the function expression for each segment and plot the function graph as follows: Figure 5 It can be seen that the probability of maneuvering is minimum at the mean, and the further away from the mean, the greater the probability of maneuvering, and the greater the slope.

[0151] like Figure 6 As shown, a trajectory anomaly detection system based on anomaly probability includes:

[0152] The data acquisition module is used to acquire satellite orbit data;

[0153] By combining the wavelet analysis module of the forecast, the satellite orbit data is stitched together and reconstructed based on the wavelet analysis method combined with the forecast, and the orbit inclination and semi-major axis information at the current moment after wavelet transformation are obtained.

[0154] The detection module is used to calculate the sub-criteria and comprehensive criteria and combine them with the limit threshold to perform anomaly detection and obtain the detection results;

[0155] The probability calculation module is used to map the anomaly probability to the comprehensive criterion to obtain the probability range of the detection result.

[0156] The content of the above method embodiments is applicable to this system embodiment. The specific functions implemented in this system embodiment are the same as those in the above method embodiments, and the beneficial effects achieved are also the same as those achieved in the above method embodiments.

[0157] A track anomaly detection device based on anomaly probability:

[0158] At least one processor;

[0159] At least one memory for storing at least one program;

[0160] When the at least one program is executed by the at least one processor, the at least one processor implements the orbital anomaly detection method based on anomaly probability as described above.

[0161] The content of the above method embodiments is applicable to the device embodiments. The specific functions implemented by the device embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those achieved by the above method embodiments.

[0162] A storage medium storing processor-executable instructions, characterized in that: the processor-executable instructions, when executed by the processor, are used to implement the orbital anomaly detection method based on anomaly probability as described above.

[0163] The content of the above method embodiments is applicable to this storage medium embodiment. The specific functions implemented in this storage medium embodiment are the same as those in the above method embodiments, and the beneficial effects achieved are also the same as those achieved in the above method embodiments.

[0164] The above is a detailed description of the preferred embodiments of the present invention. However, the present invention is not limited to the embodiments described. Those skilled in the art can make various equivalent modifications or substitutions without departing from the spirit of the present invention. All such equivalent modifications or substitutions are included within the scope defined by the claims of this application.

Claims

1. A method for detecting orbital anomalies based on anomaly probability, characterized in that, Includes the following steps: Acquire satellite orbit data; Based on the wavelet analysis method combined with forecasts, satellite orbit data is spliced ​​and reconstructed to obtain the current orbit inclination and semi-major axis information after wavelet transformation. Based on the current orbital inclination and semi-major axis information after wavelet transform, sub-criteria and comprehensive criteria are calculated and combined with the limit threshold to perform anomaly detection and obtain the detection results. By introducing anomaly probability to map the comprehensive criterion, the probability range of the detection result is obtained; The step of stitching and reconstructing satellite orbit data based on the wavelet analysis method combined with forecasts to obtain the current orbit inclination and semi-major axis information after wavelet transform specifically includes: Based on satellite orbit data, read current and historical data; Based on the data at the current moment, a forecasting model is used to make a forecast, and the predicted orbit data is obtained. The current data, historical data, and predicted orbit data are stitched together to obtain mixed data; The mixed data is decomposed using the wavelet transform method to obtain the decomposed data. After removing the high-frequency components from the decomposed data, the current orbital inclination and semi-major axis information are reconstructed using wavelet transform.

2. The orbital anomaly detection method based on anomaly probability according to claim 1, characterized in that, Also includes: Calculate the change in orbital inclination based on satellite orbital data; When the change in track inclination exceeds a preset threshold, it is determined to be track anomaly.

3. The orbital anomaly detection method based on anomaly probability according to claim 1, characterized in that, The sub-criteria include the change in semi-major axis and the change in orbital inclination.

4. The orbital anomaly detection method based on anomaly probability according to claim 3, characterized in that, The formula for calculating the comprehensive criterion is as follows: In the above formula, |Δa| k |Δi| represents the change in the semi-major axis. k This represents the change in orbital inclination, with δ1 and δ2 representing the corresponding threshold values.

5. The orbital anomaly detection method based on anomaly probability according to claim 4, characterized in that, The step of mapping the anomaly probability to the comprehensive criterion to obtain the probability interval of the detection result specifically includes: A quadratic function is used to construct the functional relationship between the probability of anomalies and the comprehensive criterion; Calculate the mean and variance of the comprehensive criterion; The function is segmented based on the mean and variance of the comprehensive criteria to obtain the probability interval of the detection result.

6. A track anomaly detection system based on anomaly probability, characterized in that, A method for performing orbital anomaly detection based on anomaly probability as described in claim 1 includes: The data acquisition module is used to acquire satellite orbit data; By combining the wavelet analysis module of the forecast, the satellite orbit data is stitched together and reconstructed based on the wavelet analysis method combined with the forecast, and the orbit inclination and semi-major axis information at the current moment after wavelet transformation are obtained. The detection module is used to calculate the sub-criteria and comprehensive criteria and combine them with the limit threshold to perform anomaly detection and obtain the detection results; The probability calculation module is used to map the anomaly probability to the comprehensive criterion to obtain the probability range of the detection result.

Images