Error analyses and calibration methods with accelerometers for optical angle encoders in rotational inertial navigation systems Fang Liu,* Wei Wang, Lei Wang, and Peide Feng School of Instrumentation Science and Optoelectronics Engineering, Beihang University, Beijing 100191, China *Corresponding author: [email protected] Received 2 July 2013; revised 27 September 2013; accepted 27 September 2013; posted 3 October 2013 (Doc. ID 193254); published 6 November 2013
By rotating a strapdown inertial navigation system (INS) over one or more axes, a number of error sources originating from the employed sensors cancel out during the integration process. Rotary angle accuracy has an effect on the performance of rotational INS (RINS). The application of existing calibration methods based on gyroscope measurements is restricted by the structure of the inertial measurement unit (IMU) and scale factor stability of the gyroscope. The multireadhead method has problems in miniaturization and cost. Hence, optical angle encoder calibration methods using accelerometers are proposed, on the basis of navigation error and accuracy requirement analyses for a single-axis RINS. The test results show that the accuracy of calibration methods proposed is higher than 4 arcsec (1σ). © 2013 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics; (120.4640) Optical instruments; (350.2770) Gratings. http://dx.doi.org/10.1364/AO.52.007724
1. Introduction
An inertial navigation system (INS) using laser gyros with rotation schemes has emerged since the 1980’s [1,2]. The position error for 15 days of a laser gyro single-axis rotational INS (RINS) for ships, introduced by [3] in 2010, is less than 2.5 nm. Research results on RINS based on fiber optics gyros (FOGs) with an advantage in reliability have been reported since the 2000’s [4,5]. Navigation errors caused by biases of inertial sensors are canceled out by rotational technology, with an acceptable cost increase. Two kinds of structure for single-axis RINSs have been developed: one is that the whole inertial measurement unit (IMU) rotates about the z axis [3]; another is that gyros x and y rotate about the z axis, but gyro z is stationary [6]. Compared with the former, the calibration method for the latter is more 1559-128X/13/327724-08$15.00/0 © 2013 Optical Society of America 7724
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
complicated, whereas the influence of the scale factor error of gyro z in the latter is less serious than that in the former [7]. Accordingly, the two kinds of structure can be extended to multiaxis RINS: one is that the IMU rotates about multiple axes; another is that there are two rotation units in the system, and each one has a group of inertial sensors with input axes that are mutually orthogonal [8]. Navigation performance is restrained by rotary angle accuracy for both kinds of IMU structures of RINSs. The influence degree depends on the attitude motion of the vehicles, which changes with different applications. Normally, the angle measurement sensor of a RINS is a high-accuracy optical angle encoder, with main errors coming from installation [9,10]. Several calibration and compensation methods for angle encoders have been developed, including methods using laser goniometers or gyroscopes in [9] and [11] and multireadheads in [12]. The laser goniometer or gyroscope method cannot be directly used in RINSs with the second IMU structure
mentioned above, due to a lack of gyroscope with input axis coinciding with rotational axis. Another disadvantage of this method is that the measurement accuracy relies heavily on the scale factor stability of the gyroscope. The effect of installation errors can be minimized efficiently by the multireadhead method, but it is unfavorable for the system designed for miniaturization and low cost. Three or more readheads and relative subdivision circuits are needed to eliminate installation errors including scale eccentricity, swash, and distortion [10,13]; the system volume will get larger, with increased cost as a result. This might be a problem in RINSs based on microelectromechanical systems (MEMS) [14]. Hence, new methods for optical angle encoder calibration with accelerometers in a single-axis RINS are proposed. Based on navigation error propagation equations, rotary angle error analyses and accuracy requirements are presented in this paper. The test results show that the estimation accuracy is higher than 4 ′′1σ. The remarkable advantage of the new methods for RINS introduced in this paper is that no additional sensor is needed, and it can be directly used in all kinds of RINS, including multiaxis RINS, no matter what the system structure is. The accuracy of the methods mentioned in this paper is not influenced by the constancy of accelerometer biases. The system needs preheating until the temperature is stabilized to eliminate the effect of the accelerometer bias changing with temperature in the process. Accelerometer scale factor stability has little effect on the measurement performance. In this paper, the IMU structure and navigation algorithm will be introduced first. The influence of rotary angle error on the navigation will be discussed in Section 3. Optical angle encoder error analyses will be presented in Section 4. Calibration methods, test results, and optimal methods will be discussed in Section 5. Conclusions will be made last. 2. Structure and Algorithm A.
IMU Structure
The IMU of a single-axis RINS includes a rotator and frame, shown in Fig. 1. The body frame (b frame, OXYZ) and the inertial sensing frame (s frame, Oxyz) are introduced. When the rotary angle is zero, the two coordinates coincide and the x axis indicates the projection of the input axis of gyro x onto
the plane perpendicular to the rotational axis. The z axis coincides with the rotation axis. The inertial sensors include three gyros and three accelerometers. Gyro z is fixed on the top of the frame, and the others are installed on the rotator. The biases of the gyros and accelerometers about the x and y axes cancel out by rotation. The angle encoder is at the bottom of the IMU. Typical performances of FOGs in the system include a 0.05°∕h bias in repeatability and stability, 0.002°∕h1∕2 angle random walk, and 50 ppm scale factor error. The performances of quartz accelerometers include a 60 μg bias in repeatability and a 20 ppm scale factor error. The optical angle encoder made by RENISHAW has an external diameter of 100 mm, with line count of 15744. The angular resolution reaches 0.8 ′′ by the subdivision technique. B. Navigation Algorithm
A continuous rotation scheme with reversal after one revolution is implemented. The navigation algorithm has two steps: coordinate transformation from the s frame to the b frame, and strapdown inertial navigation calculation. The transformation matrix relating the s frame to the b frame is expressed as 2
cos φs Cbs 4 sin φs 0
− sin φs cos φs 0
(1)
where φs is real rotary angle. The specific force vector (f) and angular rate vector (ω) of the body frame are written as f bib Cbs f sib ;
(2)
ωbib Cbs ωsib ;
(3)
where i represents the inertial reference frame. The step of coordinate transformation also includes scale factor error, and temperature and installation error compensation for inertial sensors. The conventional strapdown navigation algorithm is operated in the second step. 3. Influence of Rotary Angle Error on the Navigation
The output of the angle encoder is φs, and φ0s φs Δφs . By utilizing the small angle approximation in Eq. (1), the error of matrix Cbs due to Δφs is given by 2
− sin φs ΔCbs Δφs 4 cos φs 0
Fig. 1. IMU assembly of a single-axis RINS.
3 0 0 5; 1
− cos φs − sin φs 0
3 0 0 5; 0
(4)
where Δ is the error of variable . After the coordinate transformation, the calculated value of f bib and ωbib can be written as 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
7725
b b b b s b f 0b ib f ib Δf ib Cs ΔCs Cb f ib ;
(5)
b b s b b b ω0b ib ωib Δωib Cs ΔCs Cb ωib ;
(6)
Csb
where We define
Cbs T . 3
2
fX f bib 4 f Y 5 fZ
ωX and ωbib 4 ωY 5: ωZ
From Eqs. (5) and (6), it can be shown that " Δf bib
Δωbib
# −Δφs f Y Δφs f X ; 0
ΔCbs Csb f bib
ΔCbs Csb ωbib
" −Δφ ω # s Y Δφs ωX : 0
(7)
(8)
If there is acceleration or angular rate input, velocity error or misalignment angle error of the mathematical platform will be produced by Δφs. Phi angle error equations are introduced to analyze the propagation characteristics of the system [15]. If the duration of navigation is only a few minutes, the phi angle error equations will be simplified as Δϕ_ n −Cnb Δωbib ;
(9)
_ n Cn f b × Δϕn Cn Δf b ; ΔV b ib b ib
(10)
where n represents the local level navigation frame, ϕ represents the misalignment angle vector of the mathematical platform, V represents the velocity vector, Cnb represents the attitude matrix, and gn represents the specific force vector of the gravitational acceleration in the n frame. We define " Δϕ # " ΔV # E
Δϕn
ΔϕN ΔϕU
and
ΔV n
E
ΔV N : ΔV U
Then with Eqs. (4), (7) and (8), formulas (9) and (10) can be reformed as "
Δϕ_ E Δϕ_ N Δϕ_ U
#
" Δφ ω # s Y Cnb −Δφs ωX ; 0
3 # " −Δφs f Y ΔV_ E n b n n 4 ΔV_ 5 C f × Δϕ C Δφs f X : N b ib b 0 ΔV_ U
(11)
2
7726
"
ΔV_ E ΔV_ N ΔV_ U
# n
n
g × Δϕ
" −gΔϕ # N
gΔϕE 0
;
(13)
where g is the acceleration of gravity. Δϕ is a an important performance parameter of the INS. Formula (13) shows that ΔϕE and ΔϕN will cause a velocity error by the projection of gn. Δφs will cause ΔϕE and ΔϕN from Eq. (11), if the vehicle pitches or rolls. The error limit derived by Eq. (11) is
3
2
If the vehicle is stationary or moving with a uniform velocity, formula (12) will be simplified as
(12)
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
jΔϕE j;
ˆ jˆγ jg; jΔϕN j ≤ jΔφs j maxfjθj;
(14)
where θˆ is the change in the pitch angle, and γˆ is the change in the roll angle. According to Eq. (14), the rotary angle accuracy requirement depends on the range of attitude motion of the vehicle, and it changes with different applications. Given a misalignment angle accuracy of 3 ′′ (1σ) in pitching or rolling motion, Δφs must be less than 12 ′′ (1σ) for a ship or highway motor vehicle with a maximum pitch or roll angle change of 15°, and it is less than 4 ′′ (1σ) for a sport utility vehicle or large transport airplane with a maximum pitch or roll angle change of 45°. 4. Optical Angle Encoder Installation Error Analyses
Optical angle encoder installation errors can be compensated for by a mathematical model method that includes scale eccentricity, swash, and distortion. A. Scale Eccentricity
Scale eccentricity is the offset between the geometric center of the scale and the axis of rotation. The angle measurement error of the eccentricity, discussed in [9], is expressed as Δφse −
r sinφ0s K − sin K; R
(15)
where r is the value of eccentricity, R is the external diameter of the scale, and K is the initial phase. The geometric relationship of the parameters in Eq. (15) is shown in Fig. 2. O is the axis of rotation, O1 is the geometric center of the scale, B is the reference mark, K is the angle between OO1 and OB, the value of which is decided by installation. The real rotary angle is ∠BOX, and the output angle of the readhead is ∠BO1 X. The maximum of Δφse is related to r, R and K. Given R 50 mm, Fig. 3 shows examples of Δφse with different values of r and K. Noncoaxial machining of the scale mounting base will cause Δφse to reach scores or hundreds of arcseconds. The recommended installation accuracy of the eccentricity is 3 μm, under the condition of coaxial machining in the product manual [16]. But,
C.
Y Reference mark
Scale
O1
r
ϕs′
R
Equation (15) can be reformed as Δφse
ϕs
K
B
Readhead
O
Installation Error Analysis Summary
r r r sin K − cos K sin φ0s − sin K cos φ0s : R R R (18)
X
Comparing Eq. (18) with Eq. (16), all the optical angle encoder installation errors are expressed as Δφs a0
Fig. 2. Geometric relationship of scale eccentricity.
n X ak sin kφ0s bk cos kφ0s ;
(19)
k1
even so, Δφse reaches a maximum of 25 ′′ with a scale diameter of 50 mm, which is unacceptable in a RINS with a high accuracy requirement. B.
a0 −
Swash and Scale Distortion
The geometric center of the scale is not exactly inclined to the axis of rotation, referred to as swash. The main error is written in terms of cosines with amplitudes that depend on swash angle, scale diameter, and gringe width [9]. Scale distortion can be caused by inappropriate preloading of fastening screws on the scale. The Fourier analysis method will be used for this model. The expression combining the two kinds of encoder errors is Δφss as0
n X
ask sin kφ0s bsk cos kφ0s :
(16)
k1
The real output of the encoder refers to the zero position provided by a reference mark. Let us define Δφss 0 when φ0s 0. Substituting the initial condition in Eq. (16) yields as0 −
n X
with a constraint equation yield from the initial condition, which is
bsk :
(17)
k1
Some examples given by [10] show that the value of Δφss is about a few arcseconds and depends on the installation condition and the scale diameter.
n X
bk :
(20)
k1
The first-order term in Eq. (19) includes Δφse and the first-order term of Δφss, and the higher-order terms belong to Δφss . The coefficients of the terms in Eq. (19) get smaller as the order increases. 5. Calibration Methods and Test Results
Because of the noncoaxial machining of the scale mounting base, the eccentricity error might be very large in this single-axis RINS. A simple method to measure the eccentric value is utilizing dial test indicator, and the result is 48.0 μm. The range of Δφse is from 198 ′′ to 396 ′′ by Eq. (15), and the real value depends on K. Thus, the coefficients of the first-order terms in Eq. (19) are far more than the others. To illustrate the influence of higher-order terms clearly, a four position method is designed to diminish first-order angle error first. Then, a calibration method for second-order angle errors is introduced, and test results will be shown later. The parameters will be calibrated utilizing accelerometers along the x and y axes. The system needs preheating until the temperature is stabilized, to eliminate the effect of accelerometer biases changing with temperature in the process. A. Calibration of First-Order Errors
Eccentricity-induced angle error ( " )
100
The system is mounted with the z axis almost in the horizontal plane. By setting the rotary angle at intervals of 90° (four position), the acceleration along the x and y axes is measured at each position, denoted as axi and ayi i 0; 1; 2; 3, shown in Fig. 4. U represents the up direction. α is the small angle between the x axis and the horizontal plane when φ0s 0, φ0s of position i is i × 90°, and Δφs1i is the rotary angle error of position i. As Δφs10 0, the relationship between the measurements and Δφs1i is given by
r = 0.02µm, K = 0 50
r = 0.01µm, K = 0
0
-50 -100 r = 0.01µm, K = 50°
r = 0.02µm, K = 50° -150
0
100
200 Angle ( ° )
300
400
Fig. 3. Eccentricity-induced angle error with R 50 mm.
1 Δφs12 ax0 ax2 − ax1 − ax3 ; g 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
(21)
7727
U
U
U
α
X
ay1 O
O
ax0 ay0 (a) Position 0 (0 °)
(b) Position 1 (90 °) U
U
ay2
ax3 ax2
O
ay3
(α + ∆ϕ s13)
ϕs
ay
α ay
Fig. 5. Second-order error calibration with the consideration of installation errors.
Fig. 4. First-order error calibration (four position).
(22)
Ignoring second- and higher-order terms in Eq. (19), with substitution of the constraint equation gives Δφs1 −b1 a1 sin φ0s b1 cos φ0s :
ax
O y
(d) Position 3 (270 °)
1 Δφs11 − Δφs13 ay1 ay3 − ay0 − ay2 : g
αX
α ax
Y
O
(α + ∆ϕ s12)
(c) Position 2 (180° )
ax1
(α + ∆ϕ s11)
x
(23)
The expressions gained from the equation above are Δφs11 −b1 a1 ;
(24)
Δφs12 −2b1 ;
(25)
angle between the x axis and the up direction. The installation errors of the two accelerometers in the s frame need to be considered, and are denoted as αax and αay . The acceleration of gravity measured in the body frame is approximated as aX ≐ g; 28 aY ≐ − gαX : By a coordinate transformation with φ0s , the components of g measured by the X and Y axes accelerometers are given by ax g cosφ0s Δφs2 αax − gαX sinφ0s Δφs2 αax ≐ g1 − αX Δφs2 − αX αax cos φ0s − Δφs2 sin φ0s − αX αax sin φ0s ;
(29)
ay −gαX cosφ0s Δφs2 αay − g sinφ0s Δφs2 αay
Δφs13 −b1 − a1 :
(26)
1 − αX Δφs2 − αX αay sin φ0s
From the three equations above, the parameters in Eq. (23) can be written as
a1 b1
1 2 Δφs11 − − 12 Δφs12 :
Δφs13 ;
27
From Eqs. (22) and (27), the parameters of the first-order terms in Eq. (23) are gained from the accelerometer measurements. B.
Calibration of Second-Order Errors
After the first-order terms of the angle error have been compensated for, a curve fitting method is used to estimate the second-order terms of the angle error caused by swash and scale distortion. The system is mounted with the z axis almost in the horizontal plane and the x axis almost in the up direction, shown in Fig. 5. The system is mounted with the z axis almost in the horizontal plane. αX is the small 7728
≐ − gαX αay cos φ0s
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
Δφs2 cos φ0s :
(30)
Ignoring third- and higher-order terms in Eq. (19) gives Δφs2 a0 a1 sin φ0s b1 cos φ0s a2 sin 2φ0s b2 cos 2φ0s :
(31)
Substituting for Δφs2 following Eq. (31), Eqs. (29) and (30) are reformed as a a ax g − 1 1 − 2 − αX a0 − αX αax cos φ0s 2 2 b a 2 − a0 − αX − αax sin φ0s 1 cos 2φ0s 2 2 b1 a2 b2 0 0 0 sin 2φs cos 3φs − sin 3φs ; − 2 2 2
(32)
where αX ak and αX bk k 1; 2 are ignored because they are far less than ak and bk . The Fourier expansions for ax and ay are expressed as ax ax0
2nd- and 3rd-order terms of ax (m/s2)
b b ay −g 1 2 a0 αX αay cos φ0s 2 2 a2 b 1 − αX a0 − αX αay sin φ0s 1 cos 2φ0s 2 2 a1 b2 a2 0 0 0 sin 2φs cos 3φs sin 3φs ; (33) 2 2 2
x 10
8
-4
6 4 2 0 -2 -4 -6
0
100
200
300
400
Angle ( ° )
3 X
mxk sin kφ0s nxk cos kφ0s ;
(34)
3 X myk sin kφ0s nyk cos kφ0s :
(35)
Fig. 6. Test and fitting results of ax .
k1
k1
The parameters of the two equations above are gained from ax and ay curve fittings. Comparing Eqs. (32) and (33) with Eqs. (34) and (35), the parameters in Eq. (31) are given by
ak−1 1g nxk − myk k 2; 3 bk−1 − 1g mxk nyk
(36)
with the consideration of averaging multiple solutions. With Eq. (36) and the constraint Eq. (20), the constant parameter in Eq. (31) is written as 1 a0 −b1 − b2 mx2 ny2 mx3 ny3 : g C.
(37)
The measurements of the four position method mentioned above are shown in Table 1. The test is repeated five times with a 1 min sample time at each position. According to Eq. (27), the estimated result of the first-order errors is given by Δφs1 ′′ 157.7 − 110.4 sin φ0s − 157.7 cos φ0s : (38) Comparing Eq. (18) with Eq. (38), the scale eccentricity is calculated by ignoring the swash
ax0 a 0.203355 0.9 × 10−5 ay0 a −9.807021 1.5 × 10−5 a
Measurements of First-Order Error Calibration
ax1 a −9.808323 1.3 × 10−5 ay1 a −0.249413 1.2 × 10−5
q r R a21 b21 :
(39)
The curve fitting method mentioned above is implemented by the motor pausing for 10 s at intervals of 2° in a revolution. The test is repeated three times. Curve fitting results of second- and thirdorder terms of ax and ay are shown in Figs. 6 and 7. The dot, fork, and cross shapes represent three tests of second- and third-order terms of ax or ay respectively, and the real line is the fitting curve. According to Eqs. (31) and (36), the estimated result is given by Δφs2 ′′ 4.9 3.3 sin φ0s 4.2 cos φ0s
Calibration Results
Table 1.
and distortion, using Eq. (38), and the value is 46.7 μm. The difference between this value and that measured by a dial test indicator is 1.3 μm, which mainly results from ignoring the swash and distortion. It is more evidently indicated that the eccentricity error is far more than that of the swash and distortion in this single-axis RINS,
ax2 a −0.205914 1.1 × 10−5 ay2 a 9.790733 1.4 × 10−5
The unit of measure here is m∕s2 .
ax3 a 9.790776 1.4 × 10−5 ay3 a 0.222634 1.0 × 10−5
− 6.4 sin 2φ0s − 9.1 cos 2φ0s :
(40)
According to Eq. (31), the accuracy of Δφs2 is given by 2nd- and 3rd-order terms of ay (m/s2)
ay ay0
8
x 10
-4
6 4 2 0 -2 -4 -6
0
100
200
300
400
Angle ( ° )
Fig. 7. Test and fitting results of ay . 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
7729
12
E. Optimal Calibration Method
Angle error ( " )
8 4 0 -4 -8 -12
0
100
200
300
400
Angle ( ° )
Fig. 8. Rotary angle measurement accuracy.
σΔφs2 ≤ σa1 2σb1 σa2 2σb2
By rotating a strapdown inertial navigation system (INS) over one or more axes, a number of error sources originating from the employed sensors cancel out during the integration process. Rotary angle accuracy has an effect on the performance of rotational INS (RINS). The application of existing calibration methods based on gyroscope measurements is restricted by the structure of the inertial measurement unit (IMU) and scale factor stability of the gyroscope. The multireadhead method has problems in miniaturization and cost. Hence, optical angle encoder calibration methods using accelerometers are proposed, on the basis of navigation error and accuracy requirement analyses for a single-axis RINS. The test results show that the accuracy of calibration methods proposed is higher than 4 arcsec (1σ). © 2013 Optical Society of America OCIS codes: (000.3860) Mathematical methods in physics; (120.4640) Optical instruments; (350.2770) Gratings. http://dx.doi.org/10.1364/AO.52.007724
1. Introduction
An inertial navigation system (INS) using laser gyros with rotation schemes has emerged since the 1980’s [1,2]. The position error for 15 days of a laser gyro single-axis rotational INS (RINS) for ships, introduced by [3] in 2010, is less than 2.5 nm. Research results on RINS based on fiber optics gyros (FOGs) with an advantage in reliability have been reported since the 2000’s [4,5]. Navigation errors caused by biases of inertial sensors are canceled out by rotational technology, with an acceptable cost increase. Two kinds of structure for single-axis RINSs have been developed: one is that the whole inertial measurement unit (IMU) rotates about the z axis [3]; another is that gyros x and y rotate about the z axis, but gyro z is stationary [6]. Compared with the former, the calibration method for the latter is more 1559-128X/13/327724-08$15.00/0 © 2013 Optical Society of America 7724
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
complicated, whereas the influence of the scale factor error of gyro z in the latter is less serious than that in the former [7]. Accordingly, the two kinds of structure can be extended to multiaxis RINS: one is that the IMU rotates about multiple axes; another is that there are two rotation units in the system, and each one has a group of inertial sensors with input axes that are mutually orthogonal [8]. Navigation performance is restrained by rotary angle accuracy for both kinds of IMU structures of RINSs. The influence degree depends on the attitude motion of the vehicles, which changes with different applications. Normally, the angle measurement sensor of a RINS is a high-accuracy optical angle encoder, with main errors coming from installation [9,10]. Several calibration and compensation methods for angle encoders have been developed, including methods using laser goniometers or gyroscopes in [9] and [11] and multireadheads in [12]. The laser goniometer or gyroscope method cannot be directly used in RINSs with the second IMU structure
mentioned above, due to a lack of gyroscope with input axis coinciding with rotational axis. Another disadvantage of this method is that the measurement accuracy relies heavily on the scale factor stability of the gyroscope. The effect of installation errors can be minimized efficiently by the multireadhead method, but it is unfavorable for the system designed for miniaturization and low cost. Three or more readheads and relative subdivision circuits are needed to eliminate installation errors including scale eccentricity, swash, and distortion [10,13]; the system volume will get larger, with increased cost as a result. This might be a problem in RINSs based on microelectromechanical systems (MEMS) [14]. Hence, new methods for optical angle encoder calibration with accelerometers in a single-axis RINS are proposed. Based on navigation error propagation equations, rotary angle error analyses and accuracy requirements are presented in this paper. The test results show that the estimation accuracy is higher than 4 ′′1σ. The remarkable advantage of the new methods for RINS introduced in this paper is that no additional sensor is needed, and it can be directly used in all kinds of RINS, including multiaxis RINS, no matter what the system structure is. The accuracy of the methods mentioned in this paper is not influenced by the constancy of accelerometer biases. The system needs preheating until the temperature is stabilized to eliminate the effect of the accelerometer bias changing with temperature in the process. Accelerometer scale factor stability has little effect on the measurement performance. In this paper, the IMU structure and navigation algorithm will be introduced first. The influence of rotary angle error on the navigation will be discussed in Section 3. Optical angle encoder error analyses will be presented in Section 4. Calibration methods, test results, and optimal methods will be discussed in Section 5. Conclusions will be made last. 2. Structure and Algorithm A.
IMU Structure
The IMU of a single-axis RINS includes a rotator and frame, shown in Fig. 1. The body frame (b frame, OXYZ) and the inertial sensing frame (s frame, Oxyz) are introduced. When the rotary angle is zero, the two coordinates coincide and the x axis indicates the projection of the input axis of gyro x onto
the plane perpendicular to the rotational axis. The z axis coincides with the rotation axis. The inertial sensors include three gyros and three accelerometers. Gyro z is fixed on the top of the frame, and the others are installed on the rotator. The biases of the gyros and accelerometers about the x and y axes cancel out by rotation. The angle encoder is at the bottom of the IMU. Typical performances of FOGs in the system include a 0.05°∕h bias in repeatability and stability, 0.002°∕h1∕2 angle random walk, and 50 ppm scale factor error. The performances of quartz accelerometers include a 60 μg bias in repeatability and a 20 ppm scale factor error. The optical angle encoder made by RENISHAW has an external diameter of 100 mm, with line count of 15744. The angular resolution reaches 0.8 ′′ by the subdivision technique. B. Navigation Algorithm
A continuous rotation scheme with reversal after one revolution is implemented. The navigation algorithm has two steps: coordinate transformation from the s frame to the b frame, and strapdown inertial navigation calculation. The transformation matrix relating the s frame to the b frame is expressed as 2
cos φs Cbs 4 sin φs 0
− sin φs cos φs 0
(1)
where φs is real rotary angle. The specific force vector (f) and angular rate vector (ω) of the body frame are written as f bib Cbs f sib ;
(2)
ωbib Cbs ωsib ;
(3)
where i represents the inertial reference frame. The step of coordinate transformation also includes scale factor error, and temperature and installation error compensation for inertial sensors. The conventional strapdown navigation algorithm is operated in the second step. 3. Influence of Rotary Angle Error on the Navigation
The output of the angle encoder is φs, and φ0s φs Δφs . By utilizing the small angle approximation in Eq. (1), the error of matrix Cbs due to Δφs is given by 2
− sin φs ΔCbs Δφs 4 cos φs 0
Fig. 1. IMU assembly of a single-axis RINS.
3 0 0 5; 1
− cos φs − sin φs 0
3 0 0 5; 0
(4)
where Δ is the error of variable . After the coordinate transformation, the calculated value of f bib and ωbib can be written as 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
7725
b b b b s b f 0b ib f ib Δf ib Cs ΔCs Cb f ib ;
(5)
b b s b b b ω0b ib ωib Δωib Cs ΔCs Cb ωib ;
(6)
Csb
where We define
Cbs T . 3
2
fX f bib 4 f Y 5 fZ
ωX and ωbib 4 ωY 5: ωZ
From Eqs. (5) and (6), it can be shown that " Δf bib
Δωbib
# −Δφs f Y Δφs f X ; 0
ΔCbs Csb f bib
ΔCbs Csb ωbib
" −Δφ ω # s Y Δφs ωX : 0
(7)
(8)
If there is acceleration or angular rate input, velocity error or misalignment angle error of the mathematical platform will be produced by Δφs. Phi angle error equations are introduced to analyze the propagation characteristics of the system [15]. If the duration of navigation is only a few minutes, the phi angle error equations will be simplified as Δϕ_ n −Cnb Δωbib ;
(9)
_ n Cn f b × Δϕn Cn Δf b ; ΔV b ib b ib
(10)
where n represents the local level navigation frame, ϕ represents the misalignment angle vector of the mathematical platform, V represents the velocity vector, Cnb represents the attitude matrix, and gn represents the specific force vector of the gravitational acceleration in the n frame. We define " Δϕ # " ΔV # E
Δϕn
ΔϕN ΔϕU
and
ΔV n
E
ΔV N : ΔV U
Then with Eqs. (4), (7) and (8), formulas (9) and (10) can be reformed as "
Δϕ_ E Δϕ_ N Δϕ_ U
#
" Δφ ω # s Y Cnb −Δφs ωX ; 0
3 # " −Δφs f Y ΔV_ E n b n n 4 ΔV_ 5 C f × Δϕ C Δφs f X : N b ib b 0 ΔV_ U
(11)
2
7726
"
ΔV_ E ΔV_ N ΔV_ U
# n
n
g × Δϕ
" −gΔϕ # N
gΔϕE 0
;
(13)
where g is the acceleration of gravity. Δϕ is a an important performance parameter of the INS. Formula (13) shows that ΔϕE and ΔϕN will cause a velocity error by the projection of gn. Δφs will cause ΔϕE and ΔϕN from Eq. (11), if the vehicle pitches or rolls. The error limit derived by Eq. (11) is
3
2
If the vehicle is stationary or moving with a uniform velocity, formula (12) will be simplified as
(12)
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
jΔϕE j;
ˆ jˆγ jg; jΔϕN j ≤ jΔφs j maxfjθj;
(14)
where θˆ is the change in the pitch angle, and γˆ is the change in the roll angle. According to Eq. (14), the rotary angle accuracy requirement depends on the range of attitude motion of the vehicle, and it changes with different applications. Given a misalignment angle accuracy of 3 ′′ (1σ) in pitching or rolling motion, Δφs must be less than 12 ′′ (1σ) for a ship or highway motor vehicle with a maximum pitch or roll angle change of 15°, and it is less than 4 ′′ (1σ) for a sport utility vehicle or large transport airplane with a maximum pitch or roll angle change of 45°. 4. Optical Angle Encoder Installation Error Analyses
Optical angle encoder installation errors can be compensated for by a mathematical model method that includes scale eccentricity, swash, and distortion. A. Scale Eccentricity
Scale eccentricity is the offset between the geometric center of the scale and the axis of rotation. The angle measurement error of the eccentricity, discussed in [9], is expressed as Δφse −
r sinφ0s K − sin K; R
(15)
where r is the value of eccentricity, R is the external diameter of the scale, and K is the initial phase. The geometric relationship of the parameters in Eq. (15) is shown in Fig. 2. O is the axis of rotation, O1 is the geometric center of the scale, B is the reference mark, K is the angle between OO1 and OB, the value of which is decided by installation. The real rotary angle is ∠BOX, and the output angle of the readhead is ∠BO1 X. The maximum of Δφse is related to r, R and K. Given R 50 mm, Fig. 3 shows examples of Δφse with different values of r and K. Noncoaxial machining of the scale mounting base will cause Δφse to reach scores or hundreds of arcseconds. The recommended installation accuracy of the eccentricity is 3 μm, under the condition of coaxial machining in the product manual [16]. But,
C.
Y Reference mark
Scale
O1
r
ϕs′
R
Equation (15) can be reformed as Δφse
ϕs
K
B
Readhead
O
Installation Error Analysis Summary
r r r sin K − cos K sin φ0s − sin K cos φ0s : R R R (18)
X
Comparing Eq. (18) with Eq. (16), all the optical angle encoder installation errors are expressed as Δφs a0
Fig. 2. Geometric relationship of scale eccentricity.
n X ak sin kφ0s bk cos kφ0s ;
(19)
k1
even so, Δφse reaches a maximum of 25 ′′ with a scale diameter of 50 mm, which is unacceptable in a RINS with a high accuracy requirement. B.
a0 −
Swash and Scale Distortion
The geometric center of the scale is not exactly inclined to the axis of rotation, referred to as swash. The main error is written in terms of cosines with amplitudes that depend on swash angle, scale diameter, and gringe width [9]. Scale distortion can be caused by inappropriate preloading of fastening screws on the scale. The Fourier analysis method will be used for this model. The expression combining the two kinds of encoder errors is Δφss as0
n X
ask sin kφ0s bsk cos kφ0s :
(16)
k1
The real output of the encoder refers to the zero position provided by a reference mark. Let us define Δφss 0 when φ0s 0. Substituting the initial condition in Eq. (16) yields as0 −
n X
with a constraint equation yield from the initial condition, which is
bsk :
(17)
k1
Some examples given by [10] show that the value of Δφss is about a few arcseconds and depends on the installation condition and the scale diameter.
n X
bk :
(20)
k1
The first-order term in Eq. (19) includes Δφse and the first-order term of Δφss, and the higher-order terms belong to Δφss . The coefficients of the terms in Eq. (19) get smaller as the order increases. 5. Calibration Methods and Test Results
Because of the noncoaxial machining of the scale mounting base, the eccentricity error might be very large in this single-axis RINS. A simple method to measure the eccentric value is utilizing dial test indicator, and the result is 48.0 μm. The range of Δφse is from 198 ′′ to 396 ′′ by Eq. (15), and the real value depends on K. Thus, the coefficients of the first-order terms in Eq. (19) are far more than the others. To illustrate the influence of higher-order terms clearly, a four position method is designed to diminish first-order angle error first. Then, a calibration method for second-order angle errors is introduced, and test results will be shown later. The parameters will be calibrated utilizing accelerometers along the x and y axes. The system needs preheating until the temperature is stabilized, to eliminate the effect of accelerometer biases changing with temperature in the process. A. Calibration of First-Order Errors
Eccentricity-induced angle error ( " )
100
The system is mounted with the z axis almost in the horizontal plane. By setting the rotary angle at intervals of 90° (four position), the acceleration along the x and y axes is measured at each position, denoted as axi and ayi i 0; 1; 2; 3, shown in Fig. 4. U represents the up direction. α is the small angle between the x axis and the horizontal plane when φ0s 0, φ0s of position i is i × 90°, and Δφs1i is the rotary angle error of position i. As Δφs10 0, the relationship between the measurements and Δφs1i is given by
r = 0.02µm, K = 0 50
r = 0.01µm, K = 0
0
-50 -100 r = 0.01µm, K = 50°
r = 0.02µm, K = 50° -150
0
100
200 Angle ( ° )
300
400
Fig. 3. Eccentricity-induced angle error with R 50 mm.
1 Δφs12 ax0 ax2 − ax1 − ax3 ; g 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
(21)
7727
U
U
U
α
X
ay1 O
O
ax0 ay0 (a) Position 0 (0 °)
(b) Position 1 (90 °) U
U
ay2
ax3 ax2
O
ay3
(α + ∆ϕ s13)
ϕs
ay
α ay
Fig. 5. Second-order error calibration with the consideration of installation errors.
Fig. 4. First-order error calibration (four position).
(22)
Ignoring second- and higher-order terms in Eq. (19), with substitution of the constraint equation gives Δφs1 −b1 a1 sin φ0s b1 cos φ0s :
ax
O y
(d) Position 3 (270 °)
1 Δφs11 − Δφs13 ay1 ay3 − ay0 − ay2 : g
αX
α ax
Y
O
(α + ∆ϕ s12)
(c) Position 2 (180° )
ax1
(α + ∆ϕ s11)
x
(23)
The expressions gained from the equation above are Δφs11 −b1 a1 ;
(24)
Δφs12 −2b1 ;
(25)
angle between the x axis and the up direction. The installation errors of the two accelerometers in the s frame need to be considered, and are denoted as αax and αay . The acceleration of gravity measured in the body frame is approximated as aX ≐ g; 28 aY ≐ − gαX : By a coordinate transformation with φ0s , the components of g measured by the X and Y axes accelerometers are given by ax g cosφ0s Δφs2 αax − gαX sinφ0s Δφs2 αax ≐ g1 − αX Δφs2 − αX αax cos φ0s − Δφs2 sin φ0s − αX αax sin φ0s ;
(29)
ay −gαX cosφ0s Δφs2 αay − g sinφ0s Δφs2 αay
Δφs13 −b1 − a1 :
(26)
1 − αX Δφs2 − αX αay sin φ0s
From the three equations above, the parameters in Eq. (23) can be written as
a1 b1
1 2 Δφs11 − − 12 Δφs12 :
Δφs13 ;
27
From Eqs. (22) and (27), the parameters of the first-order terms in Eq. (23) are gained from the accelerometer measurements. B.
Calibration of Second-Order Errors
After the first-order terms of the angle error have been compensated for, a curve fitting method is used to estimate the second-order terms of the angle error caused by swash and scale distortion. The system is mounted with the z axis almost in the horizontal plane and the x axis almost in the up direction, shown in Fig. 5. The system is mounted with the z axis almost in the horizontal plane. αX is the small 7728
≐ − gαX αay cos φ0s
APPLIED OPTICS / Vol. 52, No. 32 / 10 November 2013
Δφs2 cos φ0s :
(30)
Ignoring third- and higher-order terms in Eq. (19) gives Δφs2 a0 a1 sin φ0s b1 cos φ0s a2 sin 2φ0s b2 cos 2φ0s :
(31)
Substituting for Δφs2 following Eq. (31), Eqs. (29) and (30) are reformed as a a ax g − 1 1 − 2 − αX a0 − αX αax cos φ0s 2 2 b a 2 − a0 − αX − αax sin φ0s 1 cos 2φ0s 2 2 b1 a2 b2 0 0 0 sin 2φs cos 3φs − sin 3φs ; − 2 2 2
(32)
where αX ak and αX bk k 1; 2 are ignored because they are far less than ak and bk . The Fourier expansions for ax and ay are expressed as ax ax0
2nd- and 3rd-order terms of ax (m/s2)
b b ay −g 1 2 a0 αX αay cos φ0s 2 2 a2 b 1 − αX a0 − αX αay sin φ0s 1 cos 2φ0s 2 2 a1 b2 a2 0 0 0 sin 2φs cos 3φs sin 3φs ; (33) 2 2 2
x 10
8
-4
6 4 2 0 -2 -4 -6
0
100
200
300
400
Angle ( ° )
3 X
mxk sin kφ0s nxk cos kφ0s ;
(34)
3 X myk sin kφ0s nyk cos kφ0s :
(35)
Fig. 6. Test and fitting results of ax .
k1
k1
The parameters of the two equations above are gained from ax and ay curve fittings. Comparing Eqs. (32) and (33) with Eqs. (34) and (35), the parameters in Eq. (31) are given by
ak−1 1g nxk − myk k 2; 3 bk−1 − 1g mxk nyk
(36)
with the consideration of averaging multiple solutions. With Eq. (36) and the constraint Eq. (20), the constant parameter in Eq. (31) is written as 1 a0 −b1 − b2 mx2 ny2 mx3 ny3 : g C.
(37)
The measurements of the four position method mentioned above are shown in Table 1. The test is repeated five times with a 1 min sample time at each position. According to Eq. (27), the estimated result of the first-order errors is given by Δφs1 ′′ 157.7 − 110.4 sin φ0s − 157.7 cos φ0s : (38) Comparing Eq. (18) with Eq. (38), the scale eccentricity is calculated by ignoring the swash
ax0 a 0.203355 0.9 × 10−5 ay0 a −9.807021 1.5 × 10−5 a
Measurements of First-Order Error Calibration
ax1 a −9.808323 1.3 × 10−5 ay1 a −0.249413 1.2 × 10−5
q r R a21 b21 :
(39)
The curve fitting method mentioned above is implemented by the motor pausing for 10 s at intervals of 2° in a revolution. The test is repeated three times. Curve fitting results of second- and thirdorder terms of ax and ay are shown in Figs. 6 and 7. The dot, fork, and cross shapes represent three tests of second- and third-order terms of ax or ay respectively, and the real line is the fitting curve. According to Eqs. (31) and (36), the estimated result is given by Δφs2 ′′ 4.9 3.3 sin φ0s 4.2 cos φ0s
Calibration Results
Table 1.
and distortion, using Eq. (38), and the value is 46.7 μm. The difference between this value and that measured by a dial test indicator is 1.3 μm, which mainly results from ignoring the swash and distortion. It is more evidently indicated that the eccentricity error is far more than that of the swash and distortion in this single-axis RINS,
ax2 a −0.205914 1.1 × 10−5 ay2 a 9.790733 1.4 × 10−5
The unit of measure here is m∕s2 .
ax3 a 9.790776 1.4 × 10−5 ay3 a 0.222634 1.0 × 10−5
− 6.4 sin 2φ0s − 9.1 cos 2φ0s :
(40)
According to Eq. (31), the accuracy of Δφs2 is given by 2nd- and 3rd-order terms of ay (m/s2)
ay ay0
8
x 10
-4
6 4 2 0 -2 -4 -6
0
100
200
300
400
Angle ( ° )
Fig. 7. Test and fitting results of ay . 10 November 2013 / Vol. 52, No. 32 / APPLIED OPTICS
7729
12
E. Optimal Calibration Method
Angle error ( " )
8 4 0 -4 -8 -12
0
100
200
300
400
Angle ( ° )
Fig. 8. Rotary angle measurement accuracy.
σΔφs2 ≤ σa1 2σb1 σa2 2σb2
