Research on the Shearer Positioning Method Based on SINS and LiDAR with Velocity and Absolute Position Constraints

: Accurate positioning of the shearer with a strapdown inertial navigation system (SINS) is the key technology to realize the automation of the longwall face. Unfortunately, the existing positioning methods have a strong dependence on the attitude accuracy of the SINS. The position errors gradually increase with the drift of the SINS attitude. To reduce the dependence on the SINS attitude and further increase the shearer positioning accuracy, this paper proposes a positioning method based on SINS and light detection and ranging (LiDAR) with velocity and absolute position constraints. A Kalman ﬁlter (KF) model based on these constraints was established. Simulation analysis shows that the attitude calibration between the shearer body, SINS and LiDAR, and the initial attitude alignment of the SINS are the keys to determining the shearer positioning accuracy. Even if there are small horizontal bends in the running track of the shearer and the features have small horizontal errors, an excellent positioning effect can still be obtained. In addition, four cutting processes were simulated with a reciprocating travel of 44.6 m and an advance distance of 1.2 m. Compared with the relative positioning method, the positioning accuracy of the proposed method was improved by 37%, 63%, 76%, and 69% from the ﬁrst to the fourth cutting cycle, respectively, calculated by spherical error probable (SEP) values, and positioning accuracy had a lower dependence on the installation deﬂection angles between the SINS, the LiDAR, and the SINS attitude accuracy.


Introduction
As a precious resource, coal still accounts for a large proportion of the world's primary energy consumption structure. However, the poor mining environment with high temperatures, high gas quantities, high ground stresses, etc., poses a huge threat to the personal safety of miners [1]. Pneumoconiosis is one of the diseases that threatens the health of miners which is widespread in major coal-producing countries in the world. In Poland, there was a total of 7340 cases of coal worker's pneumoconiosis in 2000-2017 [2], while in China the number was nearly 7949 in 2019 alone [3]. In addition, coal mine safety accidents also directly threaten the safety of miners. According to statistics from China Energy News, a total of 122 fatal accidents occurred in China's coal mines in 2020, with 225 deaths. Although the safety factor of coal mines in the United States, Australia and other countries is high, mine accidents also occasionally occur. These shocking numbers prompt us to explore safer and more efficient mining methods. Longwall mining technology has been widely used in European countries such as Poland [4] and the Czech Republic [5], where its advantages in mining have been fully proven. Automated mining and remote-control technology based on the longwall face can keep miners away from dangerous areas [6]. The longwall mining process could be realistically automated by achieving three goals: face alignment, creep control, and horizon control. Accurate three-dimensional positioning of the shearer is the premise for these goals [7]. Therefore, it is very important to study longwall shearer positioning for automated longwall mining. the longwall face, the higher the calibration accuracy of the required installation deflection angle. Similar to the integrated system of the SINS and OD, this relative positioning method also depends on the attitude accuracy of the SINS. To further improve the positioning accuracy, a positioning method of the velocity and absolute position constraint-aided SINS is proposed. A major novelty of this paper is to propose a positioning method that does not require the SINS attitude, especially the heading angle, to participate in the establishment of the measurement space model. Therefore, the positioning result is more stable and less dependent on the SINS attitude. The research on the proposed positioning method is also a continuation of our previous work in [30]. Thus, the purpose of this paper is to propose a higher positioning accuracy method based on a SINS/LiDAR integrated navigation system compared with [30]. Based on the velocity and absolute position constraints, a Kalman filter (KF) model was built. The comparison results and analysis were presented by simulations and experiment. The chief contributions and benefits of this study can be summarized as follows.
(1) An absolute position constraint was introduced to reduce the influence of the installation deflection angle between the SINS and LiDAR and the SINS attitude on the positioning accuracy of the shearer, compared with the relative positioning method. (2) A calibration method of the heading installation angle between the SINS and the LiDAR was proposed to improve the absolute position accuracy of the features. (3) The horizontal advancing displacement of the hydraulic support can be measured autonomously, and was obtained by the additional equipment in the traditional positioning method.
This paper is organized as follows. Section 2 describes the mining process of the longwall face. Section 3 presents the measurement models of velocity and the absolute position constraints. Section 4 describes the error model of the integrated navigation. Section 5 discusses the simulation analysis. Section 6 presents the experimental results and analysis. The conclusions are drawn in Section 7.

Longwall Mining Process
A longwall panel is a large rectangular block of coal, as shown in Figure 1 [30]. A typical panel is 1000-6000 m in length, and 150-350 m in width [24]. Before the extraction of a panel, two horizontal and permanent roadways are excavated, with some coal pillars left untouched to support the overlying strata [31]. The roadways, which are 5-6 m in width, are on both sides of the panel. The roadway along one side of the panel is called the maingate, and the roadway on the other side is called the tailgate [32]. At the end of the roadways, the major pieces of mining equipment are installed across the back of the panel, creating the longwall face. The coal is extracted by the longwall face. The goaf is the area that has been extracted. A longwall face, which includes a shearer, some hydraulic supports, and an armored face conveyor (AFC), carries out the tasks of coal breaking, coal loading, and coal transportation.
The workflow of the longwall face is shown in Figure 2. The shearer rides on the AFC and travels back and forth with the AFC as the track. The hydraulic supports, which are used to support the roof and push the AFC toward the coal, are connected to the AFC and arranged horizontally on the side of the AFC with equal spacing [33]. Figure 2a represents the shearer movement to the left, accompanied by the movement of the AFC and some hydraulic supports behind the shearer in the advancing direction. Figure 2b shows the static state of the shearer after oblique cutting. The process of the oblique cutting refers to the shearer cuts into the coal seam from the end through the curved section of the AFC, and stops after entering the straight section of the AFC again. Oblique cutting is the key process by which the shearer cuts into the coal seam. Figure 2c shows the reverse movement of the shearer, the purpose of which is to mine the coal on the left. Figure 2d depicts the movement of the shearer from the left end to the right, while the hydraulic supports lagging behind the shearer will move toward the coal to force the AFC to bend. Then, it returns to the state shown in Figure 2a. The workflow of the longwall face is shown in Figure 2. The shearer rides on the AFC and travels back and forth with the AFC as the track. The hydraulic supports, which are used to support the roof and push the AFC toward the coal, are connected to the AFC and arranged horizontally on the side of the AFC with equal spacing [33]. Figure 2a represents the shearer movement to the left, accompanied by the movement of the AFC and some hydraulic supports behind the shearer in the advancing direction. Figure 2b shows the static state of the shearer after oblique cutting. The process of the oblique cutting refers to the shearer cuts into the coal seam from the end through the curved section of the AFC, and stops after entering the straight section of the AFC again. Oblique cutting is the key process by which the shearer cuts into the coal seam. Figure 2c shows the reverse movement of the shearer, the purpose of which is to mine the coal on the left. Figure 2d depicts the movement of the shearer from the left end to the right, while the hydraulic supports lagging behind the shearer will move toward the coal to force the AFC to bend. Then, it returns to the state shown in Figure 2a.   The workflow of the longwall face is shown in Figure 2. The shearer rides on the AFC and travels back and forth with the AFC as the track. The hydraulic supports, which are used to support the roof and push the AFC toward the coal, are connected to the AFC and arranged horizontally on the side of the AFC with equal spacing [33]. Figure 2a represents the shearer movement to the left, accompanied by the movement of the AFC and some hydraulic supports behind the shearer in the advancing direction. Figure 2b shows the static state of the shearer after oblique cutting. The process of the oblique cutting refers to the shearer cuts into the coal seam from the end through the curved section of the AFC, and stops after entering the straight section of the AFC again. Oblique cutting is the key process by which the shearer cuts into the coal seam. Figure 2c shows the reverse movement of the shearer, the purpose of which is to mine the coal on the left. Figure 2d depicts the movement of the shearer from the left end to the right, while the hydraulic supports lagging behind the shearer will move toward the coal to force the AFC to bend. Then, it returns to the state shown in Figure 2a. In order to facilitate the research on the movement characteristics of the shearer, the trajectory of the shearer was sorted according to Figure 2. Ideally, the track of the shearer is mainly in the form of a straight line and a broken line, as shown in Figure 3. The shearer runs in the sequence of A-B-C-D-E-F-G-H-I-J-K-L-M-N. A cutting cycle refers to the process in which the shearer moves from one end of the longwall face to the other to complete a complete cut. Trajectories AB, EF, IJ, and MN correspond to the cutting cycle in Figure  3. The distance between two adjacent cutting cycles is defined as the cutting depth. The In order to facilitate the research on the movement characteristics of the shearer, the trajectory of the shearer was sorted according to Figure 2. Ideally, the track of the shearer is Remote Sens. 2021, 13, 3708 5 of 20 mainly in the form of a straight line and a broken line, as shown in Figure 3. The shearer runs in the sequence of A-B-C-D-E-F-G-H-I-J-K-L-M-N. A cutting cycle refers to the process in which the shearer moves from one end of the longwall face to the other to complete a complete cut. Trajectories AB, EF, IJ, and MN correspond to the cutting cycle in Figure 3. The distance between two adjacent cutting cycles is defined as the cutting depth. The cutting depth can be regarded as a constant, which is usually 0.8-1.2 m.
(d) In order to facilitate the research on the movement characteristics of the shearer, the trajectory of the shearer was sorted according to Figure 2. Ideally, the track of the shearer is mainly in the form of a straight line and a broken line, as shown in Figure 3. The shearer runs in the sequence of A-B-C-D-E-F-G-H-I-J-K-L-M-N. A cutting cycle refers to the process in which the shearer moves from one end of the longwall face to the other to complete a complete cut. Trajectories AB, EF, IJ, and MN correspond to the cutting cycle in Figure  3. The distance between two adjacent cutting cycles is defined as the cutting depth. The cutting depth can be regarded as a constant, which is usually 0.8-1.2 m    

System Description
As shown in Figure 4, both the SINS and the LiDAR are installed on the shearer, and the LiDAR faces the hydraulic supports. The SINS body frame, denoted by the b -frame, which is implicitly predefined by the calibrated sensitive axes of the inertial sensors, originates at the sensitive center of the SINS, point b O , with the axes pointing right, forward, and upward, respectively. The shearer frame, denoted by the m -frame, is fixed to the shearer with the m x -axis pointing right, the m y -axis pointing forward, and the m z -axis pointing upward. The LiDAR body frame, denoted by the s -frame, originates at the optical center, point s O , with the axes pointing right, forward, and upward, respectively.
The n -frame is the navigation frame, chosen as the local level east-north-up (ENU) coordinate.

System Description
As shown in Figure 4, both the SINS and the LiDAR are installed on the shearer, and the LiDAR faces the hydraulic supports. The SINS body frame, denoted by the b-frame, which is implicitly predefined by the calibrated sensitive axes of the inertial sensors, originates at the sensitive center of the SINS, point O b , with the axes pointing right, forward, and upward, respectively. The shearer frame, denoted by the m-frame, is fixed to the shearer with the x m -axis pointing right, the y m -axis pointing forward, and the z m -axis pointing upward. The LiDAR body frame, denoted by the s-frame, originates at the optical center, point O s , with the axes pointing right, forward, and upward, respectively. The n-frame is the navigation frame, chosen as the local level east-north-up (ENU) coordinate. As shown in Figure 4, point i is the point from the hydraulic support that is detected by the LiDAR. The coordinate s i l of point i in the s -frame can be expressed as: where ρ i represents the range between the origin s O of the LiDAR and point i , θ i represents the angle between the vector s i l and s y -axis, and ρ i and θ i are output by the LiDAR. As shown in Figure 4, point i is the point from the hydraulic support that is detected by the LiDAR. The coordinate l s i of point i in the s-frame can be expressed as: where ρ i represents the range between the origin O s of the LiDAR and point i, θ i represents the angle between the vector l s i and y s -axis, and ρ i and θ i are output by the LiDAR.

Velocity Constraint
When the shearer moves along the longwall face and is limited by the AFC, there is no motion normal to the face under ideal conditions. Here, the velocity along the z m -axis is regarded as zero. This velocity constraint can be expressed as a mathematical model as follows: where v m z denotes the shearer velocity along the z m -axis. When the measuring errors are considered, the velocity provided by the velocity constraint is: where w m z denotes the measurement noise of v m z .

Feature Description
The leg is an important part of a hydraulic support, as shown in Figure 5 [31]. The leg in the supported state includes upright and inclined forms. The upright leg is used as the research object in this paper. The green dots indicate the detection results of the leg by the LiDAR which conform to the circle model. The red center of the circle from the i-th hydraulic support is the feature i, defined in this paper. Multiple features can be obtained from a packet of LiDAR data. The absolute position of the feature i is defined, consisting of latitude L, longitude λ, and altitude h,  Figure 2 shows that with the continuous advancement of the longwall face, the hydraulic supports can move intermittently. Therefore, the absolute position of the features can also change intermittently. The process of calculating the absolute position of the features is divided into two steps: initial assignment and position update. Note that the calculation of the feature position is performed when the shearer is located at both ends of the longwall face.

Calculation of the Absolute Position of the Features
• Initial assignment of the absolute position of the features.
Ideally, the hydraulic supports are equally spaced and placed parallel to the AFC at the initial time 0 t . Hence, the feature numbered i is assigned as follows:  Figure 2 shows that with the continuous advancement of the longwall face, the hydraulic supports can move intermittently. Therefore, the absolute position of the features can also change intermittently. The process of calculating the absolute position of the features is divided into two steps: initial assignment and position update. Note that the calculation of the feature position is performed when the shearer is located at both ends of the longwall face.

•
Initial assignment of the absolute position of the features.
Ideally, the hydraulic supports are equally spaced and placed parallel to the AFC at the initial time t 0 . Hence, the feature numbered i is assigned as follows: Remote Sens. 2021, 13, 3708 where m i,0 denotes the absolute position of the feature numbered i at t 0 ; p 0 denotes the shearer position at s 0 is the direction cosine matrix from the s-frame to the b-frame at t 0 ; l s 0 q,0 ≈ [−l x 0 0] T represents the coordinate of the feature numbered q in s-frame at t 0 , which is the closest to the shearer; C b 0 m 0 is the direction cosine matrix from the m-frame to the b-frame at t 0 ; R N and R M denote the transverse and meridian radii of the curvature; and d 1 represents the distance between two adjacent hydraulic supports.
Due to the limited daily mining range of the shearer, C can be considered not to change during the mining process.

•
Update of the absolute position of the features.
According to the introduction of Section 2, the hydraulic supports lagging behind the shearer can move. Thus, when the shearer moves in the next cutting cycle, the absolute position of the features has changed. Therefore, the absolute position of the features needs to be updated. The following takes the update of the position of the feature i as an example for a detailed description.
As shown in Figure 6, t k1 and t k2 represent the two moments before and after feature i moves, respectively. The features j · · · r are the common points of t k1 and t k2 , and they are stationary at these two moments. Set L j:r,k1 = l s j,k1 , · · · , l s r,k1 and L j:r,k2 = l s j,k2 , · · · , l s r,k2 . The relationship between L j:r,k1 and L j:r,k2 is as follows where R 1 and T 1 represent the rotation and translation matrix of the s-frame between t k1 and t k2 , which can be solved using iterative closest point (ICP) [34].
Remote Sens. 2021, 13, 3708 8 of 21 where 1 R and 1 T represent the rotation and translation matrix of the s -frame between 1 k t and 2 k t , which can be solved using iterative closest point (ICP) [34].
The updated position of the feature i can be expressed as: Combining Equations (4) and (8), the position of the feature i can be re-expressed as: where , Considering the measurement error in practice, the relevant variables can be re-expressed as: , , where  0 p denotes the shearer position with error at 0 t ; δ 0 p and  The updated position of the feature i can be expressed as: where C n k1 b k1 represents the attitude matrix at t k1 . Combining Equations (4) and (8), the position of the feature i can be re-expressed as: Remote Sens. 2021, 13, 3708 8 of 20 where D s i,N denotes the displacement of feature i during the N-th advancement in the sframe; ideally, and d 2 is the cutting depth. Equations (6)-(10) represent the calculation update of the absolute position of the features. Considering the measurement error in practice, the relevant variables can be reexpressed as: where p 0 denotes the shearer position with error at t 0 ; δp 0 and the initial position error and misalignment of C Once the SINS and LiDAR are mounted on a shearer, the installation relationship between the b-frame, s-frame, and m-frame is fixed. Ideally, the three frames are reconciled.
The residual mounting angle errors δα and δβ are inevitable in actual situations, but they can be regarded as constant vectors.
Substituting (11) into (9), the feature i position with measurement error can be obtained: where l b 0 i/q,N , and w i is the measurement noise of m i . After expanding Equation (12), we can get: where The near-level longwall face is studied in this paper. Thus, the level attitude angles of the shearer can be regarded as small quantities. We can get the relationship: φ z . According to the method proposed in [22,23], δβ z + φ n 0 z and δβ x can be calibrated in advance. Hence, the position measurement model of feature i can be re-expressed as: • Shearer absolute position calculation. Analogous to the calculation of R 1 and T 1 in Equation (7), the sets L j:r,k = l s j,k , · · · , l s r,k and M j:r,k = m n j,k , · · · , m n r,k are established at t k , and the following relationship is obtained: where m n j,k = (Cˆ− 1)(m j,k − p 0 ), m n r,k = (Cˆ− 1)(m r,k − p 0 ), and R 2 and T 2 represent the rotation and translation matrix between L j:r,k and M j:r,k , respectively. The position measurement model of the shearer provided by LiDAR is: To reduce the influence of R 2 L j:r,k on T 2 , the element of set L j:r,k should be selected closer to LiDAR. After that, M j:r,k becomes the main factor affecting T 2 and p LiDAR . Equation (14) shows that the height information of feature i is affected by many error factors, so only the horizontal information is used when the position measurement is established. At the same time, considering that the error term l x (δα z + φ b 0 z ) is often small, the error-contaminated position of the shearer can be approximated as: where χ(1 : 2) denotes the first two elements of vector χ, and w LiDAR is the noise vector of p LiDAR . It is worth noting that in the actual working environment, a large shift error of the hydraulic support may lead to a decrease in the position accuracy of the corresponding feature, which in turn affects the navigation result. Therefore, in each cutting cycle, when a new feature is detected, it is necessary to determine whether the distance from the adjacent feature meets the requirement to achieve the effect of troubleshooting.

State Space Model
The state equation of the Kalman filter can be expressed as: .
where F I and G are the transition matrix and the distribution matrix of processing the white noise, respectively, as are detailed in [35]; w(t) is the white noise vector; and h x(t) is a 15-dimensional error state vector with reference to the traditional model of the KF, defined as: where δv n is the velocity error; δp is the position error; and ε b and ∇ b are the white noises of the gyros and accelerometers, respectively. In order to perform computer calculations, the continuous Kalman filter needs to be discretized. The discretization form of Equation (18) can be expressed as: where X k and X k−1 denote the discretized state vector at t k and t k−1 , respectively; Φ k/k−1 is the one-step state transition matrix, Φ k/k−1 ≈ exp(F I T S ); Γ k−1 is the system noise driving matrix, Γ k−1 ≈ (I 15 + 0.5F I T S )G; T S is the computation period; I 15 is the 15th-order unit matrix; and W k−1 is the system noise, which satisfies: where Q k is the system noise variance matrix.

Measurement Space Model
The observation equation can be expressed as: where p SI NS is the position of the shearer calculated by the SINS; ν(t) is the measurement ; H LiDAR (1 : 2, :); is the first two rows of matrix H LiDAR ; and H v (3, :) is the third row of matrix H v .
The initial position error δp 0 is the common error term of p SI NS and p LiDAR , which will be offset by the difference between them, so it is not considered in the establishment of the measurement model.
The discretized measurement equation can be expressed as: where Z k is the measurement after discretization; H k is the measurement matrix after discretization; V k is the measurement noise sequence; and R k is the variance matrix of the measurement noise sequence.

Simulation Analysis
To assess the performance of the proposed positioning method, multiple simulations are performed in this section. Affected by mining conditions and the control accuracy of the longwall equipment, the size of the longwall face required will often be different, the AFC may bend, and the hydraulic support may move inaccurately. The above situations are reflected mainly in the simulation trajectory shape and feature distribution. Therefore, we will carry out the simulations from three aspects: different trajectory lengths, different trajectory curvatures, and different feature distribution errors.

Simulations under Different Lengths of the Trajectories
In order to ensure the consistency between the simulation and experiment (mentioned in Section 6), we designed the sensor specifications with reference to the MTi-G-700 and LMS 500 used in the experiment. The sensor specifications are listed in Table 1 [36,37]. The sample rates of the SINS and LiDAR are 100 Hz and 25 Hz, respectively. The initial attitude error of the SINS was φ n T , which is set according to the literature [36].
The residual errors δβ x and δβ z after precalibration were −0.05 • and −0.98 • , respectively.  Figure 7 shows the trajectory shape and feature distribution. In Figure 7a, we designed 50 m, 100 m, and 150 m trajectories for simulation verification. The shearer moved in the direction indicated by the arrows. The operating speed of the shearer was approximately 0.12 m/s. In Figure 7b, the red dots represent the features in the first cutting cycle, which were distributed on the straight line l1 at equal intervals with d 1 = 1.8 m. The distance l x between the straight line l1 and the path of the first cutting cycle was 4 m. The blue dots were the features in the second cutting cycle, which were generated after the red features were advanced. The advancing distance d 2 was 1.2 m. Equation (13) shows that δα z has an impact on the positioning effect. For the three different simulation trajectories, let δα z take δα z 1 : δα z = −1 • , δα z 2 : δα z = 0 • , and δα z 3 : δα z = 1 • in turn, and the corresponding positioning errors are shown in Figure 8.

•
The two sets of features in the adjacent sampling period provided by the LiDAR are used as the input of the ICP algorithm, and the position increment in this period can be obtained according to the output. The process is also known as LiDAR odometry [34]. • When the shearer moves to B, the position obtained by the dead reckoning algorithm based on the SINS and LiDAR odometry can be recorded as In order to verify the effect of the above calibration method, we carried out a series of simulation tests. Assume that the shearer moved quickly on a longwall face. Set φ n 0 z to −1 • ,−0.5 • ,0 • ,0.5 • , and 1 • in turn, and record the residual error δα z after calibration. The results are listed in Table 2. Table 2 shows that the proposed calibration method has a significant effect on reducing the overall error of δα z + φ n 0 z , which is the result we expect.

Simulations under Different Trajectory Curvatures
When the hydraulic support fails to move as required, the AFC is bound to bend, and the maximum bending error is often required to be controlled within 0.1 m [19]. The movement trajectory of the shearer can often reflect the shape of AFC. Therefore, we designed three trajectories, as shown in Figure 9. The maximum bending errors corresponding to Traj 1, Traj 2, and Traj 3 were 0, 0.05, and 0.1 m, respectively. The sensor parameters, φ n 0 , δβ x , and δβ z were the same as those mentioned in Section 5.1. After calibration, δα z was −0.87 • . The feature distribution was the same as shown in Figure 7b. The results are listed in Table 2. Table 2 shows that the proposed calibration method has a significant effect on reducing the overall error of 0 n z z δα φ + , which is the result we expect.

Simulations under Different Trajectory Curvatures
When the hydraulic support fails to move as required, the AFC is bound to bend, and the maximum bending error is often required to be controlled within 0.1 m [19]. The movement trajectory of the shearer can often reflect the shape of AFC. Therefore, we designed three trajectories, as shown in Figure 9. The maximum bending errors corresponding to Traj 1, Traj 2, and Traj 3 were 0, 0.05, and 0.1 m, respectively. The sensor parameters, 0 n φ , δβ x , and δβ z were the same as those mentioned in Section 5.1. After calibration, δα z was −°0.87 . The feature distribution was the same as shown in Figure 7b.

Simulations under Different Feature Distributions
When the hydraulic supports are not moved as required, their corresponding features will not be arranged according to the ideal situation. In order to ensure the continuous advancement of the longwall face, the control accuracy of the hydraulic support is usually required to be within 0.05 m [38]. Considering that the control systems of different hydraulic supports are independent of each other, we can consider that their corresponding features are superimposed random errors at the ideal position. As shown in Figure 11, we generated a set of random sequences smaller than a , and superimposed them on the ideal feature positions to simulate the inaccurate movement of the hydraulic supports. In other words, the actual features are distributed in the square shaded area with the ideal feature as the center and 2 a as the side length. in turn to verify the influence of different feature distributions on the positioning results. The simulation trajectory was selected as Traj 1, shown in Figure 9. The sensor parameters, 0 n φ , δβ x , δβ z , and δα z were the same as those mentioned in Section

Simulations under Different Feature Distributions
When the hydraulic supports are not moved as required, their corresponding features will not be arranged according to the ideal situation. In order to ensure the continuous advancement of the longwall face, the control accuracy of the hydraulic support is usually required to be within 0.05 m [38]. Considering that the control systems of different hydraulic supports are independent of each other, we can consider that their corresponding features are superimposed random errors at the ideal position. As shown in Figure 11, we generated a set of random sequences smaller than a, and superimposed them on the ideal feature positions to simulate the inaccurate movement of the hydraulic supports. In other words, the actual features are distributed in the square shaded area with the ideal feature as the center and 2a as the side length.

Simulations under Different Feature Distributions
When the hydraulic supports are not moved as required, their corresponding features will not be arranged according to the ideal situation. In order to ensure the continuous advancement of the longwall face, the control accuracy of the hydraulic support is usually required to be within 0.05 m [38]. Considering that the control systems of different hydraulic supports are independent of each other, we can consider that their corresponding features are superimposed random errors at the ideal position. As shown in Figure 11, we generated a set of random sequences smaller than a , and superimposed them on the ideal feature positions to simulate the inaccurate movement of the hydraulic supports. In other words, the actual features are distributed in the square shaded area with the ideal feature as the center and 2 a as the side length. in turn to verify the influence of different feature distributions on the positioning results. The simulation trajectory was selected as Traj 1, shown in Figure 9. The sensor parameters, 0 n φ , δβ x , δβ z , and δα z were the same as those mentioned in Section  Figure 9. The sensor parameters, φ n 0 , δβ x , δβ z , and δα z were the same as those mentioned in Section 5.2. The positioning errors are shown in Figure 12. We can see that even though the fluctuation ranges of the errors in the east and north directions will increase with the increase of a, the maximum errors are still maintained in the small ranges. The height error is not related to the change in a. Thus, the positioning method proposed in this paper is still suitable for small hydraulic support displacement errors.
Remote Sens. 2021, 13, 3708 16 of 21 5.2. The positioning errors are shown in Figure 12. We can see that even though the fluctuation ranges of the errors in the east and north directions will increase with the increase of a , the maximum errors are still maintained in the small ranges. The height error is not related to the change in a . Thus, the positioning method proposed in this paper is still suitable for small hydraulic support displacement errors.

Experimental Composition and Design
In order to further evaluate the performance of the positioning method proposed in this paper, an experimental verification was performed. Figure 13 shows the composition of the experimental platform. In Figure 13a, the movement of a mobile carrier was controlled to simulate the operation of the shearer, and a series of cylinders was used to simulate the arrangement of the legs of the hydraulic support. As shown in Figure 13b, the proposed integrated navigation system, including the SINS (Xsens MTi-G-700) and the LiDAR (SICK LMS 500), was installed on the mobile carrier. The parameters of MTi-G-700 and LMS 500 are detailed in Section 5.1. A GPS receiver equipped with an antenna was also placed on the mobile carrier. Through network differential technology, the GPS receiver can output positions with centimeter-level positioning accuracy. The network differential results mainly provide an evaluation basis for the test. In Figure 13c, the cylinder was placed on the base, and the distance between the cylinders was 1.8 m. The advancement of the hydraulic support was simulated by moving the cylinder between two bases with a distance of 1.2 m. (a)

Experimental Composition and Design
In order to further evaluate the performance of the positioning method proposed in this paper, an experimental verification was performed. Figure 13 shows the composition of the experimental platform. In Figure 13a, the movement of a mobile carrier was controlled to simulate the operation of the shearer, and a series of cylinders was used to simulate the arrangement of the legs of the hydraulic support. As shown in Figure 13b, the proposed integrated navigation system, including the SINS (Xsens MTi-G-700) and the LiDAR (SICK LMS 500), was installed on the mobile carrier. The parameters of MTi-G-700 and LMS 500 are detailed in Section 5.1. A GPS receiver equipped with an antenna was also placed on the mobile carrier. Through network differential technology, the GPS receiver can output positions with centimeter-level positioning accuracy. The network differential results mainly provide an evaluation basis for the test. In Figure 13c, the cylinder was placed on the base, and the distance between the cylinders was 1.8 m. The advancement of the hydraulic support was simulated by moving the cylinder between two bases with a distance of 1.2 m.  Figure 12. We can see that even though the fluctuation ranges of the errors in the east and north directions will increase with the increase of a , the maximum errors are still maintained in the small ranges. The height error is not related to the change in a . Thus, the positioning method proposed in this paper is still suitable for small hydraulic support displacement errors.

Experimental Composition and Design
In order to further evaluate the performance of the positioning method proposed in this paper, an experimental verification was performed. Figure 13 shows the composition of the experimental platform. In Figure 13a, the movement of a mobile carrier was controlled to simulate the operation of the shearer, and a series of cylinders was used to simulate the arrangement of the legs of the hydraulic support. As shown in Figure 13b, the proposed integrated navigation system, including the SINS (Xsens MTi-G-700) and the LiDAR (SICK LMS 500), was installed on the mobile carrier. The parameters of MTi-G-700 and LMS 500 are detailed in Section 5.1. A GPS receiver equipped with an antenna was also placed on the mobile carrier. Through network differential technology, the GPS receiver can output positions with centimeter-level positioning accuracy. The network differential results mainly provide an evaluation basis for the test. In Figure 13c, the cylinder was placed on the base, and the distance between the cylinders was 1.8 m. The advancement of the hydraulic support was simulated by moving the cylinder between two bases with a distance of 1.2 m.
(a) As shown in Figure 14, four cutting processes were simulated. Δ and Ο represent the start and end of the trajectory, respectively. The maximum moving distance of the mobile carrier in the longitudinal direction was approximately 44.6 m, and the horizontal advancing distance was about 1.2 m. Its speed was 0.1-0.15 m/s.  As shown in Figure 14, four cutting processes were simulated. Δ and Ο represent the start and end of the trajectory, respectively. The maximum moving distance of the mobile carrier in the longitudinal direction was approximately 44.6 m, and the horizontal advancing distance was about 1.2 m. Its speed was 0.1-0.15 m/s. Figure 14. Trajectory of the mobile carrier during the experiment. Figure 14. Trajectory of the mobile carrier during the experiment.

Experimental Results and Data Analysis
To fully evaluate the advantages of the proposed positioning method, the proposed method was compared with the relative positioning method in the literature [30], and the positioning errors are shown in Figure 15. Figure 15a,b describe the positioning errors without and with calibration of α, respectively. The values of α before and after calibration were [0 0 0] T and [0.56 • 0 − 1.21 • ] T . Regardless of whether measurements were taken before or after the calibration of α, the positioning accuracy of the proposed method is better than the positioning accuracy of the relative method. It can also be seen from Figure 15 that the east and height errors of the relative positioning method without the calibration of α show an increasing-decreasing trend, while those of the proposed method are always stable. The reasons for the above phenomenon are discussed below. method was compared with the relative positioning method in the literature [30], and the positioning errors are shown in Figure 15. Figure 15a,b describe the positioning errors without and with calibration of α , respectively. The values of α before and after calibration were T [0 0 0] and °−°T [0.56 0 1.21 ] . Regardless of whether measurements were taken before or after the calibration of α , the positioning accuracy of the proposed method is better than the positioning accuracy of the relative method. It can also be seen from Figure 15 that the east and height errors of the relative positioning method without the calibration of α show an increasing-decreasing trend, while those of the proposed method are always stable. The reasons for the above phenomenon are discussed below.
The literature [30] shows that when δα z and δα x are present, positioning errors will occur in the east and height, respectively. With the reciprocating movement of the shearer, the east and height errors show an increasing-decreasing trend. The maximum errors can be expressed as δα max z S and δα max x S , where max S represents the maximum longitudinal movement distance. The positioning error caused by δα z for the proposed method is approximately δα z x l , which can be regarded as a constant, while the positioning error has nothing to do with δα x . Usually, max S is much larger than x l , so the influence of δα on the proposed method will be much smaller than on the relative method. The literature [30] shows that when δα z and δα x are present, positioning errors will occur in the east and height, respectively. With the reciprocating movement of the shearer, the east and height errors show an increasing-decreasing trend. The maximum errors can be expressed as δα z S max and δα x S max , where S max represents the maximum longitudinal movement distance. The positioning error caused by δα z for the proposed method is approximately δα z l x , which can be regarded as a constant, while the positioning error has nothing to do with δα x . Usually, S max is much larger than l x , so the influence of δα on the proposed method will be much smaller than on the relative method.
In addition, it can also be seen from Figure 15b that the east error of the relative positioning method with the calibration of α shows a gradual divergence trend. This is mainly affected by the drift of the SINS heading angle. The position errors of the proposed method are generally stable. This result fully proves that the proposed method can reduce the dependence of positioning accuracy on the SINS heading angle accuracy.
The spherical error probable (SEP) is a universal evaluation method of 3D positioning accuracy [22], which can be expressed as: where σ E , σ N , and σ H are the east, north, and height root-mean-square errors, respectively. In Figure 15, "first" through "fourth" represent the four cutting cycle time periods of the shearer, and the remaining time periods correspond to the oblique cutting process of the shearer. The SEP value in each cutting cycle was calculated, and the results are listed in Table 3. The SEP value of the proposed method in each cutting cycle was less than that of the relative method, regardless of whether or not α was calibrated. Compared with the SEP value of the relative positioning method with the calibration of α, the SEP value of the proposed method with the calibration of α was reduced by 37%, 63%, 76%, and 69% from the first to the fourth cutting cycle, respectively. It further shows that the positioning accuracy of the proposed method is better than the relative positioning method. The SEP value of the relative method with the calibration of α decreased by 84%, 69%, 69% and 61%, respectively, as compared to the relative method without the calibration of α, while the SEP value of the proposed method decreased by 46%, 30%, 21%, and 6%, respectively. The positioning accuracy of the relative method before and after calibration was improved more than the proposed method, which further verified that the positioning accuracy of the proposed method was less sensitive to the mounting angle α.

Conclusions
In this paper, we propose a positioning method based on SINS and LiDAR to solve the problem of shearer positioning. According to the motion characteristics of the longwall face, position and velocity measurement models were derived. Based on these measurement models, a 15-dimensional state Kalman filter was established. After a series of simulation analyses, we drew some meaningful conclusions. First, the pitch and azimuth residual mounting angles (δβ x and δβ z ) between the m-frame and b-frame were the main error sources causing the normal and lateral positioning errors, respectively. The maximum errors varied with the trajectory length. The azimuth residual mounting angle δα z between the s-frame and b-frame was one of the main error sources that caused the longitudinal positioning error, which is approximately constant. Therefore, accurate calibration of δβ x , δβ z and δα z is a prerequisite to ensure high positioning accuracy. Second, when compared with the horizontal curvature of the shearer running track, the inaccuracy of the feature distribution had a greater impact on the positioning effect. When the feature had a small shift error, the proposed positioning method could still achieve high accuracy. In addition, we also conducted experimental verification. The results show that the positioning errors of the proposed method are more stable and better than those of the relative positioning method. It is further verified that the proposed method can effectively improve the positioning accuracy of the shearer and has a low dependence on SINS attitude accuracy and the calibration accuracy of the mounting angle α.
The environment around the shearer is extremely complex and harsh. In the longterm mining process, it is inevitable that data loss or even sensor failure will occur. In future work, we will carry out research on sensor fault diagnosis and positioning accuracy maintenance under extreme conditions.