Abstract
It is difficult to predict tight sandstone gas in continental rift basin because of complex deep structure, fast transformation of sedimentary facies, large lateral variation of reservoir thickness, and tight and strong heterogeneity. Based on the petrophysical analysis of logging data, this paper improved the traditional low-frequency establishment method of prestack synchronous inversion which is under the constraints of fine isochronal sequence stratigraphy framework. The prestack AVO has attributed as covariate, and Bayes kriging method has applied to establish low-frequency model. Then, the relationship between density and P-wave and S-wave impedance has been introduced. The results show that the stability of density parameters has been improved, and the facies controlled prestack inversion with drilling verification has been realized which the number of layers with thickness error less than 20% accounts for 75.2%. Exploration practice shows that the facies controlled prestack inversion technology has important reference significance for other deep continental fault basins.
1. Introduction
Tight sandstone gas has become a key area of unconventional natural gas exploration which played an important role in the future oil and gas industry. There are a large number of studies on tight gas in the foreland basin depression valley slope, depression basin center, and craton slope negative structure with weak tectonic activity and relatively stable strata. Under the background of low porosity and low permeability, the reservoir of relatively higher porosity and higher permeability is called “sweet spot” of tight gas. Because of the characteristics of deep faults such as complex structure, fast sedimentary facies transformation, and strong heterogeneity, which makes lithology and gas bearing property difficult to predict. Great potential of remaining resources of tight glutenite gas reservoir in Sujiaweizi fault depression, northern Songliao Basin, but the Shahezi Formation in Xujiaweizi fault depression is in the period of strong fault depression, and the reservoir distribution is complex. Accurate identification of the spatial distribution of tight gas desserts is the key problem restricting the exploration and development of deep tight gas [1, 2].
In recent years, many scholars have made many important studies in the prediction of tight sandstone gas “sweet spot.” Goodway et al. [3] considered that the method of fluid prediction by using is more superior than the P-wave velocity ratio or model analysis method, which can obtained and to describe lithology and pore fluid by the method of seismic reflection. Goodway also pointed out that this method has the potential to distinguish gas bearing reservoir and gas water mixed layer; Gray et al. [4–6] improved the method Goodway which used lame parameters to reflect the characteristics of fluid and use them as fluid indicator for fluid identification. Russell et al. [7–9] improved the P-wave velocity equation on the basis of Biot Gassmann equation based on the previous research results, then deduced the fluid term of Gassmann equation, and applied it to reservoir fluid identification. Han et al. [8] studied the elastic impedance formula including Gassmann fluid term and proposed the prestack inversion method of pore fluid parameters. Li et al. [10, 11] combined with attribute analysis and ray elastic impedance inversion technology, and the “sweet spot” reservoir in Mahu sag of Junggar basin was predicted by using the idea of progressive control of facies, channel, physical property, and fracture, which reduced the multisolution. Summarizing the above research shows that reservoir prediction can be realized based on the correlation with a variety of elastic parameters and seismic attributes. Although the conventional prestack synchronous inversion method can transform seismic reflection information into rock elastic parameters and then quantitatively described the reservoir properties such as lithology and porosity, most seismic data only contains intermediate frequency information which low frequency and high frequency information are needed for inversion. Also, the accuracy of the low-frequency model directly affects the accuracy of the absolute elastic parameter prediction and the relative elastic property prediction [12, 13]. Therefore, it is very important to establish a reasonable low-frequency model for reservoir prediction, and the prediction of strong heterogeneous reservoir in fault lake basin requires higher accuracy of low frequency model.
In this paper, basis on the previous studies, the traditional prestack simultaneous inversion low-frequency modeling method was improved based on the constraints of fine isochronal sequence stratigraphic framework and the petrophysical characteristics of the study area. In this paper, the prestack AVO attribute has been taken as covariate [14, 15], the Bayes Kriging method to establish low-frequency model has been applied, and the relationship between P-wave and S-wave impedance and density has been introduced to improve the stability of density parameters; then, the reservoir parameters of tight glutenite in Shahezi Formation of Xujiaweizi fault depression are predicted. Based on the recognition that the characteristics of seismic data and the representation of amplitude and phase of seismic wave, in the process of low-frequency and high-frequency information fusion of well seismic combined inversion, this paper uses the data of seismic in-phase axis scale to drive the fine-layered interpretation results and prestack AVO attribute to constrain, to realize the well seismic information fusion, so as to improve the prediction accuracy of conventional prestack inversion methods [16–18].
2. Geological Setting
Xujiaweizi fault depression was formed by the coaction of the interaction in the Pacific tectonic domain and deep thermal process, and it is located in the east side of the paleo-central uplift in the northern Songliao basin. It was a Dustpan fault depression of a west fault and east overlap structure, generally NNW trending. Drilling revealed the deep fault depression from bottom to top Huoshiling Formation, Shahezi Formation, Yingcheng Formation, and Denglouku Formation is developed. The top surface of target layer Shahezi Formation is unconformity, and Shahezi Formation overlipped the Huoshiling Formation. Shahezi Formation develops fan delta and braided river delta depositional system. Sand bodies are distributed in fault depression in skirt belt along the edge of fault. Shallow lacustrine dark mudstone and swamp coal seams are high-quality source rocks (Figure 1). The thickness of the stratum is generally 500-1000 m, the central part of the fault depression can reach more than 1500 m, and the lithology includes glutenite, mudstone, sandstone, and coal [19]. Under the conditions of large-scale superimposition of source and reservoir and tight reservoir, the effective source rock controls the gas reservoir development area, the favorable facies belt controls the favorable reservoir development zone, and the reservoir development degree controls the gas reservoir enrichment.
3. Research Methods
3.1. Bayesian Global Optimization of Seismic Reflection Horizon Fine Interpretation
Fine structural interpretation is the basis of reservoir parameter identification. Through the combination of seismic geological horizon calibration and well seismic correlation, this paper tracked the third-order sequence seismic reflection horizon in the whole area, in which T41 is the top surface of Shahezi Formation, T41a is the bottom surface of Sq4 in the fourth member of Shahezi Formation, and T41b is the bottom surface of Sq3 in the third member of Shahezi Formation, T41c is the bottom surface of Sq2 in the second member of Shahezi Formation, and T42 is the bottom surface of the first member of Shahezi Formation. The comparative interpretation of seismic reflection layers gradually constrains densification and closure to realize three-dimensional interpretation from point to line and from line to plane [20, 21].
The thickness of single sand body in the study area is thin, while the thickness of fourth-order sequence is generally from 120 to 240 M. In order to clarify the distribution characteristics of single sand body or sand group and desserts, it is necessary to carry out fine interpretation of sand group level horizon. In this paper, an isochronous sequence stratigraphic model is established based on the existing fourth-order sequence stratigraphic framework and 3D seismic data volume. First, the waveforms of adjacent seismic traces are compared in the whole three-dimensional data space to identify the contact relationship of the same seismic event phase which up overtopping, top overtopping, down overtopping, and denudation, then connect the same seismic event phase with high waveform similarity. Second, the cost function of each three-dimensional connection mode of three-dimensional data volume was calculated based on Bayesian algorithm. Through global scanning, the minimum cost function that corresponds to the three-dimensional isochronal sequence stratigraphic model is obtained. The cost function is as follows [22]: where is the number of sample points in a grid node, and and are used to describe the location of sample points. The cost function takes into account the similarity and relative distance between nodes. By minimizing the cost function, the corresponding models which connect these nodes are selected, and the global sequence stratigraphic model is established. The optimization of sequence stratigraphic model can be iteratively optimized by continuously adding geological knowledge, and faults can be reasonably added to adjust a layer picked up as a standard layer. After adding these new constraints, a new three-dimensional stratigraphic model can be calculated iteratively. This paper starting from the seismic data, based on the global optimization, and fully considering the inheritance of stratum sedimentation and the overall geological rationality, the relative isochronal model is obtained, from which the interpretation horizon of each seismic codirection axis can be picked up (Figure 2), which effectively improved the interpretation accuracy of subdivision horizon [23, 24].
The subdivision layer position which is interpreted by this method is controlled by well stratification at the well point, and the direction of the same phase axis and the variation of the formation thickness are referred between wells. The variation of event direction and formation thickness is compared with each other between wells to realize seismic resolution axis interpretation of single sand layer or sand group, which lays the data foundation for the subsequent study of reservoir parameter prediction (Figure 3).
3.2. Quantitative Prediction of Reservoir Parameters by Facies Controlled Prestack Synchronous Inversion
From the petrophysical analysis results, the glutenite has the characteristics of high P-wave impedance, and the gas bearing glutenite has the characteristics of high P-wave impedance and low p-s-wave velocity ratio; from the seismic reflection characteristics, the gas bearing glutenite mainly shows class I and weak class II AVO.
In prestack synchronous inversion, the accuracy of low-frequency model directly affects the prediction accuracy of relative impedance and absolute impedance. As shown in formulas (2) and (3), under the constraint of isochronal stratigraphic model, the low-frequency model is established by using Bayesian kriging and AVO attribute constraint with double precision, so as to improve the accuracy of deterministic inversion and realize phase controlled prestack synchronous inversion. In this paper, the AVO attribute regression low-frequency model based on Bayesian kriging with double precision was established, which was mainly divided into the following steps: (1) loading the p-attribute and pseudo-Poisson’s ratio attribute of seismic prestack AVO attribute to reflect the spatial variation of wave impedance and density; (2) analyzing the correlation between AVO attribute and logging curve; (3) transforming the attribute of prestack AVO attribute; (4) taking the elastic parameters of low-frequency model as learning objectives and prestack AVO attribute and seismic transformation attribute body as constraints, then using the co-Kriging method, the model is established Based on the relationship between elastic parameters and seismic AVO attributes and seismic transformation attributes, the cross well elastic parameters are predicted; and (5) low-pass filtering is applied to the predicted elastic parameters to obtain the low-frequency trend model for inversion [25–27]. where is the estimation point, is the actual data, is the kriging coefficient, is the Lagrange multiplier, and is the variogram value between the known and unknown data.
Prestack synchronous inversion used AVA constrained sparse pulse inversion algorithm based on single channel to solve the following constrained optimization problems and minimize the objective function [24, 28]:
In the formula, the first four terms are the constrained sparse pulse inversion terms, and the last four terms are the additional constrained terms of prestack inversion. three terms are obtained by solving the minimum objective function. AVO/AVA constrained sparse pulse inversion algorithm can obtain absolute elastic parameter model. Due to the missing low-frequency components in the measured seismic data cannot be directly obtained by sparse pulse inversion, it is necessary to establish a low-frequency trend model as the constraint condition of inversion.
In order to enhance the stability of inversion algorithm and improve the accuracy of inversion results, the inversion of reflection coefficient and elastic parameters are divided into two parts. The specific implementation process are as follows: (1)First, the constrained sparse pulse reflection coefficient inversion algorithm is applied to obtain the sparse reflection coefficient sequence at each angle, and the objective function of minimization is as follows [29]: where and represent line number and track number, respectively; and are L-mode factors; is reflection coefficient; and are synthetic seismic record and original seismic data, respectively; and is balance factors.(2)Second, the AVA reflection coefficients obtained from the above inversion are weighted and superimposed to obtain the variation of elastic parameters. For each angle gather, there is an equation: , , …… , they represent the reflection coefficients of different angles, , , , they represent P-wave impedance, S-wave impedance, and density reflection coefficient, respectively. They represent P-wave impedance, S-wave impedance, and density reflection coefficient, respectively. Each term is a weight factor, which can be calculated by the coefficients of AKI Richards equation.(3)These elastic parameters are transformed into the initial estimates of P- and S-wave impedances and densities, and these initial estimates are unstable. Finally, we used the low-frequency trend information of P- and S-wave impedances and densities as constraints to minimize the objective function. Through continuous iteration, stable P- and S-wave impedances and densities can be obtained [30].
Constraints:
In the above, is the elastic parameter difference term, is the residual item of earthquake, is a low frequency trend, is a space constraint, is the time drift term, and are the minimum and maximum P-wave impedances for constraints, and are the minimum and maximum shear wave impedances for constraints, and and are the minimum and maximum densities used for constraints, respectively.
The method inputs multiple offset seismic bodies superimposed at different angles and combines logging elastic parameter data and seismic data calibration in different angle regions which can obtained seismic wavelet of different angles. Under the trend constraints, the elastic parameter data bodies such as shear wave impedance are obtained by used the above reflection coefficient inversion and elastic parameter inversion, which are suitable for reservoir prediction of complex lithology and impedance overlap [31].
The facies-controlled inversion method uses AVO attributes and isochronal sequence stratigraphic model constraints which can realize facies control, so that the vertical resolution of well logging data and the horizontal resolution of seismic data are reasonably integrated, while improving the vertical resolution need to maintain the horizontal resolution of seismic data. The results show that the prediction coincidence rate of glutenite above 10 m reaches 72% (Figure 4), and the inversion results have high coincidence rate with wells and can reflect the sedimentary characteristics of the fault depression strata, as shown in Figure 5.
3.3. Inversion Data Analysis
It can be seen from Figure 5 that the west of the study area, fan delta plain, and front subfacies are developed, while in the east, affected by the uplift and erosion of the late strata, there are few remains in the plain facies belt, and the front subfacies are mainly developed, with local interdistribution Bay microfacies.
Combined with seismic profile and sand ground ratio which is obtained from inversion calculation results of facies controlled prestack inversion, it is comprehensively judged as the west provenance fan, as shown in Figures 6 and 7; it is the development area of thin interbedding such as mud and coal seam, which is the east outer front, and the plain facies is denuded in a large area; and it is the east provenance fan, which is developed with thick glutenite. From the seismic section, we can clearly see the difference of seismic response caused by east-west provenance (Figure 6).
Because of the complexity and diversity of sedimentary structural belts, which has many types of sedimentary facies, and the lateral change of sedimentary facies belt is fast, so the sedimentary facies in different tectonic positions have different distribution characteristics and evolution rules. Based on the analysis of single well facies, combined with the data of logging, logging, core, and seismic response characteristics, in this paper, the correlation section and map of sedimentary facies is compiled, which reveals the vertical and horizontal distribution of sedimentary facies in Shahezi Formation (Figures 7 and 8).
The east-west differentiation of sedimentary facies belt in cross well seismic profile section is obvious. The eastern gentle slope area is characterized by widespread development of braided river deposits. On the horizontal plane, the sedimentary facies belt has a smooth transition, and in vertical direction, the braided river delta plain showed a deposition process of first retrogradation and then progradation. The western steep slope area was characterized by extensive development of fan delta deposits, the fan delta front deposits continue to prograde towards the center of the basin, and the sedimentary range expands. In the central area, the sedimentary water is relatively deep, and the shallow lake sediments are developed locally. The eastern and western sedimentary systems converge in the central area. Based on the provenance direction, seismic reflection characteristics, seismic attributes, and sand land ratio, the sedimentary facies map of the study area during the Sq4 period was compiled (Figure 7).
4. Discussion
This research is deployed well Songsh9p1 which is located in the NNW fault block in the south of Anda. Research found that the glutenite of fan delta developed in sequence Sq4 contacts with the mudstone of lacustrine facies in finger shape. The source and reservoir are overlapped. The overlying and surrounding dense glutenite or mudstone are the caprocks, which have favorable conditions for forming self-generating and self-storing gas reservoirs. Combined with the reservoir prediction results (Figures 7 and 9), we considered that the thickness of Sq43 glutenite in this well is large, and the reservoir is favorable. In order to evaluate the development and gas bearing capacity of glutenite reservoir in Shahezi Formation of northern fault depression, the well SongSh9p1 was deployed.
Well SongSh9p1 has completed drilling and gas testing, and the results show that the actual drilling effect of the well is good. The accumulated thickness of glutenite encountered in Sq43 is 126.6 m, the accumulated thickness of gas bearing glutenite encountered in Sq43 is 125.8 m, and the sand to ground ratio reaches 72.2%, which is consistent with the prediction. Sand gravel with fine gas bearing property was drilled in the main target layers, and after fracturing, 200300 cubic meters of industrial gas flow per day can obtained.
The results show that the accuracy of the low-frequency model directly affects the accuracy of the absolute elastic parameter prediction and also affects the accuracy of the relative elastic property prediction. The drilling results in the study area show that the facies controlled prestack inversion technology proposed in this paper can ensure the inversion accuracy of the elastic parameters inversion.
Data Availability
All data included in this study are available upon request by contact with the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
We gratefully acknowledge Northeast Petroleum University and the forward-looking basic technology research project of China National Petroleum Corporation (2021DJ0205) financial support.