[go: up one dir, main page]
Next Article in Journal
Accumulation of Potentially Toxic Metals in Ryegrass (Lolium perenne, L.) and Other Components of Lawn Vegetation in Variously Contaminated Sites of Urban Areas
Previous Article in Journal
Environmental and Economic Analysis of the Production of Oregano Oil Microparticles
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Sustainability
Volume 16
Issue 18
10.3390/su16188039
sustainability-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Deformation-Adapted Spatial Domain Filtering Algorithm for UAV Mining Subsidence Monitoring

by
Jianfeng Zha
1,*,
Penglong Miao
1,
Hukai Ling
1,
Minghui Yu
1,
Bo Sun
1,
Chongwu Zhong
2 and
Guowei Hao
2
1
School of Environment and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China
2
Yankuang Energy Group Company Limited, Jining 273500, China
*
Author to whom correspondence should be addressed.
Sustainability 2024, 16(18), 8039; https://doi.org/10.3390/su16188039
Submission received: 20 July 2024 / Revised: 11 September 2024 / Accepted: 11 September 2024 / Published: 14 September 2024
(This article belongs to the Section Pollution Prevention, Mitigation and Sustainability)
Browse Figures
Figure 1
<p>Location of the study area, working face, and observation line.</p> "> Figure 2
<p>Histogram of error distribution of conventional methods.</p> "> Figure 3
<p>Schematic diagram of surface inclination calculation. The arrows in the figure refer to the direction in which the sampling point is calculated along the x or y direction.</p> "> Figure 4
<p>Schematic diagram of mean filtering in spatial domain. The blue dots represent the four corners of the grid corresponding to each sampling point; the blue arrows point to the images representing the results after processing using the algorithm in this paper.</p> "> Figure 5
<p>Algorithm flow chart. The red dots in the graph represent the direction of the calculation along the x or y direction.</p> "> Figure 6
<p>Distribution of error intervals of sampling points after different grid treatments.</p> "> Figure 7
<p>Error proportion diagram under different grid intervals.</p> "> Figure 8
<p>Mean error of sampling points under different grid intervals.</p> "> Figure 9
<p>Road settlement map measured by UAV based on original data.</p> "> Figure 10
<p>Road settlement map obtained by spatial domain filtering.</p> "> Figure 11
<p>Comparison of path changes for different grid sizes.</p> "> Figure 12
<p>Surface settlement maps obtained through the proposed algorithm. (<b>a</b>) Phase II. (<b>b</b>) Phase III. (<b>c</b>) Phase IV. (<b>d</b>) Phase V. (<b>e</b>) Phase VI. (<b>f</b>) Phase VII.</p> "> Figure 12 Cont.
<p>Surface settlement maps obtained through the proposed algorithm. (<b>a</b>) Phase II. (<b>b</b>) Phase III. (<b>c</b>) Phase IV. (<b>d</b>) Phase V. (<b>e</b>) Phase VI. (<b>f</b>) Phase VII.</p> "> Figure 13
<p>Extracting the isoline map of surface subsidence by UAV. (<b>a</b>) Phase II isoline map of surface subsidence. (<b>b</b>) Phase III isoline map of surface subsidence. (<b>c</b>) Phase IV isoline map of surface subsidence. (<b>d</b>) Phase V isoline map of surface subsidence. (<b>e</b>) Phase VI isoline map of surface subsidence. (<b>f</b>) Phase VII isoline map of surface subsidence.</p> "> Figure 13 Cont.
<p>Extracting the isoline map of surface subsidence by UAV. (<b>a</b>) Phase II isoline map of surface subsidence. (<b>b</b>) Phase III isoline map of surface subsidence. (<b>c</b>) Phase IV isoline map of surface subsidence. (<b>d</b>) Phase V isoline map of surface subsidence. (<b>e</b>) Phase VI isoline map of surface subsidence. (<b>f</b>) Phase VII isoline map of surface subsidence.</p> "> Figure 14
<p>Accuracy assessment for different mining phases. (<b>a</b>) phase IV direction E-W. (<b>b</b>) phase IV direction N-S. (<b>c</b>) phase V direction E-W. (<b>d</b>) phase V direction N-S. (<b>e</b>) phase VI direction E-W. (<b>f</b>) phase VI direction N-S.</p> ">
Review Reports Versions Notes

Abstract

:
Underground coal mining induces surface subsidence, leading to disasters such as damage to buildings and infrastructure, landslides, and surface water accumulation. Preventing and controlling disasters in subsidence areas and reutilizing land depend on understanding subsidence regularity and obtaining surface subsidence monitoring data. These data are crucial for the reutilization of regional land resources and disaster prevention and control. Subsidence hazards are also a key constraint to mine development. Recently, with the rapid advancement of UAV technology, the use of UAV photogrammetry for surface subsidence monitoring has become a significant trend in this field. The periodic imagery data quickly acquired by UAV are used to construct DEM through point cloud filtering. Then, surface subsidence information is obtained by differencing DEM from different periods. However, due to the accuracy limitations inherent in UAV photogrammetry, the subsidence data obtained through this method are characterized by errors, making it challenging to achieve high-precision ground surface subsidence monitoring. Therefore, this paper proposes a spatial domain filtering algorithm for UAV photogrammetry combined with surface deformation caused by coal mining based on the surface subsidence induced by coal mining and combined with the characteristics of the surface change. This algorithm significantly reduces random error in the differential DEM, achieving high-precision ground subsidence monitoring using UAV. Simulation and field test results show that the surface subsidence elevation errors obtained in the simulation tests are reduced by more than 50% compared to conventional methods. In field tests, this method reduced surface subsidence elevation errors by 39%. The monitoring error for surface subsidence was as low as 8 mm compared to leveling survey data. This method offers a new technical pathway for high-precision surface subsidence monitoring in mining areas using UAV photogrammetry.

1. Introduction

China is a major coal resource country, and the increasing demand for mineral resources driven by societal development has exacerbated the adverse impacts of ground subsidence caused by coal mining on human lives [1]. Extensive and prolonged large-scale coal mining leads to the expansion of mined-out areas, creating subsidence zones that not only endanger the lives of residents in mining areas but also trigger a series of disasters, such as landslides, surface water accumulation, and power line collapses [2]. Subsidence disasters in mining areas are major factors restricting regional sustainable development, primarily reflected in infrastructure damage, water resource pollution, and limitations on land use. Therefore, accurately monitoring and providing timely early warnings of ground subsidence caused by mining are key technical approaches to addressing these constraints. Surface subsidence monitoring enables the timely detection of subsidence signs and trends, allowing for prompt action to effectively prevent and mitigate the impact of subsidence disasters. Collecting surface subsidence data provides detailed information on subsidence areas, aiding in the assessment of land reuse potential. Monitoring surface subsidence offers essential data for the sustainable development of mining areas, which is critical for developing rational land use plans and managing water resources, thereby providing robust technical support for the sustainable development of mining areas. Traditional methods for monitoring surface subsidence typically involve establishing surface movement observation stations and using instruments such as leveling instruments, total stations, and Global Navigation Satellite System (GNSS) technology to measure the coordinate changes at monitoring points during the mining process [3,4,5]. While these methods provide high monitoring accuracy, they offer displacement information for only a limited number of observation points and fail to reflect the three-dimensional deformation characteristics of the entire subsidence basin [6]. Moreover, these methods require substantial fieldwork [7].
In recent years, several advanced remote sensing technologies have been applied to monitor surface displacement in mining areas, including Interferometric Synthetic Aperture Radar (InSAR) technology, Light Detection and Ranging (LiDAR) technology, and Unmanned Aerial Vehicle (UAV) photogrammetry. InSAR technology, with its all-weather, all-day monitoring capabilities, high precision, and extensive coverage, significantly mitigates the limitations of traditional displacement monitoring methods [8,9,10]. However, surface deformation caused by underground coal mining is characterized by rapid deformation rates and small subsidence areas, which can degrade the quality of Synthetic Aperture Radar (SAR) images. The large deformation gradient at the center of subsidence basins can cause the decoherence of the interferometric image pairs, making it difficult to obtain accurate displacement information. Although many scholars have proposed methods to address this issue, high-resolution SAR data remain expensive. With the development of UAV technology, UAVs have become a promising solution due to their high flexibility, low cost, and rapid data acquisition capabilities. They can conduct high-precision aerial surveys over mined-out areas, providing extensive surface displacement data and offering a new approach to address this problem [11]. UAV LiDAR systems, which integrate laser scanners, Inertial Measurement Unit (IMU), and differential GNSS, have emerged as effective tools for quickly and accurately obtaining three-dimensional coordinates of the scanned area [12] and collecting terrain data [13,14]. The compact and lightweight nature of UAV LiDAR systems enhances operational flexibility and provides denser point clouds and higher precision [15]. However, LiDAR systems are typically expensive, and current data processing methods often struggle to effectively separate ground points from non-ground points in the point cloud, resulting in reduced monitoring accuracy.
UAV photogrammetry has become a prominent tool for monitoring subsidence in mining areas, owing to its flexibility, cost-effectiveness, and high efficiency. The primary technique for extracting subsidence information from UAV images currently involves differential analysis of the time-series digital elevation model (DEM). This approach leverages multi-temporal UAV aerial images to generate a time-series DEM, enabling the identification of three-dimensional deformation characteristics of subsidence areas through the differentiation of the DEM across different periods [16]. However, the surface settlement data obtained may have large errors due to the effects of constructing models and measurement errors. To mitigate these issues, Zha Jianfeng et al. [17] performed extensive repeated experiments to analyze the distribution regularities of the aerial triangulation point error, specifically focusing on the distribution and precision of Ground Control Points (GCPs). Their findings suggest that the digital averaging of multiple measurements can substantially enhance measurement accuracy. Paweł Cwiakała et al. [18] explored the applicability of UAV photogrammetry data for monitoring surface displacement in regions impacted by underground mining activities, highlighting that vegetation growth significantly affects elevation error. Javier Cardenal et al. [19] employed a multi-view stereo (MVS) technique to generate dense point clouds, which were subsequently utilized to create a Digital Surface Model (DSM) and high-resolution orthophotos. Monitoring points were manually extracted from the orthophotos and DSM, and their horizontal and vertical displacements over different time periods were calculated. The monitoring precision achieved an accuracy of 0.02 ± 0.01 m in the XY direction and 0.04 ± 0.02 in the Z direction. Liu et al. [20] integrated UAV photogrammetry data with LiDAR to create a dynamic subsidence basin based on the differential DEM, co-registering dense point clouds over time. This method achieved a root mean square error (RMSE) ranging from 0.2 to 0.3 m, with a minimum error of 0.17 m. Lian et al. [21] applied the Adaptive Irregular Triangulation (ATIN) algorithm to filter the overall point cloud data of the study area, resulting in a digital elevation model, and used the two-phase digital elevation model to subtract to obtain the dynamic subsidence basin of the surface, and the average RMSE was 165 mm when compared with the same period of on-site measured level data. Zhou Dawei et al. [22] utilized UAV photogrammetry for short-period monitoring of mined-out areas, deriving short-period DEM and calculating the dynamic surface subsidence basin with an accuracy of 121 mm. Yang Xuting et al. [23] investigated the impact of topographic changes on image resolution, photo overlap, and monitoring accuracy, proposing UAV altitude adjustments according to terrain variations. This approach effectively addressed the limitations of traditional UAV aerial surveys in complex terrains. Edyta Puniach et al. [24] utilized a Weighted Normalized Cross-Correlation (WNCC) algorithm for matching multi-temporal orthoimages, thereby determining horizontal ground displacements. Based on the matching results, both horizontal and vertical displacements were calculated to form an initial spatial displacement field. Subsequently, a series of filtering steps were applied to remove outliers, yielding a filtered spatial displacement field. Compared with the observatory data, the displacements were determined with an accuracy of 2–3 cm for the horizontal portion and 3–5 cm for the subsidence portion. Liu et al. [25] revealed the spatial and temporal characteristics of surface-forming subsidence by aligning multi-period UAV densely matched point clouds and then quantifying the surface deformation using the Model-to-Cloud Comparison Algorithm (M3C2), which has an average accuracy of 45 mm, 45 mm, and 64 mm in the X, Y, and Z directions, respectively. Evangelia Siafali et al. [26] conducted a comparative analysis between DSM and DEM to identify and quantify displacement on trail surfaces. ArcGIS software 10.8.1 was utilized to process the DSM data, detect and categorize displacement, and calculate the area and volume of each displacement category to assess the alteration of the trail surface. The methodology employed achieved standard deviations of 0.023 m, 0.020 m, and 0.020 m in the X, Y, and Z directions, respectively, with the corresponding estimated RMSE being 0.023 m, 0.020 m, and 0.020 m. The above studies show the widespread application of UAV photogrammetry in surface subsidence monitoring. Nonetheless, the existing methodologies exhibit limited precision in monitoring mining-induced subsidence. The limitation of measurement accuracy hinders the practical application of UAV photogrammetry technology in subsidence monitoring in mining areas.
In response to this issue, this paper proposes a method to improve measurement accuracy by reducing the resolution of UAV images. In this paper, we first design a simulation experiment to verify the feasibility of the method by verifying it at the experimental site and analyzing the accuracy improvement effect of the method. To better apply this method to subsidence monitoring in mined-out areas, UAV images were acquired from the 3308 working face (means an area where coal blasting, cutting, and transportation are carried out during the coal mining process) of a mine over seven periods between March 2021 and January 2022. The corresponding surface subsidence was calculated, and the feasibility and accuracy of the method were verified using data from surface movement observation stations.

2. Overview of the Area

2.1. Overview of the Experimental Area

To verify the feasibility of the method, a simulated experiment was conducted at a selected site within the China University of Mining and Technology (CUMT). The experimental site is located on Xingjian West Road, Nanhu Campus of CUMT, with geographical coordinates of 117°8′42″ E and 34°13′12″ N. The study area covers approximately 1200 m2, primarily consisting of asphalt roads, cement roads, and trees. The experimental area is an asphalt road with a smooth surface unobstructed by other vegetation or structures, which can serve as an ideal experimental site.

2.2. Overview of the Study Area

To validate the practical performance of the method proposed in this paper, a study area was selected in a coal mine located in Shandong Province. Most of the land above the 3308 working face is cultivated, and there are scattered buildings; the study area is shown in Figure 1. Coal mining activities can adversely affect surface soil and vegetation, leading to land subsidence, soil erosion, and reduced vegetation cover, among other issues. These problems can cause a decline in groundwater levels, resulting in drought conditions across the surrounding land and plant withering. To mitigate the impact on local agricultural production and residents’ lives, as well as to enhance the sustainable use of local land, it becomes especially important to obtain the dynamic process of surface subsidence during coal mining. This information will provide critical data support for subsequent studies on land degradation, groundwater depletion, and other related aspects.
A DJI Genie 4 RTK UAV (Shenzhen Dajiang Innovation Technology Co., Ltd., Shenzhen, China) was used to collect image data from workface 3308 in the study area. The coal mine here employs comprehensive coal mining technology, with surface elevations ranging from +43.18 to +44.12 m and working face elevations from −225.3 to −274.9 m. The strike length of the coal seam at the working face is 687 m, and the dip length is 197 m. The seam’s inclination angle varies from 2° to 7°, and the coal seam thickness ranges from 4.1 to 5.6 m, with an average thickness of 5.18 m.

3. Materials

3.1. Data Sources

3.1.1. UAV Orthophoto Data

To achieve the purpose of verifying the accuracy of the algorithms in this paper as well as monitoring the displacement generated on the ground surface by the working surface, the two regions mentioned above were acquired separately. To ensure the stability of the experimental data, all the data for this study were collected through a DJI Phantom 4 RTK small multi-rotor high-precision aerial survey UAV (Shenzhen Dajiang Innovation Technology Co., Ltd., Shenzhen, China). Detailed information is provided as shown in Table 1.
In pursuit of verifying the accuracy of the algorithm, two sets of image data were taken in the simulated experimental area, and the flight path is parallel to the observation area. The settings were adjusted to ensure an 85% directional overlap and a 75% lateral overlap amid adjacent pictures, maintaining a flight altitude of 50 m, and the Ground Sampling Distance (GSD) was 1.4 cm/pix. Primary operational information is presented in Table 2. A total of 11 measurement points were arranged within the experimental zone. From these, six served as GCPs, while the remaining points acted as Checkpoints (CPs). Both GCPs and CPs were demarcated using right-angled triangles sprayed with red paint. Each GCP was gauged using a Hi-Target ZTS-121R Total Station (Hi-Target, Canton, China), with planar coordinates obtained through traverse surveying and elevation coordinates obtained via trigonometric leveling. The planar measurement accuracy was set at ±2 mm, while the elevation accuracy up to 150 m was also set at ±2 mm.
To monitor the surface displacement caused by the extraction of the 3308 working face, this study collected imagery data from the research area over seven different periods between March 2021 and January 2022, as detailed in Table 3. The GCPs used in the experiment consisted of rigid control boards and right-angle triangles marked with red paint. The planar coordinates of the control points were measured using RTK (Shanghai Huace Navigation Technology Ltd., Shanghai, China) with an accuracy of 3 cm, and the elevation coordinates were obtained through third-order leveling with an accuracy of 1 mm.

3.1.2. Validate Data

To verify the reliability of surface subsidence data obtained from processed UAV images, this study established surface movement observation stations above the 3308 working face and measured leveling elevation for accuracy testing. Based on field research, 26 leveling points were placed at 50 m intervals in the north–south direction, labeled Z1 to Z26, and 18 leveling points at 50 m intervals in the east–west direction, labeled Q1 to Q18. The layout of the surface movement observation stations is shown in Figure 1. Due to environmental conditions in the study area, some monitoring points were damaged, with Z3 to Z9 and Q10 to Q14 being the primary points that remained intact over the long term. During the experiment, third-order leveling was used to measure the elevation of the monitoring points, achieving an accuracy of 1 mm. Periodic monitoring of these points enabled the calculation of subsidence at various stages of mining progress.

4. Principles and Methods

4.1. Data Characteristics

As a collection reflecting the surface position information of objects within the same coordinate system, the point cloud data processed from UAV orthophotos exhibit the following characteristics: (1) The accuracy of aerial triangulation from UAV imagery is significantly influenced by the precision of GCPs, meaning that GCP accuracy directly affects the precision of surface subsidence monitoring by UAV; (2) The data processed from UAV image contain substantial error, impacting measurement accuracy. Based on these characteristics, previous researchers have analyzed the error distribution in processed UAV imagery by ensuring the distribution and accuracy of GCPs and using CPs. They found that the errors in the planar coordinates and elevation coordinates of UAV measurement points follow a normal distribution but are subject to certain systematic error. The error distribution is characterized by bounded and dense features [17]. Given this characteristic, it is possible to digitally average the results of multiple measurements to improve accuracy. However, considering that multiple measurements increase labor intensity, this study designed relevant experiments to investigate this feature.

4.2. Error Distribution of Subsidence Basin Established by Conventional Method

Conventional UAV image processing methods for constructing subsidence basins are significantly affected by random error. Therefore, this paper proposes an algorithm for this problem. Currently, the conventional method is to process the image of each period to the DSM, then obtain the DEM by filtering, and then find the difference between the DEM of different periods to obtain the surface subsidence basin. To determine whether this algorithm aligns with the characteristics of random error, an analysis is conducted through the first experiment. Given the short interval between the two data collection times, the surface subsidence value can be assumed to be zero. The surface subsidence basin is then obtained using a conventional method. For point cloud filtering, the commonly used Cloth Simulation Filter (CSF) filtering method for DSM is employed.
The error distribution histogram derived from conventional methods is illustrated in Figure 2. As seen in the figure, the error does not tend toward zero, with maximum error reaching up to 70 mm. The induced error can be attributed to factors such as random error, systematic error, and filtering precision. Through the error distribution histogram, it can be seen that its distribution has a certain regularity, that is, in the same measurement mode, its error presents the characteristics of random distribution, but there is also a certain systematic error influence.
Based on the characteristics of random error, their impact can be reduced by averaging multiple measurements. However, conducting multiple observations would increase the workload of field measurements. To minimize human resource investment, this study applies the following principles to simulate multiple observations: (1) Within a specific range, if the surface displacement is minimal and falls below the UAV’s monitoring precision, it can be assumed that all monitoring points within this range show no displacement and are observations of the same precision. (2) For regions with equal displacement and monitoring precision, averaging within the spatial domain is equivalent to multiple observations, thereby reducing the influence of random error. Using these principles, this study proposes enhancing precision by reducing resolution, thus mitigating the impact of random error.

4.3. Deformation-Adapted Spatial Domain Filtering Algorithm

This paper proposes a deformation-adapted spatial domain filtering algorithm, which obtains a filtering grid by utilizing the slope values between sample points in the DEM and then digitally averaging the settlement values of the sample points within the grid to obtain a processed DEM. The regenerated DEM is then used to construct the surface subsidence basin by superimposing DEM generated at different periods. The workflow of the algorithm is as follows:
(1)
DSM acquisition. UAV orthophotos were imported into Pix4DMapper software 2.0.104, and the DSM of the study area was generated through the steps of image feature point matching, Structure from Motion (SFM), multi-view stereo (MVS), and point cloud modeling;
(2)
Point cloud filtering. To remove non-ground points from the point cloud, the CFS filtering was used to separate the non-ground points from the DSM of the study area, retaining the ground point cloud data;
(3)
Hole filling. Since some regions are left empty after point cloud filtering, the empty parts are interpolated by using Kriging interpolation.
(4)
DEM reconstruction. The filtered ground point cloud data are rasterized to regenerate a higher precision DEM;
(5)
Define the size of the grid. Let us imagine the original DEM data as S , where the subsidence value of each sampled point is S ( m , n ) and the size of the corresponding filter window centered around each sample point is S Z × Z . To achieve the conditions of adaptive mining deformation, the grid size can be determined based on the accuracy of the monitoring of drone subsidence and the magnitude of the tilt between neighboring sampling points. Using Formulas (1)–(3) to calculate the tilt value of the surface sampling points and using Formula (4) to calculate the size of the grid divided by these surface sampling points, i.e., the size of the filter grid,
i x ( m , n ) = S ( m , n ) − S ( m , n − 1 ) l m , n − ( m , n − 1 )
i y ( m , n ) = S ( m , n ) − S ( m − 1 , n ) l m , n − ( m − 1 , n )
i = i x 2 + i y 2
Z = ∆ d i
In the above formula, S represents the subsidence value of each sampling point, l m , n − ( m , n − 1 ) is the distance between two sampling points, i x and i y are the slope values in the east–west and north–south directions, respectively, i is the slope value of the sampling point, ∆ d is the precision of UAV elevation measurement, and Z is the filter window size of the sampling point. Figure 3 illustrates the calculation of the surface slope.
(6)
Spatial domain filtering treatment for adaptive mining deformation. Let us define the data output after filtering therapy in this study as g. When the grid of each sampling point is divided according to Formula (4), since the subsidence difference in all sampling points within the grid is less than the measurement error of the UAV, i.e., the difference in subsidence values of all sampling points within each filtering window is negligible, the sampling points within the grid can be averaged and used as a substitute for the center of the grid, and the averaging of all the points is equivalent to the multiple observations on the sampling points, which helps to reduce the influence of random error on UAV monitoring. Through Formula (5), the filtering process for the point S ( m , n ) can be completed, and the resultant output is g ( m , n ) , as shown in Figure 4.
g m , n = ∑ i ∈ m − Z 2 ,     m + Z 2 , j ∈ n − Z 2 ,     n + Z 2 S ( i , j ) Z × Z
When changes occur in the UAV’s monitoring accuracy and the surface slope value, the size of the filtering window corresponding to the sampling points also changes accordingly. Repeat the above steps until the filtering is completed for all the sampling points and finally output the surface subsidence data g after reducing the random error. The algorithm in this study employs an adaptive principle when calculating the filtering window size for each sampling point. It not only fully utilizes the correlation within the data itself but also takes into account its inherent physical meaning, that is, the coupling of surface deformation and UAV measurement accuracy. This allows for a more reasonable selection of the filter window size, and it has a noticeable effect on reducing the random error in the surface subsidence data.
(7)
Accuracy Validation. In order to validate the method to enhance the accuracy, the validation is carried out by calculating the RMSE by means of Bessel’s correction. The RMSE is calculated according to Formula (6).
E H = ∑ 1 n ∆ H i 2 n
In the equation, ∆ H represents the settlement values at various sampling points, while E H denotes the RMSE of the settlement values.
The specific algorithm flowchart of this paper is shown in Figure 5.

5. Results and Discussion

5.1. Simulation Experimental Results

To verify the feasibility of the algorithm presented in this study, the UAV image data obtained from the experimental area are processed using the aforementioned method. This generates a DEM after two periods of differences, and subsidence values are uniformly extracted from the different DEMs at intervals of 2 cm (the original resolution of the UAV). Because there is no actual subsidence in the selected test area, there is no inclination value between the sampling points obtained. Therefore, after considering the common range of surface slope values in the mining area, this experiment chooses to divide the filtering windows at intervals of 0.1 m, 0.5 m, and 1 m and analyzes the root mean square error size of the sampling points and the error distribution situation after filtering through different-sized windows.
Figure 6, Figure 7 and Figure 8 show that the proportion of settlement error within a 3 cm range exceeds 95%. Additionally, it can be observed that most settlement errors are concentrated within the −5 to 5 mm range. As the grid division range increases, the errors become more centralized. Figure 8 also demonstrates that as the grid division range increases, the root mean square error of the sampling points decreases, indicating that this method effectively reduces the root mean square error.
If the data processed through different grid divisions are used to regenerate DEM data, and the newly generated DEM is analyzed and compared using three-dimensional plots, results are obtained as shown in Figure 9 and Figure 10.
From the figure above, it is evident that the original data of the UAV contain numerous errors, making it challenging to ensure measurement accuracy without processing. The results obtained using the method proposed in this paper effectively reduce the impact of random error. On analyzing Figure 9 and Figure 10, it can be seen that there are a large number of errors in the subsidence basin obtained using the conventional method, and there are mainly a large number of random errors. In this paper, we propose to reduce the resolution of the UAV to minimize the impact of random error by digitally averaging the subsidence data in the grid without the need for multiple measurements to reduce the random error. On the one hand, it reduces the workload and complex data post-processing work in the field; on the other hand, it can provide more accurate monitoring information.
Because the experimental area is a road, and the planar data contain a large amount of information, the data points are densely distributed. To simplify the data and more intuitively depict the deformation characteristics and trends in the experimental area road, the central axis from different grid data sets is extracted for analysis, as shown in Figure 11. It can be observed that as the grid size increases, the profile of this path tends towards smoothness. However, it is also evident that despite the path profile becoming smoother, the subsidence value is not zero. Therefore, it can be inferred that under this measurement mode, the value approached by this curve represents the systematic error of the UAV.
The above experiments demonstrate that grid subdivision based on monitoring accuracy and tilt values, followed by mean processing, can effectively reduce random measurement error. This indicates that the spatial domain filtering method proposed in this study is feasible and can significantly improve the measurement accuracy of UAV data.

5.2. Case Application Results

To evaluate the effectiveness of the proposed method in practical scenarios, a study was conducted on the 3308 working face of a mine in Shandong Province. Given the presence of vegetation in the research area, it was necessary to remove the vegetation before extracting settlement data. After setting the corresponding color threshold according to the RGB value of the vegetation in the point cloud data of each period, the vegetation points and other noises in the ground point cloud were eliminated. To acquire more precise surface subsidence data, the gaps in the point cloud after the filtering process were filled using the Kriging method. Subsequently, the surface point cloud for each period was rasterized to re-generate the DEM with higher accuracy. By calculating the difference between the initial and final DEMs, subsidence values for points in the region were uniformly extracted. The initial subsidence point coordinates ( x ,   y , z ) of the sampling points were derived along the east–west and north–south directions, respectively, resulting in subsidence value files sorted in two different directions. By traversing the sampling points, the subsidence value between two neighboring points in the north–south direction and the distance between them were calculated to determine the north–south tilt value for each point. Similarly, the east–west tilt value for each point was computed. Once the initial tilt values for all points were determined, the grid size surrounding each point was adjusted based on the UAV’s measurement accuracy and the initial tilt value of the sampling points. The average subsidence values of all points within each grid were then taken as the actual subsidence value for each point. The processed sampling points were then used to reconstruct an elevation triangulated irregular network (TIN) to generate the DEM. The processed data were comprehensively analyzed, and the sinking basin for each phase (time period for data collection and analysis carried out at a specific point in time or during a specific period of time in the mining process) are shown Figure 12. The positive values in the Figure 12 represent the surface subsidence values.
As the above figure illustrates, the range of surface displacement above the extraction face gradually expanded due to the continuous excavation, resulting in a steady increase in the maximum surface settlement value.
To visually demonstrate the dynamic changes and impact range of the surface subsidence basin processed using the proposed algorithm, contour maps based on the surface subsidence data are plotted, as shown in Figure 13. This figure clearly illustrates the extent and magnitude of the influence on the surface above the mined-out area as the working face advances. This allows for a detailed analysis of the surface deformation characteristics due to underground coal mining and its impact on surrounding structures.

Accuracy Assessment

By comparing the subsidence basin formed using conventional methods, the superiority of the algorithm in this study is evident. To verify the accuracy improvement achieved by the proposed algorithm in practical applications, the subsidence basins established under different measurement modes are compared with settlement data acquired from surface displacement observation stations. Due to the environmental conditions in the study area, a number of monitoring sites have been damaged, of which Z3 to Z9 and Q10 to Q14 are the main long-term intact sites, which are used for accuracy assessment. The comparison of surface settlement accuracy for different phases is plotted as shown in Figure 14. Verification of the accuracy of the UAV measurement of surface subsidence using leveling data as the true value of surface subsidence is carried out. Since the above monitoring points are not monitored or not deformed in phases I to III, the data from stages IV to VI are selected for comparison in this paper, and the calculation shows that the algorithm has a significant improvement in accuracy. Under the experimental conditions of the study area selected in this paper, the results are as follows: the root mean square error for measuring surface settlement using original UAV data for the first and fourth periods is 13 mm. After processing with the proposed method, the root mean square error reduces to 8 mm, enhancing the measurement accuracy by 38.5%. Likewise, for the first and fifth periods, the root mean square error reduces from 18 mm to 11 mm—an improvement of 38.9%; and for the first and sixth periods, the error decreases from 20 mm to 12 mm—an improvement of 40.0%, as shown in Table 4. It can be noted that with an increase in the settlement value, there is a certain decline in UAV’s measurement accuracy. However, the UAV-based settlement measurement method proposed in this study consistently results in substantial accuracy improvements.
From the above results, it is shown that the method proposed in this paper has unique advantages compared with conventional methods. The method proposed in this paper can automatically adjust the grid size according to the deformation characteristics of the ground surface, adapt to different terrains, and can be applied to most of the mining areas or other ground subsidence areas. In addition, the method can obviously improve the monitoring accuracy and can better meet the requirements of surface subsidence monitoring in mining areas.
However, the method in this paper has the following shortcomings: (1) The method in this paper is based on the DEM constructed by point cloud filtering, and the error of point cloud filtering will be accumulated in the surface subsidence basin after the subsequent processing of the method in this paper. (2) The method of removing the vegetation above the cultivated land based on the RGB value of the point cloud mentioned in this paper has certain limitations, and it only applies to the case where there is a small amount of discrete vegetation above the land surface. There is a large amount of vegetation, and it is difficult to ensure the accuracy of point cloud hole interpolation after the removal. (3) For the removal of the original point cloud coarseness and noise, different experimental conditions have different choices of point cloud filtering methods, and it is difficult to generalize. (4) Although the method shows a significant improvement on the accuracy of surface subsidence monitoring, the accuracy may be limited in the area of particularly large-scale or especially intense subsidence. This is something we need to further explore, improving the algorithm in our future work to ensure that high-accuracy monitoring results can be obtained under various conditions.

6. Conclusions

(1)
In this paper, we analyze the error distribution regulation in UAV measurement in depth with the previous research. We specifically research the error distribution characteristics in UAV elevation measurement while maintaining constant factors such as UAV flight control parameters and the layout and quality of GCPs. The results of the research show that the elevation measurement errors in UAV surveys follow a normal distribution, albeit influenced by the presence of systematic errors. The error distribution is characterized by bounded and dense.
(2)
Based on the characteristics of UAV measurement error distribution, this paper proposes a spatial domain filtering algorithm for UAV with mining deformation adaptation. As the random error can be averaged over several measurements to reduce the influence of random error, the algorithm in this paper avoids the huge workload caused by several measurements and achieves the effect of averaging over several measurements by reducing the resolution of the UAV. For the grid size division of sampling point, this paper utilizes the surface changes caused by coal mining surface subsidence combined with the drone measurement accuracy to automatically change the grid. Through the experiment, it can be seen that the method of this paper can greatly reduce the influence of random error in the measurement accuracy of the drone and can better meet the requirements of surface deformation monitoring in mining areas.
(3)
In the actual case application, the algorithm of this paper and the level data, in which the RMSE can be up to 8 mm, are compared with the present conventional method to obtain the RMSE of 13 mm, and the accuracy is improved by about 40%, which greatly improves the accuracy of UAV monitoring. This method provides a new technical path for accurately obtaining the deformation information of the coal mine subsidence area.

Author Contributions

Conceptualization, J.Z. and P.M.; methodology, J.Z. and P.M.; software, P.M.; validation, P.M., H.L., M.Y. and B.S.; formal analysis, J.Z.; investigation, P.M., H.L. and M.Y.; resources, J.Z., C.Z. and G.H.; data curation, P.M.; writing—original draft preparation, P.M.; writing—review and editing, J.Z. and P.M.; project administration, J.Z.; funding acquisition, J.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by China University of Mining and Technology (Xuzhou, China) through National Key Research and Development Program of China (2023YFC3804201) and the Natural Science Foundation of China (NSFC) (U21A20109).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

Author Chongwu Zhong and Guowei Hao was employed by the Yankuang Energy Group Company Limited. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

  1. Li, D.; Deng, K.; Gao, X.; Niu, H. Monitoring and Analysis of Surface Subsidence in Mining Area Based on SBAS-InSAR. Geomat. Inf. Sci. Wuhan Univ. 2018, 43, 1531–1537. (In Chinese) [Google Scholar]
  2. Chen, Y.; Tao, Q.; Liu, G.; Wang, L.; Wang, F.; Wang, K. Detailed mining subsidence monitoring combined with InSAR and probability integral method. Chin. J. Geophys. 2021, 64, 3554–3566. (In Chinese) [Google Scholar]
  3. Yang, J.; Jiang, Y.; Zhou, J.; Huang, L.; Lu, X. Analysis on reliability and accuracy of subsidence measurement with gps technique. J. Geod. Geodyn. 2006, 26, 70–75. (In Chinese) [Google Scholar]
  4. Liu, G.; Zhang, L.; Cheng, S.; Jiang, T. Feasibility Analysis of Monitoring Mining Surface Substance Using InSAR/GPS Data Fusion. Bull. Surv. Mapp. 2005, 11, 13–16. (In Chinese) [Google Scholar]
  5. Wang, Z.; Zhang, J.; Huang, G. Precise monitoring and analysis of the land subsidence in Jining coal mining area based on InSAR technique. J. China Univ. Min. Technol. 2014, 43, 169–174. [Google Scholar]
  6. Deng, J.; Xu, Y. Three-dimensional Dynamie Monitoring and Analysis of Surface Deformation in Mining Area Based on Multi-track SAR Image. Metal Mine 2019, 10, 48–54. (In Chinese) [Google Scholar]
  7. Hu, L.; Navarro-Hernández, M.I.; Liu, X.; Tomás, R.; Tang, X.; Bru, G.; Ezquerro, P.; Zhang, Q. Analysis of regional large-gradient land subsidence in the Alto Guadalentín Basin (Spain) using open-access aerial LiDAR datasets. Remote Sens. Environ. 2022, 280, 113218. [Google Scholar] [CrossRef]
  8. Chen, B.; Deng, K.; Fan, H.; Hao, M. Large-scale deformation monitoring in mining area by D-InSAR and 3D laser scanning technology integration. Int. J. Min. Sci. Technol. 2013, 23, 555–561. [Google Scholar] [CrossRef]
  9. Chen, B.Q.; Deng, K.Z. Integration of D-InSAR technology and PSO-SVR algorithm for time series monitoring and dynamic prediction of coal mining subsidence. Surv. Rev. 2014, 46, 392–400. [Google Scholar] [CrossRef]
  10. Dong, L.; Wang, C.; Tang, Y.; Tang, F.; Zhang, H.; Wang, J.; Duan, W. Time Series InSAR Three-Dimensional Displacement Inversion Model of Coal Mining Areas Based on Symmetrical Features of Mining Subsidence. Remote Sens. 2021, 13, 2143. [Google Scholar] [CrossRef]
  11. Liu, G.S.; Lv, J.C. Research on key technologies for monitoring changes in mining areas based on unmanned aerial vehicles. Bull. Surv. Mapp. 2013, 95–98. (In Chinese) [Google Scholar]
  12. Sui, L.; Zhang, Y.; Zhang, S.; Chen, W. Filtering of airborne LiDAR point cloud data based on progressive TIN. Wuhan Daxue Xuebao (Xinxi Kexue Ban)/Geomat. Inf. Sci. Wuhan Univ. 2011, 36, 1159–1163. [Google Scholar]
  13. Polat, N.; Uysal, M. Investigating performance of Airborne LiDAR data filtering algorithms for DTM generation. Measurement 2015, 63, 61–68. [Google Scholar] [CrossRef]
  14. Cook, K.L. An evaluation of the effectiveness of low-cost UAVs and structure from motion for geomorphic change detection. Geomorphology 2017, 278, 195–208. [Google Scholar] [CrossRef]
  15. Liang, X.; Kankare, V.; Hyyppä, J.; Wang, Y.; Kukko, A.; Haggrén, H.; Yu, X.; Kaartinen, H.; Jaakkola, A.; Guan, F.; et al. Terrestrial laser scanning in forest inventories. ISPRS J. Photogramm. Remote Sens. 2016, 115, 63–77. [Google Scholar] [CrossRef]
  16. Mengxia, Z.; Xinqi, Z.; Bo, L. Dynamic Monitoring Method for Mining Area Based on UAV Images. Bull. Surv. Mapp. 2017, 43–47. (In Chinese) [Google Scholar]
  17. Jianfeng, Z.; Pengcheng, Z.; Dejun, W.; Chenyang, M.; Chongwu, Z.; Hongwu, Y.; Dong, W. UAV aerial triangulation: Point error distributions and the influencing mechanisms of ground control points on its accuracy. Coal Geol. Explor. 2023, 51, 151–161. (In Chinese) [Google Scholar]
  18. Ćwiąkała, P.; Gruszczyński, W.; Stoch, T.; Puniach, E.; Mrocheń, D.; Matwij, W.; Matwij, K.; Nędzka, M.; Sopata, P.; Wójcik, A. UAV Applications for Determination of Land Deformations Caused by Underground Mining. Remote Sens. 2020, 12, 1733. [Google Scholar] [CrossRef]
  19. Cardenal, J.; Fernández, T.; Pérez-García, J.L.; Gómez-López, J.M. Measurement of Road Surface Deformation Using Images Captured from UAVs. Remote Sens. 2019, 11, 1507. [Google Scholar] [CrossRef]
  20. Liu, J.; Liu, X.; Lv, X.; Wang, B.; Lian, X. Novel Method for Monitoring Mining Subsidence Featuring Co-Registration of UAV LiDAR Data and Photogrammetry. Appl. Sci. 2022, 12, 9374. [Google Scholar] [CrossRef]
  21. Lian, X.-g.; Liu, X.-y.; Ge, L.; Hu, H.F.; Du, Z.; Wu, Y.-r. Time-series unmanned aerial vehicle photogrammetry monitoring method without ground control points to measure mining subsidence. J. Appl. Remote Sens. 2021, 15, 024505. [Google Scholar] [CrossRef]
  22. Zhou, D.; Qi, L.; Zhang, D.; Zhou, B.; Guo, L. Unmanned Aerial Vehicle (UAV) Photogrammetry Technology for Dynamic Mining Subsidence Monitoring and Parameter Inversion: A Case Study in China. IEEE Access 2020, 8, 16372–16386. [Google Scholar]
  23. Yang, X.; Yao, W.; Zheng, J.; Ma, B.; Ma, X. UAV terrain following technology application in the mining subsidence monitoring research. Bull. Surv. Mapp. 2021, 111–115. (In Chinese) [Google Scholar]
  24. Puniach, E.; Gruszczyński, W.; Stoch, T.; Mrocheń, D.; Ćwiąkała, P.; Sopata, P.; Pastucha, E.; Matwij, W. Determination of the coefficient of proportionality between horizontal displacement and tilt change using UAV photogrammetry. Eng. Geol. 2023, 312, 106939. [Google Scholar] [CrossRef]
  25. Liu, X.; Zhu, W.; Lian, X.; Xu, X. Monitoring Mining Surface Subsidence with Multi-Temporal Three-Dimensional Unmanned Aerial Vehicle Point Cloud. Remote Sens. 2023, 15, 374. [Google Scholar] [CrossRef]
  26. Siafali, E.; Tsioras, P.A. An Innovative Approach to Surface Deformation Estimation in Forest Road and Trail Networks Using Unmanned Aerial Vehicle Real-Time Kinematic-Derived Data for Monitoring and Maintenance. Forests 2024, 15, 212. [Google Scholar] [CrossRef]
Sustainability 16 08039 g001
Figure 1. Location of the study area, working face, and observation line.
Figure 1. Location of the study area, working face, and observation line.
Sustainability 16 08039 g001
Sustainability 16 08039 g002
Figure 2. Histogram of error distribution of conventional methods.
Figure 2. Histogram of error distribution of conventional methods.
Sustainability 16 08039 g002
Sustainability 16 08039 g003
Figure 3. Schematic diagram of surface inclination calculation. The arrows in the figure refer to the direction in which the sampling point is calculated along the x or y direction.
Figure 3. Schematic diagram of surface inclination calculation. The arrows in the figure refer to the direction in which the sampling point is calculated along the x or y direction.
Sustainability 16 08039 g003
Sustainability 16 08039 g004
Figure 4. Schematic diagram of mean filtering in spatial domain. The blue dots represent the four corners of the grid corresponding to each sampling point; the blue arrows point to the images representing the results after processing using the algorithm in this paper.
Figure 4. Schematic diagram of mean filtering in spatial domain. The blue dots represent the four corners of the grid corresponding to each sampling point; the blue arrows point to the images representing the results after processing using the algorithm in this paper.
Sustainability 16 08039 g004
Sustainability 16 08039 g005
Figure 5. Algorithm flow chart. The red dots in the graph represent the direction of the calculation along the x or y direction.
Figure 5. Algorithm flow chart. The red dots in the graph represent the direction of the calculation along the x or y direction.
Sustainability 16 08039 g005
Sustainability 16 08039 g006
Figure 6. Distribution of error intervals of sampling points after different grid treatments.
Figure 6. Distribution of error intervals of sampling points after different grid treatments.
Sustainability 16 08039 g006
Sustainability 16 08039 g007
Figure 7. Error proportion diagram under different grid intervals.
Figure 7. Error proportion diagram under different grid intervals.
Sustainability 16 08039 g007
Sustainability 16 08039 g008
Figure 8. Mean error of sampling points under different grid intervals.
Figure 8. Mean error of sampling points under different grid intervals.
Sustainability 16 08039 g008
Sustainability 16 08039 g009
Figure 9. Road settlement map measured by UAV based on original data.
Figure 9. Road settlement map measured by UAV based on original data.
Sustainability 16 08039 g009
Sustainability 16 08039 g010
Figure 10. Road settlement map obtained by spatial domain filtering.
Figure 10. Road settlement map obtained by spatial domain filtering.
Sustainability 16 08039 g010
Sustainability 16 08039 g011
Figure 11. Comparison of path changes for different grid sizes.
Figure 11. Comparison of path changes for different grid sizes.
Sustainability 16 08039 g011
Sustainability 16 08039 g012aSustainability 16 08039 g012b
Figure 12. Surface settlement maps obtained through the proposed algorithm. (a) Phase II. (b) Phase III. (c) Phase IV. (d) Phase V. (e) Phase VI. (f) Phase VII.
Figure 12. Surface settlement maps obtained through the proposed algorithm. (a) Phase II. (b) Phase III. (c) Phase IV. (d) Phase V. (e) Phase VI. (f) Phase VII.
Sustainability 16 08039 g012aSustainability 16 08039 g012b
Sustainability 16 08039 g013aSustainability 16 08039 g013b
Figure 13. Extracting the isoline map of surface subsidence by UAV. (a) Phase II isoline map of surface subsidence. (b) Phase III isoline map of surface subsidence. (c) Phase IV isoline map of surface subsidence. (d) Phase V isoline map of surface subsidence. (e) Phase VI isoline map of surface subsidence. (f) Phase VII isoline map of surface subsidence.
Figure 13. Extracting the isoline map of surface subsidence by UAV. (a) Phase II isoline map of surface subsidence. (b) Phase III isoline map of surface subsidence. (c) Phase IV isoline map of surface subsidence. (d) Phase V isoline map of surface subsidence. (e) Phase VI isoline map of surface subsidence. (f) Phase VII isoline map of surface subsidence.
Sustainability 16 08039 g013aSustainability 16 08039 g013b
Sustainability 16 08039 g014
Figure 14. Accuracy assessment for different mining phases. (a) phase IV direction E-W. (b) phase IV direction N-S. (c) phase V direction E-W. (d) phase V direction N-S. (e) phase VI direction E-W. (f) phase VI direction N-S.
Figure 14. Accuracy assessment for different mining phases. (a) phase IV direction E-W. (b) phase IV direction N-S. (c) phase V direction E-W. (d) phase V direction N-S. (e) phase VI direction E-W. (f) phase VI direction N-S.
Sustainability 16 08039 g014
Table 1. UAV parameters.
Table 1. UAV parameters.
ProjectParametersProjectParameter
UAV modelDJI Phantom 4 RTKCamera nameFC6130r
Maximum flight time/min30Camera pixel20 megapixels
Camera sensorOne inch COMOS
Pixel width/px5472
Pixel high/px3648
Table 2. Flight parameters.
Table 2. Flight parameters.
TimeAltitude/mGSD/cmNumber of Images per IssueRoute OverlapLateral Overlap
23 June 2023501.410485%75%
Table 3. UAV operation information.
Table 3. UAV operation information.
DateAltitude/mGSD/cmOverlap RateNumber of ImagesWeather
4 March 2021250.780%/70%1109Sunny
1 April 2021250.780%/70%1091Sunny
18 April 2021250.780%/70%1100Cloudy
21 May 2021250.780%/70%1109Sunny
10 June 2021250.780%/70%1091Cloudy
13 July 2021250.780%/70%1100Sunny
11 January 2022250.780%/70%1100Sunny
Table 4. Measurement accuracy improvement results at different phases.
Table 4. Measurement accuracy improvement results at different phases.
PhaseRMSE in UAV Original Data Measurements (mm)RMSE in UAV Proceed Data Measurements (mm)Precision Improvement
Phase IV13838.5%
Phase V181138.9%
Phase VI201240.0%
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zha, J.; Miao, P.; Ling, H.; Yu, M.; Sun, B.; Zhong, C.; Hao, G. Deformation-Adapted Spatial Domain Filtering Algorithm for UAV Mining Subsidence Monitoring. Sustainability 2024, 16, 8039. https://doi.org/10.3390/su16188039

AMA Style

Zha J, Miao P, Ling H, Yu M, Sun B, Zhong C, Hao G. Deformation-Adapted Spatial Domain Filtering Algorithm for UAV Mining Subsidence Monitoring. Sustainability. 2024; 16(18):8039. https://doi.org/10.3390/su16188039

Chicago/Turabian Style

Zha, Jianfeng, Penglong Miao, Hukai Ling, Minghui Yu, Bo Sun, Chongwu Zhong, and Guowei Hao. 2024. "Deformation-Adapted Spatial Domain Filtering Algorithm for UAV Mining Subsidence Monitoring" Sustainability 16, no. 18: 8039. https://doi.org/10.3390/su16188039

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop