Non-cooperative spacecraft continuous low-thrust maneuver detection and inversion method

Abstract

The application discloses a non-cooperative spacecraft continuous low-thrust maneuver detection and inversion method and belongs to the technical field of aerospace. The method comprises the following steps: a dynamic basic model containing unknown low-thrust input is given; the continuous low-thrust acceleration is converted into Chebyshev coefficient form; joint estimation of the orbit state and the maneuver coefficient is completed under an extended state filtering framework; a statistical baseline is calibrated and maneuver determination is triggered based on a non-maneuver reference trajectory; only the maneuver coefficient covariance sub-block is iteratively inflated when the detection is established and the observation residual is over-limit; the thrust acceleration norm is reconstructed according to the compensated maneuver coefficient, and the final maneuver period is extracted through an output switch function. The non-cooperative spacecraft continuous low-thrust maneuver detection and inversion method can realize online detection of continuous low-thrust maneuver, thrust acceleration inversion and maneuver period output.

Description

Technical Field

[0001] This invention relates to the field of aerospace technology, and in particular to a method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft. Background Technology

[0002] Non-cooperative spacecraft are characterized by limited prior information, unknown orbital maneuver patterns, and difficulty in directly acquiring on-orbit behavior. Their maneuver detection and future trajectory assessment are key issues in space situational awareness. Compared to pulsed maneuvers, continuous low-thrust maneuvers are characterized by small thrust amplitude, long duration, slow changes in direction, and often non-periodic characteristics. This can easily lead to systematic deviations in orbit determination results based solely on nominal dynamic models, thus affecting the determination of maneuver periods and thrust parameter inversion. Existing methods for continuous low-thrust maneuver detection and inversion mainly include multi-model detection and basis function parameterization approaches. The former typically relies on a pre-built maneuver model library, achieving state estimation and maneuver determination through model switching; the latter usually expands the unknown thrust into a set of basis function coefficients and jointly estimates these coefficients with the orbital state. When processing mixed angle, range, and velocity measurement observation data, these methods often struggle to simultaneously balance model universality, computational complexity, detection timeliness, and thrust amplitude inversion accuracy. Especially under unknown continuous low-thrust maneuvering conditions, if a fixed prior model or direct parameter incremental correction method is directly adopted, problems such as unstable detection threshold, absorption of maneuvering coefficient estimation by track state, and overall divergence of filter covariance due to compensation link are likely to occur.

[0003] In existing technologies, variable-dimensional filtering and multi-model switching techniques are used to achieve maneuver detection and state estimation for the orbit estimation problem of continuous thrust spacecraft. While this approach maintains good detection performance when the model library is sufficient, the complexity of model design and computation increases significantly with the number of candidate maneuver models, and its adaptability to unknown maneuver patterns is limited. Another approach for determining the orbit of unknown continuous maneuvers uses thrust Fourier coefficients to extend the state characterization of unknown thrust and directly estimates maneuver parameters using an extended Kalman filter. However, while this approach reduces reliance on prior maneuver models, it is prone to truncation errors and coefficient leakage in the Fourier basis functions when dealing with non-periodic, locally abrupt, or short-segment continuous thrust maneuvers, making it difficult to stably recover the true continuous low thrust amplitude.

[0004] Therefore, there is a need for an integrated method that can maintain filtering stability while taking into account maneuver detection, equivalent thrust acceleration norm inversion, and maneuver period output. Summary of the Invention

[0005] The purpose of this invention is to provide a method for detecting and inverting continuous low-thrust maneuvers in non-cooperative spacecraft. This method is used to achieve online detection, thrust acceleration inversion, and output of maneuver periods for continuous low-thrust maneuvers, under the condition that only observation data and rough initial orbit values ​​are known. It utilizes Chebyshev polynomial extended state, no-maneuver reference threshold calibration, and maneuver coefficient covariance expansion compensation mechanism.

[0006] To achieve the above objectives, this invention provides a method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft, comprising the following steps: S1. Based on the rough initial orbit values ​​of the target spacecraft and the observation spacecraft, establish a spacecraft dynamics model that considers two-body gravity, second-order spherical harmonic perturbation and unknown continuous low thrust input, and obtain the basic state equations. S2. Based on the basic state equations obtained in step S1, the unknown continuous low thrust acceleration in each sampling interval is mapped to the standard interval, and the second kind of Chebyshev polynomial is used for piecewise expansion to obtain the parameterized expression of the maneuver coefficient. S3. Based on the parameterized expression of the maneuver coefficient in step S2, combine the maneuver coefficient vector with the orbital state to form a generalized state and establish an extended state transition relationship to obtain the extended system prediction model. S4. Based on the extended system prediction model, establish a multi-observer dual-angle, distance, and relative velocity hybrid observation model and perform extended Kalman filter prediction update to obtain the joint estimation results of the orbital state and maneuver coefficients of the current sampling step. S5. Based on the joint estimation results, the non-maneuvering reference trajectory is pre-run and the detection statistical baseline is calibrated. Then, the maneuvering statistics in the real-time observation branch are smoothed and thresholded to obtain the maneuvering detection results. S6. Based on the maneuver detection results and the current sampling step observation residual exceeding the threshold, perform iterative expansion compensation on the maneuver coefficient covariance sub-block and re-complete the prediction update for the current sampling step to obtain the corrected maneuver coefficient estimation result. S7. Based on the estimation results of the corrected maneuver coefficient, reconstruct the thrust acceleration norm and construct the output switching function, extract the final continuous low thrust maneuver period and output the thrust acceleration inversion results.

[0007] Preferably, step S1 specifically includes: In a geocentric inertial frame, the spacecraft's state vector is denoted as a six-dimensional vector composed of position and velocity. The unknown continuous low-thrust acceleration is introduced as an external input into the state equation, resulting in: ; in, and These represent the spacecraft's position vector and velocity vector, respectively. This represents the gravitational acceleration due to two bodies. This represents the second-order spherical harmonic perturbation acceleration. Assign a matrix to the input. The unknown continuous low-thrust acceleration is represented by the following form: Second-order spherical harmonic perturbation acceleration is expressed as follows: ; in, The gravitational constant of Earth, These are the second-order spherical harmonic coefficients. The radius of the Earth's equator. The distance from the spacecraft to the Earth's center. The component represents the position; the output is the basic dynamic model under continuous low thrust maneuvering conditions, for subsequent parameterization and filtering.

[0008] Preferably, step S2 specifically includes: Regarding the first sampling intervals First, map physical time to a standard interval. ,get: ; The second kind of Chebyshev polynomial takes the following recursive form: ; In the orbital coordinate system, the unknown acceleration is decomposed along the radial, tangential, and normal directions, and expanded using Chebyshev basis functions: ; in, , and Let the radial, tangential, and normal unit basis vectors of the orbital coordinate system be represented respectively, and their definitions be: ; in, , and These represent the Chebyshev expansion coefficients along the three axes, respectively. To expand the order, the above steps transform the unknown continuous low thrust acceleration within each sampling interval into a finite-dimensional coefficient estimation problem.

[0009] Preferably, in step S3, the generalized state vector is formed by combining the maneuver coefficient vector and the orbital state: ; in, For the maneuver coefficient vector, , , Let Chebyshev expand the coefficient subvectors. , , These represent the directions in the orbital coordinate system; Since the Chebyshev coefficients within each sampling interval are treated as constant parameters, the extended system dynamics can be expressed as: ; Linearizing the extended system dynamics at the current estimation point yields the variational equation for the state transition matrix: ; in, Describing extended dynamics with respect to generalized states Jacobian matrix, Represents the identity matrix. Indicates in The state transition matrix of each sampling step is obtained by numerical integration. This is used for subsequent filtering and prediction.

[0010] Preferably, step S4 specifically includes: For each observed spacecraft, the measurement vector is composed of observations of two angles, distance, and relative velocity. The overall observation model is represented as follows: ; The observation function for a single observer is taken as: ; Each component is determined by the relative position and relative velocity of the target spacecraft and the observation spacecraft: ; ; By performing first-order linearization on the observation function, we obtain the observation Jacobian matrix. And perform extended Kalman filter prediction updates: ; ; in, Indicates the first Each sampling step, From the state transition matrix Give, For process noise covariance, To measure the noise covariance, To measure the noise vector, The elevation angle of the target relative to the line of sight of the observing spacecraft. The azimuth angle of the target relative to the line of sight of the observation spacecraft. The relative distance between the target and the observed spacecraft. The magnitude of the relative velocity of the target with respect to the observed spacecraft. , , These are the components of the relative position vector. Let V be the relative velocity vector of the target spacecraft with respect to the observation spacecraft. For the first The posterior state covariance matrix updated after each sampling step To change the state from the first Spread to the Nonlinear state transition function for each sampling step, These are the prior generalized state estimates. For the first The updated posterior generalized state estimate after each sampling step Let be the prior state covariance matrix. Here is the Kalman gain matrix. This is the updated posterior generalized state estimate. This is the predicted measurement vector mapped from the prior state estimate. To ensure that the maneuver information in the current sampling interval is primarily driven by the current observation, the prior mean of the maneuver coefficient is reset to zero after each sampling step prediction, while maintaining a constant prior covariance.

[0011] Preferably, step S5 obtains the non-maneuvering reference trajectory based on the coarse initial value integration, and performs a pre-run on the reference observation using a model consistent with the real-time filtering. The norm of the updated maneuvering coefficients at each step is defined as: ; Before pre-run Establish statistical baselines on each sample: ; In the real-time observation branch, the maneuver statistics are exponentially smoothed, and the presence of maneuver is determined accordingly. ; in, For smoothing coefficients, To detect the threshold coefficient, This is the vector of maneuver coefficients obtained through filtering estimation. The baseline of the mean of the mobility statistics for the sample of vehicles without pre-operational mobility is given. For the first phase of non-motorized pre-operation The mobility statistic corresponding to each sample The standard deviation of the motor statistics for the pre-operation sample without motors. This is the exponentially smoothed value of the mobility statistic. The maneuver statistic is defined as the L2 norm of the maneuver coefficient estimation vector. For the first The exponential smoothing value of the motion statistic for each sampling step is used. When the smoothing statistic exceeds a threshold, the sampling interval is determined to be in a motion state. The output of this step is used both for subsequent compensation triggering and for constructing prior criteria for the output of the motion period.

[0012] Preferably, in step S6, after the detection branch gives the maneuver determination, the compensation branch decides whether to initiate covariance expansion based on the observation residual, defining the real-time observation residual and component threshold as follows: ; ; When there is a maneuver and at least one observation component satisfies When the value is greater than zero, it is considered that the maneuver coefficient covariance sub-block needs to be compensated. During the compensation iteration initialization, the orbital state covariance sub-block is not amplified; only the maneuver coefficient covariance sub-block is set to the dilated prior scale, and the prediction update of the current sampling step is re-executed. ; in, Indicates the first The covariance sub-block of the mobility coefficients in the next compensation iteration For the prior measure of the mobility coefficient, It is the expansion factor. The expansion growth rate, To measure the residual vector, This is the residual judgment threshold vector. This is the residual threshold coefficient, used to amplify the boundary of the measurement residual. After each expansion, the prediction update of the current sampling step is re-executed until all observation components do not exceed the threshold or the maximum number of iterations is reached.

[0013] Preferably, in step S7, for each sampling interval, select within the standard interval. Averaging the Chebyshev expansion results at discrete points along the three axes yields a piecewise axial acceleration estimate; further, the thrust acceleration norm is defined: ; ; Construct an output switching function to extract the maneuver period from the thrust acceleration inversion curve: ; in, Extract the proportional coefficient for the output time period. For the first Each sampling interval, along The first axial direction Estimates of the second-type Chebyshev expansion coefficients For the first Discrete points of a second-order Chebyshev polynomial in the standard interval The value at the position, where the subscript is in the formula Indicates the index of discrete points within a standard interval. For the first The thrust acceleration norm estimate for each sampling interval , , For the first The average maneuvering acceleration estimate of each component in each sampling interval The first of the thrust acceleration norm sequences in the entire sampling interval Item; based on the output switching function The non-zero interval is used to obtain the final maneuvering time period, which is used for plotting the switching function and statistical analysis of the implementation results in the attached figure.

[0014] Therefore, the present invention employs the above-mentioned method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft, which has the following beneficial effects: (1) By performing a partial expansion of the unknown continuous low thrust acceleration using the second type of Chebyshev polynomial, it is possible to characterize non-periodic and time-varying continuous low thrust maneuvers without needing to know the maneuver law in advance, thereby improving the applicability of unknown maneuver scenarios. (2) By using the pre-run calibration of the non-maneuvering reference trajectory to detect the statistical baseline, and by using the smoothing statistics of the maneuvering coefficient norm for detection, it is possible to achieve adaptive maneuvering judgment that matches the current error condition under mixed angle measurement, distance measurement and velocity measurement observation conditions, thereby reducing the risk of false detection and missed detection caused by fixed threshold. (3) By performing iterative expansion compensation only on the covariance sub-block of the maneuver coefficient, the overall over-amplification of the orbital state covariance is avoided. It can effectively recover the true thrust amplitude while maintaining the stability of the filter, and directly give the continuous low thrust maneuver period through the output switching function, thereby improving the availability of detection and inversion results.

[0015] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0016] Figure 1 This is a flowchart of a method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to the present invention; Figure 2 This is a graph showing the relationship between smooth detection statistics and judgment threshold under condition N in Embodiment 5 of the present invention; Figure 3 This is a diagram showing the thrust acceleration norm inversion result under the N-condition in Embodiment 5 of the present invention; Figure 4This is a comparison diagram of the output switching function under condition N in Embodiment 5 of the present invention and the actual maneuvering range. Detailed Implementation

[0017] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0018] based on Figure 1 The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft, as shown in this embodiment, only employs the covariance expansion compensation route to illustrate the detection, compensation, inversion, and time-period output process under 5 N continuous low-thrust conditions. Specifically, it includes: S1. Based on the initial orbital values ​​of the target spacecraft with an initial epoch of 2025-01-01 00:00:00 UTC and three observation spacecraft, a spacecraft dynamic model considering two-body gravity, second-order spherical harmonic perturbation, and unknown low thrust input is established. The actual maneuvering period of 900 s to 2700 s, thrust of 5 N, initial mass of 1000 kg, and specific impulse of 2800 s are used as the actual reference conditions for subsequent parameter inversion.

[0019] In a geocentric inertial frame, the spacecraft's state vector is denoted as a six-dimensional vector composed of position and velocity. The unknown continuous low-thrust acceleration is introduced as an external input into the state equation, resulting in: ; in, and These represent the spacecraft's position vector and velocity vector, respectively. This represents the gravitational acceleration due to two bodies. This represents the second-order spherical harmonic perturbation acceleration. Assign a matrix to the input. The unknown continuous low-thrust acceleration is represented by the following form: Second-order spherical harmonic perturbation acceleration is expressed as follows: ; in, The gravitational constant of Earth, These are the second-order spherical harmonic coefficients. The radius of the Earth's equator. The distance from the spacecraft to the Earth's center. The component represents the position; the output is the basic dynamic model under continuous low thrust maneuvering conditions, for subsequent parameterization and filtering.

[0020] In this embodiment, the initial orbital parameters of the target spacecraft are: semi-major axis 43,000 km, eccentricity 0.01, inclination 2°, right ascension of ascending node 0°, argument of perigee 0°, and mean perigee 15°; the actual thrust direction is taken as the positive direction of the orbital normal, therefore the actual peak thrust acceleration is 5.000 × 10⁻⁶. -3 m / s 2 The corresponding value is 5.000 × 10 -6 km / s 2 .

[0021] S2. Based on the basic state equations obtained in step S1, the unknown continuous low thrust acceleration in each sampling interval is mapped to the standard interval, and a piecewise expansion is performed using the second type of Chebyshev polynomial to obtain the parameterized expression of the maneuver coefficient.

[0022] Regarding the first sampling intervals First, map physical time to a standard interval. ,get: ; The second kind of Chebyshev polynomial takes the following recursive form: ; In the orbital coordinate system, the unknown acceleration is decomposed along the radial, tangential, and normal directions, and expanded using Chebyshev basis functions: ; in, , and Let the radial, tangential, and normal unit basis vectors of the orbital coordinate system be represented respectively, and their definitions be: ; in, , and These represent the Chebyshev expansion coefficients along the three axes, respectively. To expand the order, the above steps transform the unknown continuous low thrust acceleration within each sampling interval into a finite-dimensional coefficient estimation problem.

[0023] In this embodiment, the sampling interval is 60 s, the Chebyshev order is 3, and the number of sampling points in each standard interval is 10. Therefore, the number of single-axis expansion coefficients is 4, and a total of 12 maneuver coefficients are obtained in the three axes.

[0024] S3. Based on the parameterized expression of the maneuver coefficient in step S2, combine the maneuver coefficient vector with the orbital state to form a generalized state and establish an extended state transition relationship to obtain the extended system prediction model.

[0025] The generalized state vector is formed by combining the maneuver coefficient vector and the orbital state: ; in, For the maneuver coefficient vector, , , Let Chebyshev expand the coefficient subvectors. , , These represent the directions in the orbital coordinate system; Since the Chebyshev coefficients within each sampling interval are treated as constant parameters, the extended system dynamics can be expressed as: ; Linearizing the extended system dynamics at the current estimation point yields the variational equation for the state transition matrix: ; in, Describing extended dynamics with respect to generalized states Jacobian matrix, Represents the identity matrix. Indicates in The state transition matrix of each sampling step is obtained by numerical integration. This is used for subsequent filtering and prediction.

[0026] In this embodiment, the six-dimensional orbital state is combined with the 12-dimensional maneuver coefficient to obtain the 18-dimensional extended state. Therefore, the current sampling step completes the filtering prediction through the 18×18 state transition matrix.

[0027] S4. Based on the extended system prediction model, establish a multi-observer dual-angle, distance, and relative velocity hybrid observation model and perform extended Kalman filter prediction update to obtain the joint estimation results of the orbital state and maneuver coefficients for the current sampling step.

[0028] For each observed spacecraft, the measurement vector is composed of observations of two angles, distance, and relative velocity. The overall observation model is represented as follows: ; The observation function for a single observer is taken as: ; Each component is determined by the relative position and relative velocity of the target spacecraft and the observation spacecraft: ; ; By performing first-order linearization on the observation function, we obtain the observation Jacobian matrix. And perform extended Kalman filter prediction updates: ; ; in, Indicates the first Each sampling step, From the state transition matrix Give, For process noise covariance, To measure the noise covariance, To measure the noise vector, The elevation angle of the target relative to the line of sight of the observing spacecraft. The azimuth angle of the target relative to the line of sight of the observation spacecraft. The relative distance between the target and the observed spacecraft. The magnitude of the relative velocity of the target with respect to the observed spacecraft. , , These are the components of the relative position vector. Let V be the relative velocity vector of the target spacecraft with respect to the observation spacecraft. For the first The posterior state covariance matrix updated after each sampling step To change the state from the first Spread to the Nonlinear state transition function for each sampling step, These are the prior generalized state estimates. For the first The updated posterior generalized state estimate after each sampling step Let be the prior state covariance matrix. Here is the Kalman gain matrix. This is the updated posterior generalized state estimate. The predicted measurement vector is obtained by mapping the prior state estimate.

[0029] In this embodiment, three spacecraft are observed, with orbital inclinations of 80°, 3°, and 1°, respectively. The remaining major orbital parameters are on the same order of magnitude as the target spacecraft. Each observed spacecraft provides two angle measurements, one distance measurement, and one relative velocity measurement at each sampling time, resulting in a total of 12 observation components per update. The angular observation noise is set to 1 arcsec, the ranging noise to 0.2 m, and the velocity measurement noise to 1 × 10⁻⁶ m. -3 m / s.

[0030] S5. Based on the joint estimation results, the non-maneuvering reference trajectory is pre-run and the detection statistical baseline is calibrated. Then, the maneuvering statistics in the real-time observation branch are smoothed and thresholded to obtain the maneuvering detection results.

[0031] Based on the rough initial value integration, a non-maneuvering reference trajectory is obtained. A pre-run is performed on the reference observations using a model consistent with real-time filtering. The norm of the updated maneuvering coefficients at each step is defined as follows: ; Before pre-run Establish statistical baselines on each sample: ; In the real-time observation branch, the maneuver statistics are exponentially smoothed, and the presence of maneuver is determined accordingly. ; in, For smoothing coefficients, To detect the threshold coefficient, This is the vector of maneuver coefficients obtained through filtering estimation. The baseline of the mean of the mobility statistics for the sample of vehicles without pre-operational mobility is given. For the first phase of non-motorized pre-operation The mobility statistic corresponding to each sample The standard deviation of the motor statistics for the pre-operation sample without motors. This is the exponentially smoothed value of the mobility statistic. The maneuver statistic is defined as the L2 norm of the maneuver coefficient estimation vector. For the first The exponential smoothing value of the motion statistic for each sampling step is used. When the smoothing statistic exceeds a threshold, the sampling interval is determined to be in a motion state. The output of this step is used both for subsequent compensation triggering and for constructing prior criteria for the output of the motion period.

[0032] In this embodiment, the smoothing coefficient is set to 0.20, the median window length is set to 3, and the number of pre-run samples is set to 8, resulting in a baseline mean of 1.088 × 10⁻⁶. -8 Standard deviation 1.437 × 10 -8 The corresponding detection threshold is 5.398×10. -8 .

[0033] like Figure 2 As shown, the smoothing detection statistic rises rapidly after the start of the actual maneuver and remains above the judgment threshold during the main maneuver period.

[0034] S6. Based on the maneuver detection results and the current sampling step observation residual exceeding the threshold, perform iterative expansion compensation on the maneuver coefficient covariance sub-block and re-complete the prediction update for the current sampling step to obtain the corrected maneuver coefficient estimation result.

[0035] After the detection branch provides the maneuver decision, the compensation branch decides whether to initiate covariance expansion based on the observation residuals. The real-time observation residuals and component thresholds are defined as follows: ; ; When there is a maneuver and at least one observation component satisfies When the value is greater than zero, it is considered that the maneuver coefficient covariance sub-block needs to be compensated. During the compensation iteration initialization, the orbital state covariance sub-block is not amplified; only the maneuver coefficient covariance sub-block is set to the dilated prior scale, and the prediction update of the current sampling step is re-executed. ; in, Indicates the first The covariance sub-block of the mobility coefficients in the next compensation iteration For the prior measure of the mobility coefficient, It is the expansion factor. The expansion growth rate, To measure the residual vector, This is the residual judgment threshold vector. This is the residual threshold coefficient, used to amplify the boundary of the measurement residual. After each expansion, the prediction update of the current sampling step is re-executed until all observation components do not exceed the threshold or the maximum number of iterations is reached.

[0036] In this embodiment, the prior scale of the maneuver coefficient is taken as 1×10. -16 The expansion growth rate is set to 2, and the maximum number of compensation iterations is set to 3. Therefore, the covariance sub-block of the maneuver coefficient during compensation is set to 1×10. -16 2×10 -16 4×10 -16 The sequence is iteratively amplified.

[0037] S7. Based on the estimation results of the corrected maneuver coefficient, reconstruct the thrust acceleration norm and construct the output switching function, extract the final continuous low thrust maneuver period and output the thrust acceleration inversion results.

[0038] For each sampling interval, select within the standard interval Averaging the Chebyshev expansion results at discrete points along the three axes yields a piecewise axial acceleration estimate; further, the thrust acceleration norm is defined: ; ; Construct an output switching function to extract the maneuver period from the thrust acceleration inversion curve: ; in, Extract the proportional coefficient for the output time period. For the first Each sampling interval, along The first axial direction Estimates of the second-type Chebyshev expansion coefficients For the first Discrete points of a second-order Chebyshev polynomial in the standard interval The value at the position, where the subscript is in the formula Indicates the index of discrete points within a standard interval. For the first The thrust acceleration norm estimate for each sampling interval , , For the first The average maneuvering acceleration estimate of each component in each sampling interval The first of the thrust acceleration norm sequences in the entire sampling interval Item, subscript Indicates the sampling interval index; based on the output switching function The non-zero interval is used to obtain the final maneuvering time period, which is used for plotting the switching function in the attached figure and statistical analysis of the implementation results.

[0039] In this embodiment, the output time period extraction scaling factor is set to 0.35, and the peak thrust acceleration norm is 5.002 × 10⁻⁶. -6 km / s 2 The corresponding output switching threshold is 1.751 × 10⁻⁶. -6 km / s 2 The final output maneuver range is from 1020 s to 2940 s, with a starting delay of 120 s and an ending delay of 240 s.

[0040] like Figure 3 As shown, the peak value of the inverted thrust acceleration norm is 5.002 × 10⁻⁶. -6 km / s 2 , compared to the actual peak value of 5.000×10 -6 km / s 2 Basically the same; such as Figure 4 As shown, the maneuver period obtained by using the output switching function is from 1020 s to 2940 s, which is consistent with the main maneuver coverage relationship of the actual maneuver range of 900 s to 2700 s.

[0041] The above embodiments demonstrate that the present invention can complete the statistical detection, covariance expansion compensation, thrust acceleration inversion, and maneuver period output of continuous low-thrust maneuvers of non-cooperative spacecraft under the condition of only known observation data and coarse initial orbit values. The output results can be used for subsequent situation assessment and mission analysis.

[0042] Therefore, the present invention adopts the above-mentioned method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft. Under the condition that only observation data and rough initial orbit values ​​are known, it can realize online detection, thrust acceleration inversion and maneuver period output of continuous low-thrust maneuvers by means of Chebyshev polynomial extended state, no-maneuver reference threshold calibration and maneuver coefficient covariance expansion compensation mechanism.

[0043] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft, characterized in that, Includes the following steps: S1. Based on the rough initial orbit values ​​of the target spacecraft and the observation spacecraft, establish a spacecraft dynamics model that considers two-body gravity, second-order spherical harmonic perturbation and unknown continuous low thrust input, and obtain the basic state equations. S2. Based on the basic state equations obtained in step S1, the unknown continuous low thrust acceleration in each sampling interval is mapped to the standard interval, and the second kind of Chebyshev polynomial is used for piecewise expansion to obtain the parameterized expression of the maneuver coefficient. S3. Based on the parameterized expression of the maneuver coefficient in step S2, combine the maneuver coefficient vector with the orbital state to form a generalized state and establish an extended state transition relationship to obtain the extended system prediction model. S4. Based on the extended system prediction model, establish a multi-observer dual-angle, distance, and relative velocity hybrid observation model and perform extended Kalman filter prediction update to obtain the joint estimation results of the orbital state and maneuver coefficients of the current sampling step. S5. Based on the joint estimation results, the non-maneuvering reference trajectory is pre-run and the detection statistical baseline is calibrated. Then, the maneuvering statistics in the real-time observation branch are smoothed and thresholded to obtain the maneuvering detection results. S6. Based on the maneuver detection results and the current sampling step observation residual exceeding the threshold, perform iterative expansion compensation on the maneuver coefficient covariance sub-block and re-complete the prediction update for the current sampling step to obtain the corrected maneuver coefficient estimation result. S7. Based on the estimation results of the corrected maneuver coefficient, reconstruct the thrust acceleration norm and construct the output switching function, extract the final continuous low thrust maneuver period and output the thrust acceleration inversion results.

2. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 1, characterized in that, Step S1 specifically includes: In a geocentric inertial frame, the spacecraft's state vector is denoted as a six-dimensional vector composed of position and velocity. The unknown continuous low-thrust acceleration is introduced as an external input into the state equation, resulting in: ； in, and These represent the spacecraft's position vector and velocity vector, respectively. This represents the gravitational acceleration due to two bodies. This represents the second-order spherical harmonic perturbation acceleration. Assign a matrix to the input. The unknown continuous low-thrust acceleration is represented by the following form: Second-order spherical harmonic perturbation acceleration is expressed as follows: ； in, The gravitational constant of Earth, These are the second-order spherical harmonic coefficients. The radius of the Earth's equator. The distance from the spacecraft to the Earth's center. The component represents the position; the output is the basic dynamic model under continuous low thrust maneuvering conditions, for subsequent parameterization and filtering.

3. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 2, characterized in that, Step S2 is as follows: Regarding the first sampling intervals First, map physical time to a standard interval. ,get: ； The second kind of Chebyshev polynomial takes the following recursive form: ； In the orbital coordinate system, the unknown acceleration is decomposed along the radial, tangential, and normal directions, and expanded using Chebyshev basis functions: ； in, , and Let the radial, tangential, and normal unit basis vectors of the orbital coordinate system be represented respectively, and their definitions be: ； in, , and These represent the Chebyshev expansion coefficients along the three axes, respectively. To expand the order, the above steps transform the unknown continuous low thrust acceleration within each sampling interval into a finite-dimensional coefficient estimation problem.

4. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 3, characterized in that, In step S3, the maneuver coefficient vector and the orbital state jointly form the generalized state vector as follows: ； in, For the maneuver coefficient vector, , , Let Chebyshev expand the coefficient subvectors. , , These represent the directions in the orbital coordinate system; Since the Chebyshev coefficients within each sampling interval are treated as constant parameters, the extended system dynamics can be expressed as: ； Linearizing the extended system dynamics at the current estimation point yields the variational equation for the state transition matrix: ； in, Describing extended dynamics with respect to generalized states Jacobian matrix, Represents the identity matrix. Indicates in The state transition matrix of each sampling step is obtained by numerical integration. This is used for subsequent filtering and prediction.

5. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 1, characterized in that, Step S4 is as follows: For each observed spacecraft, the measurement vector is composed of observations of two angles, distance, and relative velocity. The overall observation model is represented as follows: ； The observation function for a single observer is taken as: ； Each component is determined by the relative position and relative velocity of the target spacecraft and the observation spacecraft: ； ； By performing first-order linearization on the observation function, we obtain the observation Jacobian matrix. And perform extended Kalman filter prediction updates: ； ； in, Indicates the first Each sampling step, From the state transition matrix Give, For process noise covariance, To measure the noise covariance, To measure the noise vector, The elevation angle of the target relative to the line of sight of the observing spacecraft. The azimuth angle of the target relative to the line of sight of the observation spacecraft. The relative distance between the target and the observed spacecraft. The magnitude of the relative velocity between the target and the observed spacecraft. , , These are the components of the relative position vector. Let V be the relative velocity vector of the target relative to the observing spacecraft. For the first The posterior state covariance matrix updated after each sampling step To change the state from the first Spread to the Nonlinear state transition function for each sampling step, These are the prior generalized state estimates. For the first The updated posterior generalized state estimate after each sampling step Let be the prior state covariance matrix. Here is the Kalman gain matrix. This is the updated posterior generalized state estimate. The predicted measurement vector is obtained by mapping the prior state estimate.

6. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 1, characterized in that, Step S5 obtains the non-maneuvering reference trajectory based on the coarse initial value integration. A pre-run is performed on the reference observations using a model consistent with the real-time filtering. The norm of the updated maneuvering coefficients at each step is defined as: ； Before pre-run Establish statistical baselines on each sample: ； In the real-time observation branch, the maneuver statistics are exponentially smoothed, and the presence of maneuver is determined accordingly. ； in, For smoothing coefficients, To detect the threshold coefficient, This is the vector of maneuver coefficients obtained through filtering estimation. The baseline of the mean of the mobility statistics for the sample of vehicles without pre-operational mobility is given. For the first phase of non-motorized pre-operation The mobility statistic corresponding to each sample The standard deviation of the motor statistics for the pre-operation sample without motors. This is the exponentially smoothed value of the motor statistic. The maneuver statistic is defined as the L2 norm of the maneuver coefficient estimation vector. For the first The exponential smoothing value of the motion statistic for each sampling step is used. When the smoothing statistic exceeds a threshold, the sampling interval is determined to be in a motion state. The output of this step is used both for subsequent compensation triggering and for constructing prior criteria for the output of the motion period.

7. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 5, characterized in that, In step S6, after the detection branch provides a maneuver decision, the compensation branch decides whether to initiate covariance expansion based on the observed residuals, defining the real-time observed residuals and component thresholds as follows: ； ； When there is a maneuver and at least one observation component satisfies When the value is greater than zero, it is considered that the maneuver coefficient covariance sub-block needs to be compensated. During the compensation iteration initialization, the orbital state covariance sub-block is not amplified, the maneuver coefficient covariance sub-block is set as the dilated prior scale, and the prediction update of the current sampling step is re-executed. ； in, Indicates the first The covariance sub-block of the mobility coefficients in the next compensation iteration For the prior measure of the mobility coefficient, It is the expansion factor. The expansion growth rate, To measure the residual vector, This is the residual judgment threshold vector. This is the residual threshold coefficient, used to amplify the boundary of the measurement residual. After each expansion, the prediction update of the current sampling step is re-executed until all observation components do not exceed the threshold or the maximum number of iterations is reached.

8. The method for detecting and inverting continuous low-thrust maneuvers of non-cooperative spacecraft according to claim 4, characterized in that, In step S7, for each sampling interval, select within the standard interval. Averaging the Chebyshev expansion results at discrete points along the three axes yields a piecewise axial acceleration estimate; further, the thrust acceleration norm is defined: ； ； Construct an output switching function to extract the maneuver period from the thrust acceleration inversion curve: ； in, Extract the proportional coefficient for the output time period. For the first Each sampling interval, along The first axial direction Estimates of the second-type Chebyshev expansion coefficients For the first Discrete points of a second-order Chebyshev polynomial in the standard interval The value at the position, where the subscript is in the formula. Indicates the index of discrete points within a standard interval. For the first The thrust acceleration norm estimate for each sampling interval , , For the first The average maneuvering acceleration estimate of each component in each sampling interval The first of the thrust acceleration norm sequences in the entire sampling interval Item; based on the output switching function The non-zero intervals are used to obtain the final maneuver time period.

Images