A process noise adaptive satellite autonomous continuous maneuver orbit determination method

Abstract

The application provides a satellite autonomous continuous maneuvering orbit determination method with adaptive process noise, and aims at the problem of orbit determination accuracy reduction of a satellite in a maneuvering process.The method is based on innovation, and adaptively adjusts process noise in a non-maneuvering stage and a maneuvering stage, so that the estimation accuracy reduction problem caused by fixed noise parameters is solved, and the filter can more accurately represent the dynamic state of the system.In a Gaussian noise, compared with a traditional unscented particle filter, the APNUPF exhibits better orbit determination accuracy during the maneuvering, and even if the state mutates, the APNUPF can still maintain high-precision orbit determination.In addition, in order to meet the actual engineering requirements, the robustness and reliability of the method for satellite autonomous continuous maneuvering orbit determination are comprehensively verified under non-Gaussian noise conditions.

Description

Technical Field

[0001] This invention relates to the field of spacecraft orbit determination and control technology, and in particular to an adaptive filtering method and system for high-precision orbit determination of on-orbit satellites during continuous maneuvers, belonging to the field of spacecraft simulation and control. Background Technology

[0002] With the surge in the amount of space debris in low Earth orbit and the continuous increase in the density of satellites in orbit, the frequency of orbital maneuvers performed by satellites to carry out missions such as collision avoidance, orbit transfer, and constellation maintenance has increased significantly. Real-time, high-precision orbit determination during satellite maneuvers is a core technology for ensuring maneuver control accuracy, achieving predetermined mission objectives, and guaranteeing space safety. Insufficient orbit determination accuracy during maneuvers may not only lead to mission failure but could also trigger serious space collision risks.

[0003] Currently, satellite orbit determination primarily relies on dynamic orbit determination methods, the core of which is to use Kalman filter algorithms (such as Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF)) to optimally estimate the satellite's state. However, the effectiveness of such algorithms heavily depends on an accurate system dynamic model and known noise statistics. During satellite maneuvers, the thrust generated by the thrusters introduces unknown accelerations, altering the system's dynamic characteristics. This leads to a severe "model mismatch" problem in filters based on fixed dynamic models. This mismatch causes increased state estimation bias, decreased filtering accuracy, and even filter divergence in strongly nonlinear maneuvering scenarios.

[0004] To address the model uncertainties introduced by maneuvers, numerous improvement schemes have been proposed. One approach attempts to accurately model or estimate the thrust of unknown maneuvers, such as using polynomial fitting, first-order Markov processes, or Fourier coefficients to model the thrust. However, the complexity of the space environment, the time-varying nature of the thruster's operating state, and the simplification inherent in the modeling process make it difficult to construct accurate thrust models using these methods, thus limiting their application in practical engineering. Another approach employs the Interactive Multiple Model (IMM) algorithm, which fuses the outputs of multiple models in a probabilistically weighted manner by running multiple filters corresponding to different maneuver assumptions in parallel. While the IMM algorithm improves robustness to some extent, its computational complexity increases dramatically with the number of models, making it unsuitable for satellite platforms with limited onboard computing resources.

[0005] To address nonlinear and non-Gaussian problems, the particle filter (PF) method, based on the Monte Carlo approach, has been introduced into the orbit determination domain. Particle filtering does not require the system noise to satisfy the Gaussian assumption and can approximate probability distributions of arbitrary shapes, theoretically solving the state estimation problem for strongly nonlinear systems. However, standard particle filtering suffers from severe particle degradation, where the weights of most particles approach zero, leading to a loss of particle diversity and reduced estimation accuracy. To address this issue, the unscented particle filter (UPF) algorithm is proposed. UPF utilizes the unscented transform (UT) to generate a better importance density function, effectively mitigating particle degradation and improving filtering accuracy.

[0006] Although the UPF performs well in nonlinear systems, a core challenge remains when applied to satellite maneuvering orbit determination: its process noise covariance matrix is ​​typically set to a fixed value. The uncertainty of the system model is time-varying at different stages of satellite flight. In the non-maneuvering phase, the system dynamics model is relatively accurate, and a smaller process noise setting helps achieve high-precision orbit determination. However, in the maneuvering phase, thrust introduces significant model uncertainty, requiring a larger process noise to ensure the robustness and tracking capability of the filter and prevent divergence due to model mismatch. If a smaller process noise setting is used in this phase, it will fail to effectively characterize the true dynamics of the system, resulting in the particle set not effectively covering the true state, which will also lead to deterioration in estimation accuracy or even divergence. Using a fixed process noise cannot simultaneously meet the conflicting requirements of these two phases; if the setting is too small, divergence will occur during maneuvers; if the setting is too large, the accuracy in the non-maneuvering phase will be sacrificed.

[0007] In summary, the core problem facing existing satellite maneuvering orbit determination technology lies in how to design a filtering algorithm that can adapt to dynamic changes in the system (especially the switching between maneuvering and non-maneuvering phases) and dynamically adjust the process noise level, thereby ensuring estimation accuracy under nonlinear and non-Gaussian noise environments while also considering the computational efficiency of the onboard platform. This is precisely the key technical problem that this invention aims to solve. Summary of the Invention

[0008] To address the issue of fixed noise covariance in dynamic scenarios in existing satellite maneuver orbit determination methods, this invention proposes a process noise adaptive satellite autonomous continuous maneuver orbit determination method. Based on the unscented particle filter framework, an innovation-based adaptive process noise covariance adjustment strategy is introduced, enabling the filter to sense real-time dynamic changes in the system and dynamically adjust model parameters, thereby achieving high-precision and robust orbit determination in both non-maneuver and continuous maneuver phases.

[0009] Technical solution:

[0010] A method for autonomous continuous maneuver orbit determination of a satellite with process noise adaptation includes the following steps:

[0011] Step 1: Establish the satellite orbital dynamics model, system state equations, and GNSS observation equations, and define the state vector containing position, velocity, and acceleration;

[0012] Step 2: Generate an initial particle set from the prior distribution and assign uniform weights to each particle to complete the initialization of the unscented particle filter;

[0013] Step 3: Importance sampling is performed based on Unscented Kalman Filter (UKF). Sigma points are generated based on the system state equation and GNSS observation equation described in Step 1, and particle states and state error covariance matrices are predicted. Then, UKF measurements are used to update and generate an optimized proposal distribution, from which a new particle set is sampled. The state error covariance matrix is ​​used to quantify prediction uncertainty and serves as the basis for weight update in Step 4.

[0014] Step 4: Update the new particle set and normalize the particle weights according to the importance sampling principle;

[0015] Step 5: Calculate the estimated satellite state at the current moment using the normalized weights;

[0016] Step 6: Calculate the predicted observation value based on the GNSS observation equation in Step 1 and the particle state predicted in Step 3. Use the difference between the actual observation value and the predicted observation value as the innovation. Calculate the weighted innovation norm after normalization and weighting. Adjust the process noise covariance matrix dynamically and adaptively. Feed the adjusted process noise covariance matrix back to Step 3 in the next iteration to update the noise parameters in the state prediction.

[0017] Step 7: Calculate the effective number of particles in the new particle set described in Step 3, resample when particles degenerate, and return the estimated result for the next iteration. Repeat steps 3 to 7.

[0018] Furthermore, in step 1, the satellite orbital dynamics model includes Earth's gravity, J2 term perturbation, Sun and Moon three-body perturbation, solar radiation pressure perturbation, and atmospheric drag perturbation, wherein the dynamic equations are expressed in the J2000 inertial coordinate system as follows:

[0019] ,

[0020] In the formula, The gravitational pull of Earth on the satellite, For J2 perturbation, and These are gravitational perturbations from the Sun and the Moon, respectively. Perturbed by solar radiation pressure It is atmospheric drag perturbation.

[0021] Furthermore, in step 1, the system's state equation and observation equation are expressed as follows: ,

[0022] In the formula, For state vectors, Let be the state vector from the previous time step. It is the observed vector / the actual measured value. and It is a nonlinear function. It is process noise, and its covariance matrix is ; It is measurement noise.

[0023] Furthermore, in step 2, the initialization includes: from the prior distribution N particles are generated, each particle having an initial state estimate. and the initial estimated covariance matrix Composition, and assigning initial weights to each particle. .

[0024] Furthermore, in step 3, the importance sampling includes:

[0025] Step 3.1, for each particle Generate Sigma points: In the formula, Estimate the mean of the state around the i-th particle at time k-1. A set of Sigma points were generated. Let be the state dimension. As a scale factor, Estimate the covariance matrix of the state of the i-th particle at time k-1;

[0026] Step 3.2, Predict particle state And calculate the predicted mean. With covariance :

[0027] , ,

[0028] In the formula, Let Sigma be the point of the i-th particle at time k-1. The j-th predicted Sigma point generated after nonlinear propagation, This represents the system control input at time k. and These are the mean weight and the covariance weight, respectively. It is the process noise covariance, where j and n are the indices of the Sigma points and the dimensions of the system state vector, respectively;

[0029] Step 3.3, predict the state Sigma point set Through the GNSS observation equation Nonlinear propagation is performed to obtain the Sigma point of the predicted observations. The predicted observation value corresponding to the i-th particle is calculated by weighting the observed Sigma points. ;

[0030] Step 3.4: Obtain the posterior mean using the UKF measurement update step. Covariance ;

[0031] Step 3.5, from the posterior mean Covariance Suggested distribution Mid-sampling generates new particles In the formula, Let be the state vector of the i-th particle at time k-1. This is the actual observation vector obtained at time k.

[0032] Furthermore, in step 4, the unnormalized weight of each particle is first calculated based on the importance sampling principle: Then, the weights of all particles are normalized to obtain the normalized weights. In the formula, Let be the weight of the i-th particle at time k-1. , Let be the state values ​​of the i-th particle at times k and k-1. The actual measured value at time k. , These are the probability density function and the proposal distribution function, respectively.

[0033] Furthermore, in step 5, the state estimate is the weighted average of the particle set: In the formula, To normalize the weights, Let be the state value of the i-th particle at time k.

[0034] Furthermore, in step 6, the adaptive adjustment process noise covariance matrix includes:

[0035] Step 6.1, Calculate the new information In the formula It is the information vector of the i-th particle at time k. It is the measurement value at time k. It is the predicted observation value of the i-th particle at time k;

[0036] Step 6.2, utilizing the forgetting factor And the length of the information window Calculate the weighted new interest:

[0037] In the formula, This represents the weighted information of the i-th particle at time k. The information is the historical value, and the superscript j indicates time j within the window;

[0038] Step 6.3, calculate the mean norm of the weighted innovation for all particles: , For normalized weights, N is the total number of particles;

[0039] Step 6.4, based on the mean norm of the weighted new information Adjustment process noise covariance :

[0040] ,

[0041] In the formula, and These are the up-adjustment and down-adjustment factors for the process noise covariance, respectively. and These are the upper and lower thresholds for the new information; and for The upper and lower limits.

[0042] Furthermore, in step 7, the resampling includes calculating the effective number of particles. ,like If the value is less than the threshold, resampling is performed, high-weight particles are copied, and the weights are reset to the uniform value.

[0043] Beneficial effects

[0044] (1) This invention addresses the model mismatch problem in satellite maneuvering scenarios, achieving rapid and accurate adaptive correction. Unlike traditional adaptive methods targeting measurement noise, this invention innovatively identifies the process noise covariance matrix as the core object of adaptive adjustment, directly addressing the abrupt changes in the system dynamics model caused by thruster operation during satellite maneuvers. By introducing a forgetting factor to calculate the weighted innovation norm and setting upper and lower thresholds for segmented adjustment, a step-like rapid response to sudden maneuvers is achieved. This allows the filter to rapidly increase process noise to match the model's uncertainty, preventing divergence caused by the accumulation of estimation bias. During non-maneuvering phases, it maintains relatively low process noise, thereby ensuring high-precision orbit determination. This targeted design overcomes the shortcomings of traditional methods, such as slow response or misjudgment during maneuvers.

[0045] (2) This invention combines the adaptive adjustment mechanism of novelty with unscented particle filtering, leveraging their synergistic advantages. This invention utilizes UKF to generate a proposal distribution for importance sampling for each particle. This distribution incorporates the latest observation information, significantly alleviating particle degradation problems and improving the accuracy of state estimation compared to standard particle filtering. This invention combines adaptive adjustment of process noise with particle resampling steps. While adapting to system dynamics by adjusting process noise, the resampling mechanism effectively maintains the diversity of the particle set, jointly ensuring the long-term stability of the filtering algorithm in highly nonlinear scenarios such as maneuvers. During the adaptive adjustment process, this invention sets upper and lower bound constraints on the process noise covariance matrix, effectively preventing over-adjustment due to individual abnormal observations or continuous maneuvers. Based on the particle filtering framework, this invention does not rely on the assumption that system noise satisfies a Gaussian distribution. Therefore, when facing non-Gaussian noise environments common in practical engineering, its orbit determination performance is significantly better than traditional algorithms based on Gaussian assumptions, such as Extended Kalman Filter (EKF).

[0046] (3) While ensuring excellent performance, this invention has low computational complexity, meeting the real-time requirements of resource-constrained environments such as spaceborne platforms. Compared with the Interactive Multiple Model (IMM) algorithm used to address model uncertainty, this invention does not require parallel operation of multiple filter models and probability-weighted fusion, significantly reducing computational load and single-step execution time. This makes it easier to implement real-time processing on spaceborne computers with limited computing power, memory, and power consumption, providing a feasible technical path for satellites to truly achieve autonomous continuous maneuvering orbit determination, and has significant engineering application value. Attached Figure Description

[0047] Figure 1 This is a target satellite orbit and GNSS data observation model according to an embodiment of the present invention;

[0048] Figure 2 This is a flowchart of an adaptive process noise unscented particle filter according to an embodiment of the present invention;

[0049] Figure 3 This is a diagram showing the satellite orbit position filtering estimation error of each algorithm under Gaussian noise according to an embodiment of the present invention;

[0050] Figure 4 This is a graph showing the variation of the innovation norm and the process noise covariance matrix under Gaussian noise according to an embodiment of the present invention.

[0051] Figure 5 This is a mixed Gaussian noise distribution diagram according to an embodiment of the present invention;

[0052] Figure 6 This is a quantile diagram according to an embodiment of the present invention;

[0053] Figure 7 This is a diagram illustrating the satellite orbital position filtering estimation error under mixed Gaussian noise according to an embodiment of the present invention.

[0054] Figure 8 This is a graph showing the variation of the innovation norm and process noise covariance matrix under mixed Gaussian noise according to an embodiment of the present invention. Detailed Implementation

[0055] This embodiment discloses a satellite autonomous continuous maneuver orbit determination method with process noise adaptation. To verify this method, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0056] Example 1

[0057] This embodiment aims to address the following problems existing in current satellite maneuvering orbit determination methods:

[0058] 1) PF and UPF use fixed process noise parameters, which cannot simultaneously meet the high accuracy requirements of the satellite in the non-maneuvering phase and the high robustness requirements in the maneuvering phase, resulting in model mismatch and filter performance degradation.

[0059] 2) Some existing adaptive algorithms (such as IMM) have excessively high computational complexity, making them difficult to implement on spaceborne computers.

[0060] 3) Traditional Kalman filters suffer from low orbit determination accuracy and easy divergence in non-Gaussian noise environments.

[0061] This embodiment discloses a satellite autonomous continuous maneuver orbit determination method that adapts to process noise, such as... Figures 1 to 8 As shown, it includes the following steps:

[0062] Step 1: Establish the satellite dynamics model, system state equations, and GNSS data observation equations to provide a mathematical framework for the subsequent filtering process.

[0063] In the J2000 inertial coordinate system, a high-precision orbital dynamics model of the satellite, incorporating gravitational, three-body, solar radiation pressure, and atmospheric drag perturbations, is established:

[0064] (1)

[0065] in, It is the Earth's gravity on the satellite. For J2 perturbation, and These are gravitational perturbations from the Sun and the Moon, respectively. Perturbed by solar radiation pressure It is atmospheric drag perturbation.

[0066] Solar gravitational perturbation and lunar gravitational perturbation In an inertial frame of reference, this is represented as:

[0067] (2)

[0068] in, and These are the position vectors of the Sun and the Moon, respectively. The satellite's position vector. and These represent the gravitational constants of the Sun and the Moon, respectively.

[0069] Solar pressure perturbation Simplified representation:

[0070] (3)

[0071] in, This is the light pressure coefficient. The ratio of the illuminated area to the satellite's mass. It is the solar radiation constant. This represents the average Earth-Sun distance.

[0072] Atmospheric drag perturbation Simplified representation:

[0073] (4)

[0074] In the formula, The ratio of windward area to satellite mass. This is the atmospheric drag coefficient. Atmospheric density, This is the satellite's velocity vector relative to the Earth's atmosphere.

[0075] The system's state equations and observation equations are as follows: This is the state vector.

[0076] (5)

[0077] in and It is a nonlinear function. It is process noise, and its covariance matrix is . It is measurement noise, and its covariance matrix is The core of this invention lies in solving the problem of model inaccuracy caused by abrupt changes in the system dynamics model, specifically by addressing the process noise covariance matrix. Make adaptive adjustments.

[0078] Step 2, filter initialization: generate an initial particle set and assign weights to provide initial state estimates for particle filtering: from the prior distribution. N initial particles are generated, each particle having an initial state estimate. and the initial estimated covariance matrix It is composed of particles, and each particle is assigned a uniform weight. Ensure that the filter starts operating from a reasonable prior distribution, and that the particle set of its output is directly input into step 3.

[0079] Step 3: Importance sampling is performed based on the Unscented Kalman Filter (UKF). Sigma points are generated based on the system state equation described in Step 1, and particle states are predicted. Then, UKF measurements are used to update and generate an optimized proposal distribution, from which a new particle set is sampled. Specifically, the steps include:

[0080] Step 3.1, for each particle Generate Sigma points:

[0081] (6)

[0082] In the formula, Estimate the mean of the state around the i-th particle at time k-1. A set of Sigma points were generated. Let be the state dimension. As a scale factor, Estimate the covariance matrix of the state of the i-th particle at time k-1.

[0083] Step 3.2, Predict particle state And calculate the predicted mean. With covariance :

[0084] (7)

[0085] (8)

[0086] In the formula, Let Sigma be the point of the i-th particle at time k-1. The j-th predicted Sigma point generated after nonlinear propagation, This represents the system control input at time k. and These are the mean weight and the covariance weight, respectively. It is the process noise covariance, where j and n are the indices of the Sigma points and the dimensions of the system state vector, respectively;

[0087] Step 3.3, Observation and Prediction: Observe the predicted state Sigma point set obtained in step b. GNSS observation equations through step 1 Nonlinear propagation is performed to obtain the Sigma point of the predicted observations. .

[0088] The predicted observation value corresponding to the i-th particle is calculated by weighting the observed Sigma points. , which serves as the basis for calculating new information.

[0089] Step 3.4: Obtain the posterior mean using the UKF measurement update step. Covariance .

[0090] Step 3.5, define the proposed distribution as a Gaussian distribution, derived from the posterior mean. Covariance Suggested distribution Mid-sampling generates new particles In the formula, Let be the state vector of the i-th particle at time k-1. The actual observation vector obtained at time k

[0091] Step 3 of this embodiment generates an importance sampling distribution through unscented Kalman filtering (UKF) to optimize the efficiency of particle filtering. Specifically, importance sampling is used in particle filtering to extract particles from the proposal distribution to approximate the true posterior distribution. Step 3 also utilizes UKF to generate a better proposal distribution because UKF handles nonlinear problems through unscented transformation, which is more accurate than the simple sampling of standard particle filtering, thereby alleviating the particle degeneration problem, i.e., the estimation bias caused by the concentration of particle weights. The output predicted state and covariance serve as the basis for innovation calculation.

[0092] Step 4: Update the new particle set and normalize the particle weights according to the importance sampling principle.

[0093] First, based on the importance sampling principle, update the unnormalized weights of each particle:

[0094] (9)

[0095] Then the weights of all particles are normalized:

[0096] (10)

[0097] Step 5: Calculate the estimated satellite state at the current moment using the normalized weights, which is the weighted average of all particles:

[0098] (11)

[0099] Step 6: Calculate the predicted observation value based on the GNSS observation equation in Step 1 and the particle state predicted in Step 3. Use the difference between the actual observation value and the predicted observation value as the innovation. Calculate the weighted innovation norm using the normalized weights. Based on this, dynamically and adaptively adjust the process noise covariance matrix. Feed the adjusted process noise covariance matrix back to Step 3 in the next iteration to update the noise parameters in the state prediction.

[0100] The innovation is used to correct the state estimate of the filter, thereby improving the accuracy of subsequent predictions. A large innovation indicates a decrease in prediction accuracy, suggesting a mismatch in the current system state model, i.e., entering a maneuvering state. Therefore, this invention uses innovation as an indicator for adaptively adjusting process noise.

[0101] Step 6.1, Calculate the new information :

[0102] (12)

[0103] in It is the information vector of the i-th particle at time k. It is the measurement value at time k. It is the predicted observation value of the i-th particle at time k.

[0104] Step 6.2 involves weighting the historical information to prevent over-adjustment of process noise due to abnormal measurements, thereby enhancing the system's robustness to anomalies. During the weighting process, earlier information is assigned lower weights, while newer information receives higher weights, ensuring that the adaptive adjustment of process noise is primarily based on recent changes in the system state. In this embodiment, a forgetting factor is used to weight the historical information:

[0105] (13)

[0106] in This represents the weighted information of the i-th particle at time k. The length of the information window. Forgetting factor, This is a new historical value.

[0107] Step 6.3, calculate the mean norm of the weighted innovation for all particles:

[0108] (14)

[0109] Step 6.4, based on the weighted new norm Adaptive adjustment of process noise covariance :

[0110] (15)

[0111] in, and These are the up-adjustment and down-adjustment factors for the process noise covariance, respectively. When an increase in innovation is detected, it is determined to be a model mismatch, therefore the process noise covariance is increased. This means telling the filter to now trust the observation data more. When the satellite enters the maneuvering phase, the process noise needs to increase rapidly to accommodate the system's sudden changes. The setting mechanism is to quickly increase the noise level; when the state returns to stability, the process noise is slowly reduced to dynamically maintain the current system stability. and These are the upper and lower thresholds for the new information. and Then it is The upper and lower limits of this parameter ensure that Within a reasonable range, avoid over-adjustment.

[0112] In this embodiment, a forgetting factor is introduced to perform a weighted average of the historical innovation sequence in the weighted innovation norm calculation. The forgetting factor ensures that the system pays more attention to recent dynamic changes and avoids interference from outliers. Based on the innovation threshold judgment, piecewise coefficients are used to directly intervene in process noise, achieving a step-like rapid response to sudden maneuvers. This solves the response speed problem during sudden maneuvers and significantly improves the response speed compared to traditional recursive estimation algorithms.

[0113] Step 7: Calculate the effective number of particles in the new particle set in Step 3, resample when particles degenerate, and return the state estimate for the next iteration. Repeat steps 3 to 7.

[0114] Step 7.1, calculate the effective number of particles:

[0115] (16)

[0116] Step 7.2, determine if resampling is needed. If Skip the following resampling steps; if ,right Resampling is performed, high-weight particles are copied, and low-weight particles are discarded to obtain a new particle set. And reset all particle weights Resampling can effectively alleviate particle degradation.

[0117] Output the current state estimate If resampling is performed, the resampled particle set is returned. and the weight of the reset If resampling is not performed, return the current particle set. and normalized weights Use the above results as the initial values ​​for the next iteration and return to step 3.

[0118] The following Examples 2 and 3 respectively demonstrate the performance under Gaussian noise environment and the robustness under non-Gaussian noise environment through experiments.

[0119] Example 2: Performance verification under Gaussian noise environment.

[0120] Set the following application background: Assume the total mass of the satellite is... The initial state is the semi-major axis. eccentricity Track inclination Right ascension of ascending node Perigeal argument True near point angle The specific maneuvering strategy is as follows: 1-50 seconds without inorganic power; 51-100 seconds of the first maneuver along the x-axis of the satellite orbital coordinate system; 101-150 seconds of the second consecutive maneuver along the y-axis of the satellite orbital coordinate system. The thrust for both maneuvers is 20N, and the fuel consumption rate is... .

[0121] Establish a satellite dynamics model; set up a state vector containing three axes of position, velocity, and acceleration; and as follows: Figure 1 As shown, a GNSS data observation equation containing three navigation satellites is constructed, and the measurement noise is set to Gaussian noise in this implementation.

[0122] like Figure 2 As shown, the unscented particle filtering process based on this adaptive process noise is simulated using 300 particles, and the parameters of the UPF algorithm are set. , , ; Set the parameters of the adaptive process noise strategy, a core step in this invention, including the innovation window length. Forgetting factor Up-adjustment factor of process noise covariance With the downward adjustment coefficient , The upper and lower limits are and ,in,

[0123] The upper threshold of new information With lower threshold .

[0124] Figure 3The position and velocity estimation errors of APNUPF, UPF, and IMM under Gaussian noise were compared. In the simulation, the UPF algorithm used a fixed process noise covariance matrix, APNUPF dynamically matched the system state by adjusting the process noise covariance matrix, and IMM dynamically adjusted the process noise covariance according to the model probability. The root mean square error (RMSE) of the positions of APNUPF, UPF, and IMM under the first maneuver was 1.12m, 1.87m, and 1.35m, respectively, and the RMSE of the positions under the second maneuver was 1.15m, 8.93m, and 1.54m, respectively. The single-step runtime of each particle in the three algorithms was 12.731ms, 21.620ms, and 40.625ms, respectively. In the state prediction stage of the UPF algorithm, i.e., in formula (8) The calculation shows that model mismatch leads to an underestimation of the actual state by the state covariance matrix. After multiple iterations, the error gradually accumulates, eventually causing the filter to diverge. APNUPF and IMM have better estimation accuracy. The RMSE of APNUPF for the two maneuvers was reduced by 40.11% and 87.12% compared to UPF, respectively.

[0125] Figure 4 The innovation norm and process noise covariance matrix The changes were observed. At the start of both maneuvers, the new information norm exceeded the upper threshold. At this point, adaptive adjustment is performed according to equation (15).

[0126] Example 3: Robustness verification under non-Gaussian noise environment.

[0127] The application background in this implementation is the same as that in Example 2.

[0128] Establish a satellite dynamics model; set up a state vector containing three axes of position, velocity, and acceleration; and as follows: Figure 1 As shown, a GNSS data observation equation containing three navigation satellites is constructed. In this implementation, the measurement noise is set to non-Gaussian noise.

[0129] During satellite maneuvers, the system may be subject to various noise interferences. For example, sensors during maneuvers may generate abnormal values ​​due to space radiation or electromagnetic interference, with amplitudes significantly exceeding the range of Gaussian noise. Gaussian mixture models (GMMs) simulate non-Gaussian noise by superimposing multiple Gaussian components, effectively describing the multimodal noise characteristics in complex noise environments. The distribution of mixed Gaussian noise is as follows: Figure 5 As shown, the quantile plot is as follows Figure 6As shown, this type of graph is used to determine whether the data follows a normal distribution. If the data perfectly conforms to a normal distribution, all points in the graph should fall on the diagonal (the black dashed line in the graph). However, the data at both ends of the graph deviate significantly from the normal reference line, exhibiting a shape with the right tail pointing upwards and the left tail pointing downwards, and the kurtosis is greater than the kurtosis of 3 for the normal distribution, indicating that the distribution is a heavy-tailed distribution; at the same time, the skewness of the three axes is not 0, indicating that the mixed Gaussian noise distribution used in this implementation is asymmetrical.

[0130] based on Figure 2 The adaptive process noise unscented particle filtering process shown has the same algorithm parameters and adaptive process noise strategy parameters as in Example 2.

[0131] like Figure 7 The figure shows the error curve for satellite orbital position filtering estimation under mixed Gaussian noise. Under mixed Gaussian noise conditions, the RMSE of APNUPF's position is 0.95m in the non-maneuvering phase, 1.61m in the first maneuvering phase, and 1.25m in the second maneuvering phase. After the maneuvering begins, the error curve shows some oscillations, but the error begins to converge after a period of time. Figure 8 The graphs show the variation curves of the innovation norm and the process noise covariance matrix under mixed Gaussian noise. This dynamic characteristic demonstrates that APNUPF can effectively overcome non-Gaussian noise interference and achieve stable estimation of orbital parameters under maneuvering conditions through particle weight adjustment, verifying the engineering applicability of the algorithm in satellite autonomous maneuvering orbit determination missions.

[0132] Examples demonstrate that this invention, through adaptive update of process noise, enables the filtering algorithm to better match the dynamic characteristics of the system. Regardless of whether the GNSS observation data exhibits Gaussian or non-Gaussian noise characteristics, it can achieve a position RMSE better than 2m during the satellite's autonomous maneuvering orbit determination. It can simultaneously meet the high precision requirements of the satellite's non-maneuvering phase and the high robustness requirements of the maneuvering phase. The computational complexity is approximately one-quarter of that of algorithms such as IMM, solving the problem of difficulty in implementing it on a spaceborne computer.

[0133] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A satellite autonomous continuous maneuver orbit determination method with process noise adaptive capability, characterized in that, Includes the following steps: Step 1: Establish the satellite orbital dynamics model, system state equations, and GNSS observation equations, and define the state vector containing position, velocity, and acceleration; Step 2: Generate an initial particle set from the prior distribution and assign uniform weights to each particle to complete the initialization of the unscented particle filter; Step 3: Importance sampling is performed based on unscented Kalman filter (UKF). Sigma points are generated and particle states are predicted based on the system state equation and GNSS observation equation described in Step 1. Then, UKF measurements are used to update and generate an optimized proposal distribution, and a new particle set is sampled from this distribution. Step 4: Update the new particle set and normalize the particle weights according to the importance sampling principle; Step 5: Calculate the estimated satellite state at the current moment using the normalized weights; Step 6: Calculate the predicted observation value based on the GNSS observation equation in Step 1 and the particle state predicted in Step 3. Use the difference between the actual observation value and the predicted observation value as the innovation. Calculate the weighted innovation norm using the normalized weights. Adjust the process noise covariance matrix dynamically and adaptively. Feed the adjusted process noise covariance matrix back to Step 3 in the next iteration to update the noise parameters in the state prediction. Step 7: Calculate the effective number of particles in the new particle set described in Step 3, resample when particles degenerate, and return the state estimate for the next iteration. Repeat steps 3 to 7.

2. The method according to claim 1, characterized in that, In step 1, the satellite orbital dynamics model includes Earth's gravity, J2 term perturbation, Sun and Moon three-body perturbation, solar radiation pressure perturbation, and atmospheric drag perturbation. The dynamic equations, expressed in the J2000 inertial coordinate system, are as follows: ， In the formula, The gravitational pull of Earth on the satellite, For J2 perturbation, and These are gravitational perturbations from the Sun and the Moon, respectively. Perturbed by solar radiation pressure It is atmospheric drag perturbation.

3. The method according to claim 1, characterized in that, In step 1, the system's state equation and observation equation are expressed as follows: ， In the formula, For state vectors, Let be the state vector from the previous time step. These are actual measured values. and It is a nonlinear function. It is process noise, and its covariance matrix is ; It is measurement noise.

4. The method according to claim 3, characterized in that, In step 2, the initialization includes: starting from the prior distribution N particles are generated, each particle having an initial state estimate. and the initial estimated covariance matrix Composition, and assigning initial weights to each particle. .

5. The method according to claim 1, characterized in that, In step 3, the importance sampling includes: Step 3.1, for each particle Generate Sigma points: ， In the formula, Estimate the mean of the state around the i-th particle at time k-1. A set of Sigma points were generated. Let be the state dimension. As a scale factor, Estimate the covariance matrix of the state of the i-th particle at time k-1; Step 3.2, Predict particle state And calculate the predicted mean. With covariance : ， ， In the formula, Let Sigma be the point of the i-th particle at time k-1. The j-th predicted Sigma point generated after nonlinear propagation, This represents the system control input at time k. and These are the mean weight and the covariance weight, respectively. It is the process noise covariance, where j and n are the indices of the Sigma points and the dimensions of the system state vector, respectively; Step 3.3, predict the state Sigma point set Through the GNSS observation equation Nonlinear propagation is performed to obtain the Sigma point of the predicted observations. The predicted observation value corresponding to the i-th particle is calculated by weighting the observed Sigma points. ; Step 3.4: Obtain the posterior mean using the UKF measurement update step. Covariance ; Step 3.5, from the posterior mean Covariance The proposed distribution Mid-sampling generates new particles In the formula, Let be the state vector of the i-th particle at time k-1. This is the actual observation vector obtained at time k.

6. The method according to claim 1, characterized in that, In step 4, the unnormalized weight of each particle is first calculated based on the importance sampling principle: Then, the weights of all particles are normalized to obtain the normalized weights. In the formula, Let be the weight of the i-th particle at time k-1. , Let be the state values ​​of the i-th particle at times k and k-1. The actual measured value at time k. , These are the probability density function and the proposal distribution function, respectively.

7. The method according to claim 1, characterized in that, In step 5, the state estimate is the weighted average of the particle set: ， In the formula, To normalize the weights, Let be the state value of the i-th particle at time k.

8. The method according to any one of claims 1-7, characterized in that, In step 6, the adaptive adjustment process noise covariance matrix includes: Step 6.1, Calculate the new information In the formula It is the information vector of the i-th particle at time k. It is the measurement value at time k. It is the predicted observation value of the i-th particle at time k; Step 6.2, utilizing the forgetting factor And the length of the new information window Calculate the weighted new interest: ， In the formula, This represents the weighted information of the i-th particle at time k. The information is the historical information value, where the superscript j represents the sampling sequence index within the sliding window, and k-W+j represents the corresponding historical moment within the window; Step 6.3, calculate the mean norm of the weighted innovation for all particles: , For normalized weights, N is the total number of particles; Step 6.4, based on the mean norm of the weighted new information Adjustment process noise covariance : ， In the formula, and These are the up-adjustment and down-adjustment factors for the process noise covariance, respectively. and These are the upper and lower thresholds for the new information; and for The upper and lower limits.

9. The method according to claim 1, characterized in that, In step 7, the resampling includes calculating the effective number of particles. ,like If the value is less than the threshold, resampling is performed, high-weight particles are copied, and the weights are reset to the uniform value.

Images