Investigation of Unmanned Aerial Vehicles-Based Photogrammetry for Large Mine Subsidence Monitoring

: With the exploitation of underground sources, nature receives a huge negative impact on the local environment introducing surface subsidence. A mining region needs to be observed in sequences before, during, and after coal extraction from the coal mine. Di ﬀ erent measuring methods exist to monitor subsidence, and all of them apply various instrumentation. A choice of methodology depends on access to a ﬁeld of observation and requested accuracy. Obviously, the most accurate results provide geometric leveling, but, many times, the terrain does not allow surveyors to walk over the dangerous outﬁelds. Looking for the most adequate and feasible method, this research did a comparison between observation of the same points, applying statistical analysis of di ﬀ erences between the reference points heights, and tested methods. Monitoring procedure comprised utilization of total stations (TS), global navigation satellite system (GNSS), and unmanned aerial vehicles (UAV). In this paper, the Velenje coal mine was taken as a case study, and observation data were collected during 2017.


Introduction
To avoid surface subsidence circumstances, it is, therefore, recommended before the excavation begins to make an important geological and geo-mechanical analysis of the region where a future mine will be opened. Investigations, which must be carried out, are studies of lithology, structure, stratigraphy, morphology, seismology, and hydrology of the area of interest. To assure minimal impact on the environment from mining, the methodology of mining, excavation instrumentation, and displacement monitoring planning must be analyzed and applied in advance [1]. As a part of this preparatory work, surveying procedures have to be introduced to better monitor and control behavior of subsidence.
The mining site Velenje, located in the NE region of central Slovenia in Šalek Valley, is packed with a high quantity of Pliocene sediments. Figure 1 introduces the geological profile of the Velenje basin, laying over rich lignite seam placed between Gorenje-Šoštanj block and Southern Karavanke mountain. A cross-section shows three aquifers layers Pl-1, Pl-2, and Pl-3, the Quaternary, the Pliocene-Quaternary, and the Pliocene, respectively [2]. Likewise, lithologies Triassic (T1, T2) and Paleozoic (Pz) are shown on the same profile [3]. Due to the inception of Velenje colliery, groundwater has been playing a dangerous role in mining safety. Pliocene layers are made up of diverse materials, such as marlstone, laminated mudstone, sand, and leans of gravel. At the same time, those layers correspond to natural protection penetration. Mining mine jeopardizes the safety, producing cracking of protective aquifers layers and giving a chance to adulterate lignite seam. During this process, the hydrological nature constantly changes. As groundwater tends to seep into the excavated area, the pumping wells systems are installed. Water pressure is making changes on Pl-1 and Pl-2, which are visible from the ground above [2]. Playing with natural consequences of dewatering may have important destructive effects on the regional environment. In the Šalek valley, drilling began to be explored in 1873. During this and the following year, several wells were drilled along the edge of the valley, which revealed several meters of poor-quality coal. Finally, in 1875, a 37.6-m-thick lignite layer was drilled at a depth of 101 m. The year 1875 marks the starting year of exploitation of the Velenje colliery, 145 years since a rich coal deposit was discovered by drilling over the lignite seam [4].
Today, Velenje colliery is a leader in Europe by its dimensions. The coal seam is 8.3 km long, 2.5 km wide of thickness up to 160 m, with an average thickness of 60 m, and approximately 400 m deep [5]. Several excavation methods have been used since its first exploitation: room, pillar, and block caving method. In the early 1950s, the longwall mining method was introduced; it has undergone various improvements and, today, is known as the Velenje mining method (VMM) [6]. Excavation of lignite by time formed over the mining site pull-apart basin, where actually a high subsidence rate occurs [2].
The VMM methodology, which is unique, has working principles based on excavation of coal seam layer by layer. There are two important sections of excavation, one introduced as a foot line and another as a hanging wall. The foot line section is 4-5 m high, where hydraulic shield support machinery is used to protect the lower excavation section from the upper one, which is well schematically shown in Figure 2. The hanging wall section represents a layer between two excavation neighboring sections, in the height of 5 to 17 m [5]. This section is exposed to permanent stress, crouching the coal. VMM cannot permit any coal pillar to sustain overload soil stress, and that is why subsidence can occur without any enounce.
The subsidence does not appear suddenly; its development takes time and occurs when the progress of excavation moves forward. Therefore, continuous monitoring is necessary. Velenje coal mine was chosen as a case study for this paper, with data collected in the period from January 2017 to January 2018. Figure 3 shows a 3D elevation map with positions of monitoring points G1 to G5, over the mining site. In the Šalek valley, drilling began to be explored in 1873. During this and the following year, several wells were drilled along the edge of the valley, which revealed several meters of poor-quality coal. Finally, in 1875, a 37.6-m-thick lignite layer was drilled at a depth of 101 m. The year 1875 marks the starting year of exploitation of the Velenje colliery, 145 years since a rich coal deposit was discovered by drilling over the lignite seam [4].
Today, Velenje colliery is a leader in Europe by its dimensions. The coal seam is 8.3 km long, 2.5 km wide of thickness up to 160 m, with an average thickness of 60 m, and approximately 400 m deep [5]. Several excavation methods have been used since its first exploitation: room, pillar, and block caving method. In the early 1950s, the longwall mining method was introduced; it has undergone various improvements and, today, is known as the Velenje mining method (VMM) [6]. Excavation of lignite by time formed over the mining site pull-apart basin, where actually a high subsidence rate occurs [2].
The VMM methodology, which is unique, has working principles based on excavation of coal seam layer by layer. There are two important sections of excavation, one introduced as a foot line and another as a hanging wall. The foot line section is 4-5 m high, where hydraulic shield support machinery is used to protect the lower excavation section from the upper one, which is well schematically shown in Figure 2. The hanging wall section represents a layer between two excavation neighboring sections, in the height of 5 to 17 m [5]. This section is exposed to permanent stress, crouching the coal. VMM cannot permit any coal pillar to sustain overload soil stress, and that is why subsidence can occur without any enounce.
The subsidence does not appear suddenly; its development takes time and occurs when the progress of excavation moves forward. Therefore, continuous monitoring is necessary. Velenje coal mine was chosen as a case study for this paper, with data collected in the period from January 2017 to January 2018. Figure 3 shows a 3D elevation map with positions of monitoring points G1 to G5, over the mining site.  The Active Subsidence Remediation Area (ASRA) is located between two lakes-Družmirje and Velenje lake. While the height difference between lakes is approximately 7 m, and the whole area is subjected to subsidence that amounts to 5-10 m/year, the ASRA area is constantly upgraded with earth works to maintain the desired altitude of the area between lakes [8]. Lake Družmirje is 30 m deeper than Velenje lake ( Figure 4). Depths of lakes are measured by sonar reason sound 110 in combination with the global navigation satellite system (GNSS) since 2010. This is also one of the solutions to how subsidence could be measured when the area of observation starts to be covered by water.  The Active Subsidence Remediation Area (ASRA) is located between two lakes-Družmirje and Velenje lake. While the height difference between lakes is approximately 7 m, and the whole area is subjected to subsidence that amounts to 5-10 m/year, the ASRA area is constantly upgraded with earth works to maintain the desired altitude of the area between lakes [8]. Lake Družmirje is 30 m deeper than Velenje lake ( Figure 4). Depths of lakes are measured by sonar reason sound 110 in combination with the global navigation satellite system (GNSS) since 2010. This is also one of the solutions to how subsidence could be measured when the area of observation starts to be covered by water. The Active Subsidence Remediation Area (ASRA) is located between two lakes-Družmirje and Velenje lake. While the height difference between lakes is approximately 7 m, and the whole area is subjected to subsidence that amounts to 5-10 m/year, the ASRA area is constantly upgraded with earth works to maintain the desired altitude of the area between lakes [8]. Lake Družmirje is 30 m deeper than Velenje lake ( Figure 4). Depths of lakes are measured by sonar reason sound 110 in combination with the global navigation satellite system (GNSS) since 2010. This is also one of the solutions to how subsidence could be measured when the area of observation starts to be covered by water.

Summary of Previous Research
There were several operational applied modes in monitoring mine subsidence for the longwall excavation method of mining. For instance, in Jakarta (Indonesia), a combination of three methods was used together: leveling, GNSS, and interferometric synthetic aperture radar (InSAR). Observations have been running for a period of 28 years, yielding to the conclusion that temporal and spatial variations of measurements play the main role in understanding land subsidence behavior [9].
Another comparison was done in Taiwan [10], where an integrated model of 'absolute' and "relative" modes showed a significant impact on the dependability of subsidence. However, using two sets of GNSS solutions with the consistent differential mode of geoidal information between the benchmarks provided root mean square (RMS) difference of 5 cm.
A case study of deep mining subsidence [11] provided an analysis of the combination of the probability-integral method with differential interferometric synthetic aperture radar (D-InSAR) technique. It was concluded that using the mining subsidence prediction model with D-InSAR could improve the accuracy of the quantitative estimation of deep mining subsidence rates and time threshold [12].
In China, for the Nantun mining area, the subsidence down below the Zouji highway was introduced by a new methodology named probability integration model small baseline set (PIM-SBAS). This methodology is made with a combination of probability integration model (PIM), small baseline set (SBAS), and terraSAR-X (TSX) satellite images. Unfortunately, time-series interferometric synthetic aperture radar (TS-InSAR), where SBAS and TCPInSAR (temporarily coherent point InSAR) belong, did not reach expected results due to its capability to monitoring only linear subsidence where priority is given to deformation velocity. Comparison with classical leveling

Summary of Previous Research
There were several operational applied modes in monitoring mine subsidence for the longwall excavation method of mining. For instance, in Jakarta (Indonesia), a combination of three methods was used together: leveling, GNSS, and interferometric synthetic aperture radar (InSAR). Observations have been running for a period of 28 years, yielding to the conclusion that temporal and spatial variations of measurements play the main role in understanding land subsidence behavior [9].
Another comparison was done in Taiwan [10], where an integrated model of 'absolute' and "relative" modes showed a significant impact on the dependability of subsidence. However, using two sets of GNSS solutions with the consistent differential mode of geoidal information between the benchmarks provided root mean square (RMS) difference of 5 cm.
A case study of deep mining subsidence [11] provided an analysis of the combination of the probability-integral method with differential interferometric synthetic aperture radar (D-InSAR) technique. It was concluded that using the mining subsidence prediction model with D-InSAR could improve the accuracy of the quantitative estimation of deep mining subsidence rates and time threshold [12].
In China, for the Nantun mining area, the subsidence down below the Zouji highway was introduced by a new methodology named probability integration model small baseline set (PIM-SBAS). This methodology is made with a combination of probability integration model (PIM), small baseline set (SBAS), and terraSAR-X (TSX) satellite images. Unfortunately, time-series interferometric synthetic aperture radar (TS-InSAR), where SBAS and TCPInSAR (temporarily coherent point InSAR) belong, did not reach expected results due to its capability to monitoring only linear subsidence where priority Minerals 2020, 10, 196 5 of 14 is given to deformation velocity. Comparison with classical leveling measurements showed that PIM-SBAS had better results than simple SBAS and TCPInSAR methods [13].
Coal mining subsidence was possible to monitor as well by gravimetry instead of classical leveling surveys [14]. The methodology was based on a negative linear correlation between gravity and height changes.
Airborne laser scanning (ALS) surveys were used for subsidence mapping in an Australian mine, where protentional benefits of this method were proven in the case study given in [15].
Terrestrial laser scanner (TLS) has become an all-purpose tool for applications needing huge datasets. One of such applications is mine subsidence monitoring. The research on the utilization of TLS in mining applications is presented in [16].
In Italy [17] and in New South Wales, Australia [18], subsidence was observed in three different ways, using unmanned aerial vehicles (UAV), total station, and laser imaging detection and ranging (LIDAR) measurements. Final results showed similar precision between the methods, but from experience, it was concluded that high-resolution images taken by UAV and processed by digital surface models (DSM), in comparison with two other techniques, were much more practical and faster.
For the Campine coal basin, Belgium, movement of the surface above the closed coal longwalls mine was studied for a period of 18 years using radar interferometry satellite data. The analysis was based on observation of the movement of downward and uplift of the coal basin. This effect was noticed in the majority of the closed coal mines and, in cases, where pumping water in the mine was stopped. The valuable conclusion said that even 40-60 years after closure, subsidence was still observable and occurring, of course, with a small rate. Using only satellite data for this type of analysis was absolutely adopted, but it might be supplemented by observing a larger area to better understand hydrological behavior, especially for flooding of underground zones [19].
To estimate active subsidence over the coal mine Velenje, [20] proposed a method where photogrammetric measurements obtained by an unmanned aerial vehicle (UAV) were introduced through a 3D point cloud. The observation area was divided into rectangular sectors, where also centroids of their planes were calculated. A model developed by the Faculty of Natural Sciences and Engineering (FNSE), in Ljubljana, was used to estimate the consolidation time, followed by cloud-to-cloud (C2C) analysis to identify the surface subsidence. This prognosis model gave only the prediction of subsidence from geometric and volume loss perspective.

Methods
As introduced previously through different examples, monitoring methodologies for any mining site have their own observation strategy based on environmental circumstances, availability of instrumental resources, accessibility to satellite data, the capability to manipulate by UAV, etc. Photogrammetry produces elevation models by UAV, which is a promising technology, due to its low costs regarding instrumentation and resources, as well as a huge amount of measured data. In the case study presented in this research, the overall accuracy of observing subsidence over Velenje coal mine using UAV was tested against reference methods, namely GNSS real-time kinematic (RTK) and total station (TS) tachymetry.
Our observation material consisted of measurements taken during 2017, where the test of congruence of UAV and traditional techniques were performed on five monitoring points (G1 to G5), located at the longwall site -80C over the Velenje coal mine. The heights of the monitoring points were determined by GNSS RTK and TS tachymetry. The monitoring points and the map of the terrain over the mine, as depicted in Figure 5, were located over the central mining area. Such a decision was made since this is the area where the biggest displacements are expected: by computerized model [21], where predicted final maximal vertical displacement on the surface above longwall -80C is approximately 10 m. All data were obtained from [22]. Depending on the weather conditions, availability of access to the points, and the quality and acceptance of the measurements themselves, we identified a set of valid measurements collected during 2017. The core of the research was to test the standard deviation of five control points' heights obtained by UAV against those obtained by GNSS (and TS, where possible). UAV is the most promising method for massive data collection, especially in areas with limited accessibility. It is capable of creating high definition digital elevation models (DEMs), allowing various analyses of topography and artificial objects. Limiting issues for UAV photogrammetry, besides camera characteristics, are weather conditions (for example, strong winds or snow), errors in models, or wrong identification of surface points, just to mention some of the sources. Therefore, we chose certain five, well-established points, that could be easily identified in images. Those points' coordinates were determined by GNSS and/or TS during each period when UAV campaigns were performed. The congruence of UAVobtained heights with the ones determined by GNSS and TS was tested on height differences between the method under test (UAV photogrammetry) and the reference ones (GNSS RTK and TS tachymetry). Using that approach, the subsidence itself, as well as possible differences in coordinate system definitions, became negligible because only height differences were statistically tested on equality. It means that this research did not deal with subsidence predictions, but rather only with testing the values of the control points' heights for particular measuring epochs. Therefore, there were no trends to incorporate in the statistical model. The choice of preferring GNSS RTK and TS tachymetry over geometric leveling (as the most accurate surveying method for determining heights) was a result of extremely hard environmental conditions, including wet and soft terrain with limited accessibility, where geometric leveling did not give acceptable results. There were few geometric leveling campaigns performed, but loop closures of leveling polygons, in all campaigns, did not fulfill quality assessment requirements. Therefore, this method was excluded from further investigation.
UAV measurements were performed using DJI Phantom 4 with DJI FC330 camera, in resolution 4000 px × 3000 px, focal length 3.61 mm, sensor size 6.2 mm × 4.6 mm, and pixel size 1.542 µm. The collected images were processed with 3D Survey aerial image processing software. The size of the obtained digital orthophoto (DOF) image was approximately 677 m × 743 m, with a pixel resolution of 4 cm.
Images are processed in 3DSurvey Mapping and Aerial Image Processing Software [23]. The number of ground control points (GCPs) used for fitting DOF to the resulting coordinate system was 13-17, depending on the particular campaign. Fitting of DOF image to GCPs was around 5 cm for each measuring campaign, which was the value we adopted for the standard deviation of UAVdetermined points' heights, considering, also, flight plans, disposition, and density of ground control points. The bundle adjustment of the DOF models was performed following, so-called, global mode [23], which allows for automatic checking of the produced model. The core of the research was to test the standard deviation of five control points' heights obtained by UAV against those obtained by GNSS (and TS, where possible). UAV is the most promising method for massive data collection, especially in areas with limited accessibility. It is capable of creating high definition digital elevation models (DEMs), allowing various analyses of topography and artificial objects. Limiting issues for UAV photogrammetry, besides camera characteristics, are weather conditions (for example, strong winds or snow), errors in models, or wrong identification of surface points, just to mention some of the sources. Therefore, we chose certain five, well-established points, that could be easily identified in images. Those points' coordinates were determined by GNSS and/or TS during each period when UAV campaigns were performed. The congruence of UAV-obtained heights with the ones determined by GNSS and TS was tested on height differences between the method under test (UAV photogrammetry) and the reference ones (GNSS RTK and TS tachymetry). Using that approach, the subsidence itself, as well as possible differences in coordinate system definitions, became negligible because only height differences were statistically tested on equality. It means that this research did not deal with subsidence predictions, but rather only with testing the values of the control points' heights for particular measuring epochs. Therefore, there were no trends to incorporate in the statistical model. The choice of preferring GNSS RTK and TS tachymetry over geometric leveling (as the most accurate surveying method for determining heights) was a result of extremely hard environmental conditions, including wet and soft terrain with limited accessibility, where geometric leveling did not give acceptable results. There were few geometric leveling campaigns performed, but loop closures of leveling polygons, in all campaigns, did not fulfill quality assessment requirements. Therefore, this method was excluded from further investigation.
UAV measurements were performed using DJI Phantom 4 with DJI FC330 camera, in resolution 4000 px × 3000 px, focal length 3.61 mm, sensor size 6.2 mm × 4.6 mm, and pixel size 1.542 µm. The collected images were processed with 3D Survey aerial image processing software. The size of the obtained digital orthophoto (DOF) image was approximately 677 m × 743 m, with a pixel resolution of 4 cm.
Images are processed in 3DSurvey Mapping and Aerial Image Processing Software [23]. The number of ground control points (GCPs) used for fitting DOF to the resulting coordinate system was 13-17, depending on the particular campaign. Fitting of DOF image to GCPs was around 5 cm for each measuring campaign, which was the value we adopted for the standard deviation of UAV-determined points' heights, considering, also, flight plans, disposition, and density of ground control points. The bundle adjustment of the DOF models was performed following, so-called, global mode [23], which allows for automatic checking of the produced model. Figure 6 contains the DOFs for three selected UAV campaigns (1 February, 28 June 2017, and 15 January 2018). The images showed a decrease in access opportunities to the control points, which were located in the middle of the area. The third image in row (right), from Figure 6, showed that the points G4 and G5 were covered by water, that is why data were not provided. Also, in some epochs, some of the points were under mud, which made the measurements unreliable or even impossible (for example, G3). Nevertheless, points' disposition was carefully planned and selected, having in mind the safety of the field crew, as well as expected displacements.
value was completely sufficient for subsidence monitoring since the value of subsidence was 10 m per year. If we assumed that the maximum distance from the network point to a control point G-1 to G-5 was 500 m with the maximum vertical angle of ±10°, and estimated error of instrument's height measurement of 1 mm, following the error propagation law, the resulting standard deviation of the observed point's height would be: with: • -standard deviation of a control point G-i (i = 1, 2, …, 5), • -standard deviation of a polygonometry point j, • ∆ -standard deviation of the height difference between P-j and C-i, and • -standard deviation of measured height of TS over the polygonometry point.
Since standard deviations of the measured height difference and height of the instrument were negligible to the standard deviation of the polygonometry point, the final value for the control point's height surveyed by TS tachymetry was 5 cm, which was the value used further in this research. GNSS observations were collected by LEICA Viva GS10 GNSS receiver, with a declared horizontal standard deviation of 8 mm + 1 ppm and a vertical standard deviation of 15 mm + 1 ppm, both in the RTK mode. The measuring sessions lasted 30 s per point. Having in mind measurement uncertainty of the method itself and transformation of the obtained coordinates to the state coordinate system, the standard deviation of GNSS RTK surveyed points' heights was estimated to be 5 cm.
Also, three control height measurements were made by TS tachymetry. The measurements were performed by LEICA TCRA 1201+ (Heerbrugg, Switzerland) with a standard deviation of 1" for angles and 1 mm + 1.5 ppm for distances with a single reflective prism, which was the configuration utilized in this research. Although qualitative performance for this type of TS suggested measurements of sub-centimeter precision, it had to be realized that the area where the measurements were taken was wide and rich in topography. The mine control network consisted of a three-level network: 1.
primary mine network, placed out of the influence zone, 2.
measuring points network (for accessing detail), and 3.
connecting network, which connected the networks 1 and 2 The final standard deviation of the points' position in network 2 was estimated to be 5 cm. This value was completely sufficient for subsidence monitoring since the value of subsidence was 10 m per year. If we assumed that the maximum distance from the network point to a control point G-1 to G-5 was 500 m with the maximum vertical angle of ±10 • , and estimated error of instrument's height measurement of 1 mm, following the error propagation law, the resulting standard deviation of the observed point's height would be: with: • σ 2 Since standard deviations of the measured height difference and height of the instrument were negligible to the standard deviation of the polygonometry point, the final value for the control point's height surveyed by TS tachymetry was 5 cm, which was the value used further in this research. Figure 7 depicts the time series holding heights of the points G-1 to G-5 obtained by GNSS RTK surveying. Certain gaps in measurements were noticed. They came from the inappropriate measuring conditions when the collected GNSS data did not meet the quality assurance requirements. It was shown that the appearance of subsidence was moving differently from point to point, from 2.5 m to 5.9 m.  Figure 7 depicts the time series holding heights of the points G-1 to G-5 obtained by GNSS RTK surveying. Certain gaps in measurements were noticed. They came from the inappropriate measuring conditions when the collected GNSS data did not meet the quality assurance requirements. It was shown that the appearance of subsidence was moving differently from point to point, from 2.5 m to 5.9 m. Congruence between GNSS and UAV results was tested by calculating the differences between the heights obtained by these two methods, i.e.,

Results
with: • -the height of the control point i obtained by GNSS RTK, • -the height of the control point i obtained by UAV photogrammetry, and • i = 1, 2, ...,5 is the index of the control point.
Differences in the control points' heights obtained from UAV and GNSS, following Equation (2), are given in Table 1. Locations of the control points were chosen to be over the central mine's area because the largest displacements are expected there. We processed 13 measuring campaigns in total, from 1 February to 10 October 2017. The dates of UAV flights (and, consequently, GNSS measurements) are given in column 2. Columns 3 to 7 contain differences (Equation (2)) for each control point. It could be noticed from Table 1 that the total number of the calculated differences per point varied between measuring epochs.
The main reason for inconsistency in the total number of differences (Equation (2)) per epoch was related to limited accessibility to the central area above the mine. There were more UAV campaigns during the second half of 2017, but they were not relevant to this research because almost Congruence between GNSS and UAV results was tested by calculating the differences between the heights obtained by these two methods, i.e., with: • H GNSS i -the height of the control point i obtained by GNSS RTK, • H UAV i -the height of the control point i obtained by UAV photogrammetry, and • i = 1, 2, . . . ,5 is the index of the control point.
Differences in the control points' heights obtained from UAV and GNSS, following Equation (2), are given in Table 1. Locations of the control points were chosen to be over the central mine's area because the largest displacements are expected there. We processed 13 measuring campaigns in total, Minerals 2020, 10, 196 9 of 14 from 1 February to 10 October 2017. The dates of UAV flights (and, consequently, GNSS measurements) are given in column 2. Columns 3 to 7 contain differences (Equation (2)) for each control point. It could be noticed from Table 1 that the total number of the calculated differences per point varied between measuring epochs. Table 1. Differences between GNSS and UAV-obtained heights of selected control points. The main reason for inconsistency in the total number of differences (Equation (2)) per epoch was related to limited accessibility to the central area above the mine. There were more UAV campaigns during the second half of 2017, but they were not relevant to this research because almost all of the control points were underwater or in the mud, so their identification was impossible (e.g., Figure 6, well depicts the loss of point G4 covered by water). In Figure 8, the situation directly from the observation field is presented. In the upper area of the photography, a large subsidence step crack was clearly visible, as well as the appearance of water overflow. The onset of mud and sludge could be distinguished, too.
Minerals 2020, 10, x FOR PEER REVIEW 9 of 14 all of the control points were underwater or in the mud, so their identification was impossible (e.g., Figure 6, well depicts the loss of point G4 covered by water). In Figure 8, the situation directly from the observation field is presented. In the upper area of the photography, a large subsidence step crack was clearly visible, as well as the appearance of water overflow. The onset of mud and sludge could be distinguished, too. It could be, even, seen from Table 1 that the number of usable observations of control points' heights and, consequently, the number of results used in this research decreased through the year. The N/A cells mean that no measurement on the particular point was made, while the cells shaded in grey mark the measurements that were unreliable, either due to poor identification on DOF images It could be, even, seen from Table 1 that the number of usable observations of control points' heights and, consequently, the number of results used in this research decreased through the year. The N/A cells mean that no measurement on the particular point was made, while the cells shaded in grey mark the measurements that were unreliable, either due to poor identification on DOF images or due to missed GNSS measurement. Each row in Table 1 represents one height difference (Equation (2)) per epoch. During the observation period, some parts of the test field were flooded, when GNSS measurement for an affected point could not be taken. Also, some of the points could not be approached in each epoch due to deep mud. Since the safety of the surveying crew was considered as the primary factor, some measurements had to be skipped. Therefore, the number of calculated differences (Equation (2)) corresponded to the last row of Table 1. However, since we did not test or predict the subsidence trend, the missed cells (i.e., differences) did not affect the generality of this article's conclusions. This statement means that only differences between heights obtained in the same epoch by two methodologies (UAV photogrammetry and GNSS RTK) were tested here.
The basic statistical indicators of the obtained results are presented in Table 2. For applied n m (number of statistically approved measurements), where m recalls to the particular point number, basic central tendencies were calculated: arithmetical means ∆ i , empirical standard deviations σ ∆ i , and widths of intervals: with i = 1, 2, . . . , 5. To prove the absence of outliers in the dataset from Table 1, an outliers test for measurement series was performed, using the test statistics: where σ 2 (1) is a dispersion of the series with the biggest standard deviation (in this case study, σ 2 5 ) and: with f (2) = f 1 + f 2 + f 3 + f 4 degrees of freedom. Individual degrees of freedom for particular series were f i = n i − 1, where n i = 1, 2, . . . , 5 is the number of measurements for the series i. After applying the numerical results, the test statistics (Equation (4)) gave the following conclusion: which confirmed that all outliers were correctly removed from the observation material.
To explore congruence between standard deviations σ i , i = 1, 2, . . . , 5, the Bartlett test was applied. Therefore, the null hypothesis: was tested against the alternative: H 1 : at least one σ i is different from the others.
The test statistics were: with: Table 2, Introducing numerical results into Equation (8), Bartlett test showed: After approving homogeneity of dispersions, the final step in the analysis was to test the average differences of heights obtained by two tested methods. The following test statistics were used: Standard deviation σ ∆ was developed from predicted accuracies of GNSS and UAV methods. After applying the error propagation law to (2), and considering its two components as not correlated, the standard deviation of the difference ∆ was then: Adopting the standard deviations σ GNSS = σ UAV = 50 mm, it yielded to σ ∆ = 70.7 mm. Knowing σ ∆ and introducing values of ∆ from Table 2 into Equation (10), the results presented in Table 3 were obtained. Table 3. Test statistics for mean differences of heights obtained by GNSS and UAV.

Discussion
All five test statistics from Table 3 fell below N(0, 1) = 1.96, which implied a conclusion that both methods gave congruent results and could be used mutually for monitoring mine subsidence.
To test compliance of UAV against two other methods, a total of three control height measurements were performed throughout the testing period. The results are depicted in Table 4. To test the congruence of GNSS and UAV datums with TS, we applied again the statistics (Equation (10)), but, this time, introducing the values from columns 6 and 7 from Table 4. If adopting 20 mm for the standard deviation of TS tachymetry obtained point heights, the standard deviation for the statistics was then: because we assumed earlier σ GNSS = σ UAV = 50 mm. After applying Equations (10) and (12) to data from columns 6 and 7 from Table 4, it yielded the results presented in Table 5. Again, all test statistics from Table 5 were lower than N(0, 1) = 1.96, which proved the congruence of testing methods with the referent one.
Considering the discussion on features of using UAV, GNSS, and TS for mine subsidence monitoring, one could summarize the pros and cons for each of the three methods, as given in Table 6. TS tachymetry gives the best standard deviation performance. However, it is not economic due to its high costs by means of time, human, and technical resources. Also, preparation works (the control network, for example) and rather slow and hard fieldwork, with sparse datasets that could be obtained, make this method least acceptable for utilization. The situation with GNSS is slightly better. The standard deviation of RTK mode is comparable with TS tachymetry; it needs less time in the field, as well as in the office. But, again, the amount of collected measurements is often too small for purposes of creating high-definition DEMs. On the other hand, UAVs can provide huge datasets of measurements, with a short time spent in the field. Processing of collected photos is performed using sophisticated software, and a variety of reports are available. Still, the assistance of an experienced expert is needed in all project phases, from preparing the flight plan to monitoring the flight and processing obtained photos. The standard deviation of UAV photogrammetry is lower than GNSS, but it can be adopted with correct flight height, methods of orientation, disposition of GCPs, and adjustment parameters.

Conclusions
Continuous investigation of improvements in monitoring mine subsidence is done in order to increase security and decrease costs. Therefore, various new methods and instrumentation are tested for their suitability for this purpose. Due to problems of approaching the characteristic monitoring points, weather conditions, overall accuracy, and productivity of the monitoring process, search for new methodologies are directed to massive data collection methods.
Due to the continuous increase of wet area: water coming from both lakes-Velenje and Družmirje-and the appearance of mud and sludge, monitoring of subsidence by humans is not possible for all points in a given time. That is why UAV, and GNSS to some extent, are promising methodologies/instrumentation for rapid and massive data collection. However, the usage of UAVs in monitoring is often limited due to accuracy requests. This case study was made in order to empirically and statistically approve using UAV for mine subsidence monitoring, especially in a situation when the observation field is covered by water and when human access is limited to work on site.
Five monitoring points in mining site Velenje were observed during 2017 using GNSS RTK, in order to test the congruence of UAV-obtained heights. Some extra control measurements obtained by TS tachymetry were performed, in the areas and timeframes when the right conditions were met.
Two important conclusions were drawn from our research. First, the standard deviation of UAV was congruent with GNSS, which was proven by adopting the hypothesis of equality of means (Equation (10) and Table 3). The absence of systematic effects was proven by control measurements using TS tachymetry, confirming the congruence of GNSS and UAV datums with the one used by TS. Statistical analysis of differences between the reference (TS tachymetry and GNSS RTK) and tested method (UAV photogrammetry), given in Table 5, was approved using UAV photogrammetry for monitoring mine subsidence.
The overall conclusion suggests that both GNSS RTK and UAV photogrammetry are suitable for mine subsidence monitoring. A significant advantage is given to UAV because of much larger sets of data available after measurements and, also, side products, such as high-density DEM, which can be used for other activities during and after mine excavation, likewise, to avoid human presence over the risk area of observation.
Author Contributions: D.I.S., data processing, theory, and modeling, J.R., performing measurements and providing data, M.V., theory, model definition, and development of field methods. All authors have read and agreed to the published version of the manuscript.