Research on the Shearer Positioning Method Based on the MEMS Inertial Sensors/Odometer Integrated Navigation System and RTS Smoother

Abstract

The shearer positioning method with an inertial measurement unit and the odometer is feasible in the longwall coal-mining process. However, the positioning accuracy will continue to decrease, especially for the micro-electromechanical inertial measurement unit (MIMU). In order to further improve the positioning accuracy of the shearer without adding other external sensors, the positioning method of the Rauch-Tung-Striebel (RTS) smoother-aided MIMU and odometer is proposed. A Kalman filter (KF) with the velocity and position measurements, which are provided by the odometer and closing path optimal estimation model (CPOEM), respectively, is established. The observability analysis is discussed to study the possible conditions under which the error states of KF can be estimated. A RTS smoother with the above-mentioned KF as the forward filter is built. Finally, the experiments of simulating the movement of the shearer through a mobile carrier were carried out, with a longitudinal movement distance of 44.6 m and a lateral advance distance of 1.2 m. The results show that the proposed method can effectively improve the positioning accuracy. In addition, the odometer scale factor and mounting angles can be estimated in real time.

1. Introduction

Automated mining based on a longwall face has shown significant potential to improve mining productivity, increase personnel safety, and secure environmental sustain-ability [1]. As shown in Figure 1, the longwall face equipment include a shearer, some hydraulic supports, and an armored face conveyor (AFC). The shearer rides on the AFC to cut coal back and forth. The AFC provides the running track for the shearer while transporting the coal. The hydraulic supports not only support the roof but also push the AFC towards the coal seam. The position of the shearer is directly related to the control of the AFC and hydraulic supports [2]. Hence, the positioning of the shearer is the key technology to realize automated mining. The inertial measurement unit, which contains a 3D inertial sensor, is widely used to estimate the position of the shearer due to its high reliability and autonomy [3]. The micro-electromechanical inertial measurement unit (MIMU) is especially favored by the mine engineers with the advantages of low cost and small size. However, the free inertial position error can grow quickly over time due to the drifting of the inertial devices [1], which includes repeatability biases, the slow-varying drifts, and the fast-varying drifts [4,5], and thus the integrated navigation mode with MIMU as the core component becomes a better choice. The MIMU/Global Positioning System (GPS) integrated navigation system is widely used as a conventional and low-cost positioning method. Unfortunately, GPS cannot be used in underground environments. The positioning system composed by ultra-wideband range measurements can be used in GPS-denied environments [6,7]. The susceptibility of the ultra-wideband system to occlusion will affect the stability of the integrated system. The MIMU/Doppler radar sensor has been successfully applied to the positioning of a continuous miner [8]. Unlike continuous coal mining, the surrounding environment of a longwall shearer is more complicated, which is detrimental to the accuracy of the Doppler radar sensor. The above-mentioned auxiliary methods need to exchange information with the external environment, and the positioning accuracy is inevitably affected by the environment. Exploring autonomous and robust auxiliary technology has become the first problem to be solved.

Zero-velocity update (ZUPT)-aided MIMU is a simple and robust integrated strategy [9,10]. However, it requires short and frequent stops. The motion constraint-aided MIMU ZUPT method can reduce the number of stops for simple ZUPT correction [11,12], but its accuracy cannot meet the actual demand of shearer positioning. The odometer is regarded as one of the most potentially useful autonomous speed sensors for land vehicles, and the system with the MIMU and odometer has been proven to be autonomous, robust, and feasible for shearer positioning [13,14,15]. The dead reckoning (DR) with Euler angles provided by MIMU and velocity provided by the odometer is one of the methods to estimate the position of the shearer. The inertial navigation attitude error, especially the heading angle error, is the main error source precluding further improvement of DR [1]. The motion constraint [16] and closing path optimal estimation model (CPOEM) [1] are used to improve the accuracy of DR. However, these two methods can only slow down the divergence of the position error. The reason is that DR, as an open-loop structure, cannot prevent the divergence of inertial navigation attitude errors. Establishing a Kalman filter (KF) based on MIMU and the odometer is another information fusion method, which is widely used in the field of conventional land navigation [17,18,19,20]. The advantages of this method lie in the closed-loop correction of the MIMU attitude and the real-time estimation of the inertial device parameters. The effect of correction and estimation has a lot to do with the steering maneuver of the vehicle. In the case of only MIMU and the odometer, the small and slow steering maneuver of the longwall shearer is difficult to achieve the ideal positioning effect [20,21]. Therefore, we plan to build a KF based on MIMU, odometer, motion constraint, and CPOEM to achieve better positioning.

The position error of the MIMU/odometer integrated system repeats “increase–decrease” changes along with the reciprocating operation of the shearer [22], which fits the typical case of the Rauch-Tung-Striebel (RTS) smoother. RTS smoothing is a technology that uses all observation information in a certain time interval to re-estimate the state based on the KF algorithm [23]. The RTS smoothing algorithm is one of the core technologies of the position and orientation system, which can enhance the ability of the position and orientation system to be free from external disturbance [24,25,26,27]. Since the estimation accuracy of the RTS smoothing algorithm is superior to that of filtering, it is often used as a reference for post-analysis of the integrated navigation system [28,29]. In the above applications, RTS smoothing works in off-line mode, which not only requires a large storage space, but also limits its use in scenarios with real-time requirements. An on-line smoothing method can overcome the above limitations and has been successfully used in the MIMU/GPS integrated navigation system [30] and pedestrian navigation system [31]. The core idea of on-line RTS smoothing can be summarized as smoothing is executed immediately after a certain time window to achieve a near real-time application effect. Obviously, on-line RTS smoothing technology is more suitable for the needs of longwall mining. Therefore, we propose a positioning method based on MIMU, odometer, motion constraint, CPOEM, and on-line RTS smoother, whose contributions and benefits can be summarized as follows:

(1) The RTS smoother is introduced into the MIMU and odometer integrated system, which can improve the shearer positioning accuracy without additional external sensors.

(2) Observability analysis of position measurement is added to theoretically provide the estimation conditions of the main error states on the basis of the previous work on velocity measurement.

(3) The mounting angles between the MIMU frame and the odometer frame can be estimated in real time, avoiding the tedious pre-calibration process.

2. Mathematical Models of Velocity and Position

The coordinate systems involved in this paper are shown in Figure 2, which are defined as follows: b-frame, the MIMU frame, which originates at the sensitive center of the MIMU, with the axes pointing to the AFC advance direction (right), shearer moving direction (forward), and upward; m-frame, the odometer frame, whose axes point right, forward, and upward; n-frame, the local-level east–north–up coordinate.

The aim of the strapdown inertial navigation system (SINS) alignment process is to determine the transition matrix from b-frame to n-frame, denoted by ${\mathit{C}}_b^n$. One of the purposes of MIMU and odometer joint calibration is to determine the mounting angles between m-frame and b-frame, which can be expressed by a vector, $\mathit{\alpha }$. The direction cosine matrix from b-frame to m-frame is denoted by ${\mathit{C}}_b^m$.

2.1. Measured Velocity Model

The moving speed of the shearer is denoted by $v_y^m$, which can be measured by the odometer. Taking into account the scale factor error, $\delta k_D$, and measurement noise, $w_{\mathrm{OD}}$, of the odometer, the actual output, ${\tilde{v}}_y^m$, of the odometer can be expressed as [20,21,32]:

$${\tilde{v}}_{y/\mathrm{OD}}^m=(1+\delta k_D)v_y^m+w_{\mathrm{OD}}$$

(1)

According to the motion constraint, there is no sideslip along the AFC advance direction and no motion normal to the AFC under ideal conditions, so the velocities along the $x_m$ axis and $z_m$ axis are regarded as zero [11]. The ideal velocity of the shearer in m-frame can be expressed as:

$${\mathit{v}}^m={[\mathbf{0} {\mathit{v}}_y^m \mathbf{0}]}^{\mathrm{T}}$$

(2)

Integrating (1) and (2), the measured velocity model in m-frame is:

$${\tilde{v}}_{\mathrm{OD}/\mathrm{MC}}^m=(1+\delta k_D)v^m+w_{\mathrm{OD}/\mathrm{MC}}$$

(3)

where $w_{\mathrm{OD}/\mathrm{MC}}$ is the noise vector, defined as:

$$w_{\mathrm{OD}/\mathrm{MC}}={[w_{\mathrm{MC},x} w_{\mathrm{OD}} w_{\mathrm{MC},z}]}^{\mathrm{T}}$$

(4)

where $w_{\mathrm{MC},x}$ and $w_{\mathrm{MC},z}$ are the motion constraint noise along the $x_m$ axis and the $z_m$ axis, respectively.

2.2. Measured Position Model

A typical shearer operation process is to repeat “straight cutting–oblique cutting–reverse straight cutting” to form a closed path, which is shown in Figure 3 [1]. The shearer runs in the order of A–B–C–D–E–F–G–H–I–J. Assuming that the shearer moves from the right to the left of Figure 3 during the first cutting cycle, we define the end corresponding to A, F, and I as the near end of the longwall face, and the corresponding end of the B, E, and J as the far end of the longwall face. The lengths of j − 1, j, and j + 1 cutting cycles, corresponding to the lengths of AB, EF, and IJ, respectively, are the same. The advance distance between two adjacent cutting cycles is also the same, denoted by d, which can be measured by the displacement sensor fixed in the push arm of the hydraulic support. Some points with Δ interval in each cutting cycle are selected as optimal points for information fusion of KF. The ideal advance displacement of an optimal point can be expressed as:

$${\mathit{D}}_{i,j-1|j}^m={[d 0 0]}^{\mathrm{T}},$$

(5)

where ${\mathit{D}}_{i,j-1|j}^m$ represents the advance displacement of point i from the j − 1 cutting cycle to the j cutting cycle.

According to the CPOEM principle and the longwall mining process [1], the positions of the optimal points during the next cutting cycle can be predicted through the positions of the current cutting cycle and the advance displacement. The position of point i in the j + 1 cutting cycle is expressed as:

$${\mathit{p}}_{i,j+1}={\mathit{p}}_{i,j}+\mathit{C}{\mathit{C}}_b^n{\mathit{C}}_m^b{\mathit{D}}_{i,j|j+1}^m$$

(6)

where ${\mathit{p}}_{i,j}$ and ${\mathit{p}}_{i,j+1}$ denote the positions of point i in the j and j + 1 cutting cycles, respectively, which are expressed in the form of longitude, $\lambda$, latitude, $L$, and height, $h$, and matrix $\mathit{C}$, whose function is to convert the position increment in n-frame into the form of longitude, latitude, and height, is:

$$\mathit{C}=\left[\begin{array}{ccc}\sec L/(R_N+h)& 0& 0\\ 0& 1/(R_M+h)& 0\\ 0& 0& 1\end{array}\right]$$

(7)

where $R_N$ and $R_M$ denote the transverse and meridian radius of curvature, respectively, and they are the parameters used to describe the earth ellipsoid model, which are regarded as constant values in this paper.

Figure 3. Closed path of shearer during the longwall mining.

Figure 3. Closed path of shearer during the longwall mining.

The initial positions of the optimal points are provided in the first cutting cycle. The specific process is as follows:

• The two optimal points corresponding to both ends of the longwall face can be accurately measured in advance, as mentioned in [16,33].
• To avoid additional surveying and mapping, the initial positions of the remaining optimal points can be obtained using the position estimates of the integrated system in the first cutting cycle. This initial value assignment method is mentioned in [1].

Taking into account the existence of measurement noise, the measured position, ${\tilde{\mathit{p}}}_{i/\mathrm{CPOEM}}$, of point i predicted by CPOEM can be approximated as:

$${\tilde{\mathit{p}}}_{i/\mathrm{CPOEM}}\approx {\mathit{p}}_i+{\mathit{w}}_{i/\mathrm{CPOEM}}$$

(8)

where ${\mathit{w}}_{i/\mathrm{CPOEM}}$ is the CPOEM noise vector of point i.

3. Integrated Navigation and RTS Smoother Models

3.1. Error State Equation of Integrated Navigation System

The mounting angles between m-frame and b-frame are inevitable, even if the pre-calibration process is carried out. The residual installation errors after calibration can be regarded as random constants, denoted by $\delta \mathit{\alpha }={[\delta {\alpha }_x \delta {\alpha }_y \delta {\alpha }_z]}^{\mathrm{T}}$. The scale factor error of the odometer can also be regarded as a random constant. Thus, the following equation can be obtained:

$$\{\begin{array}{c}\delta \dot{\mathit{\alpha }}={\mathbf{0}}_{3\times 1}\\ \delta {\dot{k}}_D=0\end{array}$$

(9)

where ${\mathbf{0}}_{a\times b}$ is a $a\times b$ zero matrix.

The direction cosine matrix, ${\tilde{\mathit{C}}}_b^m$, with the residual error $\delta \mathit{\alpha }$ satisfies:

$${\tilde{\mathit{C}}}_b^m={\mathit{C}}_b^m[\mathit{I}+(\delta \mathit{\alpha }\times )]$$

(10)

where $\mathit{I}$ is a third-order unit matrix and $(\delta \mathit{\alpha }\times )$ denotes the skew symmetric matrix of $\delta \mathit{\alpha }$ with 3 rows and 3 columns.

When the mounting angles, $\mathit{\alpha }$, are controlled within a small range by a precise mechanical installation, the equation ${\tilde{C}}_b^m=\mathit{I}$ can be obtained directly without the need to perform the pre-calibration process.

Some related error state vectors of the SINS satisfy the following equation:

$$\{\begin{array}{c}{\tilde{\mathit{v}}}_{\mathrm{SIN}\mathrm{S}}^n=v^n+\delta {\mathit{v}}^n\\ {\tilde{\mathit{C}}}_n^b={\mathit{C}}_n^b[\mathit{I}+({\mathit{\phi }}^n\times )]\\ {\tilde{\mathit{p}}}_{\mathrm{SIN}\mathrm{S}}=\mathit{p}+\delta \mathit{p}\end{array}$$

(11)

where ${\mathit{v}}^n$, ${\mathit{C}}_n^b$, and $\mathit{p}$ are the true velocity, attitude matrix, and position, respectively, and ${\tilde{\mathit{v}}}_{\mathrm{SIN}\mathrm{S}}^n$, ${\tilde{\mathit{C}}}_n^b$, and ${\tilde{\mathit{p}}}_{\mathrm{SIN}\mathrm{S}}$ are the error-contaminated velocity, attitude matrix, and position, respectively, calculated by the SINS. $\delta v^n={[\delta v_E \delta v_N \delta v_U]}^{\mathrm{T}}$ and $\delta \mathit{p}$ are the velocity errors and position errors of the SINS, respectively, ${\mathit{\phi }}^n={[{\phi }_E {\phi }_N {\phi }_U]}^{\mathrm{T}}$ is the misalignment angles of ${\tilde{\mathit{C}}}_n^b$ in n-frame, and the subscripts E, N, and U are the east, north, and up directions in n-frame.

A 19-dimensional error state vector is defined as:

$$\mathit{x}(t)={\left[{(\delta v^n)}^{\mathrm{T}} {({\mathit{\phi }}^n)}^{\mathrm{T}} {(\delta \mathit{p})}^{\mathrm{T}} {({\mathit{\epsilon }}^b)}^{\mathrm{T}} {({\mathbf{\nabla }}^b)}^{\mathrm{T}} {(\delta \mathit{\alpha })}^{\mathrm{T}} \delta k_D\right]}^{\mathrm{T}}$$

(12)

where ${\mathit{\epsilon }}^b$ and ${\mathbf{\nabla }}^b$ are the gyro and accelerometer biases, respectively.

Taking into account the low-speed motion characteristics of the shearer, the velocity-related terms in the SINS error equation can be ignored. Then, the SINS error equation can be simplified to:

$$\{\begin{array}{c}\begin{array}{c}\delta {\dot{\mathit{v}}}^n=-2{\mathit{\omega }}_{ie}^n\times \delta {\mathit{v}}^n+{\mathit{f}}^n\times {\mathit{\phi }}^n+{\mathit{C}}_b^n{\mathbf{\nabla }}^b\\ {\dot{\mathit{\phi }}}^n=-{\mathit{\omega }}_{ie}^n\times {\mathit{\phi }}^n-{\mathit{C}}_b^n{\mathit{\epsilon }}^b\end{array}\\ \delta \dot{\mathit{p}}=\mathit{C}\delta {\mathit{v}}^n\\ \begin{array}{c}{\dot{\mathit{\epsilon }}}^b={\mathbf{0}}_{3\times 1}\\ {\dot{\mathbf{\nabla }}}^b={\mathbf{0}}_{3\times 1}\end{array}\end{array}$$

(13)

where ${\mathit{\omega }}_{ie}^n$ denotes the rotation rate vector of the earth and ${\mathit{f}}^n={\mathit{C}}_b^n{\mathit{f}}^b$, in which ${\mathit{f}}^b$ is the specific force measured by the accelerometers.

The error state equation of the integrated navigation system can be expressed as:

$$\dot{\mathit{x}}(t)=\underset{\mathit{F}(t)}{\underbrace{\left[\begin{array}{cc}{\mathit{F}}_{\mathrm{SIN}\mathrm{S}}& {\mathbf{0}}_{15\times 4}\\ {\mathbf{0}}_{4\times 15}& {\mathbf{0}}_{4\times 4}\end{array}\right]}}\mathit{x}(t)+\mathit{w}(t)$$

(14)

where $\mathit{w}(t)$ is the noise vector of the integrated system and the 15 × 15 transition matrix ${\mathit{F}}_{\mathrm{SIN}\mathrm{S}}$ is denoted by:

$${\mathit{F}}_{\mathrm{SIN}\mathrm{S}}=\left[\begin{array}{ccccc}-2({\mathit{\omega }}_{ie}^n\times )& ({\mathit{f}}^n\times )& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathit{C}}_b^n\\ {\mathbf{0}}_{3\times 3}& -({\mathit{\omega }}_{ie}^n\times )& {\mathbf{0}}_{3\times 3}& -{\mathit{C}}_b^n& {\mathbf{0}}_{3\times 3}\\ \mathit{C}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{6\times 3}& {\mathbf{0}}_{6\times 3}& {\mathbf{0}}_{6\times 3}& {\mathbf{0}}_{6\times 3}& {\mathbf{0}}_{6\times 3}\end{array}\right]$$

(15)

The discretized error state equation corresponding to (14) is:

$${\mathit{X}}_k={\mathit{\Phi }}_{k,k-1}{\mathit{X}}_{k-1}+{\mathit{W}}_{k-1},$$

(16)

where ${\mathit{X}}_k$ and ${\mathit{X}}_{k-1}$ are the discretized error state vectors at $t_k$ and $t_{k-1}$, respectively, ${\mathit{W}}_{k-1}$ denotes the system noise matrix at $t_{k-1}$, and ${\mathit{\Phi }}_{k,k-1}$ is the discretized state transition matrix that satisfies the following equation:

$${\mathit{\Phi }}_{k,k-1}{=\mathrm{e}}^{TF(t_{k-1})}$$

(17)

where $T$ denotes the filtering period.

3.2. Measurement Equations of Velocity and Position

The velocity, ${\tilde{v}}_{\mathrm{SIN}\mathrm{S}}^m$, in m-frame calculated by the SINS can be expressed as:

$$\begin{array}{lll}{\tilde{\mathit{v}}}_{\mathrm{SIN}\mathrm{S}}^m& =& {\tilde{\mathit{C}}}_b^m{\tilde{\mathit{C}}}_n^b{\tilde{\mathit{v}}}_{\mathrm{SIN}\mathrm{S}}^n\\ & \approx & v^m+{\mathit{C}}_b^m{\mathit{C}}_n^b\delta {\mathit{v}}^n-{\mathit{C}}_b^m{\mathit{C}}_n^b({\mathit{v}}^n\times ){\mathit{\phi }}^n\\ & & -{\mathit{C}}_b^m({\mathit{v}}^b\times )\delta \mathit{\alpha }\end{array}$$

(18)

where ${\mathit{v}}^b={\mathit{C}}_n^b{\mathit{v}}^n$.

We find the difference between (3) and (18) to obtain the velocity measurement equation as:

$${\mathit{z}}_v={\tilde{\mathit{v}}}_{\mathrm{SIN}\mathrm{S}}^m-{\tilde{\mathit{v}}}_{\mathrm{OD}/\mathrm{MC}}^m={\mathit{H}}_v\mathit{x}(t)+{\mathit{w}}_{\mathrm{OD}/\mathrm{MC}}$$

(19)

where the measurement matrix, ${\mathit{H}}_v$, of velocity is expressed as:

$${\mathit{H}}_v=\left[\begin{array}{ccc}\begin{array}{cc}\begin{array}{cc}{\mathit{C}}_b^m{\mathit{C}}_n^b& -{\mathit{C}}_b^m{\mathit{C}}_n^b({\mathit{v}}^n\times )\end{array}& {\mathbf{0}}_{3\times 9}\end{array}& -{\mathit{C}}_b^m({\mathit{v}}^b\times )& -{\mathit{v}}^m\end{array}\right]$$

(20)

Similarly, the position measurement equation can be expressed as:

$${\mathit{z}}_p={\tilde{\mathit{p}}}_{\mathrm{SIN}\mathrm{S}}-{\tilde{\mathit{p}}}_{\mathrm{CPOEM}}={\mathit{H}}_p\mathit{x}(t)+{\mathit{w}}_{\mathrm{CPOEM}}$$

(21)

where the measurement matrix, ${\mathit{H}}_p$, of the position is given by:

$${\mathit{H}}_p=\left[\begin{array}{ccc}{\mathbf{0}}_{3\times 6}& \mathit{I}& {\mathbf{0}}_{3\times 10}\end{array}\right]$$

(22)

The update of the position measurement depends on the mileage calculated by the odometer, which means that the update period of the position measurement is an integer multiple of that of the velocity measurement. Therefore, the overall measurement model includes two forms: simultaneous velocity and position measurement update and separate velocity measurement update. The specific update conditions and model equations are as follows.

If the shearer is not in the first cutting cycle and the position measurement is judged to be valid by the mileage calculated by the odometer, then:

$$\{\begin{array}{c}\mathit{z}=\left[\begin{array}{c}{\mathit{z}}_v\\ {\mathit{z}}_p\end{array}\right]\\ \mathit{H}=\left[\begin{array}{c}{\mathit{H}}_v\\ {\mathit{H}}_p\end{array}\right]\end{array}$$

(23)

where $\mathit{z}$ and $\mathit{H}$, respectively, represent the overall measurement vector and measurement matrix.

Otherwise:

$$\{\begin{array}{c}\mathit{z}={\mathit{z}}_v\\ \mathit{H}={\mathit{H}}_v\end{array}$$

(24)

The discretized measurement equation is:

$${\mathit{Z}}_k={\mathit{H}}_k{\mathit{X}}_k+{\mathit{V}}_k$$

(25)

where ${\mathit{Z}}_k$, ${\mathit{H}}_k$, and ${\mathit{V}}_k$ are the discretized measurement vector, the discretized measurement matrix, and the measurement noise sequence at $t_k$, respectively.

3.3. Observability Analysis of Integrated System

The main purpose of observability analysis is to study the observability and estimatable conditions of the error state. An observability analysis method directly relies on the state and measurement equations to investigate the observability.

3.3.1. Observability Analysis Based on the Velocity Measurement

The observability analysis based on the velocity measurement has been discussed in detail in our previous work [34], so here, we directly summarize the observability conclusions about the velocity measurement with the motion characteristics of the longwall shearer.

• The velocity errors, $\delta v^n$, are observable, and the estimation accuracy is related to the estimation degree of other error states.
• The position errors, $\delta \mathit{p}$, are unobservable, but the estimation accuracy will still be improved with the effective estimation of the velocity errors.
• The acceleration and deceleration process of the shearer is the premise of exciting the error states $\delta {\alpha }_x$, $\delta {\alpha }_z$, and $\delta k_D$, which contribute to the positioning errors. The error $\delta {\alpha }_y$ is unobservable.
• The error ${\nabla }_z^b$ is observable. The separation of ${\nabla }_x^b$ and ${\nabla }_y^b$, and the distinction of ${\epsilon }_x^b$ and ${\epsilon }_y^b$, improving the estimation accuracy of ${\mathit{\phi }}^n$, depend on the turning motion of the shearer. The azimuth error, ${\phi }_U$, is directly related to the lateral positioning error, thus, restricting the estimation accuracy of error $\delta {\alpha }_z$.

It can be concluded that frequent turning of the shearer is necessary to improve the estimation accuracy of the error states. However, limited by working conditions, the longwall shearer has very few steering maneuvers. Therefore, it is difficult to achieve high shearer positioning accuracy using only the integrated system of MIMU and the odometer.

3.3.2. Observability Analysis Based on the Position Measurement

In (21), the measurement values are constructed by the position errors. Therefore, the position errors, $\delta \mathit{p}$, are observable. The estimation accuracy of $\delta \mathit{p}$ depends on the position accuracy of the optimal points predicted by CPOEM. The vectors that make up the measurement values have the same initial position errors, so the initial position errors cannot be estimated.

Taking the time derivative of (21) obtains:

$${\dot{\mathit{z}}}_p\approx \mathit{C}\mathit{\delta }{\mathit{v}}^n$$

(26)

Equation (26) indicates that the velocity errors, $\mathit{\delta }{\mathit{v}}^n$, are observable, and the estimation accuracy is determined by the position errors.

The terms related to ${\mathit{v}}^n$ in the mathematical models can be ignored due to the low-speed characteristics of the shearer. Therefore, the time derivative of (26) can be simplified as:

$$\begin{array}{l}({\mathit{C}}^{-1}){\ddot{\mathit{z}}}_p\approx ({\mathit{f}}^n\times ){\mathit{\phi }}^n-2({\mathit{\omega }}_{ie}^n\times )\mathit{\delta }{\mathit{v}}^n+{\mathbf{\nabla }}^n\\ =\left[\begin{array}{c}-f_U{\phi }_N+f_N{\phi }_U+2{\mathit{\omega }}_{ie}\sin L\delta v_N-2{\mathit{\omega }}_{ie}\cos L\delta v_U+{\nabla }_E\\ f_U{\phi }_E-f_E{\phi }_U-2{\mathit{\omega }}_{ie}\sin L\delta v_E+{\nabla }_N\\ -f_N{\phi }_E+f_E{\phi }_N+2{\mathit{\omega }}_{ie}\cos L\delta v_E+{\nabla }_U\end{array}\right]\end{array}$$

(27)

where ${\mathit{f}}^n={[f_E f_N f_U]}^{\mathrm{T}}$, ${\mathbf{\nabla }}^n={[{\nabla }_E {\nabla }_N {\nabla }_U]}^{\mathrm{T}}$, and ${\mathit{\omega }}_{ie}$ is the rotation rate of the earth.

The acceleration and deceleration process of the shearer makes at least one of $f_E$ and $f_N$ non-zero, which means that the coefficient of ${\phi }_U$ will not be 0 in (27). That is to say, the acceleration and deceleration process can improve the observability of the azimuth error, ${\phi }_U$.

Let us assume a high-precision position predicted by CPOEM is obtained. Then, the errors $\delta \mathit{p}$ and $\delta v^n$ will be accurately estimated. The estimation accuracy of the azimuth error, ${\phi }_U$, depends on the separating degree from the errors ${\phi }_E$, ${\phi }_N$, ${\nabla }_E$, and ${\nabla }_N$. Although the state errors $\delta k_D$, $\delta {\alpha }_x$, and $\delta {\alpha }_z$ are not directly observable in (21), (26), and (27), their accuracy will be improved with the estimation of error states $\delta \mathit{p}$ and ${\mathit{\phi }}^n$. Similarly, the estimation accuracy of ${\mathit{\epsilon }}^b$ will also be improved with the estimation of other state errors.

3.4. RTS Smoothing

RTS smoothing uses all the measurement values obtained in a time interval to estimate the error states at every epoch in this interval. A typical RTS smoothing process is divided into two steps: forward filtering and backward smoothing.

The forward filtering is obtained through standard KF, whose basic equation for discretization can be expressed as:

$$\{\begin{array}{c}\begin{array}{c}{\mathit{X}}_{k,k-1}={\mathit{\Phi }}_{k,k-1}{\mathit{X}}_{k-1}\\ {\mathit{P}}_{k,k-1}={\mathit{\Phi }}_{k,k-1}{\mathit{P}}_{k-1}{\mathit{\Phi }}_{k,k-1}^{\mathrm{T}}+{\mathit{Q}}_{k-1}\end{array}\\ {\mathit{K}}_k={\mathit{P}}_k{\mathit{H}}_k^{\mathrm{T}}{\mathit{R}}_k^{-1}\\ \begin{array}{c}{\mathit{X}}_k={\mathit{X}}_{k,k-1}+{\mathit{K}}_k({\mathit{Z}}_k-{\mathit{H}}_k{\mathit{X}}_{k,k-1})\\ {\mathit{P}}_k=(\mathit{I}-{\mathit{K}}_k{\mathit{H}}_k){\mathit{P}}_{k,k-1}\end{array}\end{array}$$

(28)

where ${\mathit{X}}_{k,k-1}$ and ${\mathit{P}}_{k,k-1}$ are the one-step predicted states and covariance at $t_k$, calculated from the information at $t_{k-1}$, ${\mathit{P}}_{k-1}$ and ${\mathit{P}}_k$ denote the state estimate covariance at $t_{k-1}$ and $t_k$, respectively, and ${\mathit{Q}}_{k-1}$ is the variance matrix of the system noise sequence ${\mathit{W}}_{k-1}$ at $t_{k-1}$. ${\mathit{K}}_k$ represents the filter gain matrix at $t_k$, and ${\mathit{R}}_k$ is the variance matrix of the measurement noise sequence ${\mathit{V}}_k$ at $t_k$.

The forward filtering in this paper refers to the KF based on the velocity and position measurements mentioned above. In the forward filtering process, it is required to save ${\mathit{X}}_{k,k-1}$, ${\mathit{X}}_k$, ${\mathit{\Phi }}_{k,k-1}$, ${\mathit{P}}_{k,k-1}$, and ${\mathit{P}}_k$ at every epoch in the time interval. After completing the forward filtering in the time interval, the backward smoothing is performed. The procedure of the backward smoothing is broken down into the following steps:

• Backward Smoothing Initialization

Define the time interval as [$t_j$,$t_{j+N}$], then the initialization equation is given as:

$$\{\begin{array}{c}{\mathit{X}}_{\mathrm{rts},j+N}={\mathit{X}}_{j+N}\\ {\mathit{P}}_{\mathrm{rts},j+N}={\mathit{P}}_{j+N}\end{array},$$

(29)

where ${\mathit{X}}_{j+N}$ and ${\mathit{P}}_{j+N}$ are the error state vector and the mean square error matrix of the forward filtering process at $t_{j+N}$, and ${\mathit{X}}_{\mathrm{rts},j+N}$ and ${\mathit{P}}_{\mathrm{rts},j+N}$, which are the initial values required for backward smoothing, denote the error state vector and the mean square error matrix of the backward smoothing process at $t_{j+N}$.

2.

Backward Smoothing Update

The recursive equations of the backward smoothing update are given by:

$$\{\begin{array}{c}{\mathit{X}}_{\mathrm{rts},k}={\mathit{X}}_k+{\mathit{A}}_k({\mathit{X}}_{rts,k+1}-{\mathit{X}}_{k+1,k})\\ {\mathit{P}}_{\mathrm{rts},k}={\mathit{P}}_k+{\mathit{A}}_k({\mathit{P}}_{rts,k+1}-{\mathit{P}}_{k+1,k}){\mathit{A}}_k^{\mathrm{T}}\end{array},$$

(30)

where ${\mathit{A}}_k$ is the smoother gain, determined as:

$${\mathit{A}}_k={\mathit{P}}_k{\mathit{\Phi }}_{k+1,k}^{\mathrm{T}}{\mathit{P}}_{k+1,k}^{-1},$$

(31)

Sorting out the process of KF and the RTS smoother, the flowchart is shown in Figure 4. If mod(S, Δ) = = 0, where S denotes the mileage of the shearer calculated by the odometer, the optimal point from the CPOEM data is valid and the position measurement update of KF is performed. Backward smoothing is only executed when the shearer is at both ends of the longwall face and not in the first cutting cycle. The positions of the optimal points of the next cycle will be updated immediately after the end of backward smoothing.

4. Experiments

To evaluate the performance of the proposed positioning method, experiments were carried out, as shown in Figure 5. A mobile carrier equipped with the MIMU (Xsens MTi-G-700) and the odometer simulated the movement of the shearer. The MIMU was installed on the mobile carrier through an adapter plate and the odometer was connected to its wheel. A GPS receiver with an antenna was also installed on the mobile carrier. The GPS receiver can output centimeter-level positioning results through the network differential technology. The positioning result of the GPS receiver only provides an evaluation basis for tests. The specifications of the MIMU and the initial errors are listed in

Table 1. Specifications of the MIMU and initial errors.

Initial errors	$\mathrm{Initial} \mathrm{attitude} \mathrm{errors} (°)$	[0.3;0.3;1]
	Initial velocity errors (m/s)	[0.01;0.01;0.01]
	Initial position errors (m)	[0.05;0.05;0.05]
Gyroscope	$\mathrm{Bias} \mathrm{repeatability} (°/\mathrm{h}$)	720
	$\mathrm{Random} \mathrm{walk} (°/\sqrt{\mathrm{h}}$)	0.6
Accelerometer	$\mathrm{Bias} \mathrm{repeatability} (\mathrm{mg}$)	3
	$\mathrm{Random} \mathrm{walk} (\mathrm{mg}/\sqrt{\mathrm{Hz}}$)	0.08

, which are related to the initial parameter configuration of the Kalman filter. The specifications and the initial attitude error refers to the MTI user manual [35]. The initial position error refers to the network differential positioning accuracy [36]. The initial parameters of the filter are set as Appendix A.

The mobile carrier simulated four cutting processes of the shearer with a reciprocating travel distance of 44.6 m and an advance distance of 1.2 m, as shown in Figure 6. Figure 6 is drawn by the network differential results provided by the GPS receiver. The symbols “*” and “o” denote the start and end of the trajectory, respectively.

In order to study the influence of the selection of optimal points on the estimation of error states, we set the optimal points interval, Δ, as 2, 6, and 10 m, in turn. Figure 7 shows the positioning errors of the proposed integrated system without performing RTS smoothing. A1–A4 corresponded to the time periods of the first to fourth cutting cycles, respectively. The remaining time periods corresponded to the process of the shearer processing the end coal seam. It can be seen that during the first cutting cycle, the curves of different values of Δ overlapped in the east, north, and height directions, respectively, since position measurement filtering was not performed. During the second to fourth cutting cycles, the curves of the east and north errors were smoother as Δ decreased. Figure 8 shows the positioning errors of the proposed integrated system with performing RTS smoothing. It was straightforward to see that the positioning accuracy of all three axes with performing RTS smoothing was higher than that of not performing, regardless of the value of Δ. In addition, it can be seen that there were larger burrs in the east and north directions as Δ increased. Therefore, we can choose a smaller value of Δ to obtain smoother position estimations.

Figure 9 and Figure 10, respectively, show the estimations of the accelerometer and gyro biases with and without performing RTS smoothing, while we set Δ as 2 m. It can be seen that the estimations of MIMU biases were equivalent with and without performing RTS smoothing when the estimators tended to be stable. Although the true values of MIMU biases cannot be accurately obtained, the results of the pure navigation before and after the bias compensation can reflect their estimation accuracy. We used the first 300 s of data to perform the pure navigation calculations before and after the bias compensation, and the results are shown in Figure 11. Compared with the results before compensation of the MIMU biases, the positioning accuracy was greatly improved after compensation. Therefore, the MIMU biases can be effectively estimated. The RTS smoothing technology only reduced the error fluctuations in the estimation process, and did not affect the final estimation accuracy.

The estimation results of the odometer scale factor and mounting angles are shown in Figure 12. We know that $\delta {\alpha }_y$ is always unobservable from the observability analysis. Therefore, its estimation curve is not drawn in this paper. The estimated scale factor and mounting angles, whose record order was $\delta \mathit{k}\mathit{\alpha }={[\delta k_D \delta {\alpha }_{\mathit{X}} \delta {\alpha }_z]}^{\mathrm{T}}$, with and without performing RTS smoothing, were $\delta \mathit{k}{\mathit{\alpha }}_1={[0.017 0.333°-0.307°]}^{\mathrm{T}}$ and $\delta \mathit{k}{\mathit{\alpha }}_2={[0.029 0.352° 0.093°]}^{\mathrm{T}}$, respectively. The positioning accuracy of the DR algorithm is restricted by $\delta \mathit{k}\mathit{\alpha }$. In other words, the DR results can reflect the estimation effect of these three error states. The DR navigation calculations were performed after the compensation of $\delta \mathit{k}{\mathit{\alpha }}_1$ and $\delta \mathit{k}{\mathit{\alpha }}_2$ using the data that has been compensated for the MIMU biases, and the results are shown in Figure 13. The horizontal positioning accuracy after compensation of $\delta \mathit{k}{\mathit{\alpha }}_1$ was better than that after compensation of $\delta \mathit{k}{\mathit{\alpha }}_2$, which means that the estimation accuracy of the error states $\delta k_D$ and $\delta {\alpha }_z$ with performing RTS smoothing was improved. According to the theory of observability analysis, the estimation accuracy of error states (such as $\delta k_D$ and $\delta {\alpha }_z$) can be improved by feedback when the position errors of the optimal points were reduced. Since the positioning accuracy with performing the RTS smoothing was better than that without performing (shown in Figure 7 and Figure 8), the position errors of the corresponding optimal points were also smaller. Therefore, the results of Figure 13 are consistent with the observability analysis.

The shearer is a long and narrow machine with limited space for the installation of external sensors. The small size of MEMS inertial sensors makes them a goal pursued by researchers. However, small-sized and low-cost inertial sensors often have low measurement accuracy. Therefore, we studied the impact of MIMU, which has a lower accuracy than MTi-G-700, on positioning accuracy of the mobile carrier. The idea of MIMU data generation is to superimpose errors on the original gyroscope and accelerometer data corresponding to Figure 6 to simulate lower-precision MEMS inertial sensor data. The error components of the gyroscope and accelerometer include repeatability biases, slow-varying drifts, and fast-varying drifts [4,5]. The repeatability biases can be regarded as a random constant. The slow-varying drift can be approximated as white noise due to the short correlation time of MIMU. The fast-varying drift is often abstracted as a white noise process. The white noise process is usually evaluated by random walk. In summary, the error model of the gyroscope and accelerometer can be expressed as:

$$\{\begin{array}{c}{\mathit{\epsilon }}_s^b={\mathit{\epsilon }}^b+{\mathit{\epsilon }}_w^b\\ {\mathbf{\nabla }}_s^b={\mathbf{\nabla }}^b+{\mathbf{\nabla }}_w^b\end{array}$$

(32)

where ${\mathit{\epsilon }}_s^b$ and ${\mathbf{\nabla }}_s^b$ represent the total gyroscope and accelerometer errors, respectively, ${\mathit{\epsilon }}^b$ and ${\mathbf{\nabla }}^b$ denote the constant biases of the gyroscope and accelerometer, respectively, and ${\mathit{\epsilon }}_w^b$ and ${\mathbf{\nabla }}_w^b$ are the random walk noises of the gyroscope and accelerometer, respectively.

A series of MIMU errors were designed, as listed in

Table 2. Error parameters to be accumulated.

	Gyroscope		Accelerometer
	$\mathbf{Constant} \mathbf{Bias} (°/\mathbf{h})$	$\mathbf{Random} \mathbf{Walk} (°/\sqrt{\mathbf{h}})$	$\mathbf{Constant} \mathbf{Bias} \left(\mathbf{mg}\right)$	$\mathbf{Random} \mathbf{Walk} (\mathbf{mg}/\sqrt{\mathbf{Hz}})$
Par 1	100	0	0	0
Par 2	200	0	0	0
Par 3	300	0	0	0
Par 4	0	1	0	0
Par 5	0	2	0	0
Par 6	0	3	0	0
Par 7	0	0	1	0
Par 8	0	0	2	0
Par 9	0	0	3	0
Par 10	0	0	0	1
Par 11	0	0	0	2
Par 12	0	0	0	3

. The errors in Table 2 were added to the original data corresponding to Figure 6. The position errors before and after RTS smoothing are shown in Figure 14 and Figure 15, respectively. Similar to the phenomenon in Figure 7 and Figure 8, even if the accuracy of the inertial sensors is reduced, the positioning accuracy can still be improved after performing RTS smoothing. It can be seen from Figure 15 that the east errors of Par 3 and Par 2 are much larger than those of other parameters, and the east error of Par 3 is greater than that of Par 2. This is because the constant biases of the gyroscope are not fully estimated, and the residual constant component is larger as the constant biases increase. Correspondingly, the attitude and position errors are bound to be greater. In summary, it can be concluded that the proposed method can still improve the positioning accuracy even if the low-precision inertial sensors are used, and the positioning accuracy is closely related to the sensor parameters, especially the gyroscope.

In order to further verify the performance of the method proposed in this paper, it is compared with the traditional method. The method mentioned in [1] is state-of-the-art based on the SINS and the odometer. However, this method does not have the ability to autonomously estimate the biases of inertial sensors, mounting angles between MIMU and the odometer, and the scale factor error of the odometer. Considering that the above parameters have a great influence on the traditional method, it was compensated with the parameter values obtained by the proposed method. The positioning results of the traditional method before and after the estimated parameter compensation and the method proposed in this paper are shown in Figure 16. It can be seen from Figure 16 that after the parameters are compensated, the positioning errors of the traditional method are significantly reduced. It can also be seen that the accuracy of the method proposed in this paper is better than the traditional method. In order to visually describe the accuracy improvement range of the proposed method, the spherical error probable (SEP) is calculated. The SEP is a universal evaluation method of 3D positioning accuracy [37], which is listed in

Table 3. SEP in each cutting cycle.

	Cutting Cycle	SEP (m)
Traditional method without compensation	A1	1.75
	A2	1.99
	A3	2.28
	A4	2.31
Traditional method with compensation	A1	0.21
	A2	0.29
	A3	0.33
	A4	0.17
Proposed method	A1	0.50
	A2	0.32
	A3	0.13
	A4	0.17

. A1–A4 correspond to the first to fourth cutting cycles, respectively. It can be seen from Table 3 that compared with the traditional method without compensation, the positioning accuracy of the proposed method increases by 71.43%, 83.92%, 94.30% and 92.64% in turn from A1 to A4. Since the position measurement is not performed in A1, and the optimal point position of A2 is related to that of A1, the SEP of the proposed method is larger in A1 and A2. It is significantly reduced after A3. Compared with the traditional method with compensation, the positioning accuracy of the proposed method increases by 60.61% in A3 and is equivalent in A4. The above phenomenon not only proves the superiority of the proposed method in positioning accuracy, but also further shows that the method proposed in this paper can effectively estimate some parameters, such as biases of the MIMU, odometer scale factor, etc.

5. Conclusions

This paper proposed a positioning method of the shearer based on an integrated system and RTS smoothing technology. Performing RTS smoothing on the basis of the Kalman filter is a major feature of this paper. An experiment was carried out to verify the performance of the proposed positioning method. The experimental results showed that the positioning accuracy after performing RTS smoothing was significantly improved, which was closely related to the sensor parameters, and the estimatable ability of some error states was improved. In addition, a comparison with traditional methods was also carried out. The result shows that the positioning accuracy of the proposed method can be improved by at least 60.61%.

According to the existing theories and the experimental results in this paper, it can be seen that RTS smoothing technology has a significant improvement effect on the jump phenomenon of the error states. Therefore, not only the position measurement, but other excellent measurement information that can be captured may also cause the error state to jump. At this time, RTS smoothing can still play an important role, which can be further studied.