1. Introduction
The Tarim Basin in western China is rich in oil and gas resources [
1], of which carbonate resources are about 7 × 10
9 t oil equivalent. The Ordovician carbonates in the Tarim Basin were partially eroded by the process of structural uplift, weathering and dissolution [
2]. The Middle-Lower Ordovician carbonates have strong heterogeneity, and the karst reservoirs in the strata carry most of the oil and gas resources [
3]. The target layer contains a large number of prolific oil and gas wells related to karst reservoirs [
4]. After years of exploration and development, it has been recognized that abnormal seismic reflection is the main feature of Ordovician carbonate fractured-vuggy reservoir (FVR) [
5]. The formation of FVR in the FuMan Oilfield (FMOF) is controlled by strike-slip faults [
6]. Hydrocarbon migrates along the fault zone and accumulates differentially into the FVRs, which has the characteristics of differential accumulation [
7]. The burial depth of strike-slip fault controlled reservoirs (SSFCRs) are 7000 to 10,000 m, and the height of the oil column has been proved to reach 550 m, and the potential far exceeds expectations.
The integrated seismic response of this type of SSFCRs indicates irregular reflection characteristics of waveforms, some of which are strong “beaded” reflections [
8,
9,
10,
11], and others of which are band-shaped messy reflections of different widths and lengths, with different energy levels. The strong “beaded” reflection is well-matched with high-yield oil and gas wells (oil production rate of about 100 tons/day) [
6,
7]. They are relatively easy to be identified in seismic profiles [
12,
13,
14,
15]. It is puzzling that with the deepening of exploration and development, the number of strong “beaded” reflections available for drilling is rapidly decreasing. There are in spite of no targets available for drilling. Therefore, there is an urgent need to develop new types of seismic reflection exploration. Now, researchers are exploring the SSFCRs, which are chaotic reflections with different widths and lengths. Drilling for this type of target encountered high-quality reservoirs, accompanied by drilling tool venting and drilling fluid leakage [
16,
17,
18]. The tested oil production level is higher than the strong “beading” reflection type targets, and the oil production rate of a batch of wells reached 500 tons per day [
19] and reached 1690 tons per day. However, it is difficult to effectively differentiate this type of reservoir using conventional attribute analysis based on depth migration processing data [
20], and it is a huge challenge to accurately locate the true underground location of this type of FVR. Accurately determining the true underground location of FVRs through geophysical methods is a problem that needs to be solved urgently.
Many scholars believe that this type of SSFCR is wrapped in a tight limestone formation, and each fractured-cavity body will cause diffraction when seismic waves propagate. Processors began to try to apply different migration methods to improve the resolution of heterogeneous reservoirs, such as the Kirchhoff migration [
21,
22,
23], Reverse time migration [
12,
24,
25], diffraction wave imaging [
26,
27], etc. However, no significant improvement has yet been achieved in the investigation area. Interpreters are also experimenting with root mean square amplitude attributes [
28], structural tensor attributes [
29,
30], amplitude gradient attributes, and curvature attributes [
31,
32,
33], etc. It can still effectively differentiate strong “beaded” reflection features, and it is still not obvious to effectively identify strip-shaped messy reflection features. In the investigation of the wave field anomaly identification method of ultra-deep fracture tunnel, it is considered to employ micro-seismic data for rock mass analysis of large underground caverns [
34,
35,
36], and jointly identify fracture tunnel bodies with seismic data. Unfortunately, the distance between wells in the investigation area is mostly one kilometer to five kilometers, or greater. There is no more micro-seismic data for the identification of fracture cavity bodies. Separation of characteristics of SSFCRs and reflection characteristics of stratum structure through interpretative processing methods, such as stratigraphic structure suppression and fault-solution morphological feature enhancement processing, shot domain pre-stack frequency division, and spectrum continuation [
27], can effectively highlight the characteristics of SSFCRs. This set of methods from data processing to interpretation runs through the entire process. The process control is complicated and repeated operations are time-consuming. Therefore, it is urgent to establish a set of efficient and effective methods or procedures for identifying the Seismic wave-field anomalies (SWFAs) of heterogeneous fractures and caves from the perspective of interpreters and effectively characterize them.
This paper proposes a new combination process, which is an SWFA identification method guided by structural dips of seismic data. We employed the method to highlight the SWFA characteristics of ultra-deep heterogeneous FVRs, which are called the SWFA data of heterogeneous FVRs. Empirically, this method can effectively identify the seismic response characteristics of uneven energy and irregular seismic waveforms. The band-like cluttered reflection features controlled by the fractured zone have an obvious recognition effect. Through the application in multiple actual blocks of FMOF, remarkable results have been achieved.
2. Methods
The proposed method in this paper was carried out in four major steps. The first step was to preprocess the data of seismic migration processing results to improve the signal-to-noise ratio (SNR) of the data and prepare for the following step of structural tilt and azimuth scanning. The second step was to scan the dip and azimuth of the seismic trace based on the preprocessed data to obtain interpretation data of the dip and azimuth of the formation and provide structural dip steering data for the next calculation. The third step was to perform structure-oriented transversal filtering on adjacent seismic traces on the preprocessed data, the purpose of which is to obtain the formation background data (FBD) in a local area. The fourth step was to subtract the result of the lateral filtering process from the original seismic data to obtain the residual, which was the SWFA data of heterogeneous FVRs encased in the dense formation.
2.1. Interpretative Preprocessing of Seismic Data
This step needs to be carried out according to the quality of seismic data. When the SNR is low, preprocessing, such as time-phase amplitude spectra analysis [
37], inverse Q-filter [
38], adaptive time-resampled high-resolution synchro squeezing transform [
39], band-pass filtering, three-dimensional data mixing, inter-track weighted averaging, and tilt scan stacking needs to be carried out. It can be a single pretreatment task or a combination of multiple pretreatment methods. The pre-processing method we employed is tilt-scan superimposition processing.
Inclination scan stacking enhances the coherence energy based on the similarity of the inclination and along-inclination data [
40,
41]. Normal spatio-temporal data can be decomposed into different tilt components using so-called tilt scan stacking or tilt stacking. For each inclination value, the samples are summed along the inclination angle and output as tau (time intercepts) and p inverse velocity (or inclination) values. Each track in the tau-p domain corresponds to a different inclination of the input data.
After converting to tau-p domain, the similarity of the data is determined as a measure of the similarity of the data along a certain inclination. The similarity is used to estimate and weight the p-channel samples, thereby enhancing the dominant tendency of the data. When the tau-p data is converted back to the t-x data, the principal coherent energy within a certain range of inclination will increase. After pre-processing, the SNR of the data has been significantly improved (
Figure 1).
2.2. Robust Inclination and Azimuth Scanning
Firstly, the method of estimating the vector inclination of 3D seismic data based on a multi-window is used to estimate the apparent inclination of the seismic data volume along the Inline and Cross-line [
42,
43,
44].
The quadratic surface is fitted by the similarity value of the pair of adjacent apparent inclination angles
, and the form is:
where
is the similarity,
is the apparent inline along the inline, and
is the apparent inclination along the crossline.
Solve the coefficient
using the method of least squares. The following equation is then used to obtain an estimate of the dip angle.
In this formula is the apparent inclination angle pair corresponding to the maximum value of the interpolated similar surface .
Even if this similarity surface is interpolated, the analysis window may inadvertently span a fault that divides the reflective layer into two regions with different inclination angles. In order to obtain an improved vector inclination estimated value, the multi-analysis window structure is used [
45]. In addition to the central window, a set of non-central and overlapping analysis windows are scanned. All windows contain analysis points of interest (
Figure 2). Then the amplitude variance of the
seismic trace falling in the
-th window is calculated:
where
represents the average value of
in the
-th analysis window. The window with the smallest variance will be assumed to be the best representative of the continuous reflective layer.
2.3. Structure Dip-Oriented Smoothing Filter
We constructed a dip-oriented smoothing filter, which is a transversal filter based on the scanning results of dip and azimuth. We performed horizontal seismic trace smoothing on the provided seismic data. The purpose was to apply filtering along the seismic horizon to eliminate random noise and enhance the lateral continuity of reflections through filtering. The key to constructing a guided smooth filter is to distinguish between the inclined azimuth and the noise superimposed of the reflecting layer. After estimating the dip and azimuth of the formation, a filter was used to enhance the signal along the reflector. The most commonly used filters are Mean filter and Median filter. The input data are the preprocessed seismic data, the apparent dip of the inline and the cross-line.
The Mean filter was the simplest filter to suppress random noise. On the plane, the Mean filter was a low-pass filter, which was a typical running window mean filtering technique. The output data value is the average value of all samples falling in the central analysis window. The window size is usually an odd number, such as 3 × 3, 5 × 5, or 9 × 9, or rectangular or elliptical. The Mean filter of time
is defined as:
In the formula, represents the th channel in the channel that falls into the analysis window.
The Median filter is one of the most widely used nonlinear techniques in signal and image processing. The Median filter replaces each sample on the seismic trace of the window with the median value of all samples that fall into the analysis window. It rejects dissident samples (outliers) during the execution process. The window size is usually odd. One method in the calculation is to use the ordered index (
) to sort all
samples in the analysis window:
Therefore, the median value can be denoted as:
2.4. Calculation of SWFA Information
The SWFA data of the heterogeneity SSFCRs was the difference between the seismic wave characteristics of the reservoir and the surrounding rock. Based on the original seismic data, the data generated after the structure-oriented lateral filtering was equivalent to a comprehensive response close to the seismic wave characteristics of the surrounding rock, which was regarded as the state before the formation of fractures [
23]. Therefore, the part deducting the seismic wave characteristics of the surrounding rock from the seismic data information was the seismic wave characteristics of the reservoir. The basic strategy is the original seismic data minus the structural dip to guide the laterally filtered data (Equation (4) or Equation (6)), which can keep the apparent phase of the heterogeneous attribute data volume consistent with the phase of the original seismic data. For this reason, the difference between the data after the structural dip guide filtering and the original seismic data is obtained to obtain a data volume of heterogeneous attributes (Equation (7)).
In the formula, represents the SWFA data. is original seismic data. is the laterally filtered data FBD.
2.5. Forward Modeling Analysis
In order to verify the effectiveness of the method, a comprehensive geologic model was established based on the geological characteristics of the investigation area, sonic logging data, VSP logging data, by considering the changes in lateral structures. The seismic data were acquired through numerical simulation, and the method and process were used to identify the anomalous information of the heterogeneous body wave field.
The geologic model considers the spatial variation of stratigraphic structure and the random distribution of geological bodies (
Figure 3a). The filling speed of the formation in the model varies uniformly from top to bottom as 4300 m/s, 4800 m/s, 4500 m/s, 5500 m/s, 5750 m/s~6000 m/s. The formation velocity of the geologic bodies in the model is filled randomly. Between 4600 m/s–5800 m/s, the speckled geologic bodies have the same shape, and the vertical and horizontal dimensions of the striped geological bodies are different.
A convolution algorithm was used for numerical simulation. By adjusting and adding noise parameters many times, the characteristics of the simulated profile are consistent with the actual seismic data when the factor is 0.01. The numerical simulation seismic profile is shown in
Figure 3b, and the seismic profile as indicated in
Figure 3c is the heterogeneous characteristic data calculated by this method. Compared with the model demonstrated in
Figure 3a, the SWFA specificity of the heterogeneous characteristic data, the position of the spot-shaped geological body and the strip-shaped geological body are distinguished correspondingly.
3. Actual Seismic Data Application
The investigation area is located in the northern part of the Tarim Basin in northwestern China (
Figure 4). It is a floodplain of the Tarim River. The terrain is relatively flat, and the surface elevation is 935–945 m. The underground structure is located in the Tabei Uplift-Northern Depression of the Tarim Basin. It is the slope part of the Lunnan Low Uplift that dips southward into the Northern Depression. Core slice analysis shows that the primary pores in the carbonate of the Yijianfang (TO
2y) and Yingshan (TO
1−2y) formations in the investigation area have almost disappeared. The secondary pores are chiefly secondary pores, accounting for more than 95% [
46,
47], and the reservoir space are mainly pores, holes and cracks formed by secondary dissolution and structural fracture.
The SSFCRs in the investigated area have undergone multiple faults and deep hydrothermal action, and the reservoir reformation is very strong, forming a complex karst fractured-cavity system with strong heterogeneity and complex and diverse reservoir types (
Figure 5a,b). This kind of carbonate reservoir demonstrates a reflection type characterized by low frequency and weak-strong amplitude composed of trough-peak or trough-peak-trough on the seismic profile, which is the so-called “beaded” reflection (
Figure 5c). Because karst is controlled by faults, the lateral distribution characteristics of the reservoir show the characteristics of strips spreading along the faults (
Figure 5d).
The application blocks of the actual data are selected in the SHB-YM, HD-YK and YGU blocks in the FMOF. The structural characteristics of the three blocks are different, including gentle structure, monocline and complex dome anticline. However, there is a common feature that a large number of strike-slip faults are developed, and high-quality reservoirs are concentrated near the faults. The three blocks are very suitable for carrying out the practical application and analysis of the method in this paper.
3.1. Application Example in SHB-YM Block
Taking the SHB-YM area as an example, one of the design purposes of the YM801-H6 well in the block is to drill the weak reflection controlled by the fractured zone. After entering the target formation, the vertical fracture zone is designed to drill horizontally along the formation. During the drilling of the well, 7327.50 m encountered drilling fluid loss, and the interval of 7433.09–7278.00 m encountered drilling tool venting. Subsequently, it was transferred to oil and gas test production operations, and perforation tests were carried out in the open hole section of 7278.00–7440.00 m, and finally high-yield industrial oil gas flow was obtained (
Figure 6). The maximum oil production reached 600 tons per day. After maintaining stable production for a period of time, water content was obtained. It rises in steps and then controls the smooth production of the nozzle (
Figure 7). In fact, the data profile of the best migration processing results implemented so far does not reveal the characteristics of large-scale reservoirs around the wellbore and cannot explain the characteristics of the YM801-H6 well, such as drilling and test production. That is, the seismic interpretation plan does not match the actual drilling and production.
The seismic profile along the YM801-H6 well, trajectory is characterized by band-like cluttered reflections of different widths, lengths, and energy strengths (
Figure 8a), but they are not typical. We extracted the Root mean square (RMS) attributes of the target layer in the well area along the layers and found that the attribute features have but no obvious plane distribution rules (
Figure 9a). In the results of the conventional seismic profile and attribute plane, the characteristics of a developed fault zone can only be vaguely displayed. The characteristics of weak anomalies are disturbed by the flake’s strong reflection background in the South and North. The wave group anomalies caused by these heterogeneous reservoirs are submerged in the background. This makes it impossible to effectively identify the characteristics of the FVRs.
In view of the above situation, we apply the method of this paper to carry out targeted research and calculate the heterogeneous body wave field heterogeneous information data body. From the well tie profile SWFA data section (
Figure 8b) and the target layer RMS amplitude attribute plan (
Figure 9b) extracted based on the SWFA data, the characteristics of the fault zone and the high-quality heterogeneity can be clearly distinguished. The new data reveals that develops a north-east-spanning fault zone with a width of 50–350 m (vertical principal displacement zone) in the area, which controls a large number of high-quality zonal-distributed reservoirs and provides a material basis for continuous fluid transport and energy replenishment. The SHB1-3 well, which is far away from the fault zone (about 500 m), has not drilled a good reservoir, and the production effect after reservoir reconstruction is not ideal. These clearly explain the high and stable production of the wells. The estimation results of underground volume are matched with the oil and gas output measured on the surface.
3.2. Application Example in HD-YK Block
In order to prove the effectiveness and stability of the method, it was applied to the HD-YK block of the research area. The geological characteristics of this block are more complex than those of the SHB-YM block. The target layer has a dip angle of 3–6 degrees, the maximum structure gradient is 15 degrees southeast, and the target layer has a structural drop of 450 m. The block was first discovered in 2015, and then five fractured-cavity units were depicted. The scheme design utilized seven appraisal wells and deployed 41 new development wells. Six wells and two test production wells were completed. In the past two years, researchers have been troubled by the fact that decision-makers cannot effectively compare and optimize well placement and drilling based on the characteristics of high-efficiency wells, resulting in a large gap between the number of wells and the planned total number of wells.
In the early stage, according to the plan explained based on the best migration processing data that has been implemented, there are only two fault zones with a length of 16.25 km that can deploy the effective well positions required by the plan, and no large-scale SSFCRs in the block can be demonstrated. The characteristics of the layers are the main reason for the lag in well deployment progress. From the previous interpretation results, we can conclude that the seismic profile across any line of the original seismic data can only distinguish the “beaded” reflections with strong reflected energy and multiple abnormal energies in the horizontal continuous track, and relatively clear SWFAs identified nine targets (
Figure 10a). As indicated in the RMS amplitude attribute plan, the characteristics of the fault zone are not clear. Coupled with the interference of the sheet-like strong reflection background, the diffraction wave field anomalies produced by many relatively small-scale heterogeneous reservoirs are uniform. Submerged in the background (
Figure 11a), it is difficult to delineate.
In view of the above situation, the method in this paper is used to carry out targeted research, and the data volume of the heterogeneous body wave field is obtained by calculation. The FVRs of the arbitrary profile using the heterogeneous SWFA data are clearly distinguishable, at least 14 sweet spots can be clearly identified, and the RMS amplitude attribute plan along the target layer can clearly distinguish the development of the fault zone. Distribution characteristics and heterogeneous reservoir spots (
Figure 10b and
Figure 11b), the length of the identified fault zone reached 42.66 km, and the number of well sites available for batch deployment exceeded 53.
3.3. Application Example in YGU Block
Take the YGU block as an example. The target stratum in the study area is an Ordovician TO
2y Formation–TO
12y formation. Ascribed to multi-directional compressive stress, the current structure of Ordovician stratum is a dome anticline structure (
Figure 12a), and the dip angle of TO
2y Formation varies greatly, ranging from 5°–35° (
Figure 12b). The reservoirs encountered during drilling are mainly caves, holes and fractures, with strong heterogeneity.
After years of exploration, researchers have always believed that the reservoir in the YGU block is controlled by a dome anticline structure and has a relatively unified oil-water interface. However, this view cannot solve three problems. First, the oil and gas development effect at the highest position of the structure is not the best. Second, there is no obvious relationship between the distribution of high-yield and high-efficiency wells and the height of the structure. Third, the wells located at the high position of the formation structure have rapid water cuts in production, while the wells located at the low position of the formation structure have no water production for a long time. At the same time, some researchers believe that the YGU block reservoir is controlled by strike slip faults, and the dome anticline structure only plays a role in the later adjustment and migration of oil and gas. This can solve the above three problems to a certain extent. However, this view has not been supported by direct evidence, such as earthquakes.
The method in this paper is applied to carry out targeted research, and the SWFA characteristics of the FVRs as indicated in the arbitrary profile are clearly demonstrated (
Figure 13). The RMS amplitude attribute plan along the target layer can clearly distinguish the distribution characteristics of FVRs along the fault zone (
Figure 14). Based on the research results of the heterogeneous SWFA data volume, the researchers deployed a wildcat well. This well is located in the southwest wing of the structure and is designed to drill the SSFCRs. It is 6741 m away from the main reservoir area. Vertically, 1100 m is added downward on the basis of the original unified oil-water contact, and the vertical depth reaches 7083.7 m. Drilling results: the logging interpretation of the target layer is 85.5 m/five layers of an oil layer and poor oil layer, and 9.5 m/one layer of a water bearing oil layer. Oil tested results: the open hole section 7133–7380 m, 4 mm nozzle, oil pressure is 17.35 MPa, daily oil production is 150 m
3 and natural gas is 2450 m
3. The view that the reservoir in the YGU block is controlled by faults is corroborated by this wildcat well. It is proved that the method has good applicability in complex structural areas.
This method has been applied in other blocks of FMOF, which can effectively improve the seismic wavefield abnormalities caused by fault-controlled heterogeneous fractures and caves and plays an important role in the prediction of FMOF reservoirs and guiding drilling into targets; guided by the effective improvement of the drilling success rate and the percentage of high-efficiency wells in the study area. The drilling success rate refers to the percentage of the number of industrial oil and gas wells obtained during the drilling test in the current year in the total number of wells, and the percentage of high-efficiency wells refers to the test production of the current year. The number of wells with an oil and gas equivalent of more than 50,000 tons accounted for the proportion of the total number of industrial oil and gas wells obtained from the oil test. According to the Tarim Oilfield, the drilling success rate has increased from 75% in 2017 to more than 95% in 2021, and the proportion of high-efficiency wells has increased from 20% in 2017 to 75% in 2021 (
Figure 15). Through the application of this method and process, many wells have achieved high production (
Table 1), and the current production state is stable. It supports the exploration of block reserves and promotes the deployment of 1130 km
2 high-density 3D seismic to the south of the HD-YK block (540 km
2 will be constructed first in the winter of 2021), with obvious effects.
According to the application effect of the three blocks, this method can be applied to areas with slow changes of structural characteristics or areas with a severe change of structural characteristics, which is obviously helpful to the identification of heterogeneous reservoirs. Through the method of dip scanning structure guided transverse filtering, the fault distribution characteristics in the study area can be effectively understood in the absence of structural horizon interpretation, and the initiative in time can be won for block geological understanding and well location deployment.
The method can be applied to depth domain and time domain stack data. It can also be used for reference by identifying geological abnormal bodies in gravity, electrical and magnetic exploration, and can also be applied to engineering geological disaster investigation.
4. Conclusions
The calculation method of the carbonate fault-controlled heterogeneous SWFA identification data calculation method proposed in this paper is implemented in four steps. The rapid dip angle scanning and the construction of the guided transverse filter are the key, and the more important step is to use the background removal method to strengthen the identification of wavefield anomalies caused by heterogeneous reservoirs. Of course, the application effect is closely related to the quality of seismic data, so it is essential to carry out data fidelity processing. It should be noted that the preprocessing of the seismic migration processing result data needs to be confirmed according to the seismic quality. If necessary, the preprocessing method is not limited to weighted mixing.
The application of this method effectively enhances the data’s ability to identify the seismic wavefield anomalies caused by the sub-level reflection intensity and the moderate-scale and small-scale carbonate heterogeneous reservoirs, which can further enhance the high-quality reservoirs. Effective fault zones and phenomena of high-yield of YM801-H6 and low yield of SHB1-3 can be explained (SHB-YM block). The discriminated fault zone is 42.66 km long, and the target number of well sites that can be deployed in the batch is 23 (HD-YK block). The abnormal geological bodies are distributed in strips, except for the main fracture area of the dome anticline structure (YGU block). Our application in the interpretation of real 3D seismic data has achieved very good results. It clearly explained some of the contradictory problems in geophysical prospecting, increased the number of targets available for drilling, and supported the exploration and development of reserves in the confirmed area.