CN108279010A - A kind of microsatellite attitude based on multisensor determines method - Google Patents

Abstract

A multi-sensor based micro-satellite attitude determination method uses MEMS gyroscopes, magnetometers, and sun sensors as attitude determination sensors, and realizes attitude determination of micro-satellites based on an improved federated filtering algorithm. The invention fully considers the application background of micro-satellites, and selects sensors with the advantages of low cost, small volume, and low power consumption as the attitude sensor; the MEMS gyroscope is the main attitude sensor, and integrates to obtain preliminary attitude determination information. The magnetometer and the sun sensor are used as auxiliary attitude sensors to correct the attitude determination results of the MEMS gyroscope in time. The invention proposes an improved federated filtering algorithm to realize the attitude calculation and information fusion based on the above sensors, which effectively improves the attitude determination accuracy of the micro-satellite attitude determination system, and also has the advantages of good fault tolerance and real-time performance.

Description

A multi-sensor based micro-satellite attitude determination method

technical field

The invention relates to a method for determining the attitude of a tiny satellite.

Background technique

Attitude Determination and Control System (ADCS) is one of the prerequisites for the normal operation of micro-satellites.

Currently commonly used satellite attitude sensors include gyroscopes, star sensors, sun sensors, and magnetometers. Among them, the gyroscope is the most commonly used as an inertial sensor, and the satellite attitude can be solved by integral. Traditional high-precision gyroscopes are not suitable for micro-satellites due to their high price and high power consumption; MEMS gyroscopes have advantages in terms of size, power consumption, and cost, but their accuracy is not ideal, with large drift and error accumulation over time. shortcoming. Solving the attitude determination accuracy of the micro-satellite attitude determination system based on MEMS sensors is the current primary task. The magnetometer determines the attitude by measuring the geomagnetic field vector, which is generally suitable for low- and medium-orbit satellites; but it is easily interfered by factors such as the satellite's own residual magnetism, and the attitude determination accuracy is medium. The sun sensor is used to measure the sun vector. It has the advantages of large field of view, clear outline, low power consumption, light weight, etc., and has medium accuracy; but it cannot work normally when the satellite is in the non-visible light area. Therefore, in order to ensure the attitude determination accuracy and availability of the micro-satellite attitude determination system, it is necessary to combine multiple attitude sensors for attitude determination.

The most commonly used attitude estimation algorithm is the Extended Kalman Filter (EKF). The basic idea is to linearize the nonlinear system and then perform Kalman filtering. However, from the principle of extended Kalman filtering, it can be found that it has a better filtering effect on weak nonlinear systems; when the nonlinearity of the system is relatively strong, there may be multiple extreme points, which may easily lead to slow convergence or even divergence. And the amount of calculation is large. Carlson proposed the concept of Federate Filter (FF), in which each sub-filter operates in parallel to perform state estimation, and the main filter realizes information fusion and outputs the final result. Federated Kalman filtering has better fault tolerance, and compared with traditional centralized Kalman filtering, it requires less computation. At first, it was mainly proposed for the application of fault-tolerant integrated navigation system, and then many scholars began to apply it to the field of satellite attitude determination.

For the application background of micro-satellites, the selection of micro-satellite attitude sensors needs to meet the requirements of volume, cost, power consumption, accuracy, etc. Trade-offs between various factors such as stability.

Contents of the invention

The technical solution of the present invention is to overcome the deficiencies of the prior art and provide a method for determining the attitude of micro-satellites based on multiple sensors. The attitude determination of micro-satellites with small-volume, low-power sensors has the advantages of high precision, good fault tolerance, good stability, and low calculation load.

The technical solution of the present invention is: a kind of micro-satellite attitude determining method based on multi-sensor, comprises the steps:

(1) Use MEMS gyroscope to obtain angular velocity measurement value Establish MEMS gyro error model;

(2) According to the measured value of the geomagnetic field vector obtained from the measurement and the calculated geomagnetic vector reference value B _o to establish the magnetometer error model;

(3) Measure and obtain the measured value of the sun vector Using the solar ephemeris to derive the solar vector reference value S _o , and establish the solar sensor error model;

(4) Establish the state equation according to the micro-satellite attitude kinematics equation and the MEMS gyro error model established with the quaternion as the attitude parameter;

(5) According to the magnetometer error model and the solar sensor error model, respectively establish the observation equations of sub-filter S1 and sub-filter S2;

(6) Carry out adaptive distribution to the noise variance matrix and the mean square error matrix of sub-filter S1 and sub-filter S2 according to information distribution factor β _i ;

(7) The estimated attitude quaternion is obtained by integrating the MEMS gyro angular velocity measurement value

(8) Use the state equation established in step (4) to calculate the one-step state prediction at time k and The estimated mean square error matrix P _k/k-1 ;

(9) Calculate the information distribution factor β _i , and reset the mean square error matrix P _i,k/k-1 respectively corresponding to the sub-filters S1 and S2 to the sub-filters S1 and S2;

(10) The sub-filters S1 and S2 are measured and updated according to the observation information of the magnetometer and the solar sensor, and the Kalman filter gain K _{i,k of the} sub-filter and the state estimation are respectively calculated and the mean square error P _i,k ;

(11) State estimation for sub-filters S1 and S2 Perform data fusion to obtain the final state estimate and the mean square error matrix P _k ;

(12) Estimating the state The quaternion used to correct the estimated pose with a random constant drift estimate Obtain the final attitude quaternion q _k and the random constant drift value b _k .

The concrete method of described step (1) is:

When the installation axis of the MEMS gyro is consistent with the inertial axis of the microsatellite, the output of the MEMS gyro is the measured value of the angular velocity of the body coordinate system b of the microsatellite relative to the inertial coordinate system j

The MEMS gyro error model is established as:

Among them, η _g and η _b are white noise, ω _bj is the real angular velocity, and b is a random constant value drift.

The concrete method of described step (2) is:

When the installation axis of the magnetometer is consistent with the coordinate system of the micro-satellite body, the measured value of the geomagnetic field vector in the satellite body coordinate system is obtained The magnetometer error model is established as:

where V _B is the measurement noise of the magnetometer, is the conversion matrix from the orbital coordinate system to the micro-satellite body coordinate system, and B _o is the reference value of the geomagnetic field vector in the orbital coordinate system.

The twelfth-generation international geomagnetic reference model IGRF-12 is used to calculate the reference value of the geomagnetic field vector in the Northeast geographic coordinate system at the current orbital position of the microsatellite, and converted to the orbital coordinate system to obtain the geomagnetic field vector reference value in the orbital coordinate system value B _o .

The concrete method of described step (3) is:

Select the two-axis digital sun sensor and determine the installation matrix of the two-axis digital sun sensor, align the aiming axis of the two-axis digital sun sensor with the -Y axis of the micro-satellite body coordinate system, and obtain the sun vector in the satellite body coordinate system Measurements The mounting matrix of the two-axis digital sun sensor for:

Use the solar ephemeris to calculate the unit sun vector reference value in the J2000 geocentric inertial coordinate system at the current moment, and convert it to the orbital coordinate system through the conversion matrix to obtain the sun vector reference value S _o ;

The sun sensor error model is established as:

Among them, V _s is the measurement noise of the sun sensor.

The concrete method of described step (4) is:

Establish micro-satellite attitude kinematics equation:

in, Δq=[Δq ₁ Δq ₂ Δq ₃ ] ^T ,

The error quaternion Δq is the real attitude quaternion q and the estimated attitude quaternion the error between

Among them, Δq ₁ , Δq ₂ , Δq ₃ are the vector part of the error quaternion;

The random constant drift error Δb is the true random constant drift b and the estimated random constant drift value the error between

Set the six-dimensional state variable as X=[Δq Δb] ^T , use the MEMS gyroscope error model and the attitude kinematics equation to derive the state equation:

Among them, t is the time parameter, the matrix matrix

Matrix W = [η _g η _b ] ^T ;

The state equation is discretized, and the subscripts k-1 and k represent the discretized time k-1 and k respectively, and T is the time interval, then:

X _k = Φ _k/k-1 X _k-1 +Γ _k-1 W _k-1 ,

discretization coefficient

discretization coefficient

Wherein, F _k-1 represents the matrix F(t) at time k-1 after discretization; G _k-1 represents the matrix G(t) at time k-1 after discretization.

The concrete steps of described step (5) are as follows:

The observation equation of sub-filter S1 is obtained by using the magnetometer error model and state variables:

ΔB _b = H _b X + V _B ;

Among them, the geomagnetic field vector error value ΔB _b is the geomagnetic field vector measurement value of the magnetometer and the estimated value of the geomagnetic field vector based on the IGRF-12 geomagnetic model The error between: matrix

Using the solar sensor error model and state variables, the observation equation of the sub-filter S2 is derived:

ΔS _b = H _s X + V _s ;

Among them, the sun vector error value ΔS _b is the sun vector measurement value obtained by the sun sensor and the solar vector estimate calculated based on the solar ephemeris The error between: matrix

The formula for calculating the information allocation factor β _i is as follows:

Among them, S ₁ represents the reciprocal of the sum of the absolute values of the sub-filter S1 mean square error matrix P ₁ eigenvalues; S ₂ represents the reciprocal of the sum of the absolute values of the sub-filter S2 mean square error matrix P ₂ eigenvalues.

In said step (8), one-step state prediction at time k and The calculation formula of the estimated mean square error matrix P _k/k-1 is as follows:

The mean square error matrix P _i,k/k-1 in the step (9):

The step (10) neutron filter Kalman filter gain K _i,k , state estimation And the calculation formula of mean square error P _i,k is as follows:

P _i,k = (IK _i,k H _i,k )P _i,k/k-1 ;

Final state estimation in described step (11) And the calculation formula of the mean square error matrix P _k is as follows:

The calculation formula of final attitude quaternion q _k and random constant drift value b _k in the described step (12) is as follows:

The beneficial effect of the present invention compared with prior art is:

(1) The present invention realizes the micro-satellite attitude determination based on multiple low-cost, small-volume, low-power sensors, which is easy for engineering implementation and contributes to the miniaturization, autonomy, and localization of the micro-satellite attitude determination system.

(2) Compared with the traditional concentrated Kalman filter attitude determination method, the present invention is easy to be detected and isolated without polluting the entire attitude determination system due to the parallel operation of two sub-filters when a sensor fails or does not work , which has the advantage of good fault tolerance.

(3) The improved federated filtering algorithm of the present invention is comparatively loose to initial filtering parameter selection; Even under the initial attitude error situation of large angle, also can realize fast convergence and guarantee fixed attitude accuracy;

(4) The state variable dimension of the improved federated filtering algorithm of the present invention is 6 dimensions, and the state variables are less, and the two sub-filters do not have independent state variables and only perform measurement and update, which has the advantages of small amount of calculation and good real-time performance .

Description of drawings

Fig. 1 is the basic principle diagram of the attitude determination method based on multi-sensor of the present invention;

Fig. 2 is the basic principle diagram of improved federated filtering algorithm in the present invention;

Fig. 3 is a schematic diagram of an improved federated filtering algorithm state equation in the present invention;

Fig. 4 is a principle diagram of the observation equation of the improved federated filtering algorithm in the present invention.

Detailed ways

As shown in Figure 1, it is the basic principle diagram of the attitude determination method based on multiple sensors in the present invention. The MEMS gyroscope is selected as the main attitude sensor, and the magnetometer and sun sensor are used as the auxiliary attitude sensors. An improved federation is designed. The filtering algorithm performs attitude calculation and information fusion.

As shown in Figure 2, the basic principle diagram of the improved federated filtering algorithm in the present invention, the improved federated filtering algorithm in the present invention is mainly composed of two sub-filters S1, S2 and a main filter M1, and adopts a fusion-feedback structure. The sub-filter S1 completes the measurement update of the magnetometer, and the sub-filter S2 completes the measurement update of the sun sensor; while in the main filter, the MEMS gyro angular velocity integration, state prediction, data fusion and final attitude determination results are realized. output. The two sub-filters S1 and S2 do not contain independent state variables. Therefore, although they operate in parallel, they only perform measurement updates; while the prediction process is completed in the main filter M1; this is different from the parallel filter fusion algorithm. After the main filter M1 completes the data fusion based on the measurement information of the sub-filters, it outputs the result of the attitude determination and at the same time feeds back the two sub-filters according to the principle of information distribution.

A multi-sensor based micro-satellite attitude determination method, select low-cost, low-power, small-volume sensors as attitude sensors, including MEMS gyroscopes, magnetometers and sun sensors; design an improved federated filtering algorithm to achieve The attitude determination of micro-satellites based on multi-sensors, the specific steps are as follows:

(1) The MEMS gyroscope is used as the main attitude sensor of the microsatellite, and the MEMS gyroscope error model suitable for the attitude determination algorithm of the microsatellite is established, and the estimated attitude quaternion is obtained by integrating the angular velocity measurement value; the specific steps are as follows:

(1.1) MEMS gyroscope error model mainly considers random constant value drift and first-order time correlation drift among the present invention, above-mentioned error is modeled and simplified;

(1.2) When the installation axis of the MEMS gyroscope is consistent with the inertial axis of the microsatellite, the angular velocity measurement value of the body coordinate system of the microsatellite relative to the J2000 earth-centered inertial coordinate system is obtained Through coordinate conversion and considering the orbital angular velocity ω _oj , the angular velocity of the satellite body coordinate system relative to the orbital coordinate system is obtained For attitude calculation;

(1.3) The MEMS gyro Allan variance analysis results are used as the filter parameter reference value of the improved federated filter algorithm to ensure the convergence and precision of the filter.

(2) The magnetometer is used as an auxiliary attitude sensor to obtain the geomagnetic vector measurement value, combined with the twelfth generation IGRF-12 model to obtain the geomagnetic vector reference value, and establish a magnetometer error model suitable for the improved federated filtering algorithm; specific steps as follows:

(2.1) Select the twelfth generation international geomagnetic reference model IGRF-12 to calculate the geomagnetic field vector reference value under the Northeast geographic coordinate system under the current orbital position of the microsatellite;

(2.2) When the installation axis of the magnetometer is consistent with the coordinate system of the micro-satellite body, the measured value of the geomagnetic field vector in the body coordinate system of the micro-satellite is obtained Establish a magnetometer error model suitable for the improved federated filtering algorithm.

(3) The sun sensor is used as an auxiliary attitude sensor to obtain the measured value of the sun vector under a specific installation matrix, and the reference value of the sun vector is derived by using the solar ephemeris, and the error model of the sun sensor suitable for the improved federated filtering algorithm is established; the specific steps as follows:

(3.1) Select a two-axis digital sun sensor and determine its installation matrix, align the aiming axis with the -Y axis of the microsatellite body coordinate system, and obtain the measured value of the sun vector in the microsatellite body coordinate system

(3.2) Utilize the solar ephemeris to calculate the unit solar vector reference value under the J2000 geocentric inertial coordinate system at the current moment, and convert it to the orbital coordinate system through the conversion matrix;

(3.3) Establish a sun sensor error model suitable for the improved federated filtering algorithm for measurement update.

(4) Obtain the 6-dimensional state variable by defining the error quaternion Δq and the random constant value drift error Δb, and combine the MEMS gyroscope error model and attitude kinematics equation in step (1) to establish the state equation, based on which the state equation can effectively reduce The amount of calculation of the filter and improve the stability of the filter;

(5) establish the observation equation of two sub-filters of the federated filtering algorithm in conjunction with the magnetometer in the step (2), the step (3) and the sun sensor error model;

(6) Select the reciprocal of the sum of the absolute values of the eigenvalues of the mean square error matrix as the performance evaluation standard of the sub-filter, and the main filter M1 adaptively allocates the noise variance matrix and the mean square error matrix of the sub-filters S1 and S2 ;

(7) Use the improved federated filtering algorithm to calculate and obtain the micro-satellite attitude results, and realize the micro-satellite attitude calculation and information fusion based on the above multi-sensors. The main process of the improved federated filtering algorithm includes six main parts: microsatellite quaternion attitude kinematics equation integral solution estimation quaternion, state prediction, information distribution, measurement update, data fusion and result output.

1. MEMS gyro error model

In the present invention, the J2000 geocentric inertial coordinate system is represented by j, and its coordinate origin is located at the earth's barycenter; the Z axis is perpendicular to the earth's equatorial plane, coincides with the earth's rotation axis and points to the North Pole; the X axis points to the vernal equinox in the earth's equatorial plane; the Y axis and The X-axis and Z-axis are determined according to the right-hand rule.

The coordinate system of the microsatellite body is represented by b, and the coordinate origin O is the center of mass of the microsatellite; the three axes are respectively parallel to the principal inertia axis of the microsatellite, the X axis points forward along the longitudinal axis, and the Z axis points perpendicularly to the X axis in the longitudinal symmetry plane Ground, Y, X, and Z axes are determined according to the right-hand rule, pointing to the right wing of the satellite.

In the present invention, the MEMS gyroscope is fixedly connected with the micro-satellite body coordinate system, and outputs the angular velocity measurement value of the micro-satellite body coordinate system b relative to the J2000 earth-centered inertial coordinate system j Considering that the time-dependent drift is an attenuation process; when the satellite is in a stable state of three-axis attitude, the deterministic part of the time-related drift is much smaller than the random constant drift; therefore, the random part of the time-related drift is regarded as measurement noise. In conjunction with the commonly used gyro error model, the MEMS gyro error model in the present invention is simplified as:

Among them, η _g and η _b are white noise, ω _bj is the real angular velocity, and b is a random constant value drift.

The present invention uses the actual MEMS gyro Allan variance analysis result as the filter parameter reference value of the federated filter algorithm to ensure the convergence and precision of the filter.

2. Magnetometer error model

In the present invention, the magnetometer is installed along the coordinate system of the micro-satellite body, and outputs the measured value of the geomagnetic field vector under the body coordinate system of the micro-satellite Use the twelfth-generation international geomagnetic reference field model IGRF-12 to calculate the geomagnetic field vector reference value in the Northeast geographic coordinate system at the current position; and use a series of coordinate system conversion matrices to convert it into the reference value B in the orbital coordinate system _o . Among them, the geographic coordinate system of the North East is represented by NED, the origin can be any point on the earth, the X axis points to the local north along the geographic longitude tangent, the Y axis points to the local east along the geographic latitude tangent, and the Z axis and the X axis, The Y axis points to the ground according to the right hand rule. The orbital coordinate system is denoted by o, the origin of which is the center of mass of the microsatellite, the Z-axis points to the center of the earth, the X-axis is perpendicular to the Z-axis in the orbital plane and points to the direction of travel of the satellite, and the Y-axis, X-axis, and Z-axis are determined according to the right-hand rule , parallel to the normal to the orbital plane.

In view of the complex error model will increase the calculation amount of the filtering algorithm, and easily lead to the instability of the filter; in the present invention, the magnetometer error model is simplified as:

Among them, V _B is the measurement noise of the magnetometer; is the conversion matrix from the orbital coordinate system to the microsatellite body coordinate system.

3. Sun sensor error model

In the present invention, the two-axis digital sun sensor is selected as the attitude sensor; the solar sensor's aiming axis Z axis is aligned with the microsatellite body coordinate system-Y axis, and the sun sensor's measuring axis X axis and Y axis are respectively aligned with the microsatellite body The X-axis and Z-axis of the coordinate system are consistent, and the solar sensor installation coordinate system is represented by s; therefore, the solar sensor installation matrix for:

The sun sensor can get the sun vector measurement value in the installation coordinate system Using the transformation matrix to convert it into the solar vector measurement value in the microsatellite body coordinate system The solar ephemeris can be used to calculate the solar vector reference value S _j in the J2000 geocentric inertial coordinate system at the current moment, and convert it to S _o in the orbital coordinate system. Comprehensively considering the filter accuracy, stability and calculation amount, the error model of the sun sensor in the present invention is described as:

Among them, V _s is the measurement noise of the sun sensor.

4. Equation of state

As shown in Figure 3, the present invention selects the quaternion as the attitude parameter to establish the microsatellite attitude kinematics equation (6), and defines the error quaternion Δq as the real attitude quaternion q and the estimated attitude quaternion The error between:

Δq＝[Δq ₀ Δq ₁ Δq ₂ Δq ₃ ] ^T (8)

Among them, Δq ₀ is the scalar part of the error quaternion, and Δq ₁ , Δq ₂ , and Δq ₃ are the vector parts of the error quaternion.

The attitude error of the satellite in the three-axis stable state is very small, so Δq ₀ ≈1; the dimension reduction of the error quaternion is three-dimensional Δq=[Δq ₁ Δq ₂ Δq ₃ ] ^T . At the same time, the random constant drift error Δb is defined as the true random constant drift b and the estimated random constant drift value error between. Accordingly, the 6-dimensional state variable of the improved federated filtering algorithm is set as X=[Δq Δb] ^T . Combining the MEMS gyroscope error model and the attitude kinematics equation, the state equation is derived:

Among them, t is the time parameter, the matrix matrix Matrix W = [η _g η _b ] ^T .

The state equation (9) is further discretized, and the subscripts k-1 and k represent the discretized time k-1 and k time respectively, where k is a positive integer and T is the time interval, then:

X _k ＝Φ _k/k-1 X _k-1 +Γ _k-1 W _k-1 (10)

discretization coefficient

discretization coefficient

Wherein, F _k-1 represents the matrix F(t) at time k-1 after discretization; G _k-1 represents the matrix G(t) at time k-1 after discretization.

5. Observation equation

As shown in Figure 4, the geomagnetic field vector error value ΔB _b is defined as the magnetometer geomagnetic field vector measurement value and the estimated value of the geomagnetic field vector based on the IGRF-12 geomagnetic model The error between:

The observation equation of sub-filter S1 can be obtained by combining the magnetometer error model and state variables:

ΔB _b ＝H _b X+V _B (14)

Among them, the matrix

Similarly, define the sun vector error value ΔS _b as the sun vector measurement value obtained by the sun sensor and the solar vector estimate calculated based on the solar ephemeris The error between:

In conjunction with the solar sensor error model proposed by the present invention and the state variable setting, the observation equation of the sub-filter S2 can be derived:

ΔS _b = H _s X + V _s (16)

Among them, the matrix

6. Information distribution

As shown in Figure 2, the present invention adopts an improved federated filtering algorithm with a fusion-feedback structure, adaptively allocates the noise variance matrix Q and the mean square error matrix P according to the information allocation factor, and feeds back the sub-filters S1 and S2. The adaptive calculation of the information distribution factor β _i is realized based on the eigenvalues of the mean square error matrix P. Since the eigenvalues can be positive or negative, the calculation is performed as the reciprocal of the sum of the absolute values of the eigenvalues:

i=1,2; Among them, S ₁ represents the reciprocal of the sum of the absolute values of the sub-filter S1 mean square error matrix P ₁ eigenvalues; S ₂ represents the sum of the absolute values of the sub-filter S2 mean square error matrix P ₂ eigenvalues and the reciprocal of .

7. Basic filtering process

As shown in Fig. 2, the improved federated filtering algorithm in the present invention includes six main functions of microsatellite quaternion attitude kinematics equation integral solution estimation quaternion, state prediction, information distribution, measurement update, data fusion and result output part.

The subscripts k-1 and k represent time k-1 and k respectively, the subscript i represents the sub-filter, and k-1/k represents the relevant parameters at time k estimated from time k-1. The calculation steps are as follows:

(1) Integral: Use the MEMS gyro angular velocity measurement value integration to obtain the estimated attitude quaternion at time k

(2) State prediction: According to the state equation in the main filter M1, the estimated value of the state at time k-1 Compute the one-step state prediction at time k and its estimated mean square error matrix P _k/k-1 ;

Among them, Q _k-1 is W _k-1 noise variance matrix.

(3) Information distribution: the main filter M1 calculates the information distribution factor β _i , and resets the mean square error matrix P _i,k/k-1 to the sub-filters S1 and S2:

(4) Measurement update: The sub-filters S1 and S2 perform measurement update based on the observation information of the magnetometer and solar sensor, and calculate the sub-filter Kalman filter gain K _i,k and the sub-filter state at time k respectively estimate and the mean square error P _i,k :

P _i,k ＝(IK _i,k H _i,k )P _i,k/k-1 (23)

Among them, R _i,k and Z _i,k are respectively the observation noise variance matrix and observation information of the sub-filter at time k, and H _i,k are the structural parameters of the observation equation of the sub-filter.

Due to the small error, In the present invention, the state estimation Simplified to (22), it helps to improve the stability and convergence of federated filtering.

(5) Data fusion: the state estimation of the two sub-filters by the main filter M1 Perform data fusion to obtain the final state estimate and the mean square error matrix P _k :

(6) Result output: estimate the state The quaternion used to correct the estimated pose with a random constant drift estimate Obtain the final k-time attitude quaternion q _k and random constant drift value b _k :

in, Estimated by time k-1, it is considered in the present invention that

The content that is not described in detail in the description of the present invention belongs to the known technology of those skilled in the art.

Claims (10)

1. A microsatellite attitude determination method based on multiple sensors is characterized by comprising the following steps:

(1) obtaining angular velocity measurements using MEMS gyroscopesEstablishing an MEMS gyro error model;

(2) vector measurement of earth magnetic field obtained from measurementAnd the calculated geomagnetic vector reference value B_oEstablishing a magnetometer error model;

(3) the sun vector measured value is obtained by measurementSolar vector reference value S is obtained by utilizing solar ephemeris derivation_oEstablishing an error model of the sun sensor;

(4) establishing a state equation according to a microsatellite attitude kinematics equation established by taking quaternion as an attitude parameter and an MEMS gyro error model;

(5) respectively establishing an observation equation of a sub-filter S1 and an observation equation of a sub-filter S2 according to the error model of the magnetometer and the error model of the sun sensor;

(6) factor β is distributed according to information_iAdaptively distributing the noise variance matrix and the mean square error matrix of the sub-filter S1 and the sub-filter S2;

(7) obtaining an estimated attitude quaternion by using the integral of the MEMS gyroscope angular velocity measurement value

(8) Calculating one-step state prediction at the k moment by using the state equation established in the step (4)Andis estimated mean square error matrix P_k/k-1；

(9) Calculate information distribution factor β_iAnd resetting the mean square error matrix P corresponding to the sub-filters S1 and S2 to the sub-filters S1 and S2_i,k/k-1；

(10) The sub-filters S1 and S2 carry out measurement updating according to the observation information of the magnetometer and the sun sensor, and respectively calculate the Kalman filtering gain K of the sub-filters_i,kState estimationAnd mean square error P_i,k；

(11) State estimation for sub-filters S1, S2Performing data fusion to obtain final state estimationAnd a mean square error matrix P_k；

(12) Estimating the stateFor modifying estimated attitude quaternionAnd a random constant drift estimateObtaining the final attitude quaternion q_kAnd a random constant drift value b_k。

2. The method of claim 1, wherein the method comprises:

the specific method of the step (1) comprises the following steps:

when the MEMS gyroscope installation shaft is consistent with the inertial axis of the microsatellite, the MEMS gyroscope outputs to obtain the angular velocity measured value of the microsatellite body coordinate system b relative to the inertial coordinate system j

The MEMS gyro error model is established as follows:

wherein, η_g，η_bIs white noise, omega_bjB is a random constant drift for true angular velocity.

3. A method for multi-sensor based microsatellite attitude determination according to claim 1 or 2 wherein:

the specific method of the step (2) is as follows:

when the mounting shaft of the magnetometer is consistent with the body coordinate system of the microsatellite, the vector measurement value of the geomagnetic field under the body coordinate system of the satellite is obtainedEstablishing a magnetometer error model as follows:

wherein, V_BIs the measurement noise of the magnetometer,is a transformation matrix from the orbital coordinate system to the microsatellite body coordinate system, B_oIs the earth magnetic field vector reference value in the track coordinate system.

Selecting a twelfth generation international geomagnetic reference model IGRF-12 for calculating a reference value of a geomagnetic field vector in a geographical coordinate system of the northeast of the earth at the current orbit position of the microsatellite, and converting the reference value into the orbit coordinate system to obtain a geomagnetic field vector reference value B in the orbit coordinate system_o。

4. A method for multi-sensor based microsatellite attitude determination as set forth in claim 3 wherein:

the specific method of the step (3) is as follows:

selecting a two-axis digital sun sensor, determining an installation matrix of the two-axis digital sun sensor, and enabling the aiming axis of the two-axis digital sun sensor to be consistent with the coordinate system-Y axis of the microsatellite body to obtain the measured value of the sun vector under the coordinate system of the microsatellite bodyInstallation matrix of two-axis digital sun sensorComprises the following steps:

calculating a unit sun vector reference value under a current time J2000 geocentric inertial coordinate system by using a sun ephemeris, and converting the unit sun vector reference value into an orbit coordinate system through a conversion matrix to obtain a sun vector reference value S_o；

The error model of the sun sensor is established as follows:

wherein, V_sIs the measurement noise of the sun sensor.

5. The method of claim 4, wherein the method comprises:

the specific method of the step (4) comprises the following steps:

establishing a microsatellite attitude kinematics equation:

wherein,Δq＝[Δq₁Δq₂Δq₃]^T，

the error quaternion delta q is a real attitude quaternion q and an estimated attitude quaternionThe error between;

wherein, Δ q₁，Δq₂，Δq₃A vector portion that is an error quaternion;

the random constant drift error delta b is the true random constant drift b and the estimated random constant drift valueThe error between;

setting a six-dimensional state variable to X ═ Δ q Δ b]^TAnd deducing to obtain a state equation by utilizing an MEMS gyro error model and an attitude kinematics equation:

where t is a time parameter, matrixMatrix arrayMatrix W ═ η_gη_b]^T；

Discretizing the state equation, and respectively expressing discretization k-1 time by subscripts k-1 and k, wherein the k time is T time interval, then:

X_k＝Φ_k/k-1X_k-1+Γ_k-1W_k-1，

coefficient of discretization

Coefficient of discretization

Wherein, F_k-1A matrix F (t) representing the discrete k-1 time; g_k-1The matrix G (t) representing the discrete k-1 time instant.

6. A method for multi-sensor based microsatellite attitude determination according to claim 4 or 5 wherein:

the specific steps of the step (5) are as follows:

the observation equation of the sub-filter S1 is obtained by using the error model of the magnetometer and the state variables:

ΔB_b＝H_bX+V_B；

wherein the geomagnetic vector error value Delta B_bAs vector measurement of the geomagnetic field of a magnetometerAnd a geomagnetic field vector estimation value obtained based on an IGRF-12 geomagnetic modelError between:matrix array

And (3) deducing an observation equation of the sub-filter S2 by using the sun sensor error model and the state variable:

ΔS_b＝H_sX+V_s；

wherein the sun vector error value Delta S_bSun vector measurements for sun sensorsAnd a solar vector estimated value calculated based on solar ephemerisThe error between:matrix array

7. The method of claim 6, wherein the information distribution factor β is a factor related to the attitude of the microsatellite_iThe calculation formula is as follows:

wherein S is₁Representing a sub-filter S1 mean square error matrix P₁The reciprocal of the sum of the absolute values of the eigenvalues; s₂Representing a sub-filter S2 mean square error matrix P₂The reciprocal of the sum of the absolute values of the eigenvalues.

8. A method for multi-sensor based microsatellite attitude determination as set forth in claim 7 wherein: the one-step state prediction at the k time in the step (8)Andis estimated mean square error matrix P_k/k-1The calculation formula of (a) is as follows:

the mean square error matrix P in the step (9)_i,k/k-1：

9. A method for multi-sensor based microsatellite attitude determination according to claim 7 or 8 wherein: the neutron filter Kalman filtering gain K in the step (10)_i,kState estimationAnd mean square error P_i,kThe calculation formula is as follows:

P_i,k＝(I-K_i,kH_i,k)P_i,k/k-1；

the final state estimation in said step (11)And a mean square error matrix P_kThe calculation formula of (a) is as follows:

10. the method of claim 9, wherein the method comprises determining the attitude of the microsatellite based on multiple sensorsIn the following steps: the final attitude quaternion q in the step (12)_kAnd a random constant drift value b_kThe calculation formula of (a) is as follows: