MONITORING TERRAIN DEFORMATIONS CAUSED BY UNDERGROUND MINING USING UAV DATA

Underground mining causes terrain surface deformations that lead to various threats to the environment and people, thus a systematic deformation monitoring needs to be performed. This monitoring mainly focuses only on the vertical part of the deformation and remote sensing techniques are currently very often used for this purpose. The development of Unmanned Aerial Systems (UASs) open new possibilities in this context. Most commonly, the mapping UASs are equipped with RGB cameras but also other lightweight sensors are utilized. In this work, the usefulness of UAS photogrammetry and LiDAR data is investigated in the context of detection and measurement of terrain deformations caused by underground mining. The accuracy of the methods was compared in reference to TLS data. The UAS and TLS measurements were performed in 2018 and 2019 but the subsidence was also evaluated in regards to ALS data acquired in 2011. The standard methodology based on Digital Terrain Models of Difference (DoDs) was applied to detect the subsidence. The DoD analysis was restricted to the hard surfaces. The profiles along the roads were also analysed to validate the accuracy of the data. The analysis showed that the UAS photogrammetry enables to obtain less noisy data and more accurate results of the terrain subsidence measurement than the UAS LiDAR sensors. The comparison of the DoDs showed about 33 cm subsidence between 2011 and 2018, which gives a subsidence rate of about 5 cm/year. The observed subsidence between years 2018 and 2019 was equal to about 5 to 15 cm depending on the measurement technique and investigated area.


INTRODUCTION
Underground mining causes various threats to the environment and people. One of such threats is deformation of the terrain surface, which is especially dangerous in urbanised areas as it can result in severe infrastructure damages. For this reason, it is necessary to conduct systematic monitoring of the terrain surface deformations in the vicinity of the mine.
The terrain deformation can be measured using different methods and equipment, such as inclinometers, typical surveying techniques, and remote sensing techniques, including utilization of Interferometric Synthetic Aperture Radar (InSAR), Light Detection and Ranging (LiDAR) and photogrammetric data. Each of these methods has its own advantages and limitations. The traditional surveying methods, such as levelling or total station measurements are characterized by high accuracy of the results. However, their application in terrain deformation monitoring is limited by their small spatial resolution and time-consuming field work that is necessary to acquire the data. In contrast, the remote sensing techniques allow for fast acquisition of relatively accurate data with high spatial and temporal resolution.
Recently, the InSAR technique is commonly applied to terrain deformation measurements (e.g., Pawłuszek-Filipiak and Borkowski, 2020;Zhu et al., 2020). InSAR allows for determination of vertical, relatively small deformations, with high accuracy, high time resolution, and even with no data acquisition costs. However, its application is limited to nonvegetated areas and the results have low spatial resolution in comparison to other remote sensing techniques. In contrast, airborne LiDAR data has higher spatial resolution and is suitable even for forest areas, but its vertical accuracy is much worse. Moreover, execution of flights is expensive, which results in a very limited time resolution of the performed measurements. The monitoring costs can be reduced through the use of Unmanned Aerial Systems (UASs). Most of the mapping UASs are equipped with RGB cameras, and some of them were already used in determining terrain deformations for small and non-vegetated areas. Unmanned Aerial Vehicle (UAV) photogrammetry is commonly used in many applications related to terrain deformation analyses, such as monitoring of landslides (e.g., Balek and Blahůt, 2017;Lucieer et al., 2014), glaciers (e.g., Dall'Asta et al., 2017;Kraaijenbrink et al., 2016), or open-pit mines (e.g., Esposito et al., 2017). UAV photogrammetry is less frequently used to monitor the impact of underground mining on its adjacent areas. However, this problem was investigated by e.g., Puniach et al. (2021). Recently, also light-weight LiDAR sensors utilization for deformation monitoring using UASs have been investigated (Zieher et al., 2019). However, none of the above publications analysed the application of UAV-borne Laser Scanning (ULS) data to deformation monitoring in the underground mining area. Since deformations caused by underground mining occur on larger areas and result usually in terrain subsidence, the goal of this study is to assess the potential of both types of UAS data: images and LiDAR in the monitoring of terrain vertical deformations caused by underground mining.

TEST SITE AND DATA ACQUISITION
The test site is located in the Upper Silesia region (Poland) that is affected by many closed and active coal mines. Since this region is highly populated and urbanised, even small deformations may cause severe material loses. The test site is a suburb area (~ 0.25 km 2 ) placed in the proximity of an active underground coal mine (Figure 1). Investigated area was not subjected earlier for detailed monitoring with respect to terrain deformations, however, some geodetic surveying (mainly geometrical levelling) conducted in this region mentions progressive subsidence. Selected area covers the coal seam, thus larger deformations are expected after the extraction. The test UAS data was collected during three campaigns that took place in March 2018, September 2018 and June 2019. During the first campaign, the ULS data was collected using Velodyne HDL-32E laser scanner, and images were collected with Digital Single Lens Reflex (DSLR) camera. The same camera was used to collect data in the second campaign in 2018, however, flights with laser scanner were not executed in this campaign. During the third campaign, the ULS and photogrammetric data was collected using GreenValley LiAir 50 system. During both campaigns, in which ULS data was collected, the Global Navigation Satellite System (GNSS) static data from a local base station was acquired for LiDAR data georeferencing. Ground Control Points (GCPs) for image data georeferencing were collected using GNSS real-time network. The same method was used for collecting check points. In addition, Terrestrial Laser Scanning (TLS) data was collected along the roads in July 2018 and July 2019 and Airborne Laser Scanning (ALS) point cloud collected in 2011 was acquired from the national database. The airborne LiDAR data was used as the first measurement to determine the amount of the deformation while TLS data was used to validate UAS data.

UAS photogrammetric data
UAS photogrammetric data in 2018 was collected with Nikon D800 DSLR camera equipped with 24 mm focal length lens and mounted on Aibot X6 V2 hexacopter ( Figure 2). In both UAS campaigns executed in 2018, images were collected along the same 4 parallel lines with planned image overlap equal to 90% and 45% along and across flight lines, respectively. The Ground Sampling Distance (GSD) was planned to 15 mm. The acquisition of images was executed in separate 7 flights due to several reasons. One of them was the relatively short flying time of the Aibot X6 V2 with the Nikon D800 camera. Note that maximal take-off weight of this drone is around 6.5 kg while the weight of used camera and lens is around 1.5 kg. Second reason were large differences of terrain heights (about 40 m, see Figure 1) that caused the need of performing flights at different absolute altitudes to keep overlaps and GSD close to planned values. This drone keeps constant altitude when it performs flight automatically from one to another waypoint even if these waypoints have different heights. The total number of collected images used in subsequent processing was equal to 689 and 551 for the spring and autumn campaign, respectively. Lower number of images in the second campaign was caused by smaller overlap between flights and performing two flight plans with one flight. All images were geotagged with the position determined by the drone on-board GNSS receiver.
In 2019 UAS photogrammetric data was collected with Sony a600 camera equipped with 16 mm focal length lens (24 mm equivalent of full frame sensor). Although the field of view of used cameras and lenses were similar, parameters of the flight plans in 2019 were different since the flight plan in that year was prepared to collect images and LiDAR data during the same flights. The UAS photogrammetric data was collected in 2019 for a slightly larger area in two separate flights, one with 7, second with 9 lines. Different number of lines was due to changes in terrain height. Similarly as for Aibot X6 V2, the drone used in 2019 also performs flight at constant altitude when it performs flight along a planned line. Different number of lines caused larger sidelap (~60%), but the endlap was slightly lower and equal to about 85%. The planned GSD was almost 20 mm. The total number of collected images was equal to 908.
During all campaigns, GCPs and ground check points were measured with GNSS real-time network technique that allowed to determine coordinates with 3 and 5 cm accuracy for planar and vertical components, respectively. Most of the measured points were natural and only a few were marked in the field to keep appropriate distribution. The total number of measured points in both campaigns executed in 2018 was equal to 20 and 40 GCPs and check points, respectively. However, one of measured GCPs in the second campaign was excluded from the processing due to its poor accuracy caused by the occlusion with tree foliage. In the third campaign (2019), the number of measured natural points was increased resulting in 24 and 56 GCPs and check points, respectively. In addition, 6 GCPs were created from ULS data as a geometrical centre of a planar patch fitted to the point cloud representing a small flat object, e.g., picnic table in the garden.
The photogrammetric data was processed using Agisoft Metashape (formerly Photoscan) Professional software. The accuracies of image block adjustments are shown in Table 1. Horizontal accuracies are similar for all campaigns, however, the difference in vertical accuracy is even 2.5 cm. The reason for this could be probably actual image overlap since the accuracy corresponds to the number of adjusted images.

Campaign
No of images After the image block bundle adjustments, the dense point clouds were created. Next, the point clouds were classified to extract ground points and finally, Digital Terrain Models (DTM) of 2.5 cm resolution was created. In addition also orthomosaics were created from images collected in all three campaigns.

ULS data
UAS LiDAR data was collected using mid-shelf UAV-borne LiDAR systems based on Velodyne laser scanners. In the case of the campaign executed in 2018, the UAS was equipped with the Velodyne HDL-32E laser scanner, Sensonor STIM300 Inertial Measurement Unit (IMU) and NovAtel OEM615 GNSS receiver with one antenna. The same, as in the case of RGB camera, UAV platform was used to carry this LiDAR system ( Figure 3). The GreenValley LiAir 50 LiDAR system was used to collect data in 2019. This system has the Velodyne VLP-16 laser scanner, KVH 1750 IMU and NovAtel GNSS receiver with two antennas. Both laser scanners have similar range, though HDL-32 has two times more laser diodes and larger vertical (in normal position) field of view than VLP-16. According to the specification, the navigational sensors used in the LiAir 50 system have better performance than those used in the Aibot X6 V2 system. This can be explained by dual-antenna GNSS sensor the IMU which has a Fiber Optic Gyro (FOG) instead of Micro Electro-Mechanical System (MEMS) gyroscopes. UAS LiDAR data in 2018 was collected only for the northern and central part of the investigated area according to the same flight plans as in the case of flights with RGB camera, but at 10-15 m lower flying altitude. The ULS data was georeferenced according to typical workflow. After georeferencing, the point cloud was coloured ( Figure 4a) with colours of the nearest points form dense point cloud created from UAS photogrammetric data. The average UAS LiDAR point cloud density in the central part of the investigated area was equal to 800 pts./m 2 . ULS data collected in 2019 was georeferenced according to the same workflow and was also coloured with collected RGB data, however, the vendor provided software used to colour the point cloud seems to use orthomosaic for this purpose (Figure 4b). The average point cloud density for data collected in 2019 was significantly lower and equal to 150 pts./m 2 .
ULS data was subsequently processed to extract ground points and create DTM of 5 cm GRID size. Similarly as for dense point clouds, also Agisoft software for point cloud filtering and DTM creation was used.

ALS and TLS data
The Airborne LiDAR data collected in 2011 was obtained from the national database. This data was collected with the density of 12 pts./m 2 and was classified and coloured with RGB aerial images ( Figure 5). This point cloud was used to create a DTM of 25 cm GRID size ( Figure 1) which was used as the reference surface for determining the amount of deformation. The TLS data was collected along the asphalt roads ( Figure 6) in July 2018 and July 2019 using Leica ScanStation C10 scanner. The data collected from separate stations was registered using targets placed in the field between neighbouring stations. In addition, during point cloud registration, also cloud to cloud constraints were added. These constrains were created between point clouds from all neighbouring stations which had sufficient overlap. After the registration the point clouds were georeferenced using points measured in the field with GNSS real time network technique. Detailed parameters of TLS data registration and georeferencing are shown in Table 2. The TLS point clouds were also classified to extract ground points, and the DTMs of 5 cm GRID size were created for roads and used to compare with UAS data.

METHODS
The goal of the investigation is to achieve the information about the deformation detected using UAS photogrammetric and LiDAR data. To better understand obtained deformations, the analysis of the data used and intermediate processing results should be helpful. Consequently, the analysis was executed in three steps: analysis of the data (point cloud), analysis of the DTMs, analysis of differential models which are basis for the determination of vertical deformations.

Point cloud analysis
In this step, the point clouds were investigated in terms of their differences and possible impacts to automatically created DTMs. The 65 meters long vertical cross-sections were created to validate the characteristics of the point clouds obtained by different techniques. The location of the cross section was selected to show point clouds for different objects (e.g., building, road, various vegetation). The cross-sections were created for classified data -for ALS data the fully classified point cloud was used, whereas for other data sets the point cloud was classified into three classes: ground, noise and other. The analysis of the prepared cross-sections will allow for the comparison of the noise level of the point cloud, evaluation of its completeness, and evaluation of the vegetation influence on the results.
The point cloud noise level was also investigated based on the internal error of the point clouds. The internal error was calculated as Root Mean Square Error (RMSE) of fitting planes to two selected point clouds representing planar surfaces. The first surface was a flat roof of the building. To avoid the edge effect and possible errors and noises, only the points from the central part of the roof area of about 13 m 2 were selected. The second surface was a part of a flat, concrete driveway with an area of about 22 m 2 . For each surface a plane was fitted using the least squares method and the RMSE of the distances between the points and plane were calculated.

DTM analysis
Since the investigation aims at the vertical deformations, the approach employing DTMs and differential DTMs was applied. The point cloud densities allowed to create DTMs with the GRID size of 5 cm and 2.5 cm for the UAS LiDAR and photogrammetric data, respectively. However, because of possible noises in the data, the high resolution DTMs may be very rough and negatively impact differential models. In addition, the deformations occur on larger area, thus they may be investigated using lower resolution DTMs. For that reason, DTMs of two different resolutions were created. In particular, DTMs of 5 cm GRID size were created from all UAS data and also TLS data. These models were subsequently resampled to the 0.5 m GRID size. To avoid the influence of the vegetation (missing ground points) during interpolation of DTM, the analysis of the model roughness was executed along terrain profiles in the well modelled areas, i.e., paved roads.

Analysis of differential models
The DTM of Difference (DoD) is a typical tool used for determining size of vertical deformations. However, besides the deformations caused by e.g., underground mining, the DoDs show also terrain height changes caused by other reasons, such as: • Intentional changes in the terrain surfaces (e.g., constructions, earthworks, etc.). • The data and its processing methods (e.g., data noise, points of low vegetation treated as ground points in photogrammetric models, interpolation and filtering methods). • The lack of points under the buildings and high vegetation.
• The accuracy of the acquired data, including the accuracy of georeferencing data.
Therefore, the reliable analysis of DoD in terms of deformation analysis can be performed only on not changed intentionally hard surfaces for which the DTM was created from reliable ground points. In the investigated area such surfaces are asphalt or concrete roads and driveways. Since the orthomosaics were available for all UAS campaigns as well as for ALS data acquisition, hard surfaces were identified in a pair of orthomosaics and a common polygon was extracted to evaluate differential models. Before this evaluation, the polygon was verified if the surface material shown in both orthomosaics is the same. Then, the existing DoD was cut to present only the verified surfaces.
To further analyse the subsidence, the cross-sections of DoD were created along the selected hard surfaces. Then, they were analysed in the terms of subsidence value depending on the location and the capabilities of the collected data to evaluate the deformations.

Point cloud analysis
The performed investigation of the created cross-sections ( Figure 7) enables us to draw the conclusions about the completeness and internal accuracy of the point clouds acquired using different measurement techniques. Since the TLS data was acquired along the roads, no points are expected behind obstacles or in vegetated areas near the roads. First, the analysis of the created cross-sections shows that all data sets obtained in summer are characterized by the lack of ground points under the vegetation. This is true for both ULS and ALS data and is caused by dense vegetation cover that prevents the laser beam from reaching the ground surface. In contrast, the measurements performed before growing season allowed us to measure points under the vegetation. However, it is not clear if the acquired points represent the terrain surface or maybe low vegetation that was left from the previous growing season (e.g., dry grass). What is more, in the case of ULS data collected in 2019, the points belonging to road, building walls and roof are missing. This is likely caused by the small albedo of these surfaces and higher flight altitude in 2019 than in 2018. Second, the internal error of the ULS point clouds is relatively high in comparison to the data acquired by other techniques. The high internal error manifests itself in a noticeably greater thickness of the cross-section lines. For instance, the building walls are thicker than in the case of other data sets, though ULS point clouds are much sparser than point clouds created by dense image matching.
To further analyse the internal accuracy of the point clouds, the roughness parameter was evaluated. The roughness was calculated as the RMSE of fitting planes to two selected point clouds representing planar surfaces. The resulting RMSE values describe the internal accuracy of the point cloud (Table 3).  Table 3. Roughness of point cloud measured as the RMSE of fitting plane to point cloud for part of the roof (13 m 2 ), and driveway (22 m 2 ).
The results of the experiments show that the calculated RMSE value varies from 9.2 to 77.1 mm in the case of driveway and from 9.5 to 64 mm in the case of roof depending on the data collection technique. The best accuracy was achieved by TLS data -the RMSE value was smaller than 1 cm. A slightly bigger RMSE value was obtained for UAS photogrammetric data. The resulting RMSE value was smaller than 2 cm, thus, it is considered satisfactory. The photogrammetric data collected in 2018 is characterized by better accuracy than data collected in 2019, because of the better geometric characteristics of the camera used. The lowest internal accuracy was obtained for ULS data -the RMSE value was equal to 5 -8 cm, with a lower accuracy for data collected in 2018. This can be explained by the worse IMU accuracy of the system used in 2018. The achieved ULS data accuracy confirms the results obtained in Jóźków et al. (2017).

DTM analysis
The problem of DTM interpolation can be observed on DTM cross-sections ( Figure 8). Moreover, during the DTM creation, the influence of the interpolation procedure is also enhanced by a problem of selection of appropriate resolution. High resolution of the DTM leads to high roughness of the cross-sections, especially in the case of ULS data. What is more, all the noise, which was not eliminated in the filtering procedure is clearly visible. Therefore, the models were up-sampled to 50 cm resolution using median filter. As a result, the smoothed DTMs were created. The presented cross-sections show high compatibility of the ULS models with the reference TLS data and the same shape of the DTM created from UAS photogrammetric data, though it is shifted of about 5 cm. The reason for this shift is likely the distortion of the photogrammetric model. The roughness of created 5 cm DTMs is visible the best on a horizontal hard surface, but such surface occurs only in the edge of the investigated area without GCPs. Thus, there were no constraints that could attract the photogrammetric model to the ground. The variance between models in other places may be different. More detailed discussion of the difference between models and the deformation is given in the next section.

Analysis of differential models
The DoD shows vertical changes of terrain height and is often used to determine vertical deformations. However, as mentioned earlier, beside vertical deformations, the DoD shows also changes in terrain height caused by other than deformation reasons, such as methods of data processing or the data inaccuracy. This is also visible in the DoD created for the investigated area ( Figure 9). This model shows the subsidence in the whole area, though some of the areas show terrain uplift. The reason for this is the presence of human made changes caused by earthworks (in south), or by lack of ground points in the area covered by trees resulting in wrong DTM interpolation from photogrammetric data (in north). The reliable analysis of DoD in terms of deformation analysis can be performed only on not changed intentionally hard surfaces for which the DTM was created from reliable ground points. The example of DoD for the identified hard surfaces is shown in Figure 10. This model shows that the strongest subsidence of about 40-50 cm occurred in the central part of the area. The average subsidence of the investigated area in the period 2011-2018 was evaluated to about 33 cm, that gives about 5 cm subsidence each year.
The subsidence was also investigated along 6 cross-sections created on paved roads. Figure 11 shows the cross-sections of 50 cm DoDs created for the main road crossing the investigated area from west to east. The analysis shows that the height difference between years 2011-2018 and 2011-2019 varies depending on the location. The largest differences between different data type collected in the same or similar period occur typically on the edges of the profile. This can be explained by the distortion of photogrammetric models in the edges of the area due to lack of GCPs. Note that all DoDs created using photogrammetric data are distorted in the edges of the profile similarly with respect to DoDs created from using TLS data. The smoothest DoDs were achieved by TLS data, whereas the noisiest results are presented by ULS data. This effect is visible both in 2018 and 2019. Moreover, ULS data collected in 2018 is visibly shifted from the other data sets. The shift in the same direction is present in all of the created cross-sections though its value varies between the cross-sections. This is probably a combination of the ULS data noise (low internal accuracy) and a constant shift which could be an issue related to GNSS data processing. The ULS 2018 data represents smaller section because of the issues with GNSS and IMU on-board data that could not be integrated to obtain trajectory. The differential models created using UAS data collected in 2019 are locally distorted between 200-250 m of the cross-section, because of its partial coverage by trees, and consequently lower number of ground points and less accurate interpolation of the DTM.
The cross-sections of differential models calculated for the data acquired using the same technique in two subsequent years are presented in the Figure 12. TLS data shows that during one year the subsidence along this profile was about 61 mm. The photogrammetric data shows larger subsidence equal to 80 mm on average. However, the larger subsidence can be justified by the fact that the photogrammetric data cover longer time period than the TLS data by 3 months. The subsidence along other profiles varies. The largest subsidence between 2018 and 2019 occurred along asphalt road leading from the north to the centre of the investigated area. This subsidence was on average equal to 15 cm. Unfortunately, because of the shift in ULS data acquired in 2018, the calculated differences are not reliable. Figure 11. Example cross-sections of differential models created with respects to ALS data for the main road crossing the scene from west to east.

Figure 12.
Example cross-sections of differential models created for the same data type for the main road crossing the scene from.

CONCLUSIONS
The conducted investigation showed that UAS data collected with inexpensive sensors can be used for the monitoring of terrain vertical deformations caused by underground mining.
The executed experiments showed that it is possible to detect medium scale deformations using UAV photogrammetry. However, the investigations need to be restricted to the nonvegetated areas. In addition, the data should not be collected during growing season to keep its best quality. The ULS data occurred to be less accurate, thus sufficient to monitor higher magnitude of the deformations. The investigation showed that the main issue with the ULS data is high noise present in the point cloud. This is likely the result of the performance of used sensors, especially IMU. The reduction of the IMU weigh cause the need to use MEMS technology IMUs which performance is still insufficient. On the other hand, the use of multicopters which flights are very unstable in comparison to airplanes, require to use much better IMUs. Probably the quality of ULS data could be also affected by some noise caused by Velodyne laser scanners. However, the sensor development allows to believe that the performance of relatively inexpensive lightweight scanners and IMUs will be increased. The higher accuracy of photogrammetric data was achieved thanks to large number of GCPs, however, UAS LiDAR technique does not require extensive field work to place markers and measure coordinates of GCPs.
The performed analysis showed that the smallest roughness of the data was achieved by TLS technique, whereas the largest roughness was detected for ULS data. In summer, none of the techniques, including ALS, enabled us to detect the points under the vegetation. However, outside of the growing season, the use of ULS allowed us to detect more points under the vegetation than the UAS photogrammetry. The analysis of differential models enabled us to detect the subsidence of 33 cm between years 2011 and 2018, which results in subsidence equal to about 5 cm per year. Moreover, the performed experiments allowed for detection of the subsidence of about 5 to 15 cm between the years 2018 and 2019 depending on the investigated area. The smaller subsidence values obtained by analysis of photogrammetric data can be ambiguous, because of the method accuracy.