Earthquake identification method and system for marlstone crack hole type reservoir

Abstract

The invention provides an earthquake identification method and system for a marlstone fracture hole type reservoir. The method comprises the following steps: acquiring three-dimensional earthquake data, a logging curve and a logging reservoir interpretation result; determining the position of a target layer according to the three-dimensional seismic data, the logging curve and the logging reservoir interpretation result; and performing coherent enhancement filtering based on the three-dimensional seismic data to obtain coherent enhancement filtering data. According to the method and the system, background reflection of seismic data is removed step by step, reflection characteristics of the fracture-hole type reservoir are reserved, finally, longitudinal and transverse distribution of the marlstone fracture-hole type reservoir is obtained by utilizing seismic attributes, the method and the system have few human factors in the implementation process, and the implementation efficiency is high. The method is suitable for fine prediction of the small-scale marlstone fracture-hole type reservoir, is high in operability, and provides technical support for predicting longitudinal and transverse distribution of the marlstone fracture-hole type reservoir and finding a favorable sweet spot area.

Description

Technical Field

[0001] The present invention belongs to the technical fields of oil and gas field exploration and development, reservoir characterization and reservoir description, and particularly relates to a seismic identification method and system for marl fracture-hole reservoirs. Background Art

[0002] As a low-carbon, clean energy source, natural gas is expected to contribute an increasing share of the energy mix in the future. Therefore, domestic oil and gas companies must continuously discover new natural gas reserves. Recent breakthroughs have been made in natural gas exploration within the Sichuan Basin's marls, including the discovery of unconventional marl reservoirs. Research has revealed that the marl gas reservoirs are self-generating and self-storing, with both the source rock and reservoir being organic-rich marls with low porosity and permeability. Gypsum-salt layers develop at both the top and bottom of the marls. Due to subsequent compression and rubbing, the gypsum-salt rocks have developed numerous fractures and associated cracks. These fracture and fracture systems can significantly improve the storage properties of the marl reservoirs, resulting in fracture-cavity reservoirs modified by these fractures. Exploration practice has shown that only when porous marl reservoirs are superimposed on fracture-cavity reservoirs can industrially produce gas reservoirs.

[0003] Fracture-vuggy reservoirs are currently primarily identified using seismic attributes, including root mean square amplitude, amplitude change rate, and curvature. RMS amplitude primarily reflects amplitude anomalies and is therefore highly sensitive to large amplitudes. While suitable for detecting large and medium-scale cavernous reservoirs with distinct "beaded" reflections, it struggles with small-scale fracture-vuggy reservoirs.

[0004] The amplitude change rate attribute is based on the root mean square amplitude attribute and is related to the lateral variation of seismic amplitude. When large fracture-vuggy reservoirs exist in carbonate rocks, the amplitude will vary significantly. The amplitude change rate attribute can detect the location of fracture-vuggy reservoirs. The amplitude change rate technology is mainly suitable for predicting reservoir development zones dominated by caves.

[0005] The core of coherence volume technology is to convert seismic amplitude data into correlation coefficient data through seismic trace similarity analysis, highlighting uncorrelated anomalies. The development of fracture-cavity reservoirs affects the similarity of seismic traces, causing changes in the occurrence of seismic events, resulting in significant incoherence in the coherence attributes. This technology is primarily used to detect large faults and caves.

[0006] Curvature technology reflects the degree to which seismic events deviate from a straight line or plane, describing the degree of curvature at any point on the seismic events. Attribute calculation first converts the seismic data volume into a dip data volume and an azimuth data volume. The curvature of any position in space is then calculated based on the dip and azimuth angles to form a curvature data volume. Curvature attributes are commonly used to detect faults and fractures, and can identify finer fracture information than coherence, as well as some large karst structures. Because curvature attributes are primarily influenced by the curvature of seismic events, locations with larger event dips on seismic observations are easily identified as anomalies. Furthermore, most fracture-cavity reservoirs exhibit increased amplitude on seismic observations, with relatively small changes in event curvature. Therefore, curvature attributes make it difficult to identify these fracture-cavity reservoirs using these attributes.

[0007] Other related amplitude and structural attributes are effective in identifying pore and cave reservoirs and can identify the spatial distribution characteristics of large-scale fracture-cavity systems. However, they suffer from relatively simple methods, limited practicality, and a large identification scale. Unlike the typical "beaded" reflection cave reservoirs in the Tarim Basin's carbonate rocks, the discovered marl fracture-pore reservoirs in the Sichuan Basin are smaller in scale and heavily infilled with calcite and other materials. The difference in wave impedance between the reservoir and the surrounding rock is relatively small, and there are basically no obvious beaded reflections on the seismic profile. Instead, they appear as chaotic reflections with slightly higher energy, mainly developed near faults and significantly controlled by faults. Based on existing seismic fracture-cavity prediction technology, it is difficult to achieve a precise prediction of small-scale marl fracture-pore reservoirs. Summary of the Invention

[0008] In order to solve at least one problem in the background technology, the present invention provides a seismic identification method and system for marl fracture-vuggy reservoirs.

[0009] In order to achieve the above object, the present invention adopts the following technical solutions:

[0010] A seismic identification method for marl fracture-vuggy reservoirs comprises the following steps:

[0011] Obtain 3D seismic data, well logging curves and logging reservoir interpretation results;

[0012] Determine the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results;

[0013] Perform coherence enhancement filtering based on three-dimensional seismic data to obtain coherence enhancement filtering data;

[0014] Based on the coherent enhancement filter data and the target layer, multiple seismic wavelet component data are obtained by wavelet decomposition, background component data are removed, and all remaining seismic components are reconstructed to obtain seismic data of the prominent fracture-hole type reservoir;

[0015] Based on the seismic data of the fracture-pore type reservoir, robust principal component analysis is used to highlight the longitudinal and lateral changes of the seismic amplitude, eliminate the background reflection of the surrounding rock, and obtain the seismic reflection data of the fracture-pore type reservoir;

[0016] Extracting sweet spot attributes from the data after the robust principal component analysis to obtain a seismic sweet spot data volume;

[0017] The seismic sweet spot data volume is calibrated using drilling data, a threshold of a fracture-pore type reservoir is determined, and a sweet spot attribute data volume reflecting the fracture-pore type reservoir is obtained.

[0018] Preferably, the three-dimensional seismic data includes seismic data collected in the field by arranging sources and detectors in an area manner, and three-dimensional pre-stack time migration seismic data obtained based on indoor seismic data processing;

[0019] The logging curves include acoustic logging curves, density logging curves and natural gamma ray spectrum curves;

[0020] The logging reservoir interpretation results include reservoir interfaces, thickness and physical property parameters obtained through carbonate rock optimization processing of different logging curves.

[0021] Preferably, determining the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results includes the following steps:

[0022] Generate wave impedance curve according to acoustic logging curve and density logging curve;

[0023] Extract wavelets from 3D seismic data and determine seismic wavelets;

[0024] The wave impedance curve is convolved with the seismic wavelet to determine the synthetic logging record;

[0025] Calibrate the synthetic records with the 3D seismic data to establish the corresponding relationship between the seismic records and the well data, and determine the position of the top and bottom boundaries of the target layer;

[0026] Conduct horizon interpretation on the seismic data of the target layer to obtain the top and bottom boundaries of the target layer.

[0027] Preferably, performing coherence enhancement filtering based on three-dimensional seismic data to obtain coherence enhancement filtered data includes the following steps:

[0028] Select an analysis position and obtain the inclination data volume and azimuth data volume of the analysis position;

[0029] Acquire a coherent attribute data volume based on the inclination data volume and the azimuth data volume;

[0030] A coherence enhancement filtering data volume is obtained based on the inclination data volume, the azimuth data volume and the coherence attribute data volume.

[0031] Preferably, the coherence enhancement filter equation is as follows:

[0032]

[0033] Where x, y, and t represent the three axes of the spatial coordinate system, t represents the time axis, and u(x, y, t) represents the three-dimensional seismic data volume. is the directional derivative vector of the three-dimensional seismic data volume, D is the diffusion tensor, ∈ is the continuous factor, and the continuous factor ∈ is as follows:

[0034]

[0035] Where S0 represents the initial structure tensor matrix, S ρ Represents the gradient structure tensor at the current number of iterations; Tr represents the trace of the matrix, which is the sum of the matrix diagonal elements; the value range of ∈ is [0, 1], 0≤∈≤1.

[0036] Preferably, based on the coherent enhancement filter data and the target layer, wavelet decomposition is used to obtain multiple seismic wavelet component data, background component data is removed, and all remaining seismic components are reconstructed to obtain seismic data of a prominent fracture-cavity reservoir, including the following steps:

[0037] Based on the target layer position, the seismic data at the target layer position is subjected to seismic wavelet decomposition, and each seismic data is decomposed into multiple seismic wavelet component data of different shapes and frequencies;

[0038] Screening out background seismic components that do not contain fracture-cavity reservoir information from seismic wavelet components;

[0039] The background components that do not contain fracture-cavity reservoir information are removed from the seismic wavelet components, and all the remaining seismic wavelet components are reconstructed to obtain the seismic data volume of the prominent fracture-cavity reservoir.

[0040] Preferably, among the multiple seismic wavelet component data of different shapes and frequencies, the first component represents the seismic wavelet component with the greatest commonality in the seismic data, the second component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component, the third component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component and the second component, and so on.

[0041] Preferably, based on the seismic data of the prominent fracture-pore type reservoir, robust principal component analysis is used to highlight the longitudinal and lateral changes of the seismic amplitude and eliminate the background reflection of the surrounding rock to obtain the seismic reflection data of the fracture-pore type reservoir, including:

[0042] min D,E ||A|| * +γ||E|| * ;

[0043] where ||A|| * represents the nuclear norm of matrix A; ||E|| * represents the l1 norm of the matrix E; A is a low-rank matrix, regarded as the background seismic reflection data; E is a sparse matrix, which is the seismic reflection data of the prominent fracture-pore type reservoir; γ is the balance parameter.

[0044] Preferably, extracting sweet spot attributes from the data after the robust principal component analysis to obtain a seismic sweet spot data volume includes:

[0045]

[0046] Where A(t) is the seismic signal reflection intensity; f(t) is the instantaneous frequency; and S(t) is the seismic sweet spot attribute.

[0047] Preferably, the seismic sweet spot data volume is calibrated using drilling data, a threshold of a fracture-pore type reservoir is determined, and a sweet spot attribute data volume reflecting the fracture-pore type reservoir is obtained, comprising the following steps:

[0048] Calibrate the seismic sweet spot data volume based on the fracture-cavity reservoir interpreted by well logging, and determine the minimum value of the sweet spot attribute corresponding to the fracture-cavity reservoir;

[0049] The minimum value part in the seismic sweet spot data volume is filled with invalid values, and only the part of the data greater than the minimum value is retained. The part greater than the minimum value represents the distribution of fracture-pore type reservoirs, and the fracture-pore type reservoir prediction data volume is obtained.

[0050] A seismic identification system for marl fracture-vuggy reservoirs, comprising:

[0051] Basic unit, used to obtain 3D seismic data, well logging curves and logging reservoir interpretation results;

[0052] A segmentation unit is used to determine the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results;

[0053] A calculation unit, configured to perform coherence enhancement filtering based on the three-dimensional seismic data to obtain coherence enhancement filtering data;

[0054] A component unit is used to obtain multiple seismic wavelet component data by wavelet decomposition based on the coherent enhancement filter data and the target layer, remove background component data, reconstruct all remaining seismic components, and obtain seismic data of a prominent fracture-pore type reservoir;

[0055] An analysis unit is configured to use robust principal component analysis based on the seismic data of the prominent fracture-pore type reservoir to highlight the longitudinal and lateral changes of the seismic amplitude and eliminate the background reflection of the surrounding rock to obtain the seismic reflection data of the fracture-pore type reservoir;

[0056] an extraction unit, configured to extract sweet spot attributes from the data after the robust principal component analysis to obtain a seismic sweet spot data volume;

[0057] The identification unit is used to calibrate the seismic sweet spot data volume using drilling data, determine the threshold of the fracture-pore type reservoir, and obtain the sweet spot attribute data volume reflecting the fracture-pore type reservoir.

[0058] Preferably, the segmentation unit comprises:

[0059] A curve processing module is used to generate a wave impedance curve based on the acoustic logging curve and the density logging curve;

[0060] The wavelet extraction module extracts wavelets from 3D seismic data and determines seismic wavelets;

[0061] The synthesis module is used to perform convolution processing on the wave impedance curve and the seismic wavelet to determine the synthetic logging record;

[0062] The calibration module is used to calibrate the synthetic records and 3D seismic data, establish the corresponding relationship between the seismic records and the well data, and determine the position of the top and bottom boundaries of the target layer;

[0063] The segmentation module is used to interpret the seismic data of the target layer segment and obtain the top and bottom boundary layers of the target layer.

[0064] Preferably, the calculation unit includes:

[0065] The first calculation module is used to select an analysis position and obtain a dip data volume and an azimuth data volume of the analysis position;

[0066] A second calculation module is used to obtain a coherent attribute data volume based on the inclination data volume and the azimuth data volume;

[0067] The third calculation module is used to obtain a coherence enhancement filter data volume based on the inclination data volume, the azimuth data volume and the coherence attribute data volume.

[0068] Preferably, the analysis unit comprises:

[0069] The first analysis module is used to perform seismic wavelet decomposition on the seismic data at the target layer position based on the target layer position, and decompose each seismic data into a plurality of seismic wavelet component data of different shapes and frequencies;

[0070] The second analysis module is used to filter out background seismic components that do not contain fracture-pore type reservoir information from the seismic wavelet components;

[0071] The third analysis module is used to remove the background components that do not contain fracture-cavity reservoir information from the seismic wavelet components, reconstruct all remaining seismic wavelet components, and obtain a seismic data volume that highlights the fracture-cavity reservoir.

[0072] Preferably, the identification unit includes:

[0073] A threshold determination module is used to calibrate the seismic sweet spot data volume based on the fracture-pore reservoir interpreted by well logging, and determine the minimum value of the sweet spot attribute corresponding to the fracture-pore reservoir;

[0074] The identification module is used to fill the minimum value part in the seismic sweet spot data body with invalid values, and only retain the part of the data that is greater than the minimum value. The part greater than the minimum value represents the distribution of fracture-pore type reservoirs, and obtain the fracture-pore type reservoir prediction data body.

[0075] Beneficial effects of the present invention:

[0076] The method and system of the present invention gradually remove the background reflection of seismic data, retain the reflection characteristics of the fracture-pore reservoir, and finally use seismic attributes to obtain the vertical and horizontal distribution of the marl fracture-pore reservoir. The implementation process of this method and system has few human factors, is suitable for the precise prediction of small-scale marl fracture-pore reservoirs, and has strong operability. It provides technical support for predicting the vertical and horizontal distribution of marl fracture-pore reservoirs and finding favorable sweet spots.

[0077] Other features and advantages of the present invention will be described in the following description, and in part will become apparent from the description, or will be understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures pointed out in the description and the drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0079] Figure 1 This is a flow chart of a method for predicting marl fracture-cavity reservoirs according to an embodiment of the present invention;

[0080] Figure 2 is a seismic profile according to an embodiment of the present invention;

[0081] Figure 3 is a seismic profile after coherence enhancement filtering processing according to an embodiment of the present invention;

[0082] Figure 4 is a seismic profile after wavelet decomposition and reconstruction processing according to an embodiment of the present invention;

[0083] Figure 5 is a seismic profile processed by robust principal component analysis according to an embodiment of the present invention;

[0084] Figure 6 is a sweet spot attribute profile of data processed by robust principal component analysis according to an embodiment of the present invention;

[0085] Figure 7 This is a plan view of a dessert slice along a layer according to an embodiment of the present invention;

[0086] Figure 8 This is a block diagram of a seismic identification system for marl fracture-vuggy reservoirs according to the present invention. DETAILED DESCRIPTION

[0087] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0088] A seismic identification method for marl fracture-cavity reservoirs, such as Figure 1 As shown, the following steps are included:

[0089] S1: Obtain 3D seismic data, well logging curves, and logging reservoir interpretation results as input data; S2: Determine the target layer based on the 3D seismic data, well logging curves, and logging reservoir interpretation results; S3: Perform coherence enhancement filtering on the 3D seismic data to highlight the boundary characteristics of the fracture-pore reservoir and obtain coherence enhancement filtering data; S4: Based on the coherence enhancement filtering data and the target layer, use wavelet decomposition to obtain multiple seismic wavelet component data, remove background component data, and reconstruct all remaining seismic components to obtain seismic data highlighting the fracture-pore reservoir; S5: Based on the seismic data highlighting the fracture-pore reservoir, use robust principal component analysis to perform nonlinear dimensionality reduction to highlight the longitudinal and lateral changes in seismic amplitude, further eliminate the background reflection of the surrounding rock, and obtain seismic reflection data of the fracture-pore reservoir; S6: Extract sweet spot attributes from the data after robust principal component analysis to obtain a seismic sweet spot data volume; S7: Use drilling data to calibrate the seismic sweet spot data volume, determine the threshold of the fracture-pore reservoir, and obtain a sweet spot attribute data volume reflecting the fracture-pore reservoir.

[0090] It should be noted that the main advantage of this method is that it utilizes multiple background separation methods to gradually remove background reflection interference, highlighting the reflection characteristics of marl fracture-vuggy reservoirs. This method is more targeted than traditional methods that rely on single amplitude and structural attributes. Therefore, it can predict fracture-vuggy reservoirs with smaller seismic reflection amplitude differences and produce more accurate results. This method is less prone to human factors and is suitable for the precise prediction of small-scale marl fracture-vuggy reservoirs. It also has strong operability. This method has the following four advantages: (1) Coherence enhancement filtering is used to highlight the discontinuity of seismic reflection and the boundaries of anomalies, wavelet decomposition is used to further weaken the seismic background reflection characteristics, highlight the seismic reflection of fracture-pore reservoirs, and robust principal component analysis is used to further remove the interference of reflection background, retaining only the reflection characteristics of fracture-pore bodies. Finally, the sweet spot attributes are used to clearly reflect the structural characteristics and boundary information of fracture-pore reservoirs. This method process gradually removes background reflections and obtains the vertical and horizontal distribution of fracture-pore reservoirs; (2) Wavelet decomposition is used to achieve optimal matching and effective separation of strong reflection background, weaken the shielding phenomenon of background reflection on fracture-pore reservoirs, improve the recognition accuracy of fracture-pore reservoirs, and effectively eliminate the influence of strong reflection background on fracture-pore reservoir imaging. (3) Robust principal component analysis addresses the problem of scattered seismic data features and non-Gaussian noise. Based on a basic model with sparse and low-rank constraints, the model can be converted into a convex optimization problem for accurate solution, and the solved sparse matrix is ​​used to obtain seismic data of fracture-pore type reservoirs; (4) This method is not only for the identification and prediction of marl fracture-pore type reservoirs, but can also be extended to the identification of geological anomalies such as carbonate fracture-cavity type, cave type, karst type and bioherm type, and the identification accuracy is higher than that of traditional seismic face method and single attribute method. It can detect small-scale anomalies on seismic profiles that are difficult to identify with the naked eye, thereby greatly improving the interpretation accuracy of seismic profiles.

[0091] Steps S1 to S7 are further described below.

[0092] In step S1, 3D seismic data includes field data collected by deploying sources and receivers in an area-based manner, as well as 3D prestack time migration seismic data obtained through indoor seismic data processing. Well logs include sonic logs, density logs, and natural gamma ray spectrometers. Well log reservoir interpretation results include reservoir interfaces, thicknesses, and physical properties derived from carbonate rock optimization processing of various well logs.

[0093] It should be noted that indoor seismic data processing refers to the process of processing and transforming raw seismic data collected in the field using a seismic data processing system to obtain information about underground geological structures and formation properties. Furthermore, carbonate rock optimization processing is a versatile well logging data interpretation method. It uses a structure that is independent of the model and well logging combination to establish an error model function between the measured and theoretically calculated well logging values:

[0094]

[0095] Among them LG i 、TLG i where represents the actual and theoretically calculated values ​​for the i-th well log, respectively; m represents the number of well logs used. By finding the solution that minimizes the error model function, we obtain a parameter curve related to the reservoir. This parameter curve can be used to determine the reservoir interface, thickness, and physical properties.

[0096] In step S2, the following steps are included:

[0097] S201: Generate a wave impedance curve based on the acoustic logging curve and the density logging curve.

[0098] S202: extracting wavelets from the 3D seismic data to determine seismic wavelets;

[0099] S203: Convolution processing is performed using the wave impedance curve and the seismic wavelet to determine the well logging synthetic record.

[0100] S204: Calibrate the synthetic records and the 3D seismic data to establish a corresponding relationship between the seismic records and the well data, and determine the positions of the top and bottom boundaries of the target layer;

[0101] Horizon interpretation is performed on the seismic data of the target layer segment (the segment from the top boundary to the bottom boundary of the target layer) to obtain the top boundary and bottom boundary of the target layer.

[0102] In step S3, the following steps are included:

[0103] S301: Select an analysis position and obtain a dip data volume and an azimuth data volume of the analysis position; specifically, use multi-window scanning to obtain the most similar adjacent window as the dip and azimuth of the analysis position to obtain the dip data volume and the azimuth data volume.

[0104] S302: Obtain a coherent attribute data volume based on the inclination data volume and the azimuth data volume; specifically, calculate the coherent attribute using the inclination data volume and the azimuth data volume to obtain the coherent attribute data volume. The coherent attribute is calculated as follows:

[0105]

[0106] In the formula, λ max is the maximum eigenvalue of the seismic data matrix; m is the data sample point, and the relationship between the eigenvalues ​​of the data matrix is ​​used to detect the discontinuity of the formation; λ j is the j-th eigenvalue of the earthquake data matrix.

[0107] A coherence enhancement filtering data volume is obtained based on the dip data volume, the azimuth data volume and the coherence attribute data volume; specifically, the dip data volume and the azimuth data volume are used to perform directional filtering along the stratum, and the coherence attribute results are used to judge the continuity of the stratum. Discontinuous strata are not smoothed, while continuous strata are smoothed, thereby achieving the purpose of highlighting discontinuities and obtaining a coherence enhancement filtering data volume.

[0108] The coherence enhancement filter equation is as follows:

[0109]

[0110] Among them, x, y, t represent the three axes of the spatial coordinate system, t represents the time axis, and u(x, y, t) represents the three-dimensional seismic data volume. is the directional derivative vector of the three-dimensional seismic data volume, D is the diffusion tensor, ∈ is the continuous factor, and the continuous factor ∈ is as follows:

[0111]

[0112] Where S0 represents the initial structure tensor matrix, S ρ Represents the gradient structure tensor at the current iteration number; Tr represents the trace of the matrix, which is the sum of the matrix diagonal elements; the value range of ∈ is [0, 1], 0≤∈≤1; near special structures such as faults and caves, ∈≈0, and in wide and gentle areas away from faults, ∈≈1.

[0113] In step S4, the following steps are included:

[0114] S401: Based on the target layer position, seismic wavelet decomposition is performed on the seismic data at the target layer position, and each seismic data is decomposed into multiple seismic wavelet component data of different shapes and frequencies; the first component represents the seismic wavelet component with the greatest commonality in the seismic data, the second component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component, and the third component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component and the second component, and so on, and different components reflect geological seismic characteristics of different scales;

[0115] S402: utilizing the characteristics of fracture-pore type reservoirs, which have strong amplitudes and disordered events, to filter out background seismic components that do not contain information about fracture-pore type reservoirs from seismic wavelet components;

[0116] S403: removing background components that do not contain fracture-vuggy reservoir information from the seismic wavelet components, and reconstructing all remaining seismic wavelet components to obtain a seismic data volume that highlights fracture-vuggy reservoirs.

[0117] In step S5, it includes:

[0118] S501: Assume that the seismic wavelet reconstruction data matrix is ​​D, D has a low-dimensional structure space and can be expressed as the sum of two matrices:

[0119] D=A+E;

[0120] The matrix A has a certain amount of internal structural information, which causes the rows or columns to be linearly correlated and is therefore low-rank; E contains noise and is sparse. The goal of robust principal component analysis is to find the lowest-rank matrix A and the matrix E with the fewest nonzero elements, that is, to solve the following convex optimization problem:

[0121] min D,E ||A|| * +γ||E|| * ;

[0122] Where γ is the equilibrium parameter; ||A|| * represents the nuclear norm of matrix A; ||E|| * represents the l1 norm of the matrix E; A is a low-rank matrix, which is regarded as the background seismic reflection data; E is a sparse matrix, which is the seismic reflection data of the prominent fracture-pore type reservoir.

[0123] The above-mentioned augmented Lagrange multiplier method is also called the alternating direction method. This method has the characteristics of high accuracy and fast convergence speed. It has now become the mainstream algorithm for solving sparse and low-rank models.

[0124] Step S6 includes:

[0125] The sweet spot property is the reflection intensity divided by the square root of the instantaneous frequency. Specifically, let the seismic signal be x(t), and the reflection intensity expression is as follows:

[0126]

[0127] in is the Hilbert transform of the seismic signal x(t); t is time.

[0128] The instantaneous phase expression is as follows:

[0129]

[0130] The instantaneous frequency is as follows:

[0131]

[0132] The dessert properties are calculated as follows:

[0133]

[0134] Where A(t) is the seismic signal reflection intensity; f(t) is the instantaneous frequency; and S(t) is the seismic sweet spot attribute.

[0135] In step S6, the following steps are included:

[0136] The seismic sweet spot data volume is calibrated based on the fracture-cavity reservoir interpreted by well logging, and the minimum value of the sweet spot attribute corresponding to the fracture-cavity reservoir, i.e., the threshold value, is determined;

[0137] The part of the seismic sweet spot data volume that is smaller than the threshold is filled with invalid values, and only the part of the data that is larger than the threshold is retained. The part larger than the threshold represents the distribution of fracture-pore type reservoirs, and the fracture-pore type reservoir prediction data volume is obtained.

[0138] It should be noted that the 3D seismic data of the central Sichuan Basin is used as input. Figure 1 The process shown is used to predict and analyze the marl fracture-pore type reservoir. The main marl fracture-pore type reservoirs are developed in this area. Figure 2-Figure 6 In the figure, the vertical axis is time and the horizontal axis is the track number. Figure 2 This is a seismic profile of an embodiment of the present invention; since the fracture-cavity reservoir is small in scale and difficult to distinguish from background reflections in seismic terms, it is extremely challenging to directly use raw seismic data to predict fracture-cavity bodies. Figure 3 This is a seismic profile after coherence enhancement filtering processing according to an embodiment of the present invention. Coherence enhancement filtering has optimal smoothing characteristics and edge preservation characteristics. Compared with the original data, the coherence enhancement filtered data highlights the discontinuity of the seismic phase axis and provides a richer and more detailed display of geological and structural anomalies. Figure 4 This is a seismic profile after wavelet decomposition processing according to an embodiment of the present invention. Seismic waveform decomposition is used to decompose each seismic trace into multiple seismic wavelets of different shapes and frequencies. Through the combination and optimization of wavelet components, the first principal component reflects the background reflection and does not have the reflection characteristics of the fracture-cavity reservoir, while the second principal component includes the reflection of large-scale fractures and caves. The seismic background reflection of the first principal component is removed, and the remaining seismic components are reconstructed. Compared with the original seismic profile, the background reflection of the reconstructed component is suppressed to a certain extent, which weakens the background influence and highlights the reflection characteristics of the fracture-cavity body. Figure 5 This is a seismic profile processed by robust principal component analysis according to an embodiment of the present invention, which further removes the interference of the reflection background and retains the reflection characteristics of the fracture-cavity body. The calculation results show that the fracture-cavity body features are clear and the boundaries are relatively crisp, which intuitively reflects the location of the fracture-cavity body. Figure 6This is a sweet spot property profile of an embodiment of the present invention, which can further clearly reflect the structural characteristics of the fracture-pore type reservoir. The internal information of the fracture-pore type reservoir is richer and the boundary is clearer. Figure 7 In one embodiment of the present invention, the sweet spot properties are sliced ​​along the layer, the surface property boundaries are clear, the distribution pattern of the fracture-pore type reservoir is clear, it is mainly distributed along the fault zone, and is obviously controlled by the fault, which proves that this method has good applicability to the marl fracture-pore type reservoir.

[0139] A seismic identification system for marl fracture-hole reservoirs is used for the above-mentioned seismic identification method for marl fracture-hole reservoirs, such as Figure 8 As shown, it includes: a basic unit for obtaining three-dimensional seismic data, well logging curves and logging reservoir interpretation results; a segmentation unit for determining the target layer based on the three-dimensional seismic data, well logging curves and logging reservoir interpretation results; a calculation unit for performing coherence enhancement filtering based on the three-dimensional seismic data to obtain coherence enhancement filtering data; a component unit for obtaining multiple seismic wavelet component data by wavelet decomposition based on the coherence enhancement filtering data and the target layer, removing background component data, and reconstructing all remaining seismic components to obtain seismic data of prominent fracture-pore type reservoirs; an analysis unit for highlighting the longitudinal and lateral changes of seismic amplitude and eliminating the background reflection of the surrounding rock based on the seismic data of prominent fracture-pore type reservoirs by robust principal component analysis; an extraction unit for extracting sweet spot attributes from the data after robust principal component analysis to obtain a seismic sweet spot data volume; and an identification unit for calibrating the sweet spot attribute data volume using drilling data, determining the threshold of the fracture-pore type reservoir, and obtaining the sweet spot attribute data volume reflecting the fracture-pore type reservoir.

[0140] Specifically, the segmentation unit includes: a curve processing module for generating a wave impedance curve based on the acoustic logging curve and the density logging curve; a wavelet extraction module for performing wavelet extraction on the three-dimensional seismic data to determine the seismic wavelet; a synthesis module for performing convolution processing on the wave impedance curve and the seismic wavelet to determine the well logging synthetic record; a calibration module for calibrating the synthetic record with the three-dimensional seismic data, establishing a correspondence between the seismic record and the well data, and determining the position of the top and bottom boundaries of the target layer; and a segmentation module for performing layer interpretation on the seismic data of the target layer segment to obtain the top and bottom boundaries of the target layer.

[0141] Specifically, the calculation unit includes: a first calculation module, used to select the analysis position and obtain the inclination data body and azimuth data body of the analysis position; a second calculation module, used to obtain the coherent attribute data body based on the inclination data body and the azimuth data body; a third calculation module, used to obtain the coherent enhancement filter data body based on the inclination data body, the azimuth data body and the coherent attribute data body.

[0142] Specifically, the analysis unit includes: a first analysis module, which is used to perform seismic wavelet decomposition on the seismic data at the target layer position based on the target layer position, and decompose each seismic data into multiple seismic wavelet component data of different shapes and frequencies; a second analysis module, which is used to screen out background seismic components that do not contain fracture-pore type reservoir information from the seismic wavelet components; a third analysis module, which is used to remove background components that do not contain fracture-pore type reservoir information from the seismic wavelet components, and reconstruct all remaining seismic wavelet components to obtain a seismic data body that highlights the fracture-pore type reservoir.

[0143] Specifically, the identification unit includes: a threshold determination module, which is used to calibrate the seismic sweet spot data body according to the fracture-pore type reservoir interpreted by well logging, and determine the minimum value of the sweet spot attribute corresponding to the fracture-pore type reservoir; an identification module, which is used to fill the part with the minimum value in the sweet spot data body with an invalid value, and only retain the part of the data greater than the minimum value. The part greater than the minimum value represents the distribution of the fracture-pore type reservoir, and obtains the fracture-pore type reservoir prediction data body.

[0144] It should be noted that, since the system embodiments generally correspond to the method embodiments, reference should be made to the description of the method embodiments for relevant details. The various units and modules of the seismic identification system for marl fracture-vuggy reservoirs are divided based on functional logic, but are not limited to this division; any unit that can achieve the corresponding function is sufficient. Furthermore, the specific names of the units are merely for the purpose of distinguishing them from one another and are not intended to limit the scope of protection of the present invention.

[0145] Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A seismic identification method for marl fracture-vuggy reservoirs, characterized in that: The following steps are involved: Obtain 3D seismic data, well logging curves and logging reservoir interpretation results; Determine the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results; Perform coherence enhancement filtering based on three-dimensional seismic data to obtain coherence enhancement filtering data; Based on the coherent enhancement filter data and the target layer, multiple seismic wavelet component data are obtained by wavelet decomposition, background component data are removed, and all remaining seismic components are reconstructed to obtain seismic data of the prominent fracture-hole type reservoir; Based on the seismic data of the fracture-pore type reservoir, robust principal component analysis is used to highlight the longitudinal and lateral changes of the seismic amplitude, eliminate the background reflection of the surrounding rock, and obtain the seismic reflection data of the fracture-pore type reservoir; Extracting sweet spot attributes from the data after the robust principal component analysis to obtain a seismic sweet spot data volume; The seismic sweet spot data volume is calibrated using drilling data, a threshold of a fracture-pore type reservoir is determined, and a sweet spot attribute data volume reflecting the fracture-pore type reservoir is obtained.

2. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: The three-dimensional seismic data includes seismic data collected in the field by arranging sources and detectors in an area manner, and three-dimensional pre-stack time migration seismic data obtained based on indoor seismic data processing; The logging curves include acoustic logging curves, density logging curves and natural gamma ray spectrum curves; The logging reservoir interpretation results include reservoir interfaces, thickness and physical property parameters obtained through carbonate rock optimization processing of different logging curves.

3. The seismic identification method for marl fracture-vuggy reservoirs according to claim 2, characterized in that: Determining the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results includes the following steps: Generate wave impedance curve according to acoustic logging curve and density logging curve; Extract wavelets from 3D seismic data and determine seismic wavelets; The wave impedance curve is convolved with the seismic wavelet to determine the synthetic logging record; Calibrate the synthetic records with the 3D seismic data to establish the corresponding relationship between the seismic records and the well data, and determine the position of the top and bottom boundaries of the target layer; Conduct horizon interpretation on the seismic data of the target layer to obtain the top and bottom boundaries of the target layer.

4. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: Performing coherence enhancement filtering based on 3D seismic data to obtain coherence enhancement filtering data includes the following steps: Select an analysis position and obtain the inclination data volume and azimuth data volume of the analysis position; Acquire a coherent attribute data volume based on the inclination data volume and the azimuth data volume; A coherence enhancement filtering data volume is obtained based on the inclination data volume, the azimuth data volume and the coherence attribute data volume.

5. The seismic identification method for marl fracture-vuggy reservoirs according to claim 4, characterized in that: The coherence enhancement filter equation is as follows: Where x, y, and t represent the three axes of the spatial coordinate system, t represents the time axis, and u(x, y, t) represents the three-dimensional seismic data volume. is the directional derivative vector of the three-dimensional seismic data volume, D is the diffusion tensor, ∈ is the continuous factor, and the continuous factor ∈ is as follows: Where S0 represents the initial structure tensor matrix, S ρ Represents the gradient structure tensor at the current iteration number; Tr represents the trace of the matrix, which is the sum of the diagonal elements of the matrix; the range of ∈ is [0, 1], 0≤∈≤1.

6. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: Based on the coherent enhancement filter data and the target layer, multiple seismic wavelet component data are obtained by wavelet decomposition, background component data are removed, and all remaining seismic components are reconstructed to obtain seismic data of a prominent fracture-cavity type reservoir, including the following steps: Based on the target layer position, the seismic data at the target layer position is subjected to seismic wavelet decomposition, and each seismic data is decomposed into multiple seismic wavelet component data of different shapes and frequencies; Screening out background seismic components that do not contain fracture-cavity reservoir information from seismic wavelet components; The background components that do not contain fracture-cavity reservoir information are removed from the seismic wavelet components, and all the remaining seismic wavelet components are reconstructed to obtain the seismic data volume of the prominent fracture-cavity reservoir.

7. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: Among the multiple seismic wavelet component data of different shapes and frequencies, the first component represents the seismic wavelet component with the greatest commonality in the seismic data, the second component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component, the third component is the seismic wavelet component with the greatest commonality in the seismic data after removing the first component and the second component, and so on.

8. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: Based on the seismic data of the prominent fracture-pore type reservoir, robust principal component analysis is used to highlight the vertical and horizontal changes in seismic amplitude and eliminate the background reflection of the surrounding rock to obtain seismic reflection data of the fracture-pore type reservoir, including: minutes D,E ||A|| * +γ||E|| * ; where ||A|| * represents the nuclear norm of matrix A; ||E|| * represents the l1 norm of the matrix E; A is a low-rank matrix, regarded as the background seismic reflection data; E is a sparse matrix, which is the seismic reflection data of the prominent fracture-pore type reservoir; γ is the balance parameter.

9. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: Sweet spot attributes are extracted from the data after the robust principal component analysis to obtain a seismic sweet spot data volume, including: Where A(t) is the seismic signal reflection intensity; f(t) is the instantaneous frequency; and S(t) is the seismic sweet spot attribute.

10. The seismic identification method for marl fracture-vuggy reservoirs according to claim 1, characterized in that: The method comprises the following steps: calibrating the seismic sweet spot data volume using drilling data, determining the threshold of the fracture-pore type reservoir, and obtaining the sweet spot attribute data volume reflecting the fracture-pore type reservoir. Calibrate the seismic sweet spot data volume based on the fracture-cavity reservoir interpreted by well logging, and determine the minimum value of the sweet spot attribute corresponding to the fracture-cavity reservoir; The minimum value part in the seismic sweet spot data volume is filled with invalid values, and only the part of the data greater than the minimum value is retained. The part greater than the minimum value represents the distribution of fracture-pore type reservoirs, and the fracture-pore type reservoir prediction data volume is obtained.

11. A seismic identification system for marl fracture-vuggy reservoirs, characterized in that: include: Basic unit, used to obtain 3D seismic data, well logging curves and logging reservoir interpretation results; A segmentation unit is used to determine the target layer position based on the three-dimensional seismic data, well logging curves and well logging reservoir interpretation results; A calculation unit, configured to perform coherence enhancement filtering based on the three-dimensional seismic data to obtain coherence enhancement filtering data; A component unit is used to obtain multiple seismic wavelet component data by wavelet decomposition based on the coherent enhancement filter data and the target layer, remove background component data, reconstruct all remaining seismic components, and obtain seismic data of a prominent fracture-pore type reservoir; An analysis unit is configured to use robust principal component analysis based on the seismic data of the prominent fracture-pore type reservoir to highlight the longitudinal and lateral changes of the seismic amplitude and eliminate the background reflection of the surrounding rock to obtain the seismic reflection data of the fracture-pore type reservoir; an extraction unit, configured to extract sweet spot attributes from the data after the robust principal component analysis to obtain a seismic sweet spot data volume; The identification unit is used to calibrate the seismic sweet spot data volume using drilling data, determine the threshold of the fracture-pore type reservoir, and obtain the sweet spot attribute data volume reflecting the fracture-pore type reservoir.

12. A seismic identification system for marl fracture-vuggy reservoirs according to claim 11, characterized in that: The segmentation unit comprises: A curve processing module is used to generate a wave impedance curve based on the acoustic logging curve and the density logging curve; The wavelet extraction module extracts wavelets from 3D seismic data and determines seismic wavelets; The synthesis module is used to perform convolution processing on the wave impedance curve and the seismic wavelet to determine the synthetic logging record; The calibration module is used to calibrate the synthetic records and 3D seismic data, establish the corresponding relationship between the seismic records and the well data, and determine the position of the top and bottom boundaries of the target layer; The segmentation module is used to interpret the seismic data of the target layer segment and obtain the top and bottom boundary layers of the target layer.

13. The seismic identification system for marl fracture-vuggy reservoirs according to claim 11, characterized in that: The calculation unit includes: The first calculation module is used to select an analysis position and obtain a dip data volume and an azimuth data volume of the analysis position; A second calculation module is used to obtain a coherent attribute data volume based on the inclination data volume and the azimuth data volume; The third calculation module is used to obtain a coherence enhancement filter data volume based on the inclination data volume, the azimuth data volume and the coherence attribute data volume.

14. The seismic identification system for marl fracture-vuggy reservoirs according to claim 11, characterized in that: The analysis unit comprises: The first analysis module is used to perform seismic wavelet decomposition on the seismic data at the target layer position based on the target layer position, and decompose each seismic data into a plurality of seismic wavelet component data of different shapes and frequencies; The second analysis module is used to filter out background seismic components that do not contain fracture-pore type reservoir information from the seismic wavelet components; The third analysis module is used to remove the background components that do not contain fracture-cavity reservoir information from the seismic wavelet components, reconstruct all remaining seismic wavelet components, and obtain a seismic data volume that highlights the fracture-cavity reservoir.

15. The seismic identification system for marl fracture-vuggy reservoirs according to claim 11, characterized in that: The identification unit includes: A threshold determination module is used to calibrate the seismic sweet spot data volume based on the fracture-pore reservoir interpreted by well logging, and determine the minimum value of the sweet spot attribute corresponding to the fracture-pore reservoir; The identification module is used to fill the minimum value part in the seismic sweet spot data body with invalid values, and only retain the part of the data that is greater than the minimum value. The part greater than the minimum value represents the distribution of fracture-pore type reservoirs, and obtain the fracture-pore type reservoir prediction data body.

Images