#### Abstract

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.

#### 1. 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–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–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–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 high-rate 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.

#### 2. 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.

##### 2.1. Dynamic and Kinematic Model of Satellite

The nonlinear kinematic equations of a satellite can be written as follows [18]:
where are the angular velocity components of the body frame, 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. and are the reaction wheels control torque vector and external torque inputs, respectively. is total satellite angular momentum in the body axes and defined as follows:
The angular momentum vector of the wheels is defined as follows [19]:
is moment of inertial matrix for the wheels and is wheels speeds matrix. Moreover,
where is reaction wheels angular momentum in body axes and *C* is the orientation matrix of reaction wheels.

##### 2.2. Gyroscope Model

We formulate a simplified gyro model as described in [11] as follows: where is the true rate, is the bias drift, and is the white noise accounting for angular random walk at the gyro angle level. , 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 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] The continuous noise covariance matrix is defined as where is variance of the gyro angular random walk, is the variance of the gyro rate random walk, is variance of the gyro scale factor, is elements variance of the gyro upper misalignment, and 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 star are determined by the stellar collinearity equations [20]
where *f* is the known focal length and is the vector toward the star as
where is the right ascension, is the declination of the star in the earth-centered inertial (ECI) coordinate system that are supposed to be available in a computer-accessible catalog [20]. 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 is standard deviation of the star measurement noise () with covariance matrix . The solution of this problem has been known as the q-method [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 , so the sphere containing that point can be approximated by [20]
where denotes the star measurement and the corresponding error is approximately zero-mean Gaussian noise with covariance matrix as follows [21]:
where *d* is on the order of one and is standard deviation of CCD centroiding error which is assumed to be known.

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 () strongly.

#### 3. 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 . System matrices , , and may be of different structures for different *j*. The process noise and noise measurement vectors and are white Gaussian noises of zero mean and covariance matrices and as follows:
The system mode sequence is assumed to be a first-order Markov chain with transition probabilities matrix [22]
where denotes probability, is the discrete valued modal state at time *k*, which denotes the mode in effect during the sampling period ending at , and is the transition probability from mode to mode . is the set of all possible system models. A general structure of the IMM method is shown in Figure 3.

Steps from 1 to 4 presents a complete cycle of the IMM estimator with Kalman filters. More details can be found in [22–27].

*Step 1. *Interaction/mixing of the estimates (for ):Predicted mode probability:
Transition probability:
Mixing estimate:
Mixing covariance:

*Step 2. *Model-condition filtering (for ):Predicted state:
Predicted covariance:
Measurement residual:
Residual covariance:
Filter gain:
Updated state:
Updated covariance:

*Step 3. *Mode probability update (for ):Likelihood function:
Mode probability:

*Step 4. *Combination of estimatesOverall state estimate:
Overall covariance:

##### 3.1. 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, and are process and sensor measurement noise, respectively, andThe 9-state Kalman filter contains three star attitude errors, three gyro bias errors and three scale factor errors. The 6-state Kalman filter contains three star attitude errors and three gyro bias errors. The linearization in the operating point of the system 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].

*Step 5. *Residual update (for )

*Step 6. *IRU compensated rate:

*Step 7. *Steps 1 to 4 ().

*Step 8. *Attitude, angular rate, and IRU compensation parameters update

In Steps from 5 to 8, we have considered that We will illustrate that applying this algorithm maintains the peak estimation errors less enough at the start of maneuvers. 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.

##### 3.2. 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 () 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
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 9. *Residual update (for)
where

*Step 10. *IRU compensated rate:

*Step 12. *Attitude and IRU compensation parameters update

In this structure, we have considered that where and 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).

##### 3.3. 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 () 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).

*Step 13. *Residual update IMM_CT (for )
where

*Step 14. *IRU compensated rate:

*Step 16. *Residual update IMM_L (for )

*Step 17. *PE_1 block (for )

*Step 18. *Interaction/mixing of the estimates IMM_L ( for ):Mixing estimate:
Mixing covariance:

*Step 19. *Steps 2 and 3 by replacing with , respectively.

*Step 20. *Combination of estimatesOverall state estimate:
Overall covariance:
Overall residual:
Residual covariance:

*Step 21. *PE_2 block:
where

*Step 22. *Residual of tertiary model:

*Step 24. *Attitude, angular rate, and IRU compensation parameters update

*Step 25. *Angular rate estimation (for )

*Step 26. *Overall angular rate estimate:

#### 4. 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 and associated 3*σ* boundaries over the 540 sec during the stereoimaging maneuvers (according to Table 1) are shown in Figure 9 [20].

**(a)**

**(b)**

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 rad/sec^{3/2} and 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 sec.

##### 4.1. Performance Evaluation of the IMM_CT Attitude 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
The posterior probability of EKFs () 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.

**(a)**

**(b)**

**(c)**

The gyro errors estimations are shown in Figures 11, 12, 13, and 14.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

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].

##### 4.2. 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 different models, with different star sensor measurements covariance matrix (), in order to design the IMM_L estimator in the presence of star sensor measurement noise level changes. These models are mentioned in Table 3.

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 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 ( sec).

##### 4.3. Performance Evaluation of the EIMMAE Method

This section clearly shows the effectiveness of the proposed EIMMAE structure for attitude estimation of the stereoimagery satellite. By applying this method during stereoimaging scenario (according to Table 1), the 9-or 15-state Kalman filter should be selected at the start of maneuvers ( sec), and the IMM_L structure (6-state kalman filters with different star sensor measurements noise levels) should be selected at the end of maneuvers ( sec), and then the suitable star sensor measurement noise level should be selected using the IMM_L structure in this mode. The posterior probability of EKFs () is shown in Figure 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.

**(a)**

**(b)**

**(c)**

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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

The attitude estimation errors using the EIMMAE structure will be compared with IMM_L and IMM_CT method in Figures 18 and 19 with 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 ( sec) and reduces the attitude estimation errors at the end of maneuvers while the satellite is in “uniform motion” mode ( sec) simultaneously.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

The comparison results between EIMMAE structure with IMM_L and IMM_CT for attitude estimation at the end of maneuvers ( sec) in stereoimaging scenario are summarized in Table 4. Entries of this table contain the average of absolute attitude estimation errors around each axis at the end of three main stereoimaging maneuvers during sec which have been computed around each axis using (80) as follows: According to Table 4, it is obvious that fulfillment of the EIMMAE method improves pointing accuracy of the satellite up to 25% in comparison with applying the IMM_L method and improves pointing accuracy of the satellite up to 50% in comparison with applying the IMM_CT method.

In order to show effectiveness and advantages of the proposed EIMMAE method to estimate the attitude of the satellite at the imaging moment, we will compare the 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.

#### 5. Conclusion

In this paper, we have designed a new adaptive attitude estimator (EIMMAE) based on interacting multiple models (IMM) extended Kalman filters (EKF) for attitude determination of a stereoimagery satellite. The proposed algorithm consists of two different IMM estimator structures: “IMM_CT” structure to select the suitable Kalman filter structure to maintain the peak attitude estimation errors less enough at the start of maneuvers and “IMM_L” structure to select the suitable star sensor measurement noise level to reduce attitude estimation errors while the satellite is in “uniform motion” mode. We have shown in Table 4, Figures 18 and 19 that by applying the proposed adaptive attitude estimation method the attitude estimation requirements of stereoimagery satellite will be satisfied simultaneously. In the other hand, we have also compared the obtained results of the proposed EIMMAE method with some other previous attitude estimation methods in Table 5. This table has verified the advantages and effectiveness of the proposed EIMMAE method to provide highly accurate attitude estimation for achieving stereoimaging scenario.