CN107389069B - Ground attitude processing method based on bidirectional Kalman filtering - Google Patents

Ground attitude processing method based on bidirectional Kalman filtering Download PDF

Info

Publication number
CN107389069B
CN107389069B CN201710612514.6A CN201710612514A CN107389069B CN 107389069 B CN107389069 B CN 107389069B CN 201710612514 A CN201710612514 A CN 201710612514A CN 107389069 B CN107389069 B CN 107389069B
Authority
CN
China
Prior art keywords
kalman filtering
value
attitude
quaternion
backward
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.)
Active
Application number
CN201710612514.6A
Other languages
Chinese (zh)
Other versions
CN107389069A (en
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.)
Shanghai Aerospace Control Technology Institute
Original Assignee
Shanghai Aerospace Control Technology Institute
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 Shanghai Aerospace Control Technology Institute filed Critical Shanghai Aerospace Control Technology Institute
Priority to CN201710612514.6A priority Critical patent/CN107389069B/en
Publication of CN107389069A publication Critical patent/CN107389069A/en
Application granted granted Critical
Publication of CN107389069B publication Critical patent/CN107389069B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/20—Instruments for performing navigational calculations
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/005—Navigation; 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
    • G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
    • G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)
  • Gyroscopes (AREA)

Abstract

The invention relates to a ground attitude processing method based on bidirectional Kalman filtering, which comprises the following steps: s1, establishing a gyro measurement model; s2, establishing a star sensor measurement model; s3, establishing a gyro and star sensor combined attitude determination Kalman filtering state equation and a measurement equation; s4, performing a forward Kalman filtering recursion process, and calculating an estimated attitude quaternion and a real angular velocity estimated value of the forward Kalman filtering recursion; s5, carrying out backward Kalman filtering recursion process, and calculating the estimation attitude quaternion and the real angular velocity estimation value of backward Kalman filtering in a recursion manner; and S6, obtaining an estimated attitude quaternion and a real angular velocity estimated value of the satellite by using the covariance weight filter, and determining the attitude of the satellite on the ground. According to the method, the bidirectional Kalman filtering ground attitude processing method is adopted for the data of the star sensor and the gyroscope which are transmitted in the orbit, so that the attitude determination precision is effectively improved, and a high-precision attitude reference is provided for image navigation and registration.

Description

Ground attitude processing method based on bidirectional Kalman filtering
Technical Field
The invention relates to a ground attitude processing method, in particular to a ground attitude processing method based on bidirectional Kalman filtering, and belongs to the high-precision ground attitude processing technology.
Background
With the continuous development of earth observation satellites, higher requirements are put forward on the attitude control precision and stability of the satellites. Particularly, a high-precision optical remote sensing satellite needs to have high attitude determination precision for a satellite platform, and provides a high-precision attitude reference for image navigation and registration.
In the prior art, a common method on the satellite is to adopt a high-precision gyroscope combination and a high-precision star sensor combined Kalman filtering to determine the high-precision attitude. At present, due to the limitation of satellite computing resources, a fixed gain Kalman filtering algorithm is used on a satellite, and image navigation and registration are carried out after image data and satellite ephemeris data are downloaded to the ground.
Therefore, in order to improve the attitude reference of image navigation and registration, it is necessary to research the ground high-precision attitude determination method by utilizing the advantages of ground data processing.
Disclosure of Invention
The invention aims to provide a ground attitude processing method based on bidirectional Kalman filtering, which adopts the bidirectional Kalman filtering ground attitude processing method for data of a star sensor and a gyroscope which are transmitted under an orbit, effectively improves the attitude determination precision and provides high-precision attitude reference for image navigation and registration.
In order to achieve the above object, the present invention provides a ground attitude processing method based on bi-directional kalman filtering, comprising the following steps:
s1, establishing a gyro measurement model;
s2, establishing a star sensor measurement model;
s3, establishing a Kalman filtering state equation and a measurement equation of combined attitude determination of a gyroscope and a star sensor by using a quaternion attitude kinematics equation;
s4, performing a forward Kalman filtering recursion process, and calculating an estimated attitude quaternion and a real angular velocity estimated value of the forward Kalman filtering recursion;
s5, carrying out backward Kalman filtering recursion process, and calculating the estimation attitude quaternion and the real angular velocity estimation value of backward Kalman filtering in a recursion manner;
and S6, obtaining an estimated attitude quaternion and a real angular velocity estimated value of the satellite by using the covariance weight filter, and determining the attitude of the satellite on the ground.
In S1, the gyro measurement model is:
Figure BDA0001359843490000021
where ω represents the true angular velocity; omegasRepresenting the value of the gyro measurement value converted to the satellite body coordinate system; cbgRepresenting a conversion matrix from a gyro measurement coordinate system to a satellite body coordinate system; b represents a real value of the gyro drift; n isgRepresenting white gyroscopic measurement noise; n isdRepresenting white gyro drift slope noise.
In S2, the star sensor measurement model is:
Figure BDA0001359843490000022
wherein q isscFor measuring the attitude quaternion, the output quaternion of the star sensor is converted into a value under a satellite body coordinate system; q represents the true attitude quaternion of the satellite body relative to the inertial space; n isscRepresenting the value of the measurement noise quaternion of the star sensor converted to the satellite body coordinate system; cbsRepresenting a conversion matrix from a star sensor measurement coordinate system to a satellite body coordinate system; n issIndicating that the star sensor measures the white gaussian noise.
In S3, the kalman filter state equation and the kalman filter measurement equation are respectively:
Figure BDA0001359843490000023
wherein:
Figure BDA0001359843490000024
Figure BDA0001359843490000031
Z(t)＝[qesc1qesc2qesc3]T；
V(t)＝ns；
Figure BDA0001359843490000032
Figure BDA0001359843490000033
H＝[I3×303×3]；
D＝Cbs；
and,
Figure BDA0001359843490000034
representing gyroscope constant drift true value b and gyroscope constant drift estimated value
Figure BDA0001359843490000035
The difference between the two;
Qe＝[qe1qe2qe3]and q ise＝{1 Qe}，QeIs qeThe vector part of (a), qeRepresenting an attitude deviation quaternion for converting the estimated attitude quaternion to a true attitude quaternion;
Qesc＝[qesc1qesc2qesc3]and q isesc＝{1 Qe+Cbsns}＝{1 Qesc}，QescIs qescThe vector part of (a), qescExpressing a measurement deviation quaternion for converting the estimation attitude quaternion into a measurement attitude quaternion;
Figure BDA0001359843490000036
is an estimate of the true angular velocity.
In S4, the forward kalman filtering recursion process includes:
Figure BDA0001359843490000037
Figure BDA0001359843490000038
Figure BDA0001359843490000039
Figure BDA0001359843490000041
Figure BDA0001359843490000042
Figure BDA0001359843490000043
Figure BDA0001359843490000044
wherein:
Figure BDA0001359843490000045
predicting a k-th predicted value of forward Kalman filtering for estimating the attitude quaternion;
Figure BDA0001359843490000046
Figure BDA0001359843490000047
estimating the k-1 step value of forward Kalman filtering for the estimation attitude quaternion;
t is an updating time interval of Kalman filtering;
qesc,fkmeasuring the current time value of the deviation quaternion in the forward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,fk/k-1is QePredicting a value in the k step of forward Kalman filtering;
Δbfk/k-1predicting a predicted value of delta b in the k step of forward Kalman filtering;
Qe,fk-1is QeThe value of the step k-1 of forward Kalman filtering;
Δbfk-1is delta b at the k-1 step of forward Kalman filtering;
Qe,fkis QeA current time value of forward Kalman filtering;
Δbfkis the current time value of delta b in the forward Kalman filtering;
Kfk＝Pfk/k-1·HT(HPfk/k-1HT+Rk)-1and K isfkA gain matrix for forward Kalman filtering;
Figure BDA0001359843490000048
and P isfk/k-1A prediction error covariance matrix for forward Kalman filtering;
Figure BDA0001359843490000051
Figure BDA0001359843490000052
Rk＝Cbs·ns(Cbs·ns)Tand R iskMeasuring a noise variance matrix;
Pfk＝(I-KfkH)Pfk/k-1and P isfkError covariance matrix of forward Kalman filtering;
Figure BDA0001359843490000053
and QkIs a system noise variance matrix;
Figure BDA0001359843490000054
Figure BDA0001359843490000055
estimating the current time value of the attitude quaternion in forward Kalman filtering;
Figure BDA0001359843490000056
the gyro constant drift estimation value is the current time value of forward Kalman filtering;
Figure BDA0001359843490000057
the step k-1 value of the gyro constant drift estimation value in forward Kalman filtering is obtained;
Figure BDA0001359843490000058
the current time value of the real angular velocity estimated value in forward Kalman filtering is obtained;
k＝1，2，…。
in S5, the backward kalman filtering recursion process includes:
Figure BDA0001359843490000059
Figure BDA00013598434900000510
Figure BDA00013598434900000511
Figure BDA00013598434900000512
Figure BDA00013598434900000513
Figure BDA00013598434900000514
Figure BDA00013598434900000515
wherein,
Figure BDA0001359843490000061
predicting a k-th predicted value of backward Kalman filtering for estimating the attitude quaternion;
Figure BDA0001359843490000062
Figure BDA0001359843490000063
estimating the k +1 step value of the backward Kalman filtering for the attitude quaternion;
qesc,bkmeasuring the current time value of the deviation quaternion in backward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,bk/k+1is QePredicting the k step of backward Kalman filtering;
Δbbk/k+1predicting a value delta b in the k step of backward Kalman filtering;
Qe,bk+1is QeStep k +1 of backward Kalman filtering;
Δbbk+1is delta b at the k +1 step of backward Kalman filtering;
Qe,bkis QeThe current time value of backward Kalman filtering;
Δbbkis the current time value of delta b in backward Kalman filtering;
Kbk＝Pbk/k+1·HT(HPbk/k+1HT+Rk)-1and K isbkA gain matrix that is a backward Kalman filter;
Figure BDA0001359843490000064
and P isbk/k+1A prediction error covariance matrix for backward Kalman filtering;
Figure BDA0001359843490000065
Figure BDA0001359843490000066
Pbk＝(I-KbkH)Pbk/k+1and P isbkAn error covariance matrix of backward Kalman filtering;
Figure BDA0001359843490000067
estimating a current time value of the attitude quaternion in backward Kalman filtering;
Figure BDA0001359843490000068
the gyro constant drift estimation value is the current time value of backward Kalman filtering;
Figure BDA0001359843490000071
the gyro constant drift estimation value is the value of the (k + 1) th step of backward Kalman filtering;
Figure BDA0001359843490000072
the estimated value of the true angular velocity is the current time value of backward Kalman filtering;
k＝0，1，2，…。
in S6, the specific steps are:
Figure BDA0001359843490000073
Figure BDA0001359843490000074
Figure BDA0001359843490000075
Figure BDA0001359843490000076
Figure BDA0001359843490000077
wherein,
Figure BDA0001359843490000078
a gyro constant drift estimation value of the satellite;
Figure BDA0001359843490000079
an estimated attitude quaternion for the satellite;
Figure BDA00013598434900000710
is an estimate of the true angular velocity of the satellite.
In summary, the ground attitude processing method based on the bidirectional kalman filter provided by the invention utilizes the advantages of ground data processing, and adopts the bidirectional variable gain kalman filter ground attitude processing method to effectively improve the attitude determination precision according to the characteristics that the attitude of the star sensor and the gyroscope which are transmitted under the orbit is associated with the historical attitude and the future attitude at the current moment, thereby providing a high-precision attitude reference for image navigation and registration and effectively improving the image registration precision of the optical remote sensing satellite.
Drawings
FIG. 1 is a flow chart of a ground attitude processing method based on bi-directional Kalman filtering in the present invention;
fig. 2 is a schematic diagram of a bi-directional kalman filtering process according to the present invention.
Detailed Description
A preferred embodiment of the present invention will be described in detail below with reference to fig. 1.
As shown in fig. 1, the ground attitude processing method based on bi-directional kalman filtering provided by the present invention includes the following steps:
s1, establishing a gyro measurement model;
s2, establishing a star sensor measurement model;
s3, establishing a Kalman filtering state equation and a measurement equation of combined attitude determination of a gyroscope and a star sensor by using a quaternion attitude kinematics equation;
s4, performing a forward Kalman filtering recursion process, and calculating an estimated attitude quaternion and a real angular velocity estimated value of the forward Kalman filtering recursion;
s5, carrying out backward Kalman filtering recursion process, and calculating the estimation attitude quaternion and the real angular velocity estimation value of backward Kalman filtering in a recursion manner;
and S6, obtaining an estimated attitude quaternion and a real angular velocity estimated value of the satellite by using the covariance weight filter, and determining the attitude of the satellite on the ground.
In S1, the output value of the gyro is corrected according to the calibrated installation error, and the established gyro measurement model is:
Figure BDA0001359843490000081
where ω represents the true angular velocity; omegasRepresenting the value of the gyro measurement value converted to the satellite body coordinate system; cbgA conversion matrix from a gyro measurement coordinate system to a satellite body coordinate system is represented as an installation matrix; b represents a gyro constant drift true value; n isgRepresenting white gyroscopic measurement noise; n isdRepresenting white gyro drift slope noise.
In S2, the output quaternion of the star sensor is corrected according to the calibrated installation error, and the established star sensor measurement model is as follows:
Figure BDA0001359843490000082
wherein q isscFor measuring the attitude quaternion, the output quaternion of the star sensor is converted into a value under a satellite body coordinate system; q represents a satelliteTrue attitude quaternion of the body relative to the inertial space; n isscRepresenting the value of the measurement noise quaternion of the star sensor converted to the satellite body coordinate system; cbsRepresenting a conversion matrix from a star sensor measurement coordinate system to a satellite body coordinate system; n issIndicating that the star sensor measures the white gaussian noise.
The step S3 specifically includes the following steps:
s31, defining an attitude deviation quaternion qeComprises the following steps:
Figure BDA0001359843490000091
wherein q represents a true attitude quaternion;
Figure BDA0001359843490000092
expressing an estimated attitude quaternion, namely an estimated value of a true attitude quaternion; q. q.seRepresenting an attitude deviation quaternion for converting the estimated attitude quaternion to a true attitude quaternion;
when the satellite is maneuvering at a small angle, qeIs a small quantity, which can be approximated as:
Figure BDA0001359843490000093
wherein Q iseIs qeThe vector portion of (1);
s32, obtaining the equation of quaternion attitude kinematics:
Figure BDA0001359843490000094
wherein,
Figure BDA0001359843490000095
is an estimate of true angular velocity;
when the star sensor is adopted to compensate the constant drift (null drift) of the gyro, the estimated value of the constant drift of the gyro is set as
Figure BDA0001359843490000096
Then there are:
Figure BDA0001359843490000097
wherein, delta b is the difference between the real value and the estimated value of the gyro constant drift;
s33, ignoring high-order small quantity, and obtaining a Kalman filtering state equation as follows:
Figure BDA0001359843490000098
wherein,
Figure BDA0001359843490000101
representing an oblique symmetric matrix;
Figure BDA0001359843490000102
Figure BDA0001359843490000103
s34, defining a quaternion q of the measurement deviationescComprises the following steps:
Figure BDA0001359843490000104
i.e. from the estimated attitude quaternion
Figure BDA0001359843490000105
Rotating by a measurement deviation quaternion qescCorresponding angles are obtained to obtain a quaternion q of the measurement attitudesc(ii) a Therefore, there are:
Figure BDA0001359843490000106
ignoring the high order small quantities, we get:
Figure BDA0001359843490000107
wherein Q isescIs qescThe vector portion of (1);
s35, obtaining a Kalman filtering measurement equation:
Qesc＝Qe+Cbsns；
s36, simplifying the Kalman filtering state equation and the Kalman filtering measurement equation in expression to obtain:
Figure BDA0001359843490000108
wherein, the definition:
Figure BDA0001359843490000109
Figure BDA00013598434900001010
Z(t)＝[qesc1qesc2qesc3]T；
V(t)＝ns；
Figure BDA0001359843490000111
Figure BDA0001359843490000112
H＝[I3×303×3](ii) a I is an identity matrix;
D＝Cbs；
s37, the system noise and the measurement noise satisfy:
E[W(t)]＝0,E[W(t)TW(τ)]＝Q(t)(t-τ)；
E[V(t)]＝0,E[V(t)TV(τ)]＝R(t)(t-τ)；
wherein W (t) and V (t) are uncorrelated; q (t) is a non-negative array, R (t) is a positive array, and the numerical value is determined by the gyro measurement noise, the random noise and the star sensor measurement noise respectively; (t- τ) is a pulse function.
As shown in fig. 2, the step S4 specifically includes the following steps:
s41, discretizing the Kalman filtering state equation and the measurement equation to obtain:
Figure BDA0001359843490000113
Figure BDA0001359843490000114
Figure BDA0001359843490000115
wherein T is an updating time interval of Kalman filtering; k is 1, 2, …; wkAnd VkIs a discrete value of the noise, representing the value of the noise at the time of the k-th step,
Figure BDA0001359843490000116
Vk＝ns；
s42, the forward Kalman filtering recursion process comprises the following steps:
Figure BDA0001359843490000121
Figure BDA0001359843490000122
Figure BDA0001359843490000123
Figure BDA0001359843490000124
Figure BDA0001359843490000125
Figure BDA0001359843490000126
Figure BDA0001359843490000127
wherein:
Figure BDA0001359843490000128
predicting a k-th predicted value of forward Kalman filtering for estimating the attitude quaternion;
Figure BDA0001359843490000129
Figure BDA00013598434900001210
estimating the k-1 step value of forward Kalman filtering for the estimation attitude quaternion;
qesc,fkmeasuring the current time value of the deviation quaternion in the forward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,fk/k-1is QePredicting a value in the k step of forward Kalman filtering;
Δbfk/k-1predicting a predicted value of delta b in the k step of forward Kalman filtering;
Qe,fk-1is QeThe value of the step k-1 of forward Kalman filtering;
Δbfk-1is delta b at the k-1 step of forward Kalman filtering;
Qe,fkis QeA current time value of forward Kalman filtering;
Δbfkis the current time value of delta b in the forward Kalman filtering;
Kfk＝Pfk/k-1·HT(HPfk/k-1HT+Rk)-1and K isfkA gain matrix for forward Kalman filtering;
Figure BDA0001359843490000131
and P isfk/k-1A prediction error covariance matrix for forward Kalman filtering;
Rk＝Cbs·ns(Cbs·ns)Tand R iskMeasuring a noise variance matrix;
Pfk＝(I-KfkH)Pfk/k-1and P isfkError covariance matrix of forward Kalman filtering;
Figure BDA0001359843490000132
and QkIs a system noise variance matrix;
Figure BDA0001359843490000133
estimating the current time value of the attitude quaternion in forward Kalman filtering;
Figure BDA0001359843490000134
the gyro constant drift estimation value is the current time value of forward Kalman filtering;
Figure BDA0001359843490000135
the step k-1 value of the gyro constant drift estimation value in forward Kalman filtering is obtained;
Figure BDA0001359843490000136
and the estimated value of the real angular velocity is the current time value of the forward Kalman filtering.
As shown in fig. 2, the step S5 specifically includes the following steps:
s51, discretizing the Kalman filtering state equation and the measurement equation to obtain:
Figure BDA0001359843490000137
Figure BDA0001359843490000138
wherein T is an updating time interval of Kalman filtering; k is 0, 1, 2, …;
Figure BDA0001359843490000139
Vk＝ns；
s52, the recursion process of the backward Kalman filtering is as follows:
Figure BDA00013598434900001310
Figure BDA00013598434900001311
Figure BDA00013598434900001312
Figure BDA0001359843490000141
Figure BDA0001359843490000142
Figure BDA0001359843490000143
Figure BDA0001359843490000144
wherein,
Figure BDA0001359843490000145
to estimate attitudePredicting the k-th predicted value of the quaternion in the backward Kalman filtering;
Figure BDA0001359843490000146
Figure BDA0001359843490000147
estimating the k +1 step value of the backward Kalman filtering for the attitude quaternion;
qesc,bkmeasuring the current time value of the deviation quaternion in backward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,bk/k+1is QePredicting the k step of backward Kalman filtering;
Δbbk/k+1predicting a value delta b in the k step of backward Kalman filtering;
Qe,bk+1is QeStep k +1 of backward Kalman filtering;
Δbbk+1is delta b at the k +1 step of backward Kalman filtering;
Qe,bkis QeThe current time value of backward Kalman filtering;
Δbbkis the current time value of delta b in backward Kalman filtering;
Kbk＝Pbk/k+1·HT(HPbk/k+1HT+Rk)-1and K isbkA gain matrix that is a backward Kalman filter;
Figure BDA0001359843490000148
and P isbk/k+1A prediction error covariance matrix for backward Kalman filtering;
Pbk＝(I-KbkH)Pbk/k+1and P isbkAn error covariance matrix of backward Kalman filtering;
Figure BDA0001359843490000149
estimating a current time value of the attitude quaternion in backward Kalman filtering;
Figure BDA0001359843490000151
the gyro constant drift estimation value is the current time value of backward Kalman filtering;
Figure BDA0001359843490000152
the gyro constant drift estimation value is the value of the (k + 1) th step of backward Kalman filtering;
Figure BDA0001359843490000153
and the estimated value of the real angular velocity is the current time value of backward Kalman filtering.
As shown in fig. 2, in S6, the weighted filtering is performed on each estimation value of the forward kalman filtering state and each estimation value of the backward kalman filtering state according to the weight of the covariance, specifically:
Figure BDA0001359843490000154
Figure BDA0001359843490000155
Figure BDA0001359843490000156
Figure BDA0001359843490000157
Figure BDA0001359843490000158
wherein,
Figure BDA0001359843490000159
a gyro constant drift estimation value of the satellite;
Figure BDA00013598434900001510
an estimated attitude quaternion for the satellite;
Figure BDA00013598434900001511
is an estimate of the true angular velocity of the satellite.
In summary, the ground attitude processing method based on the bidirectional kalman filter provided by the invention utilizes the advantages of ground data processing, and adopts the bidirectional variable gain kalman filter ground attitude processing method to effectively improve the attitude determination precision according to the characteristics that the attitude of the star sensor and the gyroscope which are transmitted under the orbit is associated with the historical attitude and the future attitude at the current moment, thereby providing a high-precision attitude reference for image navigation and registration and effectively improving the image registration precision of the optical remote sensing satellite.
While the present invention has been described in detail with reference to the preferred embodiments, it should be understood that the above description should not be taken as limiting the invention. Various modifications and alterations to this invention will become apparent to those skilled in the art upon reading the foregoing description. Accordingly, the scope of the invention should be determined from the following claims.

Claims (6)

1. A ground attitude processing method based on bidirectional Kalman filtering is characterized by comprising the following steps:
s1, establishing a gyro measurement model;
s2, establishing a star sensor measurement model;
s3, establishing a Kalman filtering state equation and a measurement equation of combined attitude determination of a gyroscope and a star sensor by using a quaternion attitude kinematics equation;
s4, performing a forward Kalman filtering recursion process, and calculating an estimated attitude quaternion and a real angular velocity estimated value of the forward Kalman filtering recursion;
s5, carrying out backward Kalman filtering recursion process, and calculating the estimation attitude quaternion and the real angular velocity estimation value of backward Kalman filtering in a recursion manner;
s6, obtaining an estimated attitude quaternion and a real angular velocity estimated value of the satellite by using the covariance weight filter, and determining the attitude of the satellite on the ground;
in S3, the kalman filter state equation and the kalman filter measurement equation are respectively:
Figure FDA0002574674150000011
wherein:
Figure FDA0002574674150000012
Figure FDA0002574674150000013
Z(t)＝[qesc1qesc2qesc3]T；
V(t)＝ns；
Figure FDA0002574674150000021
Figure FDA0002574674150000022
H＝[I3×303×3]；
D＝Cbs；
and,
Figure FDA0002574674150000023
representing gyroscope constant drift true value b and gyroscope constant drift estimated value
Figure FDA0002574674150000024
The difference between the two;
Qe＝[qe1qe2qe3]and q ise＝{1 Qe}，QeIs qeThe vector part of (a), qeRepresenting an attitude deviation quaternion for converting the estimated attitude quaternion to a true attitude quaternion;
Qesc＝[qesc1qesc2qesc3]and q isesc＝{1 Qe+Cbsns}＝{1 Qesc}，QescIs qescThe vector part of (a), qescExpressing a measurement deviation quaternion for converting the estimation attitude quaternion into a measurement attitude quaternion;
Figure FDA0002574674150000025
is an estimate of true angular velocity;
ngrepresenting white gyroscopic measurement noise; n isdWhite noise representing the gyro drift slope;
nsrepresenting the measurement of Gaussian white noise by the star sensor;
Cbgrepresenting a conversion matrix from a gyro measurement coordinate system to a satellite body coordinate system;
Cbsand the transformation matrix of the star sensor measurement coordinate system to the satellite body coordinate system is represented.
2. The ground attitude processing method based on bi-directional kalman filtering according to claim 1, wherein in S1, the gyro measurement model is:
Figure FDA0002574674150000026
where ω represents the true angular velocity; omegasRepresenting the value of the gyro measurement value converted to the satellite body coordinate system; and b represents the real value of the gyro drift.
3. The ground attitude processing method based on bi-directional kalman filtering of claim 2, wherein in S2, the star sensor measurement model is:
Figure FDA0002574674150000031
wherein q isscFor measuring the attitude quaternion, the output quaternion of the star sensor is converted into a value under a satellite body coordinate system; q represents the true attitude quaternion of the satellite body relative to the inertial space; n isscAnd the measured noise quaternion of the star sensor is converted into a value under a satellite body coordinate system.
4. The bi-directional kalman filter-based ground attitude processing method according to claim 3, wherein in S4, the forward kalman filter recursion process is:
Figure FDA0002574674150000032
Figure FDA0002574674150000033
Figure FDA0002574674150000034
Figure FDA0002574674150000035
Figure FDA0002574674150000036
Figure FDA0002574674150000037
Figure FDA0002574674150000038
wherein:
Figure FDA0002574674150000039
predicting a k-th predicted value of forward Kalman filtering for estimating the attitude quaternion;
Figure FDA0002574674150000041
Figure FDA0002574674150000042
estimating the k-1 step value of forward Kalman filtering for the estimation attitude quaternion;
t is an updating time interval of Kalman filtering;
qesc,fkmeasuring the current time value of the deviation quaternion in the forward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,fk/k-1is QePredicting a value in the k step of forward Kalman filtering;
Δbfk/k-1predicting a predicted value of delta b in the k step of forward Kalman filtering;
Qe,fk-1is QeThe value of the step k-1 of forward Kalman filtering;
Δbfk-1is delta b at the k-1 step of forward Kalman filtering;
Qe,fkis QeA current time value of forward Kalman filtering;
Δbfkis the current time value of delta b in the forward Kalman filtering;
Kfk＝Pfk/k-1·HT(HPfk/k-1HT+Rk)-1and K isfkA gain matrix for forward Kalman filtering;
Figure FDA0002574674150000043
and P isfk/k-1Is a forward KaerA prediction error covariance matrix of the manfiltering;
Figure FDA0002574674150000044
Figure FDA0002574674150000045
Rk＝Cbs·ns(Cbs·ns)Tand R iskMeasuring a noise variance matrix;
Pfk＝(I-KfkH)Pfk/k-1and P isfkError covariance matrix of forward Kalman filtering;
Figure FDA0002574674150000046
and QkIs a system noise variance matrix;
Figure FDA0002574674150000051
Figure FDA0002574674150000052
estimating the current time value of the attitude quaternion in forward Kalman filtering;
Figure FDA0002574674150000053
the gyro constant drift estimation value is the current time value of forward Kalman filtering;
Figure FDA0002574674150000054
the step k-1 value of the gyro constant drift estimation value in forward Kalman filtering is obtained;
Figure FDA0002574674150000055
the current time value of the real angular velocity estimated value in forward Kalman filtering is obtained;
k＝1，2，…。
5. the bi-directional kalman filter-based ground attitude processing method according to claim 4, wherein in S5, the backward kalman filter recursion process is:
Figure FDA0002574674150000056
Figure FDA0002574674150000057
Figure FDA0002574674150000058
Figure FDA0002574674150000059
Figure FDA00025746741500000510
Figure FDA00025746741500000511
Figure FDA00025746741500000512
wherein,
Figure FDA00025746741500000513
predicting a k-th predicted value of backward Kalman filtering for estimating the attitude quaternion;
Figure FDA00025746741500000514
Figure FDA00025746741500000515
estimating the k +1 step value of the backward Kalman filtering for the attitude quaternion;
qesc,bkmeasuring the current time value of the deviation quaternion in backward Kalman filtering;
qsc,kthe measured attitude quaternion at the current moment;
Qe,bk/k+1is QePredicting the k step of backward Kalman filtering;
Δbbk/k+1predicting a value delta b in the k step of backward Kalman filtering;
Qe,bk+1is QeStep k +1 of backward Kalman filtering;
Δbbk+1is delta b at the k +1 step of backward Kalman filtering;
Qe,bkis QeThe current time value of backward Kalman filtering;
Δbbkis the current time value of delta b in backward Kalman filtering;
Kbk＝Pbk/k+1·HT(HPbk/k+1HT+Rk)-1and K isbkA gain matrix that is a backward Kalman filter;
Figure FDA0002574674150000061
and P isbk/k+1A prediction error covariance matrix for backward Kalman filtering;
Figure FDA0002574674150000062
Figure FDA0002574674150000063
Pbk＝(I-KbkH)Pbk/k+1and P isbkIs a rear directionError covariance matrix of Kalman filtering;
Figure FDA0002574674150000064
estimating a current time value of the attitude quaternion in backward Kalman filtering;
Figure FDA0002574674150000065
the gyro constant drift estimation value is the current time value of backward Kalman filtering;
Figure FDA0002574674150000066
the gyro constant drift estimation value is the value of the (k + 1) th step of backward Kalman filtering;
Figure FDA0002574674150000067
the estimated value of the true angular velocity is the current time value of backward Kalman filtering;
k＝0，1，2，…。
6. the ground attitude processing method based on bi-directional kalman filtering according to claim 5, wherein in S6, specifically:
Figure FDA0002574674150000071
Figure FDA0002574674150000072
Figure FDA0002574674150000074
Figure FDA0002574674150000075
wherein,
Figure FDA0002574674150000076
a gyro constant drift estimation value of the satellite;
Figure FDA0002574674150000077
an estimated attitude quaternion for the satellite;
Figure FDA0002574674150000078
is an estimate of the true angular velocity of the satellite.
CN201710612514.6A 2017-07-25 2017-07-25 Ground attitude processing method based on bidirectional Kalman filtering Active CN107389069B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710612514.6A CN107389069B (en) 2017-07-25 2017-07-25 Ground attitude processing method based on bidirectional Kalman filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710612514.6A CN107389069B (en) 2017-07-25 2017-07-25 Ground attitude processing method based on bidirectional Kalman filtering

Publications (2)

Publication Number Publication Date
CN107389069A CN107389069A (en) 2017-11-24
CN107389069B true CN107389069B (en) 2020-08-21

Family

ID=60337161

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710612514.6A Active CN107389069B (en) 2017-07-25 2017-07-25 Ground attitude processing method based on bidirectional Kalman filtering

Country Status (1)

Country Link
CN (1) CN107389069B (en)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110319850B (en) * 2018-03-30 2021-03-16 阿里巴巴（中国）有限公司 Method and device for acquiring zero offset of gyroscope
CN109975839B (en) * 2019-04-10 2023-04-07 华砺智行(武汉)科技有限公司 Joint filtering optimization method for vehicle satellite positioning data
CN110132287B (en) * 2019-05-05 2023-05-05 西安电子科技大学 Satellite high-precision joint attitude determination method based on extreme learning machine network compensation
CN110260869B (en) * 2019-05-10 2023-04-25 哈尔滨工业大学 Improved method for reducing calculation amount of star sensor and gyroscope combined filtering
CN113686334B (en) * 2021-07-07 2023-08-04 上海航天控制技术研究所 Method for improving on-orbit combined filtering precision of star sensor and gyroscope

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH02130417A (en) * 1988-11-10 1990-05-18 Nec Corp Inter-planet navigator
US6825806B2 (en) * 2002-06-03 2004-11-30 The Boeing Company Satellite methods and structures for improved antenna pointing and wide field-of-view attitude acquisition
CN103398713A (en) * 2013-04-26 2013-11-20 哈尔滨工程大学 Method for synchronizing measured data of star sensor/optical fiber inertial equipment
CN103852082A (en) * 2012-11-30 2014-06-11 上海航天控制工程研究所 Inter-satellite measurement and gyro attitude orbit integrated smoothing estimation method
CN103957095A (en) * 2014-05-15 2014-07-30 北京航空航天大学 Time synchronization method based on intersatellite bidirectional distance measurement
CN105486312A (en) * 2016-01-30 2016-04-13 武汉大学 Star sensor and high-frequency angular displacement sensor integrated attitude determination method and system

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7761233B2 (en) * 2006-06-30 2010-07-20 International Business Machines Corporation Apparatus and method for measuring the accurate position of moving objects in an indoor environment

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH02130417A (en) * 1988-11-10 1990-05-18 Nec Corp Inter-planet navigator
US6825806B2 (en) * 2002-06-03 2004-11-30 The Boeing Company Satellite methods and structures for improved antenna pointing and wide field-of-view attitude acquisition
CN103852082A (en) * 2012-11-30 2014-06-11 上海航天控制工程研究所 Inter-satellite measurement and gyro attitude orbit integrated smoothing estimation method
CN103398713A (en) * 2013-04-26 2013-11-20 哈尔滨工程大学 Method for synchronizing measured data of star sensor/optical fiber inertial equipment
CN103957095A (en) * 2014-05-15 2014-07-30 北京航空航天大学 Time synchronization method based on intersatellite bidirectional distance measurement
CN105486312A (en) * 2016-01-30 2016-04-13 武汉大学 Star sensor and high-frequency angular displacement sensor integrated attitude determination method and system

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Enhanced forward/backward spatial filtering method for DOA estimation of narrowband coherent sources;A. Delis 等;《IEE Proceedings - Radar, Sonar and Navigation》;19960201;第143卷;第10-16页 *
一种星敏感器/陀螺地面高精度组合定姿与精度验证方法;范城城 等;《光学学报》;20161130;第36卷(第11期);第1-13页 *
空间站组装过程姿态控制方案研究;李广兴 等;《载人航天》;20120131;第18卷(第1期);第22-29页 *

Also Published As

Publication number Publication date
CN107389069A (en) 2017-11-24

Similar Documents

Publication Publication Date Title
CN107389069B (en) Ground attitude processing method based on bidirectional Kalman filtering
CN108226980B (en) Differential GNSS and INS self-adaptive tightly-coupled navigation method based on inertial measurement unit
Mu et al. A practical INS/GPS/DVL/PS integrated navigation algorithm and its application on Autonomous Underwater Vehicle
CN101788296B (en) SINS/CNS deep integrated navigation system and realization method thereof
CN104567871B (en) A kind of quaternary number Kalman filtering Attitude estimation method based on earth magnetism gradient tensor
CN105737823B (en) A kind of GPS/SINS/CNS Combinated navigation methods based on five rank CKF
CN106772524B (en) A kind of agricultural robot integrated navigation information fusion method based on order filtering
CN107728182B (en) Flexible multi-baseline measurement method and device based on camera assistance
Stančić et al. The integration of strap-down INS and GPS based on adaptive error damping
CN112798021B (en) Inertial navigation system inter-travelling initial alignment method based on laser Doppler velocimeter
Gong et al. A modified nonlinear two-filter smoothing for high-precision airborne integrated GPS and inertial navigation
CN104635251A (en) Novel INS (inertial navigation system)/ GPS (global position system) combined position and orientation method
US20230366680A1 (en) Initialization method, device, medium and electronic equipment of integrated navigation system
CN112562077B (en) Pedestrian indoor positioning method integrating PDR and priori map
CN111156986B (en) Spectrum red shift autonomous integrated navigation method based on robust adaptive UKF
CN113984044A (en) Vehicle pose acquisition method and device based on vehicle-mounted multi-perception fusion
CN110954102A (en) Magnetometer-assisted inertial navigation system and method for robot positioning
CN112146655A (en) Elastic model design method for BeiDou/SINS tight integrated navigation system
CN110929402A (en) Probabilistic terrain estimation method based on uncertain analysis
CN103884340A (en) Information fusion navigation method for detecting fixed-point soft landing process in deep space
CN110849360A (en) Distributed relative navigation method for multi-machine cooperative formation flight
Tsukerman et al. Analytic evaluation of fine alignment for velocity aided INS
CN105988129A (en) Scalar-estimation-algorithm-based INS/GNSS combined navigation method
CN115856922A (en) Loosely-coupled land combined navigation method and device, computer equipment and medium
CN111323020B (en) Independent orbit determination method based on multi-vector observation of Mars edges and centers

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
GR01 Patent grant
GR01 Patent grant