Automated Extraction of Ground Fissures Due to Coal Mining Subsidence Based on UAV Photogrammetry

Widespread ground fissures caused by coal mining subsidence are a main cause of ecological destruction in coal mining areas, and the rapid monitoring of ground fissures is essential for ecological restoration. Traditional fissure monitoring technologies are time consuming and laborious. Therefore, we developed a method to automatically extract ground fissures from high-resolution UAV images. First, a multiscale Hessian-based enhancement filter was utilized to enhance the ground fissures in grayscale images. Then, a simple single-thresholding operation was applied to segment the enhanced image to generate a binary ground fissure map. Finally, incomplete path opening was performed to eliminate the noises in the fissure extraction results. We selected the N1212 working face of the Ningtiaota Coal Mine in Shenmu County, China, as the study area. The results indicated that the ranges of correctness, completeness, and the kappa coefficient of the extracted results were 66.23–79.00%, 69.03–73.22%, and 67.91–75.88%, respectively. Image resolution is the key factor for successful fissure detection; the method proposed in this paper can extract ground fissures with a width greater than one pixel (2.64 cm), and the detection ratio for fissures with a width greater than two pixels was over 87%. Our research has solved the problem of the rapid monitoring of ground fissures to a certain extent and can act as a valuable tool for ecological restoration in mining areas.


Introduction
Coal is a valuable chemical material and energy source that has long held the second largest share of world energy consumption, with total world production exceeding 7500 Mt in 2020 [1]. Coal mining has created a series of environmental degradation problems. Coal mining can be divided into open-pit and underground mining. Unlike open-pit mining, which directly damages the ground, underground mining causes ground subsidence by destroying the underground rock strata, thus affecting the ecological environment. Largescale ground subsidence has led to the development of numerous ground fissures, especially in ecologically fragile arid areas [2].
Ground fissures created due to mining subsidence cause changes in soil properties, reduce soil quality [3], and decrease surface soil moisture [4] in the region where they are distributed. At the same time, ground fissures also result in mechanical damage to plant roots [5,6], affecting the expected growth of vegetation, leading to vegetation degradation and aggravated soil erosion in coal mining areas [7,8]. Ground fissures may even penetrate the underground mining space, seriously threatening the safety of mine production [9]. Therefore, rapid investigation and the mapping of ground fissures that have been caused by mining subsidence are essential for ecological restoration in mining areas.
Traditionally, ground fissure maps are obtained by field surveys through the use of a total station, GPS, and other means to measure the plane position of ground fissures [10,11]. Due to the dense distribution of ground fissures in mining areas, the field surveys are time consuming and labor intensive. According to actual fissure survey data, the width of ground fissures is usually at the sub-meter level [12][13][14], meaning that the resolution of current optical satellite images cannot clearly display ground fissures. Nowadays, UAV photogrammetry and remote sensing technologies are undergoing rapid development and have been widely used in many fields [15]. UAV photogrammetry is cheap and highly efficient, and it can provide images with a resolution in the centimetric range, making it possible to rapidly monitor ground fissures in mining areas.
Visual interpretation of high-resolution images, a classical method to determine linear feature extraction, cannot meet the demands of massive ground fissure detection; therefore, an automatic extraction algorithm is vital. At present, the combination of UAV photogrammetry and an automatic extraction method is frequently applied for fissure detection research in bridges, roads, and landslides [16][17][18][19][20][21]; however, relatively few studies have explored the application of automatic approaches for the ground fissure mapping in mining areas. Zhang et al. [22] proposed a fissure identification algorithm based on machine learning to detect ground fissures from UAV images of the mining subsidence area in the eolian sand area of western China. Jia et al. [23] successfully detected ground fissures via high-resolution UAV images using the matched filtering method. Zhao et al. [24] innovatively used thermal infrared images from a UAV and an edge detection operator for fissure extraction in a mining subsidence area. There are various types of ground fissures, and different mining areas have various and complex surface backgrounds; however, the existing research either has a small research scope or relies on land cover information, so it is doubtful as to whether these automatic extraction methods are universal.
Based on this context, we propose a method for the automatic extraction of mininginduced ground fissures from high-resolution UAV images. The developed method is based on a combination of multiscale Hessian-based enhancement filters and incomplete path openings and was tested in the mining subsidence area of the Ningtiaota Coal Mine. The obtained results were compared to the reference data produced by experts by means of visual interpretation. The aim of this paper is to propose an automatic ground fissure extraction method that has high generality and that can quickly monitor ground fissures in mining areas and can accurately diagnose ecological damage.

Types of Ground Fissure Observed in Mining Area
Ground fissures result from the coupling of the movement of the overlying strata and the deformation of the topsoil in the mined-out area after coal mining has taken place. The formation process is a complex mechanical one involving rocks and soil. Based on the formation mechanism, ground fissures can be classified into three types: tensile fissures (Figure 1a), sliding fissures (Figure 1b), and collapsing fissures (Figure 1c) [25]. When the tensile stress generated by the ground movement and deformation exceeds the ultimate tensile strength of the soil, it will undergo tensile failure and form a tensile fissure [26]. The movement and deformation of the ground surface cause the slippage of steep slopes in mountainous or valley terrains, and local breakage occurs to form the sliding fissure [27]. Moreover, the collapsing fissure is formed by the shear failure of the overlying rock and soil above the mined-out area, which is caused by breaking the key stratum (hard strata) that controls surface subsidence [28]. This type of fissure has the characteristics of considerable width and a large step drop.

Study Area
The study area is located above the N1212 working face of the Ningtiaota Coal Mine (38 • 58 -39 • 6 N, 110 • 9 -110 • 17 E) in the northern part of the Loess Plateau, and the southeastern edge of the Mu Us Desert, which lies within Shenmu County, Shaanxi Province, China ( Figure 2). The vegetation is sparse in the study area, and the vegetation types mainly comprise shrubs and grasslands. The study area's elevation ranges from 1220 to 1338 m, making it a typical loess gully region with complex geomorphology (Figure 3). The average annual rainfall is about 434.1 mm, which is concentrated in July, August, and September, and the annual average temperature is 8.6 • C.  The Ningtiaota Coal Mine, where the study area is located, spans east-west for 10 km, spans north-south for 19 km, has a total area of 119.77 km 2 , and a coal production capacity of more than 18.00 Mt/a and is categorized as an underground mine. High-intensity underground coal mining has caused a large number of ground fissures, collapses, and landslides on the surface of the study area, and gullies have formed after the surface runoff formed by rainfall, intensifying soil erosion.

Data Acquisition and Preprocessing
A DJI Matrice 210 RTK with a professional-grade digital camera (the Zenmuse X5S) mounted on a three-axis gimbal was used for image acquisition. The main camera parameters are listed in Table 1, and detailed information can be found on the product's official website (DJI. https://www.dji.com/ (accessed on 13 January 2022)). Sensor type CMOS UAV images were acquired on 21 April 2019, images were taken in sunny and breezy weather conditions. The front overlap and side overlap of the UAV missions were 80% and 55%, respectively, and the ground sampling distance (GSD) was 1.8 cm at a flight altitude of 80 m. A total of 377 photos covering the entire area were obtained in JPEG format. The Matrice 210 RTK was able to provide high-precision horizontal position and elevation information [29], and the ground surface of the study area had moved due to coal mining subsidence. Therefore, we did not set ground control points. Pix4D mapper software (Pix4D, Lausanne, Switzerland) was used to process the collected images to generate a digital orthophoto map (DOM) and digital surface model (DSM) of the study area with a resolution of 2.64 cm.

Methods
Although the first general edge detection algorithms were proposed in the 1980s [30,31], linear feature extraction is still a challenging task in many fields, such as medical research [32], geoscience [16,22], and signal processing [33]. For our research, the main challenges affecting automatic fissure detection are as follows: (1) the approach should be able to detect ground fissures of different widths and directions simultaneously; (2) the approach should not respond to edges but should instead enable the detection of dark curvilinear structures; (3) unlike roads, rivers, and blood vessels, ground fissures are not necessarily continuous linear structures, and broken fissures cannot be connected simply by morphological operators. Considering these challenges, we developed a processing workflow that included three main stages to extract ground fissures. First a multiscale Hessian-based enhancement filter [34,35] was used to filter the grayscale images to enhance the ground fissures. Second, a simple single-thresholding operation was used to distinguish ground fissures from the image background and to generate a binary fissure map. Finally, for the incomplete path openings [36,37], a flexible morphological operator was performed to eliminate some of the false positives.

Multiscale Hessian-Based Enhancement Filtering
The multiscale Hessian-based enhancement filter is an efficient vessel enhancement algorithm that is widely used in the field of medical image processing [32,34,38]. Some related works applied this enhancement filter for the detection of fissures in multicrystalline solar cells and pavements [39,40]. In this study, we utilized the multiscale Hessian-based fissure filter to perform fissure enhancement and background suppression operations.
The local intensity of the grayscale image I of the study area can be expressed using the Taylor expansion of the relationship between pixel x and its neighborhood, which has a granularity δx [34]: where x = (x, y) T is the two-dimensional column vector of the pixel position, ∇ x,σ and H x,σ are, respectively, the gradient vector and Hessian matrix at pixel x in the image I at scale σ. The Hessian matrix at scale σ is given as follows: Here, is the second-order partial derivative of the image I. In order to facilitate the differential operation in the image I while also reducing the effect of noise [38], the second derivative of the Gaussian function was used as the convolution kernel to process the image I. The second-order partial derivative of the image I can be described as follows: where γ is the normalization factor, and G σ (·) is the two-dimensional Gaussian function with standard deviation σ, which is written as follows: where x 2 is the square sum of the elements x and y in column vector x. Therefore, the Hessian matrix H x,σ at pixel x in the image I at scale σ is expressed as follows: The two eigenvalues λ 1,σ (x) and λ 2,σ (x) of the Hessian matrix H x,σ and their corresponding eigenvectors can be solved. Considering the gray value, the image I can be regarded as a curved surface. The eigenvalues λ 1,σ (x) and λ 2,σ (x) are the two principal curvatures at the pixel x of this curved surface, and the corresponding eigenvectors are the principal directions. Figure 4 shows the fissure morphology in the grayscale image as well as the schematic diagram of the principal curvatures of the fissure segment. The eigenvalues λ 1,σ (x) and λ 2,σ (x) of the ground fissures have the following properties [34]: In order to highlight the ground fissures in the image I so that the fissures have a sharp contrast with the complex background, our research used the ground fissure measurement function shown below: where β and c are applied to adjust the sensitivity weights of the measures R σ and S σ . Commonly, the value of β is 0.5, and the value of c is related to the overall intensity of the image. The value of v σ (x) ranges from 0 to 1, and the closer it is to 1, the more robust the ground fissure response is. The scale σ of the measurement function v σ (x) is continuously adjusted to adapt to ground fissures of different widths. In this paper, we select 70 cross-sections of ground fissures with different widths and fit their gray curves to calculate σ, which is set from 0.5 to 2.4 according to the actual situation in the study area ( Figure 5). Additionally, the maximum response value in each scale is used as the fissure measurement I e (x) at pixel x to maximize fissure enhancement:

Thresholding
After multiscale Hessian-based enhancement filtering, the ground fissures showed significant differences from the image background. Therefore, a simple single-thresholding operation was used to effectively segment the ground fissures and the image background and to generate a binary ground fissure map. It should be noted that a low threshold will cause many false positives to be present in the binary fissure map, while a high threshold will eliminate narrow-width fissures; so, the actual situation should be considered when determining the global threshold for image segmentation. In this paper, the mean and standard deviation (std) of the intensity of the enhanced image was calculated, and the global thresholds T were, respectively, set to mean, mean + 0.5 std, mean + 1 std, mean + 2 std, and mean + 3 std. The fissure detection results ( Figure 6) under different global thresholds revealed that there were a lot of noises in the binary ground fissure maps corresponding to the threshold mean and mean + 0.5 std (Figure 6b,c). Although the extraction effect of the threshold mean + 1 std was improved, there were still noises that were mixed and connected with fissures (Figure 6d), and the threshold mean + 3 std underestimated the fissure pixels, resulting in the fissures continuously breaking into short fragments (Figure 6f). The threshold mean + 2 std achieved a good balance between ground fissure extraction and noise suppression ( Figure 6e). Therefore, in our research, mean + 2 std was selected as the global threshold to segment the images that were enhanced by multiscale Hessian-based enhancement filtering.

Incomplete Path Openings
Due to visually similar objects, such as elongated shadows caused by micro-topography, the binary ground fissure map generated by the described routine still comprises numerous types of noise, which are similar in shape to the short fragments of ground fissures. Normally, noises are distributed in a disorderly manner, with the small ground fissure fragments closely arranged and having a clear tendency to extend (Figure 7). However, due to the different widths and directions of ground fissures, it is not easy to find a structural element that is suitable for the whole situation to perform morphological processing on the image. Therefore, a flexible morphological method for a path opening operation [36,37] was used to effectively suppress noise while retaining the ground fissures. When performing path opening operations, we need to build a directed graph that is based on binary ground fissure images. E is a set of points that represents all of the pixel locations in a binary image. A directed graph is defined based on these points through a binary adjacency relation ' →' (Figure 7), and ' x → y ' means that there is an edge connecting point x to point y. In this case, point y denotes the successor of point x, and point x is the predecessor of pixel y. These concepts are illustrated in Figure 7. The successor set of each point is expressed as follows: If X is a subset of set E, then the successor set of points in X can be written as follows: If there is an adjacency relationship a k −→ a k+1 , then the δ-path with length L can be defined as an L-tuple a = (a 1 , a 2 , a 3 , . . . , a L ), where a k is the element of the δ-path. The set of all δ-paths of length L is denoted as Π L . Given a path a in E, this paper denotes by σ(a) the set of its elements: σ(a) = {a 1 , a 2 , a 3 , . . . , a L } The set of δ-paths of length L in subset X of E is given as Π L (X): All δ-paths of length L in subset X are recorded as α L (X): Moreover, α L (X) represents the path openings along the length L of subset X. When path openings are used for deionizing, if a path lacks a few pixels, the entire path is divided into two paths so that the length of the divided paths cannot reach the set threshold L, eliminating all the pixels comprising the entire path. To reduce the sensitivity of the path morphology to noise, an improved incomplete path opening algorithm [37] was proposed. This algorithm allows for some of the pixels along the path to be missing, which greatly enhances the robustness of the path openings. The incomplete path openings eliminate the residual noise in the fissure extraction results by setting thresholds of the path length, L, and the number of admissible missing pixels, K p .
The images as they appear after the incomplete path opening process with different L, K p are shown in Figure 8; they show that as the path length, L, increases, the noise is better suppressed but that an excessive L leads to the false elimination of short fissures. In addition, K p is also a key factor during the incomplete path opening process. A smaller K p led to the loss of some short fissures, while a larger K p resulted in the images having a considerable amount of noise. Therefore, proper L and K p can achieve balance between fissure preservation and noise suppression. In this article, L and K p were set to 50 and 3, respectively, and were determined according to the actual situation in our study area.

Validation
The centerlines that were visible in the ground fissure maps were visually interpreted and were used as reference data. Then, the centerline of the automatically extracted result was obtained. Drawing on the previous linear element extraction research [20], the confusion matrix (Table 2) between the reference data and the extraction result was obtained to quantify the effects of the ground fissure extraction method used in this paper. According to the corresponding relationship between the fissures in the extraction data and the reference data, the following four items were calculated in the confusion matrix: (1) true positive (TP), (2) true negative (TN), (3) false positive (FP), and (4) false negative (FN). There is usually a small amount of deviation between the fissure centerline in the reference data and the extraction data in the actual implementation process. According to the actual width of the fissures, a buffer with a distance of 5 pixels (13.2 cm) was set on both sides of the fissures in reference data, and the parts of the extraction data within the buffer were considered to match the manual reference [20]. The confusion matrix was used to calculate the user's accuracy, producer's accuracy, and Cohen's kappa coefficient (κ) to evaluate the ground fissure extraction accuracy of our extraction method. In some related works, the user's accuracy and producer's accuracy are also called correctness and completeness, respectively [41,42]. The correctness, completeness, and κ can be expressed as follows: where r is the number of rows in the confusion matrix, X ii is the element in row i and column i on the major diagonal of the matrix, and X i+ and X +i denote the sum of the elements in row i and column i, respectively.

Accuracy Assessment
Due to the enormous scope of the study area, this paper selected four typical test sites that were 40 m × 40 m in size ( Figure 9) and mapped the centerline of the ground fissures by visual interpretation as the ground truth (reference data). The fissure extraction results and ground truth were used to establish a confusion matrix, and the correctness, completeness, and kappa coefficient were calculated as the evaluation indexes to determine the extraction accuracy. The extraction results for the ground fissures in the four selected sites are shown in Figure 10. In contrast with the original images, we found that some of the fissures with smaller widths were not successfully extracted, and correspondingly, the producer's accuracies (completeness) of the four selected areas were between 66.23% and 79% (Table 3). Meanwhile, some elongated shadows that were formed by micro-topography or hightexture ground were identified as fissures, leading to user's accuracy (correctness) values of 69.03% to 73.22% at the four typical sites (Table 3). According to the existing research, a kappa coefficient greater than 75% represents substantial consistency [20], and the kappa coefficients of the four sites were 69.11%, 75.88%, 68.82%, and 67.91%, respectively (Table 3). Therefore, the ground fissure extraction method proposed in this paper showed good performance when tested in the study area.

Impacts of L, K p
Although the incomplete path opening operation performed well for denoising, the shapes of the ground fissures in the study area were very complex, so a single L, K p cannot provide a comprehensive indication of fissure extraction performance. Therefore, we set L = (30,35,40,45,50,55,60,65,70) and K p = (1, 2, 3, 4, 5, 6, 7) to generate 63 ground fissure maps at each site and calculated the completeness, correctness, and kappa coefficient to evaluate the impacts of L, K p on fissure extraction. The completeness, correctness, and kappa coefficient corresponding to different L, K p at different sites are shown in Figure 11. The results indicated that the completeness decreased as the L decreased and increased as the K p increased, while the correctness changed conversely. The combination of a smaller L and larger K p made it easier to obtain higher completeness but resulted in lower correctness. Although there were differences in the completeness and correctness among the four sites, their changing trends were roughly the same. The reasons for this are simple. If a smaller L is set, the incomplete path openings would contain more fissure segments that were shorter in length, and a larger K p would allow this process to detect paths with sufficient length in more fissure segments, thus retaining the most fissure pixels but also greatly increasing the presence of non-fissure pixels in the extraction results, leading to lower correctness. Therefore, the L, K p combinations with higher kappa coefficients were mainly located near the diagonal of the coordinate system, but large L, K p combinations reduced the kappa coefficients ( Figure 11). According to the impact analysis of L, K p , 50 and 3, as employed in Section 4.3, was the combination with high kappa coefficients for all four sites, quantitatively demonstrating that this combination could achieve good ground fissure extraction results in the study area. It is worth noting that L and K p can be set flexibly according to fissure extraction demands or the need to eliminate noise while maintaining the key fissures.

Ground Fissure Extraction Result
The fissure extraction results in the study area ( Figure 12) revealed that when underground mining began in the area, fissures were generated on the ground at the starting position of the working face. As coal mining operations continued to advance, numerous C-shaped fissures were generated inside the working face and spread forward along the advancing direction, which is consistent with the results of field surveys and a theoretical analysis of the impacts of rock strata movement on fissure development [12,43]. However, the C-shape of the fissures inside the working face was incomplete, and loess gullies interrupted the fissure development process. The extracted fissures were mainly distributed in relatively flat areas, and there were few fissures in the steep slopes of the gullies (Figure 12). In addition, some areas with dense fissures were located at the top of both sides of gullies, and through field investigations, we found that these fissures were formed by landslides on steep slopes that had been disturbed by underground mining (Figure 12a,b). However, the densely distributed fissures in the extraction results of this study were not all real ground fissures. Figure 12c,d shows that noise may be caused by the erroneous extraction of the thin ground erosion gullies and the collapsed soil fragment gaps, which are very similar in shape to real ground fissures.

Discussion
Ground fissure is one of the main disasters caused by underground coal mining; so, it is necessary to develop rapid ground fissure monitoring and extraction technology. In this paper, we have successfully extracted ground fissures from UAV high-resolution images by using an automatic detection algorithm. The results of fissure extraction in the study area indicated that the proposed method can effectively avoid the edge extraction of nonfissure features under complex surface backgrounds. Meanwhile, the ground fissures with different widths and directions in the study area were extracted by using the multiscale Hessian-based enhancement filter. Furthermore, broken fissures cannot be connected simply by morphological operators, which leads to their morphological similarity with noise, but the utilization of incomplete path openings enables our method to distinguish short fissure segments and noise well. In dealing with the challenges in the extraction of ground fissures in mining area, the proposed method has achieved ideal fissure extraction accuracy (Table 3).
Automatic extraction of ground fissures in mining areas is a complicated process, and traditional edge detection operators and image segmentation methods can hardly meet the requirements. To prove this, the Canny operator [30] and object-oriented segmentation were used to extract ground fissures at the four selected sites. The object-oriented segmentation was handled in eCognition software (Trimble Germany Gmbh, Munich, Germany), using an appropriate scale parameter to segment images and then using a support vector machine (SVM) to extract ground fissures. The Canny operator is only sensitive to the edge information in the image, and unilateral edges of partial fissures were extracted, as shown in the area marked by shape A in Figure 13a. In addition, the extraction results of the Canny operator contained numerous nonfissure edges (e.g., the areas marked by shapes B in Figure 13b and C in Figure 13c), and we may not be able to distinguish which fissure the edge belongs to because of the close distance between fissures (e.g., the area marked by shape D in Figure 13c). Therefore, it is also difficult for us to evaluate the extraction accuracy of Canny operator. Although the extraction results of object-oriented segmentation can roughly reflect the shape of fissures (Figure 13e-h), this method mistakenly classified shadows in images as fissures (e.g., the area marked by shape C in Figure 13c). The kappa coefficients of the extracted results of object-oriented segmentation in the four selected sites were 48.84%, 45.77%, 54.74%, and 48.25%, respectively, which were lower than those for our method. The proposed method shows better fissure extraction performance than the Canny operator and object-oriented segmentation (e.g., the ground fissures marked by shape E in Figure 13l).
According to the fissure extraction completeness, less fluctuation was observed (66.23-79.00%) in the ratios of the real fissures at the four typical test sites, mainly because the images with the same resolutions were acquired at the same time and under the same weather conditions. However, the highest detected ratio of 79.00% also revealed that the image resolution (2.64 cm) used in this paper was still unable meet the needs of automatic small fissure extraction. The ground fissures that were detected by the automatic extraction algorithm in this study were typical linear features. Previous research has shown that the image resolution used for extraction is an important factor that affects the accuracy of linear feature extraction [44][45][46]. We randomly selected 360 locations in the study area that had ground fissures and verified whether the fissures in these locations were successfully detected in the extraction results. The results indicated that the detection ratio was 72.22%, which was close to the completeness of the four typical test sites. As shown in Figure 14, fissures can be divided into five categories based on width. All fissures with a width within one pixel were not detected, and the detected ratio of fissures with a width from one pixel to two pixels was only 45.52%. As the fissure width increased, the detection ratio increased rapidly. The detection ratio exceeded 87% when the fissure width was greater than two pixels, and all fissures were detected as long as their width was greater than four pixels. The method proposed in this paper cannot detect ground fissures with a width less than one pixel; however, it can detect fissures with widths greater than two pixels stably. In a study on automatic fissure extraction in the Super-Sauze mudslide, Stumpf et al. [16] proposed that lower resolution UAV images do not necessarily correspond to lower fissure extraction accuracy but that the spatial resolution of the image is at least equivalent to the target fissure width. Therefore, in practical applications, we should fully consider the width of ground fissures in mining areas when setting the UAV flight altitude to obtain orthoimages with the required resolution. The correctness of the fissure extraction executed at the four test areas ranged from 70.08% to 82.70%, indicating that noise is contained in the extraction results. Due to the inherent mechanism of our method proposed in this paper, the extraction results have some collapsed soil fragment gaps (Figure 15a), thin surface erosion gullies (Figure 15d), and vegetation texture features (Figure 15g), resulting in a decrease in the fissure extraction accuracy (Figure 15b,e,h). The shapes of these noise sources are very similar to the ground fissures, and their cross-sectional grayscale curves are consistent with the fissures (Figure 15c,f,i). At present, the algorithm proposed in this paper cannot distinguish these noise sources from ground fissures. Meanwhile, the low solar altitude angle leads to a reduction in the fissure extraction accuracy [16]. Our method actually uses the shadows generated by deep ground fissures to realize automatic extraction. As shown in Figure 16, the noise sources are not as deep as the ground fissures, and when the solar altitude angle is large, the shadow that affects the extraction accuracy cannot be formed. However, in the case of a lower solar altitude angle, the three noise sources are more likely to produce shadows that are falsely extracted as ground fissures. Our UAV images were collected at 16:00 on 21 April 2019, and the corresponding solar altitude angle was 26 • 45 , which reduced the fissure extraction accuracy. Furthermore, in this study, full sunlight at the surface may decrease the extraction accuracy further [16]. Although there was a great deal of noise in the images, the fissure extraction results still reflect the fissure distribution in the study area. However, the distribution characteristics of the extracted fissures are not completely consistent with underground mining theories because of the large slope areas. Through theoretical analysis and numerical simulation, existing research on slope stability has revealed that the slope is simultaneously affected by its own gravity and the additional mining stress that is imparted upon it during ground movement and deformation processes [47,48]. Disturbed by underground coal mining, slippage occurs in areas with large slopes, accompanied by the formation of sliding fissures [47]. Meanwhile, the generation of ground fissures on the slope increases horizontal displacement along the downslope direction [27], further deteriorating the stability of the slope and leading to landslides. A large number of landslides in the large slope area destroyed the fissures, and only the portion of the fissure in the top area of the slope is retained, resulting in fewer fissures being detected in large slope areas in the extraction results.
After the long-term improvement and innovation of image processing and extraction algorithms, the accuracy of the automatic extraction of linear elements has continued to improve. In some studies on the automatic extraction of road fissures, the correctness and completeness of the extraction results were both more than 70% [49,50]. Al-Rawabdeh et al. [20] obtained extraction results of 82.24% and 75.0% for correctness and completeness for the automatic extraction of landslide features. The fissure extraction accuracy in this paper may not reach the same level as the one obtained in other linear element extraction studies; however, it still solves the main problem related to the automatic extraction of ground fissures in coal mining subsidence to a certain extent. Meanwhile, the proposed method was also compared with the existing research on ground fissure extraction. Zhao et al. [24] used classical edge detection operators to extract fissures from UAV thermal infrared images, and we analyzed some defects of this type of operators in ground fissure extraction. The automatic extraction algorithm proposed by Jia et al. [23] has a fixed-scale σ; so, its adaptiveness to fissures with different widths needs further research. Moreover, denoising according to the fissure candidate area would lead to the false removal of some short fissure segments. Zhang et al. [22] divided the UAV image into four datasets according to the characteristics of background information: bright ground, dark ground, withered vegetation, and green vegetation, and then extracted ground fissures by machine learning. However, in general, the types of features in mining areas are more than the above categories. Therefore, it is noted that our proposed method has certain advantages in extracting ground fissures from complex background information.

Conclusions
In this paper, an image processing chain was developed to extract ground fissures through high-resolution UAV images, and it was tested on the working face of the Ningtiaota Coal Mine in order to realize the automatic extraction of ground fissures caused by coal mining subsidence. In our method, the multiscale Hessian-based enhancement filter was used to enhance the ground fissure information, and a simple single-thresholding operation was then employed to generate a binary ground fissure map; the incomplete path openings method was applied to eliminate the noise in the fissure extraction results. The fissure extraction results were displayed with the fissure centerlines.
Verified by manual visual interpretation, the completeness, correctness, and kappa coefficient of the fissure extraction results at the four selected test sites were 66.23-79.00%, 69.03-73.22%, and 67.91-75.88%, respectively. The path length and the number of admissible missing pixels in the incomplete path openings had a significant influence on the accuracy of the extraction results, and the L, K p of (50, 3) that was used in this study was found to be reasonable. Our research indicates that image resolution is the key factor for successful fissure detection, and the main sources of fissure extraction noise are collapsed soil fragments, thin erosion gullies, and vegetation textures. In addition, the solar altitude angle and illumination conditions also influence extraction accuracy. The fissure distribution in the study area was simultaneously affected by underground mining and topography.
This study verifies the feasibility of using high-resolution UAV images to automatically extract ground fissures caused by coal mining subsidence. Our method will greatly improve the monitoring efficiency of ground fissures and provide more valuable data for ecological restoration in mining areas.

Data Availability Statement:
The data presented in this study are available from the corresponding author upon request.