EM-EKF-based adaptive estimation method for lunar satellite formation orbit
By dynamically estimating the noise covariance matrix of lunar orbiting satellite formations using the EM-EKF algorithm, and combining it with a high-precision dynamic model and multi-source information fusion, the navigation accuracy and stability issues of lunar orbiting satellite formations under complex perturbation environments were solved, achieving high-precision autonomous navigation.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- INNOVATION ACAD FOR MICROSATELLITES OF CAS
- Filing Date
- 2026-03-13
- Publication Date
- 2026-07-10
AI Technical Summary
Existing lunar orbit satellite formation navigation systems struggle to achieve high-precision, autonomous orbit determination in complex perturbation environments. Traditional filtering methods rely on manually setting noise parameters, which cannot adapt to changes in inter-satellite distances during "breathing formations," leading to performance degradation and filter divergence.
An extended Kalman filter (EKF) algorithm based on expectation maximization (EM) is adopted to dynamically estimate the covariance matrix of process noise and observation noise through real-time observation data, thereby realizing online adaptive updating of filter parameters. This is combined with a high-precision lunar orbit dynamics model and multi-source information fusion measurement.
It improves the orbit determination accuracy and long-term stability of lunar orbit satellite formations, meets the reliability and verifiability requirements of space missions, reduces the dependence on training data and costs, and has good generalization ability.
Smart Images
Figure CN121855555B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of satellite technology, and in particular to an adaptive estimation method for lunar satellite formation orbits based on EM-EKF. Background Technology
[0002] Very low-wave radio observations are of great scientific significance in exploring the early evolution of the universe, studying high-energy celestial phenomena, and analyzing the structure of cosmic magnetic fields. However, due to the shielding of the Earth's ionosphere and the influence of ground-based electromagnetic noise, ground-based facilities struggle to achieve high-precision observations. Lunar orbit, with its electromagnetic environment far removed from Earth's interference and stable gravitational field, provides a unique natural platform for very low-wave astronomical observations. Deploying a "breathing formation" satellite system in lunar orbit, with the inter-satellite distances constantly changing according to mission requirements, exhibiting a "breathing" characteristic, can construct a space interferometric aperture array on the order of tens to hundreds of kilometers, enabling high-resolution detection of extremely low-frequency cosmic signals. To ensure the accuracy of formation observations and the smooth implementation of scientific missions, the lunar orbit breathing formation must possess highly autonomous and high-precision orbit determination capabilities, achieving formation maintenance and dynamic baseline control in complex perturbation environments. This places extremely high demands on the performance of the navigation system.
[0003] Currently, formation satellite navigation primarily employs nonlinear filtering methods such as Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF), and Particle Filter to determine orbits by fusing multi-source information including inter-satellite ranging, relative velocity measurement, and GNSS data. However, these methods generally rely on manually setting the covariance matrices of process noise and observation noise, lacking a systematic parameter determination mechanism. Under the complex perturbations and "breathing" configuration changes of lunar orbit, system dynamics and observation conditions continuously change, with inter-satellite distances spanning tens to hundreds of kilometers, causing significant fluctuations in measurement noise characteristics. Fixed parameters are difficult to adapt to different operational phases. Filtering parameters have a critical impact on accuracy, convergence, and stability. If they are mismatched with the actual system, performance degradation or even filter divergence may occur, thus limiting the long-term reliability and robustness of the formation navigation system.
[0004] In recent years, information fusion methods based on neural networks have attracted attention in the field of navigation. These methods, employing deep learning, RNNs, and LSTMs, learn system dynamics and noise characteristics from historical observations, which can be used for state estimation or to assist in adjusting parameters of traditional filters. These methods rely on the nonlinear modeling and feature extraction capabilities of neural networks, offering advantages in certain scenarios. However, their limitations are also significant: first, they heavily rely on large-scale, high-quality training data, while lunar orbital breathing formations are a novel mission, making it difficult to obtain sufficient real-world on-orbit samples; second, neural networks lack interpretability, failing to meet the reliability and verifiability requirements of space missions; and third, training costs are high, and their generalization ability is limited when facing unseen scenarios. These issues significantly restrict the application of neural network solutions in autonomous navigation for such missions. Summary of the Invention
[0005] This invention provides an adaptive estimation method for lunar satellite formation orbits based on EM-EKF. The proposed adaptive EKF navigation algorithm, which incorporates the expectation-maximization (EM) strategy, can dynamically estimate the covariance matrix of process noise and observation noise based on real-time observation data and filtering information during the filtering iteration process, thereby achieving online adaptive updating of filtering parameters.
[0006] This invention provides an adaptive estimation method for lunar satellite formation orbits based on EM-EKF, comprising:
[0007] The extended Kalman filter module of the navigation algorithm predicts the current orbital state of the lunar satellite formation based on the previous orbital state, the current observation data, and the process noise covariance matrix and observation noise covariance matrix; and
[0008] The expectation-maximization module of the navigation algorithm updates the process noise covariance matrix and the observation noise covariance matrix based on a state sequence consisting of a series of orbit states output by the extended Kalman filter module over a period of time.
[0009] Furthermore, it also includes:
[0010] A state model of the lunar satellite formation used for navigation algorithms is established in the lunar inertial coordinate system J2000; and
[0011] Establish a collaborative observation model for lunar satellite formations used in navigation algorithms.
[0012] Furthermore, establishing a state model for the lunar satellite formation used in the navigation algorithm includes:
[0013] The equation of motion for a single satellite is defined as:
[0014] ,
[0015] in, It is the satellite's three-axis inertial position. It is the satellite's three-axis inertial velocity; It is the lunar gravitational acceleration, which includes the effects of the lunar central gravity and non-spherical perturbations. The gravitational acceleration of the Earth, the Sun, and the Earth on that day;
[0016] Lunar gravitational acceleration Defined as:
[0017] ,
[0018] in, The gravitational constant of the moon, R That is the radius of the moon; n and m These are the order and degree, respectively. m = n When =0, it means that only the central gravity of the moon is considered; L It is the highest order of nonspherical perturbation. L =50; These are the normalized spherical harmonic coefficients. It is a normalized Legendre function. These are the satellite's fixed longitude and latitude;
[0019] Gravitational acceleration between the Sun and Earth Defined as:
[0020] ,
[0021] in, These are the gravitational constants of the Earth and the Sun, respectively. These are the position vectors of the Earth and the Sun relative to the Moon, respectively.
[0022] Based on the motion equations of a single satellite, a state model of the lunar satellite formation is established:
[0023] ,
[0024] in, These represent the primary star's three-axis inertial position and velocity, respectively. These represent the lunar gravitational acceleration and the Sun-Earth three-body acceleration of the host star, respectively. Let represent the three-axis inertial position and velocity of the i-th sub-star, respectively. denoted as the lunar gravitational acceleration and the Sun-Earth three-body gravitational acceleration of the i-th sub-star, respectively.
[0025] Furthermore, establishing a cooperative observation model for lunar satellite formations used in navigation algorithms includes:
[0026] The inter-satellite angle observation model is as follows:
[0027] ,
[0028] in, Representing the first i The azimuth and elevation angles between the minor star and the primary star; Indicates the three-axis inertial position of the primary star. Indicates the first i The three-axis inertial position of a star, This represents noise in inter-satellite angle measurements.
[0029] The inter-satellite distance observation model is as follows:
[0030] ,
[0031] in, Representing the i The distance between the minor star and the primary star. This indicates noise in inter-satellite distance measurements;
[0032] Based on inter-satellite angle and inter-satellite distance measurements, a collaborative observation model for satellite formations is established. Z :
[0033] ,
[0034] in, Representing the i Inter-satellite distance and angle measurement information between the minor star and the primary star. .
[0035] Furthermore, multiple satellites measure their line-of-sight angles with the main star using optical cameras and acquire distance information from the main star using X-band communication equipment. The line-of-sight angles include azimuth and elevation angles.
[0036] Furthermore, the prediction of the current orbital state of the lunar satellite formation by the extended Kalman filter module includes a prediction step and an update step, wherein:
[0037] The prediction steps include: based on the orbital state of the lunar satellite formation at the previous moment. The covariance matrix of the posterior state at the previous time step The process noise covariance matrix at the previous time step Calculate the prior estimate of the orbital state at the current moment. and the covariance matrix of the prior state at the current time The calculation formula is as follows:
[0038] ,
[0039] in, This is a nonlinear state transition function, i.e., the motion equation function of the lunar satellite formation. This is the state transition matrix;
[0040] The update steps include: prior estimation of the orbital state based on the current moment. and the prior state covariance matrix at the current moment Current observation data and the observation noise covariance matrix Perform calculations and output the posterior state estimate at the current time step. and posterior state covariance matrix The calculation formula is as follows:
[0041] ,
[0042] in, It is a nonlinear observation function, i.e., a cooperative observation model. To measure the Jacobian matrix, It is the gain matrix, the observed data. It includes the actual measured inter-star distances, azimuth angles, and elevation angles between each satellite and the main star, obtained through real-time measurements; This represents the observation residual.
[0043] Furthermore, the state sequence composed of the posterior state estimates output by the extended Kalman filter module over a period of time is regarded as a latent variable, and the process noise covariance matrix Q and the observation noise covariance matrix R are regarded as parameters to be estimated. The optimal parameter estimates are obtained by maximizing the log-likelihood function of the observed data.
[0044] Furthermore, the expectation-maximization module obtains the optimal parameter estimates through alternating iterations of the expectation step and the maximization step, including:
[0045] In the i In the next iteration, given the current parameter estimate... By calculating the log-likelihood of the complete data with respect to the posterior state distribution The expectation, construct the goal function:
[0046] ,
[0047] The complete data includes the observation matrix during the forward filtering process over a period of time, historical observation data, and the orbital state and state covariance matrix after reverse smoothing correction;
[0048] The desired outcome is approximated by the filtering result output by the extended Kalman filter module, including:
[0049] After substituting the state transition matrix and the Gaussian assumption of the cooperative observation model, the target Function decomposition into AND Q and R The two related items:
[0050] ,
[0051] The constant term const is irrelevant to the parameter to be estimated and is ignored.
[0052] The filtering result includes a state transition matrix stored for a period of time, a prior estimate of the orbital state and a prior state covariance matrix, a posterior state estimate and a posterior state covariance matrix, and an observation matrix;
[0053] In the desired step, based on the current parameter estimates... and Using observation data sequences The smoothed state estimate and smoothed covariance matrix are calculated through forward filtering and backward smoothing processes, where the superscript... i Indicates the first i iteration N The sliding time window length refers to the sequence of historical observation data; where:
[0054] The forward filtering employs a standard extended Kalman filter recursive process to obtain the posterior state estimate at each time step within the sliding time window. and posterior state covariance matrix ;
[0055] Backward smoothing works by recursively working backward from the last time step, using information from future time steps to correct historical state estimates, thus obtaining a smoothed state estimate. | and smooth covariance matrix and the cross-covariance matrix of states at adjacent time points. The recursive formula for the backward smoothing process is:
[0056] ,
[0057] in, J This is the smoothing gain matrix.
[0058] Furthermore, in the maximization step, the smoothed state estimate and smoothed covariance matrix obtained in the expectation step are used to update the estimate of the noise covariance matrix by maximizing the expectation of the log-likelihood function of the complete data.
[0059] For the target Function about Q and RBy taking the derivatives separately and setting them to zero, we obtain the closed-form update expression;
[0060] Smooth state estimate calculated based on the expected steps Define the state and predict the residual in one step:
[0061] ,
[0062] To minimize the second moment of the process error, for Maximize Q Updated formula:
[0063] ;
[0064] Smoothing covariance matrix calculated using the expected step ， Get the final Q Updated formula:
[0065] ,
[0066] in, Represents the cross-covariance between adjacent time points. F The state transition matrix is stored for the extended Kalman filter stage;
[0067] Smooth state estimate calculated based on the expected steps Observational innovation is defined as:
[0068] ;
[0069] right Maximize R Updated formula:
[0070] ;
[0071] Through a first-order Taylor expansion, we obtain the final... R Updated formula:
[0072] ,
[0073] in, H The measurement Jacobian matrix obtained by the extended Kalman filter module.
[0074] Furthermore, during the alternating iterations of the expectation step and the maximization step, a relative error criterion based on matrix norm is used to determine whether the parameter estimation has reached convergence.
[0075] After completing the first i The next iteration yields a new noise covariance matrix. and Then, calculate the relative error of parameter changes. :
[0076] ,
[0077] in, The Frobenius norm of a matrix is used to quantify the magnitude of changes in the matrix.
[0078] Set convergence threshold If satisfied If the iteration fails, it is considered to have converged, the iteration stops, and the current iteration is output. and As the optimal estimate;
[0079] Set the maximum number of iterations If the number of iterations Force the iteration to stop and output the current result.
[0080] The present invention has at least the following beneficial effects:
[0081] The lunar satellite formation orbit adaptive estimation method of the present invention integrates expectation maximization and extended Kalman filtering. During the iterative process of orbit state prediction by the extended Kalman filter module, the expectation maximization module can dynamically estimate the process noise covariance matrix and the observation noise covariance matrix based on the filtering results of the extended Kalman filter module, thereby realizing online adaptive updating of the extended Kalman filter parameters.
[0082] Updating the parameters of the extended Kalman filter using the expectation-maximization module does not require large-scale, high-quality training data, is low-cost, has strong generalization ability, and is interpretable, meeting the reliability and verifiability requirements of aerospace missions.
[0083] This invention constructs a high-precision dynamic model (state model) containing 50-order non-spherical gravitational perturbations, which can accurately describe the orbital evolution characteristics of formation satellites over long time scales, providing reliable, verifiable and engineering-usable navigation support for ultra-long-wave astronomical observation missions. Attached Figure Description
[0084] To further illustrate the above and other advantages and features of the various embodiments of the present invention, a more specific description of the embodiments of the invention will be presented with reference to the accompanying drawings. It is to be understood that these drawings depict only typical embodiments of the invention and are therefore not intended to limit its scope. In the drawings, identical or corresponding parts will be indicated by identical or similar reference numerals for clarity.
[0085] Figure 1 The flowchart of an adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to an embodiment of the present invention is shown.
[0086] Figure 2 The estimation results of the Kalman filter noise covariance parameter according to an embodiment of the present invention are shown.
[0087] Figure 3 The position error curve and velocity error curve of absolute navigation according to an embodiment of the present invention are shown. Detailed Implementation
[0088] It should be noted that the components in the accompanying drawings may be shown exaggerated for illustrative purposes and may not be to scale.
[0089] In this invention, the various embodiments are merely intended to illustrate the solutions of the invention and should not be construed as limiting.
[0090] In this invention, unless otherwise specified, the quantifiers “a” and “one” do not exclude scenarios involving multiple elements.
[0091] It should also be noted that, in the embodiments of the present invention, only a portion of the parts or components may be shown for clarity and simplicity. However, those skilled in the art will understand that, under the teachings of the present invention, the required parts or components can be added as needed for specific scenarios.
[0092] It should also be noted that within the scope of this invention, the terms "same", "equal", and "equal to" do not mean that the two values are absolutely equal, but allow for a certain reasonable error. In other words, the terms also cover "substantially the same", "substantially equal", and "substantially equal to".
[0093] It should also be noted that in the description of this invention, the terms "center," "longitudinal," "lateral," "upper," "lower," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the invention and for simplifying the description, and do not explicitly or implicitly suggest that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.
[0094] Furthermore, the embodiments of the present invention describe the process steps in a specific order. However, this is only for the convenience of distinguishing each step, and is not a limitation on the order of each step. In different embodiments of the present invention, the order of each step can be adjusted according to the process.
[0095] Addressing the urgent need for high-precision, autonomous navigation in lunar orbit breathing formations, this invention proposes an adaptive EKF navigation algorithm that integrates the expectation-maximization (EM) strategy. It also constructs a high-precision orbital dynamics model suitable for the long-term evolution of lunar orbits, achieving synergistic optimization of online estimation of navigation parameters (orbital state) and suppression of dynamic errors (errors caused by the dynamics model during on-orbit operation). This significantly improves the orbit determination accuracy and operational reliability of the formation in complex perturbation environments.
[0096] Figure 1 The flowchart of an adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to an embodiment of the present invention is shown.
[0097] like Figure 1 As shown, the adaptive estimation method for lunar satellite formation orbits based on EM-EKF includes the following steps:
[0098] Step 1: Establish a state model of the lunar satellite formation for the navigation algorithm in the lunar inertial coordinate system J2000.
[0099] A high-precision dynamic model (state model) of the satellite formation is established in the lunar inertial J2000 coordinate system, mainly considering the effects of 50th-order lunar nonspherical perturbations and Sun-Earth three-body perturbations. The motion equations of a single satellite are defined as follows:
[0100] ,
[0101] in, It is the satellite's three-axis inertial position. It is the satellite's three-axis inertial velocity; It is the lunar gravitational acceleration, which includes the effects of the lunar central gravity and non-spherical perturbations. The gravitational acceleration of the Earth's three bodies on that day.
[0102] Lunar gravitational acceleration Defined as:
[0103] ,
[0104] in, The gravitational constant of the moon, R That is the radius of the moon; n and m These are the order and degree, respectively. m = n When =0, it means that only the central gravity of the moon is considered; L It is the highest order of non-spherical perturbation, in this invention, L =50. These are the normalized spherical harmonic coefficients. It is a normalized Legendre function. These are the satellite's fixed longitude and latitude.
[0105] Gravitational acceleration between the Sun and Earth Defined as:
[0106] ,
[0107] in, These are the gravitational constants of the Earth and the Sun, respectively. These are the position vectors of the Earth and the Sun relative to the Moon, respectively.
[0108] Based on the above equations of motion for a single satellite, a state model for the lunar satellite formation is established:
[0109] ,
[0110] in, These represent the primary star's three-axis inertial position and three-axis inertial velocity, respectively. These represent the lunar gravitational acceleration and the Sun-Earth three-body acceleration of the host star, respectively. Let represent the three-axis inertial position and velocity of the i-th sub-star, respectively. denoted as the lunar gravitational acceleration and the Sun-Earth three-body gravitational acceleration of the i-th sub-star, respectively.
[0111] To address the complex perturbations in the lunar orbital environment, a high-precision orbital dynamics model was established, incorporating a 50th-order lunar non-spherical gravitational field and three-body gravitational perturbations. This model accurately characterizes the long-term orbital evolution of constellation satellites at baselines on a scale of hundreds of kilometers, providing a reliable physical basis for state prediction and significantly outperforming dynamics models that only consider central gravity or low-order perturbations.
[0112] Step 2: Establish a collaborative observation model for the lunar satellite formation used in the navigation algorithm.
[0113] In the formation measurement system, the nine satellites measure their line-of-sight angles (azimuth and elevation) with the main star using optical cameras, and acquire distance information with the main star using X-band communication equipment.
[0114] The inter-satellite angle observation model is as follows:
[0115] ,
[0116] in, Representing the first i The azimuth and elevation angles between the minor star and the primary star; Indicates the three-axis inertial position of the primary star. Indicates the first iThe three-axis inertial position of a star, This represents noise in inter-satellite angle measurements.
[0117] The inter-satellite distance observation model is as follows:
[0118] ,
[0119] in, Representing the i The distance between the minor star and the primary star. This indicates noise in inter-satellite distance measurements.
[0120] Based on the above inter-satellite angle and inter-satellite distance measurements, a collaborative observation model for satellite formations is established. Z :
[0121] ,
[0122] in, Representing the i Inter-satellite distance and angle measurement information between the minor star and the primary star. .
[0123] The collaborative observation model describes the functional relationship between observations and state variables, and is the basis for the filter to update the state estimate using observation data.
[0124] This invention employs a multi-source information fusion measurement scheme. Inter-satellite angle measurements are performed between the primary satellite and each secondary satellite using optical cameras to acquire the azimuth and elevation angles of the formation satellites. Simultaneously, inter-satellite distance measurements are performed using an X-band microwave link. By combining optical angle observations and microwave distance observations, an inter-satellite collaborative observation model is constructed, achieving a complete three-dimensional constraint on the relative positions of the formation.
[0125] Step 3: The extended Kalman filter module of the navigation algorithm predicts the orbital state of the lunar satellite formation at the current moment based on the orbital state of the lunar satellite formation at the previous moment, the observation data at the current moment, and the process noise covariance matrix and the observation noise covariance matrix; the expectation maximization module of the navigation algorithm updates the process noise covariance matrix and the observation noise covariance matrix based on the state sequence composed of a series of orbital states output by the extended Kalman filter module over a period of time.
[0126] Here, the navigation algorithm of the present invention refers to the EM-EKF (Expectation Maximization-Extended Kalman Filter) algorithm.
[0127] The navigation algorithm estimates the process noise covariance matrix online. Covariance matrix of observation noise This algorithm enables dynamic adjustment of the extended Kalman filter (EKF) parameters, significantly improving the accuracy, robustness, and long-term stability of trajectory determination. The overall algorithm framework consists of a collaborative EKF state estimation module and an EM parameter update module.
[0128] (1) EKF state estimation
[0129] Let the state vector of the formation satellites at time k be... This includes the three-dimensional inertial position and velocity of the formation satellites. Represents an n-dimensional real space; the observation vector is It is composed of the inter-star distances and inter-star angles of the formation. Let m represent the m-dimensional real space. The formation state model and observation model are respectively expressed as:
[0130] ,
[0131] in for k The state vector at any given moment contains the position and velocity information of each satellite in the formation; The nonlinear state transition function is used to perform trajectory recursion based on the established dynamic model (state model) containing 50th-order nonspherical perturbations; for k The observation vector at any given time includes the line-of-sight angle (azimuth and elevation) between the minor star and the primary star, as well as the inter-star distance measurement. It is a nonlinear observation function; and These are process noise and observation noise, assumed to be zero-mean Gaussian white noise. The covariance matrix is , The covariance matrix is , and It is usually unknown and changes with the stage of the task.
[0132] EKF linearizes the nonlinear model using a first-order Taylor expansion, recursively performing two steps: prediction and update. In other words, the extended Kalman filter module predicts the current orbital state of the lunar satellite formation, including both prediction and update steps.
[0133] The prediction step is based on the posterior estimate of the orbital state of the lunar satellite formation at the previous moment (i.e., the orbital state of the lunar satellite formation). The covariance matrix of the posterior state at the previous time step process noise covariance matrix Calculate the prior estimate of the orbital state at the current moment. and prior state covariance matrix The calculation formula is as follows:
[0134] ,
[0135] in, It is a nonlinear state transition function (the motion equation function of the lunar satellite formation). The input is the orbital state estimate of the lunar satellite formation at the previous time (time k-1), and the output is the prior estimate of the orbital state at the current time (time k). The state transition matrix is obtained by applying the nonlinear state transition function. It is obtained by performing local linearization.
[0136] The update steps are based on the prior estimate of the orbital state at the current moment. The prior state covariance matrix at the current moment Current observation data and observation noise covariance matrix Perform calculations and output the posterior state estimate of the orbit at the current time. and posterior state covariance matrix The calculation formula is as follows:
[0137] ,
[0138] in, This is a nonlinear observation function (cooperative observation model). To measure the Jacobian matrix, by adjusting the nonlinear observation function We obtain it by performing local linearization; It is the gain matrix, the observed data. It includes the actual measured inter-star distances, azimuth angles, and elevation angles between each satellite and the main star, obtained through real-time measurements; This represents the observation residual, which is the error between the actual measured observation data and the theoretical observation information calculated through a nonlinear observation function.
[0139] Traditional EKF requires manual pre-setting of the process noise covariance matrix and the observation noise covariance matrix, making it difficult to adapt to noise characteristic drift caused by large changes in inter-satellite spacing and single-machine switching operations during "breathing formation". To address this, an EM strategy is introduced to achieve online estimation of the covariance matrix.
[0140] In the EM-EKF algorithm framework, the state sequence composed of the posterior state estimates of the EKF output over a period of time is regarded as a latent variable, and the process noise covariance matrix Q and the observation noise covariance matrix R are regarded as parameters to be estimated. The optimal parameter estimates are obtained by maximizing the log-likelihood function of the observed data (inter-satellite distance and inter-satellite angle). The expectation maximization module gradually approximates the maximum likelihood estimate of the parameters (that is, the optimal parameter estimate) through alternating E-step (expectation step) and M-step (maximization step).
[0141] Because directly maximizing the log-likelihood function of the observed data requires integration over the latent variables (state sequence), it is difficult to solve analytically. The EM algorithm iterates through E-steps and M-steps, ensuring that the log-likelihood function value does not decrease in each step, thus gradually approaching its maximum point. The final convergence point is the maximum likelihood estimate of the parameters. In EM-EKF, the maximum likelihood estimate is the objective, and maximizing the log-likelihood function of the observed data is the means to achieve the objective.
[0142] In the i In the next iteration, given the current parameter estimate... By calculating the log-likelihood of the complete data with respect to the posterior state distribution The expectation, construct the goal function:
[0143] ,
[0144] The complete data includes the observation matrix during the forward filtering process over a period of time, historical observation data, and the orbital state and state covariance matrix after inverse smoothing correction. The observation matrix is composed of a nonlinear observation function. The partial derivative with respect to the state variable x is obtained.
[0145] Since the system is a Gaussian nonlinear dynamic system, and the EKF already provides approximations of the first and second moments of the state, the expected value can be approximated using the filtered output of the EKF. The first moment approximation refers to the prior state estimate output by the EKF prediction step and the posterior state estimate output by the update step. The second moment approximation refers to the prior state covariance matrix output by the prediction step and the posterior state covariance matrix output by the update step.
[0146] The filtered results output by EKF include a state transition matrix stored for a certain period of time, a prior estimate of the orbital state and a prior state covariance matrix, a posterior state estimate and a posterior state covariance matrix, and an observation matrix.
[0147] Specifically, after substituting the state transition matrix and the Gaussian assumption of the cooperative observation model, the target... The function can be decomposed into AND Q and R The two related items:
[0148] ,
[0149] The constant term const is independent of the parameter to be estimated and can be ignored.
[0150] In step E, based on the current parameter estimates... and (superscript) i Indicates the first i (second iteration), using the observed data sequence (in N The sliding window length is used to calculate the corrected state (smoothed state estimate) and state covariance matrix (smoothed covariance matrix) through forward filtering and backward smoothing processes. The observation data sequence refers to a sequence composed of historical observation data.
[0151] The sliding window length (or simply sliding window length) is a key parameter that determines how much historical data is used to update the process noise covariance matrix and the observation noise covariance matrix each time.
[0152] Forward filtering employs the standard EKF recursive process to obtain the posterior state estimate at each time step within the sliding time window. and posterior state covariance matrix .
[0153] Backward smoothing, on the other hand, recursively calculates backward from the last time step, using information from future time steps to correct historical state estimates, thus obtaining a smoothed state estimate. | and smoothed covariance matrix and the cross-covariance matrix of states at adjacent time points. .
[0154] The recursive formula for the backward smoothing process is:
[0155] ,
[0156] in, J This is the smoothed gain matrix. Smoothed state estimation has higher accuracy than EKF's posterior state estimation because it utilizes information from the entire observation time window, which is crucial for accurately estimating the noise covariance matrix.
[0157] In the M-step, using the corrected state (smoothed state estimate) and state covariance matrix (smoothed covariance matrix) obtained in the E-step, the estimate of the noise covariance matrix is updated by maximizing the expected value of the log-likelihood function of the complete data. For the target... Function about Q and RBy taking the derivatives separately and setting them to zero, we obtain the closed-form update expression.
[0158] Smooth state estimation based on E-step calculation Define the state and predict the residual in one step:
[0159] ,
[0160] To minimize the second moment of the process error, for Maximize Q Updated formula:
[0161] .
[0162] Combined with the smoothed covariance matrix calculated in the E-step To obtain the final Q Updated formula:
[0163] ,
[0164] in, Represents the cross-covariance between adjacent time points. F The state transition matrix is stored for the extended Kalman filter stage.
[0165] Smooth state estimation based on E-step calculation Observational innovation is defined as:
[0166] .
[0167] right Maximize R Updated formula:
[0168] .
[0169] Through a first-order Taylor expansion, we can obtain the final... R Updated formula:
[0170] ,
[0171] in, H This is the measurement Jacobian matrix obtained from the EKF module.
[0172] The EM algorithm is essentially an iterative optimization process. To determine whether the parameter estimation has reached convergence, this invention introduces a relative error criterion based on matrix norm. After completing the... i The next iteration yields a new noise covariance matrix. and Then, calculate the relative error of parameter changes. :
[0173] ,
[0174] in, This represents the Frobenius norm (F-norm) of the matrix, used to quantify the magnitude of matrix changes. The algorithm termination logic is as follows:
[0175] (a) Convergence determination: Set a convergence threshold If satisfied If the algorithm has converged, it stops iterating and outputs the current value. and As the optimal estimate.
[0176] (b) Maximum Iteration Protection: To prevent the algorithm from falling into an infinite loop due to failure to meet convergence accuracy under extreme observation conditions, a maximum number of iterations is set. If the number of iterations Force the iteration to stop and output the current result.
[0177] By employing the aforementioned convergence determination strategy, the algorithm can effectively control computational costs while ensuring the accuracy of parameter estimation, thus meeting the real-time and stability requirements of onboard computers.
[0178] An EM-EKF adaptive navigation architecture is adopted, embedding an Expectation-Maximization (EM) module on top of the traditional EKF framework. During algorithm execution, the EKF is first used to predict and update the orbital state. Then, based on the state residuals and observational innovations within the sliding time window, the expected value of the noise statistical characteristics is calculated through an E-step, followed by an M-step online update of the process noise covariance matrix and the observation noise covariance matrix. The updated parameters are fed back to the EKF in real time, forming a closed-loop adaptive mechanism. This design requires no manual intervention or prior noise information and can automatically adapt to the dynamic changes of the "breathing formation" at different baseline configuration stages, significantly improving navigation accuracy, convergence stability, and long-term robustness.
[0179] Figure 2 The estimation results of the Kalman filter noise covariance parameter according to an embodiment of the present invention are shown. Figure 3 The position error curve and velocity error curve of absolute navigation according to an embodiment of the present invention are shown.
[0180] Simulations were performed to verify the effectiveness of the EM-EKF-based adaptive estimation method for lunar satellite formation orbits of the present invention.
[0181] Since each satellite in the formation navigates via inter-satellite measurements with the primary satellite, the simulation employs a "one primary, one secondary" configuration (one primary satellite and one satellite) for absolute orbit determination. This configuration effectively represents the navigation scenario of the entire formation. In the navigation measurement system, the optical camera's angle measurement accuracy is 1 arcminute (3σ), and the X-band inter-satellite link's ranging accuracy is 3 meters (3σ). The initial orbital settings are as follows: position error of 3 kilometers (3σ) and velocity error of 3 meters per second (3σ).
[0182] Navigation simulation results show that the adaptive orbit estimation method of this invention has good parameter convergence characteristics. For example... Figure 2 As shown, the Kalman filter noise covariance parameter converges to a stable value after only 6 iterations. Based on this adaptive parameter, the navigation state estimation converges rapidly, as... Figure 2 As shown, the system completed the convergence process within 2 hours. Statistical analysis shows that the absolute position error after convergence is 884.56 meters (3σ), and the absolute velocity error is 0.8432 meters per second (3σ). These results verify that the proposed EM-EKF navigation method can achieve high-precision orbital state estimation, fully meeting the accuracy requirements for autonomous orbit determination in lunar orbit "breathing formations."
[0183] While some embodiments of the present invention have been described in this application, those skilled in the art will understand that these embodiments are merely illustrative. Numerous variations, alternatives, and improvements will arise in those skilled in the art under the teachings of this invention without departing from its scope. The appended claims are intended to define the scope of the invention and thereby cover methods and structures within the scope of the claims themselves and their equivalents.
Claims
1. An adaptive estimation method for lunar satellite formation orbits based on EM-EKF, characterized in that, include: The extended Kalman filter module of the navigation algorithm predicts the current orbital state of the lunar satellite formation based on the orbital state of the lunar satellite formation at the previous moment, the observation data at the current moment, and the process noise covariance matrix and the observation noise covariance matrix. as well as The expectation-maximization module of the navigation algorithm updates the process noise covariance matrix and the observation noise covariance matrix based on a state sequence composed of a series of orbital states output by the extended Kalman filter module over a period of time. This includes: treating the state sequence composed of posterior state estimates output by the extended Kalman filter module over a period of time as latent variables, treating the process noise covariance matrix Q and the observation noise covariance matrix R as parameters to be estimated, and obtaining the optimal parameter estimates by maximizing the log-likelihood function of the observed data; wherein: the expectation-maximization module obtains the optimal parameter estimates through alternating expectation and maximization steps; during the alternating iteration of expectation and maximization steps, the relative error criterion based on matrix norm is used to determine whether the parameter estimation has reached a convergence state; The expectation-maximization module obtains the optimal parameter estimates through alternating expectation and maximization steps, including: In the i In the next iteration, given the current parameter estimate... By calculating the log-likelihood of the complete data with respect to the posterior state distribution The expectation, construct the goal The function, where the complete data includes the observation matrix during the forward filtering process over a period of time, historical observation data, and the orbital state and state covariance matrix after reverse smoothing correction; The desired outcome is approximated by the filtering result output by the extended Kalman filter module, including: After substituting the state transition matrix and the Gaussian assumption of the cooperative observation model, the target Function decomposition into AND Q and R The two related items: The constant term const is irrelevant to the parameter to be estimated and is ignored. The filtering result includes a state transition matrix stored for a period of time, a prior estimate of the orbital state and a prior state covariance matrix, a posterior state estimate and a posterior state covariance matrix, and an observation matrix; In the desired step, based on the current parameter estimates... and Using observation data sequences The smoothed state estimate and smoothed covariance matrix are calculated through forward filtering and backward smoothing processes, where the superscript... i Indicates the first i iteration N The sliding time window length refers to the sequence of historical observation data; where: The forward filtering employs a standard extended Kalman filter recursive process to obtain the posterior state estimate at each time step within the sliding time window. and posterior state covariance matrix ; Backward smoothing works by recursively working backward from the last time step, using information from future time steps to correct historical state estimates, thus obtaining a smoothed state estimate. | and smoothed covariance matrix and the cross-covariance matrix of states at adjacent time points. ; In the maximization step, the smoothed state estimate and smoothed covariance matrix obtained in the expectation step are used to update the estimate of the noise covariance matrix by maximizing the expectation of the log-likelihood function of the complete data.
2. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF as described in claim 1, characterized in that, Also includes: A state model of the lunar satellite formation for navigation algorithms is established in the lunar inertial coordinate system J2000. as well as Establish a collaborative observation model for lunar satellite formations used in navigation algorithms.
3. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 2, characterized in that, Establishing a state model for the lunar satellite formation used in navigation algorithms includes: The equation of motion for a single satellite is defined as: ， in, It is the satellite's three-axis inertial position. It is the satellite's three-axis inertial velocity; It is the lunar gravitational acceleration, which includes the effects of the lunar central gravity and non-spherical perturbations. The gravitational acceleration of the Earth, the Sun, and the Earth on that day; Lunar gravitational acceleration Defined as: ， in, The gravitational constant of the moon, R That is the radius of the moon; n and m These are the order and degree, respectively. m = n When =0, it means that only the central gravity of the moon is considered; L It is the highest order of nonspherical perturbation. L =50; These are the normalized spherical harmonic coefficients. It is a normalized Legendre function. These are the satellite's fixed longitude and latitude; Gravitational acceleration between the Sun and Earth Defined as: ， in, These are the gravitational constants of the Earth and the Sun, respectively. These are the position vectors of the Earth and the Sun relative to the Moon, respectively. Based on the motion equations of a single satellite, a state model of the lunar satellite formation is established: ， in, These represent the primary star's three-axis inertial position and velocity, respectively. These represent the lunar gravitational acceleration and the Sun-Earth three-body acceleration of the host star, respectively. Let represent the three-axis inertial position and velocity of the i-th sub-star, respectively. denoted as the lunar gravitational acceleration and the Sun-Earth three-body gravitational acceleration of the i-th sub-star, respectively.
4. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 2, characterized in that, Establishing a cooperative observation model for lunar satellite formations used in navigation algorithms includes: The inter-satellite angle observation model is as follows: ， in, Representing the first i The azimuth and elevation angles between the minor star and the primary star; Indicates the three-axis inertial position of the primary star. Indicates the first i The three-axis inertial position of a star, This represents noise in inter-satellite angle measurements. The inter-satellite distance observation model is as follows: ， in, Representing the i The distance between the minor star and the primary star. This indicates noise in inter-satellite distance measurements; Based on inter-satellite angle and inter-satellite distance measurements, a collaborative observation model for satellite formations is established. Z : ， in, Representing the i Inter-satellite distance and angle measurement information between the minor star and the primary star. .
5. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 4, characterized in that, Multiple satellites measure their line-of-sight angles with the host satellite using optical cameras and acquire distance information from the host satellite using X-band communication equipment. The line-of-sight angles include azimuth and elevation angles.
6. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 1, characterized in that, The prediction of the current orbital state of the lunar satellite formation using the extended Kalman filter module includes a prediction step and an update step, wherein: The prediction steps include: based on the orbital state of the lunar satellite formation at the previous moment. The covariance matrix of the posterior state at the previous time step The process noise covariance matrix at the previous time step Calculate the prior estimate of the orbital state at the current moment. and the covariance matrix of the prior state at the current time The calculation formula is as follows: ， in, This is a nonlinear state transition function, i.e., the motion equation function of the lunar satellite formation. This is the state transition matrix; The update steps include: prior estimation of the orbital state based on the current moment. and the prior state covariance matrix at the current moment Current observation data and the observation noise covariance matrix Perform calculations and output the posterior state estimate at the current time step. and posterior state covariance matrix The calculation formula is as follows: ， in, It is a nonlinear observation function, i.e., a cooperative observation model. To measure the Jacobian matrix, It is the gain matrix, the observed data. It includes the actual measured inter-star distances, azimuth angles, and elevation angles between each satellite and the main star, obtained through real-time measurements; This represents the observation residual.
7. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 1, characterized in that, The recursive formula for the backward smoothing process is: ， in, J This is the smoothing gain matrix.
8. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 7, characterized in that, For the target Function about Q and R By taking the derivatives separately and setting them to zero, we obtain the closed-form update expression; Smooth state estimate calculated based on the expected steps Define the state and predict the residual in one step: ， To minimize the second moment of the process error, for Maximize Q Updated formula: ； Smoothing covariance matrix calculated using the expected step ， Get the final Q Updated formula: ， in, This represents the cross-covariance between adjacent time points. F The state transition matrix is stored for the extended Kalman filter stage; Smooth state estimate calculated based on the expected steps Observational innovation is defined as: ； right Maximize R Updated formula: ； Through a first-order Taylor expansion, we obtain the final... R Updated formula: ， in, H The measurement Jacobian matrix obtained by the extended Kalman filter module.
9. The adaptive estimation method for lunar satellite formation orbits based on EM-EKF according to claim 7, characterized in that, During the alternating iterations of the expectation step and the maximization step, a relative error criterion based on matrix norm is used to determine whether the parameter estimation has reached convergence: After completing the first i The next iteration yields a new noise covariance matrix. and Then, calculate the relative error of parameter changes. : ， in, The Frobenius norm of a matrix is used to quantify the magnitude of changes in the matrix. Set convergence threshold If satisfied If the iteration fails, it is considered to have converged, the iteration stops, and the current iteration is output. and As the optimal estimate; Set the maximum number of iterations If the number of iterations Force the iteration to stop and output the current result.