Design of an Extended InteractingMultiple Models Adaptive Estimator for Attitude Determination of a Stereoimagery Satellite

We will design an extended interacting multiple models adaptive estimator (EIMMAE) for attitude determination of a stereoimagery satellite. This algorithm is based on interacting multiple models (IMM) extended kalman filters (EKF) using star sensor and gyroscope data. In this method, the motion of satellite during stereoimaging manoeuvres is partitioned into two different modes: “manoeuvring motion” mode and “uniform motion” mode. The proposed algorithm will select the suitable Kalman filter structure to estimate gyro errors accurately in order to maintain the peak attitude estimation error less enough at the beginning of manoeuvres while the satellite is in “manoeuvring motion” mode and then will select the suitable star sensor measurement noise level at the end of manoeuvres while the satellite is in “uniform motion” mode to reduce attitude estimation errors. It will be shown that using the proposed algorithm, the attitude estimation accuracy of stereoimagery satellite will be increased considerably. The effectiveness of the proposed algorithm will be examined and compared with the previous proposed methods through numerical simulations.


Introduction
Data collection in stereo mode is the simplest and most convenient way for 3D topographic data acquisition to produce new and revision of old inaccurate databases and maps which has been matured over the 100 years [1][2][3][4].This methodology involves identifying and measuring targets on images of an object which have been taken from disparate viewpoints.These images are then used to compute three dimensional coordinates of the locations of the object.
In [5], it has been suggested an effective stereoimaging scenario to obtain the advantages of previous proposed methods.In this scenario, state-of-the-art solution is to control pitch and roll axes of satellite simultaneously in such a way that short "revisit period" and "repeat cycle" can be obtainable according to Figure 1.Therefore, it needs highly accurate and stable pointing maneuvers to be accomplished in a few seconds that require the satellite to rotate along a relatively large angle attitude very fast.
To achieve the above-mentioned scenario, an accurate and fast attitude determination system is crucial to provide precise attitude knowledge for successful satellite operations.So, we consider two basic modes of flight for achieving high accuracy attitude determination task during each of stereoimaging maneuvers: "uniform motion" mode with constant angular velocity while the satellite reaches the end of stereoimaging maneuvers and "maneuvering motion" mode with high angular velocity at the beginning of stereoimaging maneuvers [6].Therefore, a suitable attitude determination system for the stereoimagery satellite has to meet the following requirements properly in the two above-mentioned modes [6]: maintaining the peak attitude estimation errors less enough while the satellite is in the "maneuvering motion" mode and reducing the attitude estimation errors at the end of maneuvers while imaging takes place.
Difficulties associated with high-accuracy attitude estimation techniques to satisfy the above-mentioned requirements are due to the inherent nonlinearities of satellite dynamic model [7][8][9][10], estimation of gyroscope errors, and identifying of star sensor measurement noise levels which strongly affect the performance of the attitude estimation system at the end of maneuvers [6].
In order to deal with the above-mentioned problems, several approaches have been proposed [11][12][13][14].In [15], a traditional 6-state attitude determination filter, containing three star attitude errors and three gyro bias errors, has been proposed.But, in high-rate maneuvering operating conditions (e.g., "maneuvering motion" mode), the ability of this filter to maintain attitude knowledge performance is degraded by gyro scale factor and misalignment errors.To solve this problem, a 15-state Kalman filter (containing three star attitude errors, three gyro bias errors, and nine gyro error parameters) has been proposed in [11] to achieve attitude estimation and gyro calibration in highrate maneuvering operating conditions.Drawback associated with this structure is that in "uniform motion" mode, the 15-state filter behaves no differently from a 6-state filter in providing attitude estimation and gyro bias correction.Therefore, a multirate profile attitude estimation structure has been designed in [16,17] based on star sensor and gyroscope measurements data.This structure increases the estimation convergence rate of gyro bias and misalignment for large and fast maneuvers and will maintain the peak attitude estimation errors less enough while the satellite is in the "maneuvering motion" mode.But, this structure cannot improve attitude estimation accuracy at the end of maneuvers while the satellite is in "uniform motion" mode.
To circumvent the aforementioned problems in order to maintain the peak attitude estimation errors less enough in the "maneuvering motion" mode and to reduce the attitude estimation errors in the "uniform motion" mode, in this paper, we will design a new adaptive attitude estimation structure based on interacting multiple models (IMM) using star sensor and gyroscope measurements data.The proposed structure consists of two different IMM estimator structures which are called: "IMM CT" and "IMM L".The "IMM CT" structure will be selected automatically based on the model conditional likelihood functions to estimate the gyro errors in order to maintain the peak attitude estimation errors less enough while the satellite is in "maneuvering motion" mode.The "IMM L" structure will be selected automatically to identify the suitable star sensor measurement noise level to reduce attitude estimation errors while the satellite is in "uniform motion" mode.It will be shown that by using the proposed adaptive attitude estimation method, the attitude estimation requirements of stereoimagery satellite will be satisfied simultaneously.The effectiveness of the proposed algorithm method will be examined and compared with previous proposed methods through numerical simulations.

Mathematical Model
The satellite being studied here is assumed to be rigid and equipped with a three-axis gyro, a star sensor, and four reaction wheels.

Dynamic and Kinematic
Model of Satellite.The nonlinear kinematic equations of a satellite can be written as follows [18]: where ω = ω x ω y ω z are the angular velocity components of the body frame, ω o is orbit angular velocity, and ϕ, θ, ψ are the roll, pitch, and yaw angles, respectively.Motion dynamic equations of a rigid satellite can be written as J denotes the inertia tensor of satellite.u w and τ are the reaction wheels control torque vector and external torque inputs, respectively.h = [h x h y h z ] T is total satellite angular momentum in the body axes and defined as follows: The angular momentum vector of the wheels h a is defined as follows [19]: I w is moment of inertial matrix for the wheels and ω w is wheels speeds matrix.Moreover, where [h ωx h ωy h ωz ] T is reaction wheels angular momentum in body axes and C is the orientation matrix of reaction wheels.

Gyroscope Model.
We formulate a simplified gyro model as described in [11] as follows: where ω is the true rate, b is the bias drift, and n a is the white noise accounting for angular random walk at the gyro angle level.g sf , g ma are the gyro scale factor and misalignment errors, respectively, expressed as follows: By integrating (7), gyro rate measurement of ( 6) can be rewritten as follows [11]: where The drift bias is modeled as a random walk process as where n r is the white noise accounting for the rate random walk (RRW) noise process at the rate level.The elements of the gyro misalignment matrix are modeled as a random walk process as [11] d dt The continuous noise covariance matrix is defined as where σ 2 a is variance of the gyro angular random walk, σ 2 r is the variance of the gyro rate random walk, σ 2 sf is variance of the gyro scale factor, σ 2 usf is elements variance of the gyro upper misalignment, and σ 2 lsf is elements variance of the gyro lower misalignment.
At the end of stereoimaging maneuvers, the effect of scale factor and misalignment on the estimator is minimal.So, we can formulate a simplified gyro model while the satellite is in "uniform motion" mode as follows [11]: 2.3.Star Sensor Model.The stars are assumed to be inertially fixed neglecting the effects of proper motion and velocity aberration.It is established in that the photograph image plane coordinates of the jth star are determined by the stellar collinearity equations [20] x j = − f A 11 s x j + A 12 s y j + A 13 s z j A 31 s x j + A 32 s y j + A 33 s z j , where f is the known focal length and s j is the vector toward the jth star as s x j = cos δ j cos α j , s y j = cos δ j sin α j , where α i is the right ascension, δ i is the declination of the jth star in the earth-centered inertial (ECI) coordinate system that are supposed to be available in a computeraccessible catalog [20].A i j are elements of the satellite attitude matrix that are not known.So, if the measured stars can be identified as specific cataloged stars, then the attitude matrix and associated satellite orientation angles (ϕ, θ, ψ) will be determined from the measured stars in image coordinates and identified stars in inertial coordinates.This will be accomplished using the maximum-likelihood approach minimizing the following loss function [20]: where σ star vector j is standard deviation of the jth star measurement noise (υ j ) with covariance matrix Q star vector j .The solution of this problem has been known as the qmethod [20].Therefore, in order to determine the attitude of the satellite using star vectors, the star observation can be reconstructed in unit vector form as where N is the total number of star observations and It has been shown in [9] that nearly all the probability of the errors is concentrated on a very small area about the direction of As j , so the sphere containing that point can be approximated by [20] b j = As j + υ j , υ j As j = 0, (19) where b j denotes the jth star measurement and the corresponding error υ j is approximately zero-mean Gaussian noise with covariance matrix as follows [21]: Identified stars in ECI coordinate where d is on the order of one and σ s is standard deviation of CCD centroiding error which is assumed to be known.

Star sensor measurement
In summary, the computed orientation angles using the star vectors (19) and applying the q-method will be considered as the star sensor measurements (ϕ, θ, ψ).So, the accuracy of determined attitude (ϕ, θ, ψ) depends on the number of identified stars and their position in the photograph image plane.The block diagram of star sensor model is shown in Figure 2.
According to Figure 2, the number of identified star changes during the motion of satellite relative to ECI coordinate system.This will affect the standard deviation of star sensor output measurement noise (σ star ) strongly.

Interacting Multiple Models (IMM) Attitude Estimation Baseline
This adaptive estimator structure contains a bank of Kalman filters run in parallel at every time, each based on a particular model, to obtain model-based estimates and check the status of the operation system; the overall state estimate is a kind of combination of those model-based estimates.Therefore, each cycle of IMM algorithm consists of four major steps: (i) model conditional reinitialization (interacting or mixing of the estimates), in which the input to the filter matched to a certain mode is obtained by mixing the estimates of all filters at the previous time under the assumption that this particular mode is in effect at the present time; (ii) model-conditional filtering, performed in parallel for each mode; (iii) mode probability update, based on the model conditional likelihood functions; (iv) estimate combination, which yields the overall state estimate as the probabilistically weighted sum of the updated state estimates of all filters.The probability of a mode being in effect plays a key role in determining the weights in the combination of state estimates and covariance matrices.The dynamics of the plant for designing the Kalman filters for each mode are described as follows: The subscript j denotes quantities pertaining to model m j .System matrices F j , G j , and H j may be of different structures for different j.The process noise and noise measurement vectors ξ j and ϑ j are white Gaussian noises of zero mean and covariance matrices Q j and R j as follows: The system mode sequence is assumed to be a first-order Markov chain with transition probabilities matrix [22] where P{•} denotes probability, m(k) is the discrete valued modal state at time k, which denotes the mode in effect during the sampling period ending at t k , and π i j is the transition probability from mode m i to mode m j .S = {m 1 , m 2 , . . ., m N } is the set of all possible system models.A general structure of the IMM method is shown in Figure 3.

Process noise
Sensor noise Figure 3: General structure of the IMM method [22].

IMM Attitude Estimator in "Maneuvering Motion" Mode (IMM CT).
In this section, the IMM baseline algorithm will be adopted to design an attitude estimator while the satellite is in "maneuvering motion" mode.This estimator consists of three EKFs (N = 3) with different structures (6, 9, and 15 states) which are required in different operating conditions of the satellite during each of stereoimaging maneuvers to maintain the peak attitude estimation errors less enough while the satellite starts to maneuver ("maneuvering motion" mode).In this structure, the 15-and 9-state EKF estimators will be selected at the beginning of maneuvers while the satellite is in "maneuvering motion" mode, and the 6-state EKF estimator will be selected at the end of maneuvers while the satellite placed in the "uniform motion" mode [17].This structure is depicted in Figure 4. Using (1), (7), and (10), we can form a 15-state dynamic error model for the 15-state Kalman filter as described in a compact form of first order differential equation where C is the measurement matrix, w and v are process and sensor measurement noise, respectively, and The 9-state Kalman filter contains three star attitude errors, three gyro bias errors and three scale factor errors.The 6state Kalman filter contains three star attitude errors and three gyro bias errors.The linearization in the operating point of the system (x 0 , u 0 ) converts the system (39) to the form where Based on continuous system discretization theory, the discrete time model of (41) can be written as According to Figure 4 and using design principle described in Steps from 1 to 4, the estimate of gyro bias and misalignments, actual angular velocity of the satellite, and attitude estimation are updated according to the following steps (one cycle of the IMM CT method) [17].
But, as mentioned before, this structure cannot improve the steady-state estimation accuracy while the satellite is in "uniform motion" mode (deficiency of IMM CT structure).So, in the next section we will design the "IMM L" attitude estimator to reduce the steady-state estimation errors while the satellite is in "uniform motion" mode.

IMM Attitude Estimator in "Uniform Motion" Mode (IMM L).
As we mentioned in Section 2.3, the number of identified star changes during the motion of satellite relative to ECI coordinate system which will affect the standard deviation of star sensor measurement noise (σ star ) strongly.
Here, the IMM baseline algorithm will be adopted to (21) with different star sensor measurement noise levels as a new idea to decrease attitude estimation errors while the satellite is in "uniform motion" mode.The "IMM L" estimator will be designed based on 6-state EKF structure in order to identify the suitable star sensor measurement noise level.This will improve attitude estimation accuracy of the satellite at the end of maneuvers while imaging takes place.We will augment ( 2) and ( 5) with 6-state filter structure to consider the dynamic model of satellite and reaction wheels.Thus, the proposed IMM L structure consists of 12 states as follows: where

ω y ω z h ωx h ωy h ωz Δb
x Δb y Δb z , y = Δθ x Δθ y Δθ z ω x ω y ω z h ωx h ωy h ωz , U = ḣωx ḣωy ḣωz τ x τ y τ z , (49) where C is the measurement matrix as follows: The linearization and then discretization, convert the system (47) to the form of (43).The IMM L attitude estimator structure is shown in Figure 5.The estimation of gyro bias, actual angular velocity, and also attitude estimation of the satellite is updated according to Steps from 9 to 12 the (one cycle of the IMM L method).
Step 12. Attitude and IRU compensation parameters update x(4) x (5) x( 6) In this structure, we have considered that  where and R star j , j = 1, 2, . . .are the different covariance matrices of the star sensor measurement errors.It should be mentioned that applying this structure cannot reduce the peak attitude estimation errors at the beginning of the stereoimaging maneuvers (deficiency of IMM L structure).

EIMMAE Attitude Estimator.
In this section, we will design an extended interacting multiple models adaptive estimator (EIMMAE) to maintain the peak attitude estimation errors less enough while the satellite is in "maneuvering motion" mode and to reduce attitude estimation errors while the satellite is in "uniform motion" mode simultaneously.The proposed algorithm will select the suitable Kalman filter structure while the satellite is in "maneuvering motion" mode and then will select the suitable star sensor measurement noise level while the satellite placed in "uniform motion" mode.The proposed EIMMAE attitude estimator structure is shown in Figures 6 and 7.In order to provide combination of estimates from the IMM L and IMM CT outputs, we will consider PE1 and PE2 blocks to convert the states and covariance matrix from 15-state EKF structure to states and covariance matrix of 12-state EKF structure and vice versa.
The IMM L estimator block in Figure 6 has been elaborated in Figure 7.
In "interaction/mixing of the estimates IMM L" block, we have supposed that the probability of transition between different EKF structures of IMM L estimator is zero (μ L i| j = 0 for i / = j) and the transition probability from 15 and 9 states EKF structures to IMM L models is equal to the transition probability between EKFs of the IMM CT configuration.So, the mixing states estimations stage will be achieved as follows: According to the above assumptions and considering the EIMMAE block diagram, the estimation of gyro bias, actual angular velocity of the satellite, and attitude estimation are updated according to Steps from 13 to 26 the (one cycle of the EIMMAE method).

Simulation
A simulation is done to show effectiveness and efficiency of the EIMMAE method for attitude estimation of a stereoimagery satellite with Consider a scenario for achieving the stereoimaging mission shown in Figure 1.Fulfillment of this scenario requires four main slew maneuvers in a sequence.This sequence is shown in Figure 8 and starts from initial attitude (e.g., sun pointing) to nadir pointing which is called "Initial maneuver".Afterwards, the satellite will slew to point 1 from nadir pointing, and the first picture will be taken by a snapshot after a specific constant pointing stability and accuracy are reached.The second slew is achieved from point 1 to point 2 (e.g., Nadir pointing).The second picture may be taken by a snapshot after a specific constant pointing stability and accuracy are reached.The third slew is achieved from point 2 to target attitude at point 3.The second picture will be usually taken during slew 3 by a snapshot after a specific constant pointing stability and accuracy are reached.The last slew is to back toward nadir pointing under 3-axis control.These required maneuvers are summarized in Table 1 [5].
In this scenario, a typical star sensor data has been considered based on produced data in [20].A plot of the number of identified stars and star sensor measurements noise with a value of 0.005 deg for 3σ s and associated 3σ boundaries over the 540 sec during the stereoimaging maneuvers (according to Table 1) are shown in Figure 9 [20].
According to Figure 9(b), it should be mentioned that 3σ boundaries of the star sensor measurement errors will be changed according to the changes of the number of identified stars.In the following simulations, the noise parameters for the gyro measurements in (12) are given by σ r = 3 × 10 −10 rad/sec 3/2 and σ a = 3 × 10 −6 rad/sec 1/2 .The initial bias of the gyro for each axis is given by 0.01 deg/hr.The initial covariance of the gyro for the attitude error is set to 0.28 deg 2 , and the initial covariance for the gyro drift is set to 0.22 (deg/hr) 2 .The sampling interval for the star sensor and gyroscope is considered Δt = 0.1 sec.Estimator.In this simulation, we have supposed that the satellite will achieve a number of large angle maneuvers in a few seconds (100 seconds for each maneuver) to estimate the gyro bias, misalignments, and scale factor errors accurately, by applying this method according to Steps from 5 to 8 with N=3 (15, 9, and 6-state EKFs) and

Performance Evaluation of the IMM CT Attitude
0.999 0 0.001 0.008 0.99 0.002 0.008 0.002 0.99 The posterior probability of EKFs (P i (t)) is shown in Figure 10.This figure shows that at the beginning of the maneuvers, the probability of the 15-state and 9-state EKFs are more than 6-state EKF to estimate gyro errors accurately; but at the end of maneuvers the probability of the 6-state structure is higher than the other structures while satellite is in "uniform motion" mode.The gyro errors estimations are shown in Figures 11,12,13,and 14.These figures verify the effectiveness of IMM CT method to estimate gyro errors accurately and to maintain the peak error of attitude estimation less enough while the satellite is in the "maneuvering motion" mode as mentioned in [17].

Performance Evaluation of the IMM L Method.
According to Figure 9(b), the 3σ boundaries of the star sensor measurement errors change between at least two main levels (marked as "high" and "low") around each axis which are described in Table 2. Therefore, we should consider at least N = 2 3 different models, with different star sensor measurements covariance matrix (R star j , j = 1, 2, . . .8), in order to design the IMM L It has been shown in Figure 9(b) that 3σ boundaries of the star sensor measurement errors have the same level (high or low) around all axes during each period of time; so, the model #1 or #8 should be selected if the designed attitude estimator identifies the star sensor noise level correctly, by applying the IMM L estimator with and the transition probabilities matrix π The number of the most probable model is shown in   16.This figure shows that at the beginning of the maneuvers, the probability of the 15-state and 9-state EKFs is more than IMM L structure and probability of the IMM L structure is higher than the other structures while satellite is in "uniform motion" mode at the end of maneuvers.The posterior probability of the models for star sensor measurements noise level (Table 3) at the end of maneuvers is shown in Figure 17 while the IMM L structure selected.
The attitude estimation errors using the EIMMAE structure will be compared with IMM L and IMM CT method in Figures 18 and 19 with 0.999 0 0.001 0.008 0.99 0.002 0.008 0.002 0.99 These figures verify that the EIMMAE method will provide the benefits of both IMM CT and IMM L methods simultaneously for attitude estimation during the large and fast maneuvers of stereoimaging scenario.In the other words, the designed EIMMAE estimator maintains the peak error less enough while the satellite is in the "maneuvering motion" mode (t   obtained results with the other previous attitude estimation methods.For this purpose, we have summarized 1 sigma estimation errors of different attitude estimation methods in Table 5.
This table approves that the proposed EIMMAE algorithm improves the steady-state attitude estimation accuracy of the satellite considerably around roll and pitch axes at the end of stereoimaging maneuvers.According to this table, attitude estimation error around the yaw axis associated with EIMMAE method is more than the attitude estimation error around yaw axis associated with the other methods.The reason is that we have used only one star sensor with lower accuracy than the star sensors which have been used in the other mentioned methods.Therefore, in the similar situation, by using star sensor with the same characteristics being used in the other methods, we can provide more accurate attitude knowledge by applying the EIMMAE algorithm while imaging takes place.So, fulfillment of the EIMMAE algorithm is more effective than the other previous mentioned methods to reduce the peak attitude estimation errors at the start of the maneuvers and to increase accuracy of attitude estimation at the end of stereoimaging maneuvers.

Conclusion
In this paper, we have designed a new adaptive attitude estimator (EIMMAE) based on interacting multiple models

Figure 15 .
According to Figure 15, it is certified that the designed IMM L estimator identifies the star sensor measurement noise level correctly around all axes and recognizes the correct model (model #1 or model #8) at the end of maneuvers (t ∈ [180 200], [280 300], [380 400], . . .sec).

Figure 17 :
Figure 17: Posterior probability of the models for the star sensor measurements noise level.

Table 2 :
Standard deviation of the star sensor measurement noise levels around each axis.

Table 5 :
Comparison of attitude estimation methods.