A coal mine roadway is a longitudinally limited space with curves and branches, low illumination and high humidity, a large amount of dust, and an unstructured terrain environment. Traditional ICP algorithms have the defects of slow convergence speed and it is easy to fall into local optimums. While the NDT algorithm in the NDT + ICP algorithm has high registration efficiency, poor stability, and low registration accuracy, which are not suitable for point clouds with noise and a large amount of data. By calculating the FPFH value, the detailed description of the point cloud will be greatly increased to increase the robustness and accuracy Therefore, a feature registration method based on the FPFH + ICP algorithm is proposed to reduce the modeling error of excavation roadways and meet the requirements of intelligent excavation. First, outliers caused by dust are treated by the Euclidean clustering point cloud segmentation method, and then the calculation of the normal vector in the FPFH feature descriptor is optimized based on extracting key points from the roadway structure. The surface normal vector of each key point and its neighborhood point is estimated according to the measured point and its neighborhood point. The initial coordinate transformation matrix of a point cloud of an excavated roadway is obtained by the SAC-IA algorithm and transferred to the ICP algorithm. Finally, KD-tree is introduced into the ICP algorithm to accelerate the search speed of corresponding point pairs, and the Gauss–Newton method is used to optimize the solution of the nonlinear objective function of the algorithm to complete accurate registration of point clouds in an excavation roadway.

1. Introduction

China’s total coal consumption accounts for more than 60% of the total energy consumption for a long time, and coal will still be the pillar of energy in China for a long time. To ensure coal mine safety production, reduce costs, and increase efficiency, China vigorously promotes the construction of “mechanized, automated, and intelligent” intelligent mines and sets the development goal of an intelligent coal mine system with intelligent perception, intelligent decision-making, and automatic execution [1]. Excavation is the foundation of realizing intelligent coal mines, and the establishment of an excavation roadway model is the premise of realizing intelligent unmanned excavation [2].

It is impossible to build a complete model of the excavation roadway by measurement and other methods due to environmental characteristics of the excavation roadway such as low illumination, high humidity, and high density of particulate dust [35], reflectivity changes caused by the blockage of large particles and dust in the air, the properties of the surface of the object, the material and transparency of the surface, etc., and characteristics of excavation roadways, such as harsh working conditions, narrow spaces, and complex environments, which cannot be characterized [6, 7]. To construct the environment model around the excavation roadway, a suitable coordinate system is needed to merge the three-dimensional point clouds obtained from various perspectives into a unified coordinate system to form a complete data point cloud. A complete roadway model is constructed by point cloud registration technology [8, 9] to meet the intelligent tunneling requirements of an excavation roadway.

Point cloud registration techniques are divided into scanning-based matching methods and map feature-based matching methods at present. The purpose of scan-based matching is to further calculate the position and attitude of objects in the point cloud map by comparing the similarity between the current scanned point cloud and the reference point cloud. Common algorithms include the ICP (iterative closest points) algorithm [10] and the NDT (normal distributions transform) algorithm [11]. Tamaki et al. combined CUDA and ICP to speed up the matching speed, but it is still difficult to meet the requirements of fast registration of actual point clouds [12]. Takeuchi et al. further improved the NDT algorithm and applied it to the matching of three-dimensional point clouds to achieve a breakthrough from two-dimensional to three-dimensional [13]. Zhang et al. improved the NDT algorithm based on the SURF algorithm and applied it to ground laser scanning. The algorithm has a good matching effect on point clouds with different resolutions [14]. After that, experts began to combine the advantages of ICP and NDT to improve the matching accuracy and efficiency. Qingshan et al. combined ICP with the NDT algorithm and conducted experiments on campus to verify its reliability and accuracy in 2019 [15]; Guiyang et al. fused NDT and ICP algorithms to verify that the NDT-ICP algorithm can improve the speed and accuracy of point cloud registration [16]; Shi et al., simplified the NDT-ICP three-step method into a two-step method to improve the efficiency of the algorithm and reduce computational complexity [17]. Matching based on map features mainly realizes point cloud registration according to the feature information on the map. Chua et al. calculated the point signature descriptors of all points in the point cloud to find the corresponding point pairs [18]. Inspired by this, scholars have carried out a large amount of research on feature descriptors. Johnson et al. put forward the spin image SI descriptor, which does not change with the motion of the point cloud rigid body but has a large amount of computation and a slow computation speed [19, 20]. Rusu et al., proposed PFH and FPFH point cloud feature extraction algorithms [21, 22]. The outstanding advantages of the FPFH feature with only 33 dimensions are high computational efficiency and fewer storage resources. The sparse surface of the point cloud has excellent identification ability [23]. Tombari et al. proposed the signature of the histogram of orientation (SHOT) descriptor, which has good description performance but a complicated calculation process [24, 25]. Guo et al. summarized the existing local descriptors and analyzed the performance of each local descriptor. Experiments show that the FPFH descriptor has reasonable descriptive and computational efficiency and provides a good balance between feature-matching accuracy and computational efficiency [26].

The traditional ICP method requires a high initial position of a point cloud, so the nearest point changes when iteratively optimize the rotation matrix and translation matrix, which makes it easy to fall into a local optimum when matching point clouds in an excavation roadway with a large amount of data, destroying the roadway model structure. In the NDT + ICP algorithm, too much NDT as a parameter of the rough registration lattice leads to low accuracy of the roadway model. The closed environment and dust in the roadway make the amount of data gathered in the roadway large and there are outliers. If the lattice parameters are too small, the memory is too high, and it is difficult to remove outliers, which leads to low registration accuracy and poor stability. Based on this, the point cloud feature matching method of FPFH + ICP is proposed and optimized according to the internal environmental characteristics of the excavation roadway. By calculating the FPFH factor and adding feature constraints, a more accurate matching result is obtained without the initial value. Combined with the high-precision characteristics of the ICP algorithm, the degradation phenomenon of point clouds in the sparse roadway environment is prevented, and the problem of low matching accuracy due to the particularity of the roadway environment is solved. It provides technical support for the construction of an environmental map and underground positioning of the following tunneling roadway and is an important step toward realizing underground intelligence.

2. Constraint Conditions of Environmental Modeling of Excavation Roadway

Scientific and technical research of “digital mine” is a key research direction in the field of geology and mineral resources at present, and the dangerous area in the underground roadway is heading face. Therefore, the construction and spatial analysis of the network model of the excavation roadway is an important research topic in digital mining [27]. The model contains all the information about virtual and digital heading faces in physical space, mainly including the current position display of the road header, spatial position distribution display of obstacles, head-on shape scanning of the heading face, etc. Real-time communication and interactive feedback are carried out in physical space and virtual space, as shown in Figure 1. Based on point cloud registration technology, this paper constructs the environment model of the excavation roadway (excavation roadway) and builds the interactive foundation of virtual and real space to help constrain the movement space of the road header and control the forming quality of the excavation roadway.

2.1. Environmental Perception Mode of Excavation Roadway

The high temperature and humidity of the underground environment in coal mines are affected by the geothermal effect, heat dissipation of the human body and electromechanical equipment, and the air temperature is sometimes as high as 35°C and close to 100% relative humidity [28]. There is a lot of dust and gas in the underground environment, among which the dust mainly comes from all kinds of coal dust and rock dust produced by mining coal and rock by roadheader, processing, transportation, and loading. The total dust mass concentration can reach [29]. There is no natural light source underground, so underground lighting is generally illuminated by moving lamps and fixed illumination, which leads to low illumination of the roadway. The illumination of the heading face and transportation roadway is about 5 Lx. The visual camera works poorly in the scene of high humidity, high dust concentration, and low illumination, as shown in Figure 2 [30], so it is not suitable for modeling excavation roadways.

Based on the principle of millimeter-wave radar frequency modulated continuous wave, millimeter-wave point cloud obtains point cloud data by measuring the distance between surrounding points and radar echo intensity. Point cloud data contains a variety of target information, but the point cloud obtained by millimeter-wave radar is sparse, which brings new challenges to feature extraction, matching, and recognition of roadway structures. In addition, millimeter-wave radar for long-distance and dynamic target detection ability is extremely limited and cannot get real-time roadway environmental change information and road header position information in the construction of the excavation roadway environment model. Therefore, lidar with high precision, strong directivity, and easier real-time modeling and positioning is more suitable for building the environmental model of an excavation roadway.

2.2. Analysis of Environmental Modeling Characteristics of Excavation Roadway

Coal mine roadways can be regarded as a multifork pipe model with a rough surface, a closed structure, and a longdistance. There are many forks in the roadway, and the inner surface of the roadway is different with obstacles such as gravel, materials, electromechanical equipment, etc. [31, 32]. The internal environment of the roadway is shown in Figure 3. Its interior is narrow, damp, and dark. The ground is cut from rock along the trend of underground mountains, so it is rugged. Steps, steep slopes, tracks, ditches, puddles, and muddy land can be seen everywhere. Because of the particularity of the interior structure of the roadway, it is difficult to extract the features of the roadway based on the scanning matching method when modeling the long-distance and equal-width roadways, which leads to the lack of precision in the established roadway model and cannot well display the information of steps, steep slopes, and obstacles in the roadway [33, 34]. Therefore, it is difficult to select a suitable point cloud registration method to realize multiframe point cloud registration according to the unique characteristics of the roadway.

2.3. Analysis of Point Cloud Data Characteristics of Excavation Roadway

Roadway excavation is the first work of coal mining, which plays a decisive role in later mining. During tunneling, the road header moves forward through the walking mechanism and breaks and excavates rocks through the working mechanism to dig a passable roadway on the straight ground of the mining area. A large amount of dust will be caused by this process. There are always some outliers that deviate obviously from the main point cloud in the three-dimensional point cloud due to the interference of the large particle dust in the air, the reflectivity change caused by the surface property, the surface material, and transparency, etc. [35]. The existence of these points will cause a great deviation in the point cloud registration. Therefore, it is necessary to identify and filter outliers. In addition, the original roadway point cloud data is huge. To improve the speed of the subsequent excavation roadway model building and reduce resource occupation, the original roadway point cloud needs to be simplified and compressed.

The outlier cloud has the following mathematical characteristics:where and are the mean and variance of the average distance from the outlier to all point clouds in the neighborhood of u. The size of T depends on the number of point clouds in the U.

Therefore, for the special environment and constraint conditions of the excavation roadway, it is necessary to choose the appropriate point cloud preprocessing and point cloud registration methods.

3. Point Cloud Registration Method and Analysis of Excavation Roadway Modeling

3.1. ICP Algorithm

The ICP method uses a four-element operation model to find the nearest point in the point cloud matching relationship and iterates through the objective function of the Euclidean distance. Its principle is to find the nearest point in the target point cloud and the source point cloud to be matched according to certain constraints, and then calculate the optimal matching parameters and , so that the error function minimum error function formula is as shown in the following formula :where n is the number of nearest point pairs, is a point in the target point cloud , and is the nearest point corresponding to in the source point cloud , is the rotation matrix, and is the translation vector.

3.2. NDT Algorithm

The principle of the NDT algorithm is first gridding point clouds, then converting the data in each three-dimensional grid into a continuous differentiable probability density distribution function, and finally completing the registration of point clouds by the Hessian matrix method. Core formulas are shown in the following equations :where represents the mean vector, represents the covariance matrix, n is the total number of points in the cube, is a point in the cube, and represents a transformation matrix.

3.3. NDT + ICP Algorithm

The ICP algorithm performs point-to-point matching by finding the corresponding point pairs between two-point clouds; the NDT algorithm calculates the probability model of the point cloud for matching. For two different point cloud matching algorithms, the NDT + ICP algorithm combines the high efficiency of the NDT algorithm with the high precision of the ICP algorithm. It is mainly divided into two key parts: NDT coarse registration and ICP fine registration, as follows:(i)The point cloud to be registered is greatly narrowed to the source point cloud by the NDT algorithm, and the transformation matrix between two-point clouds is output as shown in the following equation :(ii)The point cloud transformation matrix obtained from the solution is used as the initial matrix of the ICP algorithm to ensure the consistency of registration.(iii)KD-tree is used to speed up the search for corresponding points when the ICP algorithm is used for precise matching.(iv)Gauss–Newton method is applied to solve the objective function of formula (2) to improve the computational efficiency of ICP matching.

3.4. Experimental Comparison of Algorithms

In order to verify the effectiveness of the above point cloud registration algorithm, this paper selects 3 groups of roadway point clouds. ① (roadway1-a.pcd roadway1-b.pcd), ② (roadway2-a.pcd roadway2-b.pcd), and ③ (roadway3-a.pcd roadway3-b.pcd). Among them, roadway1-a.pcd, roadway2-a.pcd, and roadway3-a.pcd are the reference point clouds, and the number of point clouds is 28370, 28484, and 27218; roadway1-b.pcd, roadway2-b.pcd, and roadway3-b.pcd are point clouds to be registered, and the number of point clouds is 28344, 28437, and 27216. The master computer operating system is Windows 10, and the running memory is an 8 G processor, Intel (R) Core (TM) i5-7400 CPU @ 3.00 GHz. The test software is Visual Studio 2017. The experimental renderings and data diagrams are shown in Table 1 of Figure 4.

The experimental results show that the convergence speed of the traditional ICP algorithm is slower, the time is the longest, and the initial position of the point cloud data set is higher; otherwise, the matching results will destroy the structural characteristics of the environment and lead to the deformation of the model. The NDT algorithm has low registration accuracy in a noisy and large data point cloud and is not suitable for a closed roadway with a dust environment. It is generally performed as a rough registration algorithm. Although the NDT + ICP algorithm improves the operation efficiency, it does not greatly improve the modeling accuracy. The FPFH point feature histogram descriptor has excellent discrimination ability on sparse point cloud surfaces, which will greatly improve the modeling accuracy of excavation roadways with a large amount of point cloud data and an unstructured terrain environment. Therefore, the FPFH + ICP registration method was selected for roadway environment modeling in this paper.

4. Modeling of Excavation Roadway Based on FPFH + ICP

4.1. Pretreatment Method of the Point Cloud in Excavation Roadway

In this paper, the Euclidean point cloud clustering segmentation method is used to remove outliers in point clouds and a Euclidean distance is defined to describe the geometric distance between two different clusters. As shown in the following equation , m and n are points that the cluster contains that can be calculated as follows:

An adaptive calculation is implemented for the threshold, and the threshold R is defined as in the following equation :, , and are coordinates of the point or center point of the point cloud to be processed. a and b are adjustment parameters whose values usually vary with the angular accuracy of the lidar. The higher the angle accuracy, the smaller their values are, and vice versa. The steps are as follows:(1)Define a point cloud set M, and set a Euclidean distance threshold R by formula (7).(2)Select another point cloud set n and find a point in the point cloud of the set N as a reference point, find K points near it through the KD-tree.(3)Calculate the Euclidean distance between each point and a certain point in the set m according to formula (6), and put all points whose distance is less than the threshold value into the set m.(4)Find another point in the set n as a new reference point and repeat the above operation until no new point in the point cloud set n is added to the point cloud set m.(5)Select another point cloud set to find the point and judge whether the belongs to the divided set; if so, abandon the point to find the next point. Otherwise, continue the above steps.

4.2. Point Cloud Registration Method for Excavation Roadway

After the point cloud preprocessing, the reference point cloud and the point cloud to be registered is obtained. Start the point’s cloud registration. In this paper, a point cloud registration method based on the FPFH + ICP algorithm is proposed. The specific flow is shown in Figure 5. First, the key points in the point cluster are extracted by the ISS algorithm to optimize the feature parameters in FPFH. Through the SAC-IA algorithm, the calculated FPFH features are used to complete the rough registration process of the point cloud; the rigid body transformation matrix obtained from rough registration is used as the initial matrix of ICP fine registration; and the KD-tree is introduced into the ICP algorithm to accelerate the search of point pairs while ensuring the consistency of registration. Finally, the Gauss–Newton method is used to optimize the objective function of the ICP algorithm iteratively to solve the optimal coordinate transformation parameters between point clouds to complete the accurate registration of point clouds.

4.2.1. Extraction of Key Points of the Point Cloud in Excavation Roadway

The key is the stability of the point cloud Compared with the original point cloud data, the number of differentiated point sets is much smaller. They are combined with local feature descriptors to form key point descriptors, which are often used to form a compact representation of point clouds without losing representativeness and descriptiveness, to accelerate the subsequent processing speed. The steps for obtaining key points by the ISS (intrinsic shape signature) algorithm are as follows:

Query all points within the radius of each point in the data obtained from the point cloud once, and calculate the weight , as shown in the following equation:(1)Calculate the variance matrix according to the weights as shown in the following equation :(2)Calculate the eigenvalue of the covariance matrix, and arrange the eigenvalues in order from large to small.(3)Set thresholds and , keep the points that meet / and / as key points, and repeat the above steps until all points are completed.

4.2.2. Calculation of FPFH Characteristics of the Point Cloud in Excavation Roadway

The Fast Point Feature Histogram (FPFH) is improved from PFH. The flow shown in Figure 6 is as follows:

For each query point , calculate the relative relationship between the point to be sought and all its domain points, which is denoted as S ();

Redefine the K field of each point, then estimate the FPFH feature as F () from the calculated SPFH feature, we getwhere is the weighted value of the SPFH feature of the i-th domain point, and represents the distance value between the point to be found and its i-th domain point for evaluating the relationship between point pairs (, ).

4.2.3. Rough Registration of Point Cloud in Excavation Roadway

To avoid the point cloud falling into the local optimum in the registration process, the initial transformation of the point cloud is carried out first. In the initial registration stage, a sample consistency initial alignment (SAC-IA) is proposed.(1)Select S sample points from the reference point cloud , ensuring that the sampled points have different FPFH characteristics as far as possible, and their pairing distance is greater than the minimum value set by the user.(2)For S sample points, the points satisfying similar FPFH characteristics are found in the point cloud to be registered, and stored in a list. Randomly select a point from these similar points as the corresponding point.(3)After a group of corresponding points is obtained, the rigid body transformation matrix between the corresponding points is calculated, the result of which is as follows:(4)The performance of the current registration transformation is judged by calculating the “distance error sum” function after the corresponding point transformation. The distance error sum function is often expressed by the Huber function and is denoted as , wherewhere is a predetermined value, and is the distance difference after the transformation of the corresponding points of group i. The ultimate goal of registration is to find a set of optimal transformations in all transformations to minimize the value of the error function. Then, the transformation is the final transformation matrix of coarse registration of the roadway point cloud, which is used as the initial matrix of the ICP precision registration algorithm.

4.2.4. Precise Registration of Point Cloud in Excavation Roadway

Find a point in the reference point cloud , locate the corresponding in the to be registered, so that  = min. Transmit the final rotation matrix and translation matrix results obtained by formula (11) to the ICP algorithm. Search the roadway point cloud via the KD-tree to accelerate the precision registration speed of the ICP algorithm. If the reference point cloud and the point cloud are registered to obtain m groups of corresponding points, the error function is formed as shown in the following formula:

The point cloud precise registration will be completed when equation (13) obtains the optimal value. Because of the nonlinearity of the objective function, it is difficult to obtain the optimal estimation parameters by the linear least square method or other multivariable function extreme solution methods. In this paper, the Gauss–Newton is used to optimize the parameters in the accurate registration of roadway point clouds to avoid the ICP algorithm falling into a local optimum.

5. Experimental Data Analysis

5.1. Point Cloud Preprocessing Experiment

Through the European point cloud clustering segmentation algorithm, the three groups of point clusters with outliers are preprocessed by the point cloud. Experimental results are shown in Figure 7 and the experimental data are shown in Table 2.

From the experimental results in Figure 7 and Table 2, it can be found that European point cloud clustering segmentation can effectively eliminate outliers and unnecessary ground points, reduce the number of point clouds without destroying the feature structure of point clouds, improve the calculation speed, and reduce the storage space for subsequent point cloud registration. Therefore, this method is suitable for the pretreatment of the point cloud in the excavation roadway.

5.2. Optimization of FPFH + ICP Point Cloud Registration Parameters

From formula (10), it can be seen that FPFH features are affected by their domain radius. That is, when the selected radius is larger, it may cause mixed mismatched point clouds, which makes the point cloud features unclear and causes a long registration time. When the selected radius is small, it is easy to be disturbed by external factors of the roadway, and it is difficult to form a perfect feature with few feature points. Therefore, it is of significance to choose appropriate parameters for registration results. In this paper, a large number of experiments are carried out, considering the factors of registration time and accuracy. That is, when the radius of the FPFH feature field is 0.1 m–1 m, the registration effect is the best. To ensure the domain radius interval, the parameters with a domain radius of 0.01 m, 0.1 m, 0.15 m, 0.5 m, and 1 m are fully selected for experimental analysis. The maximum iteration times of ICP are 50 times, and the allowable difference of the transformation matrix is 1e − 10. The experimental results of selecting a room1.pcd file as a reference point cloud (about 10,496-point clouds in the blue point cluster) and a room2.pcd file as the point cloud to be registered (about 10,486-point clouds in the green point cluster) are shown in Table 3.

According to the experimental results in Table 3, it can be seen that when the FPFH feature domain radius is set to 0.15 m, the FPFH + ICP registration algorithm takes the shortest time and has the highest registration accuracy with the best registration effect.

5.3. Comparative Experiment of Point Cloud Registration Algorithms
5.3.1. FPFH + ICP Algorithm Registration Results

The registration experiment of the same roadway data in Figure 4 was carried out by the FPFH + ICP algorithm. The scenes with concave and convex roadway wall, the scene with mining cable, and the scene with an obstacle are selected respectively. The matching results of the FPFH + ICP algorithm are shown in Figure 8.

5.3.2. Compared with the Traditional Point Cloud Registration Algorithm

To verify the practicability of the FPFH + ICP algorithm point cloud registration in the construction of environment models of excavation roadways, a comparative experiment of point cloud registration is carried out in this section. The maximum number of ICP iterations is 50, the allowable difference of the transformation matrix is 1e − 10, the resolution of the NDT cell is set to 1.5 m, and the radius of the FPFH feature neighborhood is set to 0.15 m. Translation errors of each algorithm in the X axis, Y axis, and Z axis are shown in Figure 9.

From the above experimental results and Figure 4, it can be concluded that there are a large number of ropes and electromechanical equipment in the closed and unstructured roadway environment, and the roadway wall is uneven. The matching method based on scanning cannot achieve a breakthrough in accuracy in this environment. The FPFH point feature histogram descriptor can describe the boundary, corner, plane point cloud, and other feature information inside the roadway, which greatly improves the registration accuracy to the millimeter level.

After point cloud registration, the root mean square error of two groups of point clusters is used to reflect the modeling accuracy of point cloud registration; that is, the smaller the error, the higher the accuracy. The root means the square error is shown in the following equation:where K represents the number of point cloud pairs in a point cluster, represents the Euclidean distance of two points in a pair of point clouds after point cloud registration, represents the true value of the Euclidean distance between two points in a point cloud pair after point cloud registration. Ideally, the truth value is 0. Select 20 frames of point clouds. The positioning accuracy analysis of each point cloud registration algorithm is shown in Figure 10.

According to the experimental data, the root mean square errors of the ICP algorithm are 11 cm, 12 cm, and 13 cm, and the initial position of the point cloud is required to be high. When modeling a long-closed roadway, it is easy to fall into a local optimum and lead to structural deformation of the model. The root mean square errors of the NDT + ICP algorithm are 10 cm, 9 cm, and 11 cm, but the maximum and minimum of the errors are quite different. As a rough registration algorithm, the NDT algorithm is not suitable for large amounts of noisy data, which leads to poor stability of registration results, large fluctuations of data, and low accuracy of registration results. The average error of the FPFH + ICP is about 5 cm, and the minimum accuracy reaches 2 cm, which is far better than the traditional ICP algorithm and more stable than the NDT + ICP algorithm. By comparison, the FPFH + ICP point cloud registration method is better in point cloud registration accuracy and more suitable for the coal mine roadway environment.

6. Conclusion

Aiming at the low efficiency and high error of point cloud model registration with low feature recognition of excavation roadways, this paper proposes the FPFH + ICP point cloud registration method, which reduces outliers in the roadway by the European point cloud clustering segmentation preprocessing method, simplifies the number of point clouds in the excavation roadway, and improves the registration efficiency. SAC-IA coarse registration realizes the initial rigid body transformation of the point cloud and avoids the algorithm falling into a local optimum. Finally, the point cloud of the excavation roadway is registered quickly and accurately by ICP. The above experiments show that the algorithm has a good registration effect and improves the point cloud registration accuracy of excavation roadways and also improves the registration efficiency in the fine registration process. However, the algorithm still has a slightly longer registration time and needs to be further optimized when it is used on other occasions, which is also the place where this algorithm should be improved.

Data Availability

This paper takes the tunnel point cloud of Beijing Mentougou Jingxi Coal Mine as the research object.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Authors’ Contributions

For this research article, JJ. Y. is in charge of methodology. WJ. L. is in charge of writing, reviewing, and formal analysis. YC. Z. is in charge of data curation and field inspection. BS. C. and R. Z. are in charge of testing and verification. M. W provides the funding for project. All authors have read and agreed to the published version of the manuscript.


The study was approved by the China University of Mining and Technology (Beijing). This research was funded by the National Natural Science Foundation of China (Grant no. 52104169), National Natural Science Foundation of China (Grant no. 51874308), Shanxi Special Fund for Science and Technology (Grant no. 20181102027), and National Natural Science Foundation of China (Grant no. 52121003).