SBAS-InSAR technology is developed based on InSAR technology, first proposed by Berardino in 2002 (Berardino et al., 2002). This method enhances the accuracy of deformation monitoring, and subsequently, many researchers have supplemented and further studied this technology, applying it to surface deformation monitoring and other areas (Singhroy et al., 2015; Zhao et al., 2016; Zhang et al., 2018; Yang et al., 2020; Wang et al., 2022). Due to its reliability and accuracy, this method has been widely adopted. Figure 1 illustrates the geometric representation of InSAR imaging.

After acquiring N+1 scenes of SAR imagery using this technology, M interferograms are generated based on the spatiotemporal baseline threshold, with specific constraints on their quantity satisfying Equation 1 (Li S. et al., 2022; Wang et al., 2022). N + 1 2 ≤ M ≤ N N + 1 2 (1)

This design maximizes the utilization of temporal information while controlling the risk of spatial decorrelation. For the interferogram j formed between time points t _A and t _B (where t _A < t _B ), its differential interferometric phase δϕ _j can be expressed as a superposition of multiple phases, as shown in Equation 2: δ ϕ j = 4 π λ ∆ d j + 4 π λ B ⊥ j ∆ z R s i n θ + ∆ ϕ a t m , j + ϕ n o i s e , j (2)

Where λ is the radar wavelength; B _⊥j is the vertical baseline; Δz represents the error of the Digital Elevation Model (DEM), indicating the difference between the elevation data and the actual terrain; R is the distance from the radar to the target point; and θ is the incidence angle.

Equation 3 is a system of equations constructed from the phases processed from all interference patterns: B v = δ ϕ (3)

When the coefficient matrix B is of full rank, the deformation rate v can be solved using the least squares method. In cases where B is rank-deficient, the Singular Value Decomposition (SVD) method is employed to obtain the deformation rate. Once the deformation rate is determined, the deformation displacement can be calculated based on the time series of the SAR images.

The numerical model simulated using MatDEM version 3.24 is generally composed of numerous particles. Each particle can be set with a certain moisture content. By configuring the contact model and range between particles, moisture can be assigned to each particle, ensuring that the overall moisture distribution of the model aligns with the moisture distribution patterns of the rock layers (Yu H. et al., 2023). After computing one time step, the moisture content ω of each particle is defined as follows: ω = ω 0 · 1 − β (4)

In the Equation 4, ω ₀ represents the water content at the end of the previous time step for the unit; β is the rate of water reduction.

Based on the contact relationships of each particle, the moisture content of each particle is defined sequentially. Finally, by considering the impact of the moisture field on the mechanical properties of the soil, a particle model with moisture content can be obtained. Figure 2 illustrates the transport of moisture between particles.

In Figure 2, “Particle ① → Particle ②” indicates that, driven by the hydraulic gradient, water flows from the water-rich soil particle to the water-deficient soil particle. Using the finite difference method, let the unit time be dt, and the amount of moisture transfer from Particle ① to particle ②, dM₁₂ , can be expressed by Equation 5. d M 12 = k · ω 1 − ω 2 · A · d t (5)

In the equation, k represents the permeability coefficient, and A is the equivalent area of the permeation path, taken as the average area of two particles.

The process of rainwater adhering to the surface of the rock and soil is characterized by the infiltration of rainwater, which compresses and replaces the water-gas two-phase flow within the slope’s rock and soil. Equation 6 is the infiltration formula for the soil surface in the Green-Ampt model (Ma et al., 2011; Van den Putte et al., 2013) that takes into account the ponding condition, used to calculate the rainfall attachment flow rate i_ω for slope soil particles: i ω = K S 1 + θ S − θ i U α I q (6)

In the equation, i _ω represents the rainfall attachment flow rate for soil particles onthe slope; θ _S is the saturated water content; θ _i is the initial water content; K _S denotes the saturated permeability coefficient; U _α is the water suction at the wetting front; and I _q represents the rainfall intensity.

When the water content of soil particles increases, the internal friction angle decreases. To adjust the internal friction angle of the particles, the Hoek-Brown hydraulic parameters can be introduced to weaken the rock layers (Michalowski and Park, 2020), as shown in Table 1.

The above calculation method is programmed in MatDEM, and then the water transport of rock mass under heavy rainfall conditions and the deformation rule of rock strata are studied.

The permeability of unsaturated soil is a function of volumetric water content, and the relationship between soil water saturation and matric suction is described by the Soil-Water Characteristic Curve (SWCC). According to the curve function proposed by Brooks and Corey, the relative permeability function based on effective saturation can be expressed using the soil matric suction (Brooks and Corey, 1963). Fredlund and Xing recommended the following function to fit the Soil-Water Characteristic Curve for soils with matric suction below 1 × 10⁶ kPa (Fredlund et al., 1994). The curve equation is shown in Equation 7: θ = 1 − ln 1 + ψ C r ln 1 + 10 6 C r θ s ln e + ψ α ″ n ″ m ″ (7)

To reduce parameter complexity and adapt to the requirements of discrete element simulation, Yu (Yu D. et al., 2023) simplified Equation 7, resulting in Equation 8: θ = θ s ln e + ψ a b c (8)

In this context, θ represents the volumetric water content, while θ _s denotes the saturated volumetric water content. Parameter a is related to the inflow value (kPa), and b is the soil water content parameter when the water uptake exceeds the inflow value. Parameter c refers to the residual water content, and e is the natural constant, approximately equal to 2.71828, the term ψ represents the matric suction (kPa).

Based on Childs’ research findings (Childs and Collis-George, 1950), Fredlund obtained a function that reflects the relationship between the relative permeability k_r and the volumetric water content θ of the soil, as shown in Equation 9. k r = ∫ ln ⁡ ψ ln ⁡ 10 6 θ e y − θ ψ e y θ e y d y ∫ ln ∂ e v ln ⁡ 10 6 θ ′ e y − θ s e y θ ′ e y d y (9)

In this context, y is the integral variable of matric suction, representing the logarithm of matric suction.

As shown in Figures 3a,b, the mine is located in Yincheng Town, Dexing, Jiangxi Province, and belongs to the western part of the Circum-Pacific metal mineralization belt. It is a super-large open-pit mine. The mining area consists of low mountains and hilly steep slopes, with the ore body located in a non-karst area. The surrounding rock is mainly schist, which has a relatively dense structure. The near-surface rocks are severely weathered, with a weathering zone depth of 10–30 m, and a maximum depth of 50 m. The rock fractures are well-developed and loose. Additionally, the rocks in the structural zone and its influence zone are more fragmented, have poor stability, and are prone to roof falls and collapses.

Satellite imagery in Figure 3a shows that the terrain in this area is higher in the northeast and lower in the southwest. The northern part is steep, while the southwestern part has only a narrow valley outlet, making the terrain conducive to water accumulation. The climate in the mining area is warm and humid, with an annual precipitation of up to 2,185 mm. Rainfall is concentrated in the summer months of June to August, and rainwater tends to accumulate in the study area, forming depressions. The aquifers in the mining area are mainly Quaternary porous aquifers, weathered rock, and structural fracture water-bearing zones. The fractured characteristics of the rock significantly enhance its permeability, leading to damp underground tunnels and common dripping phenomena.

Figure 3c shows the area of the mining site studied using SBAS-InSAR technology. The mining site is nearly circular, with a diameter of about 1,000 m and a maximum depth of approximately 350 m. Most of the stepped slope angles are between 37° and 53°. Additionally, the area outlined by the red dashed line in Figure 3c is the research area for subsequent numerical simulations.

Based on the regional topographical features, a total of 267 Sentinel-1A images from June 12, 2015, to September 16, 2024, were selected for analysis. The images were acquired in the IW (Interferometric Wide Swath) mode, using VV polarization and operating in the C-band, with an incidence angle of approximately 39.09°. To enhance data processing accuracy, precise orbital data were collected, and a 30 m × 30 m SRTM1 DEM dataset was utilized to eliminate interference from flat phase and topographic effects. Figure 4 shows the common area covered by the 267 images.

Figure 5 illustrates the time baseline and spatial baseline distribution. From the time baseline perspective, the images show a relatively uniform distribution of time intervals, indicating good temporal coverage, which effectively supports the analysis of long time series. Regarding the spatial baseline distribution, the Average Spatial Coherence of the image pairs ranges from 0.3 to 0.8, reflecting a high degree of closeness and stability in the matching relationships among the image set.

As shown in Figure 6, the coherence coefficient (γ) is calculated based on the interferometric phase information from multi-temporal SAR images, with coherence values ranging from 0.45 to 0.96, exhibiting a clear spatial gradient characteristic. High coherence areas (γ > 0.75) are mostly distributed on relatively stable surfaces in the mining area, where the quality of the interferometric data is excellent and can serve as reference points for ground subsidence monitoring. Moderate coherence areas (0.6≤γ ≤ 0.75) are located in the central and lower parts of the mining site. Due to the influence of the monitoring angle of the terrain, there is a moderate decrease in coherence, but it still possesses certain potential for interferometric processing. Low coherence areas (0.45≤γ ≤ 0.6) are found at the bottom of the mining site, specifically in the pit bottom area. Affected by complex terrain conditions, the data noise is relatively high, and it needs to be eliminated in subsequent processing to improve the accuracy of deformation monitoring.

As shown in Figures 7–9, the surface settlement in the study area was monitored from 2015 to 2024 using SBAS-InSAR technology, and the maps of cumulative settlement change characteristics along the line-of-sight (LOS) direction were obtained. The monitoring results show that the quarry area exhibits significant spatial and temporal differential deformation characteristics during this period. Negative values indicate subsidence deformation away from the satellite and positive values indicate uplift deformation towards the satellite.

The deformation of the mining area is primarily characterized by subsidence, with the amount of subsidence accumulating and intensifying year by year. The most significant subsidence is observed in the steep high slopes in the eastern part of the mining area, with a maximum subsidence of −168 mm, indicating a strong impact of mining activities on surface deformation. Figures 7a–i show that from 2015 to early 2017, the range of deformation in the mining area was relatively small, mainly concentrated in the bottom of the mining pit, with a maximum subsidence of about −72 mm. Starting in 2018, the range of deformation gradually expanded, and the overall subsidence significantly increased, with the maximum subsidence reaching −110 mm in 2018, and the subsidence area spreading from the bottom of the mining pit to the surrounding areas, as shown in Figures 7j–l. Figures 8d–l indicate that from 2020 to 2022, significant subsidence deformation occurred in various stepped areas of the mining site, with maximum subsidence exceeding −130 mm, and the area of subsidence further expanded. In 2023 and 2024, the deformation in the mining area reached its peak during the monitoring period, at −168 mm, primarily distributed in the steep high slope areas of the mining site, as shown in Figures 9a–f.

The analysis results based on the SBAS-InSAR method reveal the patterns of subsidence deformation and spatial distribution of the deformation areas. This provides valuable data on subsidence regions and deformation characteristics, which can support the study of rock layer instability patterns and moisture content evolution during the transition from open-pit to underground mining under heavy rainfall conditions.

The heterogeneity, anisotropy, and discontinuity of rock masses result in their strength being a complex combination of block strength and the strength of structural surfaces (Seshagiri Rao, 2020). The presence of fractures, joints, and faults in the rock layers requires a weakening treatment of the rock mechanics parameters obtained from laboratory tests (Table 2). Based on the Geological Strength Index (GSI) and the generalized Hoek-Brown equation, the reduced rock mechanics parameters are presented in Table 3.

MatDEM can transform the macroscopic mechanical properties of rock into the microscopic mechanical parameters of stacked particles, which can make the discrete element stacked materials have specific elastic modulus and strength properties (Yu H. et al., 2023). Using MatDEM built-in macro and micro mechanical properties conversion formula (Liu et al., 2013), the microscopic rock mechanical parameters are obtained as shown in Table 4.

As shown in Figure 10, the discrete element model is composed of a large number of particles, and the diameter and quantity of the particles affect the experimental results. Generally speaking, the smaller the particle diameter and the larger the quantity in the model, the higher the credibility of the test results, but this also places greater demands on computer performance. This numerical model consists of 800,000 particles, with a dispersion coefficient of 0.2 and a particle size distribution of 0.017–0.024 m. The boundary conditions of the model are set to impermeable, and the computer took a total of 96 h to run this model. Monitoring points for volumetric moisture content are established on the longitudinal, cross-sectional, and layer planes, with naming conventions as follows: longitudinal (L1-L5), cross-sectional (S1-S5), and layer (B1-B5). Table 5 presents the key parameters for conducting numerical simulation experiments.

In this experiment, the majority of the rock mass displacement is in the form of settlement deformation. To observe this intuitively, it is defined that the settlement amount is a positive value when the rock mass settles downward. The primary causes of rock deformation include: (1) Mining factors of chamber and pillar. Continuous mining causes damage to mining face and overlying rock structure. (2) Volume moisture content change factors. During the processes of seepage and erosion, the mechanical properties of the rock are weakened, leading to the formation of numerous new cracks and the development of fissures, which create new seepage pathways, thereby enhancing the weakening effect of water on the geomechanical properties (Liu et al., 2024; Wang et al., 2024).

Rock shape variable refers to the displacement difference before and after the end of mining, which represents the displacement caused by mining, and the direction is vertical downward. In Figures 11a, b, after the completion of mining in Chamber 1#, the deformation in the mining area is minimal, with a maximum deformation of 0.025 m. The slopes remain stable, with no occurrence of particle sliding or falling. In Figure 11c, after the mining of Chamber 2#, the deformation in the Chamber 1# area increases, with a maximum deformation of 0.032 m. This indicates that rock layer damage has a transmissive nature. The formation of mining Chamber 2# creates a mining area, and the surrounding rock layers need to bear more overburden weight, resulting in the continued increase of deformation in the mined-out area of Chamber 1# due to the impact of subsequent mining. As mining progresses to Chamber 3# as shown in Figure 11d, a “triangular” rock movement zone develops upward, with Chambers 1#, 2#, and 3# forming the base. In this area, deformation increases from top to bottom, reaching a maximum deformation of 0.045 m. During this stage, deformation occurs in the slopes of the mining area. Analysis using the controlled variable method indicates that the mining of Chambers 1#, 2#, and 3# contributes to surface deformation due to the transmissive nature of rock destruction. As shown in Figures 11e,f, after the mining of Chambers 4# and 5#, significant changes occur in the overburden and surface. The “triangular” rock movement zone expands further, connecting upward to the surface and downward to the mined areas and pillars. The stepwise deformation of the slopes increases gradually, with particle sliding occurring at the top of the steps.

During the pillar mining phase, significant changes occur in the overburden and surface, with increasing levels of damage becoming evident. In Figure 11g, after the completion of mining in Pillar 1#, the external shape of the rock layers exhibits minimal deformation changes, although new fractures and voids may have developed internally. Currently, these fractures and voids are insufficient to cause significant deformation of the rock layers. In Figure 11h, after the mining of Pillar 2#, the “triangular” rock movement zone experiences subsidence. The area of rock collapse covers Chambers 1#, 2#, and Pillar 1#, with a maximum deformation of 0.135 m. There is a notable subsidence of the slope steps, with the maximum deformation occurring at the top step. In Figures 11i,j, as the ongoing recovery of the pillars continues, the degree of damage to the overburden and surface increases. After the recovery of Pillar 4#, the “triangular” rock movement zone has spread horizontally to include Chambers 3#, 4#, and Pillar 3#, while vertical expansion has reached the center step, with a maximum deformation of 0.16 m. In fact, the “triangular” rock movement area is essentially a shear slip zone caused by the failure of the stress arch above the mined-out area (Zhang et al., 2021): after the extraction of the coal pillar, the stress redistribution of the overburden forms an arch structure. However, continuous mining causes the stress at the foot of the arch to exceed the limit, leading to the penetration of vertical cracks and horizontal layering, ultimately resulting in triangular shear slip and causing damage to the overburden and the ground surface.

The process of overburden collapse can be divided into three phases: (1) Chamber Mining Phase: In this stage, minor rock falls occur in the overburden mining area, and slight subsidence is observed at various slope levels on the surface. Overall, the degree of deformation in the overburden and surface is minimal, categorizing this period as the “latent deformation phase.” (2) Early Pillar Mining Phase: After the completion of mining in Pillar 1# and Pillar 2#, a “triangular” rock movement zone forms, and significant surface subsidence becomes evident. This is identified as the “explicit deformation phase.” It represents a transitional period where the lack of pillar support leads to considerable damage to internal fractures and voids within the rock layers, causing overburden collapse. (3) Late Pillar Mining Phase: After the completion of mining in Pillar 3# and Pillar 4#, the “triangular” rock movement zone continues to expand, with increasing deformation. Noticeable collapses of the overburden occur, and surface subsidence intensifies, leading to the classification of this stage as the “severe deformation phase.”

Figure 12 illustrates the displacement variation curves for the X-axis monitoring points (L1, L2, L3, L4, L5). The displacement in the mined area is relatively large, with displacements during the chamber mining phase being smaller than those during the pillar recovery phase. Displacements at L1 and L5 are minimal, primarily due to their greater distance from the mining chambers and pillars, resulting in reduced exposure to mining-related effects. In contrast, monitoring points L2, L3, and L4, located directly above the mining area, exhibit significant fluctuations in displacement, clearly affected by mining activities. After the mining of all chambers, the displacements at L2, L3, and L4 are 0.044 m, 0.028 m, and 0.025 m, respectively, with average displacement rates of 0.0088 m/stage, 0.0056 m/stage, and 0.005 m/stage for each chamber. Following the pillar mining, the displacements at L2, L3, and L4 increase to 0.095 m, 0.093 m, and 0.083 m, with average displacement rates of 0.0127 m/stage, 0.0162 m/stage, and 0.0145 m/stage for each pillar. The red dashed circle indicates that point L1 experiences a subsidence phenomenon during the mining of Pillars 3 and 4. This is likely influenced by mining activities, where regional stresses lead to localized uplift.

As shown in Figure 13, the displacement magnitudes and changes at each monitoring point are relatively consistent, demonstrating synchronous variations. This is primarily due to the alignment of the monitoring points parallel to Chamber 3#, which experiences similar impacts from the different mining phases. The displacement changes indicate that the displacements during the chamber mining phase are smaller than those during the pillar recovery phase.

Figure 14 illustrates that monitoring points S1, S2, S3, and S4 are arranged from the surface downward into the overburden. The observed pattern reveals that the displacements at monitoring points deeper in the rock layers are greater than those at the surface monitoring points.

Figure 15 presents the cloud diagrams of volumetric water content variations in the longitudinal, horizontal, and cross-sectional views. The stability of the rock layers is closely linked to changes in volumetric water content. Sudden changes in water content indicate that the internal fractures and voids within the rock layers have accumulated to a significant extent, compromising the internal structure’s ability to support the overburden, leading to a tendency for collapse. As previously analyzed, after the completion of mining in Chamber 3#, a “triangular” rock movement zone forms, with the volumetric water content in the mining area at 21.5%. Following the completion of mining in Chamber 5#, the “triangular” rock movement zone expands to the external slopes, resulting in a volumetric water content of 24.6% in the affected area. After the mining of Pillar 2#, the “triangular” rock movement zone experiences collapse, leading to significant structural damage and an increase in volumetric water content to 40.3%. Finally, after the completion of mining in Pillar 4#, a widespread collapse occurs in the “triangular” rock movement zone, with pronounced subsidence in the slopes, resulting in a volumetric water content of 44.8% within the mining area.

This is a revelation: the sudden change in volumetric water content is mainly due to the formation of fissures and subsidence cracks caused by underground mining, which accelerates water seepage and transport, firstly leading to local rock layer instability, and the unstable area can promote water seepage and transport. This is of guiding significance in preventing the formation of subsidence disasters. With the advancement of mining, the rock layer as a whole remains stable, but the deformation of the rock layer in the local area, which is often difficult to observe from the surface, can be monitored through the change of volumetric water content to monitor the deformation and the generation of unstable areas.

Figures 16a, b are scatter plots of volumetric water content θ against matric potential ψ and permeability coefficient k_r , respectively, calculated using Equations 8, 9. The red curves represent the corresponding fitted curves. Equations 10, 11 are fitting functions. Based on these functions, the matrix suction and permeability coefficient can be obtained from the volumetric water content measured on-site, providing a scientific basis for safe production.

Fitting formula: Ψ = 46.1 − 4.6 θ + 0.21 θ 2 − 0.0052 θ 3 + 0.00064 θ 4 (10) K r = 3.54 × 10 − 7 − 4.78 × 10 − 9 θ + 1.51 × 10 − 9 θ 2 − 3.82 × 10 − 12 θ 3 (11)

As shown in Figure 16a, when the volumetric water content is at 20%, the rock mass is in a dry state, and its matric suction remains at a high level, approximately 11.2 kPa, indicating good structural stability of the rock mass. As mining activities progress, the volumetric water content of the rock layer shows a significant increasing trend; for example, the water content of mining Chamber 3# at the end of mining is 21.5%, while that of mining Pillar 4# at the end of recovery is 44.8%. This corresponds to a significant decrease in matric suction from 8.2 kPa to 3 kPa. The essence of matric suction is the capillary tension of the pore liquid within the rock mass (Schmitt et al., 1994). When its value drops below the critical threshold of 3 kPa, it indicates that the capillary framework support within the rock mass is essentially lost, leading to a significant weakening of the rock mass’s shear strength. The interaction between volumetric water content and permeability coefficient exhibits completely different characteristics, as shown in Figure 16b. When the volumetric water content increases from 20% to 60%, the permeability coefficient surges from an initial value of 6.15 × 10⁻¹⁰ m/s to 3.71 × 10⁻⁶ m/s, indicating the expansion effect of increased water content on the flow pathways (Simon and Collison, 2001; Wu et al., 2021).

In summary, as the water content continues to increase, the suction of the rock mass matrix shows a significant downward trend, which directly leads to a decrease in the shear strength of the rock mass, resulting in displacement or even collapse of the overlying rock. It is noteworthy that the decay of matrix suction not only weakens the mechanical properties of the rock mass but also promotes the expansion and development of the fracture network, a process that further increases the permeability coefficient of the rock mass. The increase in the permeability coefficient accelerates the movement of moisture within the rock mass, and the intensified movement of moisture in turn leads to a continuous increase in the volumetric water content of the rock mass, thereby forming a mutually reinforcing vicious cycle.