Research on Gradient-Descent Extended Kalman Attitude Estimation Method for Low-Cost MARG

Aiming at the problem of the weak dynamic performance of the gradient descent method in the attitude and heading reference system, the susceptibility to the interference of accelerometers and magnetometers, and the complex calculation of the nonlinear Kalman Filter method, an extended Kalman filter suitable for a low-cost magnetic, angular rate, and gravity (MARG) sensor system is proposed. The method proposed in this paper is a combination of a two-stage gradient descent algorithm and the extended Kalman filter (GDEKF). First, the accelerometer and magnetometer are used to correct the attitude angle according to the two-stage gradient descent algorithm. The obtained attitude quaternion is combined with the gyroscope measurement value as the observation vector of EKF and the calculated attitude of the gyroscope and the bias of the gyroscope are corrected. The elimination of the bias of the gyroscope can further improve the stability of the attitude observation results. Finally, the MARG sensor system was designed for mathematical model simulation and hardware-in-the-loop simulation to verify the performance of the filter. The results show that compared with the gradient descent method, it has better anti-interference performance and dynamic performance, and better measurement accuracy than the extended Kalman filter.


Introduction
Navigation and guidance technologies are widely used in industrial and military fields such as unmanned aerial vehicles (UAV), autonomous underwater vehicles, mobile devices, and human motion tracking. Accurate attitude measurement is essential for the development of navigation and guidance technologies [1][2][3]. If a single type of sensor is used to measure the attitude of the carrier, different sensors have their own weak points [4,5]. For example, the accelerometer cannot separate the gravitational acceleration and the linear acceleration generated during the movement of the carrier, while the gyroscope is easily affected by temperature drift and noise signals, resulting in excessive accumulated errors, and the output of the magnetometer is easily interfered by magnetic materials near the sensor. In recent years, MARG (magnetic, angular rate, and gravity) sensor systems have received widespread attention [6]. Therefore, data from different sensors need to be integrated to provide an accurate position and attitude estimation [7,8].
Among the existing inertial and magnetic sensor attitude estimation methods, the most commonly used method is the complementary filtering method. Mahony et al., proposed a complementary filtering algorithm for UAV attitude calculation [9]. Then, Liang et al., applied the complementary filtering method to the combination of inertia and magnetometer for attitude calculation [10]. Calusdian et al., proposed a quaternion-based adaptive gain complementary filter [11]. This complementary filtering algorithm uses a proportional integral (PI) controller to estimate the gyroscope deviation and provide a decent attitude estimation. Madgwick et al., proposed a gradient descent method for human body motion posture tracking [12], which can reduce the influence of magnetic 3, the specific ideas of the proposed extended Kalman filter and the design metho different processes are introduced. The experimental results are discussed in Sect Section 5 is the conclusion.

Coordinate System Establishment
In order to express the attitude, the two coordinate systems shown in Figure   used to represent the carrier coordinate system { } b xyz and the navigation coordinat tem n . Among them, the carrier coordinate system adopts the front-right-down co nate system that conforms to the right-hand rule: the x-axis points to the front of th rier, the y-axis points to the right side of the carrier, and the z-axis points to the c vertically downward. The navigation coordinate system adopts the general North Down (NED) coordinate system.

Attitude Representation
In the three-dimensional space, the Euler angle representation method is used scribe the transformation from one coordinate system to another coordinate sy through three consecutive rotations around different coordinate axes. According definition of the Euler angle, rotate around z -axis, y -axis and x -axis in turn to o yaw angle ψ , pitch angle θ and roll angle φ . Three rotations can be mathemat expressed as an independent direction cosine matrix, as shown in Equation (1). The direction cosine matrix b n C that transformed the carrier coordinate system to the ence coordinate system can be derived as shown in Equation (2).
cos cos cos sin sin sin sin cos cos sin sin sin sin cos cos sin cos cos sin cos sin sin cos sin sin sin cos cos cos b z y x n However, due to the singularity problem in the attitude calculation of the method, this paper adopts the quaternion method to express the attitude of the ca

Attitude Representation
In the three-dimensional space, the Euler angle representation method is used to describe the transformation from one coordinate system to another coordinate system through three consecutive rotations around different coordinate axes. According to the definition of the Euler angle, rotate around z-axis, y-axis and x-axis in turn to obtain yaw angle ψ, pitch angle θ and roll angle φ. Three rotations can be mathematically expressed as an independent direction cosine matrix, as shown in Equation (1). Then, the direction cosine matrix C b n that transformed the carrier coordinate system to the reference coordinate system can be derived as shown in Equation (2).
cos θ cos ψ cos θ sin ψ − sin θ sin φ sin θ cos ψ − cos φ sin ψ sin φ sin θ sin ψ + cos θ cos ψ sin φ cos θ cos φ sin θ cos ψ + sin φ sin ψ cos φ sin θ sin ψ − sin φ cos ψ cos φ cos θ However, due to the singularity problem in the attitude calculation of the above method, this paper adopts the quaternion method to express the attitude of the carrier. The transformation from one coordinate system to another coordinate system can be achieved by making a single rotation angle α around a vector r defined in the reference coordinate system. The quaternion q is defined as: Unlike the complex rotation matrix derived from Equation (2), the directional cosine matrix of the quaternion q can be simply described as: The change angle of the Euler angle ψ, φ and θ of the carrier can be obtained by the direction cosine matrix decomposition [28]:

Attitude Determination
The accelerometer can determine the pitch and roll attitude of the carrier by measuring the acceleration of gravity under static conditions. According to the pitch and roll attitude information provided by the accelerometer, the magnetometer can determine the yaw attitude information of the carrier by measuring the geomagnetic field of the environment without magnetic interference, and then obtain the overall attitude information of the carrier.

Determination of Pitch Angle and Roll Angle
Ideally, the carrier in the static state is only affected by the acceleration of gravity g. The measured output of the accelerometer in the carrier coordinate system can be represented as b a = a x a y a z T . Then, substitute the direction cosine matrix of Equation (2) to obtain Moreover, the pitch angle and roll angle can be derived:

Determination of Yaw Angle
Since the accelerometer can only measure the angle of the carrier relative to the horizontal plane, in order to obtain the heading information of the carrier, it is necessary to use a magnetometer to measure the north component of the geomagnetism at the position of the carrier, assuming that the output of the geomagnetic field measured by the magnetometer in the carrier coordinate system According to the the Equations (7) and (8), pitch and roll angles derived from the direction cosine matrix in Equation (2), l h can be calculated: Then, the yaw angle can be obtained [29]: Among them, D represents the magnetic declination, the angle between the north direction of the geomagnetic field and the north direction of the navigation coordinate system, which varies with the location of the carrier.

Gradient Descent Kalman Filter Algorithm
This paper proposes an extended Kalman filter data fusion method. The overall block diagram is shown in Figure 2. The gyroscope is used as the state vector of the system state and input into the process model of the filter. Thus, the predicted value X k is obtained. The measurement model of the filter uses the gradient descent algorithm, which takes the measured values of the accelerometer and the magnetometer and gyroscope as input. The attitude quaternion is used as the observation value of the extended Kalman filter to update the predicted value calculated by the gyroscope.

Process Model
In the process model, the Kalman filter proposed in this paper selects the state vector composed of the quaternion form q ω of the three-axis angular rate output by the gyroscope and angular velocity bias ω b . The state vector is represented by Equation (11).

State Prediction
We assume that the angular velocity bias does not change very much from one sample to the next. When calculating the prediction of the attitude quaternion in discrete time, the derivative of the quaternion needs to be numerically integrated, as shown in the equation.
Micromachines 2022, 13, 1283 6 of 20 The angular rate measured by the gyroscope obeys a vector differential equation, which describes the attitude change rate as a quaternion derivative [30]: Therefore, the attitude quaternion b q ω,t of the state space equation arises as in Equation (14).
Taking the Jacobian matrix of the state space equation, the state transition matrix of Equation (15) is obtained.

System Noise
The system noise of the prediction process mainly comes from the gyroscope. The part of quaternion noise can be expressed by Equation (16): The measurement output of the gyroscope ω = [ω x ω y ω z ] T mainly includes two components: the ideal value ω = [ω x ω y ω z ] T and the drift value δω = [δω x δω y δω z ] T of the gyroscope in the sensor coordinate system, that is ω = ω + δω. Therefore, the state equation can be rewritten as: In a discrete system, extract the system noise w k from the formula as follows: where ∆T is the sampling time, v gk is the Gaussian white noise with a mean value of zero and is normally distributed, and its covariance matrix is ∑ g = δ 2 3×3 , and the covariance matrix of quaternion noise Q ω can be expressed as: The covariance matrix of the process system noise is constructed:

Measurement Model
In this paper, the observation vector is calculated by three sensors: the accelerometer and magnetometer. The attitude quaternion calculated by the accelerometer and magnetometer as the observation value of the Kalman filter system can be expressed as q ∇ = q ∇,0 q ∇,1 q ∇,2 q ∇, 3 T . This gives the measurement vector as shown in Equation (21).

Gradient Attitude Quaternion
Before inputting the data into the Kalman filter based on quaternion, it is very important to filter out the noise from the measurement process of the sensor. In this paper, the two-stage gradient descent algorithm is used. The output of the accelerometer and the magnetometer are compared with the horizontal components of the gravitational field and the geomagnetic field, respectively, to calculate the optimal attitude quaternion as the system measurement value. Figure 3 shows the schematic diagram of the gradient descent algorithm. When the carrier is in a static state, the measured values of the accelerometer and magnetometer in the carrier coordinate system are bâ and bm , respectively, which are converted into quaternion forms as follows: In an ideal state, the accelerometer output should be the same as the gravity vector nĝ in the navigation coordinate system after the direction cosine matrix conversion from the carrier coordinate system to the navigation coordinate system. In the navigation coordinate system, the gravity vector nĝ can be normalized as According to the measured value of the accelerometer and the reference value of the gravity vector, the target error function f g ( b n q, nĝ , bâ ) is constructed: In the gradient descent algorithm, in order to find the extreme value of the target error function, the Jacobian matrix J g ( b nq , nĝ , bâ ) of the error function needs to be calculated: From the error function and the corresponding Jacobian matrix, the gradient of the error function ∇ f ( b nq , nd , bâ ) can be obtained as The accelerometer alone cannot accurately measure the carrier's yaw angle attitude, because it cannot sense the rotational movement in the z axis. Therefore, a magnetometer is needed for further compensation. Suppose the reference vector nb of the magnetic field at the position of the carrier in the navigation coordinate system is: Assuming that after bm is transformed from the carrier coordinate system to the reference coordinate system through the rotation matrix, the output of the magnetometer in the navigation coordinate system nĥ is obtained as: In the navigation coordinate system, the projection of nĥ and nb on the plane xOy should be equal, so b 2 At the same time, the error function f m ( b n q, nb , bm ) and corresponding Jacobian matrix J m ( b n q, nb , bm ) can be obtained: Combining the Equations (25) and (31) with the corresponding Jacobian matrix Equations (26) and (32), the combined error function f ∇ and the corresponding Jacobian matrix J ∇ can be obtained, respectively, as below: By definition, the overall gradient of the combined error function ∇ f ∇ = J T ∇ f ∇ can be obtained from above.
Finally, the gradient descent algorithm is used to calculate the attitude quaternion. Bring it into the observation vector y k+1 of the system to obtain the measurement equation of the Kalman filter: Among them, q ∇/k+1 is the optimal pose quaternion estimated by the gradient descent method and q ∇/k is the optimal posture estimation value calculated last time by the proposed Kalman filter.
The general form of the gradient descent algorithm of the accelerometer and magnetometer is used in Equations (33)-(35). Since the gradient descent algorithm is a first-order iterative algorithm, in order to improve the calculation accuracy, the second derivative of Equation (33) can be used. However, this solution increases the calculation load of the overall system, which is not considered in practical applications. An alternative method is proposed here to find the best estimate of the step size µ, so that the algorithm convergence rate is greater than the carrier motion. For this reason, the step size µ is positively correlated with the system sampling time ∆T, the angular rate of the carrier motion measured by the gyroscope b ω, and the scale factor β [31]: Among them, β is the gain coefficient estimated according to the zero-mean measurement error of the screw instrument, and the angular velocity ω of the carrier movement can be calculated using the following equation: In Equation (36), when the initial state of the carrier is stationary or slowly moving, µ should take a relatively small initial value µ 0 . Finally, the adaptive step size µ can be given by the following equation.
Among them, µ 0 is the initial step size. The ideal values of the parameters µ 0 and β should enable the carrier to remain stable during static testing, and to keep fast tracking during dynamic testing without excessive overshoot. The determination of these two parameters is given in the experimental section.

Measurement Transfer
In order to obtain the measurements to align with the states the connection between measurements and states must be made. This is achieved by finding the non-linear measurement vector equation h(x) and its Jacobian. Therefore, the nonlinear equation is shown in Equation (39), where the gradient quaternion is calculated by substituting Equations (33) and (34) into Equation (35).
The covariance matrix H k is now built by the matrices, as shown in Equation (40). This matrix is the same as calculating the Jacobian of Equation (39).
where H ∇ is given by Equation (41).

Observation Noise
Firstly, define the vector composed of accelerometer, magnetometer and gyroscope: The measurement covariance matrix of vector can be shown in Equation (43).
where ∑ acc , ∑ mag and ∑ gyro represent the variance matrix of accelerometer, magnetometer and gyroscope, respectively. The corresponding covariance matrix of measurement system noise is constructed: where ∑ u is the measurement covariance matrix given by Equation (43), H k is the observation covariance matrix given by Equation (40) and H T k represents the transpose of H k .

Kalman Filter Design
The initial conditions for the calculation of the Kalman filter proposed in this paper are: The initial attitude angle information of the initial state vector quaternionΦ 0 can be calculated according to Equations (7)- (10). Φ 0 is the initial covariance matrix. In order to enable the stability of the filter, P 0 should be given a large positive value [32], Φ 0 = 10I 7×7 . I 7×7 is a seven-dimensional identity matrix. Next, the system state prediction is performed, and the state equation and covariance matrix are updated from time k to k + 1 time: Q k is the covariance matrix of process noise, which is calculated by Equation (20). Then, calculate the Kalman gain of the filter as follows: Among them, R k+1 is the covariance matrix of the measurement noise, which can be determined by Equation (44).
The last step is to compare the measured value at t k+1 with the predicted value of the measured value from the system model. According to the above algorithm, the predicted value is updated with the measured value to obtain an optimal estimate; the optimal estimation of the state variable at t k+1 is as follows: The covariance is: In this way, each new measurement value collected by the system can use Equations (49)-(51) to update the system state.

Hardware Design
In order to verify the effectiveness and accuracy of the proposed attitude measurement method, a measurement device containing MARG sensor is designed. The measurement device integrates four LSM9DS1 modules and four ICM-42688-P modules. LSM9DS1 has three-axis digital linear acceleration sensors, a three-axis digital angular velocity sensor and a three-axis digital magnetometer. ICM-42688-P also has a six-axis inertial measurement function, which together form the MARG sensor system. The MARG sensor system can be driven in different ways, which can fully meet the needs of the measurement device for the adaptability and stability of different environments. In addition, the measuring device integrates an STM32H753 microprocessor for data acquisition, transmission and calculation. At the same time, data collection is sent to the PC through the RS-422 communication serial port. MATLAB and HDNT Center can be used for the data analysis and posture calculation. The overall size of the measuring device is about 45 mm * 40 mm * 20 mm, which can be widely used in various environments. The block diagram of the hardware design is shown in Figures 4 and 5 and shows a picture of the designed measuring device.

Experiment Design and Result Analysis
The experiment is divided into a three-parts simulation experiment, hardware-inthe-loop simulation experiment and anti-interference experiment to verify the stability and performance of the system under different conditions. In the proposed method, the default value of the initial step size 0 µ and the proportional variable β in Equation

Experiment Design and Result Analysis
The experiment is divided into a three-parts simulation experiment, hardware-inthe-loop simulation experiment and anti-interference experiment to verify the stability and performance of the system under different conditions. In the proposed method, the default value of the initial step size 0 µ and the proportional variable β in Equation

Experiment Design and Result Analysis
The experiment is divided into a three-parts simulation experiment, hardware-in-theloop simulation experiment and anti-interference experiment to verify the stability and performance of the system under different conditions. In the proposed method, the default value of the initial step size µ 0 and the proportional variable β in Equation (38) is µ 0 = 0.01, β = 10. The zero bias of the self-designed MARG sensor system and random error standard deviation are shown in Table 1, which can be used to calculate the system noise and measurement noise of the proposed Kalman filter. The parameters of the proposed Kalman filter are shown in Equation (52).

Simulation Experiment
In order to verify the reliability and stability of the designed attitude calculation method, a typical sine motion model is used as the input for the simulation analysis to verify the accuracy and precision of the proposed attitude calculation method. Taking a sine motion model with a frequency of 0.2 Hz and an amplitude of 10 • as the input, the result of the attitude angle solution is shown in Figures 6-8. The left side is the comparison diagram between the Euler angle solution result and the input motion model, and the right side shows the error between the solution result and the true value. It can be seen from the figure that the proposed attitude calculation method has a maximum error of about 0.4 • for the roll angle and pitch angle, and a maximum error of about 0.5 • for the yaw angle.

Simulation Experiment
In order to verify the reliability and stability of the designed attitude calculation method, a typical sine motion model is used as the input for the simulation analysis to verify the accuracy and precision of the proposed attitude calculation method. Taking a sine motion model with a frequency of 0.2 Hz and an amplitude of 10° as the input, the result of the attitude angle solution is shown in Figures 6-8. The left side is the comparison diagram between the Euler angle solution result and the input motion model, and the right side shows the error between the solution result and the true value. It can be seen from the figure that the proposed attitude calculation method has a maximum error of about 0.4° for the roll angle and pitch angle, and a maximum error of about 0.5° for the yaw angle.  method, a typical sine motion model is used as the input for the simulation analysis to verify the accuracy and precision of the proposed attitude calculation method. Taking a sine motion model with a frequency of 0.2 Hz and an amplitude of 10° as the input, the result of the attitude angle solution is shown in Figures 6-8. The left side is the comparison diagram between the Euler angle solution result and the input motion model, and the right side shows the error between the solution result and the true value. It can be seen from the figure that the proposed attitude calculation method has a maximum error of about 0.4° for the roll angle and pitch angle, and a maximum error of about 0.5° for the yaw angle.
(a) (b)  method, a typical sine motion model is used as the input for the simulation analysis to verify the accuracy and precision of the proposed attitude calculation method. Taking a sine motion model with a frequency of 0.2 Hz and an amplitude of 10° as the input, the result of the attitude angle solution is shown in Figures 6-8. The left side is the comparison diagram between the Euler angle solution result and the input motion model, and the right side shows the error between the solution result and the true value. It can be seen from the figure that the proposed attitude calculation method has a maximum error of about 0.4° for the roll angle and pitch angle, and a maximum error of about 0.5° for the yaw angle.

Static Hardware-in-the-Loop Simulation Experiment
The vibration isolation table can isolate the vibration transmission between the outside world and the measuring device, ensuring that the device is in a static state. In the static experiment, the self-designed measuring device is placed on the precision vibration isolation table steadily, and nine-axis data are collected for 45 min in a static state. The attitude is calculated using the gyroscope angular rate integration method and the proposed method, respectively. Due to the strong ferromagnetic interference in the vibration isolation platform environment, the static z-axis yaw angle is not calculated. The experiment results are shown in Figure 9.

Static Hardware-in-the-Loop Simulation Experiment
The vibration isolation table can isolate the vibration transmission between the outside world and the measuring device, ensuring that the device is in a static state. In the static experiment, the self-designed measuring device is placed on the precision vibration isolation table steadily, and nine-axis data are collected for 45 min in a static state. The attitude is calculated using the gyroscope angular rate integration method and the proposed method, respectively. Due to the strong ferromagnetic interference in the vibration isolation platform environment, the static z -axis yaw angle is not calculated. The experiment results are shown in Figure 9. It can be obtained from Figure 9 that due to the existence of the gyroscope zero bias, the error increases and diverges with the angular rate integral. The proposed attitude calculation method can effectively suppress the error caused by the angular rate integral calculation output by the gyroscope.

Dynamic Hardware-in-the-Loop Simulation Experiment
In order to verify the dynamic performance of the proposed Kalman filter, a hardware-in-the-loop simulation platform was built according to Figure 10. A high-precision three-axis turntable was used for verification. By controlling the rotation of the inner frame and outer frame of the turntable, the pitch angle, roll angle, and yaw angle of the measuring device can be simulated. At the same time, the turntable controller can use the RS-422 serial port to realize synchronous position output. It can be obtained from Figure 9 that due to the existence of the gyroscope zero bias, the error increases and diverges with the angular rate integral. The proposed attitude calculation method can effectively suppress the error caused by the angular rate integral calculation output by the gyroscope.

Dynamic Hardware-in-the-Loop Simulation Experiment
In order to verify the dynamic performance of the proposed Kalman filter, a hardwarein-the-loop simulation platform was built according to Figure 10. A high-precision threeaxis turntable was used for verification. By controlling the rotation of the inner frame and outer frame of the turntable, the pitch angle, roll angle, and yaw angle of the measuring device can be simulated. At the same time, the turntable controller can use the RS-422 serial port to realize synchronous position output. Considering that the turntable is mainly composed of ferrous materials, the output of the magnetometer is greatly disturbed during the measurement process. Before the start of the test, the magnetometer ellipsoid fitting method based on the least square method was used to calibrate the output of the magnetometer. During the experiment, the measuring device was initially fixed in the center of the turntable and its XYZ axes were aligned with the NED Navigation Coordinate System. The experiment system was controlled to move quickly around the X and Z axes of the measuring device between 90° and −90°. In order to avoid the singularity problem in the pitch angle calculation process, the measuring device moves quickly between 80° and −80° in the Y-axis direction twice. The move- Considering that the turntable is mainly composed of ferrous materials, the output of the magnetometer is greatly disturbed during the measurement process. Before the start of the test, the magnetometer ellipsoid fitting method based on the least square method was used to calibrate the output of the magnetometer. During the experiment, the measuring device was initially fixed in the center of the turntable and its XYZ axes were aligned with the NED Navigation Coordinate System. The experiment system was controlled to move quickly around the X and Z axes of the measuring device between 90 • and −90 • . In order to avoid the singularity problem in the pitch angle calculation process, the measuring device moves quickly between 80 • and −80 • in the Y-axis direction twice. The movement speed of the turntable is set to 50 • /s, and the acceleration is set to 50 • /s 2 .
Using the solution method proposed in this paper, the extended Kalman (EKF), gradient-descent linear Kalman filter (GDLKF) and gradient-descent (GD) algorithm were used to analyze the experimental data. Where EKF uses the gyroscope output as the state vector and uses the accelerometer and magnetometer to directly calculate the attitude for a posteriori estimation, GDLKF uses gradient pose quaternions as observations for Kalman filtering. The experimental results are shown in Figures 11-13. Considering that the turntable is mainly composed of ferrous materials, the output of the magnetometer is greatly disturbed during the measurement process. Before the start of the test, the magnetometer ellipsoid fitting method based on the least square method was used to calibrate the output of the magnetometer. During the experiment, the measuring device was initially fixed in the center of the turntable and its XYZ axes were aligned with the NED Navigation Coordinate System. The experiment system was controlled to move quickly around the X and Z axes of the measuring device between 90° and −90°. In order to avoid the singularity problem in the pitch angle calculation process, the measuring device moves quickly between 80° and −80° in the Y-axis direction twice. The movement speed of the turntable is set to 50°/s, and the acceleration is set to 50°/s 2 .
Using the solution method proposed in this paper, the extended Kalman (EKF), gradient-descent linear Kalman filter (GDLKF) and gradient-descent (GD) algorithm were used to analyze the experimental data. Where EKF uses the gyroscope output as the state vector and uses the accelerometer and magnetometer to directly calculate the attitude for a posteriori estimation, GDLKF uses gradient pose quaternions as observations for Kalman filtering. The experimental results are shown in Figures 11-13.     It can be seen from Figures 11-13 that the proposed method and EKF can accurately estimate the roll and pitch angles during the test. However, during the rotation of the turntable, the output of the accelerometer is affected by both the acceleration of gravity and the acceleration of the external motion, which result in a relatively large error when the motion state changes suddenly. Relatively, the proposed method in this paper has better dynamic performance than the GDLKF method and GD method.
Since the turntable is anisotropic iron equipment, it is subject to non-uniform magnetic interference during the rotation in the yaw angle calculation process. Although calibrated by the magnetometer, the four methods all produce errors. The proposed method update is significantly better than the other method. In order to more accurately reflect the discrepancy between different calculation methods, the root mean square error (RMSE) of the results of the simulation test and the hardware-in-the-loop simulation experiment is shown in Table 2. It can be seen from the table that the proposed method in It can be seen from Figures 11-13 that the proposed method and EKF can accurately estimate the roll and pitch angles during the test. However, during the rotation of the turntable, the output of the accelerometer is affected by both the acceleration of gravity and the acceleration of the external motion, which result in a relatively large error when the motion state changes suddenly. Relatively, the proposed method in this paper has better dynamic performance than the GDLKF method and GD method.
Since the turntable is anisotropic iron equipment, it is subject to non-uniform magnetic interference during the rotation in the yaw angle calculation process. Although calibrated by the magnetometer, the four methods all produce errors. The proposed method update is significantly better than the other method. In order to more accurately reflect the discrepancy between different calculation methods, the root mean square error (RMSE) of the results of the simulation test and the hardware-in-the-loop simulation experiment is shown in Table 2. It can be seen from the table that the proposed method in this paper has better performance in a dynamic environment and can effectively track and estimate the attitude. In practical use, the environment is complex and there are many interferences. In order to verify the anti-interference performance and operational effect of the proposed attitude calculation method in a complex environment, the measurement device is fixed inside the test vehicle for practical application tests. The measurement device collects the MARG sensor system data and stores them in the data logger. After the experiment, the information in the data logger is read through the reserved interface, and then the calculation is performed in the PC. In this experiment, the third-generation Ellipse-N product of the French SBG company was used for measurement. The product has a built-in dual-frequency four-constellation GNSS module. Its integrated navigation output attitude angle is used as a reference value. The test environment and the vehicle trajectory recorded by GNSS are shown in Figure 14.

Anti-Interference Experiment
In practical use, the environment is complex and there are many interferences. In order to verify the anti-interference performance and operational effect of the proposed attitude calculation method in a complex environment, the measurement device is fixed inside the test vehicle for practical application tests. The measurement device collects the MARG sensor system data and stores them in the data logger. After the experiment, the information in the data logger is read through the reserved interface, and then the calculation is performed in the PC. In this experiment, the third-generation Ellipse-N product of the French SBG company was used for measurement. The product has a built-in dualfrequency four-constellation GNSS module. Its integrated navigation output attitude angle is used as a reference value. The test environment and the vehicle trajectory recorded by GNSS are shown in Figure 14. It can be seen from the satellite trajectory in Figure 14b that the experiment vehicle returns to the starting point after a week of driving, and the attitude of the measuring device is basically the same at the beginning and the end of the experiment, which can be used as an evaluation index for the attitude calculation effect. According to the calculation of the collected data, the result is shown in Figures 15-17. The test results show that the root mean square error of the roll angle and the pitch angle is less than 0.7°, and the root mean square error of the yaw angle is about 1.3°. It can be seen from the satellite trajectory in Figure 14b that the experiment vehicle returns to the starting point after a week of driving, and the attitude of the measuring device is basically the same at the beginning and the end of the experiment, which can be used as an evaluation index for the attitude calculation effect. According to the calculation of the collected data, the result is shown in Figures 15-17. The test results show that the root mean square error of the roll angle and the pitch angle is less than 0.7 • , and the root mean square error of the yaw angle is about 1.3 • .

Conclusions
In this paper, through the fusion of MARG sensor data, a new attitude calculation method combining a gradient descent algorithm and extended Kalman filter is proposed. The accelerometer and magnetometer data are processed through the two-stage gradient descent algorithm to correct the attitude angle, which effectively corrects gyroscope bias errors of the state vector. Meanwhile, compared with traditional external quaternion estimation methods, the proposed method can better eliminate the influence of magnetometer errors on the roll and pitch angles of the carrier. The proposed Kalman filter can provide relatively faster and more accurate attitude measurement results under different working conditions than using gradient descent and the linear Kalman filter alone.