A micro-nano satellite attitude and orbit determination method based on adaptive weighted sliding window factor graph optimization

Abstract

This invention discloses a method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization, belonging to the field of micro / nano satellite autonomous navigation technology. It involves simultaneously acquiring panoramic infrared Earth images and MEMS gyroscope angular velocity data, extracting geocentric vector and visual feature observations; constructing a factor graph containing orbital dynamics factors, gyroscope pre-integration factors, zero-bias factors, geocentric vector factors, visual feature factors, and prior factors, and defining state variables; determining dynamic weight coefficients based on the residual consistency evaluation index of geocentric vector factors and visual feature factors; solving the nonlinear least squares objective function of each weighted factor to obtain the optimal state estimate; and iterating through a sliding window to achieve continuous integrated attitude and orbit determination. This invention expands the application scope of traditional Earth sensors, realizing autonomous orbit and attitude determination relying solely on panoramic infrared Earth sensors and MEMS gyroscopes, improving the accuracy and robustness of autonomous satellite orbit and attitude determination.

Description

Technical Field

[0001] This invention relates to the field of autonomous navigation technology for micro and nano satellites, and in particular to a method for determining the attitude and orbit of micro and nano satellites based on adaptive weighted sliding window factor graph optimization. Background Technology

[0002] Orbit and attitude information is a crucial foundation for satellites to successfully complete their space missions. Existing satellite attitude and orbit determination systems typically rely on multiple types of sensors or external navigation sources, resulting in complex system composition, high payload resource consumption, and high power consumption, which is detrimental to the miniaturization and lightweight design of micro and nano satellite platforms.

[0003] Panoramic infrared Earth sensors are primarily used for horizon detection and attitude measurement, extracting the Earth's edge to calculate the geocentric vector and providing attitude reference information for satellites. However, current applications typically only utilize their horizon observation capabilities, failing to fully exploit the surface information contained in their images. Meanwhile, MEMS gyroscopes can provide continuous angular velocity measurements, but they suffer from zero-bias drift, making it difficult to maintain high-precision estimations over long periods relying solely on inertial measurements. Therefore, it is necessary to fuse MEMS gyroscope data with other observational information to improve the accuracy and long-term stability of satellite state estimation.

[0004] Existing methods for fusing gyroscope and vision-based sensors mainly include filtering-based and optimization-based methods. Filtering-based methods typically estimate the system state recursively based on the current and previous states, limiting the utilization of historical observation information. Factor graph optimization-based methods can fuse multi-epoch, multi-source heterogeneous observation information within a unified framework, making fuller use of historical constraints and suitable for handling highly nonlinear attitude and trajectory estimation problems.

[0005] However, in complex space environments, panoramic infrared Earth images may be affected by factors such as cloud cover, texture degradation, feature mismatch, and edge extraction errors. If a fixed weight is applied to all observation factors, when local observations are abnormal or degraded, the abnormal observations may have an adverse effect on the overall optimization results, thereby reducing the accuracy and robustness of autonomous orbit and attitude determination.

[0006] Therefore, the purpose of this invention is to overcome the problems in the prior art, such as insufficient utilization of panoramic infrared Earth sensor information, long-term estimation accuracy decline caused by MEMS gyroscope zero-bias drift, and abnormal observations easily affecting optimization results under fixed observation weights. The invention proposes a micro-nano satellite autonomous orbit determination and attitude determination method based on adaptive weighted sliding window factor graph optimization. Summary of the Invention

[0007] The purpose of this invention is to provide a method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization. This method can achieve integrated joint estimation of satellite orbit and attitude using only the measurement information of panoramic infrared Earth sensor and MEMS gyroscope.

[0008] The technical solution adopted in this invention is as follows:

[0009] A method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization includes:

[0010] S1. Simultaneously acquire panoramic infrared Earth images collected by the panoramic infrared Earth sensor and angular velocity measurement values ​​collected by the MEMS gyroscope; extract geocentric vector measurement values ​​and visual feature observation values ​​from the panoramic infrared Earth images respectively;

[0011] S2. Define state variables including satellite position, velocity, attitude quaternions and gyroscope zero bias, and construct an initial factor graph within a sliding window; the factor graph includes orbital dynamics factors, gyroscope pre-integration factors, zero bias factors, geocentric vector factors, visual feature factors, and prior factors;

[0012] S3. Calculate the residual consistency evaluation index of geocentric vector factor and visual feature factor in the current window respectively, and determine the dynamic weight coefficient of geocentric vector factor and visual feature factor based on the evaluation index.

[0013] S4. Construct a nonlinear least squares objective function that includes each weighted factor, solve the objective function, obtain the optimal estimate of the state variables within the sliding window, and realize the integrated determination of satellite attitude and orbit.

[0014] S5. Advance the sliding window one epoch, return to steps S1 to S4, and continue until the task ends.

[0015] Furthermore, the prior factors in the factor graph are updated as the sliding window moves forward, and the update process is as follows:

[0016] As the sliding window slides forward, the historical state variables of the sliding window are marginalized, and the historical state variables are eliminated by the Schur complement operation to obtain the marginalized prior information matrix and prior information vector.

[0017] The prior information matrix is ​​decomposed by square root to obtain the Jacobian matrix and reference residual vector of the marginalization factor. The prior factor is then updated using the Jacobian matrix and reference residual vector.

[0018] Furthermore, the marginalization of the historical state variables of the slide-out window specifically includes:

[0019] Linearize all factors involving historical state variables in the current factor graph at the current iteration point, and construct a model containing the increments of state variables to be marginalized. and retain state variable increment The normal equation:

[0020] ;

[0021] in, , , and This represents the block division corresponding to the Hessian matrix. and This indicates that the information vector corresponds to a block;

[0022] Use the Schul complement operation to eliminate the relevant information in the normal equation. The part that obtains the marginalized prior information matrix and prior information vector :

[0023] .

[0024] Furthermore, updating the prior factors using the Jacobian matrix and the reference residual vector specifically includes:

[0025] Marginalized prior information matrix Perform square root decomposition to obtain the condition that satisfies Jacobian matrix and satisfy Reference residual vector ,in The prior information vector;

[0026] Constructing marginalized prior factors Its expression is:

[0027] ;

[0028] in, This represents the increment of the retained state variable relative to the linearization point;

[0029] The marginalized prior factors Replace the existing prior factors in the factor graph to complete the update.

[0030] Furthermore, the factors in the factor graph are specifically as follows:

[0031] Orbital dynamics factor The expression is:

[0032] ;

[0033] in, and Let be the state variables of the (k-1)th epoch and the kth epoch, respectively. Let be the orbital dynamics propagation function. For time intervals;

[0034] Gyroscope preintegration factor The expression is:

[0035] ;

[0036] in, The increment of the attitude matrix between the (k-1)th epoch and the kth epoch. The gyroscope has zero bias in the MEMS gyroscope at the (k-1)th epoch. This is the gyroscope bias error, which is the difference between the actual bias and the estimated bias of the gyroscope at the k-th epoch. These are the attitude matrices of the satellite relative to the inertial frame at the (k-1)th and kth epochs, respectively;

[0037] The zero-bias factor refers to the difference between the zero-bias estimate of the gyroscope in the current optimization iteration of the k-th epoch and the zero-bias estimate of the gyroscope after optimization convergence in the (k-1)-th epoch.

[0038] The calculation process of the geocentric vector factor includes: extracting Earth horizon pixels from panoramic infrared Earth images, and obtaining the geocentric vector measurement value in the machine system based on Earth ellipsoid model fitting; calculating the geocentric direction in the inertial frame based on the satellite position vector in the current iteration, and rotating the geocentric direction to the machine system using the attitude quaternion in the current iteration to obtain the geocentric vector prediction value; calculating the vector difference between the geocentric vector measurement value and the geocentric vector prediction value as the geocentric vector factor.

[0039] The calculation process of the visual feature factor includes: transforming the three-dimensional geographic coordinate points in the pre-stored reference map into the camera coordinate system based on the satellite position and attitude in the current iteration to obtain the three-dimensional spatial point coordinates; projecting the three-dimensional spatial point coordinates onto the two-dimensional image plane according to the intrinsic parameter model of the panoramic infrared earth sensor to obtain the predicted pixel coordinates of the feature points; calculating the coordinate difference between the predicted pixel coordinates and the measured pixel coordinates of the feature points actually matched in the image, and using the coordinate difference as the visual feature factor.

[0040] Furthermore, the geocentric vector measurement values ​​obtained based on the Earth ellipsoid model fitting specifically include:

[0041] Edge detection algorithms are used to extract the edge pixels of the Earth's horizon in panoramic infrared Earth images;

[0042] The edge pixels are mapped to the unit gaze vector using a camera inverse projection model;

[0043] Based on the quadratic constraint of the Earth ellipsoid, the least squares method is used to fit and solve all unit line-of-sight vectors to obtain the geocentric vector measurement value.

[0044] Furthermore, the residual consistency evaluation index in S3 is the Mahalanobis distance, and the calculation process includes:

[0045] Geocentric vector factor based on the k-th epoch or visual feature factors Calculate its corresponding covariance matrix. ;in, It is an observed factor variable. Time represents the covariance matrix corresponding to the geocentric vector factor. Time represents the covariance matrix corresponding to the visual feature factors;

[0046] Solve for the normalized Mahalanobis distance As an indicator for evaluating residual consistency; among which... hour, Represents the geocentric vector factor. hour, Represents visual feature factors.

[0047] Furthermore, the covariance matrix The expression is:

[0048] ;

[0049] in, This represents the observation matrix of the j-th class of observation factors in the k-th epoch. This represents the covariance matrix of the j-th type of factors estimated in the current state. Let represent the observation noise covariance matrix of the j-th class of observation factors in the k-th epoch.

[0050] Furthermore, the piecewise function is:

[0051] when At that time, weighting coefficient ;

[0052] when At that time, weighting coefficient Follow Increases while decreasing smoothly;

[0053] when At that time, weighting coefficient The corresponding factor will no longer participate in the current optimization;

[0054] in, and These are the lower threshold and the upper threshold, respectively.

[0055] Furthermore, the objective function expression in S4 is as follows:

[0056] ;

[0057] in, This represents the optimal estimate of the state variables within the sliding window of the k-th epoch. Indicates the first orbital dynamics factor of epoch The covariance matrix, Indicates the first Epochal gyro preintegral factor The covariance matrix, Indicates the first Zero bias factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. Geocentric vector factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. The i-th visual feature factor of the epoch The covariance matrix, Denotes the prior factor covariance matrix. Indicates the first The number of visual feature points involved in the optimization of each epoch.

[0058] The beneficial effects of this invention are:

[0059] (1) This invention fully utilizes the dual observation capabilities of the panoramic infrared Earth sensor. It not only uses the edge of the horizon to obtain geocentric vector constraints, but also uses the visual features of the Earth surface to obtain visual observation constraints, enabling the panoramic infrared Earth sensor to have both attitude observation and orbit observation capabilities. This breaks through the traditional application mode of Earth sensors that only provide single attitude reference information.

[0060] (2) The present invention adopts the factor graph optimization method to fuse orbital dynamics constraints, MEMS gyroscope measurement information and panoramic infrared Earth sensor observation information under a unified framework, which can make full use of the spatiotemporal correlation between multi-epoch and multi-source heterogeneous observations and improve the joint estimation accuracy of position, velocity and attitude.

[0061] (3) The present invention introduces a dynamic adaptive observation weighting mechanism based on residual consistency, which can identify, reduce weight or remove abnormal observations and degraded observations in real time according to the statistical characteristics of observation residuals, thereby effectively improving the robustness and stability of autonomous orbit and attitude determination algorithm in complex environments.

[0062] (4) The present invention introduces a sliding window edge-out and prior factor transfer mechanism, which effectively controls the scale of optimization variables while retaining historical information, and avoids the continuous increase of batch optimization computation over time, making it more suitable for spaceborne real-time processing scenarios.

[0063] (5) This invention can achieve satellite autonomous orbit determination and attitude determination by relying only on panoramic infrared earth sensor and MEMS gyroscope, which reduces the system’s dependence on global navigation satellite system and ground telemetry and control resources. It has good autonomy and engineering application value, and is especially suitable for resource-constrained micro and nano satellite platforms. Attached Figure Description

[0064] Figure 1 This is a flowchart illustrating a method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization.

[0065] Figure 2 This is a schematic diagram of a fusion estimation framework based on sliding window factor graph optimization;

[0066] Figure 3 This is a position estimation error diagram of the present invention and the extended Kalman filter algorithm;

[0067] Figure 4 This is a speed estimation error diagram of the present invention and the extended Kalman filter algorithm;

[0068] Figure 5 This is a diagram showing the attitude estimation error of the present invention and the extended Kalman filter algorithm. Detailed Implementation

[0069] The present invention will be further described and illustrated below with reference to specific embodiments. The embodiments described are merely examples of the content of this disclosure and do not limit the scope of the invention. The technical features of each embodiment in the present invention can be combined accordingly, provided that there is no mutual conflict.

[0070] This embodiment provides a method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization. This method utilizes a panoramic infrared Earth sensor to acquire panoramic infrared Earth images, extracts geocentric vector measurements and visual feature measurements, and combines these with angular velocity measurements output from a MEMS gyroscope to construct a sliding window factor graph. The optimized solution achieves joint estimation of satellite position, velocity, attitude, and gyroscope zero bias. The specific process is as follows: Figure 1 As shown, it includes the following steps:

[0071] S1, Data Acquisition

[0072] This invention simultaneously acquires raw data from two types of sensors:

[0073] (1.1) Acquire panoramic infrared Earth image I at the current epoch k using a panoramic infrared Earth sensor. k In this embodiment, the panoramic infrared Earth sensor adopts a panoramic isometric imaging model.

[0074] From panoramic infrared Earth image I k The geocentric vector measurement value and visual feature measurement value are extracted. The extraction process of the geocentric vector measurement value is as follows:

[0075] For the panoramic infrared Earth image I at the k-th epoch k First, the Sobel operator is used for edge detection to extract the pixels at the edge of the Earth's horizon. Let the coordinates of the extracted i-th edge pixel be... Based on the inverse projection model of the panoramic infrared camera, edge pixels are inversely mapped to unit line-of-sight vectors. .

[0076] In this embodiment, the Earth is described using an ellipsoidal model, and the quadratic matrix describing the geometry of the Earth ellipsoid is denoted as A. At the k-th epoch, the unit line-of-sight vector... Geocentric vector in the machine system Satisfy mathematical constraints:

[0077]

[0078]

[0079] in, The geocentric vector to be determined The function; extracts all horizon edge pixels at the k-th epoch. Substituting the above constraints, a least-squares optimization objective function for the geocentric vector can be constructed, and then the geocentric vector measurement value at this epoch can be obtained. .

[0080] The process of extracting visual feature measurements is as follows:

[0081] For image I k The SIFT (Simultaneous Scale Invariant Feature Transform) algorithm is used to extract visual feature points on the land surface, which are then matched with a pre-stored reference map containing geographic coordinate information. Let the pixel coordinates of the i-th matched visual feature point in the k-th epoch be . The corresponding known 3D geographic points in the reference map are After feature matching is completed, a set of visual feature observations consisting of image feature points and known geographic points can be obtained. .

[0082] (1.2) The high-frequency angular velocity measurement value of the current epoch k is acquired by the MEMS gyroscope. .

[0083] S2. Establish the geocentric vector theoretical model, the visual feature theoretical model, and the angular velocity observation model.

[0084] (2.1) Under the k-th epoch system, the geocentric vector theory model is:

[0085]

[0086] in, Let represent the attitude matrix of the satellite relative to the inertial frame at the k-th epoch. Let represent the attitude quaternion of the satellite relative to the inertial frame at the k-th epoch. This represents the coordinate vector of the satellite in the inertial frame at the k-th epoch.

[0087] (2.2) Based on the panoramic camera projection model, the ground truth values ​​of successfully matched pixel feature points. Corresponding 3D geographic points The following visual feature theoretical model is satisfied:

[0088]

[0089] Where θ is the incident angle and f is the focal length of the panoramic infrared camera. The principal point coordinates of the camera. These are the spatial coordinates of a 3D geographic point in the camera coordinate system.

[0090] (2.3) Angular velocity observation model:

[0091] Establish an angular velocity observation model:

[0092]

[0093] in, This represents the angular velocity measurement value of the MEMS gyroscope at the k-th epoch. This is the actual angular velocity value. The gyroscope bias of the MEMS gyroscope is zero at the k-th epoch. It is white noise.

[0094] S3. Construct a dynamic window factor graph

[0095] (3.1) Define the state variable of the k-th epoch as ,in This represents the coordinates of the satellite's position at epoch k. This represents the velocity vector of the satellite in the k-th epoch. This represents the attitude quaternion of the satellite at epoch k. This indicates that the satellite's gyroscope has zero bias at the k-th epoch.

[0096] Let the sliding window length be s, then the set of state variables to be optimized within the window at time k is represented as follows: .

[0097] (3.2) Constructing orbital dynamic factors based on orbital dynamics model

[0098] In this embodiment, the orbital dynamics model adopts the perturbed two-body orbital dynamics model, whose continuous form is as follows:

[0099]

[0100] in, The gravitational constant of Earth, The total perturbation acceleration includes non-spherical gravitational acceleration, atmospheric drag, solar radiation pressure, and third-body gravitational perturbation.

[0101] Discretize the above continuous form using the orbital dynamics propagation function. The discretized orbital dynamics model can be obtained:

[0102]

[0103] in, This represents the modeling and discretization error.

[0104] Based on the discretized orbital dynamics model, the orbital dynamics factor is defined as:

[0105]

[0106] In practical implementation, the orbital dynamics propagation function This can be achieved using the Runge-Kutta numerical integration method.

[0107] (3.3) Constructing the gyroscope pre-integration factor based on MEMS gyroscope angular velocity measurements

[0108] Based on the attitude kinematics model and MEMS gyroscope angular velocity measurements, the attitude matrix increment from epoch k-1 to epoch k can be obtained. :

[0109]

[0110] in, For the time interval from epoch k-1 to epoch k, Exp(·) is an exponential mapping function.

[0111] The gyroscope preintegration factor is defined as:

[0112]

[0113] Log(·) is a logarithmic mapping function. This is the gyroscope's zero-bias prediction value. It is the gyroscope zero-bias estimation error.

[0114] (3.4) Constructing a zero-bias factor based on a gyroscope zero-bias random walk model

[0115] The zero bias of a MEMS gyroscope is described using a random walk model, and its state evolution relationship is as follows:

[0116]

[0117] in, This represents the zero-biased random walk noise between the k-th epoch and the (k+1)-th epoch.

[0118] The zero bias factor is defined as:

[0119]

[0120] (3.5) Constructing geocentric vector factors based on geocentric vector measurements

[0121] Based on the geocentric vector theoretical model, the predicted geocentric vector value is obtained. ;

[0122] The geocentric vector factor is defined as:

[0123]

[0124] (3.6) Constructing visual feature factors based on visual feature observation

[0125] Based on the aforementioned visual feature theory model, the predicted values ​​of visual feature points are obtained. ;

[0126] Visual feature factors are defined as:

[0127]

[0128] (3.7) Constructing prior factors based on sliding window and marginalization

[0129] As the window slides forward, the oldest state variable is... As a state to be marginalized, it is removed from the set of state variables to be optimized, but its information is retained as a state variable. Prior constraints.

[0130] In this embodiment, an optional implementation process is as follows:

[0131] Let the increment of the state to be edged out in the current window be... The state increment is retained. Then the linearized normal equation can be written as:

[0132]

[0133] in, , , and This represents the block division corresponding to the Hessian matrix. and This indicates that the information vector corresponds to a block.

[0134] After eliminating the states to be marginalized using Schur complement, the prior information matrix after marginalization can be obtained. With prior information vector for:

[0135]

[0136] Furthermore, regarding the prior information matrix Perform square root decomposition to obtain the Jacobian matrix of the marginalization factor. and reference residual vector ,satisfy:

[0137]

[0138] The marginalization prior factor is defined as:

[0139]

[0140] in, This represents the increment of the retained state variable relative to the linearization point.

[0141] S4, Adaptive Observation Weighting

[0142] In this embodiment, dynamic weights are introduced for the two observation factors, geocentric vector factor and visual feature factor, and a dynamic adaptive weighting mechanism based on residual consistency is adopted.

[0143] Observation factors for the k-th epoch Construct its linearized residual covariance matrix :

[0144]

[0145] in, It is an observed factor variable. hour, Represents the geocentric vector factor. hour, Represents visual feature factors, This represents the observation matrix of the j-th class of observation factors in the k-th epoch. This represents the covariance matrix of the j-th type of factors estimated in the current state. Let represent the observation noise covariance matrix of the j-th class of observation factors in the k-th epoch.

[0146] Based on the linearized residual covariance matrix, a normalized Mahalanobis distance is constructed as a residual consistency evaluation index:

[0147]

[0148] For different types of observation factors, set piecewise dynamic weight functions:

[0149]

[0150] in and Let represent the lower and upper threshold values ​​of the adaptive error function for the sensor corresponding to the j-th type of observation factor, respectively. This is an adjustment factor that can be set according to the accuracy and reliability of different sensors. When When the observation is consistent with the current state prediction, the original weights are retained; when... When this occurs, it indicates that the observation has some degradation, and it should be downweighted; when If the observation is deemed a severe anomaly, the corresponding factor is removed and will no longer participate in the current window optimization.

[0151] Weights are applied by adjusting the factor covariance. The covariance matrix of the j-th class observation factor at the k-th epoch is dynamically weighted and then used as the covariance matrix. It can be represented as:

[0152]

[0153] in, This represents the covariance matrix of observation factors of type j before dynamic weighting.

[0154] S5, Solve

[0155] (5.1) Constructing the objective function

[0156] The constructed objective function includes orbital dynamics residuals, gyro pre-integration residuals, zero-bias residuals, geocentric vector residuals, visual feature reprojection residuals, and marginalized prior residuals. All weighted factors within the window are combined into a nonlinear least squares problem, specifically expressed as:

[0157]

[0158] in, This represents the optimal estimate of the state variables within the sliding window of the k-th epoch. Indicates the first orbital dynamics factor of epoch The covariance matrix, Indicates the first Epochal gyro preintegral factor The covariance matrix, Indicates the first Zero bias factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. Geocentric vector factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. The i-th visual feature factor of the epoch The covariance matrix, Denotes the prior factor covariance matrix. Indicates the first The number of visual feature points involved in the optimization of each epoch.

[0159] The present invention provides a fusion estimation framework based on sliding window factor graph optimization, as follows: Figure 2 As shown, circular nodes represent satellite state variables over time. The system employs a sliding window and edge-based control of computational scale to ensure real-time operation onboard. Constraints are established between state nodes through five types of factors: orbital dynamics factors, gyroscope pre-integration factors, and zero-bias factors provide intrinsic motion state predictions; visual feature factors and geocentric vector factors provide external absolute observation corrections.

[0160] (5.2) Iterative solution

[0161] In this embodiment, the Levenburg-Marquardt algorithm is preferably used to solve the objective function. When the preset convergence condition is met, the state estimation result of the current window's last epoch is output, including the satellite position, velocity, attitude, and gyroscope zero bias.

[0162] Subsequently, the sliding window advances one epoch and repeats the above steps to achieve integrated autonomous orbit and attitude determination estimation in continuous time series.

[0163] This embodiment proposes a method for orbit and attitude determination using a panoramic infrared Earth sensor and a MEMS gyroscope based on factor graph optimization. In a specific embodiment, the method is verified using simulation. The initial simulation time is set to 00:58:31 on January 4, 2023. The true orbit value is generated by STK software, and the true attitude value is derived from the attitude kinematics model. It is assumed that the simulated satellite maintains a ground-pointing attitude throughout the mission. The simulation parameters are set as follows: the gyroscope has a fixed bias of [value missing]. The random walk of the angle is 0.03° / h. 0.5 The bias stability is 0.15° / h, and the rate random walk is 0.3° / h. 1.5 The panoramic infrared earth sensor has a geocentric vector measurement accuracy of 0.2°, a feature extraction error standard deviation of 1 pixel, a cloud occlusion probability of 50%, and a sampling frequency of 1 Hz for both the panoramic infrared earth sensor and the MEMS gyroscope.

[0164] Figure 3 , Figure 4 and Figure 5 The comparison results between the sliding window factor graph optimization method and the extended Kalman filter method described in this invention regarding position, velocity, and attitude estimation errors are presented. The results show that, compared to the extended Kalman filter method, the method described in this invention can effectively improve the estimation accuracy of position, velocity, and attitude. The main reason for this is that this invention can more fully utilize the multi-epoch correlation information between MEMS gyroscope measurements, geocentric vector observations, and visual feature observations within a unified optimization framework, thereby improving the overall accuracy of state estimation. From the root mean square error results, the position, velocity, and attitude estimation errors of the extended Kalman filter method are 126.8725 m, 2.0395 m / s, and 0.0207°, respectively; while the position, velocity, and attitude estimation errors of the sliding window factor graph optimization method proposed in this invention are 94.3923 m, 1.8026 m / s, and 0.0199°, respectively. The simulation results demonstrate that the autonomous orbit and attitude determination method proposed in this invention, which integrates a panoramic infrared Earth sensor and a MEMS gyroscope, can effectively improve the joint estimation accuracy of satellite position, velocity, and attitude, thus verifying the effectiveness and feasibility of the technical solution of this invention.

[0165] The above are merely preferred embodiments of the present invention. It should be noted that, without departing from the core idea of ​​the present invention, the visual feature extraction algorithm, camera imaging model, dynamics propagation method, optimization solution algorithm, and edge detection implementation can all be replaced or adjusted accordingly. For example, the visual feature extraction algorithm can be replaced with other locally invariant feature extraction algorithms, and the nonlinear optimization algorithm can be replaced with the Gauss-Newton algorithm or other algorithms suitable for solving nonlinear least squares problems. All such replacements and adjustments should be considered to fall within the protection scope of the present invention.

Claims

1. A method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization, characterized in that, include: S1. Simultaneously acquire panoramic infrared Earth images collected by the panoramic infrared Earth sensor and angular velocity measurements collected by the MEMS gyroscope. Geocentric vector measurements and visual feature observations are extracted from the panoramic infrared Earth image, respectively. S2. Define state variables including satellite position, velocity, attitude quaternions and gyroscope zero bias, and construct an initial factor graph within a sliding window; the factor graph includes orbital dynamics factors, gyroscope pre-integration factors, zero bias factors, geocentric vector factors, visual feature factors, and prior factors; S3. Calculate the residual consistency evaluation index of geocentric vector factor and visual feature factor in the current window respectively, and determine the dynamic weight coefficient of geocentric vector factor and visual feature factor based on the evaluation index. S4. Construct a nonlinear least squares objective function that includes each weighted factor, solve the objective function, obtain the optimal estimate of the state variables within the sliding window, and realize the integrated determination of satellite attitude and orbit. S5. Advance the sliding window one epoch, return to steps S1 to S4, and continue until the task ends.

2. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 1, characterized in that, The prior factors in the factor graph are updated as the sliding window moves forward. The update process is as follows: As the sliding window slides forward, the historical state variables of the sliding window are marginalized, and the historical state variables are eliminated by the Schur complement operation to obtain the marginalized prior information matrix and prior information vector. The prior information matrix is ​​decomposed by square root to obtain the Jacobian matrix and reference residual vector of the marginalization factor. The prior factor is then updated using the Jacobian matrix and reference residual vector.

3. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 2, characterized in that, The marginalization of the historical state variables of the slide-out window specifically includes: Linearize all factors involving historical state variables in the current factor graph at the current iteration point, and construct a model containing the increments of state variables to be marginalized. and retain state variable increment The normal equation: ； in, , , and This represents the block division corresponding to the Hessian matrix. and This indicates that the information vector corresponds to a block; Use the Schul complement operation to eliminate the relevant information in the normal equation. The part that obtains the marginalized prior information matrix and prior information vector : 。 4. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 2 or 3, characterized in that, The step of updating the prior factors using the Jacobian matrix and the reference residual vector specifically includes: Marginalized prior information matrix Perform square root decomposition to obtain the condition that satisfies Jacobian matrix and satisfy Reference residual vector ,in The prior information vector; Constructing marginalized prior factors Its expression is: ； in, This represents the increment of the retained state variable relative to the linearization point; The marginalized prior factors Replace the existing prior factors in the factor graph to complete the update.

5. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 1, characterized in that, The factors in the factor graph are specifically as follows: Orbital dynamics factor The expression is: ； in, and Let be the state variables of the (k-1)th epoch and the kth epoch, respectively. Let be the orbital dynamics propagation function. For time intervals; Gyroscope preintegration factor The expression is: ； in, The increment of the attitude matrix between the (k-1)th epoch and the kth epoch. The gyroscope has zero bias in the MEMS gyroscope at the (k-1)th epoch. This is the gyroscope bias error, which is the difference between the actual bias and the estimated bias of the gyroscope at the k-th epoch. These are the attitude matrices of the satellite relative to the inertial frame at the (k-1)th and kth epochs, respectively; The zero-bias factor refers to the difference between the zero-bias estimate of the gyroscope in the current optimization iteration of the k-th epoch and the zero-bias estimate of the gyroscope after optimization convergence in the (k-1)-th epoch. The calculation process of the geocentric vector factor includes: extracting Earth horizon pixels from panoramic infrared Earth images, and obtaining the geocentric vector measurement value in the machine system based on Earth ellipsoid model fitting; calculating the geocentric direction in the inertial frame based on the satellite position vector in the current iteration, and rotating the geocentric direction to the machine system using the attitude quaternion in the current iteration to obtain the geocentric vector prediction value; calculating the vector difference between the geocentric vector measurement value and the geocentric vector prediction value as the geocentric vector factor. The calculation process of the visual feature factor includes: transforming the three-dimensional geographic coordinate points in the pre-stored reference map into the camera coordinate system based on the satellite position and attitude in the current iteration to obtain the three-dimensional spatial point coordinates; projecting the three-dimensional spatial point coordinates onto the two-dimensional image plane according to the intrinsic parameter model of the panoramic infrared earth sensor to obtain the predicted pixel coordinates of the feature points; calculating the coordinate difference between the predicted pixel coordinates and the measured pixel coordinates of the feature points actually matched in the image, and using the coordinate difference as the visual feature factor.

6. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 5, characterized in that, The geocentric vector measurement values ​​obtained based on the Earth ellipsoid model fitting under the mechanical system specifically include: Edge detection algorithms are used to extract the edge pixels of the Earth's horizon in panoramic infrared Earth images; The edge pixels are mapped to the unit gaze vector using a camera inverse projection model; Based on the quadratic constraint of the Earth ellipsoid, the least squares method is used to fit and solve all unit line-of-sight vectors to obtain the geocentric vector measurement value.

7. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 1, characterized in that, The residual consistency evaluation index in S3 is Mahalanobis distance, and the calculation process includes: Geocentric vector factor based on the k-th epoch or visual feature factors Calculate its corresponding covariance matrix. ;in, It is an observed factor variable. Time represents the covariance matrix corresponding to the geocentric vector factor. Time represents the covariance matrix corresponding to the visual feature factors; Solve for the normalized Mahalanobis distance As an indicator for evaluating residual consistency; among which... hour, Represents the geocentric vector factor. hour, Represents visual feature factors.

8. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 7, characterized in that, covariance matrix The expression is: ； in, This represents the observation matrix of the j-th class of observation factors in the k-th epoch. This represents the covariance matrix of the j-th type of factors estimated in the current state. Let represent the observation noise covariance matrix of the j-th class of observation factors in the k-th epoch.

9. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 7, characterized in that, Based on the residual consistency evaluation index, the dynamic weight coefficients of the geocentric vector factor and the visual feature factor are determined using a piecewise function, which is: when At that time, weighting coefficient ; when At that time, weighting coefficient Follow Increases while decreasing smoothly; when At that time, weighting coefficient The corresponding factor will no longer participate in the current optimization; in, and These are the lower threshold and the upper threshold, respectively.

10. The method for determining the attitude and orbit of micro / nano satellites based on adaptive weighted sliding window factor graph optimization according to claim 1, characterized in that, The objective function expression in S4 is as follows: ； in, This represents the optimal estimate of the state variables within the sliding window of the k-th epoch. Indicates the first orbital dynamics factor of epoch The covariance matrix, Indicates the first Epochal gyro preintegration factor The covariance matrix, Indicates the first Zero bias factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. Geocentric vector factor of epoch The covariance matrix, This indicates the th element after dynamic weighting. The i-th visual feature factor of the epoch The covariance matrix, Denotes the prior factor covariance matrix. Indicates the first The number of visual feature points involved in the optimization of each epoch.

Images