**A Portable Support Attitude Sensing System for Accurate Attitude Estimation of Hydraulic Support Based on Unscented Kalman Filter**

#### **Xuliang Lu, Zhongbin Wang \*, Chao Tan, Haifeng Yan, Lei Si and Dong Wei**

School of Mechatronic Engineering, China University of Mining and Technology, Daxue Road, Xuzhou 221116, China; LuXuliang@cumt.edu.cn (X.L.); tccadcumt@126.com (C.T.); YanHaifeng@cumt.edu.cn (H.Y.); lei.si@cumt.edu.cn (L.S.); weidongcmee@cumt.edu.cn (D.W.)

**\*** Correspondence: wangzbpaper@126.com; Tel./Fax: +86-516-8359-0758

Received: 11 August 2020; Accepted: 15 September 2020; Published: 23 September 2020

**Abstract:** To measure the support attitude of hydraulic support, a support attitude sensing system composed of an inertial measurement unit with microelectromechanical system (MEMS) was designed in this study. Yaw angle estimation with magnetometers is disturbed by the perturbed magnetic field generated by coal rock structure and high-power equipment of shearer in automatic coal mining working face. Roll and pitch angles are estimated using the MEMS gyroscope and accelerometer, and the accuracy is not reliable with time. In order to eliminate the measurement error of the sensors and obtain the high-accuracy attitude estimation of the system, an unscented Kalman filter based on quaternion according to the characteristics of complementation of the magnetometer, accelerometer and gyroscope is applied to optimize the solution of sensor data. Then the gradient descent algorithm is used to optimize the key parameter of unscented Kalman filter, namely process noise covariance, to improve the accuracy of attitude calculation. Finally, an experiment and industrial application show that the average measurement error of yaw angle is less than 2◦ and that of pitch angle and roll angle is less than 1◦, which proves the efficiency and feasibility of the proposed system and method.

**Keywords:** support attitude; inertial measurement unit; coal mining; unscented Kalman filter; quaternion; gradient descent

#### **1. Introduction**

Hydraulic support is an important safety support equipment in the automatic coal mining working face, which is of great significance to the safety of coal miners and the normal operation of coal mining equipment. In recent years, with the continuous development of intelligent mining technology [1,2], the roboticized hydraulic support technology has been paid more and more attention by scholars and coal mining managers. In particular, the precise estimation of the support attitude of hydraulic support directly determines the intelligent ability of hydraulic support. Therefore, it is urgent to study a new attitude sensing method to obtain accurate support attitude of hydraulic support [3], especially in the special application scenarios with a large demand for the number of sensors and complex environment.

The existing studies on the support attitude sensing measurement of hydraulic support are mainly focused on the mechanism and kinematics of hydraulic support and trying to achieve the accurate estimation of support attitude by kinematic analysis of its mechanism [4]. In [5,6], Polish scholars analyzed the structural composition and empirical formula of support strength of hydraulic support, and developed a new type of support mechanism for vertical displacement of hydraulic support, but which can only approximately estimate the yaw angle of hydraulic support by experience. In [7,8], through the analysis of the four-linkage mechanism and parameters of hydraulic support, the kinematics equation is established, and then the support attitude can only be roughly estimated

due to the neglected clearance between the hinge joint and the pin shaft. In [9], a laser radar detection method is used for measuring the relative attitude angle of hydraulic support based on the position of inspection robot, but the estimation result is disturbed greatly due to the influence of robot motion vibration. In [10,11], a mechanical-electrical-hydraulic coordination technique based on the dynamic response of the load is exploited to adjust the attitude of hydraulic support, but it is easily affected by the pressure of the pump station. In [12], a virtual adjustment method of support attitude under the propulsive state of hydraulic support groups is proposed and the estimation accuracy is highly dependent on the modeling precision of the virtual model, and different hydraulic supports need different virtual models, which is difficult to establish and has poor universality. In [13], an approach for monitoring the posture of hydraulic support with fiber Bragg grating tilt sensor based on gravity is presented, but it is not applied to measure yaw angle of the canopy. However, in automatic coal mining working face with the complicated geological structures, the above methods are characterized by complex solutions, difficult modelling, large estimation disturbance and lack of universality with estimating attitude or relative attitude indirectly based on kinematic parameters of mechanism, or can only estimate a single attitude angle. In this paper, we try to explore a new approach for estimating directly the support attitude of hydraulic support with comprehensive consideration of estimation accuracy, measurement method, universality and sensor size.

Inertial measurement unit (IMU) can be independent of the existing mechanism of hydraulic support to estimate support attitude directly, which is not affected by the mechanism of hydraulic support itself. In recent years, the development of microelectromechanical system technology has accelerated leaps and bounds the development process of IMU, which makes the IMUs become low in cost, low in power consumption and small in size. Additionally, it is widely applied in robotics [14,15], navigation [16,17] and virtual reality [18]. Scientific research is not limited to sensor applications [19], but also includes the performance of IMUs [20]. The performance improvement and calibration analysis of IMU provide precise attitude estimation technologies for its application [21], such as the method of improving estimation performance based of sine rotation vector; the acceleration and magnetic field measurements are transformed into the differences between the Euler angles of the measured attitude and the predicted attitude to correct the predicted attitude [22]. However, the measurement accuracy and performance of these low-cost MEMS inertial measurement units are easily limited by the complex working environment [23,24], especially the automatic coal mining working face with complex magnetic interference produced by coal structure and high-power electromechanical equipment. Therefore, it is necessary to use a filtering algorithm to improve the estimation accuracy.

The Kalman filter (KF) and particle filter (PF) are the most studied and widely used filtering algorithms. The Kalman filter for linear filtering and the extended Kalman filter (EKF) with Jacobian matrix and the unscented Kalman filter (UKF) for nonlinear filtering are frequently used to process inertial information, including object tracking [25] and rotation estimation [26]. In order to reduce yaw drift of attitude rotation estimation in indoor scenarios, an improved EKF combining INS mechanization algorithm and zero velocity update methodology is proposed to improve the accuracy of yaw estimation [27], but the Jacobian matrix is used to update the state transition matrix and the observation matrix in the linearization steps of system equation in attitude determination, which leads to poor stability and large estimation bias. Owing to the system, equations are updated based on the sigma-point approximation generated by unscented transformation in UKF [28,29], which has been demonstrated to have more attitude solution accuracy and strong robustness than EKF in estimating attitude with highly nonlinear kinematics models. In [30,31], a quaternion-based MEMS sensors fusion algorithm and heading estimation approach based on rotation matrix and principal component analysis are proposed to improve the heading estimation accuracy of indoor pedestrians, which shows that quaternion is easier to implement than other methods in sensor fusion technology. However, the particle filter is a sequential importance sampling filtering method based on Bayesian estimation [32] and its idea is based on the Monte Carlo method to represent probability with particle set, which can be used in any state space models [33,34]. In [35], a method of filtering and predicting

time-varying signals under model uncertainty is proposed to realize dynamic estimation of the data. In [36], a particle metropolis Hastings algorithm driven by multiple parallel particle filters is proposed to make inference on dynamic and static variables. In any case, although the particle filter has good estimation performance in nonlinear filtering, the mathematical process of particle filters is more complex, and the algorithm implementation and calculation are difficult. Moreover, the resampling, particle degradation and calculation amount in the calculation process lead to the complexity of the implementation process and poor real-time performance. Therefore, the operation time is much longer than that of UKF.

Through the above research on attitude estimation of nonlinear systems, and considering the estimation accuracy, calculation speed and real-time performance of algorithm, UKF is proved to be a powerful technique for estimating attitude angle and a superior alternative to the EKF and PF in various nonlinear system filtering [37,38]. Unfortunately, the key parameter of UKF selected by experience, process noise covariance Q, has a great influence on the estimation accuracy in the complex working environment [39], so the research on tuning of parameter Q is very meaningful [40,41]. Therefore, the gradient descent algorithm [42] with the advantages of fast global convergence, faster operation speed and simple realization could be used to tune the process noise covariance Q to improve the estimation performance of UKF.

In this paper, a support attitude sensing system with a character of intrinsic safety is designed to measure the support attitude of hydraulic support in the special application scenarios with a large demand for the number of sensors and complex environment. The designed support attitude sensing system is a nine-axis low-cost MEMS measurement unit composed of a gyroscope, magnetometer and accelerometer. The proposed gradient descent algorithm-optimized quaternion-based unscented Kalman filter makes full use of the characteristics of complementation of the magnetometers without long-term drift, accelerometers and gyroscopes which are not affected by magnetic disturbance to improve the support attitude estimation accuracy of hydraulic support.

The remaining parts of this paper are organized as follows. In Section 2, the support attitude sensing system designed in our laboratory is described. Section 3 introduces the proposed optimized quaternion-based unscented Kalman filter based on the complementary characteristics of MEMS sensors in detail. In Section 4, an experiment is conducted and analyzed to validate the estimation performance of the proposed system and approach. Then, industrial applications and tests are carried out in Section 5. Section 6 summarizes our conclusions and future work.

#### **2. Support Attitude Sensing System**

#### *2.1. Establishing the Coordinate System and Attitude Angle*

The two coordinate systems and three attitude angles need to be established in the support attitude measurement of hydraulic support, as shown in Figure 1, which are the geographical coordinate system, carrier coordinate system, pitch angle, yaw angle and roll angle. The geographical coordinate system (frame n) is the North-East-Down (NED) orthogonal coordinate defined by the relative geoid. The Down axis (D axis) is perpendicular to the ellipsoid of the reference system and points to the interior of the earth. The North axis (N axis) points to true north (the velocity vector projection along the earth's rotation angle on the plane perpendicular to the D axis). Finally, the East axis (E axis) points horizontally east and completes the right-hand orthogonal coordinate system. The carrier coordinate system (frame b) is fixed to the carrier, which is the reference coordinate system. The *Zb* axis and *Xb* axis point to the bottom and forward along the longitudinal axis of the carrier, respectively. The *Yb* axis points to the right wing (side) and completes the right-hand orthogonal coordinate system with the axes *Xb* and *Zb*.

**Figure 1.** Definition of coordinate systems and carrier attitude angles: (**a**) local geographical coordinate system (*N'E'D'*) defined relative to the geoid; (**b**) carrier coordinate system (*XbYbZb*) and carrier attitude angles including roll, pitch and yaw angles.

The support attitude of hydraulic support includes yaw angle, roll angle and pitch angle. The yaw angle of rotation around *Zb* axis of the carrier is the angle between the N axis and the longitudinal axis *Xb* of the carrier, measured on the horizontal plane and clockwise is positive. The pitch angle of rotation around the *Yb* axis of the carrier is the angle between the horizontal plane and the longitudinal axis *Xb* of carrier, measured on the vertical plane and the direction of the carrier's head up is positive. The roll angle of rotation around the *Xb* axis of the carrier is the angle between the horizontal plane and the cross axis *Yb* of the carrier, measured on the cross section, and the lifting direction of the carrier on the left side is positive.

The carrier coordinate system and NED coordinate system can coordinate transformation with each other by rotating yaw angle, pitch angle and roll angle around the X axis, Y axis and Z axis in turn, and the transformation relationship from ONED to *OXbYbZb* is shown in Figure 2. The coordinate transformation matrix is also called direction cosine matrix, which is written as:

$$\begin{aligned} \mathbf{C}\_{n}^{b} &= [\phi]\_{\mathcal{N}\_{2}} [\theta]\_{E\_{1}} [\psi]\_{\mathcal{D}} \\ &= \begin{bmatrix} \cos\theta\cos\psi & \cos\theta\sin\psi & -\sin\theta \\ \cos\psi\sin\phi\sin\theta - \cos\phi\sin\psi & \sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi & \sin\phi\cos\theta \\ \cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi & \cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi & \cos\phi\cos\theta \end{bmatrix} \end{aligned} \tag{1}$$

**Figure 2.** Coordinate transformation and directional cosine matrix definition of the canopy between frame b and frame n based on hydraulic support.

On the contrary, the transformation from *OXbYbZb* to ONED is written as:

$$\mathbf{C}\_{b}^{\eta} = \begin{bmatrix} \cos\theta\cos\psi & \sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi & \cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi\\ \cos\theta\sin\psi & \sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi & \cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi\\ -\sin\theta & \sin\phi\cos\theta & \cos\phi\cos\theta \end{bmatrix} \tag{2}$$

where φ, θ and ψ are roll angle, pitch angle and yaw angle of support attitude of hydraulic support canopy, respectively.

#### *2.2. Intrinsically Safe Micro Inertial Sensor*

The support attitude sensing system based on the intrinsically safe micro inertial sensor is designed as a small strapdown inertial navigation system in our lab to meet the coal mine special environment navigation requirements of hydraulic support, mine inspection robot, support robot and driving robot. The intrinsically safe micro inertial sensor combining the IMU and magnetometer is composed of a three-axis gyroscope, three-axis magnetometer, three-axis accelerometer, an intrinsically safe microprocessor and a flameproof shell. The three single-axis gyroscopes and accelerometer use the ICM-42605 with the measuring range of ±125◦/s and ±16 g. The three single-axis magnetometers use the LSM3030C with a measurement range of ±16 gauss magnetic full scale. STM32MP157 is selected as the microprocessor of the system and it is based on the high-performance dual-core Arm® Cortex®-A7 32-bit RISC core with an operating frequency of up to 650 MHz. The support attitude sensing system uses an I2C/SPI bus for digital output.

#### *2.3. Support Attitude Estimation*

According to the mining process of automatic coal mining working face, there are two working states in the support process for the canopy of hydraulic support, namely, static working state and moving working state. An accelerometer with superior static characteristics has high measurement accuracy in static state. A gyroscope with stable static performance has high measurement accuracy in moving state. A magnetometer with the measurement stability in a long time test is easily disturbed by the external perturbed magnetic sources produced by coal rock structure and high-power shearer equipment. Nevertheless, the gyroscope is not disturbed by magnetic sources from the automatic coal mining working face. Therefore, this paper proposes a new attitude estimation method for measuring the pitch angle, yaw angle and roll angle of hydraulic support based on the characteristics of complementation of accelerometers, gyroscopes and magnetometers, which eliminates the errors derived from the MEMS sensors by the quaternion-based UKF.

In addition, the tuning parameter of UKF, the process noise covariance Q, is usually a general value set by experience, which sometimes makes UKF produce large estimation errors, especially in a complex environment. Therefore, a gradient descent algorithm is used to optimize Q to improve the estimation accuracy and performance in this paper. The flow chart of accurate estimation of support attitude of hydraulic support using the quaternion-based unscented Kalman filter through tuning process noise covariance is shown in Figure 3. The other key tuning parameter, the observation noise covariance R, is related to the deviation of the magnetometer, accelerometer and gyroscope, which can be obtained by calculating the expectation of the square of deviation of the sensors.

**Figure 3.** Block diagram of support attitude sensing system based on unscented Kalman filter for attitude measurement according to the complementary characteristics of the MEMS inertial sensor.

#### 2.3.1. Support Attitude Estimation Based on Magnetometer

The magnetometer output is the sensitive geomagnetic field intensity. On the earth's surface, the geomagnetic field always points to the north along the magnetic induction, with components in the north direction and vertical direction of that, and no component in the east direction. Therefore, magnetic field coordinate is (N, 0, D) in the geomagnetic coordinate system. The magnetometer output is *Mmag <sup>n</sup>* = *MN* 0 *MD <sup>T</sup>* when the geomagnetic coordinate system coincides with the frame b. The magnetometer output is *Mn* = *Mx My Mz <sup>T</sup>* in frame n and *mb* <sup>=</sup> *mx my mz <sup>T</sup>* in frame b. According to *mb* = *C<sup>b</sup> n* <sup>ψ</sup>=ψ*<sup>m</sup>* ·*Mmag <sup>n</sup>* , the relationship of geomagnetic field intensity between the two coordinate systems can be obtained as:

$$\begin{bmatrix} m\_x \\ m\_y \\ m\_z \end{bmatrix} = \begin{bmatrix} M\_N \cos \theta \cos \psi\_m - M\_D \sin \theta \\ M\_D \sin \phi \cos \theta + M\_N (\sin \phi \sin \theta \cos \psi\_m - \cos \phi \sin \psi\_m) \\ M\_D \cos \phi \cos \theta + M\_N (\cos \phi \sin \theta \cos \psi\_m + \sin \phi \sin \psi\_m) \end{bmatrix} \tag{3}$$

Thus, the components of geomagnetic vector in a horizontal direction are:

$$\begin{cases} \begin{aligned} m\_z \cos \phi \sin \theta + m\_y \sin \phi \sin \theta + m\_x \cos \theta &= M\_N \cos \psi\_m \\ m\_z \sin \phi - m\_y \cos \phi &= M\_N \sin \psi\_m \end{aligned} \tag{4}$$

Then, the yaw angle can be expressed as:

$$\psi\_{\rm{m}} = \arctan\left(\frac{M\_{\rm{y}}}{M\_{\rm{x}}}\right) = \arctan\left(-\frac{m\_{\rm{y}}\cos\phi - m\_{\rm{z}}\sin\phi}{m\_{\rm{x}}\cos\theta + m\_{\rm{y}}\sin\theta\sin\phi + m\_{\rm{z}}\sin\theta\cos\phi}\right) \tag{5}$$

In fact, the geomagnetic North Pole and the geographical North Pole do not coincide, and the difference between them is one magnetic declination angle α (calculated by IGRF or WMM geomagnetic model according to the longitude and latitude height of the observation point). Therefore, the real local yaw angle of the carrier obtained based on magnetometer is:

$$
\psi\_{m \neq n} = \psi\_m + \alpha \tag{6}
$$

2.3.2. Support Attitude Estimation Based on Accelerometer

The accelerometer output is the sensitive carrier acceleration, including the acceleration of the support attitude sensing system and the gravitational acceleration. For support attitude angle computation with an accelerometer, the relationship between the gravitational acceleration 0 0 *<sup>g</sup> <sup>T</sup>* in frame n and the gravitational acceleration *ab* = *ax ay az <sup>T</sup>* measured by the accelerometer in frame b could be expressed as:

$$
\begin{bmatrix} a\_x \\ a\_y \\ a\_z \end{bmatrix} = \mathbf{C}\_n^b \begin{bmatrix} 0 \\ 0 \\ g \end{bmatrix} \begin{bmatrix} -g\sin\theta \\ g\cos\theta\sin\phi \\ g\cos\theta\cos\phi \end{bmatrix} \tag{7}
$$

According to the acceleration output of Equation (7), the roll and pitch angles obtained based on accelerometer can be calculated as follows:

$$\begin{cases} \begin{aligned} \theta &= -\arctan\left(\frac{a\_x}{\sqrt{a\_y^2 + a\_z^2}}\right) \\ \phi &= \arctan\left(\frac{a\_y}{a\_z}\right) \end{aligned} \tag{8} \end{aligned} \tag{8}$$

2.3.3. Support Attitude Estimation Based on Gyroscope

The gyroscope output is the sensitive carrier angular acceleration. In this paper, quaternion is applied to calculate the support attitude of hydraulic support. Quaternion has the advantage of faster calculation speed, smooth interpolation, effectively avoiding the universal lock problem and small storage space. Quaternion is a vector composed of four elements, and matrix expression and normalized quaternion of that are defined as follows:

$$q = \begin{bmatrix} q\_0 & q\_1 & q\_2 & q\_3 \end{bmatrix}^T \begin{pmatrix} q\_0^2 + q\_1^2 + q\_2^2 + q\_3^2 = 1 \end{pmatrix} \tag{9}$$

Due to the movement of the carrier, *q* is constantly updated. The updated dynamic model based on quaternion using the angular velocity of the gyroscope is expressed as:

$$
\dot{q} = \frac{1}{2} q \otimes \omega^b \tag{10}
$$

where ω*<sup>b</sup>* = ω*<sup>x</sup>* ω*<sup>y</sup>* ω*<sup>z</sup>* is the angular velocity and measured based on the gyroscope. Then, the matrix form of Equation (10) is as follows:

.

$$
\begin{vmatrix} q\_0 \\ \dot{q}\_1 \\ \dot{q}\_2 \\ \dot{q}\_3 \end{vmatrix} = \frac{1}{2} \begin{vmatrix} 0 & -\omega\_x & -\omega\_y & -\omega\_z \\ \omega\_x & 0 & \omega\_z & -\omega\_y \\ \omega\_y & -\omega\_z & 0 & \omega\_x \\ \omega\_z & \omega\_y & -\omega\_x & 0 \end{vmatrix} \begin{vmatrix} q\_0 \\ q\_1 \\ q\_2 \\ q\_3 \end{vmatrix} \tag{11}
$$

The Picard successive approximation method of quaternion differential equation is used to solve Equation (11), and the discrete form is as follows:

$$q(t\_k) = \left[\cos(\frac{\Delta\alpha}{2})l + \sin(\frac{\Delta\alpha}{2})\frac{\Delta\Omega}{\Delta\alpha}\right]q(t\_{k-1})\tag{12}$$

where Δα and ΔΩ involve the integration of angular velocity ω*<sup>b</sup>* in the *k*th sampling period, *k* = 1, 2, ... n, Δα = α2 *<sup>x</sup>* + α<sup>2</sup> *<sup>y</sup>* + α<sup>2</sup> *<sup>z</sup>* is the angular increment of the gyroscope in the [*tk-1*,*tk*] sampling interval,

$$\alpha\_i = \int\_{t\_{k-1}}^{t\_k} \omega\_i dt, i = \mathbf{x}, y, z, \text{and } \Delta \Omega = \int\_{t\_{k-1}}^{t\_k} \begin{bmatrix} 0 & -\alpha\_x & -\alpha\_y & -\alpha\_z \\ \alpha\_x & 0 & \alpha\_z & -\alpha\_y \\ \alpha\_y & -\alpha\_z & 0 & \alpha\_x \\ \alpha\_z & \alpha\_y & -\alpha\_x & 0 \end{bmatrix} dt, \text{ $t$  is the sampling interval.}$$

*Sensors* **2020**, *20*, 5459

The coordinate transformation matrix expressed by direction cosine, Euler angle and quaternion from carrier coordinate system to the NED coordinate system is written as:

$$\begin{aligned} \mathbf{C}\_{b}^{n} &= \begin{bmatrix} \mathbf{C}\_{11} & \mathbf{C}\_{12} & \mathbf{C}\_{13} \\ \mathbf{C}\_{21} & \mathbf{C}\_{22} & \mathbf{C}\_{23} \\ \mathbf{C}\_{31} & \mathbf{C}\_{32} & \mathbf{C}\_{33} \end{bmatrix} \\ &= \begin{bmatrix} \cos\theta\cos\psi & \sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi & \cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi \\ \cos\theta\sin\psi & \sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi & \cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi \\ -\sin\theta & \sin\phi\cos\theta & \cos\phi\cos\theta \end{bmatrix} \\ &= \begin{bmatrix} q\_1^2 - q\_2^2 - q\_3^2 + q\_0^2 & 2(q\_1q\_2 - q\_0q\_3) & 2(q\_0q\_2 + q\_1q\_3) \\ 2(q\_0q\_3 + q\_1q\_2) & -q\_1^2 + q\_2^2 - q\_3^2 + q\_0^2 & 2(q\_2q\_3 - q\_0q\_1) \\ 2(q\_1q\_3 - q\_0q\_2) & 2(q\_0q\_1 + q\_2q\_3) & -q\_1^2 - q\_2^2 + q\_3^2 + q\_0^2 \end{bmatrix} \end{aligned} \tag{13}$$

Therefore, the support attitude angle can be expressed in the form of Euler angle represented by direction cosine, and is written as:

$$\begin{cases} \phi = \arctan\frac{C\_{32}}{C\_{33}}\\ \theta = \arctan\frac{-C\_{31}}{\sqrt{C\_{32}^2 + C\_{33}^2}}\\ \psi = \arctan\frac{C\_{21}}{C\_{11}} \end{cases} \tag{14}$$

The carrier yaw angle in form of Euler angle represented by quaternion is written as:

$$\psi\_{\mathcal{G}ym} = \arctan\left(\frac{2(q\_1q\_2 + q\_0q\_3)}{q\_0^2 + q\_1^2 - q\_2^2 - q\_3^2}\right) \tag{15}$$

#### **3. The Optimized Quaternion-Based Unscented Kalman Filter**

In this section, we propose an unscented Kalman filter based on quaternion for estimating the support attitude of hydraulic support. In order to improve the performance of the support attitude estimation algorithm, the gradient descent algorithm is used for tuning the key parameter of the unscented Kalman filter, i.e., process noise covariance Q. The whole procedure of accurate estimation of support attitude using the gradient descent algorithm-optimized quaternion-based unscented Kalman filter (GD-UKF) is also presented.

#### *3.1. Gradient Descent Algorithm*

Gradient descent algorithm with the advantages of simple implementation is a common method to solve unconstrained optimization problems. The basic idea of the gradient descent is to select an appropriate initial value *x*(0), update the value of *x* iteratively, and minimize the objective function until it converges. The input of gradient descent is the objective function *f*(*x*), the gradient function *g*(*x*) = ∇*f*(*x*) and the calculation accuracy ε; the output is the minimum point *x*<sup>∗</sup> of *f*(*x*). The process of the gradient descent can be summarized in detail as follows:

Step 1.1: Initial key parameters *x*(0) and *k*.

Step 1.2: Calculate the objective function *f x*(*k*) 

Step 1.3: Calculate the gradient *gk* <sup>=</sup> *<sup>g</sup>*(*x*(*k*)). If *gk* < ε, stop iteration and let *<sup>x</sup>*<sup>∗</sup> <sup>=</sup> *<sup>x</sup>*(*k*). Otherwise, let *pk* <sup>=</sup> <sup>−</sup>*g*(*x*(*k*)) and then find <sup>λ</sup>*<sup>k</sup>* for equation *<sup>f</sup> x*(*k*) + λ*kpk* = min <sup>λ</sup>≥<sup>0</sup> *<sup>f</sup>*(*x*(*k*) <sup>+</sup> <sup>λ</sup>*pk*).

.

Step 1.4: Set the variable *x*(*k*+1) = *x*(*k*) + λ*kpk* and update the function *f x*(*k*+1) . If the expression *<sup>f</sup>*(*x*(*k*+1)) <sup>−</sup> *<sup>f</sup>*(*x*(*k*)) < ε or *x*(*k*+1) <sup>−</sup> *<sup>x</sup>*(*k*) < ε, stop iteration and let *<sup>x</sup>*<sup>∗</sup> <sup>=</sup> *<sup>x</sup>*(*k*). Otherwise, let *<sup>k</sup>* <sup>=</sup> *<sup>k</sup>* <sup>+</sup> <sup>1</sup> and go back to Step 1.3 to find another better minimum.

When the objective function is a convex function, the solution of gradient descent is the global optimal solution.

#### *3.2. Unscented Transformation*

For unscented Kalman filter in a nonlinear system, unscented transformation (UT) is the core technique for propagating mean and covariance based on Cholesky decomposition, and can effectively approximate mean and covariance changes of the random variables when it undergoes a nonlinear transformation, including cross-correlation between state and measurement. The basic principle of UT can be considered in this way. We assume that the mean and covariance of a random variable *X* with n dimensions are *x* and *P*, respectively, and mean *x* and covariance *P* propagate through a nonlinear function *y* = *f*(*x*). In order to calculate the statistics of the variable y, 2n + 1 sigma points *X*(*i*) and corresponding weights *W*(*i*) are formed and calculated as follows:

#### (1) Calculate 2n + 1 sigma points

$$\begin{cases} X^{(0)} = \overline{x}, i = 0\\ X^{(i)} = \overline{x} + \left(\sqrt{(n+\lambda)P}\right)\_{i'} i = 1 \sim n\\ X^{(i)} = \overline{x} - \left(\sqrt{(n+\lambda)P}\right)\_{i'} i = n+1 \sim 2n \end{cases} \tag{16}$$

where √*p T* <sup>√</sup>*<sup>p</sup>* = *p*, √*p i* is the *i-*th column of the matrix of the square root. (2) Calculate the corresponding weight of sigma point

$$\begin{cases} \begin{aligned} \mathcal{W}\_m^{(0)} &= \frac{\lambda}{n+\lambda} \\ \mathcal{W}\_c^{(0)} &= \frac{\lambda}{n+\lambda} + \left(1 - \alpha^2 + \beta\right) \\ \mathcal{W}\_m^{(i)} = \mathcal{W}\_c^{(i)} &= \frac{\lambda}{2(n+\lambda)}, i = 1 \sim 2n \end{aligned} \tag{17}$$

where <sup>λ</sup> = <sup>α</sup>2(*<sup>n</sup>* + <sup>κ</sup>) <sup>−</sup> *<sup>n</sup>* is the scaling factor for reducing total prediction error; <sup>α</sup> is set to a small positive value to control the distribution of sigma points; κ is the parameter to be selected to ensure a positive semidefinite and is usually set to 0; and β is a nonnegative weight coefficient and is set for 2 for Gaussian distribution in this paper. The subscript *m* and *c* represent covariance and mean, respectively.

#### *3.3. Unscented Kalman Filter Design*

Generally, the dynamic system of the unscented Kalman filter can be described in two system equations: firstly, state equation; and secondly, observation equation. Considering the discrete-time nonlinear dynamic model of support attitude estimation of hydraulic support, the model state variable is *x* = *q*<sup>0</sup> *q*<sup>1</sup> *q*<sup>2</sup> *q*<sup>3</sup> *T* , which is a vector composed of four elements of quaternion; and the model observation variable is *z* = φθψ *<sup>T</sup>* , which is a vector composed of attitude angles. The state variable is unknown and is constantly updated according to the gyroscope; the observation variable is known and measured by the accelerometer and magnetometer. The details of the mathematical model are described with the following steps.

#### (1) State equation

The state prediction is based on the previous optimal estimation, and the discrete-time nonlinear dynamic state equation of unscented Kalman filter is shown as:

$$\mathbf{x}\_{k} = A\_{k-1}\mathbf{x}\_{k-1} + w\_{k-1} \tag{18}$$

where *xk* and *xk*−<sup>1</sup> are the prior estimation at time *k* and the posterior estimation at time *k* − 1, respectively; *wk* is the process noise with covariance Q and simplified as independent Gaussian white noise. *Ak* is the state transition matrix and can be obtained as:

$$A = \cos(\frac{\Delta\alpha}{2})l + \sin(\frac{\Delta\alpha}{2})\frac{\Delta\Omega}{\Delta\alpha} \tag{19}$$

#### (2) Observation equation

The correction state is the essential step for refining measurement estimation. The observation equation is expressed as

$$z\_k = H[\mathbf{x}\_k] + v\_k \tag{20}$$

where *vk* is the independent Gaussian measurement noise with noise covariance R. *Hk*[*xk*] is the output function and can be calculated as follows:

$$H[\mathbf{x}] = \begin{bmatrix} \arctan\left(2(q\_2q\_3 + q\_0q\_1) / \left(q\_0^2 - q\_1^2 - q\_2^2 + q\_3^2\right)\right) \\ \arcsin\left(2(q\_0q\_2 - q\_1q\_3)\right) \\ \arctan\left(2(q\_1q\_2 + q\_0q\_3) / \left(q\_1^2 - q\_2^2 - q\_3^2 + q\_0^2\right)\right) \end{bmatrix} \tag{21}$$

#### *3.4. Noise Covariance of Process and Observation*

The observation noise covariance and process noise covariance are the two key parameters of unscented Kalman filter, which sometimes leads to large support attitude estimation errors. However, it is essential to optimize the noise covariance of the process and observation to minimize the support attitude estimation errors of hydraulic support.

#### (1) Process noise covariance

In the process of hydraulic support moving and canopy supporting in the automatic coal mining working face, the uneven floor could cause random acceleration of the hydraulic support movement. Meanwhile the coupling effect of roof rock and hydraulic support can produce random impact on the hydraulic support in the moving process and static state, and the left and right adjacent hydraulic support could also produce impact vibration on hydraulic support. Therefore, the coupling effect between hydraulic support and the surrounding environment could lead to random vibration and acceleration on the support attitude sensing system in the static and moving process of hydraulic support, which leads to difficulty in determining the process covariance Q.

We assume that ω*<sup>x</sup>* = ω*<sup>x</sup>* + *wx*, ω*<sup>y</sup>* = ω*<sup>y</sup>* + *wy*, ω*<sup>z</sup>* = ω*<sup>z</sup>* + *wz* where ω*x*, ω*y*, ω*<sup>z</sup>* are the mean of ω*x*, ω*y*, ω*z*; *wx*, *wy*, *wz* are the process deviations of the support attitude sensing system caused by the disturbance of coupling interaction between the hydraulic support and the surrounding environment. σ<sup>2</sup> <sup>ω</sup>*<sup>x</sup>* = *E wx* · *wT x* , σ<sup>2</sup> <sup>ω</sup>*<sup>y</sup>* = *E wy* · *wT y* and σ<sup>2</sup> <sup>ω</sup>*<sup>z</sup>* = *E wz* · *wT z* are the variance of *wx*, *wy*, *wz*. Thus, Equation (11) can be rewritten as:

⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ . *q*0 . *q*1 . *q*2 . *q*3 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ = <sup>1</sup> 2 ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0 −(ω*<sup>x</sup>* + *wx*) − ω*<sup>y</sup>* + *wy* <sup>−</sup>(ω*<sup>z</sup>* <sup>+</sup> *wz*) (ω*<sup>x</sup>* + *wx*) 0 (ω*<sup>z</sup>* + *wz*) − ω*<sup>y</sup>* + *wy* ω*<sup>y</sup>* + *wy* <sup>−</sup>(ω*<sup>z</sup>* <sup>+</sup> *wz*) <sup>0</sup> (ω*<sup>x</sup>* <sup>+</sup> *wx*) (ω*<sup>z</sup>* + *wz*) ω*<sup>y</sup>* + *wy* <sup>−</sup>(ω*<sup>x</sup>* <sup>+</sup> *wx*) <sup>0</sup> ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *q*0 *q*1 *q*2 *q*3 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ = <sup>1</sup> 2 ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0 −ω*<sup>x</sup>* −ω*<sup>y</sup>* −ω*<sup>z</sup>* ω*<sup>x</sup>* 0 ω*<sup>z</sup>* −ω*<sup>y</sup>* ω*<sup>y</sup>* −ω*<sup>z</sup>* 0 ω*<sup>x</sup>* ω*<sup>z</sup>* ω*<sup>y</sup>* −ω*<sup>x</sup>* 0 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ *q*0 *q*1 *q*2 *q*3 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ + <sup>1</sup> 2 ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −*q*<sup>1</sup> −*q*<sup>2</sup> −*q*<sup>3</sup> *q*<sup>0</sup> −*q*<sup>3</sup> *q*<sup>2</sup> *q*<sup>3</sup> *q*<sup>0</sup> −*q*<sup>1</sup> −*q*<sup>2</sup> *q*<sup>1</sup> *q*<sup>0</sup> ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎣ *wx wy wz* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎦ (22)

The second term on the right of Equation (22) can be considered as the process noise *w* = <sup>1</sup> 2 ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −*q*<sup>1</sup> −*q*<sup>2</sup> −*q*<sup>3</sup> *q*<sup>0</sup> −*q*<sup>3</sup> *q*<sup>2</sup> *q*<sup>3</sup> *q*<sup>0</sup> −*q*<sup>1</sup> −*q*<sup>2</sup> *q*<sup>1</sup> *q*<sup>0</sup> ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎣ *wx wy wz* ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎦ and the process noise covariance Q is:

$$Q = E\{w \cdot w^T\} = \frac{1}{4} \begin{bmatrix} -q\_1 & -q\_2 & -q\_3 \\ q\_0 & -q\_3 & q\_2 \\ q\_3 & q\_0 & -q\_1 \\ -q\_2 & q\_1 & q\_0 \end{bmatrix} \begin{bmatrix} \sigma\_{\alpha\_x}^2 & 0 & 0 \\ 0 & \sigma\_{\alpha\_y}^2 & 0 \\ 0 & 0 & \sigma\_{\alpha\_z}^2 \end{bmatrix} \begin{bmatrix} -q\_1 & -q\_2 & -q\_3 \\ q\_0 & -q\_3 & q\_2 \\ q\_3 & q\_0 & -q\_1 \\ -q\_2 & q\_1 & q\_0 \end{bmatrix}^T \tag{23}$$

However, it can be determined from Equation (23) that the optimal Q of UKF can be obtained by optimizing variance σ<sup>2</sup> <sup>ω</sup>*<sup>x</sup>* , <sup>σ</sup><sup>2</sup> <sup>ω</sup>*<sup>y</sup>* and <sup>σ</sup><sup>2</sup> <sup>ω</sup>*<sup>z</sup>* based on the optimization algorithm.

#### (2) Observation noise covariance

The observation noise covariance R is determined by the measurement process and related to the characteristics of the measuring instrument, which can be obtained through long-term probability statistics of sensor measurement data. In this paper, the observation noise covariance R is determined by the measurement deviation of accelerometer and magnetometer.

We assume that *vax* , *vay* and *vaz* are the observation deviation of accelerometer three-axis output *ax*, *ay* and *az*. Then, the first-order Taylor series expansion of trigonometric function is applied to Equation (8), and the deviation of the roll and pitch angles measured by the accelerometer can be derived as

$$
\begin{bmatrix} v\_{0} \\ v\_{\phi} \end{bmatrix} = \begin{bmatrix} -\frac{\sqrt{a\_{y}^{2} + a\_{z}^{2}}}{a\_{x}^{2} + a\_{y}^{2} + a\_{z}^{2}} & \frac{a\_{t}a\_{y}}{\left(a\_{x}^{2} + a\_{y}^{2} + a\_{z}^{2}\right)\sqrt{a\_{y}^{2} + a\_{z}^{2}}} & \frac{a\_{t}a\_{t}}{\left(a\_{x}^{2} + a\_{y}^{2} + a\_{z}^{2}\right)\sqrt{a\_{y}^{2} + a\_{z}^{2}}} \\ 0 & \frac{a\_{z}}{a\_{y}^{2} + a\_{z}^{2}} & -\frac{a\_{y}}{a\_{y}^{2} + a\_{z}^{2}} \end{bmatrix} \begin{bmatrix} v\_{a\_{x}} \\ v\_{a\_{y}} \\ v\_{a\_{z}} \end{bmatrix} = \Phi\_{0\phi} \begin{bmatrix} v\_{a\_{x}} \\ v\_{a\_{y}} \\ v\_{a\_{z}} \end{bmatrix} \tag{24}
$$

Through Equation (24), the covariance of pitch angle and roll angle based on the accelerometer is as follows:

$$
\sigma\_{\partial \phi}^2 = E \begin{pmatrix} \begin{bmatrix} \upsilon\_{\partial} & \upsilon\_{\phi} \end{bmatrix}^T \cdot \begin{bmatrix} \upsilon\_{\partial} & \upsilon\_{\phi} \end{bmatrix} \end{pmatrix} = \Phi\_{\partial \phi} \begin{bmatrix} \sigma\_{a\_x}^2 & 0 & 0 \\ 0 & \sigma\_{a\_y}^2 & 0 \\ 0 & 0 & \sigma\_{a\_z}^2 \end{bmatrix} \Phi\_{\partial \phi}^T \tag{25}
$$

where σ<sup>2</sup> *ax* = *E vax* · *v<sup>T</sup> ax* , σ<sup>2</sup> *ay* = *E vay* · *v<sup>T</sup> ay* and σ<sup>2</sup> *az* = *E vaz* · *vT az* is the variance of accelerometer three-axis output *ax*, and *az*.

Similarly, through Equations (3) and (5) in the first-order Taylor series expansion of trigonometric function, the deviation of yaw angle measured by magnetometer can be derived as

$$\begin{aligned} \;^l \psi\_{\text{my}\text{r}} &= \begin{bmatrix} -\frac{M\_{\text{y}}}{M\_{\text{x}}^2 + M\_{\text{y}}^2} & \frac{M\_{\text{x}}}{M\_{\text{x}}^2 + M\_{\text{y}}^2} \\\\ -\frac{M\_{\text{y}}}{M\_{\text{x}}^2 + M\_{\text{y}}^2} & \frac{M\_{\text{x}}}{M\_{\text{x}}^2 + M\_{\text{y}}^2} \end{bmatrix} \begin{bmatrix} \cos\theta & \sin\theta\sin\phi & \sin\theta\cos\phi \\\\ 0 & -\cos\phi & \sin\phi \end{bmatrix} \begin{bmatrix} v\_{m\_{\text{x}}} \\ v\_{m\_{\text{y}}} \\ v\_{m\_{\text{z}}} \end{bmatrix} = \Phi\_{\psi} \begin{bmatrix} v\_{m\_{\text{x}}} \\ v\_{m\_{\text{y}}} \\ v\_{m\_{\text{z}}} \end{bmatrix} \end{aligned} \tag{26}$$

where *vmx* , *vmy* and *vmz* are the deviation of magnetometer output *mx*, *my* and *mz*.

Through Equation (26), the covariance of yaw angle based on magnetometer is as follows:

$$\boldsymbol{\sigma}\_{\psi} = \boldsymbol{E} \begin{pmatrix} \boldsymbol{\upsilon}\_{\psi\_{\text{mag}}} \cdot \boldsymbol{\upsilon}\_{\psi\_{\text{mag}}}^T \end{pmatrix} = \boldsymbol{\Phi}\_{\psi} \begin{bmatrix} \boldsymbol{\sigma}\_{\boldsymbol{m}\_x}^2 & 0 & 0 \\ 0 & \boldsymbol{\sigma}\_{\boldsymbol{m}\_y}^2 & 0 \\ 0 & 0 & \boldsymbol{\sigma}\_{\boldsymbol{m}\_z}^2 \end{bmatrix} \boldsymbol{\Phi}\_{\psi}^T \tag{27}$$

where σ<sup>2</sup> *mx* = *E vmx* · *v<sup>T</sup> mx* , σ<sup>2</sup> *my* = *E vmy* · *v<sup>T</sup> my* and σ<sup>2</sup> *mz* = *E vmz* · *vT mz* are the variances of the magnetometer output *mx*, *my* and *mz*, respectively.

Finally, it can be determined from Equations (25) and (27) that the observation noise covariance R is:

$$R = \begin{bmatrix} \sigma\_{0\phi}^2 & 0\\ 0 & \sigma\_{\psi}^2 \end{bmatrix} \tag{28}$$

#### *3.5. Unscented Kalman Filter Based on Gradient Descent*

The UKF with Kalman linear filtering framework is not the traditional method of linearizing nonlinear functions, which uses unscented transformation to propagate mean and covariance. Compared with the extended Kalman filter, the UKF without deriving Jacobian matrix can make the state variable approximate the probability density distribution of nonlinear function and is more robust and accurate. However, the process noise covariance Q affects the filtering performance to a certain extent, and gradient descent can be used to improve the unscented Kalman filter. The flow of gradient descent algorithm-optimized quaternion-based unscented Kalman filter (GD-UKF) can be elaborated as below:

Step 2.1: The key parameters of the initial setup are dimension (n = 4), initial value *x*ˆ0 = *E*(*x*0), covariance initial value *P*<sup>0</sup> = *E* (*x*<sup>0</sup> − *x*ˆ0)(*x*<sup>0</sup> − *x*ˆ0) *T* of random variable. Through Equations (16) and (17), a set of sigma points and the corresponding weights are computed by: *i* = 1, 2, ···, 2*n* + 1

$$X\_{k-1}^{(i)} = \begin{bmatrix} \hat{\mathbf{x}}\_{k-1} & \hat{\mathbf{x}}\_{k-1} + \sqrt{(n+\lambda)P\_{k-1}} & \hat{\mathbf{x}}\_{k-1} - \sqrt{(n+\lambda)P\_{k-1}} \end{bmatrix} \tag{29}$$

Step 2.2: The predicted value and covariance matrix of system state variable calculated through the state Equation (18) based on the sigma points obtained in Step 2.1 are expressed as:

$$\begin{aligned} \mathbf{X}\_k^{(i)} &= A\_{k-1} \mathbf{X}\_{k-1}^{(i)} \\ \mathbf{\hat{x}}\_k^- &= \sum\_{i=0}^{2n} \mathbf{W}\_m^{(i)} \mathbf{X}\_k^{(i)} \\ P\_k^- &= \sum\_{i=0}^{2n} \mathbf{W}\_c^{(i)} \left[ \mathbf{X}\_k^{(i)} - \mathbf{\hat{x}}\_k^- \right] \left[ \mathbf{X}\_k^{(i)} - \mathbf{\hat{x}}\_k^- \right]^T + Q \end{aligned} \tag{30}$$

Step 2.3: Through using UT again for the predicted value in Step 2.2, the new sigma points are calculated as follows: *i* = 1, 2, ···, 2*n* + 1

$$X\_k^{(i)} = \left[ \begin{array}{c} \hat{\mathbf{x}}\_k^- & \hat{\mathbf{x}}\_k^- + \sqrt{(n+\lambda)P\_k^-} \\ \end{array} \begin{array}{c} \hat{\mathbf{x}}\_k^- & \hat{\mathbf{x}}\_k^- - \sqrt{(n+\lambda)P\_k^-} \end{array} \right] \tag{31}$$

Step 2.4: The predicted value and covariance matrix of observation variable calculated through the observation in Equation (20) based on the new sigma points obtained in Step 2.3 are expressed as follows:*i* = 1, 2, ···, 2*n* + 1

$$\begin{aligned} \mathbf{Z}\_k^{(i)} &= H \begin{bmatrix} \mathbf{X}\_k^{(i)} \\ \mathbf{z}\_k^- \end{bmatrix} \\ \mathbf{Z}\_k^- &= \sum\_{i=0}^{2n} W\_{\mathbf{c}}^{(i)} \mathbf{Z}\_k^{(i)} \\ P\_{z\_k z\_k} &= \sum\_{i=0}^{2n} W\_{\mathbf{c}}^{(i)} \left[ \mathbf{Z}\_k^{(i)} - \mathbf{z}\_k^- \right] \left[ \mathbf{Z}\_k^{(i)} - \mathbf{z}\_k^- \right]^T + R \\ P\_{\mathbf{x}, z\_k} &= \sum\_{i=0}^{2n} W\_{\mathbf{c}}^{(i)} \left[ \mathbf{X}\_k^{(i)} - \mathbf{z}\_k^- \right] \left[ \mathbf{Z}\_k^{(i)} - \mathbf{z}\_k^- \right]^T \end{aligned} \tag{32}$$

Step 2.5: The unscented Kalman filter gain matrix *K* is calculated as:

$$K\_k = P\_{x\_k z\_k} P\_{z\_k z\_k}^{-1} \tag{33}$$

Step 2.6: The state update and covariance update of the unscented Kalman filter system are calculated as:

$$\begin{aligned} \mathfrak{X}\_k &= \mathfrak{X}\_k^- + K\_k [z\_k - \mathfrak{z}\_k^-] \\ P\_k &= P\_k^- - K\_k P\_{z\_k z\_k} \mathcal{K}\_k^T \end{aligned} \tag{34}$$

Step 2.7: The process noise covariance Q is updated by calculating and updating process variance σ = σ2 <sup>ω</sup>*<sup>x</sup>* <sup>σ</sup><sup>2</sup> <sup>ω</sup>*<sup>y</sup>* <sup>σ</sup><sup>2</sup> ω*z <sup>T</sup>* of the support attitude sensing system based on gradient descent. And the objective function of gradient descent is

$$f(\sigma) = \mathcal{L}(m\_{b^\*}, a\_{b^\*} \boldsymbol{\omega}^b, \sigma) = \text{RMSE}(\mathbf{x}, \mathbf{\hat{x}}\_k) = \sqrt{\mathcal{E}[(\mathbf{x} - \mathbf{\hat{x}}\_k) \cdot (\mathbf{x} - \mathbf{\hat{x}}\_k)^T]}\tag{35}$$

$$\sigma = \operatorname\*{arg\,min}\_{\sigma} f(\sigma)$$

where *L*(·) is the loss function; *RMSE*(·) is the root mean square error of estimation results; *x* is the true state value of the support attitude sensing system. The flowchart of the proposed GD-UKF is presented in Figure 4.

**Figure 4.** The flowchart of the proposed gradient descent algorithm-optimized quaternion-based unscented Kalman filter.

#### **4. Experiments and Analysis**

In this section, in order to evaluate measurement accuracy and performance of the designed support attitude sensing system based on GD-UKF, the static test and dynamic test experiment are carried out on the two-axis experiment platform. Then, the experimental results are contrasted with the Dutch XSENS high-precision attitude measurement system to prove the out-performing and practicability of the designed system.

#### *4.1. Experiment on Two-Axis Turntable*

The two-axis turntable experimental platform with U-frame structure is a high-precision electric precision machinery and is illustrated in Figure 5. The experiment platform is used for static test with an accuracy of 0.05◦ and dynamic test with an accuracy of 0.1◦. The intrinsically safe micro inertial sensor is mounted on the precision electric rotary table.

**Figure 5.** The experimental platform for testing the support attitude sensing system.

In the actual work of hydraulic support in the automatic coal mining working face, the hydraulic supports are closely arranged to support the surrounding rock. According to the coupling relationship between the hydraulic support and surrounding rock, the hydraulic support can appear yaw, pitch and tilt in the movement process, and the allowable range of the canopy pitch is ±15◦. Additionally, the roll angle of hydraulic support, that is, the inclination of the coal mining working face, is nonnegative and the maximum value is 40◦. The coupling relationship between adjacent hydraulic supports limits the allowable ranges of hydraulic support yaw to ±10◦. Therefore, in the experiment setup, the test ranges of pitch, roll and yaw of the experimental platform are set to ±15◦, 0◦ to 40◦ and ±10◦, respectively, which is consistent with the actual working state of hydraulic support. Meanwhile, the attitude angle estimation performance of the support attitude sensing system based on Kalman filter, unscented Kalman filter and the proposed gradient descent algorithm-optimized quaternion-based unscented Kalman filter are tested and compared, respectively. The dynamic test of pitch angle is set from −15◦ to +15◦ in 0 s to 200 s, and the angle is maintained at 15◦ for static test in 200 s to 500 s. The test experiment of roll angle is designed as a dynamic test from 0◦ to 40◦ in 0 s to 250 s and static test at 40◦ in 250 s to 500 s. The dynamic test from −10◦ to 10◦ in 0 s to 250 s and the static test at 10◦ in 250 s to 500 s of yaw angle are set.

#### *4.2. Result Analysis*

The experimental results of the support attitude sensing system-based KF, UKF and GD-UKF on the two-axis experimental platform are shown in Figures 6–8, in which Figures 6, 7 and 8a show three attitude angles measured by the three methods, and the attitude calculation results of the three methods are consistent with the variational trend set by the experimental platform. In this paper, the root mean square error is used to evaluate the estimation accuracy of the algorithms, and the estimation errors of the three attitude angles are shown in Figures 6b, 7 and 8. The estimation errors (shown in the red curve) of pitch angle and roll angle based on Kalman filter are relatively large, and the maximum errors are greater than 3◦, as shown in Figures 6b and 7b. Compared with the Kalman filter, the unscented Kalman filter improves the estimation accuracy to a certain extent, the blue curve is closer to the green curve than the red curve, but the estimation errors with maximum error greater than 2.4◦ are still relatively large. The reason for this is that the covariance Q, the key parameter of the algorithm, is set to the empirical value. In Figures 6 and 7, the estimation value (the pink curve) based on the proposed GD-UKF, is rather close to the theoretical value (the green curve), the dynamic and static estimation errors are less than 1◦ and 0.5◦, respectively. The reason is that Q is tuned by gradient descent, and the variances σ<sup>2</sup> <sup>ω</sup>*<sup>x</sup>* , <sup>σ</sup><sup>2</sup> <sup>ω</sup>*<sup>y</sup>* and <sup>σ</sup><sup>2</sup> <sup>ω</sup>*<sup>z</sup>* in the optimized Q matrix are infinitely close to the true process noise covariance value caused by the inherent deviation of the experimental platform, as shown in Table 1, which greatly improves the estimation performance of UKF. Therefore, the estimation errors of pitch angle and roll angle based on GD-UKF are distinctly smaller than those of UKF and KF methods, and the accuracies of dynamic estimation and static estimation are better than 1◦ and 0.5◦, respectively.

**Figure 6.** Pitch angle test of support attitude sensing system: (**a**) pitch angle estimation based on three methods; (**b**) pitch error of three methods.

**Figure 7.** Roll angle test of support attitude sensing system: (**a**) roll angle estimation based on three methods; (**b**) roll error of three methods.

**Table 1.** Process noise covariance.


**Figure 8.** Yaw angle test of support attitude sensing system: (**a**) yaw angle estimation based on three methods; (**b**) yaw error of three methods.

As shown in Figure 8a, the fluctuation of yaw angle estimation results of three methods are stable. The maximum estimation errors based on KF and UKF are 1.8◦ and 1.3◦, respectively, as shown in Figure 8b, the average estimation error of yaw angle of GD-UKF is 0.2◦ to 0.4◦ smaller than that of the other two algorithms, which is not obviously demonstrated to have the superior estimation performance of the proposed algorithm. The reason is that there is no magnetic field interference such as high-power equipment in the laboratory, and the estimations of yaw angle by the three methods are almost the same in the environment without magnetic noise interference. For the estimation performance of yaw angle based on GD-UKF, an industrial application needs to be operated in the automatic coal mining working face to validate the superiority and practicability of the proposed algorithm in yaw angle estimation.

The accuracy of our support attitude sensing system, with the static measurement accuracy better than 0.5◦ and dynamic measurement accuracy better than 1◦, is lower than that of MTi-630 AHRS based on industrial-grade MEMS sensor from XSENS [43], but the support attitude sensing system is much cheaper than MTi-630 AHRS and its size is also smaller, which can meet the support accuracy requirements of hydraulic support. Moreover, the support attitude sensing system with 200 mW power consumption has lower power consumption than MTi-630 AHRS with 345 mW power consumption and especially meets the technical requirements of an intrinsically safe system [44], which can be directly applied to the automatic coal mining working face in explosive atmospheres. Considering the support accuracy requirement, the total price and power consumption of sensors and intrinsically safe technology requirements in the automatic coal mining working face, we argue that the designed support attitude sensing system is more suitable for the special application scenarios with a large demand for the number of sensors and complex environment. Therefore, it can be directly applied to measure the support attitude of hydraulic support canopy, and is also more suitable for the coal mine special environment navigation applications of other mobile equipment like mine inspection robots, support robots, driving robots, and so on.

#### **5. Industrial Experiment and Application**

In this section, our support attitude sensing system based on the proposed GD-UKF is applied in the automatic coal mining working face to test practical performance, as shown in Figure 9. The industrial experiment and application were tested at the 13230 coal mining working face of Gengcun Mine of Yima Coal Industrial Group Co., Ltd. (Yima, China). The intrinsically safe micro inertial sensor was installed on the canopy of hydraulic support to measure support attitude of hydraulic support. The attitude angle information was transmitted to the remote monitoring center by the network switch,

and then the support attitude of hydraulic support was adjusted by the support attitude controller to satisfy the support requirements of coal mining.

**Figure 9.** Industrial application of the support attitude sensing system: (**a**) application environment; (**b**) remote monitoring center.

The initial pitch angle of the hydraulic support canopy was 3.5◦, then the angle was adjusted to −2.6◦ and remained stationary. In this process, the roll angle fluctuated and finally returned to the initial angle of 10◦ which was the inclination of the 13,230 coal mining working face, and the yaw angle was also changed from 6.9◦ to 3.9◦. The industrial application results using the support attitude sensing system are shown in Figure 10. The pitch and roll estimation errors of support sensing system base GD-UDF were less than 1◦, which could meet the requirements of the automatic coal mining working face. In order to verify the feasibility and superiority of the proposed algorithm for magnetic disturbances filtering, an industrial experiment of yaw angle was carried out, and the test results were shown in Figure 11. However, we could obviously find out that the yaw angle had a larger variation range than the other two angles and the yaw estimation error was relatively large, as shown in Figure 11a, the measurement error of original data was about 5◦ and that of yaw angle based on GD-UKF was less than 2◦. The reason is that the high-power equipment, such as a shearer in the automatic coal mining working face, can generate a certain complex external magnetic field interference, which has a great influence on the yaw estimation of hydraulic support and is difficult to eliminate.

**Figure 10.** Industrial application results of support attitude sensing system-based GD-UKF: (**a**) attitude angle estimation; (**b**) estimation error.

**Figure 11.** Estimation results of yaw angle tested in industrial experiment: (**a**) yaw angle estimation based on three methods; (**b**) estimation error of three methods.

The industrial application has verified that the support attitude estimation of hydraulic support with our support attitude sensing system using gradient descent algorithm-optimized quaternion-based unscented Kalman filter for attitude solution has high measurement accuracy, but the movement conditions of the automatic coal mining working face are complex, especially the magnetic disturbance. Therefore, the experiments for improving the yaw angle estimation algorithm need to be carried out under more complicated conditions.

#### **6. Conclusions and Future Work**

In order to tackle the problem of support attitude measurement of hydraulic support, this paper proposes a support attitude sensing system with the feature of intrinsic safety based on the low-cost MEMS-IMU and STM32MP157. The proposed gradient descent algorithm-optimized quaternion-based unscented Kalman filter makes full use of the characteristics of complementation of gyroscope, accelerometer and magnetometer to improve the accuracy of attitude estimation. In the proposed method, the gradient descent algorithm is employed to tune the process noise covariance for optimizing the state forecasting estimation. To verify the attitude estimation performance of the designed system and the proposed approach, an experiment was carried out and some analysis and comparisons were conducted. The experiment analysis results show that the support attitude sensing system based on the proposed approach has accurate attitude estimation performance with a dynamic measurement precision better than 1◦ and a static measurement precision better than 0.5◦. Finally, the industrial experiment and application for estimating support attitude of hydraulic support in the automatic coal mining working face were performed to validate the practicability and feasibility of the designed support attitude sensing system.

However, the error elimination of support attitude sensing system is still a challenge, especially eliminating the estimation error of yaw angle in the environment of complex magnetic disturbance. In future studies, the authors will focus on algorithm research on error elimination in the complex magnetic perturbation environment. These studies may include updating the covariance of the algorithm in combination with other potential interference factors affecting support attitude estimation to further improve the estimation accuracy. In addition, the tuning parameter intelligent algorithm based on the maximum likelihood estimation of mathematic statistics is also worth further research, and the covariance is approximately solved based on the statistical model trained by machine learning with existing data.

**Author Contributions:** Methodology, X.L.; software, D.W.; formal analysis, X.L.; data curation, C.T. and H.Y.; writing—original draft preparation, X.L.; writing—review and editing, Z.W. and L.S.; visualization, L.S.; supervision, Z.W. and C.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** The supports of China Postdoctoral Science Foundation (No.2019M661974) and the Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions in carrying out this research are gratefully acknowledged.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
