Present-Day Three-Dimensional Deformation across the Ordos Block, China, Derived from InSAR, GPS, and Leveling Observations

Abstract

The Ordos Block in China experiences tectonic activity and frequent earthquakes due to compression from the Tibetan Plateau and extension from the North China Block. This has prompted the construction of a high-resolution three-dimensional (3D) deformation field to better understand the region’s crustal movement. Considering the limitations of the existing geodetic observations, we used InSAR, GPS, and leveling observations to create a high-precision 3D deformation field for the Ordos Block. Spherical wavelet decomposition was used to separate tectonic and non-tectonic deformation signals. Short-wavelength non-tectonic deformation fields revealed complex surface deformation patterns caused by groundwater, oil, gas extraction, and coal mining. Long-wavelength tectonic deformation fields showed subsidence in the southern margin of the block, while the interior and northeastern margins were uplifted. By combining imaging results from the seismic velocity structure and magnetotellurics, we infer that the upwelling of deep materials beneath the northeastern margin leads to surface uplift with tensile strain rates. The crustal uplift in the area south of 38°N matches the thickening of the lower crust. The weak subsidence and eastward horizontal movement disappearing near 108°E at the southern margin support the existence of asthenosphere flow beneath the Qinling orogenic belt.

1. Introduction

The Ordos Block is situated at the intersection of three sub-plates in North China, South China, and the Tibetan Plateau (Figure 1). The internal faults in the block are not developed, and the frequency and intensity of present-day seismic activity are low, rendering it a well-known stable ancient craton in China. The active fault system around the Ordos Block is one of the main active faults in China and is a well-known seismic activity zone. Historically, five strong earthquakes of magnitude eight or above have occurred, accounting for approximately one-fourth of the recorded earthquakes of magnitude eight or above in China [1] (Figure 1). The area surrounding the Ordos Block includes a series of fault basins, including the Shanxi Basin on the eastern margin, the Weihe Basin on the southern margin, the Yinchuan Basin on the western margin, and the Hetao Basin on the northern margin. Due to the unique geographical location and seismic activity of the Ordos Block, its formation, evolution, and tectonic characteristics have been of great interest to scientists. The long-term three-dimensional (3D) differential motion of the Earth’s crust leads to the accumulation and abrupt release of strain energy, which is one of the primary causes of earthquakes [2,3]. Utilizing modern geodetic techniques to monitor and obtain regional 3D crustal motion velocity field images and their dynamic changes is the basis for deepening our understanding of continental deformation dynamics. Additionally, they play a crucial role in predicting strong earthquakes in block boundary zones [4].

The advancements in space geodetic techniques, such as the Global Positioning System (GPS), Interferometric Synthetic Aperture Radar (InSAR), and precise leveling measurements, have been widely employed to monitor the crustal deformation of the Ordos Block. Previous studies have indicated that within the 420,000 km² area of the Ordos Block, only approximately 600 GPS observation points were available for estimating the horizontal crustal deformation [5,6,7]. InSAR technology, known for providing high-spatial-resolution surface deformation fields, has yielded favorable results in monitoring deformation along the southern and eastern margins of the Ordos Block [8,9,10,11]. However, existing InSAR studies have yet to cover the northern, western, and internal regions of the block, resulting in a dearth of comprehensive InSAR deformation data that encompasses the entirety of the Ordos Block. Previous researchers have utilized multi-epoch precise leveling data spanning from 1970 to 2014 to obtain the vertical velocity field of the peripheral regions of the Ordos Block. However, leveling lines in the internal region of the block, as well as the northern Yinshan orogenic belt and the southern Qinling orogenic belt, are sparsely distributed (Figure 1) [12].

Figure 1. Tectonic background and distribution of InSAR, GPS, and leveling data. The black line represents active faults, and thick solid lines represent major fault zones. Green arrows represent GPS’s horizontal crustal movement velocity field (1991–2019) [7]. Blue/red arrows represent the vertical crustal movement velocity field from leveling measurements (1970–2015) [12]. The purple box represents Sentinel-1 satellite SAR image coverage. Blue/red circles represent M ≥ 6.0 earthquakes (400 BC-present). The inset represents the Ordos Block location and surrounding blocks. The blue polygon represents active fault zones at level I [13]. OR represents Ordos Block/Tibetan Plateau; SC represents South China Block; YM represents Yunnan-Myanmar Block; NC represents North China Block; WC represents West China; and NEA represents Northeast Asia Block.

Figure 1. Tectonic background and distribution of InSAR, GPS, and leveling data. The black line represents active faults, and thick solid lines represent major fault zones. Green arrows represent GPS’s horizontal crustal movement velocity field (1991–2019) [7]. Blue/red arrows represent the vertical crustal movement velocity field from leveling measurements (1970–2015) [12]. The purple box represents Sentinel-1 satellite SAR image coverage. Blue/red circles represent M ≥ 6.0 earthquakes (400 BC-present). The inset represents the Ordos Block location and surrounding blocks. The blue polygon represents active fault zones at level I [13]. OR represents Ordos Block/Tibetan Plateau; SC represents South China Block; YM represents Yunnan-Myanmar Block; NC represents North China Block; WC represents West China; and NEA represents Northeast Asia Block.

InSAR, GPS, and precise leveling measurement techniques have been widely used for crustal deformation monitoring; however, the observational data provided by these techniques have significantly different spatiotemporal characteristics. InSAR can provide a high-resolution surface deformation field. However, because of the geometric characteristics of InSAR side-looking imaging, it can only provide one-dimensional surface deformation results along the radar line-of-sight (LOS). GPS is currently the most widely used method for monitoring 3D surface deformation, with horizontal accuracy reaching the submillimeter level but with relatively low vertical accuracy. Precise leveling measurements can only provide discrete vertical deformation information along leveling routes, owing to the constraints of terrain conditions and observation costs. Owing to natural conditions and observation costs, existing GPS and precise leveling monitoring stations in the Ordos Block and its surrounding areas are sparsely distributed and insufficient for inverting surface deformation and strain accumulation at the kilometer scale (Figure 1). Therefore, InSAR, GPS, and precise leveling measurements have complementary advantages. Consequently, this study integrated the three methods to obtain a unified reference frame and a high-density 3D deformation field to study present-day crustal deformation in the Ordos Block and its surrounding areas. This 3D deformation data plays a crucial role in understanding the mechanisms of crustal deformation within continental regions and assessing the seismic hazard risk in the Ordos Block.

2. Data and Methods

2.1. InSAR Data Processing

The Sentinel-1A/B satellites, which are two C-band synthetic aperture radar (SAR) satellites launched by the European Space Agency in 2014 and 2016, operate within a predefined orbit tube with a radius of 50 m to ensure the coherence of SAR images through strict orbit control technology. As the Sentinel-1 satellites only provide ascending mode data over the Ordos Block and its surrounding areas, we collected five ascending mode Sentinel-1 single-look complex (SLC) images from October 2015 to April 2020, totaling 9786 scenes, to obtain InSAR deformation fields in the study area (Figure 1). Given the wide distribution of active faults around the Ordos Block, we adopted a long-strip map InSAR processing strategy that stitches consecutive frames from the same orbit into SLC data to avoid losing long-wavelength deformation signals. This approach enabled us to obtain 466 long-strip-map SLCs.

Before the SAR data were registered, the Precise Orbit Data (POD) product was used to refine the precise orbit parameters of the Sentinel-1 satellite. The Copernicus DEM data, with a resolution of 30 m, were utilized as the reference DEM for both SAR data registration and the simulation of the differential terrain phase. The elevation reference surface of the Copernicus DEM was transformed from EGM2008 to a WGS84 ellipsoidal height system to eliminate the impact of geodetic reference frame differences. To increase the coherence of the interferogram and avoid errors caused by an excessive number of views, a two-pass multi-looking approach was applied to the interferogram. During the first pass, an (8 × 2) small multi-looking ratio was utilized to produce a high-resolution differential interferogram that fully utilizes the topographic data provided by the DEM. The high-resolution differential interferogram was multi-looked again during the second pass (4 × 4) to yield a spatial resolution of approximately 110 m. To further enhance the signal-to-noise ratio of the interferogram and coherence estimation accuracy, a homogenous point filter was used to filter the interferogram [14]. The minimum cost flow (MCF) algorithm based on Delaunay triangulation performed differential interferogram phase unwrapping. Finally, the InSAR Generic Atmospheric Correction Online Service (GACOS) estimated the atmospheric phase delay [15].

This study adopted a combination of interferograms with long temporal and short spatial baselines, as shown in

Table 1. Threshold list of temporal–spatial baselines.

Number	Temporal Baselines ΔT	Spatial Baselines ΔB
1	275 d < ΔT < 455 d	ΔB < 100 m
2	640 d < ΔT < 820 d	ΔB < 200 m

and Figure 2. We set two thresholds for the temporal and spatial baselines to ensure that the acquisition times of the reference and secondary images were in the same season. Consequently, 7266 interferograms were obtained, of which, 1260 high-quality interferograms were selected for InSAR data processing. The InSAR time series analysis produced high-resolution InSAR deformation maps of the Ordos Block and its surrounding areas (Figure 2).

The accuracy of the InSAR deformation field was assessed through internal and external consistency evaluations. An internal consistency assessment typically involves computing the standard deviation of the deformation values in stable areas within the study region. However, the large size of the study area renders the identification of a stable, representative area challenging. On the contrary, an external consistency assessment involves validating the InSAR results using ground-based monitoring data, such as GPS and leveling measurements. The study area had abundant GPS monitoring and leveling measurement data, as depicted in Figure 1, which enabled an external consistency assessment. The observation periods for the external ground-based monitoring data were separate from those of the InSAR observation period (2015–2021). For instance, we utilized the latest GPS horizontal velocity field data published by Hao et al. in 2021, covering the observation period from 1999 to 2019 [7]. For leveling measurements, we chose data published by Hao et al. in 2016, covering the period from 1970 to 2014 [12]. This dataset comprises the most extensive and densest GPS and leveling observations available for the Ordos region. Figure 1 shows the point information for the two types of external ground-based monitoring data. Although these data represent long-term average velocity fields, they could still be used to reference InSAR monitoring results.

The InSAR deformation field is shown in Figure 2g. To evaluate the accuracy of the InSAR monitoring, the LOS deformation of InSAR was compared with the projection of the GPS horizontal deformation and GPS horizontal + leveling vertical deformation in the LOS direction (Figure 2f). To this end, the average InSAR value within a 1 km radius around the GPS station and the average leveling measurement value within a 2 km range were used. The external consistency assessment results indicated no significant systematic bias between InSAR and GPS horizontal deformation; the RMSE between the two was 1.37 mm/yr, with a linear regression correlation coefficient of 0.94. Similarly, there was no apparent systematic bias between InSAR and GPS horizontal + leveling deformation; the RMSE between the two was 1.78 mm/yr, with a linear regression correlation coefficient of 0.87.

The above comparative analysis demonstrates the reliable measurement accuracy of the InSAR deformation field of the Ordos Block. The higher consistency between the InSAR and GPS horizontal observations, as indicated by the root mean square error and correlation coefficient, suggests good agreement between the two methods. However, the addition of the leveling measurement results weakened these indicators, indicating differences between the vertical deformation component in the InSAR deformation field and the vertical deformation in the leveling data. This inconsistency could be attributed to several factors. The vertical deformation rates in the Ordos Block and its periphery have changed over the past few decades, whereas the long-term horizontal deformation rates exhibit a linear trend. Secondly, the vertical deformation of the crust in the Ordos Block and its periphery is complex, with significant influences from groundwater, oil, gas, and coal mining activities. Hence, differences in the vertical deformation within a 2 km radius around GPS stations may exist.

2.2. Three-Dimensional Deformation Field Construction and Multiscale Decomposition

InSAR LOS observations are known to be insensitive to deformation in the north–south direction. However, it is possible to obtain high-resolution vertical and east–west deformation fields by combining InSAR observations from both the ascending and descending orbits. Nevertheless, this study utilized only InSAR observations from the ascending orbit, which limited the ability to obtain a comprehensive 3D crustal deformation field of the Ordos Block and its surrounding areas. However, GPS and leveling measurements provided valuable discrete deformation information in the horizontal and vertical directions. The technical characteristics of the InSAR, GPS, and leveling measurements are complementary and could be combined to generate a unified, high-density 3D deformation field. This approach is particularly important for studying the crustal deformation of the Ordos Block and its surrounding areas.

To obtain a high-resolution 3D deformation field of the Ordos Block and its surrounding areas, this study collected and organized GPS and leveling measurement data [7,12] and fused them with the ascending orbit InSAR deformation field. To address the rank deficiency or ill-conditioning of the fusion function model due to the limited and uneven spatial distribution of the GPS and leveling observation data, spatial interpolation needed to be applied to these two datasets. Kriging interpolation, a statistical-based spatial interpolation method that utilizes covariance functions to interpolate a random field [16], was employed in this study. This method is known for its high accuracy and capability to measure estimation errors, making it widely applicable in geoscience fields. First, the study area was divided into fine triangular grids with a resolution of 1 km. Subsequently, based on the assumption of a constant strain rate inside each triangle using elastic theory, the 3D velocity field of any observation within the triangle was represented using the 3D velocity fields of the three vertices of the triangle. Afterward, a velocity function model (Equation (1)) was established between the grid nodes and the ground survey observation points (InSAR, GPS, and leveling measurements) inside the triangles [17,18,19,20,21,22].

$$\left[\begin{array}{cc}{G}_{sar}\hfill & G_{orb}\hfill \\ G_{gps}\hfill & \hfill 0\hfill \\ G_{level}\hfill & \hfill 0\hfill \\ {\kappa }_s{}^2{\nabla }_s{}^2\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{c}{m}_{vel}\\ m_{orb}\end{array}\right]=\left[\begin{array}{l}{d}_{sar}\\ d_{gps}\\ d_{level}\\ 0\end{array}\right]$$

(1)

Equation (1) comprises the design matrices $G_{sar}$, $G_{gps}$, and $G_{level}$ for the InSAR, GPS, and leveling measurements, respectively. In addition, ${\nabla }_s^2$ represents the finite-difference approximation of the Laplacian smoothing operator and ${\kappa }_s^2$ determines the smoothing weights. Additionally, $G_{orb}$ is the orbital design matrix to solve for any remaining orbital errors in the InSAR data, $m_{vel}$ is a 3D crustal deformation, and $m_{orb}$ represents a first-order polynomial fitting of the trend surface, which characterizes the reference difference between the InSAR and GPS, leveling frames. The InSAR data utilized in this study underwent orbital error correction. However, there is a certain disparity between the reference frame of the InSAR LOS deformation field and the Eurasian reference frame of the GPS and leveling data. Therefore, in our functional model, $m_{orb}$ represents the difference between the reference frames; $d_{sar}$, $d_{gps}$, and $d_{level}$ represent the observation values of the InSAR, GPS, and leveling measurements, respectively.

The Helmert variance component estimation (VCE) method is an algorithm utilized for the posterior estimation of variances. Initially, it assigns initial weights to different types of observations and conducts a preliminary adjustment using the least squares method. Subsequently, based on the obtained adjustment corrections, the method iteratively estimates the variances of the observations according to specific principles until the unit weight standard deviations of all observation types are equal. In this study, the Helmert VCE method was used to determine the weight ratio between the InSAR, GPS horizontal velocity field, and leveling data as 1:4.3:0.5. Finally, the 3D velocity field of the nodes was inverted using the least squares method, resulting in a high-density and high-precision 3D deformation field for the Ordos Block and its surrounding areas (Figure 3).

The spherical wavelet multi-scale analysis technique extends wavelet theory from an infinite planar space to a finite spherical space and has found wide applications in various fields, including GNSS velocity fields, gravity fields, magnetic fields, and temperature fields [23]. Spherical wavelet multi-scale analysis enables the separation of large-scale low-frequency deformation in the 3D velocity field and the characterization of small-scale high-frequency local deformation features. In this study, a Gaussian spherical wavelet multi-scale analysis method was adopted to obtain the 3D crustal deformation field of the Ordos Block and its surrounding area at various spatial scales [23,24].

The global Euler vector method initially separated the 3D deformation trend from the original 3D deformation field. Subsequently, multi-scale spherical wavelet analysis was conducted on the residual deformation field. To obtain the 3D deformation field at different wavelengths ranging from 6.9 km to 1763.4 km in the study area, we chose the maximum scale factor of the spherical model (qmin = 2) and the minimum scale factor (qmax = 10) (

Table 2. Geometric characteristics of spherical triangle meshes [23].

Scale Factor q	Faces	Vertices	Side Length (km)
0	20	12	7053.6 km
1	80	42	3526.8 km
2	320	162	1763.4 km
3	1280	642	881.7 km
4	5120	2562	440.8 km
5	20,480	10,242	220.4 km
6	81,920	40,962	110.2 km
7	327,680	163,842	55.1 km
8	1,310,720	655,362	27.5 km
9	5,242,880	2,621,442	13.8 km
10	20,971,520	10,485,762	6.9 km
11	83,886,080	41,943,042	3.4 km
12	335,544,320	167,772,162	1.7 km

). We discovered that at a wavelength of 55.1 km (q = 7), the small- and medium-scale deformation resulting from underground resource exploitation, such as groundwater, oil and gas, and coal mining, could be effectively separated, thereby extracting the large-scale tectonic deformation field caused by tectonic movement (Figure 3).

3. Three-Dimensional Deformation Characteristics

3.1. Horizontal Deformation Characteristics

Considering the relationship between the 3D deformation field and reference frame, it is crucial to note that the motion states derived from different reference frames can be significantly different. Conversely, the strain-rate field is independent of the frame. Therefore, to accurately image the strain rate field of the Ordos Block and its surrounding areas, we utilized the Median Estimation of Local Deformation (MELD) method, which is capable of robustly imaging the strain rate field from a high-resolution 3D deformation field while effectively suppressing anomalous 3D deformation [25] (Figure 4). Our results show that the strain rate within the Ordos Block was relatively low, whereas high strain rates were mainly concentrated near the peripheral faults of the block. The southwestern margin of the block, which is the foreland area where the Tibetan Plateau extends northeastward, was subjected to northward and eastward squeezing forces from the Tibetan Plateau, indicating that this region was under intense squeezing stress (Figure 4b).

Additionally, the entire eastern margin of the Ordos Block is under extensional stress due to the southeastward extensional stress from the North China Block. The shear strain rate field (Figure 4a) shows that the southwestern margin of the block exhibits intense squeezing and strike-slip characteristics, with the highest shear strain rate occurring in the middle section of the Haiyuan Fault Zone, reaching up to 50 nanostrains/yr. The nearby Liupan Mountain Fault, adjacent to the Haiyuan Fault, exhibited comparatively small shear strain rates, indicating that the left-lateral strike-slip motion of the eastern section of the Haiyuan Fault was transformed into compressional motion on the Liupan Mountain Fault. The strain rate of another fault adjacent to the Haiyuan fault, the Xiangshan-Tianjingshan fault, is also relatively high. To quantitatively estimate the slip rates of these two faults, we selected a profile crossing both faults and transformed the deformation direction parallel to the fault. We then used the screw dislocation model to calculate the slip rate of the Xiangshan-Tianjingshan fault section, which was approximately 2.58 mm/yr with a locking depth of 21.04 km. The slip rate of the Haiyuan Fault was approximately 2.93 mm/yr with a locking depth of 15.62 km.

A multi-scale analysis of the 3D velocity field (Figure 3) revealed that the short-wavelength horizontal deformation field (with wavelengths less than 55 km) primarily reflects near-fault deformation and some non-tectonic deformation. In contrast, the long-wavelength horizontal deformation field (with wavelengths greater than 55 km) exhibited solid compressive and dextral strike-slip features along the southwestern margin of the Ordos Block, which is mainly bounded by the Haiyuan-Liupanshan fault zone. This fault zone transformed the northeastward movement from the Tibetan Plateau into northward and eastward movements along the western and southern margins of the block, respectively. Consequently, dextral strike-slip and extensional belts developed on the eastern and western sides, respectively, and sinistral strike-slip and extensional belts developed on the northern and southern sides of the Ordos Block.

3.2. Vertical Deformation Characteristics

Vertical deformation in the Ordos Block and its surrounding areas is highly complex, exhibiting various surface deformation features at different spatial scales (Figure 3). The long-wavelength vertical deformation map (Figure 3g) depicts an uplift at the northeastern margin and inside of the Ordos Block and crustal sinking at the south margin. Large-scale structural deformation is mainly attributed to the tectonic movements of significant active faults along the perimeter of the block, forming the overall framework of the Ordos Block.

In the short-wavelength deformation map (Figure 3d), it is evident that several subsidence funnels developed in the northern and eastern faulted basins of the block, exhibiting subsidence rates exceeding 40 mm/yr, including the Hetao, Taiyuan, and Linfen-Yuncheng basins. Previous studies have also demonstrated the existence of rapidly subsiding areas in the Datong and Weihe basins [11]. However, the extent of subsidence has begun to narrow in recent years and uplifting has commenced in some regions. Moreover, certain areas within the block experience surface deformations owing to coal, petroleum, and natural gas mining, which exhibit small-scale and discrete distribution characteristics.

4. Discussion

4.1. Ground Deformation Caused by Groundwater Extraction

Ground subsidence is driven primarily by groundwater exploitation and active faults (Figure 3a). Most subsidence funnels are found on the hanging walls of normal faults and are controlled by two normal faults with opposite dip angles. A majority of subsidence funnels did not intersect the active faults; the ground subsidence on both sides of the faults differed significantly. The spatial distribution of ground subsidence zones runs parallel to the local fault trend, indicating the controlling effect of faults on the spatial extent of ground subsidence [9,10,11].

Ground subsidence is predominantly concentrated in areas where groundwater is extensively exploited in the peripheral faulted basins surrounding the Ordos Block. Cumulative ground subsidence can cause irreversible damage to infrastructure, such as buildings, pipelines, and roads, and poses serious threats to people’s production and safety. Many cities actively address the impact of ground subsidence by reducing groundwater extraction and implementing artificial groundwater recharge in areas with severe ground subsidence to increase groundwater replenishment. The rise in groundwater levels led to the shrinkage of ground subsidence areas and uplift in some regions. Areas of significant uplift in the peripheral fault basins of the Ordos Block include the Yungang, Taiyuan, Jinzhong, Fenyang, Xi’an, Pucheng, Fuping, Longxian, Fufeng, and Shizuishan areas in the Datong, Taiyuan, Weihe, and Yinchuan basins, respectively.

The Yuhuazhai area of Xi’an City, one of the largest urban villages in Xi’an and the area with the most severe ground subsidence in recent years, was selected as the region for underground water extraction and management in this study. Groundwater monitoring data were collected at Yuhuazhai (monitoring point 610113210027; https://geocloud.cgs.gov.cn, accessed on 27 December 2022) to extract the InSAR deformation time series (Figure 5i). Following the demolition of the Yuhuazhai urban village in 2018 and the closure of many self-provided wells, the ground deformation changed from subsidence to stability from May to October in 2018. On 29 October 2018, the groundwater in Yuhuazhai began to recharge and the groundwater level quickly recovered and stabilized at a depth of −15 m by 2020. The InSAR deformation time series revealed that the ground began to rebound between 14 October and 7 November 2018; the deformation tended to stabilize in 2020. The InSAR ground deformation and groundwater level changed synchronously, with no observed time lag effect.

4.2. Ground Deformation Caused by Oil, Gas, and Coal Mining

The Ordos Block, also known as the Ordos Basin, is the second-largest sedimentary basin in China and contains vast amounts of energy and mineral resources, including oil, natural gas, coal, coalbed methane, and sandstone-type uranium mines [26]. The Jiyuan, Wuqi, Zhidan, Huachi, and Xifeng oilfields (Figure 5g) within the Ordos Block are primarily owned by Shaanxi Yanchang Petroleum and China Petroleum Changqing Oilfields. Oilfield development began in Wuqi and Zhidan during the 1970s and 1980s, Jiyuan in the 1990s, and Huachi-Xifeng after 2000. Underground crude oil can be extracted using natural pressure during the initial oil-production phase. As oil reserves are depleted, extraction becomes more challenging and requires techniques, such as water injection, to maintain pressure and stability in the oilfield and facilitate crude oil extraction. However, during the later stages of water flooding development, the effectiveness of water injection can decrease due to declining oil well production and increased water injection, resulting in a more incredible injection amount than the output. The water flooding development can increase the formation pressure and cause surface uplift [27]. In Figure 5g, the red envelope delineates the extent of oilfield water flooding development. Significant surface uplift is observed in the water injection areas of the previously developed oilfields (Jiyuan, Wuqi, and Zhidan), with the maximum uplift rate observed in Zhidian reaching 15 mm/yr. In contrast, minor surface uplift is seen in the water flooding development areas of the recently developed oilfields (Huachi and Xifeng). This also suggests that the extensive uplift in the western part of the Ordos Basin is closely linked to oil–water flooding development stages.

The Shenfu-Dongsheng Coalfield (Figure 5h), located northwest of the Ordos Basin, is not only China’s largest coalfield with proven reserves but also one of the largest in the world. The coal seams in the coalfield are relatively shallow, making it an easily accessible mining area. The InSAR monitoring results revealed widespread ground deformation in the Shendong Coalfield, with multiple prominent deformation centers. The deformation areas were primarily rectangular or elliptical and exhibited high spatial consistency with the mining area. The InSAR monitoring results were highly effective in reflecting the spatial distribution characteristics of the influence range of mining subsidence and subsidence areas, thereby providing strong technical support for determining the boundaries of coal mining subsidence areas and for the early identification and prevention of geological hazards in subsidence areas. However, because of the large deformation gradient caused by coal mining, particularly the surface deformation caused by open-pit mining, the surface change often exceeds the limit of InSAR deformation monitoring. Therefore, complete restoration of the maximum deformation often requires small multi-views or high-frequency revisit data. Nevertheless, this study revealed the distribution characteristics of the deformation field, which effectively detected the mining area of coal mines and partially located open-pit mining areas. This demonstrates that InSAR technology can play a specific role in identifying the deformation range of coal mining subsidence areas and providing data support for post-mining land reclamation.

4.3. Ground Deformation Caused by Crustal Movement

After performing spherical wavelet decomposition, we removed the short-wavelength deformation caused by non-tectonic factors, such as groundwater, oil, natural gas, and coal mining. The long-wavelength vertical deformation field (Figure 3g) reflects the vertical tectonic movement characteristics of the Ordos Block, which are observed as crustal uplift at the northeastern margin of the block at a rate of approximately 5 mm/yr, crustal uplift within the block with a maximum uplift point between 37°N and 38°N at a rate of approximately 6 mm/yr, and subsidence in the Qinling Mountains orogenic belt at the southern margin of the block at a rate of approximately −3 mm/yr.

Seismic velocity structure and magnetotelluric imaging are effective techniques for detecting intricate lithospheric structures. The northeastern margin of the Ordos Block is a well-known distribution area for Datong Quaternary volcanoes. Previous studies on the crustal velocity structure and magnetotelluric imaging have suggested that mantle upwelling at the northeastern margin of the Ordos Block is still ongoing [28,29,30,31,32]. Our monitoring results indicate that the ground in the area is experiencing uplift and the strain rate shows tension, which further supports the upwelling of deep material beneath the northeastern margin of the Ordos Block (Figure 6).

The crustal velocity structure exhibits variations between the northern and southern regions of the Ordos Block, which are separated by 38°N. Specifically. The lower crustal thickness in the southern part of the Ordos Block was observed to be thicker than that in the northern part [34,35,36,37]. The long-wavelength vertical deformation field shows crustal uplift south of 38°N within the Ordos Block, which is consistent with the deep crustal structure and conforms to a new model of lower crustal thickening.

Previous studies have suggested that the Qinling Orogen accommodated the flow of low-viscosity lower crust from the Tibetan Plateau northward around the Sichuan Basin during the Cenozoic Himalayan-Tibetan orogeny [38,39]. However, our 3D deformation field results demonstrate that the Qinling Orogenic Belt was in a state of compression and exhibited weak subsidence in the long-wavelength vertical deformation field (Figure 3g and Figure 4b). Additionally, the east–west deformation field shows a clear eastward movement on the west side of the Qinling Orogenic Belt; however, this movement disappears east of 108°E (Figure 3h). Thus, we contend that the crustal deformation data do not support the occurrence of processes associated with the lower crustal flow beneath the Qinling Orogen. Other studies, based on velocity structure analysis, suggest the existence of an asthenospheric mantle flow channel with an eastward flow of material beneath the Qinling Orogenic Belt, located between the rigid blocks of the Ordos and Sichuan basins [28,40,41,42,43,44]. These studies propose that the direction of the asthenospheric mantle flow may deviate from that of the overlying lithosphere. Specifically, the asthenosphere does not passively drag the overlying lithosphere, nor can it induce the overlying plate to mimic its subsurface flow [44] (Figure 6). Our 3D deformation results are more consistent with asthenosphere flow, where the asthenosphere flows from west to east, causing subsidence, and decoupling between the asthenosphere and lithosphere leads to incomplete horizontal deformation, showing cessation of the eastward motion west of 108°E.

5. Conclusions

In this study, we employed a high-precision time-series InSAR technique to process Sentinel-1 SAR data from 2014 to 2020, which yielded a high-resolution InSAR deformation field for the Ordos Block and its adjacent regions with a deformation accuracy of 1.8 mm/yr. The obtained InSAR deformation field was subsequently integrated with GPS and leveling data, resulting in a high-resolution three-dimensional crustal deformation field. The multi-scale surface deformation field of the region was obtained by applying spherical wavelet decomposition and the deformation field was divided into tectonic and non-tectonic deformations using wavelength separation methods. This study found that non-tectonic deformation was mainly related to underground material exploitation, such as widespread subsidence and uplift of the surrounding block due to groundwater extraction and injection. These surface deformations are controlled by surrounding active faults and have a significant spatial correlation with them. Coal mining caused numerous discrete deformations within the block, whereas oil-water flooding resulted in local uplift. Furthermore, this study suggests that long-wavelength deformation around a block is associated with regional tectonic deformation. For instance, in the northeastern margin of the Ordos Block, the upwelling of deep materials causes surface uplift that manifests as extensional strain rates. In the area south of 38°N within the Ordos Block, the crustal uplift is consistent with a deep crustal structure, supporting a new model of lower crustal thickening. The disappearance of weak subsidence and eastward horizontal movement near 108°E at the southern margin of the Ordos Block supports the hypothesis of asthenospheric flow in the mantle beneath the Qinling Orogenic Belt.

While our study incorporates the latest multi-source geodetic observations, there are still certain limitations that need to be addressed. One such limitation arises from the disparate observation periods of InSAR, GPS, and leveling measurements, introducing uncertainties during the integration of the 3D deformation field. To enhance the robustness of the 3D deformation analysis, it is preferable to utilize both ascending and descending InSAR observations. Unfortunately, in our study area, only ascending Sentinel-1 images are available for analysis. Moreover, our high-resolution 3D deformation field can be further employed to investigate additional geoscientific issues in the Ordos region. For instance, the 3D deformation can be used to build models of the active faults surrounding the block. These models provide insights into the fault’s geometry, slip rates, and potential for future earthquakes. Additionally, the deformation field associated with oil, gas, coal mines, and groundwater extraction can be inverted to determine the precise locations and extent of underground resource exploitation. Such information would not only facilitate scientific productivity but also mitigate potential damage to surface infrastructure.