Skip to main content

ORIGINAL RESEARCH article

Front. Environ. Sci., 06 December 2021
Sec. Environmental Informatics and Remote Sensing
Volume 9 - 2021 | https://doi.org/10.3389/fenvs.2021.787053

Study of Ground Movement in a Mining Area with Geological Faults Using FDM Analysis and a Stacking InSAR Method

www.frontiersin.orgZhengyuan Qin1* www.frontiersin.orgVivek Agarwal1 www.frontiersin.orgDavid Gee1,2 www.frontiersin.orgStuart Marsh1 www.frontiersin.orgStephen Grebby1 www.frontiersin.orgYong Chen3* www.frontiersin.orgNingkang Meng3
  • 1Nottingham Geospatial Institute, Faculty of Engineering, University of Nottingham, Nottingham, United Kingdom
  • 2Terra Motion Limited, Nottingham, United Kingdom
  • 3State Key Laboratory of Coal Resources and Safe Mining, China University of Mining and Technology, Xuzhou, China

Underground coal mining activities and ground movement are directly correlated, and coal mining-induced ground movement can cause damage to property and resources, thus its monitoring is essential for the safety and economics of a city. Fangezhuang coal mine is one of the largest coalfields in operation in Tangshan, China. The enormous amount of coal extraction has resulted in significant ground movement over the years. These phenomena have produced severe damages to the local infrastructure. This paper uses the finite difference method (FDM) 3D model and the stacking interferometric synthetic aperture radar (InSAR) method to monitor the ground movement in Fangezhuang coalfield during 2016. The FDM 3D model used calibrated Fangezhuang geological parameters and the satellite InSAR analysis involved the use of ascending C-band Sentinel-1A interferometric wide (IW) data for 2016. The results show that the most prominent subsidence signal occurs in mining panel 2553N and the area between panel 2553N and fault F0 with subsidence up to 57 cm. The subsidence observed for the FDM 3D model and stacking InSAR to monitor land deformation under the influence of fault are in close agreement and were verified using a two-sample t-test. It was observed that the maximum subsidence point shifted towards the fault location from the centre of the mining panel. The tectonic fault F0 was found to be reactivated by the coal mining and controls the spatial extent of the observed ground movement. The impact of dominant geological faults on local subsidence boundaries is investigated in details. It is concluded that ground movement in the study area was mainly induced by mining activities, with its spatial pattern being controlled by geological faults. These results highlight that the two methods are capable of measuring mining induced ground movement in fault dominated areas. The study will improve the understanding of subsidence control, and aid in developing preventive measures in Fangezhuang coalfield with fault reactivation.

1 Introduction

As China’s primary source of energy, coal it is critical to the country’s social and economic growth. Large-scale coal mining, however, has the potential to substantially harm the environment in mining areas, as well as causing a variety of geological issues and societal problems (Diao et al., 2019). As one of the main effects, a thorough understanding of mining-induced subsidence is essential for preventing or mitigating such issues.

With regards to measuring surface movement, traditional methods like levelling, Global Navigation Satellite System (GNSS), 3D laser scanning, and similar provide accurate information, but can be time-consuming, expensive, and have low spatial extent and are therefore unsuited to surveying large (basin level) areas. A better alternative is to use the satellite Interferometric Synthetic Aperture Radar (InSAR) technique, which provides weather independence, sunlight independence (active sensor), high (basin-level or greater) spatial coverage, and is less tedious and more economical. Therefore, InSAR offers a spatial resolution comparable to other traditional methods of land surveying. Many researchers have used InSAR to study land movement resulting from earthquakes (Fialko, Simons, and Agnew 2001), glacial movements (Wang et al., 2015), landslides (Sun et al., 2015; Kang et al., 2017; Zhang et al., 2020), volcanic bulging (Fournier et al., 2010; Albino et al., 2020), groundwater extraction (Bell et al., 2008; Motagh et al., 2017; Castellazzi et al., 2018) and coal mining (Zhang et al., 2015; Yang et al., 2017; Diao et al., 2018; Dong et al., 2021).

Extensive research has been carried out to study ground movement induced by underground mining activities (Park et al., 2012; Yang et al., 2019). In recent years, several research tools, such as numerical analysis (Deck and Anirudh 2010; Shabanimashcool and Li 2012; Shi M et al., 2021) and InSAR time-series analysis (Sowter et al., 2016; Gee et al., 2017; Grebby et al., 2019; Ghayournajarkar and Fukushima 2020) have been widely employed to analyse the mechanism of ground movement. Furthermore, several studies have also focused on ground movement associated with fault activation. Bell et al. (2005) reported mining subsidence from Great Britain, Germany and Colombia, and stated that mining area experiencing reactivation of faults should be surveyed properly before construction and a safe gap of at least 10 m should be maintained between the fault zone’s edge and any structures. Moreover, Mohammady et al. (2019) employed Random Forest theory to analyse subsidence susceptibility and found that gap from the fault, elevation, slope angle and water table had the largest influence on ground deformation. Gumilar et al. (2015) and Pacheco-Martínez et al. (2013) discovered a direct link between fault position and high land deformation rates.

In terms of numerical analysis, the finite difference method (FDM) is commonly utilised within the FLAC3D software because of its efficacy as a tool for solving rock mechanics and geo-mechanical issues. It can handle material heterogeneity, nonlinearity, complicated boundary conditions, ground condition pressures and gravity. The FDM model idealises the rock mass as a continuous medium that deforms according to a given constitutive law, satisfying compatibility and equilibrium criteria, and yields approximate partial differential equation solutions. Due to the discontinuities (e.g., geological fault) in the study area, a discontinuous numerical method was chose to study the ground movement based on the FDM analysis with FLAC3D software. In FLAC3D, the ground movement value of any point in the model can be monitored; it is typically more applicable for continuous and uniform subsidence. FDM analysis is widely used to study the ground movement in coal fields (Cheng et al., 2019; Parmar et al., 2019; Sikora and Wesołowski 2021; Yan et al., 2021).

Only a few attempts have been made to study 3D model simulation of ground movement in conjunction with InSAR- derived ground movement to understand the effect of a geological fault (Jeanne et al., 2019; Perry et al., 2020; Francioni et al., 2021; Shi Y et al., 2021). Although these previous studies suggest that faults do affect the ground movement, the nature and magnitude of effects have not been fully explored. The current approach uses the hybrid method, which is the combination of the numerical method and stacking InSAR analysis to study how the fault affects the ground movement induced by mining. To study the ground movement associated with subsurface mineral extraction at Fangezhuang coal mine, a fully 3-D elastoplastic FDM model was constructed. The results of the developed FDM model were verified by comparing the outputs of the model with the ground movement revealed by stacking InSAR analysis. Finally, the ground movement under the influence of fault in Fangezhuang coal mine is discussed.

This paper is divided in five sections. The section 1 provides an Introduction, and the section 2 describes the study area. Section 3 describes the methodology for studying the ground movement with the 3D FDM model and stacking InSAR technique. In section 4, we reveal the ground movement in Fangezhuang coal mine using FDM model and InSAR technique and perform a comparison between FDM model and InSAR result and the ground movement trend in Fangezhuang coal mine is revealed. Finally, conclusions and future scope of the study are discussed in section 5.

2 Study Area

Fangezhuang coal mine lies in the Kailuan area, Hebei Province, China, approximately 23 km east of Tangshan city. While most coal lies in plain, the ground surface consists of farmland, villages, mining factories and ground sunk fissure caused by mining subsidence. Fangezhuang coal mine has an area of 32.33 km2 in which the north-south direction is 12.25 km, and the highest east-west inclination is 3.92 km. The coal seam is nearly horizontal, with a mean thickness of 3 m. This region is characterised by large ground movements, with thick loose layers and fast advancement of the working panel. The study area is marked in Figure 1.

FIGURE 1
www.frontiersin.org

FIGURE 1. Location of the study area. Blue spot is Fangezhuang coal mine in Tangshan city. Purple rectangle area is mining panel 2553N. Red line shows the F0 geological fault.

2.1 Geological Structure

Within the Fangezhuang coal mine, one major fault exists: F0 (Figure 1). F0 fault is a normal fault that runs through the Fangezhuang coal mine in a north-south direction. According to the drilling data, the fault extends about 3,000 m, and the drop ranges from 14 m to 37 m. The fault is a high-angle normal fault, inclined to SWW, with a dip angle of 70°–84°. F0 fault runs through the whole monocline structural area and has a significant impact on mining production. Meanwhile, there are associated faults of a specific scale with a significant drop on both sides of the F0 fault, and most of the associated faults here develop along the strike of the F0 fault.

2.2 Coal Seam and Rock Strata

Figure 2 shows the cross-section along H-H’. The strata revealed through exploration of Fangezhuang coal mine are Ordovician (O), Carboniferous (C), Permian (P) and Quaternary (Q), as described below.

1) The O stratum, 500–900 m thick, comprises fine-grained dolomite, stratified limestone.

2) The C stratum, 140–290 m thick, comprises siltstone, mudstone and fine sandstone. This formation also contains 1-3 layers of unstable thin coal seams.

3) The P stratum has a mean thickness of 550 m. From bottom to top, the sequence is split into the four forms showed below.

• The Damiaozhuang formation mainly contains siltstone, mudstone and medium sandstone. Four layers of the coal seam can be mined, namely No. 5 coal seam, No. 7 coal seam, No. 8 coal seam and No. 9 coal seam.

• The Tangjiazhuang formation is mainly composed of coarse-to-medium sandstone, followed by fine sandstone. The lower strata are interbedded with 1–4 layers of thin coal lines.

• The Guye formation mainly comprises of medium-coarse sandstone with a small amount of mudstone and siltstone.

• The Wali formation is mainly composed of medium-coarse sandstone, fine sandstone and siltstone. At the bottom, a layer of aluminum mudstone with a thickness of about 4–5 m is developed.

4) The Q stratum is mainly composed of clay, sand and gravel layer. The Quaternary alluvium covers the whole Fangezhuang mine field. The thickness of alluvium varies from 54 m to 424 m and gradually thickens from north to south.

FIGURE 2
www.frontiersin.org

FIGURE 2. Cross section view along H-H’ (marked in Figure 1).

In this research, No. 5 coal seam, located at around 560 m deep in Permian strata, is the main workable seam. The mining panel 2553N which started in January 2016 and terminated in December 2016 along the northwest-southeast direction, is chosen to analyse ground movement. For the geological structure, F0 fault cut No. 5 coal seam along the southeast direction, resulting in a decrease of the depth of No. 5 coal seam at the east side of F0.

3 Methodology and Data Used

In this study, FLAC3D was adopted for numerically predicting ground movement caused due by coal mining. FLAC3D is a 3-D finite-difference computer program to solve geological problems (Kumar et al., 2016; Shi M et al., 2021). In FLAC3D, the initial model is the geomechanical model with the engineering scale, which is able to simulate the real geological condition in the mining coalfield. The null model plugged in FLAC3D is used to delete the elements to simulate the mining extraction activities. When the elements representing the coal are deleted to leave a void, the overlying elements will have a free displacement boundary and the stress field will redistribute, resulting in the overlying elements caving. The caving is spread upwards to the top surface, which represent the ground movement. The InSAR analysis used 25 Sentinel-1 ascending, Interferometric Wide (IW) Single Look Complex (SLC) images between January 2, 2016 and December 21, 2016. The C-Band (wavelength of 5.6 cm) data have a pixel spacing of 2.3 m in range and 13.9 m in azimuth, and a revisit period of 12 days (Torres et al., 2012).

3.1 Establishment of the 3D FDM Model

The mining panel 2553N shown in Figure 1 is selected as the computational model domain on the basis of spatial distribution of the considered mining area. Based on the geological properties of the Fangezhuang coalfield, a 3D FDM simulation model of panel 2553N and F0 fault is set up, which is shown in Figure 3. The parameters of the FDM model are set as follows: the size of panel 2553N is 900 m × 180 m, No.5 coal thickness is 3 m, mining depth is 560 m, angle of the fault F0 is 70°, and step excavation distance is 75 m. The coal mineral extraction was simulated as long-wall mining with step by step excavation. The ground surface deformation map was obtained for the study area from January 2016 to December 2016 when the ground surface was stable (where the unbalanced force is less than 10–5 of the maximum unbalanced force) (Du et al., 2019). The dimensions of the model are 1800 m long, 900 m wide, and 590 m high. The top surface does not show any geomorphic feature and is assumed horizontal in the model. The shortest horizontal distance between the mining boundary and the model boundary is 360 m to remove the boundary effect.

FIGURE 3
www.frontiersin.org

FIGURE 3. FLAC3D FDM model.

3.1.1 Generalisation of Strata and Faults

Mining panel 2553N lies in the Permian strata, and No. 5 coal seam is the main workable coal seam. As a result, two geological strata Quaternary and Permian, as well as No. 5 coal seam are simulated in a FDM model. The structure of each strata is made with reference to the geological report of Fangezhuang coal mine and borehole data. The fault could be designed with interface elements in FLAC3D if the thickness of the fault is small (Cai et al., 2021). However, when the thickness of a fault is more than 10 m, it is appropriate to use a layer with a certain thickness to simulate a fault (Xu et al., 2013). In this paper, we choose a layer with two boundary surfaces to model the F0 fault around 10–20 m thick.

3.1.2 Computational Mesh

FLAC3D supports numerous element shapes, including hexahedrons, tetrahedrons, pyramids, triangular prisms, and others (Abbasi et al., 2013). With a consideration of achieving a balance between processing time and simulation accuracy, a set of octahedral elements was adopted to mesh the simulation model. The meshing approach is as follows. It can be seen in Figure 3, from top to bottom, the model is divided into 13 layers which are loose layer, fine sandstone, siltstone, medium coarse sandstone, mudstone, medium coarse sandstone, siltstone, fine sandstone, siltstone, No. 5 coal, siltstone, No. 7 coal and siltstone. Each layer has a thickness of 3–160 m. The closer a layer is to the coal seam roof, the thinner the layer is set; the farther it is from the roof, the thicker it is set.

3.1.3 The Constitutive Model and Boundary Conditions

The Mohr-Coulomb model is used as the constitutive model and Mohr–Coulomb yield criterion are used in this study, along with displacement boundary conditions. No vertical displacement is allowed at the bottom of the simulation model, and no displacement is allowed in the direction perpendicular to the lateral boundaries. Trapezoidal distributed load is applied in the horizontal direction. Gravity force is applied to the model. A numerical model is obtained from these results, which reflects the mining and geological condition of the region after long-wall panel exploitation.

3.1.4 Monitor Points

The 3D FDM model simulates the whole No. 5 coal extraction of panel 2553N which commenced in January 2016 and stopped in December 2016. Based on the mesh of the FDM model, on the top boundary of the model we set up 224 monitoring points (16 × 14) to obtain the ground movement during the coal mining process.

3.1.5 Mechanical Parameters of Rock Strata

Intact rock and discontinuities make up rock layers. The mechanical properties of complete rock samples acquired through laboratory testing differ substantially from the mechanical properties of rock strata (Hoek and Brown 1997). However, to some extent, the reliability of numerical simulation results depends on the choice of rock mechanical parameters. The methods for determining mechanical parameters mainly include empirical reduction method, engineering rock mass classification method and displacement back analysis method. This research adopts the back analysis with orthogonal test and numerical simulation to determine the rock mechanical parameters.

The detailed procedure is as follows (Xu et al., 2013):

1) Based on the geological report of Fangezhuang coal mine, the initial rock parameters are shown in Table 2 (Ren 2017). The average value of the deformation modulus  E¯, Poisson’s ratio μ¯, cohesion  C¯ and internal friction angle ¯ are chosen as the experimental factors. The initial mechanical parameter divided by the corresponding average value is calculated for each mechanical parameter of each strata. These values are listed as KEi ,  Kµi , KCi, and KØi in Table 1. Therefore, each strata mechanical parameters could be given as KEiE¯ ,  Kµiµ¯ , KCiC¯ and KØiد. The maximum ground movement Wmax is chosen as the test indicator.

2) The orthogonal test with four factors and five levels L25 (54) was set up to test the rock parameters (Table 2). All the testing schemes are listed in an orthogonal table (Table 3).

3) 25 numerical simulations were conducted using FLAC3D and the results for each scheme were listed in Table 3. According to the result in Table 3, it was found that the maximum ground movement of the 4th scheme is closest to the measured ground movement of 16.2 cm (geological report of Fangezhuang coal mine). Therefore, the parameters of the four factors are  E¯=3.6GPa, μ¯=0.31,  C¯=6.3MPa and ¯=37°. The final parameters for each layer shown in Table 4.

TABLE 1
www.frontiersin.org

TABLE 1. Initial mechanical parameters of rock formations.

TABLE 2
www.frontiersin.org

TABLE 2. Levels of test factors.

TABLE 3
www.frontiersin.org

TABLE 3. Orthogonal test design and results.

TABLE 4
www.frontiersin.org

TABLE 4. Calibrated mechanical parameters used in 3D model.

3.2 Stacking InSAR Analysis

An InSAR stacking method was employed to map the ground deformation due to the decimetre scale rates of deformation expected to occur following the long- wall extraction of panel 2553N. The processing chain utilised is summarized in Figure 4. The Sentinel-1 SLC data were initially deburst and merged before co-registration to the slant-range coordinate system of the master image (January 2, 2016). Phase ramps attributed to orbital errors were subtracted using the precise orbit determination and topographic phase was simulated and removed using a DLR digital elevation model (DEM) from the TanDEM-X mission (Rizzoli et al., 2017). The data were multi-looked by a factor of 4 in range and 1 in azimuth and interferograms, with approximately a 10 m resolution, were generated between consecutive SAR acquisitions irrespective of the perpendicular baseline. The interferometric fringes (phase cycles that correspond to a displacement of half of the sensor wavelength) are related to the surface deformation. Generating interferograms only over the shortest epochs the minimises like likelihood of quantitatively underestimating the deformation due to the ambiguity effect when deformation gradients are high, such as over active mining sites. Further, an accurate estimation of the deformation depends upon high coherence (or phase correlation) between the two forming SAR images, hence, utilizing consecutive image pairs helps maintain coherence. In addition, the interferograms were filtered using a modified Goldstein filter to further improve coherence and the quality of phase. The interferograms were unwrapped from modulo-2π phase to relative deformation using a statistical-cost network-flow algorithm (Chen and Zebker 2001) with respect to a reference point located at Tangshan (39.6309° N, 118.1802° E). The unwrapped interferograms were subsequently stacked and an average rate of motion was derived from a least squares covariance analysis of the unwrapped phase. Once linear velocities had been generated, the relative height change for each image acquisition was calculated in accordance with that of Berardino et al. (2002). Finally, the line-of-sight time-series were projected into the vertical, by means of dividing by the cosine of the incidence angle (∼ 0.639 radians), to facilitate an appropriate comparison with the model.

FIGURE 4
www.frontiersin.org

FIGURE 4. Stacking InSAR processing chain used for the study.

4 Results and Discussion

4.1 Comparison Between FDM Model and Stacking InSAR Result

Figure 5 shows the vertical surface deformation over the study period as generated by (a) InSAR stacking, (b) 3D FDM model and (c) Residual of the InSAR movement minus FDM model. Surface uplift is represented by a positive value and subsidence by a negative value. Both the InSAR and 3D FDM model results show the severe subsidence pattern nearby the mining panel 2553N. For the FDM model, the ground movement ranges from 34.2 to 57.6 cm, while for the stacking InSAR result, the vertical deformation ranges from 15.2 to 58.3 cm. Based on the mining ground movement theory, if the ground movement is greater than ±10 mm/ year then the area is considered unstable (Zheng et al., 2018). From the perspective of ground movement, both the FDM model and InSAR measurements confirm that the area is unstable due to the coal mineral extractions. The ground movement basin is located within the mining panel 2553N and the area between the panel and fault F0.

FIGURE 5
www.frontiersin.org

FIGURE 5. Three-dimensional deformation generated by (A) InSAR result (B) 3D FDM model result (C) Residual result (InSAR—3D FDM model).

It can be seen from Figure 5A that the shape of the subsidence basin is not symmetrical over the mining area (mining panel 2553N), and the maximum subsidence point is located at point B, which is about 75 m from point A (centre of panel 2553N). The subsidence map shows a clear boundary nearby Fault F0. The boundary corresponds to the fault F0 and clear ground movement differences can be seen either side of fault F0. Figure 5B shows the similar subsidence basin revealed by FDM modelling. The maximum subsidence point is also located at point B but with a marginally smaller magnitude compared with InSAR. The relative severe ground movement with a maximum value of 58 cm is located in the northeast part of panel 2553N. In the east of the fault F0, the maximum ground movement is approximately 15 cm. The east of panel 2553N presents relatively uniform and homogeneous subsidence, with a magnitude of up to 20 cm. Figure 5C shows the residual map of the FDM model and the InSAR result. The red colour represents where the ground movement of FDM model is less than the InSAR result. In other words, negative residual measurements represent underestimation of the FDM model with respect to the InSAR result. It can be seen that the average ground movement obtained by InSAR is smaller, especially in the area of severe ground movement between panel 2553N and fault F0. Nonetheless, the surface subsidence caused by long-wall mining is significant. InSAR only detects the medium-scale ground movement. As a result, the 3D FDM model proposed in this study is a valuable tool for ground movement analysis.

Since the numerical model does not take the rock fissure and other geological properties of rock mass into account, the simulated results appear relatively smoother than the InSAR result. In addition, the InSAR measurements are relative to a reference point and are therefore not absolute. Thus, when comparing the InSAR results directly to the results from the FDM-model, it is worth noting that there may be some offset if the InSAR reference point was not entirely stable (i.e., it has a non-zero displacement) during the observed period. Furthermore, InSAR results may be affected by noise such as atmospheric, ionospheric (Liao et al., 2018) and unwrapping correlation error (Yunjun et al., 2019). Thus, any differences between the measurements could be due to several aspects relating to the InSAR processing, e.g., noise and the fact that the InSAR measurements are relative to a reference point that may be moving slightly.

The ground movement pattern may result from underground mining, fault reactivation and other factors. Figure 6 shows the time series surface deformation of the FDM model and InSAR analysis at points A, B, C and D. It can be seen that the maximum subsidence values gradually increase from 14th January to 20th September. After 20th September, both the FDM model and InSAR result show that the rate of increase of the maximum ground movement value decreases. Although the mining panel is advancing from September to December, the ground movement is increasing slowly at a relatively steady pace. According to mining subsidence theory, when the advancing distance reaches 1.2–1.4 times the average mining depth, the mining activity will reach full mining condition (Wang et al., 2020). This suggests that after 9 months of mineral extraction (panel 2553N advancing 675 m along the strike direction), the mineral extraction reached supercritical mining. From January to March, uplift is observed in the InSAR time series, which might be the noise effect in the stacking analysis. Compared with the ground movement value at point A, point B experiences less ground movement, confirming that the maximum ground movement is not located at the centre of panel 2553N. There is a clear difference in ground movement between C and D, which are located on different sides of fault F0.

FIGURE 6
www.frontiersin.org

FIGURE 6. Time series comparison for FDM model and InSAR subsidence and at point A, B, C and D.

Table 5 summarises the two-sample t-test results for points A, B, C and D (marked in Figure 5) to determine, quantitatively, whether the differences between the InSAR and FDM model series at these points are within the acceptable range (Agarwal et al., 2020). The mean difference between the InSAR and the FDM model is not the same at each of these points. As a result, it is necessary to investigate if the difference is significant. For the test, the Null Hypothesis (H0) is taken as the average difference is zero (µ = 0) and alternate Hypothesis (HA) as the average difference is not zero (µ ≠ 0).

TABLE 5
www.frontiersin.org

TABLE 5. Two-paired t-test for FDM model and InSAR result.

The t-value computed for the comparison of the InSAR and FDM modelling outputs at point A is 1.32. This means that these values are 1.32 standard deviations off of the mean. A t-value less than the critical t-value (2.069 for point A) is necessary to accept the null hypothesis. The t-value for point A is less than the crucial t-value, meaining that the null hypothesis is accepted, indicating that there is no significant difference between the InSAR result and the FDM model at point A. Another approach to confirm this is to look at the p-value (0.25) and alpha (0.05 for the 95 percent confidence limit), where a p-value greater than alpha means the null hypothesis is accepted. Acceptance of the null hypothesis in both cases demonstrates that the InSAR-derived subsidence exists, and that the FDM model is consistent at point A. From Table 5, similar results can be seen for points B, C, and D. As a result, we have sufficient evidence that InSAR-based subsidence is consistent and agrees well with FDM model values for all four points A, B, C, and D.

In order to examine the ground movement trends, the transect lines M–M′ along the coal seam strike direction, N–N′ along the coal seam dip direction and F-F′ vertically across the fault F0 were extracted (as shown in Figure 5). Figure 7 shows the land subsidence profiles of observation lines M–M′, N–N′ obtained from the FDM model and stacking InSAR analysis. It can be seen that the ground movement trends of the two lines are almost identical. The root-mean-square error along M–M′ in the dip direction is 0.113 m, with a maximum difference 0.141 m. The root-mean-square error along N–N′ in the strike direction is 0.105 m, with the largest difference of 0.210 m. The blue and orange lines show InSAR and FDM-model subsidence, respectively. The grey shading shows the standard deviation of the displacements observed at each location for both InSAR and FDM subsidence. The overlapping darker area for both the curves highlights that the deviation between the methods is small and within the standard deviation error limits. It can be seen that the subsidence curves along with M–M′ and N–N’ do not remain horizontal at the bottom of the basin but appear to be fluctuating somewhat.

FIGURE 7
www.frontiersin.org

FIGURE 7. Land subsidence profiles of FDM model and InSAR result along (A) M–M′ and (B) N-N′.

From the subsidence profile of N–N′ in Figure 7B, it can be seen that the InSAR result shows the ground movement in the southern part of N–N’ is slightly greater than that in the northern part. This is likely because the ground movement obtained by InSAR includes some residual subsidence caused by adjacent goaf (Fan et al., 2021). According to the geological report of Fangezhuang coal mine, mining along a panel to the south of panel 2553N stopped in 2015 and may cause the residual ground movement in 2016.

4.2 Influence of Fault F0 on Ground Movement

Figure 8 shows the land subsidence profiles of observation lines F–F’ obtained from the FDM model and stacking InSAR analysis. The grey shading shows the standard deviation of the displacements observed at each location for both InSAR and FDM model subsidence. The overlapping darker area for both the curves highlight that the deviation obtained from both the methods is small and within the standard deviation error limits.

FIGURE 8
www.frontiersin.org

FIGURE 8. Land subsidence profiles of FDM model and InSAR result along F-F′.

For the profile along F-F′, it can be observed that at the distance of around 155 m there is a distinct change in slope in the subsidence profile. Interestingly, the fault F0 also lies at this location, which appears to be controlling this subsidence pattern. At this location, the fault is attributable for subsidence of 100–300 mm/ year on its right, and a subsidence rate of 300–500 mm/ year on its left. It can also be seen from Figure 5 that fault F0 has a noticeable effect on the spatial distribution of ground movement. This is most likely due to differential vertical compaction of various thicknesses of compressible soil deposited on both sides of the faults, which is greater in the hanging wall in normal faults (Burbey 2002). This observed structural control of land subsidence causes differential subsidence rates on either side of the fault line, which could potentially cause damage to villages, trains, and other infrastructure.

Fault F0 crosses the Fangezhuang coal mine in a north-south direction. The fault trace divides the subsidence rates and the arrows in Figure 8 indicate the relative displacement of the fault. To the west of F0, along F-F’, the relative severe subsidence has a maximum value of 425 mm obtained for the FDM model result and 433 mm measured by InSAR. The east of F0 exhibits relatively stable subsidence, with a value up to 225 and 152 mm obtained with the FDM model and InSAR, respectively. In addition, Figure 5 shows the maximum subsidence point (point A) is not above the centre of panel 2553N, but it approaches F0. It can be seen that fault F0 has a blocking effect on overburden movement. The surface subsidence morphology under the influence of fault is not symmetric about the centre of goaf. The maximum subsidence is shifted towards the fault. The ground movement patterns show a good agreement with the results derived from the recent study (Yu 2020). The precise different ground movement on both sides of F0 indicates that the fault is reactivated by panel 2553N coal mining. These results suggest that the subsidence in the study area is a dynamic process. Overall, the deviation obtained from the FDM model and InSAR method is small and within the standard deviation error limits. The InSAR- derived surface subsidence curve coincides with the FDM model simulation result.

5 Conclusion

This paper studies the ground movement in Fangezhuang coalfield in Tangshan city using a complete 3D FDM analysis in conjunction with stacking InSAR analysis. In order to obtain a reliable model, an orthogonal test was applied to back analyse the strata parameters, and the subsequent FDM model result was compared with the ground movement measured using InSAR. The differences in the maximum ground movement obtained from the FDM model analysis and stacking InSAR analysis were within 25%. Both results show that when the fault is activated by underground mining, the maximum subsidence value will not be located precisely above the centre of the mining panels. Instead, it is shifted towards the fault, and the geological fault clearly affects the spatial distribution of ground movement. The most dominant subsidence occurs in mining panel 2553N and the area between panel 2553N and fault F0, and has subsided by up to 57 cm. Overall, the ground movement patterns and magnitudes obtained by two different methods have a relatively good consistency, which attests the efficacy of the FDM numerical model in simulating the impact of the mining activity. The ground movement pattern in the Fangezhuang coal mine is spatially controlled by the geological fault to some degree, causing differential subsidence that could affect infrastructures. The method proposed in this paper can help to improve the understanding of subsidence control, develop preventive measures in Fangezhuang coalfield that consider fault reactivation, and forecast the impact of future underground mining activities in Fangezhuang coal mine. Further numerical analyses are required at a more realistic mine scale to study the ground movement trends in more detail.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding authors.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

The project was supported by the University of Nottingham Faculty of Engineering Research Excellence PhD scholarship. The work was supported financially by the Fundamental Research Funds for the Central Universities (2019XKQYMS63).

Conflict of Interest

Author DG was employed by the company Terra Motion 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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors would like to thank the European Space Agency (ESA) for providing Sentinel-1 data for the analysis. The authors would also like to thank Fangezhuang coal mine for providing the geological report with levelling data and strata parameters.

References

Abbasi, B., Russell, D., and Taghavi, R. (2013). “FLAC3D Mesh and Zone Quality,” in FLAC/DEM Symposium (China: Itasca Consulting Group).

Google Scholar

Agarwal, V., Kumar, A., L. Gomes, R., and Marsh, S. (2020). Monitoring of Ground Movement and Groundwater Changes in London Using InSAR and GRACE. Appl. Sci. 10, 8599. doi:10.3390/app10238599

CrossRef Full Text | Google Scholar

Albino, F., J Biggs, C. Y., and Li, Z. (2020). 'Automated Methods for Detecting Volcanic Deformation Using Sentinel‐1 InSAR Time Series Illustrated by the 2017–2018 Unrest at Agung, Indonesia. J. Geophys. Res. Solid Earth 125, e2019JB017908. doi:10.1029/2019jb017908

CrossRef Full Text | Google Scholar

Bell, F. G., Donnelly, L. J., Genske, D. D., and Ojeda, J. (2005). Unusual Cases of Mining Subsidence from Great Britain, Germany and Colombia. Environ. Geol. 47, 620–631. doi:10.1007/s00254-004-1187-9

CrossRef Full Text | Google Scholar

Bell, J. W., Falk, A., Ferretti, A., Bianchi, M., and Novali, F. (2008). Permanent Scatterer InSAR Reveals Seasonal and Long‐term Aquifer‐system Response to Groundwater Pumping and Artificial Recharge, 44. Water Resources Research. doi:10.1029/2007wr006152

CrossRef Full Text | Google Scholar

Berardino, P., Fornaro, G., Lanari, R., and Sansosti, E. (2002). A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms. IEEE Trans. Geosci. Remote Sensing 40, 2375–2383. doi:10.1109/tgrs.2002.803792

CrossRef Full Text | Google Scholar

Burbey, T. (2002). The Influence of Faults in basin-fill Deposits on Land Subsidence, Las Vegas Valley, Nevada, USA. Hydrogeology J. 10, 525–538. doi:10.1007/s10040-002-0215-7

CrossRef Full Text | Google Scholar

Cai, W., Dou, L., Si, G., and Hu, Y. (2021). 'Fault-induced Coal Burst Mechanism under Mining-Induced Static and Dynamic Stresses. Engineering 7, 687–700.

Google Scholar

Castellazzi, P., Longuevergne, L., Martel, R., Rivera, A., Brouard, C., and Chaussard, E. (2018). Quantitative Mapping of Groundwater Depletion at the Water Management Scale Using a Combined GRACE/InSAR Approach. Remote Sensing Environ. 205, 408–418. doi:10.1016/j.rse.2017.11.025

CrossRef Full Text | Google Scholar

Chen, C. W., and Zebker, H. A. (2001). Two-dimensional Phase Unwrapping with Use of Statistical Models for Cost Functions in Nonlinear Optimization. J. Opt. Soc. Am. A. 18, 338–351. doi:10.1364/josaa.18.000338

CrossRef Full Text | Google Scholar

Cheng, Q., Shi, Y., and Zuo, L. (2019). Numerical Simulation and Analysis of Surface and Surrounding Rock Failure in Deep High-Dip Coal Seam Mining. Geotech Geol. Eng. 37, 4285–4299. doi:10.1007/s10706-019-00907-z

CrossRef Full Text | Google Scholar

Deck, O., and Anirudh, H. (2010). Numerical Study of the Soil-Structure Interaction within Mining Subsidence Areas. Comput. Geotechnics 37, 802–816. doi:10.1016/j.compgeo.2010.07.001

CrossRef Full Text | Google Scholar

Diao, X., Bai, Z., Wu, K., Zhou, D., and Li, Z. (2018). Assessment of Mining-Induced Damage to Structures Using InSAR Time Series Analysis: a Case Study of Jiulong Mine, China. Environ. Earth Sci. 77, 1–14. doi:10.1007/s12665-018-7353-2

CrossRef Full Text | Google Scholar

Diao, X., Wu, K., Chen, R., and Yang, J. (2019). Identifying the Cause of Abnormal Building Damage in Mining Subsidence Areas Using InSAR Technology. IEEE Access 7, 172296–172304. doi:10.1109/access.2019.2956094

CrossRef Full Text | Google Scholar

Dong, L., Wang, C., Tang, Y., Tang, F., Zhang, H., Wang, J., et al. (2021). Time Series InSAR Three-Dimensional Displacement Inversion Model of Coal Mining Areas Based on Symmetrical Features of Mining Subsidence. Remote Sensing 13, 2143. doi:10.3390/rs13112143

CrossRef Full Text | Google Scholar

Du, S., Wang, Y., Zheng, M., Zhou, D., and Xia, Y. (2019). Goaf Locating Based on InSAR and Probability Integration Method. Remote Sensing 11, 812. doi:10.3390/rs11070812

CrossRef Full Text | Google Scholar

Fan, H., Wang, L., Wen, B., and Du, S. (2021). A New Model for Three-Dimensional Deformation Extraction with Single-Track InSAR Based on Mining Subsidence Characteristics. Int. J. Appl. Earth Observation Geoinformation 94, 102223. doi:10.1016/j.jag.2020.102223

CrossRef Full Text | Google Scholar

Fialko, Y., Simons, M., and Agnew, D. (2001). The Complete (3-D) Surface Displacement Field in the Epicentral Area of the 1999MW7.1 Hector Mine Earthquake, California, from Space Geodetic Observations. Geophys. Res. Lett. 28, 3063–3066. doi:10.1029/2001gl013174

CrossRef Full Text | Google Scholar

Fournier, T. J., Pritchard, M. E., and Riddick, S. N. (2010). 'Duration, Magnitude, and Frequency of Subaerial Volcano Deformation Events: New Results from Latin America Using InSAR and a Global Synthesis. Geochem. Geophys. Geosystems 11. doi:10.1029/2009gc002558

CrossRef Full Text | Google Scholar

Francioni, M., Stead, D., Sharma, J., Clague, J. J., and Brideau, M. A. (2021). An Integrated InSAR-Borehole Inclinometer-Numerical Modeling Approach to the Assessment of a Slow-Moving Landslide. Environ. Eng. Geosci. 27 (3), 287–305. doi:10.2113/EEG-D-20-00109

CrossRef Full Text | Google Scholar

Gee, D., Bateson, L., Sowter, A., Grebby, S., Novellino, A., Cigna, F., et al. (2017). Ground Motion in Areas of Abandoned Mining: Application of the Intermittent SBAS (ISBAS) to the Northumberland and Durham Coalfield, UK. Geosciences 7, 85. doi:10.3390/geosciences7030085

CrossRef Full Text | Google Scholar

Ghayournajarkar, N., and Fukushima, Y. (2020). 'Determination of the Dipping Direction of a Blind Reverse Fault from InSAR: Case Study on the 2017 Sefid Sang Earthquake, Northeastern Iran. Earth, Planets and Space 72, 1–15. doi:10.1186/s40623-020-01190-6

CrossRef Full Text | Google Scholar

Grebby, S., Orynbassarova, E., Sowter, A., Gee, D., and Athab, A. (2019). Delineating Ground Deformation over the Tengiz Oil Field, Kazakhstan, Using the Intermittent SBAS (ISBAS) DInSAR Algorithm. Int. J. Appl. earth observation geoinformation 81, 37–46. doi:10.1016/j.jag.2019.05.001

CrossRef Full Text | Google Scholar

Gumilar, I., AbidinAbidin, H. Z., HutasoitHutasoit, L. M., Hakim, D. M., Sidiq, T. P., and Andreas, H. (2015). Land Subsidence in Bandung Basin and its Possible Caused Factors. Proced. Earth Planet. Sci. 12, 47–62. doi:10.1016/j.proeps.2015.03.026

CrossRef Full Text | Google Scholar

Hoek, E., and Brown, E. T. (1997). Practical Estimates of Rock Mass Strength. Int. J. rock Mech. mining Sci. 34, 1165–1186. doi:10.1016/s1365-1609(97)80069-x

CrossRef Full Text | Google Scholar

Jeanne, P., Farr, T. G., Rutqvist, J., and Vasco, D. W. (2019). Role of Agricultural Activity on Land Subsidence in the San Joaquin Valley, California. J. Hydrol. 569, 462–469. doi:10.1016/j.jhydrol.2018.11.077

CrossRef Full Text | Google Scholar

Kang, Y., Zhao, C., Zhang, Q., Lu, Z., and Li, B. (2017). Application of InSAR Techniques to an Analysis of the Guanling Landslide. Remote Sensing 9, 1046. doi:10.3390/rs9101046

CrossRef Full Text | Google Scholar

Kumar, R., Choudhury, D., and Bhargava, K. (2016). Simulation of Rock Subjected to Underground Blast Using FLAC3D. JGS Spec. Publ. 2, 508–511. doi:10.3208/jgssp.ind-27

CrossRef Full Text | Google Scholar

Liao, H., Meyer, F. J., Scheuchl, B., Mouginot, J., Joughin, I., and Rignot, E. (2018). Ionospheric Correction of InSAR Data for Accurate Ice Velocity Measurement at Polar Regions. Remote Sensing Environ. 209, 166–180. doi:10.1016/j.rse.2018.02.048

CrossRef Full Text | Google Scholar

Mohammady, M., Pourghasemi, H. R., and Amiri, M. (2019). 'Land Subsidence Susceptibility Assessment Using Random forest Machine Learning Algorithm. Environ. Earth Sci. 78, 1–12. doi:10.1007/s12665-019-8518-3

CrossRef Full Text | Google Scholar

Motagh, M., Shamshiri, R., Haghshenas HaghighiHaghighi, M., Wetzel, H.-U., Akbari, B., Nahavandchi, H., et al. (2017). Quantifying Groundwater Exploitation Induced Subsidence in the Rafsanjan plain, southeastern Iran, Using InSAR Time-Series and In Situ Measurements. Eng. Geology. 218, 134–151. doi:10.1016/j.enggeo.2017.01.011

CrossRef Full Text | Google Scholar

Pacheco-Martínez, J., Hernandez-Marín, M., Burbey, T. J., González-Cervantes, N., Ortíz-LozanoZermeño-De-Leon, J. Á., Zermeño-De-Leon, M. E., et al. (2013). Land Subsidence and Ground Failure Associated to Groundwater Exploitation in the Aguascalientes Valley, México. Eng. Geology. 164, 172–186. doi:10.1016/j.enggeo.2013.06.015

CrossRef Full Text | Google Scholar

Park, I., Choi, J., Jin Lee, M., and Lee, S. (2012). Application of an Adaptive Neuro-Fuzzy Inference System to Ground Subsidence hazard Mapping. Comput. Geosciences 48, 228–238. doi:10.1016/j.cageo.2012.01.005

CrossRef Full Text | Google Scholar

Parmar, H., Bafghi, A. Y., and Najafi, M. (2019). Impact of Ground Surface Subsidence Due to Underground Mining on Surface Infrastructure: the Case of the Anomaly No. 12 Sechahun, Iran. Environ. Earth Sci. 78, 1–14. doi:10.1007/s12665-019-8424-8

CrossRef Full Text | Google Scholar

Perry, A. P., Kalenchuk, K. S., and McKinnon, S. D. (2020). “A Numerical Modeling Investigation of Coupled-Mechanism Subsidence over an Underground Carlin Trend Mine,” in 54th US Rock Mechanics/Geomechanics Symposium (OnePetro).

Google Scholar

Ren, X. (2017). 'Law of Underground Pressure and Deformation Control in Fangezhaung Mining Panel. North China University of Science and Technology.

Google Scholar

Rizzoli, P., Martone, M., Gonzalez, C., Wecklich, C., Borla Tridon, D., Bräutigam, B., et al. (2017). Generation and Performance Assessment of the Global TanDEM-X Digital Elevation Model. ISPRS J. Photogrammetry Remote Sensing 132, 119–139. doi:10.1016/j.isprsjprs.2017.08.008

CrossRef Full Text | Google Scholar

Shabanimashcool, M., and Li, C. (2012). Numerical Modelling of Longwall Mining and Stability Analysis of the gates in a Coal Mine. Int. J. rock Mech. mining Sci. 51, 24–34. doi:10.1016/j.ijrmms.2012.02.002

CrossRef Full Text | Google Scholar

Shi, M., Yang, H., Wang, B., Peng, J., Gao, Z., and Zhang, B. (2021). Improving Boundary Constraint of Probability Integral Method in SBAS-InSAR for Deformation Monitoring in Mining Areas. Remote Sensing 13, 1497. doi:10.3390/rs13081497

CrossRef Full Text | Google Scholar

Shi, Y., Zhao, Minmin., and Jian, Hao. (2021). Study on Numerical Models in Predicting Surface Deformation Caused by Underground Coal Mining. Geotechnical Geol. Eng. 39 (6), 1–17.

CrossRef Full Text | Google Scholar

Sikora, P., and Wesołowski, M. (2021). Numerical Assessment of the Influence of Former Mining Activities and Plasticity of Rock Mass on Deformations of Terrain Surface. Int. J. Mining Sci. Tech. 31, 209–214. doi:10.1016/j.ijmst.2020.11.001

CrossRef Full Text | Google Scholar

Sowter, A., Bin Che Amat, M., Cigna, F., Marsh, S., Athab, A., and Alshammari, L. (2016). Mexico City Land Subsidence in 2014-2015 with Sentinel-1 IW TOPS: Results Using the Intermittent SBAS (ISBAS) Technique. Int. J. Appl. earth observation geoinformation 52, 230–242. doi:10.1016/j.jag.2016.06.015

CrossRef Full Text | Google Scholar

Sun, Q., Zhang, L., Ding, X. L., Hu, J., Li, Z. W., and Zhu, J. J. (2015). Slope Deformation Prior to Zhouqu, China Landslide from InSAR Time Series Analysis. Remote Sensing Environ. 156, 45–57. doi:10.1016/j.rse.2014.09.029

CrossRef Full Text | Google Scholar

Torres, R., Snoeij, P., Geudtner, D., Bibby, D., Davidson, M., Attema, E., et al. (2012). GMES Sentinel-1 mission. Remote Sensing Environ. 120, 9–24. doi:10.1016/j.rse.2011.05.028

CrossRef Full Text | Google Scholar

Wang, L., Deng, K., and Zheng, M. (2020). Research on Ground Deformation Monitoring Method in Mining Areas Using the Probability Integral Model Fusion D-InSAR, Sub-band InSAR and Offset-Tracking. Int. J. Appl. Earth Observation Geoinformation 85, 101981. doi:10.1016/j.jag.2019.101981

CrossRef Full Text | Google Scholar

Wang, Y., Wang, L., Zhang, Y., and Yang, T. (2015). “Investigation of Snow Cover Change Using Multi-Temporal PALSAR InSAR Data at Dagu Glacier, China,” in 2015 IEEE International Geoscience and Remote Sensing Symposium (IGARSS) (IEEE), 747–750. doi:10.1109/igarss.2015.7325872

CrossRef Full Text | Google Scholar

Xu, N., Kulatilake, P. H. S. W., Tian, H., Wu, X., Nan, Y., and Wei, T. (2013). Surface Subsidence Prediction for the WUTONG Mine Using a 3-D Finite Difference Method. Comput. Geotechnics 48, 134–145. doi:10.1016/j.compgeo.2012.09.014

CrossRef Full Text | Google Scholar

Yan, Y., Yan, W., Liu, J., and Guo, J. (20212021). The Prediction Model of Super Large Subsidence in High Water Table Coal Mining Areas Covered with Thick Unconsolidated Layer. Geofluids.

Google Scholar

Yang, X., Wen, G., Dai, L., Sun, H., and Li, X. (2019). Ground Subsidence and Surface Cracks Evolution from Shallow-Buried Close-Distance Multi-Seam Mining: a Case Study in Bulianta Coal Mine. Rock Mech. Rock Eng. 52, 2835–2852. doi:10.1007/s00603-018-1726-4

CrossRef Full Text | Google Scholar

Yang, Z., Li, Z., Zhu, J., Yi, H., Hu, J., and Feng, G. (2017). Deriving Dynamic Subsidence of Coal Mining Areas Using InSAR and Logistic Model. Remote Sensing 9, 125. doi:10.3390/rs9020125

CrossRef Full Text | Google Scholar

Yu, Q. (2020). 'Study on Mechanism of Mininig Subsidence Affected by Fault'. China Coal Research Institute.

Google Scholar

Yunjun, Z., Fattahi, H., and Amelung, F. (2019). Small Baseline InSAR Time Series Analysis: Unwrapping Error Correction and Noise Reduction. Comput. Geosciences 133, 104331. doi:10.1016/j.cageo.2019.104331

CrossRef Full Text | Google Scholar

Zhang, Y., Meng, X. M., Dijkstra, T. A., Jordan, C. J., Chen, G., Zeng, R. Q., et al. (2020). Forecasting the Magnitude of Potential Landslides Based on InSAR Techniques. Remote Sensing Environ. 241, 111738. doi:10.1016/j.rse.2020.111738

CrossRef Full Text | Google Scholar

Zhang, Z., Wang, C., Tang, Y., Zhang, H., and Fu, Q. (2015). Analysis of Ground Subsidence at a Coal-Mining Area in Huainan Using Time-Series InSAR. Int. J. Remote Sensing 36, 5790–5810. doi:10.1080/01431161.2015.1109725

CrossRef Full Text | Google Scholar

Zheng, M., Deng, K., Fan, H., and Du, S. (2018). Monitoring and Analysis of Surface Deformation in Mining Area Based on InSAR and GRACE. Remote Sensing 10, 1392. doi:10.3390/rs10091392

CrossRef Full Text | Google Scholar

Keywords: finite difference method, ground movement, Interferometric SAR, stacking, geological fault, 3D modelling

Citation: Qin Z, Agarwal V, Gee D, Marsh S, Grebby S, Chen Y and Meng N (2021) Study of Ground Movement in a Mining Area with Geological Faults Using FDM Analysis and a Stacking InSAR Method. Front. Environ. Sci. 9:787053. doi: 10.3389/fenvs.2021.787053

Received: 30 September 2021; Accepted: 16 November 2021;
Published: 06 December 2021.

Edited by:

Junshi Xia, Riken, Japan

Reviewed by:

Bruno Adriano, Riken Center for Advanced Intelligence Project (AIP), Japan
Yosuke Aoki, The University of Tokyo, Japan

Copyright © 2021 Qin, Agarwal, Gee, Marsh, Grebby, Chen and Meng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhengyuan Qin, zhengyuan.qin@nottingham.ac.uk; Yong Chen, chenyong@cumt.edu.cn

Download