# A LiDAR-Aided Inertial Positioning Approach for a Longwall Shearer in Underground Coal Mining

Authors: Jiangtao Zheng; Sihai Li; Nan Li; Qiangwen Fu; Shiming Liu; Gongmin Yan
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The absolute three-dimensional position of a longwall shearer is fundamental to longwall mining automation. The positioning of the longwall shearer is usually realized by the inertial navigation system (INS) and odometer (OD). However, the position accuracy of this positioning approach gradually decreases over time due to the gyro drift. To further increase the positioning accuracy of the shearer, this paper proposes a positioning approach based on the INS and light detection and ranging (LiDAR). A Kalman filter (KF) model based on the observation provided by detecting hydraulic supports which are part of the longwall face, using the LiDAR, is established. The selection scheme of the point features is studied through a set of simulations. In addition, compared with that of the approach based on the INS and OD, the shearer positioning accuracy obtained using the proposed approach is higher. When the shearer moves along a 350 m track for 6 cutting cycles and lasts about 7.1 h, both east and north position errors can be maintained within 0.2 m and the height error within 0.1 m.


---

## Body

## 1. Introduction

Coal is a basic energy source and valuable raw material in China, and it plays a key role in the national economic domain [1]. However, with the increase in the mining scale and depth, the underground mining risks increase, which necessitates the development of remote and automation solutions that can enhance the personnel safety and increase the mining productivity. A primary solution is to realize underground coal mining based on the longwall face. The main mechanical equipment of the longwall face includes a shearer, certain hydraulic supports, and an armored face conveyor (AFC). The positioning of the longwall shearer is fundamental to the face alignment, creep control, and horizon control, which are the three focus areas in the automated longwall mining process [2]. Therefore, it is critical to examine the longwall shearer positioning.Because of its high autonomy and continuous three-dimensional positioning ability, the inertial navigation system (INS) is widely used to realize the positioning of mining equipment [3–5]. Because the pure inertial positioning error can increase rapidly over time due to the drifting of the inertial devices, integrated navigation approaches involving the INS and auxiliary sensors, which are more stable, are widely applied. The integrated navigation approach based on the INS and global navigation satellite system (GNSS) is a perfect way to achieve the positioning of the mobile platform [6, 7]. However, the GNSS signal is invalid in underground environments. Fan et al. [8] established a shearer positioning method based on the INS and wireless sensor networks (WSNs). In this method, the WSN contained a mobile node and several anchor nodes. The mobile node was installed on the shearer, and the anchor nodes were deployed under the beam of the hydraulic supports. Results showed that the positioning precision based on INS/WSN met the demand of actual working condition. However, the shearer positioning performance is affected by the accuracy of the anchor nodes, which reduces as the hydraulic supports advance. Yang et al. [9] proposed an integrated strategy involving the INS and the zero-velocity update (ZUPT) and motion constraint (MC) techniques. In this case, the INS/ZUPT and INS/MC strategies were adopted when the shearer was stationary and in motion, respectively. According to the MC, the vertical and lateral velocities of the shearer are regarded as zero when the shearer moves normally. However, the lack of a longitudinal velocity constraint in the MC leads to a significant increase in the positioning errors based on the INS/MC when the shearer is maintained in the same orientation for a large period. The shearer positioning method with the INS and odometer (OD) is recognized as a feasible solution [2, 10, 11]. The system uses the dead reckoning (DR) algorithm, with the Euler angles and velocity provided by the INS and OD, respectively, to estimate the shearer position. The positioning accuracy of the DR is related to the OD scale factor error, installation noncoincidence between the INS and the vehicle body, and attitude accuracy of the INS [12]. To improve the shearer positioning accuracy of the DR, the precalibration algorithm can effectively estimate the OD scale factor error and installation noncoincidence [13, 14]. However, the influence of the attitude drift of the DR on the shearer positioning accuracy is inevitable. Wang et al. [15] proposed a closing path optimal estimation model to improve the positioning accuracy calculated by DR. It can improve the positioning accuracy without adding additional sensors, compared with the traditional DR. However, the positioning accuracy of this approach depends on the measurement accuracy of the hydraulic support advancing, and the position accuracy provided by the closing path optimal estimation model can rapidly decrease over time. To further improve the performance and robustness of the shearer positioning system, new technologies must be introduced.The light detection and ranging (LiDAR) technique is being increasingly used in the mining industry because of its excellent environment perception. Li et al. [16] proposed the SLAM algorithm based on LiDAR, which could realize the remote control and automatic positioning of a mine rescue robot. Wu et al. [17] examined the navigation during unmanned driving in a mine and realized the stable output of the tunnel map by using the SLAM algorithm based on LiDAR. W. Charles Kerfoot used the LiDAR approach to examine the issue of mining impacts in coastal regions [18]. Ralston Jonathon C et al. [19] introduced a LiDAR into the longwall face and successfully controlled the creep of the longwall face by installing the LiDAR at the end of the AFC. The surface of the laser equipment is easily condensed by dust, which can affect the equipment performance. The independent dust removal system [20] and air curtain dust removal technology [21] can ensure the cleanliness of the surface of the equipment. With the continuous improvement of dust removal technology, the application of the LiDAR in the mine is bound to be more extensive.Although LiDAR has been used in the longwall face, it is mainly installed on the AFC for creep measurement. DR based on the INS and OD is a universally recognized and widely studied shearer positioning approach. The proposed integrated approach in this paper is to install both INS and LiDAR on the shearer. The proposed approach has two advantages: (1) the observation provided by LiDAR can correct the INS results during the entire mining process, and the positioning accuracy is better than that attained using the DR; (2) the proposed approach still retains the function of measuring the creep of the longwall face. The remaining paper is organized as follows. Section2 describes the working principle of the longwall face. Section 3 presents the mathematical model of the proposed approach. Section 4 describes the error model of the integrated navigation. Section 5 discusses the conducted demonstration tests. Section 6 presents the concluding remarks.


## 2. Working Principle of the Longwall Face

The longwall panel is shown in Figure1. A typical longwall panel [22, 23] has a length and width of 1000–6000 m and 150–350 m, respectively. The width of the roadway is 5–6 m. Before the extraction of a longwall panel, two roadways are excavated with certain coal pillars retained to support the overlying strata. The longwall panel between the two roadways is the coal seam to be extracted by the longwall face, and the goaf is the area that has been extracted. In general, the intersection of the longwall face and the roadway is the area in which the equipment and personnel are concentrated, and the roof support is a necessary measure to improve the safety [24]. Coal is mined by the longwall face, which includes the shearer, certain hydraulic supports, and the AFC, also known as a “three-machine” configuration. The relationship between the three machines and their respective functionalities is as follows: the shearer rides on the AFC and travels back and forth with the AFC as the track. The shearer operation speed is approximately 0.1–0.2 m/s. The AFC, which also assists in the coal transportation, is placed along the coal seam. The hydraulic supports, which are used to support the roof and drive the AFC toward the coal seam, are connected with the AFC and arranged horizontally on the side of the AFC in an equidistant manner. Usually, a longwall face has dozens or even hundreds of sets of hydraulic supports, which are independent of one another. The action of the hydraulic supports can be simply summarized as the static support and dynamic movement. According to the actions of hydraulic supports, we define zone I and zone II. The hydraulic supports in zone I can move to drive the AFC to the coal seam through the jacks, while the hydraulic supports in zone II can static support to provide a safe space for the shearer (as shown in Figure 1). The division rule of zones I and II is that zone I lags behind the rear drum of the shearer about 2-3 hydraulic supports, and zone II covers the remaining hydraulic supports [25].Figure 1
Typical top view of a longwall panel and composition of a longwall face.The trajectory of the shearer can be illustrated considering a sample operation of cutting a triangle coal with an oblique cut at the end [22]. The simplified shearer trajectory is shown in Figure 2(a). To distinguish the two ends of the longwall face, the end at which the shearer is located at the initial moment is termed as end 1, and the other end is termed as end 2. Roadway 1 corresponds to the roadway close to end 1 of the longwall face, and the other roadway is roadway 2. The shearer runs in the sequence A-B-C-D-E-F-G-H-I-J between roadways 1 and 2, and the shearer movement direction is perpendicular to the roadway. A cutting cycle corresponds to the shearer completing one coal cutting between two ends. The movement distance of the shearer along the advancing direction is equidistant, which is called the cutting depth. The cutting depth of each cycle is constant, whose size is about 0.6–1 m. The movements of the shearer from end 1 to end 2 and from end 2 to end 1 are defined as forward and backward movements, respectively. In particular, the shearer position is affected by the movement of the hydraulic supports and the AFC. When the displacement of the hydraulic supports is not equal to the default cutting depth, the trajectory of the shearer changes from straight to curved, as indicated by the dotted line in Figure 2(b). The creep corresponds to the longitudinal displacement of the longwall face, which occurs when the AFC advances. The creep of the longwall face causes the shearer to move toward a roadway, as indicated by the dotted line in Figure 2(c).Figure 2
The simplified shearer trajectory. (a) Ideal shearer trajectory. (b) Shearer trajectory with bend error. (c) Shearer trajectory with creep.
(a)(b)(c)


## 3. Mathematical Model Analysis

### 3.1. System Description

As shown in Figure3, both INS and LiDAR are installed on the shearer body, and the LiDAR faces the hydraulic supports. The related coordinate systems are defined as follows:n-frame: the navigation frame, chosen as the local-level east-north-up (ENU) coordinateb-frame: the INS body frame, implicitly predefined by the calibrated sensitive axes of the inertial sensors, with the origin located at the sensitive center of INS (point OI) and axes pointing right, forward, and upwardm-frame: the shearer body frame, rigidly fixed to the shearer with the xm-axis pointing right, the ym-axis directing forward, and the zm-axis upwardl-frame: the LiDAR body frame, originating at point OL with the xl-axis pointing forward, the yl-axis directing left, and the zl-axis upwardFigure 3
The installation relationship between the INS and LiDAR.In order to facilitate subsequent thesis research, thes-frame is established by rotating the l-frame clockwise by 90°. The s-frame has the same origin as that of the l-frame, and its axes are consistent with the yl-axis, xl-axis, and zl-axis, respectively.It is a common way to express the position of the shearer in the east-north-up coordinate, which is also used to define the position of the shearer in this paper. Since the displacement of the shearer is much smaller than the radius of the earth in one day, then-frame of all positions within the mining range can be considered to be consistent.The pointi is the one that is detected by LiDAR from the surrounding environment. The coordinate Lil of point i in the l-frame can be expressed as(1)Lil=xlylzlT=ρicosθiρisinθi0T,where ρi represents the range between origin OL of the LiDAR and point i and θi represents the angle between the vector Lil and the xl-axis of the l-frame.The vectorLis of point i in the s-frame can be expressed as(2)Lis=ClsLil=−ρisinθiρicosθi0T,where Cls is the transformation matrix from the l-frame to the s-frame. Cls can be expressed as(3)Cls=0−10100001.Note that the LiDAR data in the following analyses are expressed in thes-frame.The measurement error of the LiDAR can be divided into systematic and statistical errors. The LiDAR range and angle involving systematic errors can be expressed as [26, 27](4)ρ˜i=1+δKρi+δb,θ˜i=4a1θi+a22−4a1a3−a22a1,where ρ˜i and θ˜i are the range and angle with errors measured using the LiDAR, δK and δb are the range systematic error model coefficients, and a1−a3 are the angle systematic error model coefficients.The vectorL˜i/LiDARs of point i measured using the LiDAR in the s-frame can be expressed as(5)L˜i/LiDARs=−ρ˜isinθ˜iρ˜icosθ˜i0T≈Lis+δKLis+δb−sinθicosθi0+ρiδθi−cosθi−sinθi0,where δθi=θ˜i−θi.The mounting angles between thes-frame and b-frame are denoted by a vector, α, while the mounting angles between the m-frame and b-frame can be expressed by another vector, β. In general, α and β are inevitable but can be controlled within small enough ranges through precise installation. In order to ensure the positioning accuracy of the shearer, the precalibration processes pertaining to α and β are necessary, whose residual errors after precalibrations can be recorded as δα and δβ, respectively. The direction cosine matrices from the s-frame to the b-frame and from the m-frame to the b-frame can be expressed as(6)C˜sb=I−δα×Csb,C˜mb=I−δβ×Cmb,where Csb and Cmb denote the true direction cosine matrices; C˜sb and C˜mb denote the error-contaminated direction cosine matrices; I is the third-order unit matrix; and δα and δβ are unknown but fixed after precalibrations.


### 3.2. Point Features’ Definition

A hydraulic support consists of four major components, namely, the canopy, gob shield, legs, and base plate, as shown in Figure4 [23]. The legs tilt about 10° to realize the optimal support. The red dots indicate the detection results of the leg by the LiDAR, and the red dots conform to the ellipse model. Because the outer diameter of the leg is usually of the order of decimeters, the lengths of the major and minor axes of the ellipse are similar, which can be approximated as the outer diameter of the leg. The ellipse can be treated as a cycle with a known diameter. The center coordinates of the circle can be calculated by fitting the red dots. The center of the circle is the point feature required in this paper. We can obtain multiple point features from a packet of LiDAR data.Figure 4
Schematic diagram of the hydraulic support and point feature.The principle of using the point features to correct the shearer position is to establish the relationship among the point feature position, shearer position, and relative position between the point feature and the shearer. Subsequently, the Kalman filter (KF) is applied to estimate the position errors of the point features and shearer simultaneously. This algorithm requires the considered point features to be in a static state. It is known from Figure1 that the hydraulic supports in zone I have moved, while the hydraulic supports in zone II are stationary. It means that the point features corresponding to the hydraulic supports in zone II are static, which can be used in the proposed algorithm. The relationship between point feature i, which is selected to participate in the KF, and INS satisfies the following equation:(7)Li/INSs=CbsCnbmin−Pn,where Li/INSs is the coordinate of the point feature i in the s-frame calculated using the INS, Cnb denotes the transformation matrix from the n-frame to the b-frame, min is the position of the point feature i in the n-frame, and Pn is the shearer position calculated using the INS in the n-frame.Furthermore, how to select the point features participated in the KF must be considered. The spatial relationship between the point features and LiDAR is shown in Figure5. The shearer speed is consistent with the ys-axis of the s-frame in the forward movement, while it is opposite in the backward movement. We divide zone II into zone A and zone B with the xs-axis of the s-frame as the boundary, as shown in Figure 5. Zone I and zone II are defined in Section 2. Finally, we think that the selection of point features participated in the KF can be uniquely determined by the following three aspects:The number of point features: suppose thatN adjacent point features are selected to participate in the KF and are numbered as i…i+N−1.The direction relationship between the point features and thexs-axis of the s-frame: when the shearer moves forward as shown in Figure 5(a), the direction relationship between the point features and the xs-axis includes three modes: point features only in zone A, only in zone B, and both in zone A and B. There are also three modes when the shearer moves backward in the same way, as shown in Figure 5(b). Therefore, considering the forward and backward movement, the direction relationship modes can be summarized as listed in Table 1.The positional relationship between the point features and thexs-axis: as shown in Figure 5(a), Li/ys indicates the ordinate of the point feature i in the s-frame. U=Lj/ys denotes the set composed of the ordinates of all point features participated in the KF, and j=i,…,i+N−1. d=minU denotes the minimum value of U. Therefore, as along as d is determined, the positional relationship between the point features and the xs-axis can be fixed.Figure 5
The spatial relationship between the point features and LiDAR. (a) The spatial relationship in the forward movement. (b) The spatial relationship in the backward movement.
(a)(b)Table 1
The modes of direction relationship between the point features and thexs-axis.
Direction relationshipForward movingBackward movingMod 1Only in zone BOnly in zone AMod 2Only in zone BOnly in zone BMod 3Only in zone BBoth in zone A and BMod 4Both in zone A and BOnly in zone AMod 5Both in zone A and BOnly in zone BMod 6Both in zone A and BBoth in zone A and BMod 7Only in zone AOnly in zone AMod 8Only in zone AOnly in zone BMod 9Only in zone ABoth in zone A and B


### 3.3. Calculation of the Creep

When the shearer moves to end 1 of the longwall face two times sequentially, corresponding totk and tk1, the sets of the point features extracted from roadway 1 are denoted as L1 and L2, respectively. The rotation matrix Rr and translation matrix Tr can be solved using the iterative closest point (ICP) algorithm [28, 29], and they satisfy the following relationship:(8)L1=RrL2+Tr.The creep refers to the displacement of the longwall face moving towards the roadway during one advancing process. In order to visually express the creep, we choose them-frame at the initial time t0, which can also be recorded as mt0, to describe it. The creep, denoted by Pcr, can be expressed as(9)Pcr=C˜nmt0CbnCsbTry,where χy denotes the second component of the vector χ and C˜nmt0 denotes the error-contaminated transformation matrix from the n-frame to the mt0-frame.


## 3.1. System Description

As shown in Figure3, both INS and LiDAR are installed on the shearer body, and the LiDAR faces the hydraulic supports. The related coordinate systems are defined as follows:n-frame: the navigation frame, chosen as the local-level east-north-up (ENU) coordinateb-frame: the INS body frame, implicitly predefined by the calibrated sensitive axes of the inertial sensors, with the origin located at the sensitive center of INS (point OI) and axes pointing right, forward, and upwardm-frame: the shearer body frame, rigidly fixed to the shearer with the xm-axis pointing right, the ym-axis directing forward, and the zm-axis upwardl-frame: the LiDAR body frame, originating at point OL with the xl-axis pointing forward, the yl-axis directing left, and the zl-axis upwardFigure 3
The installation relationship between the INS and LiDAR.In order to facilitate subsequent thesis research, thes-frame is established by rotating the l-frame clockwise by 90°. The s-frame has the same origin as that of the l-frame, and its axes are consistent with the yl-axis, xl-axis, and zl-axis, respectively.It is a common way to express the position of the shearer in the east-north-up coordinate, which is also used to define the position of the shearer in this paper. Since the displacement of the shearer is much smaller than the radius of the earth in one day, then-frame of all positions within the mining range can be considered to be consistent.The pointi is the one that is detected by LiDAR from the surrounding environment. The coordinate Lil of point i in the l-frame can be expressed as(1)Lil=xlylzlT=ρicosθiρisinθi0T,where ρi represents the range between origin OL of the LiDAR and point i and θi represents the angle between the vector Lil and the xl-axis of the l-frame.The vectorLis of point i in the s-frame can be expressed as(2)Lis=ClsLil=−ρisinθiρicosθi0T,where Cls is the transformation matrix from the l-frame to the s-frame. Cls can be expressed as(3)Cls=0−10100001.Note that the LiDAR data in the following analyses are expressed in thes-frame.The measurement error of the LiDAR can be divided into systematic and statistical errors. The LiDAR range and angle involving systematic errors can be expressed as [26, 27](4)ρ˜i=1+δKρi+δb,θ˜i=4a1θi+a22−4a1a3−a22a1,where ρ˜i and θ˜i are the range and angle with errors measured using the LiDAR, δK and δb are the range systematic error model coefficients, and a1−a3 are the angle systematic error model coefficients.The vectorL˜i/LiDARs of point i measured using the LiDAR in the s-frame can be expressed as(5)L˜i/LiDARs=−ρ˜isinθ˜iρ˜icosθ˜i0T≈Lis+δKLis+δb−sinθicosθi0+ρiδθi−cosθi−sinθi0,where δθi=θ˜i−θi.The mounting angles between thes-frame and b-frame are denoted by a vector, α, while the mounting angles between the m-frame and b-frame can be expressed by another vector, β. In general, α and β are inevitable but can be controlled within small enough ranges through precise installation. In order to ensure the positioning accuracy of the shearer, the precalibration processes pertaining to α and β are necessary, whose residual errors after precalibrations can be recorded as δα and δβ, respectively. The direction cosine matrices from the s-frame to the b-frame and from the m-frame to the b-frame can be expressed as(6)C˜sb=I−δα×Csb,C˜mb=I−δβ×Cmb,where Csb and Cmb denote the true direction cosine matrices; C˜sb and C˜mb denote the error-contaminated direction cosine matrices; I is the third-order unit matrix; and δα and δβ are unknown but fixed after precalibrations.


## 3.2. Point Features’ Definition

A hydraulic support consists of four major components, namely, the canopy, gob shield, legs, and base plate, as shown in Figure4 [23]. The legs tilt about 10° to realize the optimal support. The red dots indicate the detection results of the leg by the LiDAR, and the red dots conform to the ellipse model. Because the outer diameter of the leg is usually of the order of decimeters, the lengths of the major and minor axes of the ellipse are similar, which can be approximated as the outer diameter of the leg. The ellipse can be treated as a cycle with a known diameter. The center coordinates of the circle can be calculated by fitting the red dots. The center of the circle is the point feature required in this paper. We can obtain multiple point features from a packet of LiDAR data.Figure 4
Schematic diagram of the hydraulic support and point feature.The principle of using the point features to correct the shearer position is to establish the relationship among the point feature position, shearer position, and relative position between the point feature and the shearer. Subsequently, the Kalman filter (KF) is applied to estimate the position errors of the point features and shearer simultaneously. This algorithm requires the considered point features to be in a static state. It is known from Figure1 that the hydraulic supports in zone I have moved, while the hydraulic supports in zone II are stationary. It means that the point features corresponding to the hydraulic supports in zone II are static, which can be used in the proposed algorithm. The relationship between point feature i, which is selected to participate in the KF, and INS satisfies the following equation:(7)Li/INSs=CbsCnbmin−Pn,where Li/INSs is the coordinate of the point feature i in the s-frame calculated using the INS, Cnb denotes the transformation matrix from the n-frame to the b-frame, min is the position of the point feature i in the n-frame, and Pn is the shearer position calculated using the INS in the n-frame.Furthermore, how to select the point features participated in the KF must be considered. The spatial relationship between the point features and LiDAR is shown in Figure5. The shearer speed is consistent with the ys-axis of the s-frame in the forward movement, while it is opposite in the backward movement. We divide zone II into zone A and zone B with the xs-axis of the s-frame as the boundary, as shown in Figure 5. Zone I and zone II are defined in Section 2. Finally, we think that the selection of point features participated in the KF can be uniquely determined by the following three aspects:The number of point features: suppose thatN adjacent point features are selected to participate in the KF and are numbered as i…i+N−1.The direction relationship between the point features and thexs-axis of the s-frame: when the shearer moves forward as shown in Figure 5(a), the direction relationship between the point features and the xs-axis includes three modes: point features only in zone A, only in zone B, and both in zone A and B. There are also three modes when the shearer moves backward in the same way, as shown in Figure 5(b). Therefore, considering the forward and backward movement, the direction relationship modes can be summarized as listed in Table 1.The positional relationship between the point features and thexs-axis: as shown in Figure 5(a), Li/ys indicates the ordinate of the point feature i in the s-frame. U=Lj/ys denotes the set composed of the ordinates of all point features participated in the KF, and j=i,…,i+N−1. d=minU denotes the minimum value of U. Therefore, as along as d is determined, the positional relationship between the point features and the xs-axis can be fixed.Figure 5
The spatial relationship between the point features and LiDAR. (a) The spatial relationship in the forward movement. (b) The spatial relationship in the backward movement.
(a)(b)Table 1
The modes of direction relationship between the point features and thexs-axis.
Direction relationshipForward movingBackward movingMod 1Only in zone BOnly in zone AMod 2Only in zone BOnly in zone BMod 3Only in zone BBoth in zone A and BMod 4Both in zone A and BOnly in zone AMod 5Both in zone A and BOnly in zone BMod 6Both in zone A and BBoth in zone A and BMod 7Only in zone AOnly in zone AMod 8Only in zone AOnly in zone BMod 9Only in zone ABoth in zone A and B


## 3.3. Calculation of the Creep

When the shearer moves to end 1 of the longwall face two times sequentially, corresponding totk and tk1, the sets of the point features extracted from roadway 1 are denoted as L1 and L2, respectively. The rotation matrix Rr and translation matrix Tr can be solved using the iterative closest point (ICP) algorithm [28, 29], and they satisfy the following relationship:(8)L1=RrL2+Tr.The creep refers to the displacement of the longwall face moving towards the roadway during one advancing process. In order to visually express the creep, we choose them-frame at the initial time t0, which can also be recorded as mt0, to describe it. The creep, denoted by Pcr, can be expressed as(9)Pcr=C˜nmt0CbnCsbTry,where χy denotes the second component of the vector χ and C˜nmt0 denotes the error-contaminated transformation matrix from the n-frame to the mt0-frame.


## 4. Integrated Navigation Model

Considering the measuring errors, the navigation system variables are redefined as follows:(10)C˜nb=CnbI+ϕn×,V˜n=Vn+δVn,P˜n=Pn+δPn,m˜in=min+δmin,where Vn denotes the shearer velocity; C˜nb, V˜n, P˜n, and m˜in denote the error-contaminated attitude matrix, velocity, shearer position, and point feature i position, respectively; φn denotes the misalignment angles of C˜nb; φn× is the skew symmetric matrix of φn; and δVn, δPn, and δmin denote the velocity error, shearer position error, and point feature i position error, respectively.Because the positionm˜in of the point feature i can change in different cutting cycles, we need to assign an initial value to it in each cutting cycle. Assume that when the LiDAR detects the point feature i for the first time in the current cutting cycle at tk, the point feature i can be assigned by the following equation:(11)m˜i/kn=P˜kn+C˜bntkC˜sbL˜i/LiDARstk,where m˜i/kn, P˜kn, C˜bntk, and L˜i/LiDARstk denote the error-contaminated position of point feature i, shearer position, attitude matrix, and measurement of the point feature i obtained using the LiDAR, at tk, respectively.


### 4.1. State Space Model

The state vector of the integrated navigation system is defined as(12)X=φnTδVnTδPnTεbT∇bTδminT…δmi+N−1nTT,where εb and ∇b denote the gyro and accelerometer biases in the b-frame, respectively, and δmin…δmi+N−1n denote the position errors of the point features i∼i+N−1 participated in the integrated navigation.The system state equation is(13)X˙=FX+Gw,where F is the state transition matrix and w is the white noise process whose mean and variance are as follows:(14)Ewt=0,EwtwTτ=qδt−τ,where q is the variance intensity matrix.The state transition matrixF can be expressed as(15)F=FINS015×3N03N×1503N×3N,where FINS is the 15×15 state transition matrix based on the classic INS error model [30] and 015×3N, 03N×15, and 03N×3N denote the 15 × 3N, 3N × 15, and 3N × 3N zero matrices, respectively.According to Qin et al. [31], the discretized state equation can be expressed as(16)Xk+1=Φk+1,kXk+Wk.The discretized state transition matrixΦk+1,k and system noise matrix Wk satisfy equations (17) and (18), respectively:(17)Φk+1,k=eTsFtk,where Ts represents the filtering period.(18)EWk=0,EWkWjT=Qkδkj,Qk=∫tktk+1Φk+1,kGkqGkTΦk+1,kTdt.


### 4.2. Measurement Space Model

By substituting equation (10) into (7), the variable L˜i/INSs can be expressed as(19)L˜i/INSs=C˜bsC˜nbm˜in−P˜n≈Lis−CbsCnbmin−Pn×φn−CbsCnbδPn+CbsCnbδmin−Cbsmib−Pb×δα,where min−Pn× is the skew symmetric matrix of the vector min−Pn and mib−Pb=Cnbmin−Pn.After the difference between equations (19) and (5), the measurement equation of the point feature i is(20)Zi/k=L˜i/INSs−L˜i/LiDARs=Hi/kXk+νi/k,where νi/k denotes the measurement noise whose mean and variance are expressed as equation (21) and the measurement matrix Hi/k is expressed as equation (22).(21)Eνi/k=0,Eνi/kνi/jT=Ri/kδkj,where Ri/k is the variance intensity matrix of νi/k.(22)Hi/k=HINSHm,where HINS is the measurement matrix of the 15d INS state vector and Hm is the measurement matrix of the N point features’ state vector.(23)HINS=−CbsCnbmin−Pn×03×3−CbsCnb03×6,Hm=03×3k−1CbsCnb03×3N−k,where k represents the sequence number of the current point feature in all the participated filtering point features.The total measurement matrix ofN point features can be expressed as(24)Zk=Zi/k⋮Zi+N−1/k,Hk=Hi/k⋮Hi+N−1/k.


### 4.3. Theoretical Analysis

The purpose of theoretical analysis is to find the error sources that affect the navigation results and discuss the influence mechanism of the error sources on positioning results. Theoretical analysis is the premise of experimental design.According to equation (11), the shearer positions at adjacent times tk and tk+1 can be expressed as follows:(25)P˜kn=m˜i/kn−C˜bntkC˜sbL˜i/LiDARstk,P˜k+1n=m˜i/k+1n−C˜bntk+1C˜sbL˜i/LiDARstk+1.Due to the low-speed linear running of the shearer and high sampling frequency of the LiDAR, we can obtain equationsC˜bntk≈C˜bntk+1, ρi/k≈ρi/k+1, and δθi/k≈δθi/k+1. The displacement increment of the shearer in a sampling period by substituting equations (5), (6), and (10) into equation (25) can be expressed as(26)ΔP˜kn=P˜k+1n−P˜kn=C˜bnC˜sbL˜i/LiDARstk−L˜i/LiDARstk+1≈ΔPkn−Cbnφb+δα×ΔPkb+CbnδKΔPkb+CbnCsbδbΔθi/kcosθi/ksinθi/k0+CbnCsbδθi/kρi/kΔθi/k−sinθi/kcosθi/k0,where ΔPkb=CnbΔPkn, φb=Cnbφn, and Δθi/k=θi/k+1−θi/k.When the mounting anglesα and β are small, we can obtain the equations Cmb≈Csb≈I and ΔPkb≈ΔPks≈ΔPkm. The displacement increment of the shearer in the m-frame can be expressed as(27)ΔPkm=0ΔSk0T,where ΔSk is the mileage of the shearer in a sampling period.According to the law of sines,Δθi/k can be expressed as(28)Δθi/k=ΔSksinθi/kρi/k=ΔSkLxsρi/k2,where Lxs is the abscissa of the point feature i in the s-frame, which can ideally be regarded as a constant value for all point features.Substituting equations (27) and (28) into (26) yields(29)CnbP˜kn−Pkn=∑CnbΔP˜kn−ΔPkn≈∑ϕzbΔSk0−ϕxbΔSk+δαz0−δαxSk+0δK0Sk+∑δbLxsρi/k2cosθi/kΔSksinθi/kΔSk0+∑δθi/kLxsρi/k−sinθi/kΔSkcosθi/kΔSk0,where Sk represents the mileage of the shearer at tk, and Sk=∑ΔSk.It can be seen from equation (29) thatϕxb and ϕzb can cause the vertical and lateral position errors of the shearer, respectively.δαx and δαz can cause the vertical and lateral position errors of the shearer, respectively, which change linearly with the mileage of the shearer. The maximum mileage Smax of the shearer occurs when it moves to end 2 of the longwall face. The vertical and lateral position errors pertaining to δαx and δαz are also the largest at this position, which can be expressed as δαxSmax and δαzSmax, respectively.δK only causes the longitudinal position error which also changes linearly with the mileage of the shearer. When the shearer moves to end 2 of the longwall face, the maximum longitudinal position error pertaining to δK occurs, which can be expressed as δKSmax.Bothδb and δθi/k can cause the lateral and longitudinal position errors. The lateral and longitudinal position error forms caused by δb and δθi/k are opposite, and they are related to the selection of the point features mentioned in Section 3.2.


## 4.1. State Space Model

The state vector of the integrated navigation system is defined as(12)X=φnTδVnTδPnTεbT∇bTδminT…δmi+N−1nTT,where εb and ∇b denote the gyro and accelerometer biases in the b-frame, respectively, and δmin…δmi+N−1n denote the position errors of the point features i∼i+N−1 participated in the integrated navigation.The system state equation is(13)X˙=FX+Gw,where F is the state transition matrix and w is the white noise process whose mean and variance are as follows:(14)Ewt=0,EwtwTτ=qδt−τ,where q is the variance intensity matrix.The state transition matrixF can be expressed as(15)F=FINS015×3N03N×1503N×3N,where FINS is the 15×15 state transition matrix based on the classic INS error model [30] and 015×3N, 03N×15, and 03N×3N denote the 15 × 3N, 3N × 15, and 3N × 3N zero matrices, respectively.According to Qin et al. [31], the discretized state equation can be expressed as(16)Xk+1=Φk+1,kXk+Wk.The discretized state transition matrixΦk+1,k and system noise matrix Wk satisfy equations (17) and (18), respectively:(17)Φk+1,k=eTsFtk,where Ts represents the filtering period.(18)EWk=0,EWkWjT=Qkδkj,Qk=∫tktk+1Φk+1,kGkqGkTΦk+1,kTdt.


## 4.2. Measurement Space Model

By substituting equation (10) into (7), the variable L˜i/INSs can be expressed as(19)L˜i/INSs=C˜bsC˜nbm˜in−P˜n≈Lis−CbsCnbmin−Pn×φn−CbsCnbδPn+CbsCnbδmin−Cbsmib−Pb×δα,where min−Pn× is the skew symmetric matrix of the vector min−Pn and mib−Pb=Cnbmin−Pn.After the difference between equations (19) and (5), the measurement equation of the point feature i is(20)Zi/k=L˜i/INSs−L˜i/LiDARs=Hi/kXk+νi/k,where νi/k denotes the measurement noise whose mean and variance are expressed as equation (21) and the measurement matrix Hi/k is expressed as equation (22).(21)Eνi/k=0,Eνi/kνi/jT=Ri/kδkj,where Ri/k is the variance intensity matrix of νi/k.(22)Hi/k=HINSHm,where HINS is the measurement matrix of the 15d INS state vector and Hm is the measurement matrix of the N point features’ state vector.(23)HINS=−CbsCnbmin−Pn×03×3−CbsCnb03×6,Hm=03×3k−1CbsCnb03×3N−k,where k represents the sequence number of the current point feature in all the participated filtering point features.The total measurement matrix ofN point features can be expressed as(24)Zk=Zi/k⋮Zi+N−1/k,Hk=Hi/k⋮Hi+N−1/k.


## 4.3. Theoretical Analysis

The purpose of theoretical analysis is to find the error sources that affect the navigation results and discuss the influence mechanism of the error sources on positioning results. Theoretical analysis is the premise of experimental design.According to equation (11), the shearer positions at adjacent times tk and tk+1 can be expressed as follows:(25)P˜kn=m˜i/kn−C˜bntkC˜sbL˜i/LiDARstk,P˜k+1n=m˜i/k+1n−C˜bntk+1C˜sbL˜i/LiDARstk+1.Due to the low-speed linear running of the shearer and high sampling frequency of the LiDAR, we can obtain equationsC˜bntk≈C˜bntk+1, ρi/k≈ρi/k+1, and δθi/k≈δθi/k+1. The displacement increment of the shearer in a sampling period by substituting equations (5), (6), and (10) into equation (25) can be expressed as(26)ΔP˜kn=P˜k+1n−P˜kn=C˜bnC˜sbL˜i/LiDARstk−L˜i/LiDARstk+1≈ΔPkn−Cbnφb+δα×ΔPkb+CbnδKΔPkb+CbnCsbδbΔθi/kcosθi/ksinθi/k0+CbnCsbδθi/kρi/kΔθi/k−sinθi/kcosθi/k0,where ΔPkb=CnbΔPkn, φb=Cnbφn, and Δθi/k=θi/k+1−θi/k.When the mounting anglesα and β are small, we can obtain the equations Cmb≈Csb≈I and ΔPkb≈ΔPks≈ΔPkm. The displacement increment of the shearer in the m-frame can be expressed as(27)ΔPkm=0ΔSk0T,where ΔSk is the mileage of the shearer in a sampling period.According to the law of sines,Δθi/k can be expressed as(28)Δθi/k=ΔSksinθi/kρi/k=ΔSkLxsρi/k2,where Lxs is the abscissa of the point feature i in the s-frame, which can ideally be regarded as a constant value for all point features.Substituting equations (27) and (28) into (26) yields(29)CnbP˜kn−Pkn=∑CnbΔP˜kn−ΔPkn≈∑ϕzbΔSk0−ϕxbΔSk+δαz0−δαxSk+0δK0Sk+∑δbLxsρi/k2cosθi/kΔSksinθi/kΔSk0+∑δθi/kLxsρi/k−sinθi/kΔSkcosθi/kΔSk0,where Sk represents the mileage of the shearer at tk, and Sk=∑ΔSk.It can be seen from equation (29) thatϕxb and ϕzb can cause the vertical and lateral position errors of the shearer, respectively.δαx and δαz can cause the vertical and lateral position errors of the shearer, respectively, which change linearly with the mileage of the shearer. The maximum mileage Smax of the shearer occurs when it moves to end 2 of the longwall face. The vertical and lateral position errors pertaining to δαx and δαz are also the largest at this position, which can be expressed as δαxSmax and δαzSmax, respectively.δK only causes the longitudinal position error which also changes linearly with the mileage of the shearer. When the shearer moves to end 2 of the longwall face, the maximum longitudinal position error pertaining to δK occurs, which can be expressed as δKSmax.Bothδb and δθi/k can cause the lateral and longitudinal position errors. The lateral and longitudinal position error forms caused by δb and δθi/k are opposite, and they are related to the selection of the point features mentioned in Section 3.2.


## 5. Simulation Results and Discussion

To evaluate the performance of the proposed approach, the simulations are performed. First, the simulations are carried out to study the feasible point feature selection scheme, which has been proved to have an impact on the shearer positioning in the analysis of equation (29). Second, a posterior information correction model is proposed to improve the lateral position error, and the simulation results with and without the posterior information correction model are compared and analyzed. Finally, the displacement error of the hydraulic supports and creep of the longwall face are introduced in the simulation, and the proposed approach is compared with the DR based on the INS and OD.It can be seen from equation (29) that the difference of the initial attitude only affects the decomposition of the position error in the ENU coordinate system. Therefore, all simulations in this paper are conducted under the zero initial attitudes for research convenience. This means that φn=φb, and the lateral, longitudinal, and vertical direction of the shearer are consistent with the east, north, and height direction, respectively.


### 5.1. Simulation of Different Selection Methods of Point Features

In order to study the influence of different selection methods of point features on the shearer positioning, the shearer trajectory as shown in Figure6 is designed. The maximum displacement of the shearer is 350 m. The advancing distance of the longwall face is 1 m. The distance between the adjacent point features is 1.5 m. The point features are arranged parallel at a distance of 4 m related to the shearer trajectory. The operating speed of the shearer is about 0.1 m/s.Figure 6
Simulation trajectory of the shearer.The specifications of the INS and LiDAR are listed in Table2. The range and angle systematic error model coefficients are set according to the literature [27]. The sample rates of the INS and LiDAR are 100 Hz and 25 Hz, respectively. The initial attitude error of the INS after the initial alignment [30] is φ0n=−0.23′0.24′−15.60′T. The residual error δα of the mounting angle α between the s-frame and the b-frame after precalibration [32] is −0.60′013.80′T.Table 2
Specifications of the INS and LiDAR.
GyroscopeBias repeatability0.05°/hRandom walk0.01°/hAccelerometerBias repeatability50ugRandom walk10ug/HzLiDARRange systematic error model coefficientδK0.0006Range systematic error model coefficientδb0.035Angle systematic error model coefficienta10.00004Angle systematic error model coefficienta20.9925Angle systematic error model coefficienta30.0531The purpose of this simulation is to verify the influence of the selection of point features participated in the KF, mentioned in Section3.2, on the shearer positioning. First, the estimation effects of the shearer position and attitude under different direction relationship modes between point features and LiDAR are studied, which are shown in Figures 7(a) and 7(b), respectively. The direction relationship modes are set as Table 1. The number of point features participated in the KF is 3. The minimum distance d about the positional relationship between point features and LiDAR is 0 m.Figure 7
Position and attitude errors under different selection methods of point features: shearer (a) position and (b) attitude errors under different direction relationship modes; shearer (c) position and (d) attitude errors under different numbers of point features; shearer (e) position and (f) attitude errors under different minimum distances.
(a)(b)(c)(d)(e)(f)Then, the estimation effects of the shearer position and attitude under different numbers of point features participated in the KF are studied, which are shown in Figures7(c) and 7(d), respectively. The numbers of point features are set as 1–4. The direction relationship mode is set as mod 2 listed in Table 1. The minimum distance d about the positional relationship between point features and LiDAR is 0 m.Finally, the estimation effects of the shearer position and attitude under different minimum distancesd about the positional relationship between point features and LiDAR are studied, which are shown in Figures 7(e) and 7(f), respectively. The minimum distance d is set as 0–7.5 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.As shown in Figure7(a), when the mileage of the shearer increases from 0 m to 350 m for the first time at 0–3540 s, both east and north position errors of the shearer accumulate from zeros to first peak values. When the mileage of the shearer decreases from 350 m to 0 m at 4412–7952 s, the east and north position errors pertaining to mod 2 and mod 7 can decrease from the peak values to near zeros, while the others cannot decrease to near zeros. When the mileage of the shearer increases from 0 m to 350 m for the second time at 8824–12364 s, the peak values of east and north position errors pertaining to mod 2 and mod 7 are similar to the first peak values, while the others cannot be equal to the first peak values. The above phenomenon is more obvious on the north position error curve. That is to say, the east and north position errors corresponding to the modes except mod 2 and mod 7 have shown a certain divergence trend. The height errors of all modes are similar, which change in the zigzag shape. We can see from Figure 7(b) that the direction relationship modes have little effect on the INS attitude, especially the pitch and roll.It can be concluded from Figure7(a) that the direction relationship modes should be selected from mod 2 and mod 7 in order to ensure the long-tern stability of the shearer positioning accuracy.As shown in Figures7(c) and 7(e), increasing both the number of point features and the minimum distances between the selected point features and LiDAR can improve the north positioning accuracy of the shearer. Although the east position error increases first and then decreases with the increase of d, the final improvement is obviously better than that caused by increasing the number of point features. The height position errors have nothing to do with the number of point features and minimum distances. It can be seen from Figures 7(d) and 7(f) that the number of point features and minimum distances have little effect on the INS attitude.It can be concluded from Figures7(c) and 7(e) that increasing both the number of point features and the minimum distances between the selected point features and LiDAR can improve the positioning accuracy, but we should give priority to the latter.


### 5.2. Simulation with and without the Posterior Information Correction Model

Observing Figures7(a), 7(c), and 7(e), it can be found that the north position error can be quickly reduced by adjusting the selection of point features, while the east position error is difficult to achieve. This is because the north position error is mainly caused by δb and δθi/k (the error caused by δK is small), which are related to the selection of point features, and the east position error is not only related to δb and δθi/k but also ϕzb and δαz. In order to minimize the east position error without affecting the north position, we propose a posterior information correction model (PICM). It can be found from equation (29) that δαz has no effect on the longitudinal position error (corresponding to the north position error of the simulation), and the form of the east position error as shown in Figure 7 is similar to that caused by δαz. Therefore, the core idea of the PICM is to add a correction angle η to the precalibrated mounting angle α so that the east position error can be minimized. The calculation process of the correction angle η is as follows.When the shearer moves to end 2 of the longwall face for the first time, its mileage also reaches the maximum, and the shearer positionP˜End2mt0 in the mt0-frame at this time can be expressed as(30)P˜End2mt0=Cnmt0P˜End2n,where P˜End2n represents the shearer position calculated using the INS when it first moves to end 2 of the longwall face and P˜End2mt0=pxpypzT denotes the shearer position calculated using the INS in the mt0-frame.Ideally, the lateral displacement of the shearer in themt0-frame in the first cutting cycle is 0. Therefore, the lateral displacement error of the shearer is px when it first moves to end 2 of the longwall face. The calculation of the correction angle η is as follows:(31)η=00pxSmaxT,where the maximum mileage Smax of the shearer is equal to the difference between the length of the AFC and length of the shearer body.The mounting angleαc between the s-frame and the b-frame adding the correction angle becomes(32)αc=α−η,where α denotes the mounting angle between the s-frame and the b-frame after precalibration mentioned in Section 3.1.Equations (30)–(32) are the mathematical model of the PICM. In order to examine the effect of the PICM, the position and attitude errors of the shearer, with and without the PICM, are presented in Figures 8(a) and 8(b), respectively. The minimum distance d is set as 7.5 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.Figure 8
Position and attitude errors with and without the PICM: (a) shearer position error; (b) shearer attitude error.
(a)(b)We can see from Figure8(a) that the east position accuracy of the shearer is significantly improved with the PICM, which does not affect the north and height position accuracy. It can be seen from Figure 8(b) that the PICM has little effect on the INS attitude.Through the analyses of Figures7 and 8, we can summarize the following point features’ selection scheme:Step 1: fix the number of point features participated in the KF. We can set it as 1.Step 2: determine the minimum distanced about the positional relationship between point features and LiDAR. The principle of determining d is that the larger d, the better, but the hydraulic support corresponding to the selected point feature has not moved in the current cutting cycle.Step 3: determine the direction relationship modes between mod 2 and mod 7 as listed in Table1. The longitudinal position errors py−Smax are calculated separately under mod 2 and mod 7 when the shearer first moves to end 2 of the longwall face, then the direction relationship mode is the one, whose longitudinal position error is smaller between mod 2 and mod 7.Step 4: recalculate the mounting angleαc between the s-frame and the b-frame according to the PICM.


### 5.3. Comparison between the Proposed Approach and DR

The bend and creep of the longwall face are inevitable during actual advancing. Figure9 shows the trajectory with bend and creep to further verify the performance of the proposed INS-/LiDAR-integrated navigation system. The maximum bend error is usually controlled within 0.1 m according to the current control requirements of the longwall face [33]. The creep is approximately 0.16 m as the longwall face advances 1 m [34].The total time of the simulation is about 7.1 h.Figure 9
Simulation trajectory of the shearer with bend and creep: (a) horizontal trajectory with bend and creep; (b) height trajectory with bend.
(a)(b)The specifications of the INS and LiDAR are listed in Table2. The initial attitude error of the INS after the initial alignment is φ0n=−0.18′0.13′−20.9′T. The residual error δα after adding the PICM is −0.6′025.5′T. The residual error δβ of the mounting angle β after precalibration is −0.72′020.3′T [35]. The minimum distance d is set as 9 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.The comparison results between the proposed approach based on the INS and LiDAR, traditional DR based on the INS and OD, and improved DR with the closing path optimal estimation model are shown in Figure10. As shown in Figure 10(a), the position error of the proposed approach is significantly smaller and more stable than that of the traditional DR and improved DR, especially in the east and height directions. As shown in Figure 10(b), the attitude errors of the traditional DR and improved DR exhibit continual divergence, whereas the attitude error of the proposed approach remains stable. It can also be seen from Figure 10(a) that the largest position errors occur at end 2 of each cutting cycle, which are listed in Table 3. E, N, and H mean east, north, and height. Table 3 indicates that the east and north maximum position errors of the proposed approach can be maintained within 0.2 m and the height error within 0.1 m. Compared with the larger and larger maximum position errors of the traditional DR and improved DR, the proposed approach is more stable and accurate.Figure 10
Comparison of the results between the proposed algorithm and DR algorithm: (a) comparison of the shearer position errors; (b) comparison of the shearer attitude errors.
(a)(b)Table 3
Maximum position errors of each cutting cycle.
Cutting cycleTraditional DRImproved DRProposed approachE (m)N (m)H (m)E (m)N (m)H (m)E (m)N (m)H (m)First−0.2100.0870.005−0.2100.0800.005−0.0810.1250.075Second−0.2120.0730.001−0.2120.0810.001−0.0280.1660.082Third−0.5000.107−0.115−0.3320.107−0.103−0.1880.1880.074Fourth−0.4990.106−0.121−0.3300.097−0.110−0.1110.1490.071Fifth−0.6320.177−0.328−0.4530.181−0.2600.1280.1630.054Sixth−0.6360.188−0.336−0.4570.178−0.2680.1080.1330.047The creep of each cutting cycle calculated relative to the previous cycle is listed in Table4. Table 4 indicates that the proposed approach can calculate the creep of the longwall face. Although the errors show a slowly diverging trend as the longwall face advances, the calculated creep can be considered accurate, considering the limited number of times that the longwall face advances each day.Table 4
The creep of the longwall face.
CreepFirstSecondThirdFourthFifthCalculated (m)−0.172−0.223−0.174−0.224−0.176Standard (m)−0.176−0.225−0.175−0.225−0.175Errors (m)0.0040.0020.0010.001−0.001


## 5.1. Simulation of Different Selection Methods of Point Features

In order to study the influence of different selection methods of point features on the shearer positioning, the shearer trajectory as shown in Figure6 is designed. The maximum displacement of the shearer is 350 m. The advancing distance of the longwall face is 1 m. The distance between the adjacent point features is 1.5 m. The point features are arranged parallel at a distance of 4 m related to the shearer trajectory. The operating speed of the shearer is about 0.1 m/s.Figure 6
Simulation trajectory of the shearer.The specifications of the INS and LiDAR are listed in Table2. The range and angle systematic error model coefficients are set according to the literature [27]. The sample rates of the INS and LiDAR are 100 Hz and 25 Hz, respectively. The initial attitude error of the INS after the initial alignment [30] is φ0n=−0.23′0.24′−15.60′T. The residual error δα of the mounting angle α between the s-frame and the b-frame after precalibration [32] is −0.60′013.80′T.Table 2
Specifications of the INS and LiDAR.
GyroscopeBias repeatability0.05°/hRandom walk0.01°/hAccelerometerBias repeatability50ugRandom walk10ug/HzLiDARRange systematic error model coefficientδK0.0006Range systematic error model coefficientδb0.035Angle systematic error model coefficienta10.00004Angle systematic error model coefficienta20.9925Angle systematic error model coefficienta30.0531The purpose of this simulation is to verify the influence of the selection of point features participated in the KF, mentioned in Section3.2, on the shearer positioning. First, the estimation effects of the shearer position and attitude under different direction relationship modes between point features and LiDAR are studied, which are shown in Figures 7(a) and 7(b), respectively. The direction relationship modes are set as Table 1. The number of point features participated in the KF is 3. The minimum distance d about the positional relationship between point features and LiDAR is 0 m.Figure 7
Position and attitude errors under different selection methods of point features: shearer (a) position and (b) attitude errors under different direction relationship modes; shearer (c) position and (d) attitude errors under different numbers of point features; shearer (e) position and (f) attitude errors under different minimum distances.
(a)(b)(c)(d)(e)(f)Then, the estimation effects of the shearer position and attitude under different numbers of point features participated in the KF are studied, which are shown in Figures7(c) and 7(d), respectively. The numbers of point features are set as 1–4. The direction relationship mode is set as mod 2 listed in Table 1. The minimum distance d about the positional relationship between point features and LiDAR is 0 m.Finally, the estimation effects of the shearer position and attitude under different minimum distancesd about the positional relationship between point features and LiDAR are studied, which are shown in Figures 7(e) and 7(f), respectively. The minimum distance d is set as 0–7.5 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.As shown in Figure7(a), when the mileage of the shearer increases from 0 m to 350 m for the first time at 0–3540 s, both east and north position errors of the shearer accumulate from zeros to first peak values. When the mileage of the shearer decreases from 350 m to 0 m at 4412–7952 s, the east and north position errors pertaining to mod 2 and mod 7 can decrease from the peak values to near zeros, while the others cannot decrease to near zeros. When the mileage of the shearer increases from 0 m to 350 m for the second time at 8824–12364 s, the peak values of east and north position errors pertaining to mod 2 and mod 7 are similar to the first peak values, while the others cannot be equal to the first peak values. The above phenomenon is more obvious on the north position error curve. That is to say, the east and north position errors corresponding to the modes except mod 2 and mod 7 have shown a certain divergence trend. The height errors of all modes are similar, which change in the zigzag shape. We can see from Figure 7(b) that the direction relationship modes have little effect on the INS attitude, especially the pitch and roll.It can be concluded from Figure7(a) that the direction relationship modes should be selected from mod 2 and mod 7 in order to ensure the long-tern stability of the shearer positioning accuracy.As shown in Figures7(c) and 7(e), increasing both the number of point features and the minimum distances between the selected point features and LiDAR can improve the north positioning accuracy of the shearer. Although the east position error increases first and then decreases with the increase of d, the final improvement is obviously better than that caused by increasing the number of point features. The height position errors have nothing to do with the number of point features and minimum distances. It can be seen from Figures 7(d) and 7(f) that the number of point features and minimum distances have little effect on the INS attitude.It can be concluded from Figures7(c) and 7(e) that increasing both the number of point features and the minimum distances between the selected point features and LiDAR can improve the positioning accuracy, but we should give priority to the latter.


## 5.2. Simulation with and without the Posterior Information Correction Model

Observing Figures7(a), 7(c), and 7(e), it can be found that the north position error can be quickly reduced by adjusting the selection of point features, while the east position error is difficult to achieve. This is because the north position error is mainly caused by δb and δθi/k (the error caused by δK is small), which are related to the selection of point features, and the east position error is not only related to δb and δθi/k but also ϕzb and δαz. In order to minimize the east position error without affecting the north position, we propose a posterior information correction model (PICM). It can be found from equation (29) that δαz has no effect on the longitudinal position error (corresponding to the north position error of the simulation), and the form of the east position error as shown in Figure 7 is similar to that caused by δαz. Therefore, the core idea of the PICM is to add a correction angle η to the precalibrated mounting angle α so that the east position error can be minimized. The calculation process of the correction angle η is as follows.When the shearer moves to end 2 of the longwall face for the first time, its mileage also reaches the maximum, and the shearer positionP˜End2mt0 in the mt0-frame at this time can be expressed as(30)P˜End2mt0=Cnmt0P˜End2n,where P˜End2n represents the shearer position calculated using the INS when it first moves to end 2 of the longwall face and P˜End2mt0=pxpypzT denotes the shearer position calculated using the INS in the mt0-frame.Ideally, the lateral displacement of the shearer in themt0-frame in the first cutting cycle is 0. Therefore, the lateral displacement error of the shearer is px when it first moves to end 2 of the longwall face. The calculation of the correction angle η is as follows:(31)η=00pxSmaxT,where the maximum mileage Smax of the shearer is equal to the difference between the length of the AFC and length of the shearer body.The mounting angleαc between the s-frame and the b-frame adding the correction angle becomes(32)αc=α−η,where α denotes the mounting angle between the s-frame and the b-frame after precalibration mentioned in Section 3.1.Equations (30)–(32) are the mathematical model of the PICM. In order to examine the effect of the PICM, the position and attitude errors of the shearer, with and without the PICM, are presented in Figures 8(a) and 8(b), respectively. The minimum distance d is set as 7.5 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.Figure 8
Position and attitude errors with and without the PICM: (a) shearer position error; (b) shearer attitude error.
(a)(b)We can see from Figure8(a) that the east position accuracy of the shearer is significantly improved with the PICM, which does not affect the north and height position accuracy. It can be seen from Figure 8(b) that the PICM has little effect on the INS attitude.Through the analyses of Figures7 and 8, we can summarize the following point features’ selection scheme:Step 1: fix the number of point features participated in the KF. We can set it as 1.Step 2: determine the minimum distanced about the positional relationship between point features and LiDAR. The principle of determining d is that the larger d, the better, but the hydraulic support corresponding to the selected point feature has not moved in the current cutting cycle.Step 3: determine the direction relationship modes between mod 2 and mod 7 as listed in Table1. The longitudinal position errors py−Smax are calculated separately under mod 2 and mod 7 when the shearer first moves to end 2 of the longwall face, then the direction relationship mode is the one, whose longitudinal position error is smaller between mod 2 and mod 7.Step 4: recalculate the mounting angleαc between the s-frame and the b-frame according to the PICM.


## 5.3. Comparison between the Proposed Approach and DR

The bend and creep of the longwall face are inevitable during actual advancing. Figure9 shows the trajectory with bend and creep to further verify the performance of the proposed INS-/LiDAR-integrated navigation system. The maximum bend error is usually controlled within 0.1 m according to the current control requirements of the longwall face [33]. The creep is approximately 0.16 m as the longwall face advances 1 m [34].The total time of the simulation is about 7.1 h.Figure 9
Simulation trajectory of the shearer with bend and creep: (a) horizontal trajectory with bend and creep; (b) height trajectory with bend.
(a)(b)The specifications of the INS and LiDAR are listed in Table2. The initial attitude error of the INS after the initial alignment is φ0n=−0.18′0.13′−20.9′T. The residual error δα after adding the PICM is −0.6′025.5′T. The residual error δβ of the mounting angle β after precalibration is −0.72′020.3′T [35]. The minimum distance d is set as 9 m. The direction relationship mode is set as mod 2 listed in Table 1. The number of point features is 1.The comparison results between the proposed approach based on the INS and LiDAR, traditional DR based on the INS and OD, and improved DR with the closing path optimal estimation model are shown in Figure10. As shown in Figure 10(a), the position error of the proposed approach is significantly smaller and more stable than that of the traditional DR and improved DR, especially in the east and height directions. As shown in Figure 10(b), the attitude errors of the traditional DR and improved DR exhibit continual divergence, whereas the attitude error of the proposed approach remains stable. It can also be seen from Figure 10(a) that the largest position errors occur at end 2 of each cutting cycle, which are listed in Table 3. E, N, and H mean east, north, and height. Table 3 indicates that the east and north maximum position errors of the proposed approach can be maintained within 0.2 m and the height error within 0.1 m. Compared with the larger and larger maximum position errors of the traditional DR and improved DR, the proposed approach is more stable and accurate.Figure 10
Comparison of the results between the proposed algorithm and DR algorithm: (a) comparison of the shearer position errors; (b) comparison of the shearer attitude errors.
(a)(b)Table 3
Maximum position errors of each cutting cycle.
Cutting cycleTraditional DRImproved DRProposed approachE (m)N (m)H (m)E (m)N (m)H (m)E (m)N (m)H (m)First−0.2100.0870.005−0.2100.0800.005−0.0810.1250.075Second−0.2120.0730.001−0.2120.0810.001−0.0280.1660.082Third−0.5000.107−0.115−0.3320.107−0.103−0.1880.1880.074Fourth−0.4990.106−0.121−0.3300.097−0.110−0.1110.1490.071Fifth−0.6320.177−0.328−0.4530.181−0.2600.1280.1630.054Sixth−0.6360.188−0.336−0.4570.178−0.2680.1080.1330.047The creep of each cutting cycle calculated relative to the previous cycle is listed in Table4. Table 4 indicates that the proposed approach can calculate the creep of the longwall face. Although the errors show a slowly diverging trend as the longwall face advances, the calculated creep can be considered accurate, considering the limited number of times that the longwall face advances each day.Table 4
The creep of the longwall face.
CreepFirstSecondThirdFourthFifthCalculated (m)−0.172−0.223−0.174−0.224−0.176Standard (m)−0.176−0.225−0.175−0.225−0.175Errors (m)0.0040.0020.0010.001−0.001


## 6. Conclusions

This paper proposes a shearer positioning approach based on the INS and LiDAR. The shearer position accuracy is determined by the INS attitude, mounting angles between thes-frame and the b-frame, and ranging and angle accuracy of the LiDAR. By adjusting the selection of the point features, the positioning accuracy of the shearer can be improved, especially the north position accuracy. When the posterior information model is executed, the accuracy of the east position can be further improved. The test results prove that the proposed approach can meet the current demand for shearer positioning. In addition, the proposed approach can calculate the creep of the longwall face by using the ICP algorithm when the shearer moves to end 1 of the longwall face. The accurate shearer position and longwall face creep can provide reference for the automated mining of the longwall face.



---

*Source: 6616090-2021-02-25.xml*