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

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.

Figure 1 

Stereoimaging Scenario.

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]: 𝜔_𝑥=𝜙̇−𝜓̇sin𝜃−𝜔_𝑜𝜔cos𝜃sin𝜓,_𝑦=𝜃̇cos𝜑+𝜓̇cos𝜃sin𝜑−𝜔_𝑜,𝜔cos𝜑cos𝜓+sin𝜑sin𝜃sin𝜓_𝑧=𝜓̇cos𝜃cos𝜑−𝜃̇sin𝜙−𝜔_𝑜−sin𝜑cos𝜓+cos𝜑sin𝜃sin𝜓, (1) 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 𝐽𝜔̇=𝜏+𝜔×ℎ−𝑢_𝑤 (2)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: ℎ=𝐽𝜔+𝐶ℎ_𝑎. (3) The angular momentum vector of the wheels ℎ_𝑎 is defined as follows [19]: ℎ_𝑎=𝐼_𝑤𝐶^𝑇𝜔+𝐼_𝑤𝜔_𝑤 (4)𝐼_𝑤 is moment of inertial matrix for the wheels and 𝜔_𝑤 is wheels speeds matrix. Moreover, 𝑢_𝑤=ℎ̇_𝜔𝑥ℎ̇_𝜔𝑦ℎ̇_𝜔𝑧^𝑇=𝐶ℎ̇_𝑎, (5) 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: 𝜔_𝑔=𝜔+𝑏+𝑔_sf+𝑔_ma+𝑛_𝑎, (6) where 𝜔 is the true rate, 𝑏 is the bias drift, and 𝑛_𝑎 is the white noise accounting for angular random walk at the gyro angle level. 𝑔_sf, 𝑔_ma are the gyro scale factor and misalignment errors, respectively, expressed as follows: 𝑔_sf=diag𝑘_sf𝑥𝑘_sf𝑦𝑘_sf𝑧𝑔⋅𝜔,_ma𝑥𝑔_ma𝑦𝑔_ma𝑧=0𝑘_𝑥𝑦𝑘_𝑥𝑧𝑘_𝑦𝑥0𝑘_𝑦𝑧𝑘_𝑧𝑥𝑘_𝑧𝑦0⋅𝜔_𝑥𝜔_𝑦𝜔_𝑧. (7) By integrating (7), gyro rate measurement of (6) can be rewritten as follows [11]: 𝜔_𝑔=(𝐼+𝐾)𝜔+𝑏+𝑛_𝑎, (8) where 𝑘𝐾=^sf𝑥𝑘_𝑥𝑦𝑘_𝑥𝑧𝑘_𝑦𝑥𝑘_sf𝑦𝑘_𝑦𝑧𝑘_𝑧𝑥𝑘_𝑧𝑦𝑘_sf𝑧. (9) The drift bias is modeled as a random walk process as 𝑑𝑏𝑑𝑡^𝑥𝑏_𝑦𝑏_𝑧=𝑛_𝑟𝑥𝑛_𝑟𝑦𝑛_𝑟𝑧, (10) 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] 𝑑𝑘𝑑𝑡^sf𝑥𝑘_sf𝑦𝑘_sf𝑧=𝑛_sf𝑥𝑛_sf𝑦𝑛_sf𝑧,𝑑𝑘𝑑𝑡^𝑥𝑦𝑘_𝑥𝑧𝑘_𝑦𝑥𝑘_𝑦𝑧𝑘_𝑧𝑥𝑘_𝑧𝑦=𝑛_𝑔𝑥𝑦𝑛_𝑔𝑥𝑧𝑛_𝑔𝑦𝑥𝑛_𝑔𝑦𝑧𝑛_𝑔𝑧𝑥𝑛_𝑔𝑧𝑦. (11) The continuous noise covariance matrix is defined as 𝑄^Gyro=𝜎²_𝑎𝐼_3×30_3×30_3×30_3×30_3×30_3×3𝜎²_𝑟𝐼_3×30_3×30_3×30_3×30_3×30_3×3𝜎²_sf𝐼_3×30_3×30_3×30_3×30_3×30_3×3𝜎²_usf𝐼_3×30_3×30_3×30_3×30_3×30_3×3𝜎²_lsf𝐼_3×3, (12) where 𝜎²_𝑎 is variance of the gyro angular random walk, 𝜎²_𝑟 is the variance of the gyro rate random walk, 𝜎²_sf is variance of the gyro scale factor, 𝜎²_usf is elements variance of the gyro upper misalignment, and 𝜎²_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]: 𝜔_𝑔=𝜔+𝑏+𝑛_𝑎. (13)

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 𝑗th star are determined by the stellar collinearity equations [20] 𝑥_𝑗𝐴=−𝑓¹¹𝑠_𝑥𝑗+𝐴₁₂𝑠_𝑦𝑗+𝐴₁₃𝑠_𝑧𝑗𝐴₃₁𝑠_𝑥𝑗+𝐴₃₂𝑠_𝑦𝑗+𝐴₃₃𝑠_𝑧𝑗,𝑦_𝑗𝐴=−𝑓²¹𝑠_𝑥𝑗+𝐴₂₂𝑠_𝑦𝑗+𝐴₂₃𝑠_𝑧𝑗𝐴₃₁𝑠_𝑥𝑗+𝐴₃₂𝑠_𝑦𝑗+𝐴₃₃𝑠_𝑧𝑗, (14) where f is the known focal length and 𝑠_𝑗 is the vector toward the 𝑗th star as𝑠_𝑥𝑗=cos𝛿_𝑗cos𝛼_𝑗,𝑠_𝑦𝑗=cos𝛿_𝑗sin𝛼_𝑗,𝑠_𝑧𝑗=sin𝛿_𝑗, (15) where 𝛼_𝑖 is the right ascension, 𝛿_𝑖 is the declination of the 𝑗th 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]: 𝑗𝐴=^𝑁_𝑗=1𝜎⁻²_{star_vector_𝑗}𝑏̃^𝑇_𝑗𝐴𝑠_𝑗, (16) where 𝜎_{star_vector_𝑗} is standard deviation of the 𝑗th star measurement noise (𝜐_𝑗) with covariance matrix 𝑄^star_vector_𝑗. 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 𝑏_𝑗=𝐴𝑠_𝑗,𝑗=1,2,…,𝑁, (17) where N is the total number of star observations and 𝑏_𝑗=1𝑓²+𝑥²_𝑗+𝑦²_𝑗−𝑥_𝑗−𝑦_𝑗𝑓,𝑠_𝑗=𝑠_𝑥𝑗𝑠_𝑦𝑗𝑠_𝑧𝑗^𝑇. (18) 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] 𝑏̃_𝑗=𝐴𝑠_𝑗+𝜐_𝑗,𝜐_𝑗𝐴𝑠_𝑗=0, (19) where 𝑏̃_𝑗 denotes the 𝑗th star measurement and the corresponding error 𝜐_𝑗 is approximately zero-mean Gaussian noise with covariance matrix as follows [21]: 𝑄^star_vector_𝑗=𝜎²_𝑠1+𝑑𝑥²_𝑖+𝑦²_𝑖1+𝑑𝑥²_𝑖²𝑑𝑥_𝑖𝑦_𝑖²𝑑𝑥_𝑖𝑦_𝑖²1+𝑑𝑦²_𝑖², (20) 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.

Figure 2 

Star sensor model.

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.

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: 𝑥(𝑘+1)=𝐹_𝑗(𝑘)⋅𝑥(𝑘)+𝐺_𝑗(𝑘)⋅𝑈(𝑘)+𝜉_𝑗(𝑘),𝑧(𝑘)=𝐻_𝑗(𝑘)⋅𝑥(𝑘)+𝜗_𝑗(𝑘). (21) 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: 𝐸𝜉_𝑗(𝑘)𝜉^𝑇_𝑗(𝑘)=𝑄_𝑗,𝐸𝜗_𝑗(𝑘)𝜗^𝑇_𝑗(𝑘)=𝑅_𝑗. (22) The system mode sequence is assumed to be a first-order Markov chain with transition probabilities matrix [22] 𝑃𝑚_𝑗(𝑘+1)∣𝑚_𝑖(𝑘)=𝜋_𝑖𝑗(𝑘)∀𝑚_𝑖,𝑚_𝑗∈𝑆,_𝑗𝜋_𝑖𝑗(𝑘)=1,𝑖=1,2,…,𝑁, (23) 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.

Figure 3 

General structure of the IMM method [22].

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 𝑗=1,2,…,𝑁):Predicted mode probability: 𝜇_𝑗(𝑘+1∣𝑘)=^𝑁_𝑖=1𝜋_𝑖𝑗𝜇_𝑖(𝑘). (24)Transition probability: 𝜇_𝑖∣𝑗=^𝑁_𝑖=1𝜋_𝑖𝑗𝜇_𝑖(𝑘)𝜇_𝑗(𝑘+1∣𝑘). (25)Mixing estimate: 𝑥⁰_𝑗(𝑘∣𝑘)=^𝑁_𝑖=1𝑥_𝑖(𝑘∣𝑘)𝜇_𝑖∣𝑗(𝑘). (26)Mixing covariance: 𝑃⁰_𝑗(𝑘∣𝑘)=^𝑁_𝑖=1𝑃_𝑗(𝑘∣𝑘)+𝑥⁰_𝑗(𝑘∣𝑘)−𝑥_𝑖(×𝑥𝑘∣𝑘)_0𝑗(𝑘∣𝑘)−𝑥^𝑇_𝑖(𝑘∣𝑘)𝜇_𝑖∣𝑗(𝑘). (27)

Step 2. Model-condition filtering (for 𝑗=1,2,…,𝑁):Predicted state: 𝑥_𝑗(𝑘+1∣𝑘+1)=𝐹_𝑗(𝑘)𝑥⁰_𝑗(𝑘∣𝑘)+𝐺_𝑗(𝑘)𝑈(𝑘)+𝜉_𝑗. (28)Predicted covariance: 𝑃_𝑗(𝑘+1∣𝑘)=𝐹_𝑗(𝑘)𝑃⁰_𝑗(𝑘∣𝑘)𝐹^𝑇_𝑗(𝑘)+𝑄_𝑗(𝑘). (29)Measurement residual: 𝑟_𝑗=𝑧(𝑘+1)−𝐻_𝑗(𝑘+1)𝑥_𝑗(𝑘+1∣𝑘). (30)Residual covariance: 𝑆_𝑗=𝐻_𝑗(𝑘+1)𝑃_𝑗(𝑘+1∣𝑘)𝐻^𝑇_𝑗(𝑘+1)+𝑅_𝑗(𝑘+1). (31)Filter gain: 𝐾_𝑗=𝑃_𝑗(𝑘+1∣𝑘)𝐻^𝑇_𝑗(𝑘+1)𝑆_𝑗⁻¹. (32)Updated state: 𝑥_𝑗(𝑘+1∣𝑘+1)=𝑥_𝑗(𝑘+1∣𝑘)+𝐾_𝑗𝑟_𝑗. (33)Updated covariance: 𝑃_𝑗(𝑘+1∣𝑘+1)=𝑃_𝑗(𝑘+1∣𝑘)+𝐾_𝑗𝑆_𝑗𝐾^𝑇_𝑗. (34)

Step 3. Mode probability update (for 𝑗=1,2,…,𝑁):Likelihood function: 𝐿_𝑗()=1𝑘+12𝜋𝑆_𝑗−1exp2𝑟^𝑇_𝑗𝑆_𝑗⁻¹(𝑘+1)𝑟_𝑗. (35)Mode probability: 𝜇_𝑗()=𝜇𝑘+1^𝑗(𝑘+1∣𝑘)𝐿_𝑗(𝑘+1)∑^𝑁_𝑖=1𝜇_𝑗(𝑘+1∣𝑘)𝐿_𝑗(𝑘+1). (36)

Step 4. Combination of estimatesOverall state estimate: 𝑥(𝑘+1∣𝑘+1)=^𝑁_𝑖=1𝑥_𝑖(𝑘+1∣𝑘+1)𝜇_𝑖(𝑘+1). (37)Overall covariance: 𝑃()=𝑘+1∣𝑘+1_𝑁𝑖=1{𝑃()+𝑥()−𝑥𝑘+1∣𝑘+1𝑘+1∣𝑘+1_𝑖(×𝑥()−𝑥𝑘+1∣𝑘+1)𝑘+1∣𝑘+1_𝑖(𝑘+1∣𝑘+1)^𝑇×𝜇_𝑖(𝑘+1). (38)

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.

Figure 4 

IMM_CT attitude estimator structure.

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 𝑥̇=𝑓_15×1(𝑥,𝑢,𝑡)+𝑤,𝑦=𝐶𝑥+𝑣, (39) where C is the measurement matrix, 𝑤 and 𝑣 are process and sensor measurement noise, respectively, and𝑥=Δ𝜃_𝑥Δ𝜃_𝑦Δ𝜃_𝑧Δ𝑏_𝑥Δ𝑏_𝑦Δ𝑏_𝑧Δ𝑘_sf𝑥Δ𝑘_sf𝑦Δ𝑘_sf𝑧Δ𝑘_𝑥𝑦Δ𝑘_𝑥𝑧Δ𝑘_𝑦𝑥Δ𝑘_𝑦𝑧Δ𝑘_𝑧𝑥Δ𝑘_𝑧𝑦^𝑇,𝑦=Δ𝜃_𝑥Δ𝜃_𝑦Δ𝜃_𝑧^𝑇,𝑢=𝜔_𝑥𝜔_𝑦𝜔_𝑧^𝑇. (40)The 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 𝑥̇=𝐴𝑥+𝐵𝑢+𝑤,𝑦=𝐶𝑥+𝑣, (41) where 𝐶=𝐼_3×30_3×123×15𝐼,𝐵=^3×30_12×3. (42) Based on continuous system discretization theory, the discrete time model of (41) can be written as 𝑥()(𝑘)(𝑘)(𝑘)(𝑘)(𝑘),𝑦(𝑘)(𝑘)(𝑘)(𝑘).𝑘+1=𝐹⋅𝑥+𝐺⋅𝑈+𝜉=𝐻⋅𝑥+𝜗 (43) 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 𝑗=1,2,3) 𝑥́_𝑗(4∶15)=𝑥_𝑗()−𝑥(),́́́4∶154∶15𝑥(1)𝑥(2)𝑥(3)_𝑛=10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛×𝜑𝜃𝜓₀−𝜑𝜃𝜓_𝑛,𝜑𝜃𝜓₀=𝜑𝜃𝜓_𝑛−1+10−sin(𝜃_𝑛−1)0cos(𝜑_𝑛−1)cos𝜃_𝑛−1sin(𝜑_𝑛−1)0−sin𝜑_𝑛−1cos𝜑_𝑛−1cos𝜃_𝑛−1⁻¹×𝑥_𝑗(1)𝑥_𝑗(2)𝑥_𝑗(3)_𝑛. (44)

Step 6. IRU compensated rate: 𝜔=𝐼−𝐾𝜔_𝑔−𝑏̂. (45)

Step 7. Steps 1 to 4 (𝑁=3).

Step 8. Attitude, angular rate, and IRU compensation parameters update 𝜑𝜃𝜓_new=𝜑𝜃𝜓_𝑛+10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛⁻¹⋅𝑥(1)𝑥(2)𝑥(3)_𝑛+1,𝑏̂(𝑛+1)=𝑏̂(𝑛)+𝑥_4∶6(),𝑔𝑛+1_𝑖𝑗(𝑛+1)=𝑔_𝑖𝑗(𝑛)+𝑥_7∶15(),𝜔𝑛+1_𝑛=𝜔_𝑛−𝑥(4)𝑥(5)𝑥(6)_𝑛+1−𝜔_𝑥000𝜔_𝑦000𝜔_𝑧⋅𝑥(7)𝑥(8)𝑥(9)_𝑛+1−𝜔_𝑦000𝜔_𝑧000𝜔_𝑥⋅𝑥()𝑥()𝑥()101112_𝑛+1−𝜔_𝑧000𝜔_𝑥000𝜔_𝑦⋅𝑥()𝑥()𝑥()131415_𝑛+1. (46)

In Steps from 5 to 8, we have considered that 𝑅_𝑗=max𝜎^star2𝑄,𝑗=1,2,3,₁=𝑄^Gyro,𝑄₂=𝜎²_𝑎𝐼_3×30_3×30_3×30_3×3𝜎²_𝑟𝐼_3×30_3×30_3×30_3×3𝜎²_sf𝐼_3×3,𝑄₃=𝜎²_𝑎𝐼_3×30_3×30_3×3𝜎²_𝑟𝐼_3×3. (47) 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 (𝜎^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: 𝑥̇=𝑓_12×1(𝑥,𝑢,𝑡)+𝑤,𝑦=𝐶𝑥+𝑣, (48) where 𝑥=Δ𝜃_𝑥Δ𝜃_𝑦Δ𝜃_𝑧𝜔_𝑥𝜔_𝑦𝜔_𝑧ℎ_𝜔𝑥ℎ_𝜔𝑦ℎ_𝜔𝑧Δ𝑏_𝑥Δ𝑏_𝑦Δ𝑏_𝑧,𝑦=Δ𝜃_𝑥Δ𝜃_𝑦Δ𝜃_𝑧𝜔_𝑥𝜔_𝑦𝜔_𝑧ℎ_𝜔𝑥ℎ_𝜔𝑦ℎ_𝜔𝑧,𝑈=ℎ̇_𝜔𝑥ℎ̇_𝜔𝑦ℎ̇_𝜔𝑧𝜏_𝑥𝜏_𝑦𝜏_𝑧, (49) where C is the measurement matrix as follows: 𝐶=𝐼_9×90_3×3𝐼_3×30_3×39×12. (50) 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).

Figure 5 

IMM_L attitude estimator structure.

Step 9. Residual update (for𝑗=1,2,…,8) 𝑥_𝑗́(4∶9)=𝑥_𝑗(),𝑥́4∶9_𝑗(10∶12)=𝑥_𝑗()−𝑥(),́́́10∶1210∶12𝑥(1)𝑥(2)𝑥(3)_𝑛=10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛×𝜑𝜃𝜓₀−𝜑𝜃𝜓_𝑛, (51) where 𝜑𝜃𝜓₀=𝜑𝜃𝜓_𝑛−1+10−sin(𝜃_𝑛−1)0cos(𝜑_𝑛−1)cos𝜃_𝑛−1sin(𝜑_𝑛−1)0−sin𝜑_𝑛−1cos𝜑_𝑛−1cos(𝜃_𝑛−1)⁻¹×𝑥_𝑗(1)𝑥_𝑗(2)𝑥_𝑗(3)_𝑛. (52)

Step 10. IRU compensated rate: 𝜔=𝜔_𝑔−𝑏̂. (53)

Step 11. Steps 1 to 4 (𝑁=8)

Step 12. Attitude and IRU compensation parameters update 𝜑𝜃𝜓_new=𝜑𝜃𝜓_𝑛+10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛⁻¹⋅𝑥(1)𝑥(2)𝑥(3)_𝑛+1,𝜔_𝑛=𝑥(4)𝑥(5)𝑥(6)_𝑛+1,𝑏̂(𝑛+1)=𝑏̂(𝑛)+𝑥_10∶12(𝑛+1). (54)

In this structure, we have considered that 𝑅_𝑗=𝑅^star_𝑗0_3×30_3×30_3×3𝜎²_𝑟𝐼_3×30_3×30_3×30_3×3𝜎²_wheel𝐼_3×39×9,𝑄_𝑗=𝑄=𝑄_{Satellite_Kinematic3×3}0_3×30_3×30_3×30_3×3𝑄^{Satellite_Dynamic}_3×30_3×30_3×30_3×30_3×3𝑄^ℎ𝑤_3×30_3×30_15×30_15×30_3×3𝜎²_𝑎𝐼_3×312×12, (55) where 𝑄^ℎ𝑤_3×3=𝜎²_ℎ𝑤𝐼_3×3,𝑄^{Satellite_Dynamic}_3×3=𝜎²_SD𝐼_3×3,𝑄^{Satellite_Kinematic}_3×3=𝜎²_SK𝐼_3×3, (56) and 𝑅^star_𝑗,𝑗=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).

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.

Figure 6 

EIMMAE structure.

Figure 7 

IMM_L estimator block.

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 (𝜇_{𝐿𝑖∣𝑗}=0for𝑖≠𝑗) 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: 𝑥⁰_𝑗(𝑘∣𝑘)=𝑥́_𝑗(𝑘∣𝑘)𝜇₃(𝑘)+𝑥₁(𝑘∣𝑘)𝜇_1∣3(𝑘)+𝑥₂(𝑘∣𝑘)𝜇_2∣3(𝑘). (57) 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 𝑗=1,2) 𝑥́_𝑗(4∶15)=𝑥_𝑗()−𝑥(),́́́4∶54∶15𝑥(1)𝑥(2)𝑥(3)_𝑛=10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛×𝜑𝜃𝜓₀−𝜑𝜃𝜓_𝑛, (58) where 𝜑𝜃𝜓₀=𝜑𝜃𝜓_𝑛−1+10−sin(𝜃_𝑛−1)0cos(𝜑_𝑛−1)cos𝜃_𝑛−1sin(𝜑_𝑛−1)0−sin𝜑_𝑛−1cos𝜑_𝑛−1cos(𝜃_𝑛−1)⁻¹×𝑥_𝑗(1)𝑥_𝑗(2)𝑥_𝑗(3)_𝑛. (59)

Step 14. IRU compensated rate: 𝜔=𝐼−𝐾𝜔_𝑔−𝑏̂. (60)

Step 15. Steps 1 and 2 (𝑁=2).

Step 16. Residual update IMM_L (for 𝑗=1,2,…,8) 𝑥́_𝑙𝑗(4∶9)=𝑥_𝑙𝑗(),𝑥́4∶9_𝑙𝑗(10∶12)=𝑥_𝑙𝑗(10∶12)−𝑥_𝑙(),́́́10∶12𝑥(1)𝑥(2)𝑥(3)_𝑙𝑗=10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛×𝜑𝜃𝜓₀−𝜑𝜃𝜓_𝑛,𝜑𝜃𝜓₀=𝜑𝜃𝜓_𝑛−1+10−sin(𝜃_𝑛−1)0cos(𝜑_𝑛−1)cos𝜃_𝑛−1sin(𝜑_𝑛−1)0−sin𝜑_𝑛−1cos𝜑_𝑛−1cos(𝜃_𝑛−1)⁻¹×𝑥_𝑙𝑗(1)𝑥_𝑙𝑗(2)𝑥_𝑙𝑗(3). (61)

Step 17. PE_1 block (for 𝑗=1,2) 𝑥_𝑗(1∶3)=𝑥_𝑗́(),𝑥1∶3_𝑗(4∶6)=𝜔_𝑛−1−𝜔_𝑥000𝜔_𝑦000𝜔_𝑧𝑛−1⋅𝑥_𝑗(7)𝑥_𝑗(8)𝑥_𝑗(9)−𝜔_𝑦000𝜔_𝑧000𝜔_𝑥𝑛−1⋅𝑥_𝑗()𝑥10_𝑗()𝑥11_𝑗()−𝜔12_𝑧000𝜔_𝑥000𝜔_𝑦⋅𝑥_𝑗()𝑥13_𝑗()𝑥́14_𝑗(),𝑥15_𝑗(7∶9)=wheelsmomentum,𝑥_𝑗(10∶12)=𝑥́_𝑗(4∶6). (62)

Step 18. Interaction/mixing of the estimates IMM_L ( for 𝑗=1,2,…,8):Mixing estimate: 𝑥⁰_𝑙𝑗(𝑘∣𝑘)=𝑥́_𝑙𝑗(𝑘∣𝑘)𝜇₃(𝑘)+𝑥₁(𝑘∣𝑘)𝜇_1∣3(𝑘)+𝑥₂(𝑘∣𝑘)𝜇_2∣3(𝑘). (63)Mixing covariance: 𝑃⁰_𝑙𝑗(𝑘∣𝑘)=²_𝑖=1𝑃_𝑗(𝑘∣𝑘)+𝑥⁰_𝑙𝑗(𝑘∣𝑘)−𝑥_𝑖(×𝑥𝑘∣𝑘)_0𝑙𝑗(𝑘∣𝑘)−𝑥^𝑇_𝑖(𝑘∣𝑘)𝜇_𝑖∣3(𝑘)+𝑃_𝑙𝑗(𝑘∣𝑘)+𝑥⁰_𝑙𝑗(𝑘∣𝑘)−𝑥́_𝑙𝑗(×𝑥𝑘∣𝑘)_0𝑙𝑗(𝑘∣𝑘)−𝑥́^𝑇_𝑙𝑗(𝑘∣𝑘)𝜇_3∣3(𝑘). (64)

Step 19. Steps 2 and 3 by replacing 𝑥,𝐹_𝑗,𝑥⁰_𝑗,𝐺_𝑗 with 𝑥_𝑙,𝐹_𝑙𝑗,𝑥⁰_𝑙𝑗,𝐺_𝑙𝑗, respectively.

Step 20. Combination of estimatesOverall state estimate: 𝑥_𝑙(𝑘+1∣𝑘+1)=⁸_𝑖=1𝑥_𝑙𝑖(𝑘+1∣𝑘+1)𝜇_𝑙𝑖(𝑘+1). (65)Overall covariance: 𝑃_𝑙()=𝑘+1∣𝑘+1_8𝑖=1𝑃_𝑙𝑖()+𝑥𝑘+1∣𝑘+1_𝑙(𝑘+1∣𝑘+1)−𝑥_𝑙𝑖(×𝑥()−𝑥𝑘+1∣𝑘+1)𝑘+1∣𝑘+1_𝑙𝑖(𝑘+1∣𝑘+1)^𝑇×𝜇_𝑙𝑖(𝑘+1). (66)Overall residual: 𝑟̂=𝑧(𝑘+1)−𝐻_𝑙(𝑘+1)𝑥_𝑙(𝑘+1∣𝑘). (67)Residual covariance: 𝑆_𝑙()=𝑘+1∣𝑘+1_8𝑖=1𝑆_𝑙𝑖()+𝑟̂𝑘+1∣𝑘+1_𝑙(𝑘+1∣𝑘+1)−𝑟̂_𝑙𝑖(×𝑟̂𝑘+1∣𝑘+1)_𝑙(𝑘+1∣𝑘+1)−𝑟̂_𝑙𝑖(𝑘+1∣𝑘+1)^𝑇×𝜇_𝑙𝑖(𝑘+1). (68)

Step 21. PE_2 block: 𝑥₃=𝐴⋅𝑥_𝑙(1∶3), (69) where 𝐼𝐴=^3×30_3×60_3×30_3×30_3×6𝐼_3×30_9×30_9×60_9×3. (70)

Step 22. Residual of tertiary model: 𝑟₃̂=𝑧(𝑘+1)−𝐻(𝑘+1)𝑥_𝑙(),𝑆𝑘+1∣𝑘₃=𝐴⋅𝑆_𝑙⋅𝐴^𝑇. (71)

Step 23. Steps 3 and 4

Step 24. Attitude, angular rate, and IRU compensation parameters update 𝜑𝜃𝜓_new=𝜑𝜃𝜓_𝑛+10−sin𝜃_𝑛0cos𝜑_𝑛cos𝜃_𝑛sin𝜑_𝑛0−sin𝜑_𝑛cos𝜑_𝑛cos𝜃_𝑛⁻¹⋅𝑥(1)𝑥(2)𝑥(3)_𝑛+1𝑏̂(𝑛+1)=𝑏̂(𝑛)+𝑥_4∶6(),𝑔𝑛+1_𝑖𝑗(𝑛+1)=𝑔_𝑖𝑗(𝑛)+𝑥_7∶15(𝑛+1). (72)

Step 25. Angular rate estimation (for 𝑗=1∶2) 𝜔_𝑗=𝜔_𝑗−𝑥_𝑗(4)𝑥_𝑗(5)𝑥_𝑗(6)_𝑘+1−𝜔_𝑥000𝜔_𝑦000𝜔_𝑧⋅𝑥_𝑗(7)𝑥_𝑗(8)𝑥_𝑗(9)_𝑘+1−𝜔_𝑦000𝜔_𝑧000𝜔_𝑥⋅𝑥_𝑗()𝑥10_𝑗()𝑥11_𝑗(12)_𝑘+1−𝜔_𝑧000𝜔_𝑥000𝜔_𝑦⋅𝑥_𝑗()𝑥13_𝑗()𝑥14_𝑗(15)_𝑘+1. (73)

Step 26. Overall angular rate estimate: 𝜔=𝜇₁𝜔₁+𝜇₂𝜔₂+𝜇₃𝜔_𝑙. (74)

4. Simulation

A simulation is done to show effectiveness and efficiency of the EIMMAE method for attitude estimation of a stereoimagery satellite with 𝐽=389.99−3.28−11.57−3.28391.83−7.42−11.57−7.42176.58Kgm². (75) 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].

Table 1 

Mission sequence for stereoimaging.

Figure 8 

Stereoimaging Scenario.

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

(a)

(b)

Figure 9 

(a) Number of identified stars (b) Star sensor measurement errors and boundaries.

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 𝜎_𝑟=3×10⁻¹⁰ rad/sec^3/2 and 𝜎_𝑎=3×10⁻⁶ 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², and the initial covariance for the gyro drift is set to 0.22 (deg/hr)². The sampling interval for the star sensor and gyroscope is considered Δ𝑡=0.1 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 𝜋_𝑖𝑗=0.99900.0010.0080.990.0020.0080.0020.99. (76) 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)

Figure 10 

Posterior probability of the EKFs structures in IMM_CT method.

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

(a)

(b)

(c)

Figure 11 

Gyro bias estimation.

(a)

(b)

(c)

Figure 12 

Scale factor estimation.

(a)

(b)

(c)

Figure 13 

Upper misalignment error estimation.

(a)

(b)

(c)

Figure 14 

Lower misalignment error estimation.

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.

Table 2 

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

Therefore, we should consider at least 𝑁=2³ different models, with different star sensor measurements covariance matrix (𝑅^star_𝑗,𝑗=1,2,…8), 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.

Table 3 

Different model for star sensor measurement noise levels.

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 𝜎_SK=1×10⁻⁸,𝜎²_ℎ𝑤=3.2×10⁻⁴,𝜎_SD=3.1×10⁻⁶,𝜎_wheel=2.1×10⁻⁶, (77) and the transition probabilities matrix 𝜋𝜋=0.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.930.010.010.010.010.010.010.010.010.93. (78) The number of the most probable model is shown in Figure 15.

Figure 15 

Number of correct model in IMM_L structure.

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 (𝑡∈[180200],[280300],[380400],… 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 (𝑡∈[100150],[200250],[300350],[400450] 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 (𝑡∈[180200],[280300],[380400] 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)

Figure 16 

Posterior probability of the EKFs structures in EIMMAE structure.

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)

Figure 17 

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

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.99900.0010.0080.990.0020.0080.0020.99. (79) 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 (𝑡∈[100150],[200250],[300350],[400450] sec) and reduces the attitude estimation errors at the end of maneuvers while the satellite is in “uniform motion” mode (𝑡∈[180200],[280300],[380400] sec) simultaneously.

(a)

(b)

(c)

Figure 18 

Estimation errors using IMM_L and EIMMAE.

(a)

(b)

(c)

Figure 19 

Estimation errors using IMM_CT and EIMMAE.

The comparison results between EIMMAE structure with IMM_L and IMM_CT for attitude estimation at the end of maneuvers (𝑡∈[180200],[280300],[380400] 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 𝑡∈[180200],[280300],[380400] sec which have been computed around each axis using (80) as follows: 1Δ𝜑=31200−180²⁰⁰_𝑡=180𝑒_𝜑(𝑡+1)300−280³⁰⁰_𝑡=280𝑒_𝜑(𝑡+1)400−380⁴⁰⁰_𝑡=380𝑒_𝜑(𝑡). (80) 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.

Table 4 

Comparing results of EIMMAE structure with IMM_L and IMM_CT methods.

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.

Table 5 

Comparison of attitude estimation methods.

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.