A Fast Weakly-Coupled Double-Layer ESKF Attitude Estimation Algorithm and Application

: Aimed at the problem of small unmanned aerial vehicle (UAV) attitude solution accuracy and real-time performance in short-range navigation ﬂight, in this paper, we propose a fast weakly-coupled double-layer error-state Kalman ﬁlter (DL-ESKF) attitude estimation algorithm. Considering the application of short-range navigation, we designed an improved attitude error model for low-cost gyroscope/accelerometer/magnetometer devices. In addition, we reasonably simpliﬁed certain factors that affect the attitude solution to reduce the ﬁltering calculation burden. For the data coupling phenomenon caused by the different sampling frequencies of the attitude sensor data in the ﬁltering process, we designed a new attitude algorithm combined with the ESKF and hierarchical ﬁlter. The ﬁrst layer of ﬁlters used an accelerometer and the second layer used a magnetometer to correct the attitude error. We also built an ofﬂine and real-time test platform to verify the performance of the proposed algorithm in a simulation and ﬂight test environment compared with the classic attitude algorithms. The experimental results demonstrated that the proposed algorithm not only improved the attitude solution accuracy and stability but also reduced the ﬁlter running time.


Introduction
In the autonomous navigation and control of small unmanned aerial vehicles (UAVs) [1], accurate and reliable attitude estimation [2,3] is indispensable. In recent years, the microelectromechanical systems inertial measurement unit (MEMS IMU) [4,5] has been developed rapidly. Due to the advantages [6,7] of being lightweight, small size, low power consumption, and low cost, the MEMS IMU has become an excellent selection for small UAVs in attitude estimation. The UAV attitude estimation system [8,9] is composed of a gyroscope, accelerometer, magnetometer, and microprocessor, and uses information fusion algorithms to estimate the optimal attitude vector, which includes the roll angle, pitch angle, and yaw angle. This can still output steady attitude information without the data from the other sensors.
There are two main fusion algorithms to fuse the attitude sensors data-constant gain observer (complementary filter) and Kalman filter. The complementary filter is a multisensor fusion algorithm that integrates the gyroscope angular rate, corrects the gyroscope bias [10], and then uses the accelerometer measurement value to correct the roll angle and pitch angle and the magnetometer measurement value to correct the yaw angle [11]. The state variable of the complementary filter can usually be defined in the form of an attitude vector, such as a Euler angle or quaternion and gyro bias [12].
The complementary filter is a pretty simple algorithm, and this has very little computation burden and is easy to implement in the application. These advantages make it more suitable for low-cost embedded systems. Euston M [13] proposed a complementary filter attitude algorithm, which combined low-cost IMU [14,15] and a magnetometer using the first-order dynamic model to provide good attitude estimation for the UAV. Wu J [16] proposed a fast complementary filter fusion algorithm based on quaternion attitude estimation and designed a step-by-step filter to eliminate the effect of magnetometer outliers to improve the solution accuracy and real-time performance.
Another attitude fusion algorithm is the Kalman filter [17,18] that can model and optimize the attitude sensors system to estimate the UAV attitude. The Kalman filter can be divided into two parts-time update and measurement update when performing attitude solutions [19]. The angular rate of the gyroscope is used to calculate the attitude prediction value at a faster data sample frequency in the time update section. The acceleration of the accelerometer and the magnetic field strength of the magnetometer correct the attitude prediction at a different data sample frequency in the measurement update section and, meanwhile, constantly adjust the attitude error [20] through the Kalman filter gain.
Song Y [21] proposed a quaternion attitude Kalman filter algorithm and used the adaptive filter to modify the measurement noise covariance matrix. This solved the problem of large errors of MEMS IMU [22,23] and reduced the impact of gyroscope random errors in attitude estimation. Sabatelli S [24] proposed a quaternion double-level Kalman filter of the attitude estimation with a 9-axis MEMS IMU. The first level used the accelerometer measured value, and the second level used the magnetometer measured value, which can reduce the coupling phenomenon between the sensors measured values, and make the filtering solution more flexible and safe.
Unlike the complementary filter, the feedback gain of the Kalman filter is real-time and always reflects the state variable changes of the system. At the same time, the Kalman filter can generate the optimal estimated state for the linear system. The corresponding state covariance changes indicate the uncertainty of the estimated state variable, which is very important for analyzing the system stability. The Kalman filter also considers the statistical characteristics of the process noise and observation noise when modeling the system.
In this paper, we studied small UAV low-cost attitude sensors based on the analysis and summary of common attitude fusion filtering algorithms. For the problem of the solution accuracy and attitude real-time performance of UAV, we proposed a new method combined with the Kalman filter and hierarchical filter, namely the weakly-coupled double-layer error-state Kalman filter (DL-ESKF) attitude estimation algorithm. The proposed algorithm was performed and verified in simulation and flight tests. Compared with other classic attitude estimation algorithms, this algorithm has the following advantages.

1.
Considering the small UAV short-range autonomous navigation flight environment, we designed a low-cost attitude IMU error model, which reasonably simplified the impact of the cone error and earth radius on the attitude solution accuracy.

2.
For the data coupling phenomenon caused by the inconsistency of measurement sensor update frequency during the attitude filtering, a weakly-coupled double-layer ESKF attitude filtering algorithm was presented to improve the attitude solution accuracy and robustness.

3.
To verify the effectiveness of the proposed algorithm, we analyzed the simulation and flight test compared with other attitude algorithms. It was more comprehensive to study the performance of the proposed algorithm through the mathematical statistics methods [25] and the running time of the algorithm.
The remainder of this paper is organized as follows-the mathematical attitude model is introduced and designed in Section 2. The weakly-coupled DL-ESKF attitude estimation is proposed in Section 3. Then, the experimental simulation and flight test in Section 4 demonstrates the performance of the proposed algorithm. Finally, the conclusion and future work are given in Section 5.

Quaternion Attitude Determination Model
In general, UAV strapdown inertial attitude representation methods include the Euler angle, DCM, and quaternion. However, the Euler angle has a "Gimbal Lock" problem when the pitch angle is 90 degrees, which can occur with the wrong attitude solution. The DCM is a matrix containing three DOF, and its matrix multiplication and inversion operations can make the real-time performance decrease due to the matrix computation burden. Compared with the previous two methods, the quaternion consists of the four elements, which can quickly calculate the rotation of the rigid carrier and the attitude transformation [19].
where Q Q Q = [q 0 , q 1 , q 2 , q 3 ] T is the quaternion vector. The q 0 and the [q 1 , q 2 , q 3 ] T are the real part and the imaginary part of the quaternion, respectively. The [i i i, j j j, k k k] T is the unit vector and orthogonal of each. The quaternion continuous multiplication operation can be defined by Equation (2).
where ⊗ is the quaternion multiplication, and the quaternion must satisfy the norm principle.
In this paper, the quaternion attitude differential equation is established to update the UAV attitude.
where the Q Q Q * = q 0 − q 1 i i i − q 2 j j j − q 3 k k k is the conjugate complex operation of the Q Q Q, ω ω ω = [ω x , ω y , ω z ] T is three-axis angular rate of the gyroscope, and the Ω Ω Ω ω can be represented as follows.
The first-order Picard iteration method [21] is used to solve the quaternion attitude differential equation, and the discrete model is established by Equation (6).
where the Q Q Q(t k ) and Q Q Q(t k−1 ) denote the unit quaternion attitude at time t k and t k−1 , respectively. The q q q(∆t k )(∆t k = t k − t k−1 ) can be derived by Equation (7).
The UAV attitude is usually composed of the roll angle φ(t k ), pitch angle θ(t k ), and yaw angle ψ(t k ). The the updated quaternion needs to transform the attitude by the Equation (9).
where the range of the roll and yaw angle are the (−π, π), the range of the pitch angle is the (− π 2 , π 2 ).

Attitude Sensors Model
Low-cost attitude sensors mainly include the IMU, which consists of a gyroscope, accelerometer, and magnetometer as shown in Figure 1. The IMU is designed of the size 4 cm × 3 cm and consists of the microcontroller unit (MCU) , a three-axis magnetometer (Mag), and a six-axis IMU used for the attitude sensor data sampling and processing. In general, the three axis of the attitude sensors are orthogonal to each other and fixed on the center of gravity of the UAV in ideal conditions.
whereω ω ω(t k ),â a a(t k ), andm m m(t k ) are the measured three-axis angular rate, three-axis acceleration, and three-axis magnetic field strength value in the body frame [26], respectively. We assumed calibration [27] before using the attitude sensors, and the scale factors k ω , k a , and k m are all 1.
denote the attitude sensor bias, but they change slowly over time. The n n n ω (t k ), n n n a (t k ), and n n n m (t k ) are the attitude sensor noise and are assumed to be the zero-mean Gaussian white noise. The C C C b n (t k ) is the rotation matrix from the navigation frame to the body frame, and g g g is the gravity vector in the navigation frame.

Attitude Filtering System Model
The different attitude system models have different selections to the attitude state estimation parameters, and this is because the kinds of state estimation parameters can affect the establishment and design of the system model. There are two kinds of system models-direct state Kalman filter (DSKF) and ESKF. However, compared to DSKF, the ESKF has the following advantages [28]. On the one hand, the magnitude of the state estimation value is different during the filtering due to the selection of the error state. On the other hand, the error state filtering method can make a large estimated error and the initial error value converge to the steady attitude value. If the corresponding parameter settings are adjusted properly, it will be more robust than the DSKF. Given these, this paper presents an error state quaternion attitude update model based on the Kalman filter algorithm.
The main principle of ESKF is to treat the real state vectorX X X as a combination of the nominal state vector X X X and error state vector δ X X X; the mathematical representation between them is expressed as follows: whereX X X ∈ 3 denotes the state vector, and ⊕ is the typical addition operation. For the quaternion, the mathematical representation can be expressed by Equation (13).
In this paper, the true state vectorX X X and the nominal state vector X X X are both 4×1 unit quaternion vectors, and the error state vector δ X X X is the 6×1 vector. In addition, the b b b ω is the gyroscope random bias in the body frame. where Given the above analysis, the continuous Gaussian error-state system equation is established by Equations (15) and (16).
where δX X X is the error-state vector, f f f (t) is the state transition matrix function, g g g(t) ∼ N(0, q q q(t)) is the process noise. δZ Z Z(t) is the measurement vector, h h h(t) is the measurement matrix function, v v v(t) ∼ N(0, r r r(t)) is the measurement noise. g g g(t) and v v v(t) are assumed to be the zero-mean Gaussian white noise and uncorrelated each other. The δ(t) is the Dirac function.

Attitude Error Dynamic Equation
When the UAV is static, zero-bias is the main factor affecting the gyroscope's output accuracy. Eliminating the zero-bias can improve the three-axis angular rate measurement accuracy of the gyroscope. However, when the UAV is flying, the gyroscope random drift error will constantly change over time and, therefore, must be designed and estimated [30].
where the b b b 0 is the gyroscope zero-bias, the b b b r is the gyroscope random bias, and the w w w b is the gyroscope error white noise. This is often a constant after the gyroscope is started, and thus Equation (19) can be derived as follows.
Based on the strapdown inertial attitude error vector equation [30] and the gyroscope error model [31], the attitude error dynamic model can be established and derived as follows.
where the δϕ ϕ ϕ is the attitude error in the body frame and C C C n b is the rotation matrix from the body frame to the navigation frame. The ω ω ω n in is the angular rate vector of the navigation frame relative to the inertial frame, and b b b n ω is the gyroscope bias in the navigation frame. The τ g is a first-order Markov correlated time constant, and n n n ω is the gyroscope measured value noise.
In this paper, we mainly used low-cost MEMS IMU. Compared with high-precision inertial devices [5], the low-cost sensor noise is large. In addition, small UAVs often perform tasks in the short-range navigation scenarios. The flight distance is generally several kilometers to dozens of kilometers. The impact of cone error and the earth rotation angular rate can be approximately ignored for a short flight distance. Based on these conditions, this paper reasonably established and simplified the attitude error dynamic equation.

Attitude Error Measurement Equation
The error state was selected as the filtering parameters estimated in this paper. According to the system equation of Equation (15), the corresponding measurement vector was the error quantity. The measurement sensor was the low-cost accelerometer and magnetometer shown in Figure 1, and the sensor sampling rate was different. Given practical application considerations, it was necessary to establish an attitude error measurement equation.
where the φ a (t k ) and θ a (t k ) are the roll angle and pitch angle of the accelerometer calculated by Equation (24), the ψ m (t k ) is the yaw angle of magnetometer calculated by Equation (25). Theφ(t k ), θ(t k ), andψ(t k ) are the attitude estimation of quaternion updating solution. The v v v a (t k ) is the measurement noise of the roll angle error and pitch angle error, v v v m (t k ) is the measurement noise of the yaw angle error.
where the a a a(t k ) = [a x (t k ), a y (t k ), a z (t k )] T is the three-axis acceleration in the body frame, m m m(t k ) = [m x (t k ), m y (t k ), m z (t k )] T is the magnetic field strength in the body frame.

Attitude Filter Algorithm
The attitude angle error and gyro bias were used as the state estimation of the Kalman filter. Based on the short-range UAV application scenarios, the error state dynamic equation and error measurement equation were derived and simplified reasonably. However, the sampling frequency of the IMU and the magnetometer measurement values were different in actual applications. If all the measured values were filtered and updated uniformly, not only do they lose the data of the high sampling frequency but also the coupling problem occurs between different sensor data resulting in filtering divergence. The weakly-coupled hierarchical filtering method [22] can update the sensor data of different sampling frequencies, which solves the coupling problem between the sensor measured data and makes the filtering more flexible and safe.
We combined the ESKF and weakly-coupled hierarchical filter, and used two different sampling frequency measurement sensors. The flow-chart of computations performed by the proposed DL-ESKF attitude algorithm is summarized in Figure 2. This method allowed the sensor values to enter the corresponding filtering measurement update to correct the attitude error when the different data update rates were triggered and could reduce the mutual interference between the sensor data. The weakly-coupled double-layer ESKF algorithm is as follows: (1) attitude filtering time update δ X X X k/k−1 and P P P k/k−1 to the attitude error dynamic equation: where δ X X X = [δφ ϕ ϕ, b b b ω ] T is the attitude error state vector, P P P is the error state covariance with the principle of symmetric positive matrix. Φ Φ Φ is the state transform matrix and can be derived by the discretization of the state transition matrix funtion f f f (t k ) and obtained by Equations (27)- (29). Q Q Q is the process noise covariance matrix by Equations (30) and (31).
where the dt is the gyroscope sampling update time.
The the first-layer attitude filtering measurement update δ X X X a k1 and P P P a k1 to the attitude error measurement Equation (22): where the δZ Z Z a k1 is the measurement error value about the roll angle and pitch angle, R R R a k1 is the measurement noise covariance matrix about the roll angle and pitch angle. K K K a k1 is the first-layer filtering gain matrix, and H H H a k1 is the measurement matrix by Equation (33).
The second-layer attitude filtering measurement update δ X X X m k2 and P P P m k2 to the attitude error measurement Equation (23) where the δZ Z Z m k2 is the measurement error value about the yaw angle, R R R m k2 is the measurement noise covariance matrix about the yaw angle. K K K m k2 is the second-layer filtering gain matrix, and H H H m k2 is the measurement matrix by Equation (35). The pseudo code of the proposed attitude algorithm (Algorithm 1) is given to elaborate the principle in more detail as follows.

Algorithm 1 Proposed weakly-coupled double-layer ESKF.
1: Input Input Input: three-axis angular rate ω ω ω, acceleration a a a, and magnetic field strength m m m; 2: Initialization Initialization Initialization:X X X t 0 , X X X t 0 , δ X X X t 0 , P P P t 0 , Q Q Q, R R R a , R R R m , τ g ; 3: Output Output Output: UAVs attitude (roll φ, pitch θ, and yaw ψ); Filter Update Filter Update Filter Update: 4: while t k > 0(k = 1, 2, ..., ∞) do 5: t k+1 = t k + dt; 6: Update gyroscope data ω ω ω t k+1 ; 7: Gyroscope quaternion attitude update Equation (8); 8: Transform quaternion attitude into Euler angle Equation (9); time update time update time update: 9: Calculate the error state and covariance prediction matrix by the Equation (26); First-layer measurement update First-layer measurement update First-layer measurement update: 10: if The acceleration is obtained (t k < t a a a t k < t k+1 ) then 11: Update accelerometer data a a a t k ; 12: Calculate the roll angle and pitch measurement angle by the Equation (24); 13: Obtain the roll angle and pitch measurement angle error by the Equation (22); 14: Use K K K a k1 to continuously correct the error state; 15: update the error state δ X X X a k1 and error covariance matrix P P P a k1 by the Equation (32); 16: end if Second-layer measurement update Second-layer measurement update Second-layer measurement update: 17: if The manetic field strength is obtained (t k < t m m m t k < t k+1 ) then 18: Update magnetometer data m m m t k ; 19: Calculate the yaw measurement angle by the Equation (25); 20: Obtain the yaw measurement angle error by the Equation (23); 21: Use K K K m k2 to continuously correct the error state; 22: update the error state δ X X X m k2 and error covariance matrix P P P m k2 by the Equation (34); 23: end if Correct UAV Attitude Correct UAV Attitude Correct UAV Attitude: 24: Perform quaternion attitude correction by the Equation (13), and quaternion normalization (3); 25: Transform corrected quaternion attitude into Euler angle Equation (9); the error state resets the error state resets the error state resets 26: The update estimated error state is setted 0; 27: end while

Simulation and Flight Verification Platform
The algorithm simulation of the attitude sensor data was collected using the commercial autopilot platform [32] of Figure 3, and the corresponding UAV specifications are listed in Table 1. The sensors are the IMU, magnetometer, barometer, and GPS module, and the detailed parameters are given in the Table 2.  The self-designed flight control autopilot was developed as shown in Figure 4, and was mainly used to verify the applied performance of the proposed algorithm for the UAV and its specifications are the same as Table 1. The core processing unit of the autopilot was an STM32F407VGT6 microprocessor with a 168 MHz clock frequency. The sensors are the IMU (MPU6500, in, magnetometer (HMC5883), barometer, and GPS module. In addition, the performance regarding the self-designed flight control autopilot is listed in the Table 3.

Static Conditions
The experimental environment in which the platform was placed was an outdoor flat place to collect the attitude sensor data under static conditions. Figure 5-7 is the attitude angle error under static conditions. The attitude error curve of the Mahony-CF is volatile, and the error is also larger than the Song-KF and the DL-ESKF. The reasons are as follows-(1) The UAV body vibration caused by the rotation of the blade can be transmitted to the attitude sensors, which affects the accuracy of the measurement data under the static condition. (2) This paper's research object was the low-cost sensors of the small UAV attitude estimation, and the low-cost sensors have large accuracy errors compared with high-precision IMU sensors. (3) The Song-CF only used the complementary frequency characteristics of the sensors to calculate the attitude, and did not consider the influence of the sensors noise to the attitude solution accuracy in the attitude system modeling process, which damaged the reliability and stability of the calculation.   In addition, the solution accuracy of the Song-KF was better than Mahony-CF. However, the Song-KF constructed the different sensors data update frequency as one measurement equation for the accelerometer and magnetometer. This method occurred as the coupling phenomenon between the acceleration and magnetic field strength, enlarged the dimension of the measurement matrix, and decreased the attitude solution accuracy. Therefore, the proposed algorithm not only considered the impact of low-cost sensors noise on the solution accuracy to the system model but also included a weakly-coupled double-layer measurement update for the problem of different update frequencies of the accelerometer and the magnetometer. The attitude angle error of the proposed algorithm was lower than the previous two methods under static conditions. Figures 8-10 show the attitude angle under static conditions. We regarded the attitude angle of the autopilot [32] as the truth, which was obtained from the EKF attitude estimator. We found that the attitude estimated values of the proposed DL-ESKF algorithm were closer to the truth. The stability and convergence of the proposed algorithm during the 5000 s static conditions were better than the other two methods.
To more intuitively compare the attitude solution accuracy and running time computing performance of the three filtering algorithms, Tables 7 and 8 provide the RMSE and running time of the algorithm. The proposed algorithm was smaller than the other two algorithms under static conditions in terms of the attitude RMES calculation. The running time of the Machony-CF was the smallest, and this is because it does not calculate the state covariance matrix compared with the KF and, thus, reduces the computation burden. The proposed algorithm was faster than the Song-KF, and slightly slower with the Mahony-CF; this is attributed to the simplification of the model and design of algorithm architecture.

Dynamic Condition
According to pre-setting waypoints and UAV flight path as shown in Figure 11, the UAV was tested by the mode of the autonomous positioning and navigation, and, at the same time, the attitude sensor data were collected. Figures 12-14 show the attitude angle error under dynamic conditions. The Mahony-CF and Song-KF have a larger attitude solution error than the DL-ESKF. This is because the body jitter and wind speed increase the uncertainty of the attitude sensor noise value during the flight, and the Mahony-CF cannot suppress the sensor noise. This disadvantage amplifies the effect on the attitude solution accuracy and increases the influence of the attitude sensor data on each other when the UAV performs certain maneuvers, and damages the attitude solution accuracy of the Song-KF. The proposed algorithm considers these factors so that the attitude angle error under dynamic conditions is smaller than the other two algorithms, and the solved curve is relatively smoother.     Particularly for the yaw angle, the difference between the DL-ESKF and the truth was the smallest. Table 9 denotes that the attitude RMSE changed under dynamic conditions; the RMSE of the proposed algorithm was smaller than the other two algorithms. Table 10 shows that the proposed algorithm was faster than the Song-CF, and closer to the Mahony-CF in terms of the running time performance, mainly because the weakly-coupled double-layer measurement update method reduced the coupling problem between the attitude sensor data. This also decreased the dimension of the measurement matrix, and made the multiplication and reversible operations of the matrix quicker.

Flight Real-Time Test
The effectiveness of the proposed algorithm in this paper was verified by simulation, and then the actual flight verification was performed on the self-developed flight control hardware as shown in Figure 18. The attitude changes of the UAV were sent to the ground station in real-time through digital transmission to observe the correctness of the solution values. There are many parts of UAVs in the flight control system (FCS), and these include sensors of perception, remote control, ground station, and PWM rotors. The UAV FCS sets a task list of different calling frequencies, as shown in Figure 19, to ensure steady operation.

Dynamic Test
If the multi-rotor UAV maneuvering amplitude is too large, this will cause a loss of control during the flight so that the dynamics of the proposed algorithm cannot be tested. Therefore, the UAV is first stationary in the direction of the south horizontal plane, then it is held for a large-scale operation, and finally, it is stationary in the direction of the south horizontal plane.
It can be seen from Figure 20 that after a large-scale operation, and then placing the UAV in the south horizontal plane, the roll angle, pitch angle, and yaw angle can quickly converge to the initial attitude values. We demonstrated that the proposed algorithm can meet the needs of UAV flying at a large-scale operation.

Autonomous Flight Test
The UAV will fly autonomously after setting the waypoint and flight speed on the UAV ground station. The flight path is shown in Figure 21.
In this paper, the multirotor UAV was used in the autonomous flight. The flight speed was slower than the fixed-wing UAV, and the changes in attitude were also smaller than the fixed-wing UAV. The changes in attitude solution were relatively stable in the Figure 22. The reason why the yaw angle continued changing from 200 s to 600 s in Figure 23 is that when the UAV was flying from one waypoint to another, it adjusted the nose direction of the UAV. At this time, the calculated yaw angle value constantly changed in real time.

Conclusions and Future Work
In this paper, we presented a new weakly-coupled double-layer error Kalman filter attitude estimation algorithm to provide the reliable and steady attitude solution during small UAV real-time flight.
For the small UAV short-range navigation application, a low-cost gyroscope/accelerometer/magnetometer attitude error mathematical model was designed. A double-layer ESKF attitude algorithm was proposed to solve the coupling phenomenon of the attitude sensor data due to the different sampling frequencies. The proposed DL-ESKF attitude algorithm was verified in simulation compared with the other attitude algorithms under the static and dynamic conditions of the UAV. In addition, we tested on large-scale operations and autonomous flight, which provided good attitude solution values. Future work will study the impact of harmful acceleration and wind speed on the attitude solution accuracy when the flight velocity is faster for the fixed-wing UAV, to design a more robust and reliable attitude algorithm.

Conflicts of Interest:
The authors declare no conflict of interest.