Reference / Paper · 2010

An Efficient Orientation Filter for Inertial and Inertial/Magnetic Sensor Arrays

The 30 April 2010 internal report that introduced the Madgwick gradient-descent orientation filter.

Author: Sebastian O. H. Madgwick
Year: 2010
Type: Reference / Paper
Credit: Sebastian Madgwick. Author's copy via x-io Technologies (x-io.co.uk)

orientation filter
imu
quaternion
madgwick filter

An efficient orientation filter for inertial and inertial/magnetic sensor arrays Sebastian O.H. Madgwick April 30, 2010 Abstract This report presents a novel orientation filter applicable to IMUs consisting of tri-axis gyroscopes and accelerometers, and MARG sensor arrays that also include tri-axis magnetometers. The MARG implementation incorporates magnetic distortion and gyroscope bias drift compensation. The filter uses a quaternion representation, allowing accelerometer and magnetometer data to be used in an analytically derived and optimised gradient-descent algorithm to compute the direction of the gyroscope measurement error as a quaternion derivative. The benefits of the filter include: (1) computationally inexpensive; requiring 109 (IMU) or 277 (MARG) scalar arithmetic operations each filter update, (2) effective at low sampling rates; e.g. 10 Hz, and (3) contains 1 (IMU) or 2 (MARG) adjustable parameters defined by observable system characteristics. Performance was evaluated empirically using a commercially available orientation sensor and reference measurements of orientation obtained using an optical measurement system. A simple calibration method is presented for the use of the optical measurement equipment in this application. Performance was also benchmarked against the propriety Kalman-based algorithm of orientation sensor. Results indicate the filter achieves levels of accuracy exceeding that of the Kalman-based algorithm; < 0.6◦ static RMS error, < 0.8◦ dynamic RMS error. The implications of the low computational load and ability to operate at low sampling rates open new opportunities for the use of IMU and MARG sensor arrays in real-time applications of limited power or processing resources or applications that demand extremely high sampling rates. 1 Contents 1 Introduction 3 2 Quaternion representation 4 3 Filter derivation 3.1 Orientation from angular rate . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Orientation from vector observations . . . . . . . . . . . . . . . . . . . . . . 3.3 Filter fusion algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Magnetic distortion compensation . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Gyroscope bias drift compensation . . . . . . . . . . . . . . . . . . . . . . . 3.6 Filter gains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 6 6 10 11 12 13 4 Experimentation 4.1 Equipment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Orientation from optical measurements . . . . . . . . . . . . . . . . . . . . . 4.3 Calibration of frame alignments . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Experimental proceedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 14 14 16 18 5 Results 5.1 Typical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Static and dynamic performance . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Filter gain vs. performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Sampling rate vs. performance . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Gyroscope bias drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 19 21 22 23 23 6 Discussion 24 7 Conclusions 25 A IMU filter implementation optimised in C 29 B MARG filter implementation optimised in C 30 2 1 Introduction The accurate measurement of orientation plays a critical role in a range of fields including: aerospace [1, 2, 3], robotics [4, 5], navigation [6, 7] and human motion analysis [8, 9] and machine interaction [10]. Whilst a variety of technologies enable the measurement of orientation, inertial based sensory systems have the advantage of being completely self contained such that the measurement entity is constrained neither in motion nor to any specific environment or location. An IMU (Inertial Measurement Unit) consists of gyroscopes and accelerometers enabling the tracking of rotational and translational movements. In order to measure in three dimensions, tri-axis sensors consisting of 3 mutually orthogonal sensitive axes are required. A MARG (Magnetic, Angular Rate, and Gravity) sensor is a hybrid IMU which incorporates a tri-axis magnetometer. An IMU alone can only measure an attitude relative to the direction of gravity which is sufficient for many applications [4, 2, 8, 1]. MARG systems, also known as AHRS (Attitude and Heading Reference Systems) are able to provide a complete measurement of orientation relative to the direction of gravity and the earth’s magnetic field. A gyroscope measures angular velocity which, if initial conditions are known, may be integrated over time to compute the sensor’s orientation [11, 12]. Precision gyroscopes, ring laser for example, are too expensive and bulky for most applications and so less accurate MEMS (Micro Electrical Mechanical System) devices are used in a majority of applications [13]. The integration of gyroscope measurement errors will lead to an accumulating error in the calculated orientation. Therefore, gyroscopes alone cannot provide an absolute measurement of orientation. An accelerometer and magnetometer will measure the earth’s gravitational and magnetic fields respectively and so provide an absolute reference of orientation. However, they are likely to be subject to high levels of noise; for example, accelerations due to motion will corrupt measured direction of gravity. The task of an orientation filter is to compute a single estimate of orientation through the optimal fusion of gyroscope, accelerometer and magnetometer measurements. The Kalman filter [14] has become the accepted basis for the majority of orientation filter algorithms [4, 15, 16, 17] and commercial inertial orientation sensors; xsens [18], micro-strain [19], VectorNav [20], Intersense [21], PNI [22] and Crossbow [23] all produce systems founded on its use. The widespread use of Kalman-based solutions are a testament to their accuracy and effectiveness, however, they have a number of disadvantages. They can be complicated to implement which is reflected by the numerous solutions seen in the subject literature [3, 4, 15, 16, 17, 24, 25, 26, 27, 28, 29, 30, 31, 32]. The linear regression iterations, fundamental to the Kalman process, demand sampling rates far exceeding the subject bandwidth; for example, a sampling rate between 512 Hz [18] and 30 kHz [19] may be used for a human motion caption application. The state relationships describing rotational kinematics in three-dimensions typically require large state vectors and an extended Kalman filter implementation [4, 17, 24] to linearise the problem. These challenges demand a large computational load for implementation of Kalmanbased solutions and provide a clear motivation for alternative approaches. Many previous approaches to address these issues have implemented either fuzzy processing [2, 5] or fixed filters [33] to favour accelerometer measurements of orientation at low angular velocities and the integrated gyroscope measurements at high angular velocities. Such an approach is simple 3 but may only be effective under limited operating conditions. Bachman et al [34] proposed an alternative approach where the filter achieves an optimal fusion of measurements data at all angular velocities. However, the process requires a least squares regression, which also brings in an associated computational load. Mahony et al [35] developed the complementary filter which is shown to be an efficient and effective solution; however, performance is only validated for an IMU. This report introduces novel orientation filter that is applicable to both IMUs and MARG sensor arrays addressing issues of computational load and parameter tuning associated with Kalman-based approaches. The filter employs a quaternion representation of orientation (as in: [34, 17, 24, 30, 32]) to describe the coupled nature of orientations in three-dimensions and is not subject to the problematic singularities associated with an Euler angle representation1 . A complete derivation and empirical evaluation of the new filter is presented. Its performance is benchmarked against an existing commercial filter and verified with optical measurement system. Innovative aspects of the proposed filter include: a single adjustable parameter defined by observable systems characteristics; an analytically derived and optimised gradientdescent algorithm enabling performance at low sampling rates; an on-line magnetic distortion compensation algorithm; and gyroscope bias drift compensation. 2 Quaternion representation A quaternion is a four-dimensional complex number that can be used to represent the orientation of a ridged body or coordinate frame in three-dimensional space. An arbitrary orientation of frame B relative to frame A can be achieved through a rotation of angle θ around an axis A r̂ defined in frame A. This is represented graphically in figure 1 where the mutually orthogonal unit vectors x̂A , ŷA and ẑA , and x̂B , ŷB and ẑB define the principle axis of coordinate frames A and B respectively. The quaternion describing this orientation, A B q̂, is defined by equation (1) where rx , ry and rz define the components of the unit vector A r̂ in the x, y and z axes of frame A respectively. A notation system of leading super-scripts and sub-scripts adopted from Craig [37] is used to denote the relative frames of orientations and vectors. A leading sub-script denotes the frame being described and a leading super-script denotes the frame this is with reference to. For example, A B q̂ describes the orientation of A frame B relative to frame A and r̂ is a vector described in frame A. Quaternion arithmetic often requires that a quaternion describing an orientation is first normalised. It is therefore conventional for all quaternions describing an orientation to be of unit length. A B q̂ = q1 q2 q3 q4 = cos 2θ −rx sin 2θ −ry sin 2θ −rz sin 2θ (1) The quaternion conjugate, denoted by ∗ , can be used to swap the relative frames described A by an orientation. For example, B A q̂ is the conjugate of B q̂ and describes the orientation of A frame A relative to frame B. The conjugate of B q̂ is defined by equation (2). B A ∗ B q̂ = A q̂ = 1 q1 −q2 −q3 −q4 Kuipers [36] offers a comprehensive introduction to this use of quaternions. 4 (2) ẑ B ẑ A A θ r̂ ŷ B ŷ A x̂A x̂B Figure 1: The orientation of frame B is achieved by a rotation, from alignment with frame A, of angle θ around the axis A r. The quaternion product, denoted by ⊗, can be used to define compound orientations. B A For example, for two orientations described by A B q̂ and C q̂, the compounded orientation C q̂ can be defined by equation (3). A B A C q̂ = C q̂ ⊗ B q̂ (3) a ⊗ b = a1 a2 a3 a4 ⊗ b 1 b 2 b 3 b 4  T a1 b 1 − a2 b 2 − a3 b 3 − a4 b 4 a1 b2 + a2 b1 + a3 b4 − a4 b3   =  a1 b 3 − a2 b 4 + a3 b 1 + a4 b 2  a1 b 4 + a2 b 3 − a3 b 2 + a4 b 1 (4) A A ∗ v=A B q̂ ⊗ v ⊗ B q̂ (5) For two quaternions, a and b, the quaternion product can be determined using the Hamilton rule and defined as equation (4). A quaternion product is not commutative; that is, a ⊗ b 6= b ⊗ a. A three dimensional vector can be rotated by a quaternion using the relationship described in equation (5) [36]. A v and B v are the same vector described in frame A and frame B respectively where each vector contains a 0 inserted as the first element to make them 4 element row vectors. B A The orientation described by A B q̂ can be represented as the rotation matrix B R defined by equation (6) [36].   2 2q1 − 1 + 2q22 2(q2 q3 + q1 q4 ) 2(q2 q4 − q1 q3 ) A 2 2  2(q3 q4 + q1 q2 ) (6) B R = 2(q2 q3 − q1 q4 ) 2q1 − 1 + 2q3 2 2 2(q2 q4 + q1 q3 ) 2(q3 q4 − q1 q2 ) 2q1 − 1 + 2q4 Euler angles ψ, θ and φ in the so called aerospace sequence [36] describe an orientation of frame B achieved by the sequential rotations, from alignment with frame A, of ψ around 5 ẑB , θ around ŷB , and φ around x̂B . This Euler angle representation of A B q̂ is defined by equations (7), (8) and (9). ψ = Atan2 2q2 q3 − 2q1 q4 , 2q12 + 2q22 − 1 θ = −sin−1 (2q2 q4 + 2q1 q3 ) φ = Atan2 2q3 q4 − 2q1 q2 , 2q12 + 2q42 − 1 3 Filter derivation 3.1 Orientation from angular rate (7) (8) (9) A tri-axis gyroscope will measure the angular rate about the x, y and z axes of the senor frame, termed ωx , ωy and ωz respectively. If these parameters (in rads−1 ) are arranged into the vector S ω defined by equation (10), the quaternion derivative describing the rate of change of orientation of the earth frame relative to the sensor frame SE q̇ can be calculated [38] as equation (11). S ω = 0 ωx ωy ωz (10) 1S S S E q̇ = E q̂ ⊗ ω (11) 2 The orientation of the earth frame relative to the sensor frame at time t, E S qω,t , can S be computed by numerically integrating the quaternion derivative E q̇ω,t as described by equations (12) and (13) provided that initial conditions are known. In these equations, S ωt is the angular rate measured at time t, ∆t is the sampling period and SE q̂est,t-1 is the previous estimate of orientation. The sub-script ω indicates that the quaternion is calculated from angular rates. 1S S S E q̇ ω,t = E q̂est,t-1 ⊗ ωt (12) S S S E qω,t = E q̂est,t-1 + E q̇ ω,t ∆t (13) 2 3.2 Orientation from vector observations A tri-axis accelerometer will measure the magnitude and direction of the field of gravity in the sensor frame compounded with linear accelerations due to motion of the sensor. Similarly, a tri-axis magnetometer will measure the magnitude and direction of the earth’s magnetic field in the sensor frame compounded with local magnetic flux and distortions. In the context of an orientation filter, it will initially be assumed that an accelerometer will measure only gravity, a magnetometer will measure only the earth’s magnetic field. 6 If the direction of an earth’s field is known in the earth frame, a measurement of the field’s direction within the sensor frame will allow an orientation of the sensor frame relative to the earth frame to be calculated. However, for any given measurement there will not be a unique sensor orientation solution, instead there will infinite solutions represented by all those orientations achieved by the rotation the true orientation around an axis parallel with the field. In some applications it may be acceptable to use an Euler angle representation allowing an incomplete solution to be found as two known Euler angles and one unknown [5]; the unknown angle being the rotation around an axis parallel with the direction of the field. A quaternion representation requires a complete solution to be found. This may be achieved through the formulation of an optimisation problem where an orientation of the sensor, SE q̂, ˆ is that which aligns a predefined reference direction of the field in the earth frame, E d, S with the measured direction of the field in the sensor frame, ŝ, using the rotation operation described by equation (5). Therefore SE q̂ may be found as the solution to (14) where equation (15) defines the objective function. The components of each vector are defined in equations (16) to (18). ˆ S ŝ) min f (SE q̂, E d, (14) ˆ S ŝ) = S q̂ ∗ ⊗ E dˆ ⊗ S q̂ − S ŝ f (SE q̂, E d, E E (15) S q̂∈<4 E S E q̂ = q1 q2 q3 q4 d = 0 dx dy dz E ˆ S ŝ = 0 sx sy sz (16) (17) (18) Many optimisation algorithms exist but the gradient descent algorithm is one of the simplest to both implement and compute. Equations (19) describes the gradient descent algorithm for n iterations resulting in an orientation estimation of SE q̂n+1 based on an ‘initial guess’ orientation SE q̂0 and a step-size µ. Equation (20) computes the gradient of the solution surface defined by the objective function and its Jacobian; simplified to the 3 row vectors defined by equations (21) and (22) respectively. S S E qk+1 = E q̂ k − µ ˆ S ŝ) ∇f (SE q̂ k , E d, , k = 0, 1, 2...n S E S ˆ ∇f ( q̂ , d, ŝ) E k ˆ S ŝ) = J T (S q̂ , E d)f ˆ (S q̂ , E d, ˆ S ŝ) ∇f (SE q̂ k , E d, E k E k  (19) 2dx ( 12 − q32 − q42 ) + 2dy (q1 q4 + q2 q3 )+ S E ˆ S f (E q̂ k , d, ŝ) = 2dx (q2 q3 − q1 q4 ) + 2dy ( 21 − q22 − q42 )+ 2dx (q1 q3 + q2 q4 ) + 2dy (q3 q4 − q1 q2 )+  2dz (q2 q4 − q1 q3 ) − sx 2dz (q1 q2 + q3 q4 ) − sy  2dz ( 21 − q22 − q32 ) − sz 7 (20) (21)  2dy q4 − 2dz q3 S E ˆ  J (E q̂k , d) = −2dx q4 + 2dz q2 2dx q3 − 2dy q2 2dy q3 + 2dz q4 2dx q3 − 4dy q2 + 2dz q1 2dx q4 − 2dy q1 − 4dz q2  −4dx q3 + 2dy q2 − 2dz q1 −4dx q4 + 2dy q1 + 2dz q2 2dx q2 + 2dz q4 −2dx q1 − 4dy q4 + 2dz q3  2dx q1 + 2dy q4 − 4dz q3 2dx q2 + 2dy q3 (22) Equations (19) to (22) describe the general form of the algorithm applicable to a field predefined in any direction. However, if the direction of the field can be assumed to only have components within 1 or 2 of the principle axis of the global coordinate frame then the equations simplify. An appropriate convention would be to assume that the direction of gravity defines the vertical, z axis as shown in equation (23). Substituting E ĝ and normalised accelerometer measurement S â for E dˆ and S ŝ respectively, in equations (21) and (22) yields equations (25) and (26). E S ĝ = 0 0 0 1 â = 0 ax ay az  2(q2 q4 − q1 q3 ) − ax fg (SE q̂, S â) =  2(q1 q2 + q3 q4 ) − ay  2( 21 − q22 − q32 ) − az   −2q3 2q4 −2q1 2q2 2q1 2q4 2q3  Jg (SE q̂) =  2q2 0 −4q2 −4q3 0  (23) (24) (25) (26) The earth’s magnetic field can be considered to have components in one horizontal axis and the vertical axis; the vertical component due to the inclination of the field which is between 65◦ and 70◦ to the horizontal in the UK [39]. This can be represented by equation (27). Substituting E b̂ and normalised magnetometer measurement S m̂ for E dˆ and S ŝ respectively, in equations (21) and (22) yields equations (29) and (30). E S b̂ = 0 bx 0 bz m̂ = 0 mx my mz   2bx (0.5 − q32 − q42 ) + 2bz (q2 q4 − q1 q3 ) − mx fb (SE q̂, E b̂, S m̂) =  2bx (q2 q3 − q1 q4 ) + 2bz (q1 q2 + q3 q4 ) − my  2bx (q1 q3 + q2 q4 ) + 2bz (0.5 − q22 − q32 ) − mz 8 (27) (28) (29)  −2bz q3 S E  Jb (E q̂, b̂) = −2bx q4 + 2bz q2 2bz q4 −4bx q3 − 2bz q1 2bx q3 + 2bz q1 2bx q2 + 2bz q4 2bx q4 − 4bz q2 2bx q1 − 4bz q3  −4bx q4 + 2bz q2 −2bx q1 + 2bz q3  2bx q2 2bx q3 (30) As has already been discussed, the measurement of gravity or the earth’s magnetic field alone will not provide a unique orientation of the sensor. To do so, the measurements and reference directions of both fields may be combined as described by equations (31) and (32). Whereas the solution surface created by the objective functions in equations (25) and (29) have a minimum defined by a line, the solution surface define by equation (31) has a minimum define by a single point, provided that bx 6= 0. fg (SE q̂, S â) S S E S (31) fg,b (E q̂, â, b̂, m̂) = fb (SE q̂, E b̂, S m̂) JgT (SE q̂) S E (32) Jg,b (E q̂, b̂) = JbT (SE q̂, E b̂) A conventional approach to optimisation would require multiple iterations of equation (19) to be computed for each new orientation and corresponding senor measurements. Efficient algorithms would also require the step-size µ to be adjusted each iteration to an optimal value; usually obtained based on the second derivative of the objective function, the Hessian. However, these requirements considerably increase the computational load of the algorithm and are not necessary in this application. It is acceptable to compute one iteration per time sample provided that the convergence rate governed by µt is equal or greater than the physical rate of change of orientation. Equation (33) calculates the estimated orientation SE q ∇,t computed at time t based on a previous estimate of orientation SE q̂ est,t-1 and the objective function gradient ∇f defined by sensor measurements S ât and S m̂t sampled at time t. The form of ∇f is chosen according to the sensors in use, as shown in equation (34). The sub-script ∇ indicates that the quaternion is calculated using the gradient descent algorithm. S S E q ∇,t = E q̂ est,t-1 − µt ∇f = ∇f k∇f k JgT (SE q̂est,t-1 )fg (SE q̂est,t-1 , S ât ) T S Jg,b (E q̂est,t-1 , E b̂)fg,b (SE q̂est,t-1 , S â, E b̂, S m̂) (33) (34) An optimal value of µt can be defined as that which ensures the convergence rate of SE q ∇,t is limited to the physical orientation rate as this avoids overshooting due an unnecessarily large step size. Therefore µt can be calculated as equation (35) where ∆t is the sampling period and SE q̇ω,t is the physical orientation rate measured by gyroscopes and α is an augmentation of µ to account for noise in accelerometer and magnetometer measurements. 9 µt = α 3.3 S E q̇ω,t ∆t, α > 1 (35) Filter fusion algorithm An estimated orientation of the sensor frame relative to the earth frame, SE qest,t , is obtained through the fusion of the orientation calculations, SE qω,t and SE q∇,t ; calculated using equations (13) and (33) respectively. The fusion of SE q̂ω,t and SE q∇,t is described by equation (36) where γt and (1 − γt ) are weights applied to each orientation calculation. S S S E qest,t = γt E q∇,t + (1 − γt )E qω,t , 0 ≤ γt ≤ 1 (36) An optimal value of γt can be defined as that which ensures the weighted divergence of SE qω is equal to the weighted convergence of SE q∇ . This is represented by equation (37) µt is the convergence rate of SE q∇ and β is the divergence rate of SE qω expressed as where ∆t the magnitude of a quaternion derivative corresponding to the gyroscope measurement error. Equation (37) can be rearranged to define γt as equation (38). (1 − γt )β = γt µt ∆t β +β ∆t γt = µt (37) (38) Equations (36) and (38) ensure the optimal fusion of SE qω,t and SE q∇,t assuming that the convergence rate of SE q ∇ governed by α is equal or greater than the physical rate of change of orientation. Therefore α has no upper bound. If α is assumed to be very large then µt , defined by equation (35), also becomes very large and the orientation filter equations simplify. A large value of µt used in equation (33) means that SE q̂est,t-1 becomes negligible and the equation can be re-written as equation (39). S E q ∇,t ≈ −µt ∇f k∇f k (39) The definition of γt in equation (38) also simplifies as the β term in the denominator becomes negligible and the equation can be rewritten as equation (40). It is possible from equation (40) to also assume that γt ≈ 0. γt ≈ β∆t µt (40) Substituting equations (13), (39) and (40) into equation (36) directly yields equation (41). It is important to note that in equation (41), γt has been substituted as both as equation (39) and 0. β∆t ∇f S S S q̇ ∆t (41) q = −µ + (1 − 0) q̂ + t E est,t E est,t-1 E ω,t µt k∇f k 10 Equation (41) can be simplified to equation (42) where SE q̇est,t is the estimated rate of change of orientation defined by equation (43) and SE q̂˙ ,t is the direction of the error of SE q̇est,t defined by equation (44). S S S E qest,t = E q̂est,t-1 + E q̇est,t ∆t (42) S S S ˙ E q̇est,t = E q̇ ω,t − β E q̂,t (43) S ˙ E q̂,t = ∇f k∇f k (44) It can be seen from equations (42) to (44) that the filter calculates the orientation SE qest by numerically integrating the estimated orientation rate SE q̇ est . The filter computes SE q̇est as the rate of change of orientation measured by the gyroscopes, SE q̇ ω , with the magnitude of the gyroscope measurement error, β, removed in the direction of the estimated error, SE q̂˙ , computed from accelerometer and magnetometer measurements. Figure 2 shows a block diagram representation of the complete orientation filter implementation for an IMU. Accelerometer S ât J Tg (SE q̂ est,t−1)f g (SE q̂ est,t−1, S ât ) ∇f k∇f k z−1 β Gyroscope S ω t Z 1S q̂ ⊗ S ωt 2 E est,t−1 S E q̇ est,t .dt q kqk S E q̂ est,t z−1 Figure 2: Block diagram representation of the complete orientation filter for an IMU implementation 3.4 Magnetic distortion compensation Measurements of the earth’s magnetic field will be distorted by the presence of ferromagnetic elements in the vicinity of the magnetometer. Investigations into the effect of magnetic distortions on an orientation sensor’s performance have shown that substantial errors may be introduced by sources including electrical appliances, metal furniture and metal structures within a buildings construction [40, 41]. Sources of interference fixed in the sensor frame, termed hard iron biases, can be removed through calibration [42, 43, 44, 45]. Sources of interference in the earth frame, termed soft iron, cause errors in the measured direction of 11 the earth’s magnetic field. Declination errors, those in the horizontal plane relative to the earth’s surface, cannot be corrected without an additional reference of heading. Inclination errors, those in the vertical plane relative to the earth’s surface, may be compensated for as the accelerometer provides an additional measurement of the sensor’s attitude. The measured direction of the earth’s magnetic field in the earth frame at time t, E ĥt , can be computed as the normalised magnetometer measurement, S m̂t , rotated by the estimated orientation of the sensor provided by the filter, SE q̂est,t-1 ; as described by equation (45). The effect of an erroneous inclination of the measured direction earth’s magnetic field, E ĥt , can be corrected if the filter’s reference direction of the earth’s magnetic field, E b̂t , is of the same inclination. This is achieved by computing E b̂t as E ĥt normalised to have only components in the earth frame x and z axes; as described by equation (46). E ∗ ĥt = 0 hx hy hz = SE q̂est,t-1 ⊗ S m̂t ⊗ SE q̂est,t -1 E p h2x + h2y 0 hz b̂t = 0 (45) (46) Compensating for magnetic distortions in this way ensures that magnetic disturbances are limited to only affect the estimated heading component of orientation. The approach also eliminates the need for the reference direction of the earth’s magnetic field to be predefined; a potential disadvantage of other orientation filter designs [17, 24]. 3.5 Gyroscope bias drift compensation The gyroscope zero bias will drift over time, with temperature and with motion. Any practical implementation of an IMU or MARG sensor array must account for this. An advantage of Kalman-based approaches is that they are able to estimate the gyroscope bias as an additional state within the system model [26, 30, 15, 24]. However, Mahony et al [35] showed that gyroscope bias drift may also be compensated for by simpler orientation filters through the integral feedback of the error in the rate of change of orientation. A similar approach will be used here. The normalised direction of the estimated error in the rate of change of orientation, SE q̂˙ , may be expressed as the angular error in each gyroscope axis using equation (47); derived as the inverse to the relationship defined in equation (11). The gyroscope bias, S ωb , is represented by the DC component of S ω and so may removed as the integral of S ω weighted by an appropriate gain, ζ. This would yield the compensated gyroscope measurements S ωc , as shown in equations (48) and (49). The first element of S ωc is always assumed to be 0. S ∗ S ˙ ω,t = 2SE q̂est,t -1 ⊗ E q̂,t S ωb,t = ζ X S (47) ω,t ∆t (48) ωc,t = S ωt − S ωb,t (49) t S The compensated gyroscope measurements, S ωc , may then be used in place of the of the gyroscope measurements, S ω, in equation (11). The magnitude of the angular error in 12 each axis, S ω is equal to a quaternion derivative of unit length. Therefore the integral gain ζ directly defines the rate of convergence of the estimated gyroscope bias, S ωb , expressed as the magnitude of a quaternion derivative. As this process requires the use of the filter estimate of a complete orientation, SE q̂est , it is only applicable to a MARG implementation of the filter. Figure 3 shows a block diagram representation of the complete filter implementation for a MARG sensor array, including the magnetic distortion and gyroscope bias drift compensation. S ∗ S S E q̂ est,t-1 ⊗ m̂t ⊗ E q̂ est,t-1 h Magnetometer S m̂t Accelerometer S ât Z 0 q h2x + h2y 0 hz i Group 1 J Tg,b (SE q̂ est,t-1, E b̂t )f g,b (SE q̂ est,t-1, S â, E b̂t , S m̂) .dt ∇f k∇f k 2SE q̂ ∗est,t-1 ⊗ SE q̂˙ ²,t ζ z -1 β Group 2 S Gyroscope ω t 1S q̂ ⊗ S ω c,t 2 E est,t-1 S E q̇ est,t Z .dt q kqk S E q̂ est,t z -1 Figure 3: Block diagram representation of the complete orientation filter for an MARG implementation including magnetic distortion (Group 1) and gyroscope drift (Group 2) compensation 3.6 Filter gains The filter gain β represents all mean zero gyroscope measurement errors, expressed as the magnitude of a quaternion derivative. The sources of error include: sensor noise, signal aliasing, quantisation errors, calibration errors, sensor miss-alignment, sensor axis nonorthogonality and frequency response characteristics. The filter gain ζ represents the rate of convergence to remove gyroscope measurement errors which are not mean zero, also expressed as the magnitude of a quaternion derivative. These errors represent the gyroscope bias. It is convenient to define β and ζ using the angular quantities ωβ and ω̇ζ respectively, where ω̃β represents the estimated mean zero gyroscope measurement error of each axis and ω̇ζ represents the estimated rate of gyroscope bias drift in each axis. Using the relationship described by equation (11), β may be defined by equation (50) where q̂ is any unit quaternion. Similarly, ζ may be described by equation (51). 13 1 β= q̂ ⊗ 0 ω̃β ω̃β ω̃β = 2 r 3˜ ζ= ω̇ζ 4 4 Experimentation 4.1 Equipment r 3 ω̃β 4 (50) (51) The filter was tested using the xsens MTx orientation sensor [18] containing 16 bit resolution tri-axis gyroscopes, accelerometers and magnetometers. The device and accompanying software offer a mode of operations where raw sensor data may be logged at a rate of 512 Hz and then post-processed to provide calibrated sensor measurements. The calibrated sensor measurements could then be processed by the proposed filter to provide the estimated orientation of the sensor. The software also incorporates a propriety Kalman-based orientation filter to provide an additional estimate orientation. As both the Kalman-based algorithm and proposed filter’s outputs could be computed using identical sensor data, the performance of each algorithm could be evaluated relative to one-another, independent of sensor performance. A Vicon system, consisting of 8 MX3+ cameras connected to an MXultranet server [46] and Nexus [47] software was used to provide reference measurement of the orientation sensor’s actual orientation. The system is an array of IR (Infrared) sensitive cameras with incorporated IR flood lights. The cameras are fixed at calibrated positions and orientations so that the measurement subject is within the field of view of multiple cameras. The Cartesian positions of IR reflective optical markers fixed to the measurement subject may then be computed in the coordinate frame of the camera array. Cameras were fixed at a height of approximately 2.5 m, evenly distributed around the perimeter of a 4 m by 4 m enclosure. Each camera was orientated to face toward the centre of the room, approximately 30◦ to 60◦ to horizontal. Experiments were conducted with the measurement subject in the centre of the room at a height of approximately 1 m. To measure the orientation of the sensor, it was fixed to an optical orientation measurement platform specifically designed for this application. The system was used to log the positions of optical markers at a rate of 120 Hz. 4.2 Orientation from optical measurements The orientation measurement platform is comprised of 3 500 mm, mutually orthogonal rods, rigidly connected at the central position along each length. Optical markers were positioned at both ends of each rod and the orientation sensor fixed to a platform at the point where the rods coincide. The platform was constructed from an aluminium central hub, carbon fibre rods and assembled using adhesives to ensure that it had no magnetic properties that may interfere with the orientation sensor’s magnetometer. Additional optical markers were placed at arbitrary but dissimilar positions along the lengths of the rods to break the rotational symmetry and aid the identification of each rod within the measurement data. Figure 4 14 shows an annotated photograph of the orientation measurement platform where C istart , C iend , C jstart , C jend , C kstart and C kend define the measured position of each marker within the camera frame. These positions may be used to define 3 mutually orthogonal unit vectors, C x̂M , C ŷM , and C ẑM within the camera frame representing the direction the x, y and z axes of the orientation measurement platform coordinate frame; as described by equations (52), (53) and (54). These vectors define the rotation matrix describing the orientation of the measurement platform in the camera frame; as shown in equation (55). C C z end j end Sensor C xstart C C C xend j start z start Figure 4: Photograph of the orientation measurement platform C C x̂M = C ŷM = iend − C istart kC iend − C istart k jend − C jstart kC jend − C jstart k (52) C kend − C kstart kC kend − C kstart k C C C C x̂ ŷ ẑ R = M M M M (53) C C ẑM = (54) (55) Due to measurement errors and tolerances in the marker frame construction, the rotation matrix defined by equation (55) cannot be considered orthogonal and so does not represent a pure rotation. Bar-Itzhack provides a method [48] where by an optimal ‘best fit’ quaternion may be extracted from an imprecise and non-orthogonal rotation matrix. The method requires the construction of the symmetric 4 by 4 matrix, K, defined by equation (56), where rmn corresponds to the element of the mth row and nth column of C M R. The optimal C quaternion M q̂ is found as the normalised Eigen vector corresponding to the maximum Eigen value of K. Equation (57) defines the optimal quaternion accounting the alternative 15 quaternion element order convention assumed by the method, where v1 , v2 , v3 and v4 define elements of the normalised Eigen vector.   r11 − r22 − r33 r21 + r12 r31 + r13 r23 − r32 1  r21 + r12 r22 − r11 − r33 r32 + r23 r31 − r13   K=  (56) r32 + r23 r33 − r11 − r22 r12 − r21  3  r31 + r13 r23 − r32 r31 − r13 r12 − r21 r11 + r22 + r33 C v v v v q̂ = (57) 4 1 2 3 M 4.3 Calibration of frame alignments In order to compare the optical measurement of the platform orientation in the camera frame, C S M q̂, and the orientation filter’s estimated orientation of the earth in the sensor frame, E q̂est , it is necessary to know the alignment of the earth frame relative to the camera frame, C E q̂, S and the alignment of the measurement platform relative to the sensor frame, M q̂. Once these quantities are found the optical measurement of the sensor frame orientation in the earth frame, SE q̂meas , may be defined by equation (58). Although the use of optical measurement equipment in this application is documentaed [26, 24, 41], little discussion is offered on the calibration of these two quantities. S C M S E q̂meas = E q̂ ⊗ C q̂ ⊗ M q̂ (58) The earth frame x and z axes are defined by the earth’s magnetic and gravitational fields respectively. Measurements of these fields in the camera frame can be used to define the alignment C E q̂. The direction of gravity was measured using a pendulum constructed from a 1 m length of cotton thread with a small weight fixed to one end. Optical markers where fixed at either end of the length of cotton and the pendulum left to come to rest. Additional optical markers were required at arbitrary fixed positions relative to the static pendulum to break the rotational symmetry of the optical marker constellation. Figure 5 shows an annotated photograph of the pendulum where C pstart and C pend define the position of the optical markers in the camera frame. The mean position of each marker over a period of time defines the direction of the pendulum in the camera frame, C p̂. This directly defines the earth frame z axis in the camera frame, C ẑE ; as shown in equation (59).   z1 C C p̄end − p̄start C C  = ẑE = z2  (59) p̂ = C k p̄end − C p̄start k z3 The direction of earth’s magnetic field was measured using a magnetic compass constructed from a 1 m carbon-fibre rod with neodymium magnets fixed to each end; polarising each end to either magnetic north or south. The compass was hung from a 1 m length of cotton thread and left to come to rest. Optical markers were fixed to either end of the rod as well as to an arbitrary position along the length, but off-set form the rod’s axis, to beak the rotational symmetry optical marker constellation. Figure 60 shows an annotated photograph of the compass where C cstart and C cend define the position of the optical markers 16 C C pend pstart Figure 5: Photograph of the pendulum and optical markers used to measure the direction of gravity in the camera frame in the camera frame. The mean position of each marker over a period of time define the direction of the compass in the camera frame, C ĉ; as shown in equation (60). C cstart C cend Figure 6: Photograph of the magnetic compass and optical markers used to measure the direction of the earth’s magnetic field in the camera frame C C c̄end − C c̄start ĉ = C k c̄end − C c̄start k (60) Due to measurement errors and the imbalance of the suspended magnetic compass, C ĉ cannot be assumed to be orthogonal to the direction of gravity defined by C p̂ and so cannot be used to directly define the earth x axis. The non-orthogonal component of C ĉ can be computed as the vector projection of C ĉ on C p̂. This can then be removed from C ĉ to define 17 the the earth x axis direction in the camera frame, C xE ; as shown in equation (61). Once normalised, this defines the earth x axis; C x̂E , as shown in equation (62). C ĉ · C p̂ C p̂ kC p̂k2   x1 C xE C  = x2  x̂E = C k xE k x3 C xE = C ĉ − (61) (62) The earth frame y axis in the camera frame, C ŷE , can be calculated as the vector that is orthogonal to both C x̂E and C ẑE and so be defined by equation (63) where signs are chosen to satisfy the relative axis direction direction convention. The alignment of the earth frame C may be defined as the rotation matrix C x̂E , C ŷE , and C ẑE . The E R, constructed from quaternion representation, C E q̂, can then be extracted form this using Bar-Itzhack’s method [48].   p ±p1 − x21 − z12 C (63) ŷE = ±p1 − x22 − z22  , C x̂E · C ŷE = 0, C ẑE · C ŷE = 0 2 2 ± 1 − x3 − z3 C C x̂E C ŷE C ẑE (64) ER = In order to find the alignment SM q̂ it is assumed that the static error of the orientation filter’s Kalman-based algorithm is mean zero. The mean algorithm’s output, SE q̄ˆKalman , was computed for the measurement platform held stationary for a period of approximate 10 C seconds. This was used with the alignment E C q̂ and optical measurement M q̂ to define the alignment of the measurement platform in the sensor frame, SM q̂, as equation (65). S C E Sˆ M q̂ = M q̂ ⊗ C q̂ ⊗ E q̄Kalman 4.4 (65) Experimental proceedure The optical measurement data and raw orientation sensor data were logged simultaneously. The raw orientation sensor data was then processed by the accompanying software to provide the calibrated sensor data and Kalman-based algorithm output. This data then synchronised the optical measurement data, with the optical measurement data interpolated to match the higher sampling rate of the orientation sensor data. The calibrated sensor data was then processed through both the IMU and MARG implementations of the proposed orientation filter and calibrated orientation measurements