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

A kind of microsatellite attitude based on multisensor determines method Download PDF

Info

Publication number
CN108279010A
CN108279010A CN201711362902.XA CN201711362902A CN108279010A CN 108279010 A CN108279010 A CN 108279010A CN 201711362902 A CN201711362902 A CN 201711362902A CN 108279010 A CN108279010 A CN 108279010A
Authority
CN
China
Prior art keywords
microsatellite
attitude
matrix
error
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201711362902.XA
Other languages
Chinese (zh)
Inventor
倪枫
韩逸飞
李文杰
刘肖姬
陈路
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Microelectronic Technology Institute
Mxtronics Corp
Original Assignee
Beijing Microelectronic Technology Institute
Mxtronics Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Microelectronic Technology Institute, Mxtronics Corp filed Critical Beijing Microelectronic Technology Institute
Priority to CN201711362902.XA priority Critical patent/CN108279010A/en
Publication of CN108279010A publication Critical patent/CN108279010A/en
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01CMEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/005Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 with correlation of navigation data from several sources, e.g. map or contour matching

Abstract

A kind of microsatellite attitude based on multisensor determines method, and sensor is determined as posture using MEMS gyro, magnetometer, sun sensor, realizes that the posture of microsatellite determines based on a kind of modified Federated Filtering.The present invention has fully considered the application background of microsatellite, select to have many advantages, such as low cost, small size, low-power consumption sensor as attitude sensor;MEMS gyro is main attitude sensor, and integral obtains preliminary posture and determines information.Magnetometer and sun sensor are as assisting what attitude sensor corrected MEMS gyro in time to determine appearance result.The present invention proposes a kind of modified Federated Filtering and realizes that the attitude algorithm based on the sensor is merged with information, effectively increases the accuracy of attitude determination that microsatellite attitude determines system, also has many advantages, such as that fault-tolerance is good, real-time is good.

Description

A kind of microsatellite attitude based on multisensor determines method
Technical field
The present invention relates to a kind of microsatellite attitudes to determine method.
Background technology
Pico-satellite (Attitude Determination and Control System, ADCS) is One of the premise of microsatellite normal work.
Currently used satellite attitude sensor has gyro, star sensor, sun sensor, magnetometer etc..Wherein, top Spiral shell is the most commonly used as inertial sensor, and the attitude of satellite can be solved using integral.Traditional high accuracy gyroscope instrument is high due to price It is expensive, power consumption is big etc., and reasons are not particularly suited for microsatellite;MEMS gyroscope takes advantage in volume, power consumption and cost etc., But precision is not ideal enough, has the shortcomings that drift is big, error is accumulated at any time.Solve the microsatellite based on MEMS sensor Attitude and heading reference system accuracy of attitude determination is current top priority.Magnetometer determines posture by measuring geomagnetic fieldvector, general to be applicable in In LEO satellite;But it is easily interfered by factors such as satellite itself remanent magnetism, accuracy of attitude determination is medium.Sun sensor is for measuring Solar vector has many advantages, such as that visual field is big, clear-cut, low in energy consumption, light weight, and precision is medium;But when satellite be in it is non-can Cisco unity malfunction when light-exposed area.Thus, to ensure that microsatellite attitude determines the accuracy of attitude determination and availability of system, need Joint many attitude sensor carries out determining appearance.
The most commonly used Attitude estimation algorithm be extended Kalman filter (Extended Kalman Filter, EKF), basic thought is that nonlinear system linearization then is carried out Kalman filtering.But from the original of Extended Kalman filter It can be found that it has preferable filter effect to weakly non-linear system in reason;When system it is non-linear more strong when, can There can be multiple extreme points, be easy to cause convergence slowly, or even diverging, and calculation amount is larger.Carlson proposes federal filter The concept of wave (Federate Filter, FF), each subfilter concurrent operation carry out state estimation, and senior filter is realized Information fusion output final result.Federated Kalman Filtering has preferable fault-tolerance, with traditional centralized Kalman filter phase Than calculation amount is smaller.Originally primarily directed in fault-tolerance combined navigation systematic difference propose, later many scholars start by It determines field applied to the attitude of satellite.
For the application background of microsatellite, need to meet in the selection of microsatellite attitude sensor volume, at Various demands such as sheet, power consumption, precision;And based on the posture of multisensor determine algorithm must precision, calculation amount, Weighed between the various aspects factor such as filter stability.
Invention content
The technology of the present invention solves the problems, such as:In place of overcome the deficiencies in the prior art, provide a kind of based on multisensor Microsatellite attitude determines that method, this method are suitable for microsatellite attitude determination and control system, realize based on multiple low Cost, small size, low-power consumption sensor microsatellite attitude determine, and with high-precision, zmodem, good stability, The advantages such as calculation amount is low.
Technical solution of the invention is:A kind of microsatellite attitude based on multisensor determines method, including such as Lower step:
(1) MEMS gyro is utilized to obtain angular velocity measurement valueEstablish MEMS gyro error model;
(2) according to the geomagnetic fieldvector measured value for measuring acquisitionAnd the ground magnetic vector reference value B being calculatedo, build Vertical magnetometer error model;
(3) it measures and obtains sun vector measurement valueIt is derived by solar vector reference value S using sun ephemeriso, build Vertical sun sensor error model;
(4) it is missed according to the microsatellite attitude kinematical equation and MEMS gyro established using quaternary number as attitude parameter Poor model foundation state equation;
(5) according to magnetometer error model and sun sensor error model, subfilter S1 is established respectively and son filters The observational equation of device S2;
(6) according to information distribution factor βiTo the noise variance matrix and mean square error of subfilter S1 and subfilter S2 Poor matrix carries out self-adjusted block;
(7) it integrates to obtain estimation attitude quaternion using MEMS gyro angular velocity measurement value
(8) state equation established in step (4) is utilized to calculate one step status predication of k momentAndEstimation Square Error matrix Pk/k-1
(9) information distribution factor β is calculatedi, and it is corresponding to subfilter S1, S2 resetting subfilter S1, S2 Square error matrix Pi,k/k-1
(10) subfilter S1, S2 carries out measurement update according to the observation information of magnetometer and sun sensor, counts respectively Operator filtering device Kalman filtering gain Ki,k, state estimationAnd mean square error Pi,k
(11) to the state estimation of subfilter S1, S2Data fusion is carried out, final state estimation is obtainedWith And Square Error matrix Pk
(12) by state estimationEstimate attitude quaternion for correctingWith Random Constant Drift estimated valueIt obtains Final attitude quaternion qkWith Random Constant Drift value bk
The specific method of the step (1) is:
When MEMS gyro installation axle is consistent with the microsatellite axes of inertia, MEMS gyro exports to obtain microsatellite ontology seat Angular velocity measurement values of the mark system b relative to inertial coodinate system j
Establishing MEMS gyro error model is:
Wherein, ηg, ηbFor white noise, ωbjFor true angular velocity, b is Random Constant Drift.
The specific method of the step (2) is:
When magnetometer installation axle is consistent with microsatellite body coordinate system, obtains earth's magnetic field under satellite body coordinate system and swear MeasurementEstablishing magnetometer error model is:
Wherein, VBFor the measurement noise of magnetometer,For the conversion square of orbital coordinate system to microsatellite body coordinate system Battle array, BoFor geomagnetic fieldvector reference value under orbital coordinate system.
Select the 12nd generation world earth magnetism reference model IGRF-12 for calculating under microsatellite current orbit position The reference value of magnetic vector in east northeast geographic coordinate system, and convert to orbital coordinate system, obtain earth magnetism under orbital coordinate system Field vector reference value Bo
The specific method of the step (3) is:
It selects two axis digital sun sensors and determines the installation matrix of two axis digital sun sensors, by two number of axle words The guidance axis of sun sensor is consistent with microsatellite body coordinate system-Y-axis, obtains solar vector in satellite body coordinate system Under measured valueThe installation matrix of the two axis digital sun sensorFor:
Using unit solar vector reference value under sun ephemeris computation current time J2000 geocentric inertial coordinate systems, and lead to It crosses transition matrix to convert to orbital coordinate system, obtains solar vector reference value So
Establishing sun sensor error model is:
Wherein, VsFor the measurement noise of sun sensor.
The specific method of the step (4) is:
Establish microsatellite attitude kinematical equation:
Wherein,Δ q=[Δ q1 Δq2 Δq3]T,
Error quaternion Δ q is true attitude quaternion q and estimation attitude quaternionBetween error;
Wherein, Δ q1, Δ q2, Δ q3For the vector section of error quaternion;
Random Constant Drift error delta b is true Random Constant Drift b and estimation Random Constant Drift valueBetween mistake Difference;
Sextuple state variable is set as X=[Δ q Δs b]T, utilize MEMS gyro error model and attitude kinematics side Journey is derived by state equation:
Wherein, t is time parameter, matrixMatrix
Matrix W=[ηg ηb]T
By state equation discretization, with subscript k-1, k indicates the k-1 moment of discretization respectively, and at the k moment, T is between the time Every then having:
Xkk/k-1Xk-1k-1Wk-1,
Discretization coefficient
Discretization coefficient
Wherein, Fk-1The matrix F (t) at the k-1 moment after indicating discrete;Gk-1The matrix G at the k-1 moment after indicating discrete (t)。
The step (5) is as follows:
The observational equation of subfilter S1 is obtained using magnetometer error model and state variable:
ΔBb=HbX+VB
Wherein, geomagnetic fieldvector error amount Δ BbFor magnetometer geomagnetic fieldvector measured valueWith based on IGRF-12 earth magnetism The geomagnetic fieldvector estimated value that model obtainsBetween error:Matrix
Using sun sensor error model and state variable, it is derived by the observational equation of subfilter S2:
ΔSb=HsX+Vs
Wherein, solar vector error amount Δ SbThe sun vector measurement value obtained for sun sensorWith based on the sun The solar vector estimated value that ephemeris is calculatedBetween error:Matrix
Described information distribution factor βiCalculation formula is as follows:
Wherein, S1Indicate subfilter S1 Square Error matrix P1The inverse of the sum of the absolute value of characteristic value;S2Indicate son filter Wave device S2 Square Error matrix P2The inverse of the sum of the absolute value of characteristic value.
One step status predication of k moment in the step (8)AndEstimation Square Error matrix Pk/k-1Meter It is as follows to calculate formula:
Square Error matrix P in the step (9)i,k/k-1
Subfilter Kalman filtering gain K in the step (10)i,k, state estimationAnd mean square error Pi,kMeter It is as follows to calculate formula:
Pi,k=(I-Ki,kHi,k)Pi,k/k-1
Final state estimation in the step (11)And Square Error matrix PkCalculation formula it is as follows:
Final attitude quaternion q in the step (12)kWith Random Constant Drift value bkCalculation formula it is as follows:
Compared with the prior art, the invention has the advantages that:
(1) present invention realizes the microsatellite attitude determination based on multiple low costs, small size, low-power consumption sensor, It is easy to Project Realization, microsatellite attitude is contributed to determine miniaturization, autonomy-oriented, the production domesticization development of system.
(2) present invention is compared with traditional concentrated Kalman filter method for determining posture, due to two sub- filter parallel operations, When a sensor failure or when not working, easily it is detected and is isolated without causing entire attitude determination system to pollute, Have the advantages that fault-tolerance is good.
(3) modified Federated Filtering of the invention is more loose to initial filter parameter selection;Even if there are big angles Under the initial attitude error condition of degree, it can also realize Fast Convergent and ensure accuracy of attitude determination;
(4) the state variable dimension of modified Federated Filtering of the invention is 6 dimensions, and state variable is less, and two Subfilter does not have separate state variable only to carry out measurement update, has calculation amount small, the good advantage of real-time.
Description of the drawings
Fig. 1 is the attitude determination method basic principle figure the present invention is based on multisensor;
Fig. 2 is modified Federated Filtering basic principle figure in the present invention;
Fig. 3 is modified Federated Filtering state equation schematic diagram in the present invention;
Fig. 4 is modified Federated Filtering observational equation schematic diagram in the present invention.
Specific implementation mode
As shown in Figure 1, for the present invention is based on the attitude determination method basic principle figures of multiple sensors, MEMS tops are selected Spiral shell devises a kind of modified connection as main attitude sensor, magnetometer, sun sensor as auxiliary attitude sensor Nation's filtering algorithm carries out attitude algorithm and is merged with information.
Such as modified Federated Filtering basic principle figure in Fig. 2 present invention it is found that the modified federation in the present invention filters Wave algorithm is mainly made of two subfilters S1, S2 and a senior filter M1, and takes fusion-feedback arrangement.Son filter Wave device S1 completes the measurement update of magnetometer, and subfilter S2 completes the measurement update of sun sensor;And in senior filter In, realize MEMS gyro angular speed integral, status predication, data fusion and the final output for determining appearance result.Two son filtering Device S1, S2 do not contain independent state variable, thus, although their concurrent operations, only carry out measuring update;And it is pre- Survey process is then completed in senior filter M1;This point is different from parallel filter blending algorithm.Senior filter M1 is according to son After the measurement information of filter completes data fusion, while appearance result is determined in output, two are fed back to according to information sharing principle Subfilter.
A kind of microsatellite attitude based on multisensor determines method, select low cost, low-power consumption, small size biography Sensor is as attitude sensor, including MEMS gyro, magnetometer and sun sensor;Design a kind of modified federated filter calculation Method realizes that the microsatellite attitude based on multisensor determines, is as follows:
(1) main attitude sensor of the MEMS gyroscope as microsatellite is established and is suitable for microsatellite attitude determination The MEMS gyro error model of algorithm, and integrated by angular velocity measured value and obtain estimation attitude quaternion;Specific steps are such as Under:
(1.1) MEMS gyro error model mainly considers that Random Constant Drift is drifted about with single order time correlation in the present invention, Model simultaneously abbreviation to above-mentioned error;
(1.2) it is opposite that microsatellite body coordinate system is obtained when MEMS gyro installation axle is consistent with the microsatellite axes of inertia In the angular velocity measurement value of J2000 geocentric inertial coordinate systemsIt is converted by coordinate and considers orbit angular velocity ωojIt obtains The angular speed of satellite body coordinate system relative orbit coordinate systemFor attitude algorithm;
(1.3) it is referred to using MEMS gyro Allan the results of analysis of variance as the filtering parameter of modified Federated Filtering Value, ensures the convergence and precision of filter.
(2) magnetometer obtains geomagnetic fieldvector measured value as auxiliary attitude sensor, in conjunction with the 12nd generation IGRF-12 Model obtains ground magnetic vector reference value, establishes the magnetometer error model suitable for modified Federated Filtering;Specific steps It is as follows:
(2.1) select the 12nd generation world earth magnetism reference model IGRF-12 for calculating microsatellite current orbit position Geomagnetic fieldvector reference value of the lower east northeast ground under geographic coordinate system;
(2.2) it is obtained when magnetometer installation axle is consistent with microsatellite body coordinate system under microsatellite body coordinate system Geomagnetic fieldvector measured valueEstablish the magnetometer error model suitable for modified Federated Filtering.
(3) sun sensor obtains the sun vector measurement value under specific installation matrix, profit as auxiliary attitude sensor It is derived by solar vector reference value with sun ephemeris, establishes and is missed suitable for the sun sensor of modified Federated Filtering Differential mode type;It is as follows:
(3.1) two axis digital sun sensors are selected and determine that it installs matrix, guidance axis and microsatellite ontology are sat It is consistent to mark system-Y-axis, obtains the sun vector measurement value under microsatellite body coordinate system
(3.2) unit solar vector under sun ephemeris computation current time J2000 geocentric inertial coordinate systems is utilized to refer to Value, and converted to orbital coordinate system by transition matrix;
(3.3) the sun sensor error model suitable for modified Federated Filtering is established, for measuring update.
(4) 6 dimension state variables are obtained by definition error quaternion Δ q, Random Constant Drift error delta b, in conjunction with step (1) the MEMS gyro error model in establishes state equation with attitude kinematics equations, can effectively be dropped based on the state equation The calculation amount of low filter and the stability for improving filter;
(5) Federated Filtering is established in conjunction with step (2), the magnetometer in step (3) and sun sensor error model The observational equation of two subfilters;
(6) select the inverse of the sum of absolute value of Square Error matrix characteristic value as the performance evaluation mark of subfilter Standard, senior filter M1 carry out self-adjusted block to the noise variance matrix and Square Error matrix of subfilter S1, S2;
(7) it utilizes modified Federated Filtering to calculate and obtains microsatellite attitude as a result, realizing based on above-mentioned more sensings The microsatellite attitude resolving of device is merged with information.Modified Federated Filtering main process includes microsatellite quaternary number Attitude kinematics equations Integration Solving estimates quaternary number, status predication, information distribution, measures update, data fusion and result Export six major parts.
1, MEMS gyro error model
J2000 geocentric inertial coordinate systems are indicated with j in the present invention, and coordinate origin is located at earth centroid;Z axis perpendicular to Earth equatorial plane, and earth rotation overlapping of axles and is directed toward the arctic;X-axis is directed toward the first point of Aries in earth equatorial plane;Y-axis and X Axis, Z axis are determined according to the right-hand rule.
Microsatellite body coordinate system indicates that coordinate origin O is the barycenter of microsatellite with b;Three axis are defended with small respectively The principal moments axis of star is parallel, and X-axis is directing forwardly along the longitudinal axis, and Z axis is perpendicularly oriented to ground, Y and X-axis, Z axis in the vertical plane of symmetry with X-axis It is determined according to the right-hand rule, is directed toward satellite right flank.
In the present invention, MEMS gyro is connected with microsatellite body coordinate system, exports microsatellite body coordinate system b phases To the angular velocity measurement value of J2000 geocentric inertial coordinate systems jIn view of time correlation drift is an attenuation process;When When satellite is in three-axis attitude stabilization state, the certainty part in time correlation drift is more much smaller than Random Constant Drift;Cause And time correlation drift random partial is considered as measurement noise.In conjunction with common gyroscope error model, the MEMS tops in the present invention Spiral shell error model abbreviation is:
Wherein, ηg, ηbFor white noise, ωbjFor true angular velocity, b is Random Constant Drift.
Filtering parameter reference of the present invention using practical MEMS gyro Allan the results of analysis of variance as Federated Filtering Value, ensures the convergence and precision of filter.
2, magnetometer error model
Magnetometer is installed along microsatellite body coordinate system in the present invention, is exported under microsatellite body coordinate system Geomagnetic fieldvector measured valueCurrent location north is calculated using the 12nd generation international geomagnetic reference field model IGRF-12 Geomagnetic fieldvector reference value under eastern ground geographic coordinate system;And it is converted to track using a series of Conversion Matrix of Coordinate and sits Reference value B under mark systemo.Wherein, east northeast geographic coordinate system indicate that origin can be tellurian any point, X with NED Axis is directed toward local north orientation along geographic meridian tangent line, and Y-axis is directed toward local east orientation, Z axis and X-axis, Y-axis along geographical weft tangent line It is directed toward ground according to the right-hand rule.Orbital coordinate system indicates that coordinate origin is the barycenter of microsatellite with o, and Z axis is directed toward the earth's core, X-axis is perpendicularly oriented to satellite direction of travel in orbit plane with Z axis, and Y-axis is determined with X-axis, Z axis according to the right-hand rule, with rail The normal parallel of road plane.
It will increase filtering algorithm calculation amount in view of complicated error model, and be easy to cause the unstable of filter;This hair It is bright middle to be by magnetometer error model abbreviation:
Wherein, VBFor the measurement noise of magnetometer;For the conversion square of orbital coordinate system to microsatellite body coordinate system Battle array.
3, sun sensor error model
It is attitude sensor that two axis digital sun sensors are selected in the present invention;By sun sensor guidance axis Z axis with it is micro- The alignment of moonlet body coordinate system-Y-axis, sun sensor measure axis X-axis, Y-axis respectively with microsatellite body coordinate system X-axis, Z axis is consistent, and sun sensor mounting coordinate system is indicated with s;Therefore sun sensor installs matrixFor:
Sun sensor can obtain sun vector measurement value under mounting coordinate systemIt is converted using transition matrix At sun vector measurement value under microsatellite body coordinate systemCurrent time J2000 can be calculated using sun ephemeris Solar vector reference value S under geocentric inertial coordinate systemj, it is S under conversion to orbital coordinate systemo.Consider filter accuracies, Stability and calculation amount, the error model of sun sensor is described as in the present invention:
Wherein, VsFor the measurement noise of sun sensor.
4, state equation
As shown in figure 3, the present invention selects quaternary number to establish microsatellite attitude kinematical equation (6) as attitude parameter, And it is true attitude quaternion q and estimation attitude quaternion to define error quaternion Δ qBetween error:
Δ q=[Δ q0 Δq1 Δq2 Δq3]T (8)
Wherein, Δ q0For the scalar component of error quaternion, Δ q1, Δ q2, Δ q3For the vector section of error quaternion.
Attitude error very little under satellite three-axis stabilization state, thus Δ q0≈1;Error quaternion dimensionality reduction is three-dimensional Δ q= [Δq1 Δq2 Δq3]T.Meanwhile it defining Random Constant Drift error delta b and being true Random Constant Drift b and estimate random normal Value drift valueBetween error.Accordingly, 6 dimension state variables of modified Federated Filtering are set as X=[Δ q Δs b]T。 In conjunction with MEMS gyro error model and attitude kinematics equations, it is derived by state equation:
Wherein, t is time parameter, matrixMatrix Matrix W=[ηg ηb]T
By state equation (9) further discretization, with subscript k-1, k indicates the k-1 moment after discretization respectively, when k It carves, k is positive integer, and T is time interval, then has:
Xkk/k-1Xk-1k-1Wk-1 (10)
Discretization coefficient
Discretization coefficient
Wherein, Fk-1The matrix F (t) at the k-1 moment after indicating discrete;Gk-1The matrix G at the k-1 moment after indicating discrete (t)。
5, observational equation
As shown in figure 4, defining geomagnetic fieldvector error amount Δ BbFor magnetometer geomagnetic fieldvector measured valueWith based on The geomagnetic fieldvector estimated value that IGRF-12 geomagnetic models obtainBetween error:
The observational equation of subfilter S1 can be obtained in conjunction with magnetometer error model and state variable:
ΔBb=HbX+VB (14)
Wherein, matrix
Similarly, solar vector error amount Δ S is definedbThe sun vector measurement value obtained for sun sensorWith based on The solar vector estimated value that sun ephemeris is calculatedBetween error:
It is arranged in conjunction with sun sensor error model proposed by the invention and state variable, sub- filter can be derived by The observational equation of wave device S2:
ΔSb=HsX+Vs (16)
Wherein, matrix
6, information is distributed
As shown in Fig. 2, the present invention takes the modified Federated Filtering of fusion-feedback arrangement, according to information distribution because Son carries out self-adjusted block to noise variance matrix Q and Square Error matrix P, feeds back to subfilter S1 and S2.Information is distributed Factor-betaiAdaptive polo placement is realized according to the characteristic value of Square Error matrix P.Be divided into view of characteristic value it is positive and negative, with feature The inverse of the sum of value absolute value is calculated:
I=1,2;Wherein, S1Indicate subfilter S1 Square Error matrix P1The inverse of the sum of the absolute value of characteristic value;S2 Indicate subfilter S2 Square Error matrix P2The inverse of the sum of the absolute value of characteristic value.
7, basic filtering
As shown in Figure 2, the modified Federated Filtering in the present invention includes microsatellite quaternary number attitude kinematics side Journey Integration Solving estimates that quaternary number, status predication, information distribution, measurement update, data fusion and result output six are main Part.
Indicate that k-1, k moment, subscript i indicate that subfilter, k-1/k expressions are estimated by the k-1 moment respectively with subscript k-1, k Meter obtains k moment relevant parameters.Process of solution is as follows:
(1) it integrates:It integrates to obtain the k moment using MEMS gyro angular velocity measurement value and estimates attitude quaternion
(2) status predication:According to state equation in senior filter M1, by k-1 moment state estimationsCalculate the k moment A step status predicationAnd its estimation Square Error matrix Pk/k-1
Wherein, Qk-1For Wk-1Noise variance matrix.
(3) information is distributed:Senior filter M1 is into row information distribution factor βiCalculating, and to subfilter S1, S2 reset Square Error matrix Pi,k/k-1
(4) update is measured:Subfilter S1, S2 is measured more according to the observation information of magnetometer and sun sensor Newly, the subfilter Kalman filtering gain K at k moment is calculated separatelyi,k, subfilter state estimationAnd mean square error Pi,k
Pi,k=(I-Ki,kHi,k)Pi,k/k-1 (23)
Wherein, Ri,k、Zi,kThe respectively observation noise variance matrix and observation information of k moment subfilters, Hi,kFor son The structural parameters of filter observational equation.
Due to error very little,By state estimation in the present inventionChange Letter is (22) formula, contributes to the stability and convergence that promote federated filter.
(5) data fusion:The state estimation of M1 pairs of two subfilters of senior filterData fusion is carried out, is obtained most Whole state estimationAnd Square Error matrix Pk
(6) result exports:By state estimationEstimate attitude quaternion for correctingWith Random Constant Drift estimated valueObtain final k moment attitude quaternions qkWith Random Constant Drift value bk
Wherein,Estimate to obtain by the k-1 moment, think in the present invention
The content that description in the present invention is not described in detail belongs to the known technology of professional and technical personnel in the field.

Claims (10)

1. a kind of microsatellite attitude based on multisensor determines method, which is characterized in that include the following steps:
(1) MEMS gyro is utilized to obtain angular velocity measurement valueEstablish MEMS gyro error model;
(2) according to the geomagnetic fieldvector measured value for measuring acquisitionAnd the ground magnetic vector reference value B being calculatedo, establish magnetic strength Count error model;
(3) it measures and obtains sun vector measurement valueIt is derived by solar vector reference value S using sun ephemeriso, establish the sun Sensor error model;
(4) according to the microsatellite attitude kinematical equation and MEMS gyro error model established using quaternary number as attitude parameter Establish state equation;
(5) according to magnetometer error model and sun sensor error model, subfilter S1 and subfilter S2 are established respectively Observational equation;
(6) according to information distribution factor βiTo the noise variance matrix and Square Error matrix of subfilter S1 and subfilter S2 Carry out self-adjusted block;
(7) it integrates to obtain estimation attitude quaternion using MEMS gyro angular velocity measurement value
(8) state equation established in step (4) is utilized to calculate one step status predication of k momentAndEstimation it is square Error matrix Pk/k-1
(9) information distribution factor β is calculatedi, and reset the corresponding mean square error of subfilter S1, S2 to subfilter S1, S2 Matrix Pi,k/k-1
(10) subfilter S1, S2 carries out measurement update according to the observation information of magnetometer and sun sensor, calculates separately son Filter Kalman filtering gain Ki,k, state estimationAnd mean square error Pi,k
(11) to the state estimation of subfilter S1, S2Data fusion is carried out, final state estimation is obtainedAnd it is square Error matrix Pk
(12) by state estimationEstimate attitude quaternion for correctingWith Random Constant Drift estimated valueIt obtains finally Attitude quaternion qkWith Random Constant Drift value bk
2. a kind of microsatellite attitude based on multisensor according to claim 1 determines method, it is characterised in that:
The specific method of the step (1) is:
When MEMS gyro installation axle is consistent with the microsatellite axes of inertia, MEMS gyro exports to obtain microsatellite body coordinate system b Angular velocity measurement value relative to inertial coodinate system j
Establishing MEMS gyro error model is:
Wherein, ηg, ηbFor white noise, ωbjFor true angular velocity, b is Random Constant Drift.
3. a kind of microsatellite attitude based on multisensor according to claim 1 or 2 determines that method, feature exist In:
The specific method of the step (2) is:
When magnetometer installation axle is consistent with microsatellite body coordinate system, obtains geomagnetic fieldvector under satellite body coordinate system and measure ValueEstablishing magnetometer error model is:
Wherein, VBFor the measurement noise of magnetometer,For orbital coordinate system to the transition matrix of microsatellite body coordinate system, Bo For geomagnetic fieldvector reference value under orbital coordinate system.
The 12nd generation world earth magnetism reference model IGRF-12 is selected to be sweared for calculating earth's magnetic field under microsatellite current orbit position The reference value of in east northeast geographic coordinate system is measured, and is converted to orbital coordinate system, geomagnetic fieldvector under orbital coordinate system is obtained Reference value Bo
4. a kind of microsatellite attitude based on multisensor according to claim 3 determines method, it is characterised in that:
The specific method of the step (3) is:
It selects two axis digital sun sensors and determines the installation matrix of two axis digital sun sensors, two number of axle word sun are quick The guidance axis of sensor is consistent with microsatellite body coordinate system-Y-axis, obtains measurement of the solar vector under satellite body coordinate system ValueThe installation matrix of the two axis digital sun sensorFor:
Using unit solar vector reference value under sun ephemeris computation current time J2000 geocentric inertial coordinate systems, and by turning It changes under matrix conversion to orbital coordinate system, obtains solar vector reference value So
Establishing sun sensor error model is:
Wherein, VsFor the measurement noise of sun sensor.
5. a kind of microsatellite attitude based on multisensor according to claim 4 determines method, it is characterised in that:
The specific method of the step (4) is:
Establish microsatellite attitude kinematical equation:
Wherein,Δ q=[Δ q1 Δq2 Δq3]T,
Error quaternion Δ q is true attitude quaternion q and estimation attitude quaternionBetween error;
Wherein, Δ q1, Δ q2, Δ q3For the vector section of error quaternion;
Random Constant Drift error delta b is true Random Constant Drift b and estimation Random Constant Drift valueBetween error;
Sextuple state variable is set as X=[Δ q Δs b]T, using MEMS gyro error model and attitude kinematics equations, derive Obtain state equation:
Wherein, t is time parameter, matrixMatrixSquare Battle array W=[ηg ηb]T
By state equation discretization, with subscript k-1, k indicates the k-1 moment of discretization respectively, and at the k moment, T is time interval, then Have:
Xkk/k-1Xk-1k-1Wk-1,
Discretization coefficient
Discretization coefficient
Wherein, Fk-1The matrix F (t) at the k-1 moment after indicating discrete;Gk-1The matrix G (t) at the k-1 moment after indicating discrete.
6. a kind of microsatellite attitude based on multisensor according to claim 4 or 5 determines that method, feature exist In:
The step (5) is as follows:
The observational equation of subfilter S1 is obtained using magnetometer error model and state variable:
ΔBb=HbX+VB
Wherein, geomagnetic fieldvector error amount Δ BbFor magnetometer geomagnetic fieldvector measured valueWith based on IGRF-12 geomagnetic models Obtained geomagnetic fieldvector estimated valueBetween error:Matrix
Using sun sensor error model and state variable, it is derived by the observational equation of subfilter S2:
ΔSb=HsX+Vs
Wherein, solar vector error amount Δ SbThe sun vector measurement value obtained for sun sensorWith based on sun ephemeris institute The solar vector estimated value being calculatedBetween error:Matrix
7. a kind of microsatellite attitude based on multisensor according to claim 6 determines method, it is characterised in that:Institute State information distribution factor βiCalculation formula is as follows:
Wherein, S1Indicate subfilter S1 Square Error matrix P1The inverse of the sum of the absolute value of characteristic value;S2Indicate subfilter S2 Square Error matrix P2The inverse of the sum of the absolute value of characteristic value.
8. a kind of microsatellite attitude based on multisensor according to claim 7 determines method, it is characterised in that:Institute State one step status predication of k moment in step (8)AndEstimation Square Error matrix Pk/k-1Calculation formula it is as follows:
Square Error matrix P in the step (9)i,k/k-1
9. a kind of microsatellite attitude based on multisensor according to claim 7 or 8 determines that method, feature exist In:Subfilter Kalman filtering gain K in the step (10)i,k, state estimationAnd mean square error Pi,kCalculation formula It is as follows:
Pi,k=(I-Ki,kHi,k)Pi,k/k-1
Final state estimation in the step (11)And Square Error matrix PkCalculation formula it is as follows:
10. a kind of microsatellite attitude based on multisensor according to claim 9 determines method, it is characterised in that: Final attitude quaternion q in the step (12)kWith Random Constant Drift value bkCalculation formula it is as follows:
CN201711362902.XA 2017-12-18 2017-12-18 A kind of microsatellite attitude based on multisensor determines method Pending CN108279010A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711362902.XA CN108279010A (en) 2017-12-18 2017-12-18 A kind of microsatellite attitude based on multisensor determines method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711362902.XA CN108279010A (en) 2017-12-18 2017-12-18 A kind of microsatellite attitude based on multisensor determines method

Publications (1)

Publication Number Publication Date
CN108279010A true CN108279010A (en) 2018-07-13

Family

ID=62801701

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711362902.XA Pending CN108279010A (en) 2017-12-18 2017-12-18 A kind of microsatellite attitude based on multisensor determines method

Country Status (1)

Country Link
CN (1) CN108279010A (en)

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109655070A (en) * 2018-12-28 2019-04-19 清华大学 A kind of multi-mode attitude determination method of remote sensing micro-nano satellite
CN110017837A (en) * 2019-04-26 2019-07-16 沈阳航空航天大学 A kind of Combinated navigation method of the diamagnetic interference of posture
CN110260869A (en) * 2019-05-10 2019-09-20 哈尔滨工业大学 A kind of improved method reducing star sensor and gyro Federated filter calculation amount
CN110285815A (en) * 2019-05-28 2019-09-27 山东航天电子技术研究所 It is a kind of can in-orbit whole-process application micro-nano satellite multi-source information attitude determination method
CN110736468A (en) * 2019-09-23 2020-01-31 中国矿业大学 Spacecraft attitude estimation method assisted by self-adaptive kinematics model
CN111207772A (en) * 2020-01-14 2020-05-29 上海卫星工程研究所 Method for testing light path and polarity of multi-head star sensor
CN111641456A (en) * 2018-11-07 2020-09-08 长沙天仪空间科技研究院有限公司 Laser communication method based on satellite
CN111854764A (en) * 2020-07-20 2020-10-30 中国科学院微小卫星创新研究院 Spacecraft attitude determination method and system based on inter-satellite measurement information
CN112212860A (en) * 2020-08-28 2021-01-12 山东航天电子技术研究所 Distributed filtering micro-nano satellite attitude determination method with fault tolerance
CN112278329A (en) * 2020-10-30 2021-01-29 长光卫星技术有限公司 Nonlinear filtering method for remote sensing satellite attitude determination
CN113291493A (en) * 2021-05-13 2021-08-24 航天科工空间工程发展有限公司 Method and system for determining fusion attitude of satellite multi-sensor
CN113483765A (en) * 2021-05-24 2021-10-08 航天科工空间工程发展有限公司 Satellite autonomous attitude determination method

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060074558A1 (en) * 2003-11-26 2006-04-06 Williamson Walton R Fault-tolerant system, apparatus and method
CN102607564A (en) * 2012-03-09 2012-07-25 北京航空航天大学 Small satellite autonomous navigation system based on starlight/ geomagnetism integrated information and navigation method thereof
CN102914308A (en) * 2012-10-24 2013-02-06 南京航空航天大学 Anti-outlier federated filtering method based on innovation orthogonality
CN103697894A (en) * 2013-12-27 2014-04-02 南京航空航天大学 Multi-source information unequal interval federated filtering method based on filter variance matrix correction
CN104655152A (en) * 2015-02-11 2015-05-27 北京航空航天大学 Onboard distributed type POS real-time transmission alignment method based on federal filtering
CN104913781A (en) * 2015-06-04 2015-09-16 南京航空航天大学 Unequal interval federated filter method based on dynamic information distribution
CN105758401A (en) * 2016-05-14 2016-07-13 中卫物联成都科技有限公司 Integrated navigation method and equipment based on multisource information fusion

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060074558A1 (en) * 2003-11-26 2006-04-06 Williamson Walton R Fault-tolerant system, apparatus and method
CN102607564A (en) * 2012-03-09 2012-07-25 北京航空航天大学 Small satellite autonomous navigation system based on starlight/ geomagnetism integrated information and navigation method thereof
CN102914308A (en) * 2012-10-24 2013-02-06 南京航空航天大学 Anti-outlier federated filtering method based on innovation orthogonality
CN103697894A (en) * 2013-12-27 2014-04-02 南京航空航天大学 Multi-source information unequal interval federated filtering method based on filter variance matrix correction
CN104655152A (en) * 2015-02-11 2015-05-27 北京航空航天大学 Onboard distributed type POS real-time transmission alignment method based on federal filtering
CN104913781A (en) * 2015-06-04 2015-09-16 南京航空航天大学 Unequal interval federated filter method based on dynamic information distribution
CN105758401A (en) * 2016-05-14 2016-07-13 中卫物联成都科技有限公司 Integrated navigation method and equipment based on multisource information fusion

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
屠善澄: "《卫星姿态动力学与控制 2》", 30 September 1998, 宇航出版社 *
张国良等: "《组合导航原理与技术》", 31 May 2008, 西安交通大学出版社 *
梅昌明等: "基于矢量观测的微小卫星姿态确定研究", 《2015年小卫星技术交流会论文集》 *
王鹏: "基于星载敏感器的卫星自主导航及姿态确定方法研究", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 *
秦永元等: "《卡尔曼滤波与组合导航原理》", 30 June 2015, 西北工业大学出版社 *
郁丰等: "微小卫星姿态确定系统多信息融合滤波技术", 《上海交通大学学报》 *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111641456A (en) * 2018-11-07 2020-09-08 长沙天仪空间科技研究院有限公司 Laser communication method based on satellite
CN109655070B (en) * 2018-12-28 2022-05-17 清华大学 Multi-mode attitude determination method for remote sensing micro-nano satellite
CN109655070A (en) * 2018-12-28 2019-04-19 清华大学 A kind of multi-mode attitude determination method of remote sensing micro-nano satellite
CN110017837A (en) * 2019-04-26 2019-07-16 沈阳航空航天大学 A kind of Combinated navigation method of the diamagnetic interference of posture
CN110017837B (en) * 2019-04-26 2022-11-25 沈阳航空航天大学 Attitude anti-magnetic interference combined navigation method
CN110260869A (en) * 2019-05-10 2019-09-20 哈尔滨工业大学 A kind of improved method reducing star sensor and gyro Federated filter calculation amount
CN110285815A (en) * 2019-05-28 2019-09-27 山东航天电子技术研究所 It is a kind of can in-orbit whole-process application micro-nano satellite multi-source information attitude determination method
CN110285815B (en) * 2019-05-28 2023-10-20 山东航天电子技术研究所 Micro-nano satellite multi-source information attitude determination method capable of being applied in whole orbit
CN110736468A (en) * 2019-09-23 2020-01-31 中国矿业大学 Spacecraft attitude estimation method assisted by self-adaptive kinematics model
CN111207772A (en) * 2020-01-14 2020-05-29 上海卫星工程研究所 Method for testing light path and polarity of multi-head star sensor
CN111207772B (en) * 2020-01-14 2021-07-13 上海卫星工程研究所 Method for testing light path and polarity of multi-head star sensor
CN111854764A (en) * 2020-07-20 2020-10-30 中国科学院微小卫星创新研究院 Spacecraft attitude determination method and system based on inter-satellite measurement information
CN112212860A (en) * 2020-08-28 2021-01-12 山东航天电子技术研究所 Distributed filtering micro-nano satellite attitude determination method with fault tolerance
CN112212860B (en) * 2020-08-28 2023-03-03 山东航天电子技术研究所 Distributed filtering micro-nano satellite attitude determination method with fault tolerance
CN112278329A (en) * 2020-10-30 2021-01-29 长光卫星技术有限公司 Nonlinear filtering method for remote sensing satellite attitude determination
CN113291493B (en) * 2021-05-13 2022-09-23 航天科工空间工程发展有限公司 Method and system for determining fusion attitude of multiple sensors of satellite
CN113291493A (en) * 2021-05-13 2021-08-24 航天科工空间工程发展有限公司 Method and system for determining fusion attitude of satellite multi-sensor
CN113483765A (en) * 2021-05-24 2021-10-08 航天科工空间工程发展有限公司 Satellite autonomous attitude determination method

Similar Documents

Publication Publication Date Title
CN108279010A (en) A kind of microsatellite attitude based on multisensor determines method
CN102809377B (en) Aircraft inertia/pneumatic model Combinated navigation method
CN105486312B (en) A kind of star sensor and high frequency angular displacement sensor integrated attitude determination method and system
CN112697138B (en) Bionic polarization synchronous positioning and composition method based on factor graph optimization
CN100356139C (en) Miniature assembled gesture measuring system for mini-satellite
CN108051866B (en) Based on strap down inertial navigation/GPS combination subsidiary level angular movement isolation Gravimetric Method
CN104698485B (en) Integrated navigation system and air navigation aid based on BD, GPS and MEMS
CN104655152B (en) A kind of real-time Transfer Alignments of airborne distributed POS based on federated filter
CN101556155A (en) Small satellite attitude determination system and method thereof
CN107529376B (en) The method of the microsatellite non-cooperative target Relative Navigation of multimodality fusion
CN104697520B (en) Integrated gyro free strap down inertial navigation system and gps system Combinated navigation method
CN107421550A (en) A kind of earth Lagrange joint constellation autonomous orbit determination methods based on H_2O maser
CN103697894B (en) Multi-source information unequal interval federated filter method based on the correction of wave filter variance battle array
CN101949703A (en) Strapdown inertial/satellite combined navigation filtering method
CN105737858A (en) Attitude parameter calibration method and attitude parameter calibration device of airborne inertial navigation system
CN108181916A (en) The control method and device of moonlet relative attitude
CN110285815A (en) It is a kind of can in-orbit whole-process application micro-nano satellite multi-source information attitude determination method
CN110304279A (en) A kind of mass center on-orbit calibration compensation method of electric propulsion satellite
CN103712598A (en) Attitude determination system and method of small unmanned aerial vehicle
CN113291493B (en) Method and system for determining fusion attitude of multiple sensors of satellite
CN111207773B (en) Attitude unconstrained optimization solving method for bionic polarized light navigation
CN113834483A (en) Inertial/polarization/geomagnetic fault-tolerant navigation method based on observability degree
CN115900770B (en) Online correction method and system for magnetic sensor in airborne environment
CN112461262A (en) Device and method for correcting errors of three-axis magnetometer
CN112525203A (en) Spacecraft autonomous astronomical navigation method based on angle constraint auxiliary measurement

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20180713