In the Shunbei 1 and 5 zones of the study area, the Yijianfang Formation and the Yingshan Formation have 3 wells with conventional logging curves, such as density, P-wave velocity, gamma and resistivity curves, which can be used to analyze reservoir properties. Moreover, 31 wells provide information on the location of drilling fluid leakage or lost drilling tool, a critical factor for effective reservoir identification (Duan et al., 2020).

The observed seismic data in the study area has a bin size of 25 m × 25 m, a sample interval of 1 ms, and a dominant seismic frequency of 19 Hz. The seismic frequency is mainly distributed between 7 and 40 Hz. In order to more accurately depict the distribution of the reservoirs, this study applied a deep learning algorithm to interpret the faults. The main geological horizons are T₇ ⁴, T₇ ⁶, T₇ ⁸, T₈ ⁰, and T₈ ¹ (Figures 1B, D). In order to effectively identify the Shunbei fault-controlled karst reservoirs, a variety of reservoir identification methods have been improved and applied with satisfying results. In this section, the optimization methods are elaborated.

Generally, large caves form at the end of the dissolution process, as a result of the dissolution of tight limestone by surface water or deep hydrothermal fluids along faults. The longer the period of dissolution, the greater the potential for the formation of large caves. During tectonic evolution, some of these cavern reservoirs may be subsequently filled. Depending on the amount of filling, caverns can be further categorized as either unfilled or filled. Unfilled caves provide a high-quality reservoir space. When drilling this type of reservoir, drilling fluid leakage and drilling tools loss are possible. Under the influence of flowing water and gravity, certain caverns frequently fill with sediments, typically a complex mixture of sand and mud. The intergranular pore is the primary reservoir of the cave.

There are generally two types of seismic identification methods for large karst caves. The first is an attribute calculation-based identification method (Cheng et al., 2020), which directly extracts attributes from observed seismic data. This method is simple and convenient, and it fully respects the seismic reflection. However, it is also prone to errors in the location, number, and volume of identified caves due to the influence of the wavelet. The second method is an inversion of geological parameters, which can eliminate the wavelet effect to restore the cave position and achieve a more accurate cave volume (Hamid, et al., 2018; Li et al., 2020; Liu et al., 2020). However, for the Shunbei karst reservoirs with low dissolution degree, the conventional inversion of the karst cave reflection is easily confused with the strong stratigraphic reflection, resulting in low imaging accuracy (Figure 3). The inverted P-impedance in Figure 3B is obtained from the seismic data in Figure 3A. In Figure 3B, the bright yellow color indicates low impedance. Because karst caves are more porous than the surrounding rocks, their velocity and density are lower, resulting in a drop in P-impedance. The arrows in Figure 3B show that the caves are submerged due to the low impedance interference from the stratigraphic reflection, making it difficult to distinguish their shapes.

To solve the problem that conventional inversion methods are difficult to accurately characterize large caves, a new cave identification method is proposed in this study. The improvements cover three aspects. First, a low-pass filter term is constructed using the Hanning window function and Fourier transform. The new term can optimize the objective function of the inversion, thus improving the accuracy of the inversion results. Second, based on the joint constraint of well logging and seismic attributes, an improved low-frequency model is constructed to ensure accuracy and characterize the lateral heterogeneity features. This provides an optimization constraint for the improved inversion method. Third, a background interference reflection suppression calculation method is established to realize the extraction of background interference reflection based on the dynamic window, on which the nonlinear transformation equation is constructed to obtain high accuracy identification results of large caves.

Around caves and faults, it can be seen from the outcrop that there are associated dissolved pore and fracture zones. These zones are formed primarily by the dissolution and widening of small-to medium-sized fractures on either side of strike-slip faults. This type of reservoir consists of dissolved pore and fracture networks. Its distribution is controlled by strike-slip faults, which serve as both reservoir space and communication channel (Cheng et al., 2020). Therefore, for a comprehensive analysis of fault-controlled reservoirs, it is crucial to correctly identify this reservoir type. This section focuses on the identification of this type of reservoir and proposes a modified method to improve the accuracy of the identification.

Due to the influence of pores and fractures, the seismic reflection of dissolved pore and fracture zones typically displays anomalies. These anomalous responses can be amplified by the energy envelope ( F e n v e l o p e ). The equation is as follows: F e n v e l o p e = f 2 D + g 2 D (5)

In Eq. 5, f ∙ is the real part obtained by Hilbert transformation of seismic data, and g ∙ is the imaginary part.

Due to the strong influence of the sedimentary strata and the unconformity surface, it is impossible to precisely define the boundary of the dissolved pore and fracture zones. This can lead to interference responses, which may result in overestimation of the volume of this type of reservoir. In order to diminish the interference caused by the sedimentary strata and the unconformity surface, the background trend suppression method proposed in Section 4.2 is still utilized to optimize the seismic data and eliminate the interference. The phase shift transformation is then applied to restore the exact location of the characteristic reservoir. Based on Eq. 5, the energy envelope ( F e n v e l o p e ) reflecting the dissolved pore and fracture zone is then calculated. The exact expression is: F e n v e l o p e = f 2 s h i f t D − D t r e n d + g 2 s h i f t D − D t r e n d (6)

In Eq. 6, D t r e n d is the seismic background trend derived from seismic data D . s h i f t ∙ stands for the phase transformation calculation. f ∙ and g ∙ represent the real and imaginary components of the seismic data after background suppression, which is attained using Hilbert transform.

During the formation of strike-slip faults, tectonic stress and fluid dissolution result in the formation of fault-fracture zones consisting of micro-fractures. These percolation-capable zones form primarily around the main faults. The fault-fracture zones serve as vital flow channels for reservoirs. Their distribution characteristics have a significant impact on the extent of karst reservoirs. The seismic events of the fault-fracture zones are often chaotic and irregular (Figure 2). By precisely characterizing the fault-fracture zone, one can not only determine the extent of the fault-controlled reservoir, but also effectively constrain its interior description and reservoir modeling.

Interior fault characteristics of fault-controlled reservoirs are typically described using conventional fault attributes (e.g., coherence, curvature, etc.). For instance, the eigenvalue coherent is sensitive to discontinuous seismic events, and its recognition results are often linearly distributed, which is insufficient to accurately describe the distribution of fault-fracture zones. The structural tensor algorithm is extensively used to address the problem of fault-fracture zone identification in the Shunbei field (Li et al., 2020; Liu et al., 2020). In the case of high degree of dissolution, the detection performance of the gradient structure tensor is better. However, the region of the Shunbei field with low degree of dissolution is strongly influenced by strong seismic reflection events. The structural tensor data can only partially describe large karst caves with low dissolution, as demonstrated by the practical application. To address the problem of fault-fracture zone identification, this study proposes a method that improves the detection accuracy of fault-fracture zones and is applicable to fault-fracture bodies with different degrees of dissolution.

Figure 2 clearly demonstrates that the seismic reflection of the fault-fracture zone is disorganized. Texture is a measure that reflects the uniformity or smoothness of the target. Texture roughness can characterize disorganized features in terms of texture feature analysis. The coarser the texture, the more developed the fracture zone. In characterizing the spatial relationship and amplitude pattern of seismic data, the gray level co-occurrence matrix algorithm (Gao, 2011) can extract texture features from seismic data. In this research, this algorithm has been applied to the detection of fault-fracture zones.

The gray level co-occurrence matrix texture feature algorithm considers the seismic amplitude as the gray level of the image, and various amplitude values correspond to pixels of different gray levels (Gao, 2003). If the amplitude of seismic data with a given window range is divided into N different grades, it can be considered to have N different gray levels. The selection of N is primarily determined by the numerical distribution of the seismic amplitude and the required resolution. Consequently, each amplitude unit in 3D seismic data D a , b , c ( a , b , c are the spatial positions of the data point in the 3D seismic amplitude matrix) can be mapped to one of N different gray levels, each of which has a unique digital representation. For example, the gray level ranges from 0 to 15 when N is 16, allowing digitization of the grayscale matrix and generation of the 3D grayscale texture matrix G T a , b , c .

Based on the grayscale texture matrix G T a , b , c , the gray level co-occurrence matrix G i m , n i = x , y , z can be obtained along the x , y , and z axes. The size of the gray level co-occurrence matrix is N × N . For example, in the x-axis direction, the data in the m th row and n th column of G x m , n represent the number of times in the x-axis direction that a grayscale sample point m and a grayscale sample point n in the G T are adjacent. In other words, the gray level co-occurrence matrix records the adjacent frequency of each pair of grayscale data. The gray co-occurrence matrix G x m , n along the x-axis can be calculated using this method. Since the data points are adjacent, the calculated gray level co-occurrence matrix is axisymmetric, and the matrices G y m , n and G z m , n along the y-axis and z-axis can be obtained in the same way.

When the difference between m and n is larger, it implies that the heterogeneity of the seismic is more heterogeneous and the fault-fracture zone is more developed. In the principal diagonal of G i m , n i = x , y , z , the ability to characterize heterogeneity is weakest. The ability to characterize heterogeneity increases progressively from the principal diagonal to both.

In order to further optimize the ability of gray-level co-occurrence matrix to characterize heterogeneity, a grayscale texture contrast algorithm is presented in this research. On the basis of the calculated gray level co-occurrence matrix G i m , n i = x , y , z , the contrast weight coefficient ω = m − n is multiplied by each item at different positions in the matrix. In particular, at the principal diagonal position m − n = 0 , which represents adjacent sample points with the same gray level, the contrast weight coefficient is 0, indicating that there is no contrast. If the difference between m and n is 1, the contrast weight is 1. If the difference between m and n is 2, the contrast weight is 4. As the value of m − n increases, the contrast weight increases exponentially. The specific equations of the grayscale texture contrast g c o n t r a s t are: g c o n t r a s t = ∑ m , n = 0 N − 1 G i m , n R m − n 2 i = x , y , z (7) R x = 2 N y N z N x − 1 (8) R y = 2 N x N z N y − 1 (9) R z = 2 N x N y N z − 1 (10)

In above equations, N x , N y , and N z are the numbers of seismic data sampling points along the x , y , and z axes, respectively. R denotes the maximum probability of obtaining a specific data pair in the 3D grayscale texture matrix. The grayscale texture contrast value at the center point of the 3D seismic data can be obtained using Eq. 7.

There are still issues, although the fact that the grayscale texture contrast method can effectively characterize the distribution of fault-fracture zones. The interference is caused by the strong seismic reflection shielding effect of the T₇ ⁴ unconformity surface. Typically, the observed seismic reflection events consist of a group of compound harmonics. The effective reflection in the lower part of the reservoir is obscured by the strong seismic reflection of the unconformity interface. Consequently, it is challenging to identify the fault-fracture zone. In order to extract effective data from seismic data and suppress the shielding effect, we explored a set of methods for optimizing the observed seismic. Wavelet decomposition is a seismic attribute analysis technique based on waveform characteristics. It can transform seismic data into the frequency domain and decompose the original seismic data into a series of wavelets with varying shapes and dominant frequencies (Zhu et al., 2016). These series of wavelets can be considered as atoms of seismic data. Their equations are as follows: S e i s m i c t = ∑ A i W i f , α , φ + n o i s e (11)

In Eq. 11, A i denotes the amplitude of the i th atom W i ; f , α , φ are the frequency, central position and phase of each atom, respectively; n o i s e represents the interference noise.

Atoms can be extracted using Morlet wavelet, Gaussian function, Hamming window function, Hanning window function, etc. This paper employs the matching pursuit method (Sinha, et al., 2005; Yang et al., 2021). The matching pursuit method not only improves the resolution in the frequency domain, but also in the time domain. It is able to perform the atomic decomposition of seismic waves and extract the effective information from them, thereby reducing the shielding effect of strong reflection events. Texture contrast attribute for fault-fracture zone identification can be calculated based on strong seismic reflection suppression. However, due to a certain amount of noise interference after the strong seismic reflection suppression, the subsequent predictions still contain anomalous responses. The reason is that these interference noises are confused with fractures when calculating the fault-fracture zone, which leads to errors. To solve this problem, the interference noise was reduced using the Principal Component Analysis (PCA) algorithm after strong seismic reflection suppression. The fault-fracture zone was then calculated using the data that after strong seismic reflection suppression, PCA denoising and texture contract calculation.

To evaluate the effectiveness of the P-impedance inversion method proposed in Section 4.2.1, data from a research area in the Shunbei 1 Zone were selected (Figure 4A). Figures 4B, C compare P-impedance inversion results before and after improvement. The well point in Figure 4 is the measured P-impedance curve. In comparison, the conventional inversion results show obvious low-frequency banding effects (black circles in Figure 4B). However, for fault-controlled reservoirs with strong heterogeneity and few well data, this method reduces the ability to characterize reservoir heterogeneity and easily obscures details of lateral changes. If the reservoirs are identified based on these inversion results, there will be a significant error in reserve estimation. In contrast, the new method can more accurately reflect the characteristics of caves (Figure 4C). The new inversion results are more consistent with geological understanding. To quantify the new inversion method, a cross plot analysis of the measured P-impedance (y-axis data in Figures 4D, E) and the inverted P-impedance (x-axis data in Figures 4D, E) was performed. The measured data in Figures 4D, E are the logging data after filtering and resampling to the seismic frequency range. The inverted P-impedance is extracted from the pseudo wells. It can be seen that the distribution of the two is basically in the range of the 45° slope in Figure 4E. The correlation coefficient between measured and inverted data based on the new method is 0.86971 (Figure 4E), while the correlation coefficient based on the conventional method is 0.582894 (Figure 4D). The new inversion method provides reliable data for subsequent characterization of karst caves.

To verify the effectiveness of the proposed cave identification method, we constructed a 3D karst cave model based on the observed data from the Shunbei field (Figure 5A). The model includes feathers of karst caves and sedimentary strata. The arrows in Figure 5A indicate the locations of the karst caves. To create the synthetic seismic records, the reflection coefficient was calculated using this geological model and convolved with a wavelet with a dominant frequency of 25 Hz (Figure 5B). According to the synthetic records, caves developing at the top of the target layer are easily obscured by sedimentary strata (inside the black ellipse). Figure 5C illustrates the inverted P-impedance derived from the seismic in Figure 5B. The inverted results improve the resolution of the caves compared to the seismic reflection, but are still affected by the low impedance of the sedimentary strata. For instance, as shown in Figure 5C, the shallow low impedance stripe response interferes with the reflection of the caves within the red ellipse, making it difficult to discern the contour of the caves. When carving the cave with the threshold value, the response of the cave inside the yellow ellipse is more susceptible to removal than that of the shallow low-impedance strip. This problem may cause the following cave identification to be incorrect. Figure 5D presents the I P t r e n d as determined by the method described in Section 4.2.2. Figure 5E is the difference between the inverted I P (Figure 5C) and I P t r e n d (Figure 5D). The region of low-impedance in the figure indicates the karst caves. Figure 5F shows the caves ( I P M ) predicted by Eq. 4. Comparing Figures 5C, E, it is clear that the accuracy of cave identification in Figure 5E has been significantly improved, especially for the cave zones marked in the figure. Comparing Figures 5A, F, it is clear that the method proposed in this paper can accurately characterize the location of caves.

To further verify the applicability of the proposed method, an application area in the Shunbei field was selected. The Shunbei field has a lower degree of dissolution than the Tahe field, and its fault-controlled reservoirs are more developed. In addition, the cave reflection in the Shunbei field is weaker than that in the Tahe field, making it more susceptible to being confused with the reflection of sedimentary strata. Figure 6 shows a cross section of the Shunbei 5 Zone through a well. In the figure, the black curve represents the well trajectory. The green markers indicate the location of drilling fluid leakage, indicating the presence of reservoirs. Figure 6A shows the seismic reflection amplitude data. Figure 3B shows the conventional inverted results. Figure 6B shows the P-impedance obtained with the new inversion method. Comparing Figure 6B to Figure 3B, the new method provides more information and can more accurately represent horizontal heterogeneity. The bright red color in Figure 6B indicates a low-impedance value, revealing that the response of the caves is obscured by the low-impedance of the sedimentary strata (as indicated by the arrows in Figure 6B). And the exterior contours of the caves are difficult to identify. Figure 6C shows the cave results predicted by the new identification method. The results are in excellent agreement with the reservoir position revealed by the drilling data. In addition, the caves are concentrated around the main faults, and their size gradually decreases from the shallow to the deep target layer, which is consistent with geological theory.

The theoretical model test and the field case demonstrate that the novel method effectively suppresses background interference and more accurately describes the boundaries of caves. The characteristics of the cave distribution are consistent with geological knowledge, thus validating the feasibility and effectiveness of the proposed method.

A 3D dissolved pore and fracture equivalent model was added to the model in Figure 5A with reference to the measured data (Figure 7A). The red dashed line indicates the extent of the dissolved pore and fracture zone. Figure 7B shows the synthetic seismic. A comparison of Figure 7B and Figure 5B shows that when the dissolved pore and fracture zones are added, the seismic events become more chaotic and the amplitude energy varies to different degrees. Figure 7C shows the energy envelope derived from the seismic data in Figure 7B. Figure 7C shows that although the non-optimized energy envelope can identify the spatial location of the dissolved pore and fracture zones to some extent, there are still errors in the identification results (white arrows in Figure 7). For instance, the reservoir near 4,600 ms is affected by the upper lithologic interface and its boundary cannot be accurately characterized. In addition, the response of the reservoir near 4,650 ms is obscured by the interference of the sedimentary strata. In addition, the reservoir response near 4,700 ms is comparable to the energy level of the sedimentary strata near 4,650 ms. If a single threshold is used, these errors would lead to an inaccurate description of the distribution of dissolved pore and fracture zones in subsequent modeling. Figure 7D shows the results of the improved energy envelope. The range of dissolved pore and fracture zones derived by the new method is consistent with the theoretical model in Figure 7A, effectively reducing the interference of the lithologic interface and sedimentary strata. Thus, it is possible to characterize the spatial distribution range of 3D dissolved pore and fracture zones using a single threshold.

Figure 8A shows the seismic data of the Shunbei field, and the black part of the figure is the strike-slip fault identified based on the deep learning algorithm. Figure 8B shows the non-improved energy envelope derived from the seismic data shown in Figure 8A. Figure 8C shows the improved energy envelope. As shown in Figures 8B, C, it is obvious that the improved method can more accurately identify the distribution of cavity belts. In addition, the upper part of the reservoir has a higher degree of dissolution than the lower part. The longitudinal zonal distribution indicates that faults control the distribution of dissolved pore and fracture zones. Two leakage points are marked on the trajectory of Well C in Figure 8 (green lines), which is consistent with the predication results. The comparison shows that there is a good match between the dissolved pore and fracture zones and the location of the strike-slip fault in the study area, which is consistent with geological knowledge and validates the applicability of the proposed identification method.

Figure 9A shows the seismic amplitude reflection, while Figure 9B shows the identification results of the fault-fracture zone based on Eq. 7. Comparing the two figures, it can be seen that the proposed method can accurately delineate the fracture-damaged zone around the main faults. In the figure, the bright red color indicates areas where fractures are developed, while the blue color indicates areas where the fractures are not developed. Figure 9 demonstrates that the fracture-damaged zone has a typical funnel shape and that the intensity of the fractures progressively decreases from the main fault to the sides. However, there is till a problem, as shown in Figure 9B, interference is caused by the strong seismic reflection shielding effect of the unconformity surface (inside the black dashed circle in Figure 9B).

Based on the matched pursuit algorithm, we implemented a wave decomposition optimization process for seismic data, which is used to eliminate the interference of the strong seismic reflection. Figure 10 compares seismic data before and after processing to suppress the strong seismic reflection. An enlarge view of the unconformity interface is shown at the bottom of the figure. The comparison shows that the strong seismic reflection at the unconformity interface has been effectively suppressed and the precision of the target reservoir details has been effectively improved (at the yellow arrows in Figure 10).

The anomalous information resulting from the strong seismic reflection can be effectively suppressed (Figure 11A). However, the subsequent predictions still contain anomalous responses due to some interference noise after the strong seismic reflection suppression (indicated by yellow arrows in Figure 11B). The reason is that these interference noises are confused with fractures in the calculation of the fault-fracture zone, which leads to errors. To address this issue, the interference noise was reduced using the Principal Component Analysis (PCA) algorithm after strong seismic reflection suppression (Figure 11C). The fault-fracture zone (Figure 11D) was then calculated using the data that after strong seismic reflection suppression, PCA denoising and texture contract calculation.

Comparing Figures 11A, C, it is clear that PCA processing has successfully suppressed the noise. Similarly, Figures 11B, D show that the artifacts have been significantly reduced after PCA optimization (as indicated by the yellow arrows in the figures). In addition, a comparison of Figures 9B, 11D reveals that the method proposed in this study not only effectively weakens the shielding effect of strong seismic reflections, thereby improving the accuracy of the target characterization, but also preserves the true response of the fault-fracture zone.

To further verify the effectiveness of the fault-fracture zone identification method, a new model was created by adding the fault-fracture zone information (Figure 12A) to the model in Figure 7A, as shown in Figure 12B. Figure 12C shows the synthetic seismic. Comparing Figures 7B, 12C, it can be seen that the addition of fracture information makes the seismic events appear more chaotic. Figure 12D is a superimposed profile of the faults, predicted caves and fault-fracture zone. The seismic data in Figure 12C are used to predict fault-fracture zone. Figure 12E shows the synthetic seismic with optimization of the strong reflection suppression and PCA processing. By comparing Figures 12C, E, it can be seen that the optimized seismic effectively suppresses the effect of the strong reflection, thereby improving the accuracy of reservoir characterization (as indicated by the arrows in Figures 12C, E). Figure 12F is an overlay profile of the faults, identified caves and fault-fracture zone. The predicted fault-fracture zone is derived from the optimized seismic data of Figure 12E. The light purple color indicates the developed part of the fault-fracture zone, while the blue color indicates the undeveloped part. Comparative analysis reveals that the optimized method eliminates the influence of the interface (in the red dashed boxes in Figure 12). The detection results show a strong correlation with faults and karst caves, demonstrating that the new method accurately describes the extent of fault-fracture zones.

To further support the applicability of the fault-fracture zone identification method presented in this study, we select another field data set for testing. Figure 13 illustrates the results in a cross section that extends over 6 wells in a northeast to southwest orientation. Figure 13A shows the original seismic data. Figure 13B shows the predicted fault-fracture zone based on the modified method. The results can accurately capture the spatial distribution range of fault-fracture zones. The predicted fracture body is consistent with the geological comprehension and reservoirs indicated by the drilling data.

The fine identification of multi-level fracture-controlled karst reservoirs can be facilitated by the proposed integrated methods. Using multi-level characterization methods, the spatial distribution of the fracture-controlled karst reservoirs in the Shunbei field can be effectively analyzed in three dimensions. A series of optimization processing, inversion and interpretation techniques are applied to the seismic data (Figure 14A) to identify different levels of fault-controlled karst reservoir zones, such as large cave zones, dissolved pore and fracture zones, and fault-fracture zones (Figure 14B). As shown in the profile of the characterization results (Figure 14B), the reservoirs are mainly developed around the main fault zones (Figure 14C), which are significantly controlled by strike-slip faults. The shallow parts of the reservoirs are more developed in terms of faults and fractures.

It is clear from the map diagram (Figure 15) that the karst caves have a relatively limited extent and were formed mainly around the main faults. The dissolution pore and fracture zones are mainly present in the area surrounding karst caves and extent in the direction of the faults. The fault-fracture zones are the most developed, allowing the main faults and their associated fractures to be clearly characterized, as well as the extent of the reservoirs.