- Email: [email protected]
Nonlinear empirical mode predictive drift extraction on fiber optical gyroscope
Optik 122 (2011) 415–421
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Nonlinear empirical mode predictive drift extraction on fiber optical gyroscope Haiting Tian a,∗ , Ruiming Yuan a,1 , Chunxi Zhang b,2 , Jing Jin b,3 , Ningfang Song b,4 , Song Lin b,4 a b
North China Electric Power Research Institute Co., Ltd, No. 1, Dizang’an Nanxiang, Fuxingmenwai Dajie, Xicheng District, Beijing, China Institute of Opto-Electronic Technology, Beihang University, No. 37 Xueyuan Road, Haidian Block, Beijing 100191, PR China
a r t i c l e
i n f o
Article history: Received 4 August 2009 Accepted 20 February 2010
Keywords: Attitude determination Stellar sensor Gyroscope Model error filter
a b s t r a c t The integration drift of inertial unit is a fatal error on navigation and attitude determination system. Many approaches have been developed to extract and eliminate it. A novel fiber optical gyroscope drift extraction algorithm was proposed in this paper to ameliorate the performance of fiber optical gyroscope by extracting the trend component as the compensation. In this algorithm, the attitude quaternion of stellar sensor and the output of fiber optical gyroscope were imported into a nonlinear empirical mode predictive algorithm to extract the drift component of gyroscope. Combing the advantages of predictive filter and empirical mode decomposition in nonlinear signal processing, our proposal is more precise in drift extraction and data approximation. Software simulation and hardware verifications on gyroscope were launched, in which the results had proven the capability of the algorithm. © 2010 Elsevier GmbH. All rights reserved.
1. Introduction With the development of human technique in aerospace exploration, many kinds of spacecrafts have been constructed, which include satellites, space shuttles, space ships and space stations [1]. The related technology is being the most attractive realm of human science. In the voyage of a spacecraft, Gyroscope, which works as a key device for angular velocity detection, had played an important role in inertial navigation system. Because of gyroscope’s important responsibility, many modifications have been undertaken to improve its performance [2,3]. However, the drift of gyroscope, which accumulates with time, is hard to be eliminated due to the inherent flaw of inertial sensors themselves [4]. Since the gyroscope drift is a main error source in inertial system, many research works had been deployed for exploring an alternative solution of decreasing the drift of gyroscope [5,6]. As a gyroscope is a self-contained and no radiating navigation system, its error propagation cannot be separated by merely improving the precision of the sensor itself [7]. The common approach is to compensate the integration drift of gyroscope by importing the external reference information obtained from other navigation or attitude determination sensors. The Inertial Rate Sen-
∗ Corresponding author. Tel.: +86 10 8807 2458; fax: +86 10 8807 2501. E-mail addresses: [email protected] (H. Tian), [email protected] (R. Yuan), [email protected] (C. Zhang), [email protected] (J. Jin), [email protected] (N. Song), [email protected] (S. Lin). 1 Tel.: +86 10 8807 1503; fax: +86 10 8807 2501. 2 Tel.: +86 10 82316904; fax: +86 10 82316887 818. 3 Tel.: +86 10 82316887 821; fax: +86 10 82316887 818. 4 Tel.: +86 10 82316887 820; fax: +86 10 82316887 818. 0030-4026/$ – see front matter © 2010 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2010.02.024
sor (IRS)/Celestial Vector Sensor (CVS) or Global Position System (GPS)/IRS joint attitude determination algorithms have been proposed. With the advantage of no drift accumulation in CVS and GPS, the integration drift could be impaired in these joint attitude determination systems [8,9]. Amongst these joint attitude determination approaches, the attitude determination based on stellar sensor (ADSS) is enchanting as the merits of stellar sensor in the application of inter-planetary spacecraft voyage are obvious. Compared with the other CVS devices, the stellar sensor is advantageous for its high precision, solid state, high stable and wide workable range, which makes it suitable for deep space exploration [10]. In this paper, a novel fiber optical gyroscope (FOG) drift extraction algorithm is proposed. This proposal imports the stellar sensor output as the reference, and extract the drift component from FOG through a Nonlinear Empirical Mode Predictive (NEMP) process. This NEMP architecture had blended the advantage of nonlinear predictive algorithm in nonlinear data approximation and the merits of empirical mode calculation in trend component analysis. This paper has consisted of five sections. In Section 2, the related works are introduced. The NEMP algorithm on FOG drift extraction is discussed in Section 3. The verification and the simulation are listed in Section 4, and conclusion is drawn in Section 5. 2. Related works The related works are introduced in this section. In the traditional attitude determination, the Kalman filter has been proven to be extremely useful for attitude estimation using vector attitude measurements and gyro measurements. For spacecraft attitude
416
H. Tian et al. / Optik 122 (2011) 415–421
estimation, the Kalman filter is most applicable to spacecraft attitude determination with three axis gyroscopes [11]. In the application of Gyroless SAD, the analytical models of rate motion was used [12]. This approach has been successfully used in a RealTime Sequential Filter (RTSF) algorithm which propagates state estimates and error covariance using a dynamic model [12,13]. However, the RTSF is essentially a Kalman filter in which the gyro bias is modeled as a Gaussian process with known covariance. Also, fairly accurate models of angular momentum are required in order to obtain accurate estimates. Subsequently, the design process for choosing the model error covariance becomes difficult [13,14]. A new method of performing optimal state estimation in the presence of significant model error has been proposed by Mook and Junkins [15]. This method, called the Minimum Model Error (MME) estimator, is different from most filter and smoother algorithms, does not presume that the model error is represented by a Gaussian process. Instead, the model error is determined during the MME estimation process [16]. Therefore, accurate state estimates can be determined without the use of precise system representations in the assumed model. This algorithm has been successfully used to estimate the attitude of an actual spacecraft without the utilization of gyro measurements [16]. However, the MME estimator is a batch (off-line) estimator which must utilize post experiment measurements [16]. Another approach in Gyroless application SAD is the nonlinear predictive filter, which combines the good qualities of both the Kalman filter and the MME estimator [14]. The new algorithm is based on a predictive tracking scheme introduced by Prof. Lu in his publication [17]. The predictive filter algorithm developed in Prof. J. Crassidis’s paper is reformulated as a filter and estimator with a stochastic measurement process [14]. As the drift of gyroscope is a nonlinear, non-stationary random process, it is hard to be extracted by stationary separation approaches [5]. For the non-stationary signal extraction, the common used approach is regression algorithm, which builds the trend model by analyzing many sub-components’ trend. These characters make the regression algorithm complex and imprecise [18]. With the development in wavelet study, some novel gyroscope drift extraction algorithms based on wavelet analysis were proposed [19]. The wavelet transform related drift extraction algorithms need multi-size decomposition and low frequency component verification in drift extraction process, which also make it complicated in process [20]. To resolve these problems, a kind of Empirical Mode Decomposition (EMD) algorithm, which was proposed by Dr. Norden E. Huang, is adopted in gyroscope drift extraction [21,22]. Unlike other signal decomposition techniques, which map the signal space onto a space spanned by a predefined basis, the idea behind this method is to decompose a general data set into a number of basis functions termed intrinsic mode functions (IMFs), which are derived directly from the data, in a natural way [21,23].
3. Algorithm architecture The architecture of stellar sensor based FOG drift extraction system is introduced in this section, which includes an angular velocity generation block and a FOG drift extraction block. The angular velocity generation block is based on a regression process, which relies on the attitude quaternion from stellar sensor, and the FOG initial output. The generated spacecraft angular velocity was sent to the FOG drift extraction block to extract the drift trend component from FOG output. This drift extraction block is based on the NEMP algorithm to complete the calculation of error separation and drift trend extraction. The architecture of the system is shown in Fig. 1.
Fig. 1. The system architecture.
3.1. Stellar sensor The stellar sensor we used was a mature production. The field of view is 8.5◦ × 8.5◦ ; the line of sight Accuracy is 10 arcsec; and the Update rate is 10 Frames/s. The output of the stellar sensor is a quaternion sequence, which represent the attitude variation from the stellar sensor body coordinates system to the ECICS. 3.2. Reference coordinates system The RCS mentioned in the paper is J2000 Reference Coordinates, which is a sub-type of ECICS and defined as the X and Z axes point toward the mean vernal equinox and mean rotation axes of the Earth on January 1, 2000 at 12:00:00.00 Universal Time Coordinated (UTC). J2000.0 = 2000 January 1.5 = JD 2451545.0 Terrestrial Dynamical Time (TDT) [24]. 3.3. Angular velocity generation algorithm As the satellite attitude update is described by the satellite rotation angular velocity under Satellite Body Coordinates System (SBCS), which contains the components in three dimensions and should be measured by at least three gyroscopes. Based on this principle, the reference rotation angular velocity revived from stellar sensor should be mirrored on the reference frame of FOG. The stellar sensor, which exported the data under the control of synchronous trigger, outputted the satellite attitude as the quaterT nion vector q(i), q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] . Similarly, the FOG outputted the vehicle rotation angular velocity in its reference T frame, expressed as ω(i) = [ωx (i), ωy (i), ωz (i)] . Firstly, an installation error compensation for both the FOG and the stellar sensor needed to be conducted to guarantee the data obtained from these two sensors are in high precision mode. Express the FOG output data at initial time t0 as the vector ω(0) = [ωx (0), ωy (0), ωz (0)]T , where x, y, z represent the three axes of FOG output. The generated angular velocity of spacecraft is ω(i), ˆ ω(i) ˆ = T [ω ˆ x (i), ω ˆ y (i), ω ˆ z (i)] , which was calculated in this block. Before the start of the angular velocity reconstruction regression process, some initial calculations needed to be completed by using the initial measurement parameters. The initial reference angular velocity Modulus |ω(0)| was calculated as follows:
|ω(0)| =
ωx (0)2 + ωy (0)2 + ωz (0)2
(1)
The initial cosine parameter and the sine parameter were also obtained:
⎧ |ω(0)| ⎪ ⎨ c(0) = cos 2
⎪ ⎩ s(0) = sin(|ω(0)|/2) |ω(0)|
(2)
H. Tian et al. / Optik 122 (2011) 415–421
The initial exchange matrix T1 was obtained at this time:
⎡ c(0)q (1) −s(0)q (1) −s(0)q (1) −s(0)q (1) ⎤ 0 1 2 3 ⎢ c(0)q1 (1) s(0)q0 (1) −s(0)q3 (1) s(0)q2 (1) ⎥ T1 = ⎣ ⎦ c(0)q2 (1) s(0)q3 (1) s(0)q0 (1) −s(0)q1 (1) c(0)q3 (1)
−s(0)q2 (1)
s(0)q1 (1)
The function h[·] was as follows:
h(ω(i), ˆ q(i − 1)) =
s(0)q0 (1)
(4)
y
ω ˆ z (1) where the constant Nu was a redundant parameter and was abandoned in calculation. Involved the estimated angular velocity ω(1) ˆ into the regression process and the rest estimated angular velocity ω(i) ˆ with i = 1, 2, . . . , n were calculated. The iteration polynomial was shown as follows: ⎧
This drift extraction block was based on the NEMP algorithm to complete the calculation of error separation and drift trend extraction. The first launched process was the FOG error estimation, which was based on a nonlinear predictive architecture [14]. The system state function was built as follows: eω (i + 1) = f [ω(i), q(i + 1), q(i + 2), i] + I · d(i)
(6)
where eω (i) was the state estimate FOG error vector, d(i) was the model error vector, and f [·] was the quaternion to angular velocity converter algorithm, shown as: −1 (i + 1) · q(i + 2) f [ω(i), q(i + 1), q(i + 2), i] = ω(i) − Tqw
(7)
where Tqw (i) was the converter matrix, which was defined as: c(i − 1)q1 (i) c(i − 1)q2 (i) c(i − 1)q3 (i)
Tqw (i) =
−s(i − 1)q3 (i) s(i − 1)q2 (i) −s(i − 1)q1 (i) s(i − 1)q0 (i)
(9)
where v(i) was the system measurement noise, which was a GUASS white noise with zero mean and known variance R, the characteristics of v(i) was as follows:
E v(i)
=0
E v(i)vT (i)
= Rıi1i2
q(i − 1)
0 ⎢ −ωz (i)
˝(i) = ⎣
ωy (i) −ωx (i)
ωz (i)
−ωy (i)
0 −ωx (i)
ωx (i) 0
−ωy (i)
−ωz (i)
ωx (i)
⎤
ωy (i) ⎥ ⎦ ωz (i)
(12)
0
Made a Taylor series expansion on the measurement vector qˆ (i) of measurement function, the polynomial was generated as follows: qˆ (i) = q(i − 1) + Z[ω(i)] ˆ + S[ω(i)]d(i) ˆ Z[ω(i)] ˆ =
pi 1
k!
Lfk (hi ),
i = 1, 2, . . . , m
(13)
k=1
=
∂Lfk−1 (hi ) ∂ω ˆ Hi k = 0
f
k≥1
(14)
was a m × m tri-diagonal matrix, where the elements on its diagonal were ii = (1/pi !). S[ω(i)] ˆ was a m × l matrix, whose line elements were as follows: Si = Lg1 |Lfpi −1 (hi )|, . . . , Lgj |Lfpi −1 (hi )|
(15)
In the polynomial above, gj was the jth column of the model error matrix G(i), which was a unit matrix in our proposal, G(i) = I, the Lee-derivation definition of gj was: Lgj |Lfpi −1 (hi )| =
∂Lfpi −1 (hi ) ∂ω ˆ
gj ,
j = 1, 2, . . . , l
(16)
Based on the system model, a cost function could be built as follow: J=
1 1 T T [q(i) − qˆ (i)] R−1 [q(i) − qˆ (i)] + d(i) Wd(i) 2 2
(17)
where W was a positive semi-definite matrix, and the balance between two cost function relied on the selection of W. When W satisfy the covariance constraint conditions defined as follows, the balance between two cost functions was achieved. 1 T y(tk ) − h(ˆxtk , tk )y(tk ) − h(ˆxtk , tk ) R M M
(18)
k=1
(8)
The vectors c(i) and s(i) are the cosine basic function and the sine basic function, defined as above. Defined the measurement function as: qˆ (i) = h[ω(i), ˆ q(i − 1)] + v(i)
⎡
Lfk (hi )
3.4. FOG drift extraction algorithm
−s(i − 1)q2 (i) −s(i − 1)q3 (i) s(i − 1)q0 (i) s(i − 1)q1 (i)
2
(11)
(5)
An estimated spacecraft angular velocity, ω(i), ˆ could be obtained from the regression above. After then, the generated spacecraft angular velocity was sent to the FOG drift extraction block to continue the drift trend component extraction process.
−s(i − 1)q1 (i) s(i − 1)q0 (i) s(i − 1)q3 (i) −s(i − 1)q2 (i)
2
|ω(i)|
where pi was the lowest step of d(i) that appeared in ω(i) ˆ when the vector qˆ (i) was made a derivative. pi was also called as related step. Lfk (hi ) was the k step Lee-derivative of hi on f (·), which was defined as follows:
ω ˆ z (i)
c(i − 1)q0 (i)
cos
˝(i) ·I+ · sin m
where ˝ was an angular velocity matrix, represented as:
⎡ Nu ⎤ Nu ⎢ ωˆ x (1) ⎥ =⎣ ⎦ = T1−1 · q(2) ω(c) ˆ ω ˆ (1)
−s(i − 1)q3 (i) s(i − 1)q2 (i) −s(i − 1)q1 (i) s(i − 1)q0 (i)
|ω(i)|
(3)
The estimated angular velocity ω(1) ˆ could be generated through the attitude quaternion vector q(2) and the T1−1 , which was the inversion of the exchange matrix T1 .
2 2 2 |ω(i ˆ − 1)| = ω ˆ x (i − 1) + ω ˆ y (i − 1) + ω ˆ z (i − 1) ⎪ ⎪ ⎪ ⎪ |ω(i ˆ − 1)| ⎪ ⎪ c(i − 1) = cos ⎪ 2 ⎪ ⎪ ⎪ ⎪ sin(|ω(i ˆ − 1)|/2) ⎪ ⎪ s(i − 1) = ⎪ |ω(i ˆ − 1)| ⎨ c(i − 1)q0 (i) −s(i − 1)q1 (i) −s(i − 1)q2 (i) ⎪ Ti = c(i − 1)q1 (i) s(i − 1)q0 (i) −s(i − 1)q3 (i) ⎪ ⎪ c(i − 1)q2 (i) s(i − 1)q3 (i) s(i − 1)q0 (i) ⎪ ⎪ ⎪ c(i − 1)q3 (i) −s(i − 1)q2 (i) s(i − 1)q1 (i) ⎪ x ⎪ ⎪ ⎪ ⎪ ω ˆ x (i) ⎪ ⎪ = Ti−1 · q(i + 1) + d(i) ⎪ ⎩ ωˆ y (i)
417
(10)
where M was the total measurement point number. Based on the model error minimization criterion, minimized the model error cost function J, we could get the representation of model error as follow: d(i) =
T
[S[ω(i)]] ˆ R−1 [S[ω(i)]] ˆ +W T
−1
·[S[ω(i)]] ˆ R−1 [Z[ω(i)] ˆ − e(i) − b(i)]
(19)
Imported the model error d(i) back in to the system state function to modify the system error state vectoreω (i + 1) and continue the regression process. The system error state vectoreω (i + 1) represented the error component of FOG.
418
H. Tian et al. / Optik 122 (2011) 415–421
At second step, a drift trend extraction process was deployed to filter out the noise and short time interference from the separated error component of FOG. This extraction calculation was based on the empirical model decomposition (EMD) [21]. The target of drift trend extraction was to decompose an original gyroscope output signal x(i) as a intrinsic mode component h(i) and a final residual r(i). The FOG error components with the different time step-size were represented as the different Intrinsic Mode Functions (IMFs), K h(i) = h (i), where the hk (i) was a IMFs sequence. As the EMD k=1 k algorithm was based on the time local trait, it could be used in nonlinear and non-stationary process as a sifting function [22]. For a real-valued signal, the EMD performed the mapping in equation below. x(i) = h(i) + r(i) =
K
hk (i) + r(i),
(20)
k=1
where dk (i) was a set of IMFs, and K was the number of IMF. The r(i) was the residual. The first IMF was calculated in envelope calculation block and IMFs separation block. The drift trend extraction algorithm contained three blocks, which were envelope calculation block, IMFs separation block and condition judgment block. The ANFIS interpolation unit was embedded in the envelope calculation block to approximate the maximum and minimum envelope by the limited sample data. The condition judgment block was to judge whether the output data from the IMFs separation block contained another IMF or the final residual data. The error component of FOG, which was separated before, was sent to the envelope calculation block as the input signal. The envelope calculation block responded to identify all extrema of x(i) and then Interpolate between maxima and minima of x(i). This block ended up with the maximum envelope emax (i) and minimum envelope emin (i). The envelope interpolation method in our paper was developed from ANFIS algorithm in the ANFIS interpolation unit. With the inference capability of fuzzy logic and the learning capability of neural networks, the ANFIS, which is based on the utilization of first order SUGENO model, is more flexible and has a preferable capability in approximating a linear or nonlinear function [25,26]. These advantages make it very suitable in this nonlinear envelope interpolation. Input signals of the ANFIS model were time sequence value t(i). The output signals were the ANFIS approximations of maximum envelope emax (i) and minimum envelope emin (i), defined as eˆmax (i) and eˆmin (i). Considered the approximation, the time sequence value t(i) could be divided into 4 parts on average and a bell shaped function was used as the membership function of a network fuzzy set, which was defined as follows: (x) =
1
1+|
x−ci 2b , ai |
(21)
where x was the input variable, and a, b and c were premise parameters. There were 5 layer nodes and 16 fuzzy if-then rules in the structure of ANFIS, in which 16 membership functions and 24 premise parameters were for layer 1, and 48 consequent parameters were for layer 4. There were totally 864 fitting parameters needed in the conventional ANFIS structure. Compared with the conventional cubic spline interpolation, the ANFIS approach with nonlinear multi layers structure was better in nonlinear data approximation. The ANFIS structure needed to be trained before approximation [27]. The training data was the sampled maxima data xmax (i) and minima data xmin (i) and the responded time point t(i), which could be expressed as a vector array, [xmax (i), xmin (i), t(i)]. In ANFIS training process, the time point t(i) was set as the input and the data xmax (i) and xmin (i) were treated as the output.
A mixture algorithm has been used to train parameters of ANFIS for quick parameters identification [27,28]. The training data was brought to the inputs, the premise parameters were assumed to be fixed and the optimal consequent parameters were estimated by an iterative least mean square procedure in the forward pass. In the reverse pass, the patterns were propagated again, but this time the consequent parameters were assumed to be fixed, and the back-propagation was used to modify the premise parameters [26]. The iteration training number of the ANFIS model was 1000; the minimum square estimation between the practical system output and the model output was used as a network training index. After the training, the ANFIS approximation network would be used to approximate the maximum envelope emax (i) and the minimum envelope emin (i). We treated a serial time sequence, t(i), where i was a continuous time value, as the approximation network input. The output of the ANFIS network would be the approximation of the maximum envelope, which was defined as eˆmax (i), and the minimum envelope, which was defined as eˆmin (i). The IMFs separation block, which computed the mean m(i) = (emax (i) + emin (i))/2 and Subtracted from the signal to obtain h(i) = x(i) − m(i). There were two conditions needed to be satisfied if h(t) was IMF, which was introduced in paper [29]. Ideally, h(t) should be an IMF since the construction of h(t) described above has made it satisfy all the requirements of IMF. If h(t) did not meet the definition of a true IMF, the sifting process had to be repeated many times to correct this until the first IMF imf1 was obtained. After the first IMF was separated, it was need to judge whether the output data from the IMFs separation block contains another IMF [22]. This process was completed in the condition judgment block. The first IMF, imf1 , should contain the finest scale or the shortest period component of the signal. We could separate imf1 from the rest of the data by: h(i) − imf1 (i) = r1 (i),
(22)
Since the residue, r1 (i), still contained long period components, it should be treated as the new data and subject to the same sifting process as described above. The sifting process stopped by any of the following predetermined criteria: either when the component, imfn , or the residue, rn (i), becomed so small and less than the pre-determined value of substantial consequence, or when the residue, rn (i), becomed a monotonic function and no another IMF can be extracted from it. After the decomposition, the gyroscope signal had been separated in to a monotonic residual and the IMFs components [21]. The final residual rN (i) was a zero mean noise process, which did not contain drift trend information. However, the IMFs components separated in this section contained the trend characters that had disguised in FOG error. All of these trend characters could be utilized by navigation system to calibrate or compensate the performance of FOG. As these trend characters in IMFS had been decomposed into different time step-size domains, it had become easy to select the notable component as the main drift trend of FOG.
4. Verification and simulation To verify the performance of our FOG drift extraction algorithm, software simulation and FOG testing experiments were launched. The related simulation and experiments were grouped into three types, which were angular velocity generation algorithm simulation, FOG trend extraction verification and system experiment on FOG.
H. Tian et al. / Optik 122 (2011) 415–421
419
Table 1 The deviations of the recovered angular velocities of circular orbit. Axis
Mean deviation (rad/s)
X axis Y axis Z axis
−7
2.282 × 10 2.261 × 10−7 2.248 × 10−7
Standard deviation (rad/s) 2.846 × 10−7 2.820 × 10−7 2.841 × 10−7
Table 2 The deviations of the recovered angular velocities of circular orbit using conventional algorithm.
Fig. 2. The spacecraft attitude quaternion of circular orbit.
4.1. Angular velocity generation algorithm simulation The angular velocity generation algorithm was simulated on ground. The vehicle attitude information was generated by the spacecraft orbit simulation software Satellite Tools Kit (STK), which was a powerful satellite simulation soft developed by Analytical Graphics, Inc. (AGI). The reference time format mentioned in the simulation was J2000 Reference Coordinates, which was a subtype of ECICS and defined as the X and Z axes point toward the mean vernal equinox and mean rotation axes of the Earth on January 1, 2000 at 12:00:00.00 Universal Time Coordinated (UTC). J2000.0 = 2000 January 1.5 = JD 2451545.0 Terrestrial Dynamical Time (TDT) Two individual satellite orbits were built to verify the performance of our algorithm in different orbit and attitude characteristics. 1. Circular orbit verification The first algorithm simulation was in a circular orbit condition, which was generated from the High-Precision Orbit Propagator of the STK simulation platform. The parameters of this orbit were as follows: • • • • • • •
Semi major axis: 6878.137 km Eccentricity: 1.85671◦ × 10−16 Inclination: 45◦ Period: 5676.978029 s Argument of Perigee: 0◦ Mean anomaly: 30◦ RAAN: 0◦
Axis
Mean deviation (rad/s)
Standard deviation (rad/s)
X axis Y axis Z axis
6.642 × 10−7 7.302 × 10−7 6.665 × 10−7
8.873 × 10−7 9.691 × 10−7 8.635 × 10−7
software platform, are shown in Fig. 2. The estimated angular velocities generated by our angular velocity generation algorithm are shown in Fig. 3. The deviations of the recovered angular velocities, which included the mean error and the standard deviation results, are listed in Table 1. Shown from the verification results, the attitude recovering deviations were about 1 percent of the input angular velocity range. For these deviations, the data acquisition error of STK software, which had a round-off error in small data calculation, was the most contribution. As this phenomenon might also appears in the data acquisition process of stellar sensor, an external filter should be utilized to smooth the data. As a comparison, another group of simulation was also launched using the algorithm that introduced in article [30]. The deviations of the recovered angular velocities in circular orbit using conventional algorithm are listed in Table 2. 2. Critically inclined orbit verification The second algorithm simulation was in a critically inclined orbit condition, and the parameters of this orbit are as follows: • • • • • • •
Semi major axis: 12,578.137 km Eccentricity: 0.461118◦ Inclination: 63.434949◦ Period: 14038.976969 s Argument of Perigee: 0◦ RAAN: 0◦ Mean anomaly: 329.517505◦
The spacecraft attitude vector, expressed as the quaternion T q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] , which were generated from STK
The spacecraft attitude vector, expressed as the quaternion T q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] , which were generated from STK software platform, are shown in Fig. 4.
Fig. 3. The estimated angular velocities of circular orbit.
Fig. 4. The spacecraft attitude quaternion of critically inclined orbit.
420
H. Tian et al. / Optik 122 (2011) 415–421 Table 5 FOG drift extraction results by the trend extraction algorithm. The final residuals after the NFL-EMD Angular velocity
Mean
Standard deviation
Variance
0◦ /s 1◦ /s 10/s 100◦ /s
0.3875 −3.3 −1.3373 0.4303
5.5843 0.0369 2.7572 2.2696
31.1846 1.3625 75.0559 51.5117
Table 6 FOG drift extraction results by the conventional EMD algorithm. The final residuals after the conventional EMD Fig. 5. The estimated angular velocity of critically inclined orbit.
Angular velocity ◦
The estimated angular velocities of critically inclined orbit generated by our angular velocity generation algorithm are shown in Fig. 5. The deviations of the recovered angular velocities, which included the mean error and the standard deviation results, are listed in Table 3. The deviations of the recovered angular velocities in critically inclined orbit using conventional algorithm are listed in Table 4.
0 /s 1◦ /s 10◦ /s 100◦ /s
Mean
Standard deviation
Variance
0.705 −12.1 −2.81 0.6825
14.9735 0.1375 13.3541 9.381
64.7358 12.276 35.4732 46.6724
4.2. FOG trend extraction verification To verify the performance of our FOG trend extraction algorithm, we had launched a specified verification on FOG. The FOG we used in experiment was in middle precision (with the original bias stability of 1◦ /s). After being used and frayed for a long time, the output signal of this sensor was blended with some drift components and the precision has decreased. The bias stability has become as 3◦ /s. Four test angular velocities were set as the input parameter, they were 0, 1, 10 and 100◦ /s. The output data of the FOG were sent to the trend extraction program. The ANFIS interpolation approach used in the trend extraction algorithm was the same as the architecture mentioned above, in which the membership functions number for each input was 10 and the membership function type was generalized bell membership function. The intrinsic mode data, which was extracted from the FOG error data, could represent the drift components in different time step-size. An analysis had been made to examine the statistic characteristics (Mean, Standard deviation and Variance) of the final residuals in four different input angular velocities, and the analysis results are listed in Table 5. As the comparison, we used the conventional EMD algorithm to complete another drift extraction experiment. The architecture of
Fig. 6. The import attitude quaternion of FOG.
Table 3 The deviations of the recovered rotation velocities of critically inclined orbit. Axis X axis Y axis Z axis
Mean deviation (rad/s) −7
2.337 × 10 3.632 × 10−7 2.230 × 10−7
Standard deviation (rad/s)
Fig. 7. The FOG output data.
−7
2.943 × 10 4.881 × 10−7 2.806 × 10−7
Table 4 The deviations of the recovered rotation velocities of critically inclined orbit by conventional algorithm. Axis
Mean deviation (rad/s)
Standard deviation (rad/s)
X axis Y axis Z axis
6.621 × 10−7 7.387 × 10−7 6.720 × 10−7
9.071 × 10−7 8.672 × 10−7 9.607 × 10−7
this EMD algorithm was introduced in paper [21] and the interpolation function was the cubic spline function. The analysis results of this verification data are listed in Table 6. Compared with the conventional EMD algorithm, our trend extraction algorithm based on nonlinear empirical model predictive architecture could achieve higher performance and the results had shown that the bias stability of the final residual in our drift extraction was 0.62 times as the residual in conventional EMD extraction.
H. Tian et al. / Optik 122 (2011) 415–421
Fig. 8. The trend components in three axes.
4.3. System performance testing To test the performance of our proposal in real FOG operation, FOG drift extraction verification was launched and introduced in this section. It was known that the angular velocity of orbit spacecrafts are about −1◦ /s → +1◦ /s, which is a small velocity importation. We set the geostationary orbit satellite angular velocity as the import angular velocity of the FOG. The FOG was fixed on a revolving table for simulating the rotation of the satellite. There were three axes contained in this FOG, the X axis of it pointed to the sky, the Y axis pointed to the north direction, the Z axis pointed to the east direction. The import attitude quaternion of FOG are shown in Fig. 6. The FOG output data are shown in Fig. 7. Extracting the drift trend components by our algorithm, the results of the trend components in three axes are shown in Fig. 8. 5. Conclusion An FOG drift extraction algorithm was exhibited in the paper and the related simulations were discussed. The algorithm, which could extract the drift trend from the output data of FOG, made it possible to eliminate the integration drift of gyroscope. The utilization of this approach could be on spacecraft, satellite and space-shuttle. The algorithm verification was held on ground, which had verified the performance of the algorithm and the results were described. Some modifications have been deployed to gain a more approximation of the actual angular velocity. Besides, the ground verification of algorithm was in good condition, which had not considered the external noise and potential interference so much. The experiment in space orbit environment will be hold in next step to verify the algorithm more deeply. The observation of large points, which appear at the simulation result sporadically, is a potential threat to the application and need to be eliminated. The related modification is under development currently. Acknowledgment The author gratefully acknowledges Mrs. Li Min, who has made remarkable contribution for this paper. References [1] M. Kuritsky, M. Goldstein, I. Greenwood, H. Lerman, J. McCarthy, T. Shanahan, M. Silver, J. Simpson, Inertial navigation, Proc. IEEE 71 (October (10)) (1983) 1156–1176.
421
[2] R. Zhu, Y. Zhang, Q. Bao, A novel intelligent strategy for improving measurement precision of fog, IEEE Trans. Instrum. Meas. 49 (December (6)) (2000) 1183–1188. [3] L.G. Richard, Initial navigation technology from 1970–1995, J. Inst. Navigat. 42 (January (1)) (1995) 165–185. [4] Y. Akcayir, Y. Ozkazanc, Gyroscope drift estimation analysis in land navigation systems Control Applications, in: CCA 2003. Proceedings of 2003 IEEE Conference, vol. 2, June 2003, pp. 1488–1491. [5] L. Miao, J.F. Zhang, S.J. Shen, Data analysis and modeling of fiber optic gyroscope drift, J. Beij. Inst. Technol. 11 (January (1)) (2002) 50–55. [6] J. Li, H. Wang, A. Liu, Trend extraction of gyro’s drift based on modified empirical mode decomposition, Syst. Eng. Electron. 27 (June (6)) (2005) 1080–1082. [7] R. Zhu, Y. Zhang, Q. Bao, A novel intelligent strategy for improving measurement precision of fog, IEEE Trans. Instrum. Meas. 49 (December (6)) (2000) 1183–1188. [8] H.S. Ahn, S.H. Lee, Gyroless attitude estimation of sun-pointing satellites using magnetometers, Geosci. Rem. Sens. Lett. IEEE 2 (January (1)) (2005) 8–12. [9] M. Grewal, M. Shiva, Application of kalman filtering to gyroless attitude determination and control system for environmental satellites, in: Decision and Control, 1995, Proceedings of the 34th IEEE Conference, vol. 2, December 1995. [10] J.C. Yoon, B.S. Lee, K.H. Choi, Spacecraft orbit determination using gps navigation solutions, Aerospace Sci. Technol. 4 (April (3)) (2000). [11] G.-S. Huang, J.-C. Juang, Application of nonlinear kalman filter approach in dynamic gps-based attitude determination, in: Circuits and Systems, 1997. Proceedings of the 40th Midwest Symposium, vol. 2, August 1997, pp. 1440–1444. [12] M. Grewal, M. Shiva, Application of kalman filtering to gyroless attitude determination and control system for environmental satellites, in: Decision and Control, 1995. Proceedings of the 34th IEEE Conference, vol. 2, December 1995, pp. 1544–1552. [13] M. Challa Solar, Anomalous, and Magnetospheric Particle Explorer (SAMPEX) Real Time Sequential Filter (RTSF), Evaluation Report, NASA-Goddard Space Flight Center, Greenbelt, MD, USA, April 1993. [14] J. Crassidis, F. Markley, Predictive filtering for attitude estimation without rate sensors, Journal of Guidance, Contr. Dynam. 20 (May–June (3)) (1997) 522–527. [15] D. Mook, J. Junkins, Minimum model error estimation for poorly modeled dynamic systems, J. Guidance Contr. Dynam. 3 (January–February (4)) (1988) 367–375. [16] V. Parameswaran, J. Raol, Estimation of model error for nonlinear system identification, in: Control Theory and Applications, IEE Proceedings, vol. 141, no. 6, November 1994, pp. 403–408. [17] P. Lu, Nonlinear predictive controllers for continuous systems, J. Guidance Contr. Dynam. 17 (May–June (3)) (1994) 553–560. [18] H. Wang, J. Chen, W. Huang, Y. Zhang, J. Tang, Detection and analysis of the fiber optic gyroscope’s random drift, Opt. Tech. 30 (September(5)) (2004) 623–627. [19] W.W. Dai, Z. Yang, S.L. Lim, O. Mikhailova, J. Chee, Processing and analysis of ecg signal using nonorthogonal wavelet transform, in: Engineering in Medicine and Biology Society, 1998. Proceedings of the 20th Annual International Conference of the IEEE, vol. 1, October–November 1998, pp. 139–142. [20] Y. Zheng, D.B.H. Tay, L. Li, Signal extraction and power spectrum estimation using wavelet transform scale space filtering and bayes shrinkage, Signal Process. 80 (August (8)) (2000) 1535–1549. [21] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. A: Math. Phys. Eng. Sci. 454 (March (1971)) (1998) 903–995. [22] H. Zhang, Q. Gai, Research on properties of empirical mode decomposition method, Intelligent Control and Automation, 2006. WCICA 2006. The Sixth World Congress, vol. 2, June 2006, pp. 10001–10004. [23] T. Tanaka, D.P. Mandic, Complex empirical mode decomposition, Signal Process. Lett. IEEE 14 (February (2)) (2007) 101–104. [24] A. Brown, G. Moy, Long duration strapdown stellar-inertial navigation using satellite tracking, in: Position Location and Navigation Symposium, 1992. Record. 500 Years After Columbus—Navigation Challenges of Tomorrow, IEEE PLANS’92,IEEE, March 1992, pp. 194–201. [25] C. Lee, J.-S. Wang, Self-adaptive neuro-fuzzy systems: structure and learning Intelligent Robots and Systems, (IROS 2000), in: Proceedings. 2000 IEEE/RSJ International Conference, vol. 1, 2000, pp. 52–57. [26] M. Sugeno, G.T. Kang, Structure identification of fuzzy model, Fuzzy Sets Syst. 28 (October (1)) (1988) 15–33. [27] J.-S. Jang, Anfis: adaptive-network-based fuzzy inference system, IEEE Trans. Syst. Man Cybernet. 23 (May/June (3)) (1993) 665–685. [28] T. Amaral, M. Crisostomo, V. Pires, Adaptive neuro-fuzzy inference system for modelling and control Intelligent Systems, in: Proceedings. 2002 First International IEEE Symposium, vol. 1, 2002, pp. 67–72. [29] P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, Signal Process. Lett. IEEE 11 (February (2)) (2004) 112–114. [30] X.Q. Chen, Y.H. Geng, On orbit calibration algorithm with star sensors of gyros, Syst. Eng. Electron. 27 (December (12)) (2005) 2112–2116.
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.de/ijleo
Nonlinear empirical mode predictive drift extraction on fiber optical gyroscope Haiting Tian a,∗ , Ruiming Yuan a,1 , Chunxi Zhang b,2 , Jing Jin b,3 , Ningfang Song b,4 , Song Lin b,4 a b
North China Electric Power Research Institute Co., Ltd, No. 1, Dizang’an Nanxiang, Fuxingmenwai Dajie, Xicheng District, Beijing, China Institute of Opto-Electronic Technology, Beihang University, No. 37 Xueyuan Road, Haidian Block, Beijing 100191, PR China
a r t i c l e
i n f o
Article history: Received 4 August 2009 Accepted 20 February 2010
Keywords: Attitude determination Stellar sensor Gyroscope Model error filter
a b s t r a c t The integration drift of inertial unit is a fatal error on navigation and attitude determination system. Many approaches have been developed to extract and eliminate it. A novel fiber optical gyroscope drift extraction algorithm was proposed in this paper to ameliorate the performance of fiber optical gyroscope by extracting the trend component as the compensation. In this algorithm, the attitude quaternion of stellar sensor and the output of fiber optical gyroscope were imported into a nonlinear empirical mode predictive algorithm to extract the drift component of gyroscope. Combing the advantages of predictive filter and empirical mode decomposition in nonlinear signal processing, our proposal is more precise in drift extraction and data approximation. Software simulation and hardware verifications on gyroscope were launched, in which the results had proven the capability of the algorithm. © 2010 Elsevier GmbH. All rights reserved.
1. Introduction With the development of human technique in aerospace exploration, many kinds of spacecrafts have been constructed, which include satellites, space shuttles, space ships and space stations [1]. The related technology is being the most attractive realm of human science. In the voyage of a spacecraft, Gyroscope, which works as a key device for angular velocity detection, had played an important role in inertial navigation system. Because of gyroscope’s important responsibility, many modifications have been undertaken to improve its performance [2,3]. However, the drift of gyroscope, which accumulates with time, is hard to be eliminated due to the inherent flaw of inertial sensors themselves [4]. Since the gyroscope drift is a main error source in inertial system, many research works had been deployed for exploring an alternative solution of decreasing the drift of gyroscope [5,6]. As a gyroscope is a self-contained and no radiating navigation system, its error propagation cannot be separated by merely improving the precision of the sensor itself [7]. The common approach is to compensate the integration drift of gyroscope by importing the external reference information obtained from other navigation or attitude determination sensors. The Inertial Rate Sen-
∗ Corresponding author. Tel.: +86 10 8807 2458; fax: +86 10 8807 2501. E-mail addresses: [email protected] (H. Tian), [email protected] (R. Yuan), [email protected] (C. Zhang), [email protected] (J. Jin), [email protected] (N. Song), [email protected] (S. Lin). 1 Tel.: +86 10 8807 1503; fax: +86 10 8807 2501. 2 Tel.: +86 10 82316904; fax: +86 10 82316887 818. 3 Tel.: +86 10 82316887 821; fax: +86 10 82316887 818. 4 Tel.: +86 10 82316887 820; fax: +86 10 82316887 818. 0030-4026/$ – see front matter © 2010 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2010.02.024
sor (IRS)/Celestial Vector Sensor (CVS) or Global Position System (GPS)/IRS joint attitude determination algorithms have been proposed. With the advantage of no drift accumulation in CVS and GPS, the integration drift could be impaired in these joint attitude determination systems [8,9]. Amongst these joint attitude determination approaches, the attitude determination based on stellar sensor (ADSS) is enchanting as the merits of stellar sensor in the application of inter-planetary spacecraft voyage are obvious. Compared with the other CVS devices, the stellar sensor is advantageous for its high precision, solid state, high stable and wide workable range, which makes it suitable for deep space exploration [10]. In this paper, a novel fiber optical gyroscope (FOG) drift extraction algorithm is proposed. This proposal imports the stellar sensor output as the reference, and extract the drift component from FOG through a Nonlinear Empirical Mode Predictive (NEMP) process. This NEMP architecture had blended the advantage of nonlinear predictive algorithm in nonlinear data approximation and the merits of empirical mode calculation in trend component analysis. This paper has consisted of five sections. In Section 2, the related works are introduced. The NEMP algorithm on FOG drift extraction is discussed in Section 3. The verification and the simulation are listed in Section 4, and conclusion is drawn in Section 5. 2. Related works The related works are introduced in this section. In the traditional attitude determination, the Kalman filter has been proven to be extremely useful for attitude estimation using vector attitude measurements and gyro measurements. For spacecraft attitude
416
H. Tian et al. / Optik 122 (2011) 415–421
estimation, the Kalman filter is most applicable to spacecraft attitude determination with three axis gyroscopes [11]. In the application of Gyroless SAD, the analytical models of rate motion was used [12]. This approach has been successfully used in a RealTime Sequential Filter (RTSF) algorithm which propagates state estimates and error covariance using a dynamic model [12,13]. However, the RTSF is essentially a Kalman filter in which the gyro bias is modeled as a Gaussian process with known covariance. Also, fairly accurate models of angular momentum are required in order to obtain accurate estimates. Subsequently, the design process for choosing the model error covariance becomes difficult [13,14]. A new method of performing optimal state estimation in the presence of significant model error has been proposed by Mook and Junkins [15]. This method, called the Minimum Model Error (MME) estimator, is different from most filter and smoother algorithms, does not presume that the model error is represented by a Gaussian process. Instead, the model error is determined during the MME estimation process [16]. Therefore, accurate state estimates can be determined without the use of precise system representations in the assumed model. This algorithm has been successfully used to estimate the attitude of an actual spacecraft without the utilization of gyro measurements [16]. However, the MME estimator is a batch (off-line) estimator which must utilize post experiment measurements [16]. Another approach in Gyroless application SAD is the nonlinear predictive filter, which combines the good qualities of both the Kalman filter and the MME estimator [14]. The new algorithm is based on a predictive tracking scheme introduced by Prof. Lu in his publication [17]. The predictive filter algorithm developed in Prof. J. Crassidis’s paper is reformulated as a filter and estimator with a stochastic measurement process [14]. As the drift of gyroscope is a nonlinear, non-stationary random process, it is hard to be extracted by stationary separation approaches [5]. For the non-stationary signal extraction, the common used approach is regression algorithm, which builds the trend model by analyzing many sub-components’ trend. These characters make the regression algorithm complex and imprecise [18]. With the development in wavelet study, some novel gyroscope drift extraction algorithms based on wavelet analysis were proposed [19]. The wavelet transform related drift extraction algorithms need multi-size decomposition and low frequency component verification in drift extraction process, which also make it complicated in process [20]. To resolve these problems, a kind of Empirical Mode Decomposition (EMD) algorithm, which was proposed by Dr. Norden E. Huang, is adopted in gyroscope drift extraction [21,22]. Unlike other signal decomposition techniques, which map the signal space onto a space spanned by a predefined basis, the idea behind this method is to decompose a general data set into a number of basis functions termed intrinsic mode functions (IMFs), which are derived directly from the data, in a natural way [21,23].
3. Algorithm architecture The architecture of stellar sensor based FOG drift extraction system is introduced in this section, which includes an angular velocity generation block and a FOG drift extraction block. The angular velocity generation block is based on a regression process, which relies on the attitude quaternion from stellar sensor, and the FOG initial output. The generated spacecraft angular velocity was sent to the FOG drift extraction block to extract the drift trend component from FOG output. This drift extraction block is based on the NEMP algorithm to complete the calculation of error separation and drift trend extraction. The architecture of the system is shown in Fig. 1.
Fig. 1. The system architecture.
3.1. Stellar sensor The stellar sensor we used was a mature production. The field of view is 8.5◦ × 8.5◦ ; the line of sight Accuracy is 10 arcsec; and the Update rate is 10 Frames/s. The output of the stellar sensor is a quaternion sequence, which represent the attitude variation from the stellar sensor body coordinates system to the ECICS. 3.2. Reference coordinates system The RCS mentioned in the paper is J2000 Reference Coordinates, which is a sub-type of ECICS and defined as the X and Z axes point toward the mean vernal equinox and mean rotation axes of the Earth on January 1, 2000 at 12:00:00.00 Universal Time Coordinated (UTC). J2000.0 = 2000 January 1.5 = JD 2451545.0 Terrestrial Dynamical Time (TDT) [24]. 3.3. Angular velocity generation algorithm As the satellite attitude update is described by the satellite rotation angular velocity under Satellite Body Coordinates System (SBCS), which contains the components in three dimensions and should be measured by at least three gyroscopes. Based on this principle, the reference rotation angular velocity revived from stellar sensor should be mirrored on the reference frame of FOG. The stellar sensor, which exported the data under the control of synchronous trigger, outputted the satellite attitude as the quaterT nion vector q(i), q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] . Similarly, the FOG outputted the vehicle rotation angular velocity in its reference T frame, expressed as ω(i) = [ωx (i), ωy (i), ωz (i)] . Firstly, an installation error compensation for both the FOG and the stellar sensor needed to be conducted to guarantee the data obtained from these two sensors are in high precision mode. Express the FOG output data at initial time t0 as the vector ω(0) = [ωx (0), ωy (0), ωz (0)]T , where x, y, z represent the three axes of FOG output. The generated angular velocity of spacecraft is ω(i), ˆ ω(i) ˆ = T [ω ˆ x (i), ω ˆ y (i), ω ˆ z (i)] , which was calculated in this block. Before the start of the angular velocity reconstruction regression process, some initial calculations needed to be completed by using the initial measurement parameters. The initial reference angular velocity Modulus |ω(0)| was calculated as follows:
|ω(0)| =
ωx (0)2 + ωy (0)2 + ωz (0)2
(1)
The initial cosine parameter and the sine parameter were also obtained:
⎧ |ω(0)| ⎪ ⎨ c(0) = cos 2
⎪ ⎩ s(0) = sin(|ω(0)|/2) |ω(0)|
(2)
H. Tian et al. / Optik 122 (2011) 415–421
The initial exchange matrix T1 was obtained at this time:
⎡ c(0)q (1) −s(0)q (1) −s(0)q (1) −s(0)q (1) ⎤ 0 1 2 3 ⎢ c(0)q1 (1) s(0)q0 (1) −s(0)q3 (1) s(0)q2 (1) ⎥ T1 = ⎣ ⎦ c(0)q2 (1) s(0)q3 (1) s(0)q0 (1) −s(0)q1 (1) c(0)q3 (1)
−s(0)q2 (1)
s(0)q1 (1)
The function h[·] was as follows:
h(ω(i), ˆ q(i − 1)) =
s(0)q0 (1)
(4)
y
ω ˆ z (1) where the constant Nu was a redundant parameter and was abandoned in calculation. Involved the estimated angular velocity ω(1) ˆ into the regression process and the rest estimated angular velocity ω(i) ˆ with i = 1, 2, . . . , n were calculated. The iteration polynomial was shown as follows: ⎧
This drift extraction block was based on the NEMP algorithm to complete the calculation of error separation and drift trend extraction. The first launched process was the FOG error estimation, which was based on a nonlinear predictive architecture [14]. The system state function was built as follows: eω (i + 1) = f [ω(i), q(i + 1), q(i + 2), i] + I · d(i)
(6)
where eω (i) was the state estimate FOG error vector, d(i) was the model error vector, and f [·] was the quaternion to angular velocity converter algorithm, shown as: −1 (i + 1) · q(i + 2) f [ω(i), q(i + 1), q(i + 2), i] = ω(i) − Tqw
(7)
where Tqw (i) was the converter matrix, which was defined as: c(i − 1)q1 (i) c(i − 1)q2 (i) c(i − 1)q3 (i)
Tqw (i) =
−s(i − 1)q3 (i) s(i − 1)q2 (i) −s(i − 1)q1 (i) s(i − 1)q0 (i)
(9)
where v(i) was the system measurement noise, which was a GUASS white noise with zero mean and known variance R, the characteristics of v(i) was as follows:
E v(i)
=0
E v(i)vT (i)
= Rıi1i2
q(i − 1)
0 ⎢ −ωz (i)
˝(i) = ⎣
ωy (i) −ωx (i)
ωz (i)
−ωy (i)
0 −ωx (i)
ωx (i) 0
−ωy (i)
−ωz (i)
ωx (i)
⎤
ωy (i) ⎥ ⎦ ωz (i)
(12)
0
Made a Taylor series expansion on the measurement vector qˆ (i) of measurement function, the polynomial was generated as follows: qˆ (i) = q(i − 1) + Z[ω(i)] ˆ + S[ω(i)]d(i) ˆ Z[ω(i)] ˆ =
pi 1
k!
Lfk (hi ),
i = 1, 2, . . . , m
(13)
k=1
=
∂Lfk−1 (hi ) ∂ω ˆ Hi k = 0
f
k≥1
(14)
was a m × m tri-diagonal matrix, where the elements on its diagonal were ii = (1/pi !). S[ω(i)] ˆ was a m × l matrix, whose line elements were as follows: Si = Lg1 |Lfpi −1 (hi )|, . . . , Lgj |Lfpi −1 (hi )|
(15)
In the polynomial above, gj was the jth column of the model error matrix G(i), which was a unit matrix in our proposal, G(i) = I, the Lee-derivation definition of gj was: Lgj |Lfpi −1 (hi )| =
∂Lfpi −1 (hi ) ∂ω ˆ
gj ,
j = 1, 2, . . . , l
(16)
Based on the system model, a cost function could be built as follow: J=
1 1 T T [q(i) − qˆ (i)] R−1 [q(i) − qˆ (i)] + d(i) Wd(i) 2 2
(17)
where W was a positive semi-definite matrix, and the balance between two cost function relied on the selection of W. When W satisfy the covariance constraint conditions defined as follows, the balance between two cost functions was achieved. 1 T y(tk ) − h(ˆxtk , tk )y(tk ) − h(ˆxtk , tk ) R M M
(18)
k=1
(8)
The vectors c(i) and s(i) are the cosine basic function and the sine basic function, defined as above. Defined the measurement function as: qˆ (i) = h[ω(i), ˆ q(i − 1)] + v(i)
⎡
Lfk (hi )
3.4. FOG drift extraction algorithm
−s(i − 1)q2 (i) −s(i − 1)q3 (i) s(i − 1)q0 (i) s(i − 1)q1 (i)
2
(11)
(5)
An estimated spacecraft angular velocity, ω(i), ˆ could be obtained from the regression above. After then, the generated spacecraft angular velocity was sent to the FOG drift extraction block to continue the drift trend component extraction process.
−s(i − 1)q1 (i) s(i − 1)q0 (i) s(i − 1)q3 (i) −s(i − 1)q2 (i)
2
|ω(i)|
where pi was the lowest step of d(i) that appeared in ω(i) ˆ when the vector qˆ (i) was made a derivative. pi was also called as related step. Lfk (hi ) was the k step Lee-derivative of hi on f (·), which was defined as follows:
ω ˆ z (i)
c(i − 1)q0 (i)
cos
˝(i) ·I+ · sin m
where ˝ was an angular velocity matrix, represented as:
⎡ Nu ⎤ Nu ⎢ ωˆ x (1) ⎥ =⎣ ⎦ = T1−1 · q(2) ω(c) ˆ ω ˆ (1)
−s(i − 1)q3 (i) s(i − 1)q2 (i) −s(i − 1)q1 (i) s(i − 1)q0 (i)
|ω(i)|
(3)
The estimated angular velocity ω(1) ˆ could be generated through the attitude quaternion vector q(2) and the T1−1 , which was the inversion of the exchange matrix T1 .
2 2 2 |ω(i ˆ − 1)| = ω ˆ x (i − 1) + ω ˆ y (i − 1) + ω ˆ z (i − 1) ⎪ ⎪ ⎪ ⎪ |ω(i ˆ − 1)| ⎪ ⎪ c(i − 1) = cos ⎪ 2 ⎪ ⎪ ⎪ ⎪ sin(|ω(i ˆ − 1)|/2) ⎪ ⎪ s(i − 1) = ⎪ |ω(i ˆ − 1)| ⎨ c(i − 1)q0 (i) −s(i − 1)q1 (i) −s(i − 1)q2 (i) ⎪ Ti = c(i − 1)q1 (i) s(i − 1)q0 (i) −s(i − 1)q3 (i) ⎪ ⎪ c(i − 1)q2 (i) s(i − 1)q3 (i) s(i − 1)q0 (i) ⎪ ⎪ ⎪ c(i − 1)q3 (i) −s(i − 1)q2 (i) s(i − 1)q1 (i) ⎪ x ⎪ ⎪ ⎪ ⎪ ω ˆ x (i) ⎪ ⎪ = Ti−1 · q(i + 1) + d(i) ⎪ ⎩ ωˆ y (i)
417
(10)
where M was the total measurement point number. Based on the model error minimization criterion, minimized the model error cost function J, we could get the representation of model error as follow: d(i) =
T
[S[ω(i)]] ˆ R−1 [S[ω(i)]] ˆ +W T
−1
·[S[ω(i)]] ˆ R−1 [Z[ω(i)] ˆ − e(i) − b(i)]
(19)
Imported the model error d(i) back in to the system state function to modify the system error state vectoreω (i + 1) and continue the regression process. The system error state vectoreω (i + 1) represented the error component of FOG.
418
H. Tian et al. / Optik 122 (2011) 415–421
At second step, a drift trend extraction process was deployed to filter out the noise and short time interference from the separated error component of FOG. This extraction calculation was based on the empirical model decomposition (EMD) [21]. The target of drift trend extraction was to decompose an original gyroscope output signal x(i) as a intrinsic mode component h(i) and a final residual r(i). The FOG error components with the different time step-size were represented as the different Intrinsic Mode Functions (IMFs), K h(i) = h (i), where the hk (i) was a IMFs sequence. As the EMD k=1 k algorithm was based on the time local trait, it could be used in nonlinear and non-stationary process as a sifting function [22]. For a real-valued signal, the EMD performed the mapping in equation below. x(i) = h(i) + r(i) =
K
hk (i) + r(i),
(20)
k=1
where dk (i) was a set of IMFs, and K was the number of IMF. The r(i) was the residual. The first IMF was calculated in envelope calculation block and IMFs separation block. The drift trend extraction algorithm contained three blocks, which were envelope calculation block, IMFs separation block and condition judgment block. The ANFIS interpolation unit was embedded in the envelope calculation block to approximate the maximum and minimum envelope by the limited sample data. The condition judgment block was to judge whether the output data from the IMFs separation block contained another IMF or the final residual data. The error component of FOG, which was separated before, was sent to the envelope calculation block as the input signal. The envelope calculation block responded to identify all extrema of x(i) and then Interpolate between maxima and minima of x(i). This block ended up with the maximum envelope emax (i) and minimum envelope emin (i). The envelope interpolation method in our paper was developed from ANFIS algorithm in the ANFIS interpolation unit. With the inference capability of fuzzy logic and the learning capability of neural networks, the ANFIS, which is based on the utilization of first order SUGENO model, is more flexible and has a preferable capability in approximating a linear or nonlinear function [25,26]. These advantages make it very suitable in this nonlinear envelope interpolation. Input signals of the ANFIS model were time sequence value t(i). The output signals were the ANFIS approximations of maximum envelope emax (i) and minimum envelope emin (i), defined as eˆmax (i) and eˆmin (i). Considered the approximation, the time sequence value t(i) could be divided into 4 parts on average and a bell shaped function was used as the membership function of a network fuzzy set, which was defined as follows: (x) =
1
1+|
x−ci 2b , ai |
(21)
where x was the input variable, and a, b and c were premise parameters. There were 5 layer nodes and 16 fuzzy if-then rules in the structure of ANFIS, in which 16 membership functions and 24 premise parameters were for layer 1, and 48 consequent parameters were for layer 4. There were totally 864 fitting parameters needed in the conventional ANFIS structure. Compared with the conventional cubic spline interpolation, the ANFIS approach with nonlinear multi layers structure was better in nonlinear data approximation. The ANFIS structure needed to be trained before approximation [27]. The training data was the sampled maxima data xmax (i) and minima data xmin (i) and the responded time point t(i), which could be expressed as a vector array, [xmax (i), xmin (i), t(i)]. In ANFIS training process, the time point t(i) was set as the input and the data xmax (i) and xmin (i) were treated as the output.
A mixture algorithm has been used to train parameters of ANFIS for quick parameters identification [27,28]. The training data was brought to the inputs, the premise parameters were assumed to be fixed and the optimal consequent parameters were estimated by an iterative least mean square procedure in the forward pass. In the reverse pass, the patterns were propagated again, but this time the consequent parameters were assumed to be fixed, and the back-propagation was used to modify the premise parameters [26]. The iteration training number of the ANFIS model was 1000; the minimum square estimation between the practical system output and the model output was used as a network training index. After the training, the ANFIS approximation network would be used to approximate the maximum envelope emax (i) and the minimum envelope emin (i). We treated a serial time sequence, t(i), where i was a continuous time value, as the approximation network input. The output of the ANFIS network would be the approximation of the maximum envelope, which was defined as eˆmax (i), and the minimum envelope, which was defined as eˆmin (i). The IMFs separation block, which computed the mean m(i) = (emax (i) + emin (i))/2 and Subtracted from the signal to obtain h(i) = x(i) − m(i). There were two conditions needed to be satisfied if h(t) was IMF, which was introduced in paper [29]. Ideally, h(t) should be an IMF since the construction of h(t) described above has made it satisfy all the requirements of IMF. If h(t) did not meet the definition of a true IMF, the sifting process had to be repeated many times to correct this until the first IMF imf1 was obtained. After the first IMF was separated, it was need to judge whether the output data from the IMFs separation block contains another IMF [22]. This process was completed in the condition judgment block. The first IMF, imf1 , should contain the finest scale or the shortest period component of the signal. We could separate imf1 from the rest of the data by: h(i) − imf1 (i) = r1 (i),
(22)
Since the residue, r1 (i), still contained long period components, it should be treated as the new data and subject to the same sifting process as described above. The sifting process stopped by any of the following predetermined criteria: either when the component, imfn , or the residue, rn (i), becomed so small and less than the pre-determined value of substantial consequence, or when the residue, rn (i), becomed a monotonic function and no another IMF can be extracted from it. After the decomposition, the gyroscope signal had been separated in to a monotonic residual and the IMFs components [21]. The final residual rN (i) was a zero mean noise process, which did not contain drift trend information. However, the IMFs components separated in this section contained the trend characters that had disguised in FOG error. All of these trend characters could be utilized by navigation system to calibrate or compensate the performance of FOG. As these trend characters in IMFS had been decomposed into different time step-size domains, it had become easy to select the notable component as the main drift trend of FOG.
4. Verification and simulation To verify the performance of our FOG drift extraction algorithm, software simulation and FOG testing experiments were launched. The related simulation and experiments were grouped into three types, which were angular velocity generation algorithm simulation, FOG trend extraction verification and system experiment on FOG.
H. Tian et al. / Optik 122 (2011) 415–421
419
Table 1 The deviations of the recovered angular velocities of circular orbit. Axis
Mean deviation (rad/s)
X axis Y axis Z axis
−7
2.282 × 10 2.261 × 10−7 2.248 × 10−7
Standard deviation (rad/s) 2.846 × 10−7 2.820 × 10−7 2.841 × 10−7
Table 2 The deviations of the recovered angular velocities of circular orbit using conventional algorithm.
Fig. 2. The spacecraft attitude quaternion of circular orbit.
4.1. Angular velocity generation algorithm simulation The angular velocity generation algorithm was simulated on ground. The vehicle attitude information was generated by the spacecraft orbit simulation software Satellite Tools Kit (STK), which was a powerful satellite simulation soft developed by Analytical Graphics, Inc. (AGI). The reference time format mentioned in the simulation was J2000 Reference Coordinates, which was a subtype of ECICS and defined as the X and Z axes point toward the mean vernal equinox and mean rotation axes of the Earth on January 1, 2000 at 12:00:00.00 Universal Time Coordinated (UTC). J2000.0 = 2000 January 1.5 = JD 2451545.0 Terrestrial Dynamical Time (TDT) Two individual satellite orbits were built to verify the performance of our algorithm in different orbit and attitude characteristics. 1. Circular orbit verification The first algorithm simulation was in a circular orbit condition, which was generated from the High-Precision Orbit Propagator of the STK simulation platform. The parameters of this orbit were as follows: • • • • • • •
Semi major axis: 6878.137 km Eccentricity: 1.85671◦ × 10−16 Inclination: 45◦ Period: 5676.978029 s Argument of Perigee: 0◦ Mean anomaly: 30◦ RAAN: 0◦
Axis
Mean deviation (rad/s)
Standard deviation (rad/s)
X axis Y axis Z axis
6.642 × 10−7 7.302 × 10−7 6.665 × 10−7
8.873 × 10−7 9.691 × 10−7 8.635 × 10−7
software platform, are shown in Fig. 2. The estimated angular velocities generated by our angular velocity generation algorithm are shown in Fig. 3. The deviations of the recovered angular velocities, which included the mean error and the standard deviation results, are listed in Table 1. Shown from the verification results, the attitude recovering deviations were about 1 percent of the input angular velocity range. For these deviations, the data acquisition error of STK software, which had a round-off error in small data calculation, was the most contribution. As this phenomenon might also appears in the data acquisition process of stellar sensor, an external filter should be utilized to smooth the data. As a comparison, another group of simulation was also launched using the algorithm that introduced in article [30]. The deviations of the recovered angular velocities in circular orbit using conventional algorithm are listed in Table 2. 2. Critically inclined orbit verification The second algorithm simulation was in a critically inclined orbit condition, and the parameters of this orbit are as follows: • • • • • • •
Semi major axis: 12,578.137 km Eccentricity: 0.461118◦ Inclination: 63.434949◦ Period: 14038.976969 s Argument of Perigee: 0◦ RAAN: 0◦ Mean anomaly: 329.517505◦
The spacecraft attitude vector, expressed as the quaternion T q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] , which were generated from STK
The spacecraft attitude vector, expressed as the quaternion T q(i) = [q0 (i), q1 (i), q2 (i), q3 (i)] , which were generated from STK software platform, are shown in Fig. 4.
Fig. 3. The estimated angular velocities of circular orbit.
Fig. 4. The spacecraft attitude quaternion of critically inclined orbit.
420
H. Tian et al. / Optik 122 (2011) 415–421 Table 5 FOG drift extraction results by the trend extraction algorithm. The final residuals after the NFL-EMD Angular velocity
Mean
Standard deviation
Variance
0◦ /s 1◦ /s 10/s 100◦ /s
0.3875 −3.3 −1.3373 0.4303
5.5843 0.0369 2.7572 2.2696
31.1846 1.3625 75.0559 51.5117
Table 6 FOG drift extraction results by the conventional EMD algorithm. The final residuals after the conventional EMD Fig. 5. The estimated angular velocity of critically inclined orbit.
Angular velocity ◦
The estimated angular velocities of critically inclined orbit generated by our angular velocity generation algorithm are shown in Fig. 5. The deviations of the recovered angular velocities, which included the mean error and the standard deviation results, are listed in Table 3. The deviations of the recovered angular velocities in critically inclined orbit using conventional algorithm are listed in Table 4.
0 /s 1◦ /s 10◦ /s 100◦ /s
Mean
Standard deviation
Variance
0.705 −12.1 −2.81 0.6825
14.9735 0.1375 13.3541 9.381
64.7358 12.276 35.4732 46.6724
4.2. FOG trend extraction verification To verify the performance of our FOG trend extraction algorithm, we had launched a specified verification on FOG. The FOG we used in experiment was in middle precision (with the original bias stability of 1◦ /s). After being used and frayed for a long time, the output signal of this sensor was blended with some drift components and the precision has decreased. The bias stability has become as 3◦ /s. Four test angular velocities were set as the input parameter, they were 0, 1, 10 and 100◦ /s. The output data of the FOG were sent to the trend extraction program. The ANFIS interpolation approach used in the trend extraction algorithm was the same as the architecture mentioned above, in which the membership functions number for each input was 10 and the membership function type was generalized bell membership function. The intrinsic mode data, which was extracted from the FOG error data, could represent the drift components in different time step-size. An analysis had been made to examine the statistic characteristics (Mean, Standard deviation and Variance) of the final residuals in four different input angular velocities, and the analysis results are listed in Table 5. As the comparison, we used the conventional EMD algorithm to complete another drift extraction experiment. The architecture of
Fig. 6. The import attitude quaternion of FOG.
Table 3 The deviations of the recovered rotation velocities of critically inclined orbit. Axis X axis Y axis Z axis
Mean deviation (rad/s) −7
2.337 × 10 3.632 × 10−7 2.230 × 10−7
Standard deviation (rad/s)
Fig. 7. The FOG output data.
−7
2.943 × 10 4.881 × 10−7 2.806 × 10−7
Table 4 The deviations of the recovered rotation velocities of critically inclined orbit by conventional algorithm. Axis
Mean deviation (rad/s)
Standard deviation (rad/s)
X axis Y axis Z axis
6.621 × 10−7 7.387 × 10−7 6.720 × 10−7
9.071 × 10−7 8.672 × 10−7 9.607 × 10−7
this EMD algorithm was introduced in paper [21] and the interpolation function was the cubic spline function. The analysis results of this verification data are listed in Table 6. Compared with the conventional EMD algorithm, our trend extraction algorithm based on nonlinear empirical model predictive architecture could achieve higher performance and the results had shown that the bias stability of the final residual in our drift extraction was 0.62 times as the residual in conventional EMD extraction.
H. Tian et al. / Optik 122 (2011) 415–421
Fig. 8. The trend components in three axes.
4.3. System performance testing To test the performance of our proposal in real FOG operation, FOG drift extraction verification was launched and introduced in this section. It was known that the angular velocity of orbit spacecrafts are about −1◦ /s → +1◦ /s, which is a small velocity importation. We set the geostationary orbit satellite angular velocity as the import angular velocity of the FOG. The FOG was fixed on a revolving table for simulating the rotation of the satellite. There were three axes contained in this FOG, the X axis of it pointed to the sky, the Y axis pointed to the north direction, the Z axis pointed to the east direction. The import attitude quaternion of FOG are shown in Fig. 6. The FOG output data are shown in Fig. 7. Extracting the drift trend components by our algorithm, the results of the trend components in three axes are shown in Fig. 8. 5. Conclusion An FOG drift extraction algorithm was exhibited in the paper and the related simulations were discussed. The algorithm, which could extract the drift trend from the output data of FOG, made it possible to eliminate the integration drift of gyroscope. The utilization of this approach could be on spacecraft, satellite and space-shuttle. The algorithm verification was held on ground, which had verified the performance of the algorithm and the results were described. Some modifications have been deployed to gain a more approximation of the actual angular velocity. Besides, the ground verification of algorithm was in good condition, which had not considered the external noise and potential interference so much. The experiment in space orbit environment will be hold in next step to verify the algorithm more deeply. The observation of large points, which appear at the simulation result sporadically, is a potential threat to the application and need to be eliminated. The related modification is under development currently. Acknowledgment The author gratefully acknowledges Mrs. Li Min, who has made remarkable contribution for this paper. References [1] M. Kuritsky, M. Goldstein, I. Greenwood, H. Lerman, J. McCarthy, T. Shanahan, M. Silver, J. Simpson, Inertial navigation, Proc. IEEE 71 (October (10)) (1983) 1156–1176.
421
[2] R. Zhu, Y. Zhang, Q. Bao, A novel intelligent strategy for improving measurement precision of fog, IEEE Trans. Instrum. Meas. 49 (December (6)) (2000) 1183–1188. [3] L.G. Richard, Initial navigation technology from 1970–1995, J. Inst. Navigat. 42 (January (1)) (1995) 165–185. [4] Y. Akcayir, Y. Ozkazanc, Gyroscope drift estimation analysis in land navigation systems Control Applications, in: CCA 2003. Proceedings of 2003 IEEE Conference, vol. 2, June 2003, pp. 1488–1491. [5] L. Miao, J.F. Zhang, S.J. Shen, Data analysis and modeling of fiber optic gyroscope drift, J. Beij. Inst. Technol. 11 (January (1)) (2002) 50–55. [6] J. Li, H. Wang, A. Liu, Trend extraction of gyro’s drift based on modified empirical mode decomposition, Syst. Eng. Electron. 27 (June (6)) (2005) 1080–1082. [7] R. Zhu, Y. Zhang, Q. Bao, A novel intelligent strategy for improving measurement precision of fog, IEEE Trans. Instrum. Meas. 49 (December (6)) (2000) 1183–1188. [8] H.S. Ahn, S.H. Lee, Gyroless attitude estimation of sun-pointing satellites using magnetometers, Geosci. Rem. Sens. Lett. IEEE 2 (January (1)) (2005) 8–12. [9] M. Grewal, M. Shiva, Application of kalman filtering to gyroless attitude determination and control system for environmental satellites, in: Decision and Control, 1995, Proceedings of the 34th IEEE Conference, vol. 2, December 1995. [10] J.C. Yoon, B.S. Lee, K.H. Choi, Spacecraft orbit determination using gps navigation solutions, Aerospace Sci. Technol. 4 (April (3)) (2000). [11] G.-S. Huang, J.-C. Juang, Application of nonlinear kalman filter approach in dynamic gps-based attitude determination, in: Circuits and Systems, 1997. Proceedings of the 40th Midwest Symposium, vol. 2, August 1997, pp. 1440–1444. [12] M. Grewal, M. Shiva, Application of kalman filtering to gyroless attitude determination and control system for environmental satellites, in: Decision and Control, 1995. Proceedings of the 34th IEEE Conference, vol. 2, December 1995, pp. 1544–1552. [13] M. Challa Solar, Anomalous, and Magnetospheric Particle Explorer (SAMPEX) Real Time Sequential Filter (RTSF), Evaluation Report, NASA-Goddard Space Flight Center, Greenbelt, MD, USA, April 1993. [14] J. Crassidis, F. Markley, Predictive filtering for attitude estimation without rate sensors, Journal of Guidance, Contr. Dynam. 20 (May–June (3)) (1997) 522–527. [15] D. Mook, J. Junkins, Minimum model error estimation for poorly modeled dynamic systems, J. Guidance Contr. Dynam. 3 (January–February (4)) (1988) 367–375. [16] V. Parameswaran, J. Raol, Estimation of model error for nonlinear system identification, in: Control Theory and Applications, IEE Proceedings, vol. 141, no. 6, November 1994, pp. 403–408. [17] P. Lu, Nonlinear predictive controllers for continuous systems, J. Guidance Contr. Dynam. 17 (May–June (3)) (1994) 553–560. [18] H. Wang, J. Chen, W. Huang, Y. Zhang, J. Tang, Detection and analysis of the fiber optic gyroscope’s random drift, Opt. Tech. 30 (September(5)) (2004) 623–627. [19] W.W. Dai, Z. Yang, S.L. Lim, O. Mikhailova, J. Chee, Processing and analysis of ecg signal using nonorthogonal wavelet transform, in: Engineering in Medicine and Biology Society, 1998. Proceedings of the 20th Annual International Conference of the IEEE, vol. 1, October–November 1998, pp. 139–142. [20] Y. Zheng, D.B.H. Tay, L. Li, Signal extraction and power spectrum estimation using wavelet transform scale space filtering and bayes shrinkage, Signal Process. 80 (August (8)) (2000) 1535–1549. [21] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. A: Math. Phys. Eng. Sci. 454 (March (1971)) (1998) 903–995. [22] H. Zhang, Q. Gai, Research on properties of empirical mode decomposition method, Intelligent Control and Automation, 2006. WCICA 2006. The Sixth World Congress, vol. 2, June 2006, pp. 10001–10004. [23] T. Tanaka, D.P. Mandic, Complex empirical mode decomposition, Signal Process. Lett. IEEE 14 (February (2)) (2007) 101–104. [24] A. Brown, G. Moy, Long duration strapdown stellar-inertial navigation using satellite tracking, in: Position Location and Navigation Symposium, 1992. Record. 500 Years After Columbus—Navigation Challenges of Tomorrow, IEEE PLANS’92,IEEE, March 1992, pp. 194–201. [25] C. Lee, J.-S. Wang, Self-adaptive neuro-fuzzy systems: structure and learning Intelligent Robots and Systems, (IROS 2000), in: Proceedings. 2000 IEEE/RSJ International Conference, vol. 1, 2000, pp. 52–57. [26] M. Sugeno, G.T. Kang, Structure identification of fuzzy model, Fuzzy Sets Syst. 28 (October (1)) (1988) 15–33. [27] J.-S. Jang, Anfis: adaptive-network-based fuzzy inference system, IEEE Trans. Syst. Man Cybernet. 23 (May/June (3)) (1993) 665–685. [28] T. Amaral, M. Crisostomo, V. Pires, Adaptive neuro-fuzzy inference system for modelling and control Intelligent Systems, in: Proceedings. 2002 First International IEEE Symposium, vol. 1, 2002, pp. 67–72. [29] P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, Signal Process. Lett. IEEE 11 (February (2)) (2004) 112–114. [30] X.Q. Chen, Y.H. Geng, On orbit calibration algorithm with star sensors of gyros, Syst. Eng. Electron. 27 (December (12)) (2005) 2112–2116.