Oilfield Reservoir Parameter Inversion Based on Ground Deformation 2D Measurements Acquired by Msbas Time Series InSAR Method

1. Introduction

oil, as the largest energy source in the world, its extraction and storage are of great significance to national economic development. However, it also brings environmental problems that cannot be ignored [1,2,3]. Oil extraction can cause not only subsidence but also uplift. This is because fluid injection into the reservoir is often carried out to enhance oil recovery, resulting in excessive underground reservoir pressure, which manifests as a ground uplift [4,5]. Whether it is subsidence or uplift, severe ground deformation not only damages oil well production facilities but also affects the safety of surrounding infrastructure such as buildings, railways, and roads. Therefore, ground deformation monitoring and reservoir parameter inversion in oil extraction areas are essential for ensuring the safety of oilfields. The Qaidam Basin, often referred to as China’s “treasure trove, is not only abundant in various mineral resources but also hosts important oil and gas accumulation zones in its western region. By the year 2020, the Qaidam Basin boasts 32 identified oil and gas fields, with a combined proven reserve exceeding 6.4 billion tons [6]. Previous studies on this area mostly described the reservoir characteristics of oilfields from the perspectives of geology, sedimentary petrology, and geochemistry [7,8], there hasn’t been any research on utilizing Interferometric Synthetic Aperture Radar (InSAR) technology for deformation characteristics monitoring and reservoir parameter inversion. Therefore, using InSAR technology to monitor the deformation of the oilfields in the western basin not only provides guidance for sustainable oil extraction but also offers a new approach to understanding the reservoir characteristics in this region.

Conventional monitoring techniques require significant time and effort. However, InSAR technology, due to its benefits of extensive coverage, all-weather capabilities, precision, and efficiency, has found broad application in monitoring urban ground deformation [9,10,11], as well as in the monitoring and early warning of geological disasters such as earthquakes, landslides, and others [12,13,14]. Currently, the most commonly employed InSAR techniques include Permanent Scatterer InSAR (PS-InSAR), Small Baseline Subset InSAR (SBAS-InSAR), and Distributed Scatterer InSAR (DS-InSAR), which can not only effectively overcome the shortcomings of phase decoherence, terrain residual and atmospheric delay, but also quickly obtain the deformation time series in line-of-sight (LOS) direction [15]. Meanwhile, some scholars have also applied InSAR technology to monitor ground deformation in oilfields, geothermal fields, and mines, which is caused by underground reservoir changes, thus exposing the intimate connection between ground deformation characteristics and the state of underground reservoirs [16,17,18]. Sun et al. [19] used multi-track PS-InSAR technology to quantitatively analyze three types of ground subsidence, including oilfields, in the lower Liaohe plain of China. Juncu et al. [20] extracted deformations in a geothermal field by SBAS InSAR technology and found that the primary factor contributing to the subsidence was the observed pressure drawdown. Chen et al. [21] used Differential InSAR (D-InSAR) and SBAS technology to acquire cumulative deformation data within a mining area in China and compared it with the progress of mining operations. However, one-dimensional deformation measurement can not only fail to provide more comprehensive and intuitive information for detecting the ground deformation characters of oilfields but also struggle to reflect the real ground deformation conditions, often leading to deformation interpretation deviations, which is called the LOS fuzzy problem.

In addition, for oilfields, changes in reservoir pore pressure caused by fluid extraction and injection, as well as geometric shapes, positions, volumes, and other physical parameters of the reservoir, greatly influence the ground deformation characteristics. Therefore, using InSAR deformation results to invert reservoir parameters in oilfields can provide a quick and timely understanding of the reservoir state, which is of great significance for oilfield stability evaluation. Klemm et al. [22] monitored an oilfield in the Middle East with PS-InSAR and inverted ground deformations using a geomechanical model. They found that monitoring ground deformations and applying geomechanical inversion can yield valuable insights into the dynamic characteristics of reservoirs. Yang et al. [23] obtained subsidence information of a hydrocarbon reservoir in West Texas, USA, using the SBAS method and simulated ground subsidence based on reservoir parameters. By comparing the InSAR observed deformations with the modeled deformations, they found a high degree of agreement between the two. Ji et al. [24] conducted subsidence monitoring in the Karamay oilfield in Xinjiang by the SBAS method and then inverted its reservoir geometric parameters based on the Okada model. Nevertheless, there is still limited research on applying InSAR deformation results to reservoir parameter inversion in oilfields. The inversion theory is relatively lacking, and most scholars use a single model for the inversion under the assumption that there exists only one source of deformation in the subsurface. There is a lack of comparative studies on joint inversion of multiple-source models for complex deformations. Besides, some scholars have jointly inverted the seismic source parameters with the LOS and azimuthal deformation fields of ascending and descending to improve the inversion efficiency and accuracy by increasing the constraints of nonlinear inversion [25,26], but there are no similar attempts currently in the field of subsurface reservoir parameters inversion, and they all directly utilize the one-dimensional deformation in the inversion.

Based on the 113 ascending and descending Sentinel-1A SAR images from January 2021 to December 2022, Multidimensional Small Baseline Subset (MSBAS) InSAR technology is employed in this study to obtain the vertical and east-west ground deformation of the oilfields in the western Qaidam Basin and then analyze its temporal and spatial characteristics and causes. Afterwards, we combine the obtained two-dimensional deformation measurements for constraints and use a single-source model and a dual-source model to invert the reservoir parameters of the two complex deformation areas in the Youshashan oilfield within the study area, and then comparative analysis and evaluation of the inversion results are conducted. The related results can provide references for the safety maintenance of oilfields in the western and the selection of reservoir parameter inversion models and methods in other areas.

2. Study Area and Datasets

2.1. Background of the Study Area

Located in the northern part of the Qinghai-Tibetan Plateau, the Qaidam Basin is one of four major basins in China. It is surrounded by the Kunlun Mountains, Qilian Mountains, Altyn Mountains, and other mountains, covering an expanse of about 240,000 square kilometers. Within the basin, there are not only extensive salt lakes and marshes but also abundant reserves of oil, coal, and various metallic minerals.

The study area is located in the western part of China’s Qaidam Basin. Its eastern region is characterized by hills and mountains and its western region features relatively flat topography, with an average altitude of approximately 3000 meters. The local climate is a typical plateau continental climate, with strong sunshine, cold winters and cool summers, and less rainfall. In this region, the Paleocene and Neocene systems are mainly developed, and the Quaternary is uplifted and thinly deposited. The reservoirs can be broadly categorized into two types: porous clastic rocks and fractured, pore-vug diamictites. In the vertical direction, these two reservoir types are alternately stacked, forming a reservoir distribution characteristic that develops continuously in time and alternately overlaps in space. The unique geological and geographical conditions have fostered abundant oil resources. Following the commencement of development in the Youshashan (YSH) oilfield in 1957, the Huatugou (HTG) oilfield, the Youyuangou (YYG) oilfield, the Shizigou (SZG) oilfield, and the Ganchaigou (GCG) oilfield within the study area were subsequently developed. At the same time, the fluid injection technique was employed to increase oilfield production. By 2014, the total oil production in this region had exceeded 6.5 million tons [27].

Figure 1. (a) Overview of the western part of Qaidam Basin overlayed on the topography. (b) Detailed satellite image of the study area. The oilfields’ boundary is depicted by the white polygon.

Figure 1. (a) Overview of the western part of Qaidam Basin overlayed on the topography. (b) Detailed satellite image of the study area. The oilfields’ boundary is depicted by the white polygon.

2.2. Data and Processing

The dataset used in this study is 58 ascending and 56 descending Sentinel-1A SAR images covering the period from January 2021 to December 2022. Sentinel-1A is a satellite equipped with C band synthetic aperture radar with an orbital height of 690 km and a revisit period of 12 days. Strict orbital control technology is adopted to ensure the space baseline is small enough to improve the SAR image interference effect. GAMMA software [28] was used for D-InSAR processing in this study. At the same time, terrain errors are removed by the digital elevation model (DEM) with 30 m resolution, and to remove the atmospheric delay, the Generic Atmospheric Correction Online Service (GACOS) data is also introduced. Finally, to improve the quality of the InSAR result, the interference pairs with spatio-temporal baselines less than 200 meters and 49 days are selected to perform the phase unwrapping by the minimum-cost flow (MCF) method.

Figure 2. Spatio-temporal baseline of the ascending (a) and descending (b) interfering pairs.

Figure 2. Spatio-temporal baseline of the ascending (a) and descending (b) interfering pairs.

Table 1. Main parameters of the SAR data.

Sensor	Wavelength	Azimuth/Range pixel spacing	Orbit direction	Path	Temporal coverage
Sentinel-1A	5.6cm	13.99cm/2.33cm	Ascending	143	2021.01.11-2022.12.20
			Descending	48	2021.01.05-2022.12.26

Table 1. Main parameters of the SAR data.

Sensor	Wavelength	Azimuth/Range pixel spacing	Orbit direction	Path	Temporal coverage
Sentinel-1A	5.6cm	13.99cm/2.33cm	Ascending	143	2021.01.11-2022.12.20
			Descending	48	2021.01.05-2022.12.26

3. Methodology

3.1. MSBAS-InSAR Method

MSBAS technology is improved on the basis of SBAS technology. The result generated by SBAS is the time series of LOS displacement for each pixel in a SAR image. In the case of a single acquisition geometry without supplementary data, it is not possible to fully decompose the LOS solution obtained from SBAS into its three components of ground displacement (vertical, east-west, and north-south directions). However, when the InSAR data come from multiple acquisition geometries, the SBAS method can be adjusted to generate an approximate solution that includes time series of multiple components [27]. Therefore, by simultaneously processing multiple orbit D-InSAR data, MSBAS can obtain the vertical and horizontal east-west deformation time series of the overlapping area [28]. The main steps of the MSBAS method are as follows: firstly, D-InSAR processing is performed on one or more sets of ascending and descending orbit data, followed by resampling the unwrapped and geocoded differential interferograms to a consistent area and grid size. Secondly, a temporal matrix is created and decoherent pixels are removed. Finally, each pixel is solved by singular value decomposition (SVD) [29] to acquire the two- dimensional deformation rate component, and then the deformation time series can be reconstructed by numerically integrating the deformation rate. The use of Tikhonov regularization during the matrix solving process helps solve ill-posed problems in the observation equation. The inversion matrix is:

$$\left(\begin{array}{cc}-\cos \theta \sin \phi \Delta t& \cos \phi \Delta t\\ \lambda L& \end{array}\right)\left(\begin{array}{c}{V}_E\\ V{}_U\end{array}\right)=\left(\begin{array}{c}\widehat{\varphi }\\ 0\end{array}\right)$$

(1)

where $\phi , \theta$represent the incident and azimuth angle, $\Delta t$represents the time interval between adjacent SAR images, $\lambda$is a regularization parameter, L represents the identity matrix, $V_E$ and $V_U$ are the east-west and vertical deformation rate, respectively, $\widehat{\varphi }$ is the obtained differential interferometry measurements. If we denote the coefficient matrix by C:

$$C=\left(\begin{array}{cc}-\cos \theta \sin \phi \Delta t& \cos \phi \Delta t\\ \lambda L& \end{array}\right)$$

(2)

Then, we can obtain the east-west and vertical deformation rate by the following formulas:

$$\left(\begin{array}{c}{V}_E\\ V_U\end{array}\right)={(C^{T}C)}^{-1}{C}^T\left(\begin{array}{c}\widehat{\varphi }\\ 0\end{array}\right)$$

(3)

3.2. The Source Model and Inversion Method

In 1977, Matsu’Ura [30] inverted the fault properties of the Tango earthquake based on geodetic data and clearly proposed the concept of “geodetic inversion” for the first time. Geodetic inversion mainly consists of three components: geodetic data, geophysical inversion model, and inversion algorithm. In this study, the deformation data obtained through InSAR technology is used as geodetic data, and appropriate inversion models and algorithms will be employed to obtain the underground reservoir parameters of the deformation area in the oilfields. Next, the geophysical inversion models and inversion algorithms used in this study will be introduced.

3.2.1. Finite Prolate Spheroidal Model

In 1988, Yang et al. [31] derived and calculated the analytical formula for the arbitrary-oriented prolate spheroidal cavity model within a finite-dimensional elastic half-space. They proposed the finite prolate spheroidal model and found that this model has relatively more parameter degrees of freedom compared to the ellipsoidal point source model. The Finite Prolate Spheroidal Model is composed of eight parameters, including the three-dimensional coordinates of the ellipsoid’s center, the long axis, the short axis, pressure variation, as well as the strike and dip of the long axis. A space rectangular coordinate system with O as the origin is established, as shown in Figure 3. M represents the coordinates $\left(x_0, y_0,-d\right)$ of the ellipsoid source center, a and b are the major and minor semi-axis of the ellipsoid, respectively, $\phi$ is the strike, and $\theta$ is the dip angle.

3.2.2. Dipping Dike with Uniform Opening Model

In 1985, Japanese scholar Yoshimitsu Okada [32] analyzed existing research results on ground deformation caused by elastic half-space faults and proposed a general expression for ground deformation caused by fault dislocations of point sources and finite rectangular surface sources. In 1992, Okada [33] further improved the theory of this fault dislocation model. According to the elastic half-space isotropic dislocation theory, the displacement of a point on the ground due to the dislocation of a rectangular surface within an elastic medium is proportional to the amount of displacement of that surface. The proportionality coefficient is uniquely determined by the relative positions, the geometric dimensions, the inclination, the depth, and the elastic medium of the dislocation surface. If there are multiple dislocation surfaces underground, then the displacement of the point on the ground is the sum of the displacement vectors caused by the respective dislocations of the multiple dislocation surfaces.

Establish a space rectangular coordinate system with O as the origin, as shown in Figure 4. There are seven main parameters of the dipping dike with uniform opening model, including the length along the strike (L), the width along the inclination (W), the depth of the geometric center of the dislocation surface (d), the inclination angle ($\theta$), the strike ($\phi$), and the projected coordinates (x₀, y₀) of the geometric center point M on the ground surface.

3.2.3. The Nonlinear Bayesian Inversion Method

Geodetic inversion is a discipline that studies the evolutionary characteristics and laws of ground deformation based on geodetic observation data combined with prior information, and infers the internal and physical parameters of the earth, thereby revealing the internal dynamics of the earth [34]. In the process of geodetic inversion, the functional relationship between observed data and model parameters is as follows:

$$d=G(m)+\epsilon$$

(4)

where d denotes the observed data, m represents the model parameters, G is the Green’s function that connects the model parameters and ground deformation, and $\epsilon$ is the model error.

Compared with the traditional deterministic inversion algorithm that determines the optimal model parameters by minimizing the error between observation data and model simulation data, the nonlinear Bayesian inversion algorithm in the stochastic inversion algorithm used in this paper [35]. By obtaining the posterior probability density function (PDF) of all model parameters, this algorithm not only can determine the optimal parameter values but also can evaluate the uncertainty of the parameter values well, which is more practical than the traditional algorithms. In addition, this paper proposes to utilize the 2D deformation field measured by MSBAS to constrain and invert the vertical and east-west deformations of the oilfields, so as to obtain more reliable reservoir parameters of the oilfield. At the same time, quadtree sampling is applied to improve inversion efficiency. The likelihood function p(d|m) can be expressed as follows:

$$p(d|m)={(2\pi )}^{-N/2}{\left|{\sum }_d\right|}^{-\frac{1}{2}}\times \exp \left[-\frac{1}{2}{(d-Gm)}^T{\sum }_d^{-1}(d-Gm)\right]$$

(5)

where $d={\left[d_u, d_e\right]}^T$ represents the vertical and east-west observation data, m represents the model parameters, N refers to the number of data points, ${\Sigma }_d$ refers to the variance-covariance matrix of the data, and $Gm={\left[G_{u}m,G_{e}m\right]}^T$is the 2D simulated data calculated from the model parameters m through the Green’s function G.

After obtaining the prior information and likelihood function, the PDF $p(m|d)$ can be computed by the Bayesian formula. It refers to the probability of the current parameter effectively accounting for the observed data when considering the prior information.

In addition, the Markov Chain Monte Carlo (MCMC) sampling method and the Metropolis-Hastings rule [36,37] will be used to randomly modify the model parameters and calculate new likelihood function values. The newly calculated likelihood function value will be compared with the previously computed value to decide whether the current model parameters are approved or rejected. This process is iterated until the termination condition is met.

4. Results

4.1. 2D Deformation Monitoring and Analysis

The MSBAS method is employed in this research to simultaneously process ascending and descending SAR images, leading to the generation of a vertical average annual deformation rate map and an east-west average annual deformation rate map of the study area, as depicted in Figure 5a and 5b, respectively. It can be observed that from January 2021 to December 2022, the GCG oilfield, HTG oilfield, and YSS oilfield in the study area exhibited varying degrees of ground deformation. Specifically, the GCG oilfield and YSS oilfield experienced notable ground uplift and east-west deformation, with the highest uplift rate reaching 48 mm/year. Only uneven ground subsidence occurred without significant horizontal displacement in the HTG oilfield, with a maximum subsidence rate of about 12 mm/year. Previous studies have indicated that the injection of fluid and gas during oil extraction can cause ground uplift[16]. The study area is characterized by multiple faults and complex structures, posing challenges in oil production. To overcome these challenges and achieve increased and sustained production, water injection methods have been widely adopted in the oilfields. Therefore, the ground uplift in the GCG oilfield and YSS oilfield may be associated with the oilfield’s water injection activities. However, although the HTG oilfield employed the same method to enhance production, its ground didn’t rise but subsided. This may be due to the fact that the HTG oilfield is located in the urban area of Mangya City, where residents predominantly depend on groundwater extraction for their daily water supply. The substantial groundwater and oil extraction as well as city construction has resulted in the compaction of the ground in this area. Considering the complexity of the ground subsidence causes in the HTG oilfield, this study will focus on discussing and analyzing the deformation characteristics of the uplift areas (area A, area B, area C) within the GCG oilfield and the YSS oilfield.

As shown in Figure 5, area A exhibits a maximum uplift rate of 16 mm/year, along with a maximum east-west deformation rate of 8 mm/year. Most of the area B moves westward, and there is an eastward displacement in a small local area. The maximum uplift rate reaches 30 mm/year, and the maximum east-west deformation rate reaches 7 mm/year. The maximum uplift rate of area C is 48 mm/year, and it shows horizontal deformation characteristics with east-west opposing movement centered on the uplift center. The eastern flank of the uplift center is shifting towards the east with a velocity of approximately 15 mm/year, while the western flank of the uplift center is moving in the westward direction at a rate of about 6 mm/year. Moreover, the ground deformation in area C is near China National Highway 315, potentially posing risks to the road and its surrounding infrastructure. Therefore, it is necessary to enhance the supervision and upkeep of this area. The YSS oilfield is located at the high point of the Youshagou in the northwest of the Youshashan structure in the Qaidam Basin. It is a typical abnormally low-pressure reservoir so relying solely on natural energy for development is difficult to achieve satisfactory results. Therefore, water injection development has been conducted in the YSS oilfield since 1996, and it has currently transitioned into the late production stage. The long history of water injection and the large amount of water injection required for production maintenance coincide with the serious ground deformation in area B and area C. However, the GCG oilfield, which just started trial production in 2020, boasts substantial oil and gas resources and has a small amount of water injection, so the ground deformation there is also relatively small.

To investigate the spatiotemporal evolution patterns of each uplift in the research area, Figure 6a,b respectively show the time series cumulative vertical and east-west deformation from January 2021 to December 2022. According to Figure 6, area A has been uplifting throughout the observation period, and the uplifted area has gradually expanded southward. As of December 20, 2022, the cumulative uplift reached 34 mm. In the horizontal direction, it moved westward significantly from May 11, 2021, to June 28, 2021, and then resumed its slow movement. Area B exhibited a small relative uplift before September 2021, and a significant increase in relative uplift with time after September 2021, while the east-west deformation always remained relatively stable and slowly increasing. Both the vertical and horizontal deformation rates in area C increased slowly before February 2022 but significantly accelerated thereafter. The maximum vertical displacement was up to 106 mm, and the maximum horizontal displacement was up to 42 mm. Ground deformation is an intuitive reflection of changes in underground reservoir conditions. The greater the intensity of fluid injection into the oil layer, the faster the pore pressure of the reservoir increases. Therefore, we speculate that the variation in ground deformation rate is related to the variation in injection intensity.

Furthermore, we choose the six points’ (P1-P6) (Figure 5) time series InSAR results to specifically analyze the temporal changes in the oilfield ground, as illustrated in Figure 7. P1 and P2 are both located in area A. From their deformation time series, it is observed that although the surface in area A has been uplifted compared to January 2021, there is a brief subsidence followed by a rebound every year from May to August. This may be due to a decrease in reservoir pressure caused by oil extraction during this time period, leading to surface subsidence, but it quickly rebounds due to the injection of fluid. P3 and P4 are situated on the west and east sides of area B, respectively. Unlike the roughly linear deformation trend of P3, the P4 point remained relatively stable before August 2021 but rapidly experienced uplift after August 2021, accompanied by a minor horizontal displacement of less than 5 mm. This may be because the intensity of water injection activities increased in the west side of area B since about August 2021. P5 and P6 are situated on the north and south sides of area C, respectively. They both exhibit a gradual deformation trend before February 2022 and an accelerated deformation trend after February 2022, which is consistent with the increasing injection efforts required for the late-stage production of the oilfield.

We also plotted three profile lines, L1, L2, and L3, in area A, B, and C respectively. The two-dimensional cumulative deformation along these profile lines is shown in Figure 8. The vertical axis represents the cumulative displacement, and the horizontal axis represents the distance along the direction indicated by the arrows in Figure 5. L1 exhibits an uplift in the vertical direction and reaches its maximum deformation at approximately 1.3 km, after which it gradually decreases. In the horizontal direction, there is a cumulative displacement towards the west before 1.4 km, while there is a cumulative displacement towards the east after 1.4 km. From the analysis above, it can be inferred that there may be a deformation center in area A in the vicinity of 1.3 km to 1.4 km. Unlike the single peak feature of L1, L2 and L3 exhibit two uplift peaks. Combined with their horizontal deformation characteristics, we speculate that the deformation in areas B and C is relatively complex and there is more than one deformation center.

4.2. Reservoir Modeling

Inverting the reservoir parameters based on the deformation information obtained through InSAR technology, combined with corresponding geophysical models allows us to have a deeper understanding of the underground process of injecting fluids and the reservoir’s dynamic behavior. Currently, commonly used geophysical models include the Mogi model [38], the Okada model [33], and the ellipsoid model [39]. The point source model approximates the deformation as a circle, which differs significantly from the observed deformation, while the Okada model and finite ellipsoid model use more adjustable parameters to simulate complex deformation, which is closer to the observed deformation. Therefore, we choose the dipping dike with uniform opening model in the Okada model and the finite prolate spheroidal model to invert the reservoir parameters in area B and area C in this study. Additionally, considering the complex deformation features exhibited in these two areas, this study not only uses the most common single-source model for inversion but also introduces a dual-source model combining a rectangular dislocation surface and a dipping ellipsoid for inversion.

Due to the ambiguity of one-dimensional LOS deformation, it may be difficult for nonlinear inversion to converge quickly. Therefore, we combine InSAR vertical and east-west cumulative ground deformation variables to constrain the inversion, which can increase the robustness of nonlinear inversion. Subsequently, the Bayesian inversion algorithm is used to search for the optimal parameters, and predefined search ranges are set for each parameter. For example, the rectangular length is set from 100 to 1000 m, and the width is set from 10 to 500 m. After 1,000,000 iterations, the parameters reach convergence. In all inversion cases, we typically configure Poisson’s ratio to a value of 0.25.

4.2.1. Inversion Result of Area B

We assume that the ground deformation in area B is the result of a rectangular dislocation source in the shallow reservoir and use the dipping dike with uniform opening model to invert the oilfield reservoir parameters. The optimal parameters acquired through inversion are shown in Table 2. It can be observed that the length and width of the rectangular dislocation surface source are 980.668 m and 399.642 m, respectively. The depth range of the reservoir is 171-227 m, trending from northwest to southeast with an inclination angle of approximately 3°, indicating it is nearly horizontal. The positive value of the opening suggests reservoir expansion. Figure 10 gives the 3D schematic of a rectangular dislocation surface source with optimal parameters in area B.

Figure 9. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T1 corresponds to the profile line in Figure 13; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

Figure 9. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T1 corresponds to the profile line in Figure 13; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

We assume that the ground deformation in area B arises from the combined effect of a dipping ellipsoid source and a rectangular dislocation surface source in the shallow reservoir. Accordingly, we use the finite prolate spheroidal model and the dipping dike with uniform opening model to invert the reservoir parameters of the oilfield. The optimal parameters acquired through inversion are shown in Table 3. For the ellipsoid source, its burial depth is about 218 m. The length of its major and minor semi-axis indicates its morphology as a flat ellipsoid, and its minor axis (30 m) reflects the thickness of the reservoir. The ellipsoid’s strike angle relative to the north direction is 65° (clockwise direction), and the dip angle is 17°. DP/mu stands for the ratio of the pressure change to the shear modulus. For the rectangular surface source, its burial depth is about 232 m and its length and width are 728 m and 297 m, respectively. The strike angle is 39° (counterclockwise) and the dip angle is nearly horizontal. The opening of the dislocation surface is 0.063m. Figure 12 presents the 3D schematic of the two sources in area B.

Figure 11. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T1 corresponds to the profile line in Figure 13; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

Figure 11. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T1 corresponds to the profile line in Figure 13; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

Figure 12. 3D schematic of the two deformation sources in area B.

Figure 12. 3D schematic of the two deformation sources in area B.

4.2.2. Comparative Analysis of Inversion Results for Area B

The residuals between the observed and simulated deformation fields serve as a reliable indicator of the inversion results quality. From the residual plots in Figure 9 and Figure 11, it can be observed that the residuals in the east-west direction for the single-source model are large and widely distributed, whereas the residuals for the dual-source model are smaller in both two directions and mostly exist in the far-field region, which may be due to the influence of the random noise and topography. Additionally, the models used in this study are based on the assumption that they are modeled in a homogeneous elastic half-space, whereas the area B geology is somewhat inhomogeneous. Considering the inherent limitations of the models themselves (e.g., lack of complete rotational freedom), the presence of some residuals, even in the near-field region, is expected.

Furthermore, to compare and evaluate the inversion effects of the two models more intuitively and quantitatively, we drew the residual histograms and the deformation profiles along T1 between their two-dimensional simulated deformation fields and observed deformation fields, as shown in Figure 13. From the distribution of residuals, it can be seen that the mean and root mean square values of vertical and east-west residuals of the dual-source model are smaller than those of the single-source model. However, compared with the vertical direction, the east-west residuals of both models are smaller and more concentrated towards the zero. This may be due to the fact that the ground uplift caused by the increase of pore pressure is mainly vertical uplift, and the horizontal deformation is smaller and simpler, so the simulation effect of the east-west direction is better than that of the vertical direction. On the other hand, in terms of the fit between the observed and simulated deformations along T1, the simulated curves of the dual-source model are more consistent with the trend of the observed curves, so its degree of fit is higher. In summary, whether from the residual analysis or the fitting degree of deformation, the dual-source model is more suitable for the inversion of reservoir parameters for complex deformations and yields more reliable inversion results than the single-source model.

Figure 13. (a) The distribution of vertical and horizontal residuals for the single-source and dual-source model (b) The distribution of vertical residuals for the dual-source model. (c) The distribution of horizontal residuals for the single-source model. (d) The distribution of horizontal residuals for the dual-source model. (e) The vertical deformation profile alone T1. (f) The horizontal deformation profile along T1. .

Figure 13. (a) The distribution of vertical and horizontal residuals for the single-source and dual-source model (b) The distribution of vertical residuals for the dual-source model. (c) The distribution of horizontal residuals for the single-source model. (d) The distribution of horizontal residuals for the dual-source model. (e) The vertical deformation profile alone T1. (f) The horizontal deformation profile along T1. .

4.2.3. Inversion Results for Area C

The comparative analysis of the inversion results of the single-source model and the dual-source model in area B reveals that the dual-source model is more suitable for inverting the reservoir parameters of complex deformation. Considering the extremely complex deformation in area C, which cannot be reliably inverted using the single-source model, this paper directly applies the dual-source model to invert the complex deformation in area C. The optimal parameters acquired through inversion are shown in Table 4. It can be seen from the inversion results that the ground deformation of area C is caused by a rectangular dislocation surface source on the north side of the shallow reservoir and an inclined ellipsoid source on the south side. For the rectangular fault surface source, its length and width are approximately 667 m and 74 m, respectively, with a burial depth of around 451 m. As for the dipping ellipsoid source, its major semi-axis is approximately 728 m, the minor semi-axis is around 17 m, and the burial depth is 546 m. A 3D schematic of two sources with optimal parameters for area C is given in Figure 15. In Figure 16, the root-mean-square values of the vertical and east-west residuals are 6.45 mm and 5.95 mm, respectively, and the fit between the simulated deformations and the observed deformations along the profile T2 is also high, which indicates that this inversion is effective and the reliability of the obtained parameters is high, and further validates the feasibility of the dual-source model.

Figure 14. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T2 corresponds to the profile line in Figure 16 ; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

Figure 14. The observed deformation field in vertical (a) and horizontal (d) directions, the black arrow T2 corresponds to the profile line in Figure 16 ; The modeled vertical deformation field in vertical (b) and horizontal (c) directions; The residuals in vertical (c) and horizontal (e) directions.

Figure 15. 3D schematic of the two deformation sources in area C.

Figure 15. 3D schematic of the two deformation sources in area C.

Figure 16. (a) The distribution of vertical residuals for the dual-source model. (b) The distribution of horizontal residuals for the dual-source model. (c) The vertical deformation profile alone T2. (e) The horizontal deformation profile along T2.

Figure 16. (a) The distribution of vertical residuals for the dual-source model. (b) The distribution of horizontal residuals for the dual-source model. (c) The vertical deformation profile alone T2. (e) The horizontal deformation profile along T2.

5. Discussion

In this paper, we use the MSBAS method to monitor the ground deformation of the oilfields in the western Qaidam Basin based on Sentinel-1A satellite data. Derived from the SBAS technology, the MSBAS technology inherits the technical advantages of SBAS, which can overcome the influence of interferometric phase decorrelation and correct multiple errors, including terrain residual, orbital errors, and atmospheric delay errors. In addition, regularization or time filtering can remove high-frequency noise and improve the temporal resolution.

Different from only obtaining the one-dimensional deformation field in the line-of-sight direction, we acquired the two-dimensional time series deformation field in vertical and east-west directions by MSBAS technology. A series of results indicate that besides vertical deformation, there is also a significant amount of deformation in the horizontal east-west direction at the ground of the oilfields, with a maximum east-west deformation of 16 mm/year. Unlike the LOS deformation, which is a mixture of vertical and horizontal deformation, independent vertical and east-west deformations can more comprehensively and intuitively reflect the real deformation effects of the oilfields. This provides more reliable support and guidance for oil field development planning and land subsidence management.

Subsequently, based on the acquired two-dimensional cumulative deformation, we perform reservoir modeling for two complex deformation areas using geophysical models. In real work scenarios, the injection and production activities in oilfields are often carried out simultaneously at multiple places, resulting in complex ground deformation. By comparing and analyzing the simulation effects of the single-source model and the dual-source model, this study found that there are often multiple deformation sources underground in areas with complex deformation. The deformation field simulated by the single-source model is relatively simple and cannot effectively reflect the observed deformation field. The dual-source model not only simulates smaller residual mean values and root mean square values than the single-source model but also has a higher fit with the observed deformation field. Therefore, for areas with complex deformation, the reservoir parameters inverted using the dual-source model are more reliable and better represent the mapping relationship between reservoir changes and ground deformation in oilfields. Moreover, the two-dimensional deformation variable used as the model input variable in this study adds nonlinear inversion constraints compared to most current one-dimensional inversion methods, so the inversion efficiency and accuracy are further improved. Unfortunately, due to the lack of water injection and oil production data for the oilfields in the study area, further validation of the inversion results cannot be conducted. In future research on reservoir parameter inversion, efforts should be made to strengthen the collection of such data and the validation of models.

6. Conclusions

This study conducted vertical and east-west deformation monitoring of the oilfields in the western Qaidam Basin by simultaneously processing ascending and descending Sentinel-1A SAR images from 2021 to 2022. Based on this, the reservoir parameters of two deformation areas were inverted using single-source and dual-source models. The main results are as follows:

In the study area, the HTG oilfield exhibited uneven ground subsidence, with a maximum subsidence rate of 12 mm/year. At the same time, the GCG oilfield and YSS oilfield experienced substantial ground uplift, with a maximum rate of 48 mm/year. Along with this uplift process, the east-west deformation rate in this area also reached 16 mm/year. The cause of the ground uplift in the oilfield may be related to increased reservoir pore pressure resulting from water injection. Changes in water injection intensity may lead to changes in deformation rates.

Combining the analysis of inversion results, it can be concluded that the introduction of a two-dimensional deformation field helps improve the non-uniqueness of the inversion results, enhance the robustness of the inversion process, and consequently obtain more reliable oil reservoir parameters. Furthermore, the dual-source model is more suitable for inverting complex deformation reservoir parameters compared to the single-source model.