An Accurate Digital Subsidence Model for Deformation Detection of Coal Mining Areas Using a UAV-Based LiDAR

: Coal mine surface subsidence detection determines the damage degree of coal mining, which is of great importance for the mitigation of hazards and property loss. Therefore, it is very important to detect deformation during coal mining. Currently, there are many methods used to detect deformations in coal mining areas. However, with most of them, the accuracy is difﬁcult to guarantee in mountainous areas, especially for shallow seam mining, which has the characteristics of active, rapid, and high-intensity surface subsidence. In response to these problems, we made a digital subsidence model (DSuM) for deformation detection in coal mining areas based on airborne light detection and ranging (LiDAR). First, the entire point cloud of the study area was obtained by coarse to ﬁne registration. Second, noise points were removed by multi-scale morphological ﬁltering, and the progressive triangulation ﬁltering classiﬁcation (PTFC) algorithm was used to obtain the ground point cloud. Third, the DEM was generated from the clean ground point cloud, and an accurate DSuM was obtained through multiple periods of DEM difference calculations. Then, data mining was conducted based on the DSuM to obtain parameters such as the maximum surface subsidence value, a subsidence contour map, the subsidence area, and the subsidence boundary angle. Finally, the accuracy of the DSuM was analyzed through a comparison with ground checkpoints (GCPs). The results show that the proposed method can achieve centimeter-level accuracy, which makes the data a good reference for mining safety considerations and subsequent restoration of the ecological environment.


Introduction
Deformation detection has been defined as the identification of geometric state differences based on multiple periods of data capturing. The detection of the surface subsidence of coal mining areas is a part of deformation detection and has become a hot topic to mitigate hazards and property loss.
Generally, the surface subsidence caused by coal mining is the main source of danger for the destruction of buildings and structures, inevitably causing surface collapse and environmental damage [1][2][3]. Therefore, scholars have adopted various methods to observe the surface subsidence of coal mining areas. The traditional geodetic method detects surface subsidence with fixed points on the ground [4,5]. Although high-precision data can be obtained, it is still point-to-point acquisition, which is inefficient and expensive. Moreover, this method only measures local subsidence, and full coverage of a mining area cannot be obtained. Recently, Shi et al. [6][7][8][9] tried to detect surface deformation with the method of interferometry synthetic aperture radar (InSAR), which can obtain accurate vertical displacement measurements. However, the speed of the surface subsidence of a shallow coal seam is relatively fast, and a long period of SAR satellite observation easily causes incoherence of SAR images [8,9]. UAV oblique photogrammetry can obtain of wavy dunes. The study area has a mid-temperate, continental, semi-arid climate. There are regular droughts. There is little rain (mean annual precipitation of 435.7 mm; the annual average evaporation is 1774.1 mm, which is 4−5 times the precipitation), sparse vegetation, and serious soil erosion. The ecological environment is very fragile [34]. Thus, it is of great significance to monitor the deformation of the coal mining areas. In this study, a local part of a working face was selected as the study area. UAV−based LiDAR was used to obtain the surface deformation data, and checkpoints were obtained by leveling.

Mining and Geological Conditions
We took the 431,301 working face of Liangshuijing Coal Mine as the study area. The working face, which mined 3−4 coal seams, is stable and has flat seams. The coal seam thickness is 1.17−1.43 m, the average thickness is 1.3 m, and the average depth of coal seam is 138 m. The working face adopts the longwall, fully mechanized, and full−seam mining method, the roof of which is managed with the all fall method. The ground is covered by loose sandy soil with a strong flow characteristic. The entire study area is 858 m long from east to west and 466 m long from north to south, giving an area of 399,431 m 2 . The location of the study area and the working face is shown in Figure 1.

Reference Data
Our data include airborne laser scanning data and geodetic data. The geodetic data were obtained by laying observation piles on the ground. Then, static and dynamic GNSS [35], total station observation, and precision leveling were used to obtain geodetic data. In the study, the checkpoints include plane (RTK observation) and elevation (leveling observation) data.

LiDAR Data
UAV−based LiDAR was used to collect initial data, and the DSuM was generated by the data for detecting the deformation. The endurance time of UAV in each flight was 35 min, and the effective working time was about 20 min. The laser scanner was a RIEGL miniVUX-1UAV with scanning range of 250 m. The same aerial survey parameters were used to ensure each data acquisition would have the same system error. The parameters included flying height, flying speed, laser scanning speed, pulse emission frequency, and Yushen Coal Field is located in the transition zone of the Muus Desert and the Loess Plateau. The eastern part consists of a loess ridge and a valley; the western part consists of wavy dunes. The study area has a mid-temperate, continental, semi-arid climate. There are regular droughts. There is little rain (mean annual precipitation of 435.7 mm; the annual average evaporation is 1774.1 mm, which is 4-5 times the precipitation), sparse vegetation, and serious soil erosion. The ecological environment is very fragile [34]. Thus, it is of great significance to monitor the deformation of the coal mining areas. In this study, a local part of a working face was selected as the study area. UAV-based LiDAR was used to obtain the surface deformation data, and checkpoints were obtained by leveling.

Mining and Geological Conditions
We took the 431,301 working face of Liangshuijing Coal Mine as the study area. The working face, which mined 3-4 coal seams, is stable and has flat seams. The coal seam thickness is 1.17-1.43 m, the average thickness is 1.3 m, and the average depth of coal seam is 138 m. The working face adopts the longwall, fully mechanized, and full-seam mining method, the roof of which is managed with the all fall method. The ground is covered by loose sandy soil with a strong flow characteristic. The entire study area is 858 m long from east to west and 466 m long from north to south, giving an area of 399,431 m 2 . The location of the study area and the working face is shown in Figure 1.

Reference Data
Our data include airborne laser scanning data and geodetic data. The geodetic data were obtained by laying observation piles on the ground. Then, static and dynamic GNSS [35], total station observation, and precision leveling were used to obtain geodetic data. In the study, the checkpoints include plane (RTK observation) and elevation (leveling observation) data.

LiDAR Data
UAV-based LiDAR was used to collect initial data, and the DSuM was generated by the data for detecting the deformation. The endurance time of UAV in each flight was 35 min, and the effective working time was about 20 min. The laser scanner was a RIEGL miniVUX-1UAV with scanning range of 250 m. The same aerial survey parameters were used to ensure each data acquisition would have the same system error. The parameters included flying height, flying speed, laser scanning speed, pulse emission frequency, and weather conditions. The flying height was 50 m, and the flying speed was 8 m/s. The laser scanning speed was 100 lines/s, and the pulse emission frequency was 100 kHz. Each laser beam had five echoes, and data collection was performed under the same conditions on different dates. Two groups of point clouds were obtained by UAV-based LiDAR.
The diagram of the UAV-based LiDAR data collection and the mining process of the underground working face is shown in Figure 2. The nth LiDAR data collection was carried out to obtain the initial shape of the surface. After a period of time, the coal seams were mined, which caused surface deformation. The m-th LiDAR data collection was carried out to obtain the current surface morphology. With the calculation of two phases of LiDAR data, the surface shape change caused by underground working face mining was obtained, which is called the subsidence basin. weather conditions. The flying height was 50 m, and the flying speed was 8 m/s. The laser scanning speed was 100 lines/s, and the pulse emission frequency was 100 kHz. Each laser beam had five echoes, and data collection was performed under the same conditions on different dates. Two groups of point clouds were obtained by UAV-based LiDAR. The diagram of the UAV-based LiDAR data collection and the mining process of the underground working face is shown in Figure 2. The nth LiDAR data collection was carried out to obtain the initial shape of the surface. After a period of time, the coal seams were mined, which caused surface deformation. The m−th LiDAR data collection was carried out to obtain the current surface morphology. With the calculation of two phases of LiDAR data, the surface shape change caused by underground working face mining was obtained, which is called the subsidence basin.  Figure 2 shows the process of mining subsidence. The layered structure on the right of Figure 2 represents the geological structure, which includes the ground surface, sandy soil layer, rock layer, coal seam, and other strata. The main part of Figure 2 shows the relationship between coal seam mining and surface deformation. The surface deformation was different at different mining locations. The black arrow indicates the mining direction of the coal seam. Numbers 1 and 2 indicate the locations of the shearer on different dates. The continuous mining produced a goaf in the coal seam, caused surface deformation, and formed a subsidence basin on the surface.
The data collection dates were 7 November 2020 and 19 May 2021, and the corresponding coal seam mining locations are shown in Figure 3. The original point cloud statistics are shown in Table 1.   Figure 2 shows the process of mining subsidence. The layered structure on the right of Figure 2 represents the geological structure, which includes the ground surface, sandy soil layer, rock layer, coal seam, and other strata. The main part of Figure 2 shows the relationship between coal seam mining and surface deformation. The surface deformation was different at different mining locations. The black arrow indicates the mining direction of the coal seam. Numbers 1 and 2 indicate the locations of the shearer on different dates. The continuous mining produced a goaf in the coal seam, caused surface deformation, and formed a subsidence basin on the surface.
The data collection dates were 7 November 2020 and 19 May 2021, and the corresponding coal seam mining locations are shown in Figure 3. The original point cloud statistics are shown in Table 1.

Ground Checkpoints
The fixed observation piles were used as the objects to obtain the partial subsidence value of the ground surface, because there were few fixed markers in the study area. The observation piles were made by precasting concrete, whose structure is shown in Figure  4b, and the locations of the observation piles are shown as black triangles in Figure 3; they are also called ground checkpoints. The processing of GCPs included the design of ground checkpoints' locations, observation pile burying, geodetic surveying, measurement calculations, and subsidence polyline drawing, as shown on the right-hand side of Figure 5. In order to verify the measurement accuracy of the point cloud, DEM, and DSuM, the observation piles were measured on 7 November 2020 and 19 May 2021, the same dates as the UAV−based LiDAR scanning. The results of GCPs were also used to calculate the max subsidence value and boundary angle.

Ground Checkpoints
The fixed observation piles were used as the objects to obtain the partial subsidence value of the ground surface, because there were few fixed markers in the study area. The observation piles were made by precasting concrete, whose structure is shown in Figure 4b, and the locations of the observation piles are shown as black triangles in Figure 3; they are also called ground checkpoints. The processing of GCPs included the design of ground checkpoints' locations, observation pile burying, geodetic surveying, measurement calculations, and subsidence polyline drawing, as shown on the right-hand side of Figure 5. In order to verify the measurement accuracy of the point cloud, DEM, and DSuM, the observation piles were measured on 7 November 2020 and 19 May 2021, the same dates as the UAV-based LiDAR scanning. The results of GCPs were also used to calculate the max subsidence value and boundary angle.

Data Processing
When the shearer mined different positions on the working face, the two periods of UAV-based LiDAR surveying were carried out to obtain the surface morphology of the working face in different periods, as shown in Figures 2 and 3. The GCPs were mainly used to verify the accuracy of the point cloud, DEM, and DSuM, which were all generated directly or indirectly based on LiDAR. Similarly to previous related studies [7,17,[35][36][37], we used geodetic data as a reference to analyze the accuracy of the point cloud, DEM, and DSuM. The multi-scale morphological algorithm [38] and progressive triangulation filter algorithm [39,40] were chosen for filtering the point cloud, and the Kriging interpolation algorithm [41] was chosen for generating DEM. This paper also analyzes the influence of resolution on DEM accuracy and proposes a new model termed DSuM based on the optimal Remote Sens. 2022, 14, 421 6 of 19 resolution of DEM to detect the deformation of coal mining areas. The main purposes of coal mining subsidence monitoring are subsidence value acquisition and boundary angle calculation. This section introduces the calculation methods for the subsidence value and boundary angle, and the process of generating DSuM from the UAV-based LiDAR data.

Data Processing
When the shearer mined different positions on the working face, the two periods of UAV-based LiDAR surveying were carried out to obtain the surface morphology of the working face in different periods, as shown in Figures 2 and 3. The GCPs were mainly

Data Processing
When the shearer mined different positions on the working face, the two periods of UAV-based LiDAR surveying were carried out to obtain the surface morphology of the working face in different periods, as shown in Figures 2 and 3. The GCPs were mainly

Subsidence Value and Boundary Angle
Subsidence data analysis includes data preprocessing, data calculation, subsidence polyline drawing, and boundary angle calculation. The data preprocessing is used to remove error data caused by human factors; data calculation is to obtain the difference between the two periods' data. Then, subsidence polyline diagram is drawn, and boundary angles are obtained. The subsidence of the observation points is calculated by Formula (1): where ω n represents the subsidence value of point n; H 0 n and H m n represent the elevation of point n during the first and the m-th observations.
The boundary angle is used to determine the range and boundary of the subsidence caused by coal mining. Due to the observation error and the seasonal variations of the soil, an observation point with a subsidence value close to zero could not be accurately determined, so points with subsidence values of 10 and 100 mm were taken as the boundaries of the subsidence basin, as shown in Figure 4. Boundary angles calculated by different boundary thresholds are different. The boundary angle is calculated as follows: where L 1 represents the horizontal distance between the boundary of the underground goaf and the strike or incline boundary of the subsidence basin, H 0 represents the average depth of the coal seam, and δ 0 represents the boundary angle.

UAV-Based LiDAR Data Processing
The data processing of UAV-based LiDAR mainly included data checking, point cloud processing, and interpolation calculation. Data checking was to remove erroneous data; point cloud processing included registration, denoising, and filtering; the purpose of interpolation calculation was to generate the DEM. After data checking, the produced raw data were converted to point clouds. The point cloud was roughly registered by pose data generated by GNSS/INS integration, and then fine registration was carried out based on the iterative closest point (ICP) algorithm. Subsequently, we removed the noise points with the morphological method [42]. Then, the PTFC algorithm was used to obtain the ground point cloud, and the DEM was generated from the ground point cloud. Finally, an accurate DSuM was obtained through multiple periods of DEM difference calculation.
We compared the GCPs and the point cloud data at the same location to verify the point cloud accuracy, which is a direct accuracy verification method. Subsequently, we evaluated the accuracy of the DEM by comparing it with GCPs. The Kriging interpolation algorithm was used to generate the DEM, and the most suitable grid size was determined by trial and error.
There are two methods to evaluate the accuracy of a DEM. The elevation error statistical analysis method compares the DEM with the reference DEM or checkpoints; the logical analysis method is a qualitative method for overall accuracy evaluation, including a visual interpretation method, contour analysis, visual analysis, and other methods [41]. We adopted the method of comparing DEM to checkpoints. The differences in elevation between the checkpoints and the DEM with different resolutions were calculated, and the results were displayed by box plots. Meanwhile, the mean error (ME), the mean absolute error (MAE), and the root mean square error (RMSE), were calculated [19]: where R m represents the value of the DEM and Z m represents the value of a checkpoint. Since directly calculating the difference between two point clouds is difficult, we used the difference between any two DEMs created by the point cloud to represent the ground surface subsidence deformation during the period. The difference between the two DEMs is called the DSuM. Finally, the accuracy of the DSuM was analyzed by GCPs. The analysis of a point cloud, DEM, and DSuM is shown in Figure 5.

Pipeline of DSuM
The difference between DEMs obtained on any two different dates can be calculated to obtain the ground surface elevation change during this period, in our case, the ground surface subsidence value caused by underground coal mining. In this paper, the difference between the DEMs was calculated, and represented by DSuM, which is a digital model that represents the value of ground subsidence in an ordered array of values. Each pixel value of the DSuM represents the subsidence value of the pixel location caused by coal mining. The schematic diagram of DSuM generation is shown in Figure 6, which is a digital expression of whole ground surface subsidence value. A DSuM can be calculated from any two DEMs obtained on different dates, and can represent the subsidence value of any position in the ground surface during the mining processing. Finally, data mining was conducted based on the DSuM to obtain outputs such as the maximum surface subsidence value, a subsidence contour map, the subsidence range, the subsidence area, and the subsidence boundary angle, as shown in Figure 5.
where represents the value of the DEM and represents the value of a checkpoint. Since directly calculating the difference between two point clouds is difficult, we used the difference between any two DEMs created by the point cloud to represent the ground surface subsidence deformation during the period. The difference between the two DEMs is called the DSuM. Finally, the accuracy of the DSuM was analyzed by GCPs. The analysis of a point cloud, DEM, and DSuM is shown in Figure 5.

Pipeline of DSuM
The difference between DEMs obtained on any two different dates can be calculated to obtain the ground surface elevation change during this period, in our case, the ground surface subsidence value caused by underground coal mining. In this paper, the difference between the DEMs was calculated, and represented by DSuM, which is a digital model that represents the value of ground subsidence in an ordered array of values. Each pixel value of the DSuM represents the subsidence value of the pixel location caused by coal mining. The schematic diagram of DSuM generation is shown in Figure 6, which is a digital expression of whole ground surface subsidence value. A DSuM can be calculated from any two DEMs obtained on different dates, and can represent the subsidence value of any position in the ground surface during the mining processing. Finally, data mining was conducted based on the DSuM to obtain outputs such as the maximum surface subsidence value, a subsidence contour map, the subsidence range, the subsidence area, and the subsidence boundary angle, as shown in Figure 5.  After point cloud data collection, the DEMs were generated from clean ground point clouds, and the DSuM was obtained through multiple periods of DEM difference calculation. Then, coal mining subsidence was analyzed based on the DSuM. We adopted point, line, and area analysis from different perspectives. The point analysis used the method of simulating traditional monitoring to extract subsidence values from points, which involves drawing the strike and incline subsidence polyline and obtaining the max subsidence value. The line analysis easily extracted high-density points from the DSuM, which was used to extract points in the strike and incline direction in this paper. It is robust for discrete monitoring points, and it can plot the subsidence curve graph and calculate the strike and incline boundary angle. The area analysis can analyze the whole subsidence area. It was used to calculate subsidence area, maximum subsidence value, and the ratio of subsidence area.
Finally, data mining was carried out based on the DSuM to obtain the maximum subsidence value, subsidence area, subsidence distribution of surface, boundary angle, and other parameters for later analysis.

GCPs Analysis
The GCPs results were obtained by Formula (1) based on the geodetic data. The results indicate that the maximum ground subsidence value of the working face was 1826 mm, and the minimum subsidence value was 0 mm, during 7 November 2020 and 19 May 2021, respectively. Taking ground observation point A176 as the origin, the subsidence values of all GCPs are shown as blue squares in Figure 7.
strike and incline boundary angle. The area analysis can analyze the whole subsidence area. It was used to calculate subsidence area, maximum subsidence value, and the ratio of subsidence area.
Finally, data mining was carried out based on the DSuM to obtain the maximum subsidence value, subsidence area, subsidence distribution of surface, boundary angle, and other parameters for later analysis.

GCPs Analysis
The GCPs results were obtained by Formula (1) based on the geodetic data. The results indicate that the maximum ground subsidence value of the working face was 1826 mm, and the minimum subsidence value was 0 mm, during 7 November 2020 and 19 May 2021, respectively. Taking ground observation point A176 as the origin, the subsidence values of all GCPs are shown as blue squares in Figure 7. The GCPs covered the entire subsidence area, including the non−subsidence area and the maximum subsidence area, and can reflect the surface subsidence in a period of time. The blue polyline in Figure 7 shows that the subsidence boundary with the threshold of 10 mm located between A205 and A206, and the horizontal distances of A205 and A206 from the boundary of goaf were 121 and 136 m, respectively. The average mining depth of the coal seam was 138 m. According to Formula (2), the mining subsidence boundary angle with subsidence threshold of 10 mm was between 45° and 48.7°. The average was 46.9°. Similarly, the boundary angle with a subsidence threshold of 100 mm was 61.2°.

Accuracy Assessment
The ground point cloud could be obtained after filtering the original point cloud. Since a ground point cloud is composed of a series of discrete points, we took the average The GCPs covered the entire subsidence area, including the non-subsidence area and the maximum subsidence area, and can reflect the surface subsidence in a period of time. The blue polyline in Figure 7 shows that the subsidence boundary with the threshold of 10 mm located between A205 and A206, and the horizontal distances of A205 and A206 from the boundary of goaf were 121 and 136 m, respectively. The average mining depth of the coal seam was 138 m. According to Formula (2), the mining subsidence boundary angle with subsidence threshold of 10 mm was between 45 • and 48.7 • . The average was 46.9 • . Similarly, the boundary angle with a subsidence threshold of 100 mm was 61.2 • .

Accuracy Assessment
The ground point cloud could be obtained after filtering the original point cloud. Since a ground point cloud is composed of a series of discrete points, we took the average elevation values of points near the GCP to calculate the differences between ground point cloud and GCPs, and the distances from the selected points to the GCPs had to be less than a set threshold. The results indicate that the RMSE of elevation was 60.6 mm on 7 November 2020, and 59.9 mm on 19 May 2021, and the results are shown in Table 2. Next, we calculated the accuracy of the DEM and analyzed the influence of resolution on DEM accuracy. The differences between DEM and GCPs are shown in Figure 8. As resolution increases, the average error remains basically unchanged, but the error distribution becomes more discrete. The error distribution of the resolution of 0.1 m was most concentrated, and all errors were less than 100 mm.
Next, we calculated the accuracy of the DEM and analyzed the influence of resolution on DEM accuracy. The differences between DEM and GCPs are shown in Figure 8. As resolution increases, the average error remains basically unchanged, but the error distribution becomes more discrete. The error distribution of the resolution of 0.1 m was most concentrated, and all errors were less than 100 mm. The difference results of DEM from data at different resolutions and GCPs are shown in Table 3. The results include ME, MAE, and RMSE, and indicate that with different resolutions, the maximum ME is 323 mm and the minimum is -34 mm; the maximum MAE is 537 mm and the minimum is 74 mm; and the maximum RMSE is 97 mm.  The difference results of DEM from data at different resolutions and GCPs are shown in Table 3. The results include ME, MAE, and RMSE, and indicate that with different resolutions, the maximum ME is 323 mm and the minimum is -34 mm; the maximum MAE is 537 mm and the minimum is 74 mm; and the maximum RMSE is 97 mm. The ME did not change, and the optimal state was in the range of 0.05-1 m, but it increased when the grid size was larger than 1 m. The change in MAE was similar to that of the ME. Similarly, the RMSE reached an equilibrium state within 0.05-1 m. When the resolution was greater than 1 m, the RMSE gradually increased as the grid size increased, and it increased sharply when the resolution became greater than 1 m.
According to Figure 8, the data with 0.1 m grid size are the most concentrated. Therefore, we selected 0.1 m as the parameter for conducting the study to ensure research accuracy.
Finally, the accuracy of the DSuM was analyzed. The DSuM of the study area was obtained by the subtraction of DEM in 7 November 2020 from DEM in 19 May 2021, as shown in Figure 9. Different colors represent different subsidence values in the DSuM. Additionally, the main area of subsidence caused by coal mining is mainly distributed in the goaf and its surroundings.
The DSuM contains plane coordinate (X, Y) and subsidence value. The subsidence value of any position can be acquired. The subsidence monitoring accuracy of UAV-based LiDAR can be acquired by comparing the subsidence value of the DSuM with those of the GCPs. The results are shown in Figure 10.

racy.
Finally, the accuracy of the DSuM was analyzed. The DSuM of the study area was obtained by the subtraction of DEM in 7 November 2020 from DEM in 19 May 2021, as shown in Figure 9. Different colors represent different subsidence values in the DSuM. Additionally, the main area of subsidence caused by coal mining is mainly distributed in the goaf and its surroundings. The DSuM contains plane coordinate (X, Y) and subsidence value. The subsidence value of any position can be acquired. The subsidence monitoring accuracy of UAV-based LiDAR can be acquired by comparing the subsidence value of the DSuM with those of the GCPs. The results are shown in Figure 10. Figure 7 shows the maximum subsidence value and minimum subsidence value. The subsidence polyline trends obtained by DSuM and GCPs were basically consistent. The maximum subsidence value calculated by DSuM was 1860 mm, and the maximum subsidence value of GCPs was 1872 mm. The inflection points of the subsidence polyline were consistent. After removing the abnormal subsidence value of ground checkpoint (A 203), the maximum and minimum of the subsidence difference's absolute values were 82 and 3 mm, the average difference was 45 mm, and 75% of the difference was within 50 mm. All of the differences were less than 100 mm, and their distribution was relatively uniform. The absolute values of differences among the point cloud, DEM, and DSM are shown in Figure 10.

Analysis of DSuM
In order to comprehensively verify the accuracy of the DSuM, we conducted point, line, and area analysis based on the DSuM.

Point Analysis
The traditional method of coal mining subsidence monitoring is geodetic surveying by a fixed observation pile, which is a point survey method. For the point analysis of DSuM, we used the method of simulating traditional monitoring to extract subsidence values point by point. The locations of monitoring points were set according to requirements of relevant specifications. The interval between adjacent points was 15 m, the total number was 50, and the whole surface subsidence change area was covered, as shown in Figure 11, by black triangles.  Figure 7 shows the maximum subsidence value and minimum subsidence value. The subsidence polyline trends obtained by DSuM and GCPs were basically consistent. The maximum subsidence value calculated by DSuM was 1860 mm, and the maximum subsidence value of GCPs was 1872 mm. The inflection points of the subsidence polyline were consistent. After removing the abnormal subsidence value of ground checkpoint (A 203), the maximum and minimum of the subsidence difference's absolute values were 82 and 3 mm, the average difference was 45 mm, and 75% of the difference was within 50 mm. All of the differences were less than 100 mm, and their distribution was relatively uniform. The absolute values of differences among the point cloud, DEM, and DSM are shown in Figure 10.

Analysis of DSuM
In order to comprehensively verify the accuracy of the DSuM, we conducted point, line, and area analysis based on the DSuM.

Point Analysis
The traditional method of coal mining subsidence monitoring is geodetic surveying by a fixed observation pile, which is a point survey method. For the point analysis of DSuM, we used the method of simulating traditional monitoring to extract subsidence values point by point. The locations of monitoring points were set according to requirements of relevant specifications. The interval between adjacent points was 15 m, the total number was 50, and the whole surface subsidence change area was covered, as shown in Figure 11, by black triangles.
The traditional method of coal mining subsidence monitoring is geodetic surveying by a fixed observation pile, which is a point survey method. For the point analysis of DSuM, we used the method of simulating traditional monitoring to extract subsidence values point by point. The locations of monitoring points were set according to requirements of relevant specifications. The interval between adjacent points was 15 m, the total number was 50, and the whole surface subsidence change area was covered, as shown in Figure 11, by black triangles. The subsidence value of each monitoring point was extracted by DSuM and is shown in the form of a polyline graph in Figure 12. The results indicate that the average maximum value of the subsidence basin is approximately 1650 mm, and the maximum single point subsidence is 1761 mm. The shape of the strike and incline polyline graph is the same as that of the standard mining subsidence polyline graph, whose shape is like a basin or a half-basin. The results can reflect the elevation deformation of the surface caused by coal mining, and meet the requirements of mining subsidence monitoring.

Line Analysis
Line analysis is robust for discrete monitoring points. Therefore, it is necessary to obtain high-density points on the observation line. We extracted high-density points from DSuM in the strike and incline directions and plotted the subsidence curve graph in We extracted monitoring points from DSuM in an interval of 0.1 m, and drew the points on graphs, as shown in Figure 13a,b−the green points. Finally, we extracted and drew the strike and incline curves, as shown in Figure 13 as the red curves. The results indicate that the average maximum value of the subsidence basin is approximately 1650 mm, and the maximum single point subsidence is 1761 mm. The shape of the strike and incline polyline graph is the same as that of the standard mining subsidence polyline graph, whose shape is like a basin or a half-basin. The results can reflect the elevation deformation of the surface caused by coal mining, and meet the requirements of mining subsidence monitoring.

Line Analysis
Line analysis is robust for discrete monitoring points. Therefore, it is necessary to obtain high-density points on the observation line. We extracted high-density points from DSuM in the strike and incline directions and plotted the subsidence curve graph in Figure 13. The locations of strike and cline lines are shown in Figure 11  The strike curve is shaped like a basin, conforming to the characteristics of mining subsidence, and it shows that a super-full mining state in the strike direction was reached. The incline observation line is in a valley state and is symmetrically distributed due to the small width of the working face. According to the strike and incline curves, we found that the maximum subsidence value is about 1700 mm. Since the average accuracy of UAVbased LiDAR data cannot reach 10 mm, 100 mm was used as the subsidence boundary threshold to extract and calculate the boundary angle by Formula (2). In the strike direction, the horizontal distance between the subsidence boundary and the goaf boundary is 68 m; in the incline direction, the horizontal distance is 60 m, and the average coal mining depth is 138 m. Therefore, when the boundary subsidence threshold is 100 mm, the strike and cline boundary angle are 63.8° and 66.5°.

Area Analysis
In order to accurately express the subsidence values of different areas and the subsidence area with different subsidence values, it was necessary to perform area statistical analysis based on the DSuM. Therefore, we classified and counted the DSuM at the interval of 100 mm. The results are shown in Figure 14. We extracted monitoring points from DSuM in an interval of 0.1 m, and drew the points on graphs, as shown in Figure 13a,b-the green points. Finally, we extracted and drew the strike and incline curves, as shown in Figure 13 as the red curves.
The strike curve is shaped like a basin, conforming to the characteristics of mining subsidence, and it shows that a super-full mining state in the strike direction was reached. The incline observation line is in a valley state and is symmetrically distributed due to the small width of the working face. According to the strike and incline curves, we found that the maximum subsidence value is about 1700 mm. Since the average accuracy of UAV-based LiDAR data cannot reach 10 mm, 100 mm was used as the subsidence boundary threshold to extract and calculate the boundary angle by Formula (2). In the strike direction, the horizontal distance between the subsidence boundary and the goaf boundary is 68 m; in the incline direction, the horizontal distance is 60 m, and the average coal mining depth is 138 m. Therefore, when the boundary subsidence threshold is 100 mm, the strike and cline boundary angle are 63.8 • and 66.5 • .

Area Analysis
In order to accurately express the subsidence values of different areas and the subsidence area with different subsidence values, it was necessary to perform area statistical analysis based on the DSuM. Therefore, we classified and counted the DSuM at the interval of 100 mm. The results are shown in Figure 14. The upper half of Figure 14 is the subsidence classification map of the study area. It shows the locations of different subsidence values, and the different colors represent different subsidence values. The lower part of Figure 14 is a histogram which represents the number of pixels for each subsidence value. According to the classification map, there are some areas with negative subsidence values, especially the red circled area, which are caused by man-excavation, road construction, etc. Figure 14 shows that the subsidence area is distributed in an oval shape. The main subsidence area is located at the center of the study area. Due to the influence of natural rainfall, there is a large area with subsidence values less than 100 m. In particular, there are very few areas where the subsidence value is less than −100 mm, as shown in the red circle in Figure 14. In addition, according to the histogram in Figure 14, we found that the larger the subsidence value, the smaller the number of the corresponding grids. Next, we counted the number of grids with different subsidence values, and calculated the areas of different subsidence values. The results are shown in Table 4.  The upper half of Figure 14 is the subsidence classification map of the study area. It shows the locations of different subsidence values, and the different colors represent different subsidence values. The lower part of Figure 14 is a histogram which represents the number of pixels for each subsidence value. According to the classification map, there are some areas with negative subsidence values, especially the red circled area, which are caused by man-excavation, road construction, etc. Figure 14 shows that the subsidence area is distributed in an oval shape. The main subsidence area is located at the center of the study area. Due to the influence of natural rainfall, there is a large area with subsidence values less than 100 m. In particular, there are very few areas where the subsidence value is less than −100 mm, as shown in the red circle in Figure 14. In addition, according to the histogram in Figure 14, we found that the larger the subsidence value, the smaller the number of the corresponding grids. Next, we counted the number of grids with different subsidence values, and calculated the areas of different subsidence values. The results are shown in Table 4. Table 4 indicates that areas with a subsidence value of less than 100 mm occupied more than 55% of the study area, in which most of the subsidence values are caused by external factors. In order to determine the subsidence area caused by coal mining, we took 100 mm as the minimum subsidence threshold caused by coal mining. The area of subsidence value larger than 100 mm is about 180,075 m 2 . The goaf area is 49,789 m 2 . The ratio is 3.6:1. According to Table 4, the area with a subsidence value of 100-300 mm is 102,585 m 2 , which means the majority of subsidence was caused by coal mining. The maximum subsidence value is about 1700 mm. We found that the larger the subsidence value, the smaller the corresponding area, and the area ratios of different subsidence values are shown in Table 4. The isoline map is a classical graphical representation, such as a contour map, which can be used to represent elevation of terrain. In this study, we took the isoline map to express the subsidence value, named the subsidence isoline map, which can clearly show the surface subsidence values in different locations. The subsidence isoline map shows different forms when setting different interval values, as shown in Figure 15.  Table 4 indicates that areas with a subsidence value of less than 100 mm occupied more than 55% of the study area, in which most of the subsidence values are caused by external factors. In order to determine the subsidence area caused by coal mining, we took 100 mm as the minimum subsidence threshold caused by coal mining. The area of subsidence value larger than 100 mm is about 180,075 m 2 . The goaf area is 49,789 m 2 . The ratio is 3.6:1. According to Table 4, the area with a subsidence value of 100−300 mm is 102,585 m 2 , which means the majority of subsidence was caused by coal mining. The maximum subsidence value is about 1700 mm. We found that the larger the subsidence value, the smaller the corresponding area, and the area ratios of different subsidence values are shown in Table 4.
The isoline map is a classical graphical representation, such as a contour map, which can be used to represent elevation of terrain. In this study, we took the isoline map to express the subsidence value, named the subsidence isoline map, which can clearly show the surface subsidence values in different locations. The subsidence isoline map shows different forms when setting different interval values, as shown in Figure 15. It can be seen in Figure 15 that the detailed characteristics of the subsidence isoline map with different interval values are different. The smaller isoline interval values correspond to more detailed subsidence characteristics, but also cause the noise to be more pronounced. Through a comprehensive comparison, we found that the isoline map with It can be seen in Figure 15 that the detailed characteristics of the subsidence isoline map with different interval values are different. The smaller isoline interval values correspond to more detailed subsidence characteristics, but also cause the noise to be more pronounced. Through a comprehensive comparison, we found that the isoline map with a subsidence interval value of 100 mm can retain the details of subsidence characteristics of coal mining and reduce the impact of noise.
After analysis of the subsidence isoline map with an interval value of 100 mm, we took the subsidence value of 100 mm as the minimum subsidence threshold and extracted the subsidence boundary. The result is shown in Figure 16. a subsidence interval value of 100 mm can retain the details of subsidence characteristics of coal mining and reduce the impact of noise.
After analysis of the subsidence isoline map with an interval value of 100 mm, we took the subsidence value of 100 mm as the minimum subsidence threshold and extracted the subsidence boundary. The result is shown in Figure 16.  Figure 16 indicates that the area of subsidence is larger than the goaf area, and the subsidence area is offset to the right relative to the goaf. The subsidence area was 140,271 m 2 during 7 November 2020 and 19 May 2021, which is smaller than the 180,075 m 2 calculated in Table 4, and the ratio of subsidence area to goaf was 2.8:1.
A subsidence value larger than 100 mm is considered to be caused by coal mining. We divided the subsidence values into six grades, and the corresponding area of each grade was calculated, and then expressed in a specific color, as shown in Figure 17.   Figure 16 indicates that the area of subsidence is larger than the goaf area, and the subsidence area is offset to the right relative to the goaf. The subsidence area was 140,271 m 2 during 7 November 2020 and 19 May 2021, which is smaller than the 180,075 m 2 calculated in Table 4, and the ratio of subsidence area to goaf was 2.8:1.
A subsidence value larger than 100 mm is considered to be caused by coal mining. We divided the subsidence values into six grades, and the corresponding area of each grade was calculated, and then expressed in a specific color, as shown in Figure 17. a subsidence interval value of 100 mm can retain the details of subsidence characteristics of coal mining and reduce the impact of noise.
After analysis of the subsidence isoline map with an interval value of 100 mm, we took the subsidence value of 100 mm as the minimum subsidence threshold and extracted the subsidence boundary. The result is shown in Figure 16.  Figure 16 indicates that the area of subsidence is larger than the goaf area, and the subsidence area is offset to the right relative to the goaf. The subsidence area was 140,271 m 2 during 7 November 2020 and 19 May 2021, which is smaller than the 180,075 m 2 calculated in Table 4, and the ratio of subsidence area to goaf was 2.8:1.
A subsidence value larger than 100 mm is considered to be caused by coal mining. We divided the subsidence values into six grades, and the corresponding area of each grade was calculated, and then expressed in a specific color, as shown in Figure 17.  As shown in Figure 17, the areas with different subsidence grades gradually expanded outward with the goaf from the center. The subsidence can be divided into six grades: 100-300 mm, 300-600 mm, 600-900 mm, 900-1200 mm, 1200-1500 mm, and more than 1500 mm. The corresponding subsidence areas are 65,101, 29,083, 16,354, 13,178, 11,123, and 5432 m 2 . The ratio among the subsidence grades is 46:21:12:9:8:4.

Discussion
The accuracy of the DSuM was determined by point cloud and DEM. Therefore, we analyzed the accuracy of the point cloud and DEM, and the results are shown in Tables 2 and 3. We also analyzed the influence of interpolation grid size (resolution of DEM) on the accuracy of the DEM. The results indicate that the grid size of 0.1 m is suitable for the study area (Table 3, Figure 8). Compared to the results of Yu H. et al. [26], the average error of DEM generated by UAV-based LiDAR increased from 0.32 to 0.039 m. Although UAV oblique photogrammetry can achieve centimeter-level accuracy by arranging high-density ground control points [43], it is not applicable in coal mining subsidence monitoring due to the difficulty of setting ground control points in the coal mining area. The proposed method can achieve centimeter accuracy without ground control points. This paper analyzed and compared the accuracy of point cloud, DEM, and DSuM. The results indicate that the accuracy of the DSuM is better than the accuracy of point clouds and DEMs, as shown in Tables 2 and 3 and Figures 7 and 10. The reason is that the DSuM generated by two periods of UAV-based LiDAR data adopts the same flight parameters, and the data processing process is exactly the same. Therefore, it can eliminate some systematic errors and improve the accuracy.
We calculated the strike and incline boundary angles with the threshold of 100 mm, which are 63.8 • and 66.5 • , respectively. The boundary angle value is similar to the strike boundary angle of 61.2 • calculated by GCPs. The result indicates that line analysis could be applied to the calculation of boundary angle with the subsidence threshold of 100 mm, but it is difficult to detect a boundary of 10 mm. In the area analysis, there was a series of subsidence values less than zero. The reasons for this may be external factors, such as rain wash and crop farming. Therefore, to further determine the subsidence area caused by coal mining, we set the subsidence value of 100 mm as the maximum threshold value caused by non-coal mining to eliminate the influences of other factors. We calculated and drew the subsidence isoline map with four different interval values. The result indicates that the subsidence isoline map with an interval of 100 mm is suitable for subsidence monitoring; see Figure 15. It can perfectly express the characteristics of coal mining subsidence. We calculated the area of subsidence area according to the entire study area and the isoline map (Figures 14 and 16 and Table 4). The results show that the area calculated by the isoline map is smaller than the area calculated according to the entire study area, the reason for which is that the area calculated by the entire study area includes some subsidence not caused by coal mining. The final subsidence area caused by coal mining should be the result calculated by the isoline map with interval value of 100 mm.

Conclusions and Future Works
Subsidence detection is important work for coal mining safety. There are various methods used to detect the deformation of mining areas. However, the accuracy is difficult to guarantee in mountainous areas. In this study, UAV-based LiDAR was used to monitor the ground surface subsidence of the working face of a coal mining area, and a new model termed DSuM is proposed to detect the subsidence deformation of a coal mining area. The accuracy of the DSuM was verified by GCPs, and the data mining was performed based on the DSuM to obtain the parameters required for coal mining subsidence monitoring. Subsequently, point, line, and area analyses of the DSuM were conducted. The results indicate that UAV-based LiDAR can be used to monitor continuous changes of the working face in a coal mining area, which provides basic information for subsidence prediction and damage recovery in mining areas.
Some issues are still worth investigating. The accuracy of UAV-based LiDAR can reach the centimeter level, which is not suitable for areas with small subsidence values, especially when the subsidence values are less than 100 mm. On the other hand, the ground surface deformation also includes horizontal movement, but the DSuM cannot express horizontal deformation. Therefore, future work will be focused on further improving the accuracy of the DSuM and creating a model that can represent horizontal deformation.