Construction of Mining Subsidence Basin and Inversion of Predicted Subsidence Parameters Based on UAV Photogrammetry Products Considering Horizontal Displacement

Abstract

Constructing high-precision subsidence basins is of paramount importance for mining subsidence monitoring. Traditional unmanned aerial vehicle (UAV) photogrammetry techniques typically construct subsidence basins by directly differencing digital elevation models (DEMs) from different monitoring periods. However, this method often neglects the influence of horizontal displacement on the accuracy of the subsidence basin. Taking a mining area in Ordos, Inner Mongolia, as an example, this study employed the normalized cross-correlation (NCC) matching algorithm to extract horizontal displacement information between two epochs of a digital orthophoto map (DOM) and subsequently corrected the horizontal position of the second-epoch DEM. This ensured that the planar positions of ground feature points remained consistent in the DEM before and after subsidence. Based on this, the vertical displacement in the subsidence area (the subsidence basin) was obtained via DEM differencing, and the parameters of the post-correction subsidence basin were inverted using the probability integral method (PIM). The experimental results indicate that (1) the horizontal displacement was influenced by the gully topography, causing the displacement within the working face to be segmented on both sides of the gully; (2) the influence of the terrain on the subsidence basin was significantly reduced after correction; (3) the post-correction surface subsidence curve was smoother than the pre-correction curve, with abrupt error effects markedly diminished; (4) the accuracy of the post-correction subsidence basin increased by 43.12% compared with the total station data; and (5) comparing the measured horizontal displacement curve with that derived using the probability integral method revealed that the horizontal displacement on the side of an old goaf adjacent to the newly excavated working face shifted toward the advancing direction of the new working face as mining progressed. This study provides a novel approach and insights for using low-cost UAVs to construct high-precision subsidence basins.

1. Introduction

As China’s primary energy source, coal has significantly contributed to the country’s economic development [1]. However, geological hazards, such as surface subsidence, ground fissures, and landslides, caused by coal mining pose severe threats to the ecological environment of mining areas and the safety of residential zones [2]. Traditional mining area monitoring methods include leveling, trigonometric leveling, and global navigation satellite system (GNSS) measurements. These methods typically set up strike and dip observation lines above goaf areas to monitor subsidence [3]. Although these traditional techniques offer high accuracy, they have drawbacks, such as high labor intensity, high costs, challenging operational conditions, and the inability to capture surface deformations over large areas [4,5]. With the continuous advancement of surveying and mapping technologies, interferometric synthetic aperture radar (InSAR), UAV photogrammetry, and light detection and ranging (LiDAR) technologies have been successfully applied in mining subsidence monitoring [6]. InSAR, in particular, has attracted attention due to its high resolution, all-weather capability, wide coverage, and high accuracy [7,8]. However, the high gradient deformation in the center of subsidence basins often causes decorrelation in interferometric images, leading to phase unwrapping failures and the inability to obtain true deformation information [9,10]. In contrast, the UAV method has become increasingly important for mining area monitoring due to its efficiency, flexibility, and adaptability to various terrains. Although UAV LiDAR excels in vegetation penetration and high-precision data acquisition, its high cost limits its widespread application. Conversely, UAV photogrammetry offers a low cost, ease of operation, and high-resolution imagery [11,12], making it widely used in mining subsidence monitoring [13].

In recent years, with the successful application of UAV photogrammetry in mining monitoring, scholars have focused on improving the quality of UAV imagery and the accuracy of surface subsidence monitoring. Zhang [14] enhanced the quality of original images by optimizing UAV photography schemes. Additionally, Tang et al. [15] improved subsidence monitoring accuracy by eliminating subsidence data from vegetated areas and re-interpolating the data. After comparing five filtering algorithms, Lian et al. [16] found that the adaptive triangulated irregular network (TIN) filtering algorithm performed best in extracting ground points. They used this algorithm to generate a DEM of the mining area and then constructed a subsidence basin model via direct differencing. However, the current UAV photogrammetry method primarily relies on directly differencing multi-temporal DEMs to construct the subsidence basin. This method only reflects surface elevation changes at the same geographical coordinates, neglecting the impact of horizontal displacement on vertical deformation [17]. Therefore, it is essential to simultaneously obtain horizontal displacement information related to vertical deformation to construct a high-precision subsidence basin model.

Partial progress has been made in the extraction of horizontal displacements in mining areas. Zhan et al. [18] treated ground feature point clouds as rigid bodies, combined the fast point feature histogram (FPFH) and iterative closest point (ICP) algorithms to register the point clouds, and extracted 3D surface movements from the transformation matrix. He et al. [19] extracted the horizontal displacements of the main cross-section using airborne LiDAR point cloud data and improved the binary shape contextual feature description operator for the topographic features of the mining area. Yang et al. [20] extracted horizontal displacements from slope feature maps before and after subsidence using image sub-pixel correlation technology. Yang et al. [21] proposed a weighted hue-based ICP (WHICP) algorithm to estimate 3D surface displacements using multi-temporal color point clouds generated from UAV photogrammetry. Current research primarily focuses on extracting 3D deformation from point cloud data. High-resolution DOMs captured simultaneously contain more feature information, which has increasingly attracted attention from scholars. In the fields of glacier dynamics, landslides, and earthquakes, the NCC algorithm has been successfully applied to extract horizontal surface displacement from high-resolution imagery captured before and after deformation events. Scambos et al. [22] used an NCC-based image matching algorithm to achieve sub-pixel image registration of satellite imagery, generating a high-precision velocity field for glacier movement. Zhao et al. [23] applied the open-source software MicMac1.0 (which employs the NCC algorithm) to extract horizontal displacement in landslide areas before and after deformation. Similarly, Provost et al. [24] used the MicMac image matching library to perform sub-pixel image registration on Sentinel-2 satellite imagery, producing high-resolution surface displacement fields to estimate co-seismic displacement following earthquakes. These applications provide new insights for horizontal displacement extraction in mining subsidence studies. Puniach et al. [25] used high-resolution DOMs to determine the horizontal displacement induced by underground mining, and by comparing various image-matching algorithms, found that the weighted normalized cross-correlation (WNCC) algorithm achieved an accuracy of 1–2 pixels for horizontal displacement extraction. Zhu [26] extracted the horizontal displacement by matching corresponding feature points between two DOMs using the scale-invariant feature transform (SIFT) algorithm. Zhu et al. [27] proposed the probability integral method-based prior-guided COSI-Corr (PIM-CC) method, which connects horizontal displacement with subsidence, thereby effectively improving the accuracy of horizontal displacement extraction in mining subsidence monitoring. However, surface horizontal displacement and subsidence in mining areas are caused by underground mining-induced deformation, making it difficult to fully extract the 3D surface deformation solely from high-resolution DOMs.

Based on this, this study proposes an innovative approach to constructing mining subsidence basins, incorporating both horizontal displacement and vertical deformation using UAV-derived DOM and DEM data. This method consists of two key steps: (1) extracting the horizontal displacement between two epochs of DOMs using the NCC method, and (2) correcting the horizontal position of the second-epoch DEM based on the extracted horizontal displacement and constructing the subsidence basin by differencing the corrected second-epoch DEM and the first-epoch DEM. Taking the 2S201 working face of a mining area in Ordos, Inner Mongolia, as the study area, the effectiveness and accuracy of the proposed method were validated by comparing it with traditional methods for constructing subsidence basins. Additionally, the surface movement rule caused by multi-working face mining in a loess gully region was analyzed, and subsidence prediction parameters were inverted using the probability integral method based on data from the entire basin.

2. Research Areas and Data Introduction

2.1. Overview of Research Areas

Ordos City is located in the southwest of the Inner Mongolia Autonomous Region of China. Its underground space is rich in coal resources and is one of China’s important energy bases [28]. This study focused on the 2S201 working face of a mining area in Ordos City, Inner Mongolia, as the research area. Figure 1 illustrates the geographical location of the study area. The mining area is situated in the northeastern part of the Ordos Basin, belonging to the Ordos structural erosion plateau, with an overall topography that is higher in the north and west and lower in the south and east. The elevation ranges from 1000 m to 2000 m. The mining area features a complex and diverse terrain, surrounded by the Yellow River on three sides in the northeast and west and connecting with the Loess Plateau in the south. The landscape includes a variety of geomorphological types, such as grasslands, plateaus, sandy lands, and hills. The study area comprises various geomorphological units, including eroded depositional landforms, accumulative landforms, and fluvial terrains, reflecting the region’s diverse landscape. Overall, the area exhibits an undulating topography, with the southern part characterized by well-developed dendritic gullies, forming a typical structural erosion plateau landscape.

The 2S201 working face has a strike length of 1252 m, a dip length of 260 m, an average mining thickness of 3.3 m, an average depth of 200 m, and a dip angle of 2°. Mining began in mid-July 2018 and concluded in late October 2018. The surface primarily comprises arsenic sandstone, characterized by weak rock formations, poor cementation between sand grains, and low structural strength. The sandstone turns into mud upon contact with water and disintegrates into sand when exposed to wind. Prolonged, large-scale, and high-intensity mining activities have led to varying degrees of ground subsidence, resulting in structural damage, such as fissured houses and collapsed roads. This region is also highly susceptible to landslides and mudslides, further exacerbating the area’s ecological vulnerability.

2.2. Data

In this study, high-resolution imaging data of the mining area were collected using a Trimble UX5 UAV equipped with a SONY A5100 DSLR camera. The specific parameters of the UAV imaging acquisition are listed in

Table 1. Parameters for UAV photography collection.

No.	Acquisition Date	UAV	Camera	Course Overlap (%)	Lateral Overlap (%)	Flight Altitude (m)
1	9 June 2018	Trimble UX5	SONY A5100	80	80	230
2	16 April 2019	Trimble UX55	SONY A5100	80	80	230

. The ground control points (GCPs) coordinate system is based on the Beijing 1954 coordinate system, and GCP coordinates were obtained using GPS RTK measurements. To verify the accuracy of 3D surface deformation measurements obtained from UAV photogrammetry, 11 ground monitoring points were established within the study area, with their locations shown in Figure 1c. Ground data were collected on 9 June 2018 and 16 April 2019, using a Leica TSO9 robotic total station with an angular accuracy of 2″, ensuring that the deformation monitoring accuracy met the requirements for mining engineering surveys.

2.3. Data Preprocessing

The aerial images captured with the UAV were imported into the UASMaster8.0 software for data processing. The main workflow included image distortion correction, image pyramid construction, relative orientation, absolute orientation, and aerial triangulation to obtain point clouds and the DOM of the monitoring area. The point cloud data generated with the UAV included ground points, various buildings in the mining area, and vegetation. The point cloud data were denoised and filtered to reduce noise and the influence of non-ground points on the DEM accuracy. This study used an improved progressive triangulated irregular network densification filtering algorithm [29,30] to remove non-ground points such as vegetation and buildings, thereby extracting ground points. The IDW algorithm was then used to interpolate these ground points to generate the required DEM.

The obtained DOM and DEM were then registered to ensure precise alignment within the same coordinate system and achieve a consistent ground resolution. Following this, the common monitoring region of the DOM and DEM was cropped to ensure consistent coverage. This process provided accurate data for subsequent horizontal displacement extraction and the construction of subsidence basins that considered horizontal displacement.

3. Methodology

The surface imagery captured via UAV photogrammetry contains rich three-dimensional coordinate information. High-resolution DOMs and DEMs of the monitoring area can be generated by applying advanced data-processing techniques. The NCC algorithm was introduced in this study to extract the horizontal displacement information between two epochs of high-resolution DOMs. This displacement information was used to correct the horizontal position of the second DEM to ensure consistency in the planar coordinates of ground objects before and after subsidence. Subsequently, a subsidence basin model was constructed via DEM differencing to obtain the vertical subsidence values of each ground object. This study’s technical workflow is shown in Figure 2.

3.1. Horizontal Displacement Extraction Method

The extraction of the surface horizontal displacement from high-resolution imagery before and after deformation using image sub-pixel correlation techniques has been widely applied in fields such as glacier movement, landslides, and earthquakes [31,32,33,34]. Monitoring the surface horizontal displacement in mining areas is critical for protecting surface structures and the early warning of geological hazards, such as surface cracks and landslides. In this study, we primarily utilized the MM2DPosSism interface in the open-source software MicMac1.0 to extract the horizontal displacement between the before and after deformation DOMs of the mining area. MicMac applies regularization techniques to ensure reliable results, even when using small correlation windows. It utilizes a multi-resolution approach, starting with low-resolution calculations and progressively increasing the resolution at each matching level. The NCC algorithm was additionally utilized, offering options for non-linear cost functions and adjustable correlation thresholds, thereby preventing low-correlation areas from impacting the accuracy of high-correlation regions [35].

The NCC algorithm is widely used in template matching. In this process, a matching window of size $n\times n$ is constructed by selecting a neighborhood around the pixel location $(p_x,p_y)$ in the reference image. Similarly, an identical matching window is constructed at the target pixel location $(p_x+u,p_y+v)$ in the target image, where variables $u$ and $v$ represent the horizontal and vertical displacements between two images, specifically of $I_2$ relative to $I_1$. The correlation between these two matching windows is then computed using the following:

$$NCC(u,v)={\displaystyle \frac{{\displaystyle \sum _{(x,y)\in W_p}\left(I_1\left(p_x,p_y\right)-{\overline{I}}_1\left(p_x,p_y\right)\right)}\ast \left(I_2\left(p_x+u,p_y+v\right)-{\overline{I}}_2\left(p_x+u,p_y+v\right)\right)}{\sqrt{{\displaystyle \sum _{(x,y)\in W_p}{\left(I_1\left(p_x,p_y\right)-{\overline{I}}_1\left(p_x,p_y\right)\right)}^2}\ast {\left(I_2\left(p_x+u,p_y+v\right)-{\overline{I}}_2\left(p_x+u,p_y+v\right)\right)}^2}}}$$

(1)

In Equation (1), $NCC(u,v)$ represents the normalized cross-correlation coefficient, with a range of values from [−1, 1]. $W_p$ denotes the matching window; $I_1(p_x,p_y)$ represents the pixel value within the matching window in the reference image; and ${\overline{I}}_1(p_x,p_y)$ is the average pixel value within the matching window in the reference image. $I_2(p_x+u,p_y+v)$ represents the pixel value within the matching window of the target image, and ${\overline{I}}_2(p_x+u,p_y+v)$ is the average pixel value in the target image’s matching window.

By default, MicMac uses the following linear correlation cost formula: $Cost=1-cor$. A non-linear cost formula is applied to reduce the impact of poorly correlated areas—such as those with snow or water features lacking clear contours—on the calculation.

$$C_1=\mathit{Max}\left(Cor,C^{min}\right)$$

(2)

$$C_2={\displaystyle \frac{C_1-C^{min}}{1-C^{min}}}$$

(3)

$$C_3=C_2^{\gamma }$$

(4)

$$Cost=\left(1-C_3\right)\ast \left(1-C^{min}\right)$$

(5)

where $C\mathit{or}$ is the normalized cross-correlation coefficient, and $C^{min}$ is the correlation threshold. With a higher $\gamma$, the influence of the correlation close to 1 is higher. The parameters selected for this study are shown in

Table 2. Calculation parameters for horizontal displacement.

Parameter	Value
$C^{min}$	0.5
$\gamma$	2
Size of window	9 × 9

.

3.2. Subsidence Basin Construction Method

Currently, subsidence basin models constructed by differencing two-epoch DEMs do not consider the impact of horizontal displacement on subsidence values; as such, the obtained surface subsidence values only reflect the change in ground elevation at the same geographic coordinates, rather than the true vertical displacement of ground feature points [17]. We propose an improved method for constructing the subsidence basin model to address this limitation. Using the planar positions of ground feature points from the first-epoch DEM as the reference, this method first corrects the positions of ground feature points in the second DEM to align with their corresponding positions in the first DEM by using the east–west and north–south horizontal displacement information extracted from two epochs of a DOM. Subsequently, the subsidence basin model is constructed by differencing the two-epoch DEMs. The specific correction process is as follows: In Figure 3, points A and A′ represent the same surface object before and after horizontal displacement, with ${\mathsf{\Delta }}_{EW}$ and ${\mathsf{\Delta }}_{NS}$ representing the horizontal displacement in the east–west and north–south directions, respectively. To accurately calculate the vertical displacement of point A, we first correct the position of point A′ in the second DEM to align with the position of point A in the first DEM using the horizontal displacement information. Then, the elevation value corresponding to point A′ is interpolated at the location of point A in the second DEM using the inverse distance weighting (IDW) interpolation method. The other ground feature points are processed likewise to construct the first-epoch DEM corresponding to each feature point in the second-epoch DEM. Finally, the two-epoch DEMs are differenced to construct a subsidence basin model that eliminates the effects of horizontal displacement.

3.3. Parameter Inversion Using the Probability Integral Method

Common methods for predicting mining subsidence include the typical curve, profile function, and probability integral methods. The probability integral method, based on stochastic medium theory, is effective for predicting mining subsidence parameters and has been widely applied to predict mining subsidence in coal mines across China [36,37].

(1)

The mathematical modeling of the subsidence value $W(x,y)$ at any point $A(x,y)$ on the surface is expressed as follows:

$$W(x,y)={\displaystyle \frac{1}{W_0}}{W}^0(x)W^0(y)$$

(6)

$$W^0(x)={\displaystyle \frac{W_0}{2}}\left[erf({\displaystyle \frac{\sqrt{\pi }}{r}}x)-erf({\displaystyle \frac{\sqrt{\pi }}{r}}(x-l)\right]$$

(7)

$$W^0(y)={\displaystyle \frac{W_0}{2}}\left[erf({\displaystyle \frac{\sqrt{\pi }}{r_1}}y)-erf({\displaystyle \frac{\sqrt{\pi }}{r_2}}(y-L)\right]$$

(8)

where $W_0$ represents the maximum subsidence value of the ground surface under full-mining conditions ($W_0=mq\cos \alpha$); $m$ is the mining thickness; $q$ is the subsidence coefficient; $\alpha$ is the coal seam dip angle; and $W^0(x)$ and $W^0(y)$ represent the surface subsidence values in the main section under limited mining conditions in the strike and dip directions, respectively. $l=D_3-S_3-S_4$ is the calculated length of the strike working face, where $D_3$ is the strike length of the working face, and $S_3$ and $S_4$ are the left and right deviations of the inflection point, respectively. $L=(D_1-S_1-S_2){\displaystyle \frac{\sin (\theta +\alpha )}{\sin \theta }}$ is the calculated length of the dip working face, where $D_1$ is the dip length of the working face, and $S_1$ and $S_2$ are the upper and lower deviations of the inflection point, respectively. Here, $\theta$ represents the Propagation Angle, while $\alpha$ is the dip angle of coal seam. $r$, $r_1$, and $r_2$ represent the primary influence radii for the strike, lower, and upper directions, respectively. $erf(x)$ is the Gaussian error function, where $erf(x)={\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle {\int }_0^x{e}^{-x^2}}dx$.

(2)

The mathematical modeling of the horizontal displacement $u(x,y,\phi )$ in the $\phi$ direction at any surface point $A(x,y)$ is

$$U(x,y,\phi )={\displaystyle \frac{1}{W_0}}[U^0(x)W^0(y)\cos \phi +U^0(y)W^0(x)\sin \phi ]$$

(9)

$$U^0(x)=b{W}_0\left(e^{-\pi {\left({\displaystyle \frac{x}{r}}\right)}^2}-e^{-\pi {\left({\displaystyle \frac{\left(x-l\right)}{r}}\right)}^2}\right)$$

(10)

$$U^0(y)=W_0\left(b_1{e}^{-\pi {\left({\displaystyle \frac{y}{r_1}}\right)}^2}-b_2{e}^{-\pi {\left({\displaystyle \frac{\left(y-L\right)}{r_2}}\right)}^2}\right)+\cot \theta \cdot W^0(y)$$

(11)

where $b$, $b_1$, and $b_2$ are the horizontal displacement coefficients.

This study used the least-squares method for surface fitting, iteratively solving for the parameters of the 2S201 working face via multiple iterations. The objective was to minimize the sum of squared deviations between the fitted surface and the observed subsidence basin. The mathematical model is as follows [38]:

$$\left.\begin{array}{l}{V}_1{V}_1^T={\displaystyle \sum _{k=1}^n{\left[W_m-W(x,y)\right]}^2}=\min \\ V_2{V}_2^T={\displaystyle \sum _{k=1}^n{\left[U_m-U(x,y)\right]}^2}=\min \end{array}\right\}$$

(12)

where $V_1{V}_1^T$ and $V_2{V}_2^T$ represent the deviations between the observed subsidence and horizontal displacement values and the corresponding least-squares fitted values, respectively. $W_m$ and $U_m$ represent the observed subsidence and horizontal displacement values, respectively. $W(x,y)$ and $U(x,y)$ represent the least-squares fitted subsidence and horizontal displacement values, respectively.

3.4. Accuracy Assessment Methods

3.4.1. Accuracy Assessment Method of Horizontal Displacement

The accuracy of the extracted horizontal displacement was evaluated using the ground monitoring point data obtained by the total station as a reference. The root-mean-square error (RMSE) was calculated using the following formula:

$$M=\pm \sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n{\mathsf{\Delta }}_i^2}}{n}}}$$

(13)

where $M$ represents the RMSE; ${\mathsf{\Delta }}_i$ is the displacement difference in the east–west or north–south directions; and $n$ is the number of monitoring points.

3.4.2. Accuracy Assessment Method of Subsidence Basin

(1)

Verification of Internal Coincidence Accuracy

The PIM is described as “the most commonly used method” in the 2017 edition of the “Code for Coal Pillar Retention and Coal Mining in Buildings, Water Bodies, Railways and Main Shafts and Lanes” [39]. Therefore, this study evaluated the internal coincidence accuracy of the subsidence basin models before and after correction by analyzing the residuals between the pre-correction and post-correction subsidence basin models and the parameter-inverted subsidence basin model.

(2)

Accuracy Verification Based on Total Station Monitoring Points

In this study, subsidence values at ground monitoring points obtained from total station measurements serve as benchmarks to evaluate the accuracy of subsidence values in subsidence basins before and after correction. The correction accuracy is assessed by calculating the RMSE.

4. Results Analysis

4.1. Horizontal Displacement Extraction

Disparity maps containing the east–west and north–south horizontal displacement information were generated using the proposed method with the precisely registered and cropped DOM from two different periods, as shown in Figure 4. Figure 4a and Figure 4b present the east–west and north–south horizontal displacements in the study area, respectively, with the maximum horizontal displacement reaching 1.1 m. The horizontal displacement was mainly distributed along the boundaries of the working face and in areas with significant topographic variation. The direction of the horizontal displacement at the boundaries of the working face pointed toward the center of the working face. The influence of the newly excavated working face 2S202 at the intersection of the two working faces caused the horizontal displacement on the side adjacent to working face 2S202 to point toward the center of 2S202. Cross-sections A-A′ and B-B′ were drawn to further study the relationship between the horizontal displacement and topography in areas with significant elevation changes within the 2S201 working face. The horizontal displacement, slope, and slope of the slope (SOS) along these cross-sections were extracted to observe their correlations. Figure 5a–d show that the horizontal displacement suddenly changed when the slope or the slope of the slope was large.

The study area is a typical hilly and gully region, including various terrain elements such as flat land, slopes, and steep cliffs [40]. The analysis revealed that horizontal displacements are distributed not only along the boundary of the subsidence basin but also within the working face. A combined analysis of the horizontal displacement in the east–west and north–south directions, along with the topography, revealed that the gully topography significantly affected the horizontal displacement. As shown in Figure 6a, the A-A′ profile intersected with the gully at the 385 m mark, with horizontal displacement primarily in the east–west direction, reaching a maximum displacement value of −0.688 m, until it reached the gully bottom at 455 m, where the east–west displacement trended toward 0. Between 265 m and 385 m, the two adjacent gullies caused the east–west horizontal displacement to shift toward the bottom of the gully. As shown in Figure 6b, when the profile crossed the bottom of the gully and intersected with the other side (530 m), the horizontal displacement primarily occurred in the north–south direction, progressively increasing from the gully bottom to the top, where it reached the maximum displacement of 0.873 m. In the section of the A-A′ profile between 1000 m and 1600 m, the vertical profile displacement (EW) gradually approached zero as the gully disappeared, similar to Wang et al.‘s findings [41]. In summary, the gully topography influenced the horizontal displacement, causing the displacement within the working face to be segmented on both sides of the gully. The horizontal displacement reached its maximum at the top of the gully and gradually decreased toward the bottom, wherein it eventually trended toward zero.

4.2. Subsidence Basin Construction

Conventional UAV photogrammetry techniques typically construct subsidence basins by directly differencing two periods of a DEM, as shown in Figure 7a. However, this method does not account for the impact of horizontal displacement on the accuracy of the subsidence basin. The first step in constructing a subsidence basin that considers horizontal displacement is to precisely register the generated DEM and DOM and crop them to ensure they have the same coverage, guaranteeing data consistency. Then, the planar position of the second-epoch DEM is corrected to align with the same feature points in both DEMs using the extracted horizontal displacement information. Based on this, the vertical displacement of each feature point in the monitored area is calculated with the difference between the corrected second-epoch DEM and the first-epoch DEM, constructing a more accurate subsidence basin model for the mining area, as shown in Figure 7b. The results indicate that the maximum subsidence in the basin’s center was approximately 2.95 m.

Figure 7a,b and the local maps of areas I and II in Figure 8 demonstrate that the subsidence values in areas I and II in the uncorrected subsidence basin were significantly influenced by the topography; meanwhile, the topographic influence on the subsidence values in these areas was significantly reduced in the corrected subsidence basin. Cross-sections A-A′ and C-C′ were plotted for further analysis, and the subsidence curves of the main section of the subsidence basin before and after correction in area I were obtained, as shown in Figure 9. Figure 9 reveals the following findings: (1) The post-correction A-A′ and C-C′ subsidence curves exhibit greater stability than the pre-correction curves, with less fluctuation, further demonstrating that the post-correction subsidence basin effectively mitigated the influence of the topography. (2) In the gully terrain region (Zone I) along profile A-A′ (385 m to 530 m), the pre-correction subsidence curve exhibits significant fluctuations along with the gully topography, while the post-correction subsidence curve shows a more gradual trend. As shown in Figure 9b, the gully topography in region I between 150 m and 300 m of the C-C′ profile caused the pre-correction subsidence curve to fluctuate more in the gully area than in the flat topography. However, the subsidence curve trends were consistent in both the gully and flat topography areas after correction. In summary, the proposed method for constructing subsidence basins demonstrates better applicability than traditional methods.

4.3. Parameter Estimation of Subsidence Prediction Using the Probability Integral Method

The traditional PIM for parameter estimation typically relies on monitoring point data along the strike and dip profiles. However, arranging the required strike and dip observation lines in the complex topography of mining areas is challenging, making parameter estimation difficult. UAV photogrammetry can address this issue by providing planar data. Previous studies [42,43,44] have verified the feasibility of estimating subsidence parameters using the PIM across the entire basin. This study extracted estimated subsidence parameters for the pre-correction and post-correction subsidence basins. The results shown in

Table 3. Parameter inversion results of probability integral method.

Probability Integral Method Parameters	Subsidence Coefficient $\mathit{q}$	Tangent of Major Influence Angle $\mathbf{tan}\mathit{\beta }$	Propagation Angle ${\mathit{\theta }}_0$ (°)	Horizontal Displacement Coefficient $\mathit{b}$	Deviation of Inflection Point (m)				Ratio of Fitting Mean Square Error to Maximum Measured Subsidence Value (%)
					S1	S2	S3	S4
Pre-correction	0.92	1.7	89	0.27	33	42	36	26	8.0
Post-correction	0.89	1.8	89	0.27	33	42	35	25	6.2

indicate the fitting mean square errors of the post-correction and uncorrected subsidence basins accounted for 6.2% and 8.0% of the maximum subsidence value, respectively. This indicates that the parameters obtained post-correction had significantly higher accuracy than those obtained pre-correction, further validating the effectiveness of horizontal displacement correction in improving the accuracy of subsidence basin models.

These parameters were used to simulate the subsidence basin in the mining area to further validate the reliability of the predicted parameters obtained using the PIM, as shown in Figure 10. This simulated basin was compared with the measured subsidence basin constructed using UAV photogrammetry, as shown in Figure 11. The strike and dip profiles of both the measured and simulated basins were extracted for further analysis, as depicted in Figure 12. The subsidence curves derived from the parameter fitting and the UAV-measured subsidence curves shown in Figure 12a,b exhibit a similar trend, further confirming the reliability of the estimated parameters.

The predicted parameters were used to invert the horizontal displacements in the east–west and north–south directions. The strike and dip profiles of the measured and inverted horizontal displacements were extracted for analysis, as shown in Figure 13 and Figure 14, respectively, wherein the overall trend of the measured and predicted horizontal displacements is similar. However, large horizontal displacements were observed in the center of the working face due to the influence of the topography [45], which deviated from the predicted displacement, as highlighted in the circled area in Figure 13. The comparison of the predicted horizontal displacement curves using the PIM and UAV-measured horizontal displacement curves revealed that the horizontal displacement curves inverted using the superposition method for multiple working faces significantly differed from the UAV-measured displacement curves at the junction of the two working faces, as shown in Figure 14. As indicated in Figure 14b, due to the mining of the 2S202 working face, the horizontal displacement on the side adjacent to the 2S202 working face shifted toward the center of the 2S202 working face, consistent with Yan’s findings [46]. In conclusion, the surface points on the side adjacent to the newly excavated working face tended to shift toward the advancing direction of the new face as it progressed. However, the horizontal displacement model derived by simply superimposing multiple working faces struggled to represent this characteristic accurately. Further research is needed on the horizontal displacement inversion of multiple working faces using the PIM.

4.4. Accuracy Assessment

4.4.1. Accuracy Assessment of Horizontal Displacement

The horizontal displacement error calculation results are shown in Figure 15. The minimum error in the north–south direction was 0.052 m, the maximum error was 0.302 m, the average error was 0.162 m, and the RMSE was 0.172 m. In the east–west direction, the minimum error was −0.023 m, the maximum error was −0.254 m, the average error was −0.167 m, and the RMSE was 0.178 m.

4.4.2. Accuracy Assessment of Subsidence Basin

(1)

Verification of Internal Coincidence Accuracy

In this study, horizontal displacements measured at ground monitoring points using a total station serve as benchmarks to assess the accuracy of horizontal displacements extracted via the NCC algorithm. According to the residual statistics for the entire basin, the strike profile, and the dip profile (shown in Figure 16, Figure 17 and Figure 18, respectively), the internal coincidence accuracy of the post-correction data was significantly better than that of the pre-correction data. The standard deviations were reduced by 89.2 mm, 115.1 mm, and 86 mm, respectively. This indicates that the correction process effectively reduced random errors in the measurement data, significantly enhancing the data consistency.

(2)

Accuracy Verification Based on Total Station Monitoring Points

The subsidence values derived from subtracting the elevations of ground monitoring points observed by the total station over two periods were compared with the subsidence values from the pre-correction and post-correction subsidence basins to assess the data precision more accurately before and after correction, respectively. The RMSE was calculated, and the results are shown in Figure 19. The RMSE of the subsidence values pre-correction was 0.218 m, which significantly decreased to 0.124 m post-correction. This result indicates that the correction process effectively reduced measurement errors and improved the data accuracy. When the influence of horizontal displacement was not considered, the monitoring points with larger discrepancies compared with the total station measurements were points 2, 3, 6, 10, and 11, with residuals of −0.398 m, −0.296 m, 0.260 m, −0.260 m, and −0.258 m, respectively. After accounting for horizontal displacement, these discrepancies were reduced to −0.216 m, 0.012 m, 0.216 m, −0.170 m, and −0.061 m, respectively, showing a significant improvement in accuracy and eliminating some errors caused by topography fluctuations. This demonstrates that the corrected data aligned more closely with the total station measurements, further validating the correction method’s effectiveness.

5. Discussion

In recent years, research on surface deformation caused by underground mining has primarily focused on vertical subsidence [6,47,48], while studies on horizontal displacement have been relatively scarce. However, horizontal displacement can easily lead to surface cracks, which may trigger landslide-related geological hazards [23,49]. Therefore, monitoring horizontal displacement and analyzing its deformation patterns is of great significance for the protection of surface buildings and structures. Currently, some researchers have proposed utilizing the abundant ground feature points captured in UAV-generated DOMs to extract horizontal displacement information in mining areas [25,50]. However, studies that use these horizontal displacement data, in conjunction with terrain features and mining information, to conduct in-depth analyses of horizontal displacement patterns are still lacking. This study introduces the NCC algorithm to extract horizontal displacement information between high-resolution DOMs captured before and after deformation. Based on the study area’s typical hilly and gully topography, an in-depth analysis of the extracted horizontal displacement information revealed the following: (1) the direction of horizontal displacement is influenced by the gully topography, with horizontal displacements on both sides of the gully moving toward the gully bottom; and (2) horizontal displacement is affected by the proximity of adjacent working faces, causing horizontal displacement on the side closer to the newly excavated working face to shift in the advancing direction of the new face.

Furthermore, from the analysis of the subsidence basin, it can be observed that traditional methods for constructing subsidence basins [5,42,51,52] are significantly affected by topography when obtaining subsidence values, leading to rapidly increasing errors in subsidence values in gully topography areas and possibly resulting in completely inaccurate subsidence values. In contrast, the method proposed in this study for constructing the subsidence basin effectively reduces the topography’s influence, as shown in Figure 7 and Figure 8, thereby achieving more accurate subsidence values. This provides more precise information to support the protection of the local ecological environment and human life and production.

Despite the advantages demonstrated by this study, several issues still need further improvement in future research: (1) there is uncertainty in extracting horizontal displacement in mining areas with large vegetation coverage; (2) the horizontal displacement patterns at the junctions of multiple working faces could be analyzed in greater depth using time-series data; and (3) further validation and optimization are required to assess this method’s applicability under different mining conditions.

6. Conclusions

When using DEM differencing to obtain subsidence information in mining areas, the traditional UAV photogrammetry technique does not strictly ensure that the obtained subsidence values represent the vertical displacement of the same feature points, leading to errors caused by horizontal displacement. Therefore, this paper proposed an improved method in response to this issue. This method utilizes sub-pixel image correlation technology to extract the horizontal displacement information between two epochs of a DOM. The planar position of the second-epoch DEM is corrected based on the extracted horizontal displacement information, ensuring that the same feature points in both DEMs are aligned in the vertical direction. This approach allows for the construction of a more accurate subsidence basin model. This study was demonstrated using a case in Inner Mongolia. The key findings are as follows:

(1)

The NCC algorithm was used to calculate the correlation between two epochs of a DOM, thereby extracting the horizontal displacement between the two DOMs. The results show that the maximum horizontal displacement was 1.1 m, with an RMSE of 0.172 m in the east–west direction and 0.178 m in the north–south direction. The horizontal displacement in the monitoring area conformed to the surface movement pattern of the mining subsidence.

(2)

Analyzing the correlation between the horizontal displacement profile and the slope and slope of the slope revealed that the topography significantly affected the horizontal displacement, especially in areas with steep slopes featuring a pronounced gradient change.

(3)

Using a case from a mine in Ordos, Inner Mongolia, the horizontal displacement extracted from two epochs of a DOM was used to correct the planar position of the second-epoch DEM. On this basis, the constructed subsidence basin showed a significant reduction in the influence of the topography. A comparative analysis of the subsidence basin profiles before and after correction revealed that the post-correction surface subsidence curve exhibited a smoother characteristic than the pre-correction curve, effectively reducing the impact of sudden error. The accuracy of the post-correction subsidence basin was improved by 43.12% compared with the total station monitoring data, providing a new method and perspective for constructing high-precision subsidence basins.

(4)

The results of utilizing the “planar” data of the entire basin for the inversion analysis of the parameters for mining subsidence prediction using the PIM indicate that the fitting mean square error of the inverted parameters for the post-correction subsidence basin accounted for 6.2% of the maximum subsidence value. Additionally, comparing the UAV-measured and inverted subsidence curves along the strike and dip profiles revealed that the subsidence trends were largely consistent, further validating the inverted parameters’ reliability.