Fracture prediction method fusing prestack and poststack seismic attributes

Abstract

The application provides a fracture prediction method fusing prestack and poststack seismic attributes, which comprises the following steps: step 1, prestack fracture development direction and fracture development density prediction is performed; step 2, a fine description of buried hill reservoir fracture system is extracted and analyzed through coherent scanning based on eigenvalues of gradient structure tensor, fracture zone detection attribute and chaotic zone detection attribute; step 3, poststack fracture and hole detection and enhancement technology is used to effectively detect fracture development of the buried hill reservoir; and step 4, the minimum distribution range and the maximum distribution range of the fracture are predicted to realize buried hill reservoir fracture development area prediction. The fracture prediction method fusing prestack and poststack seismic attributes realizes buried hill reservoir fracture development area prediction, improves the reliability of fracture development prediction, accurately describes the fracture distribution characteristics, and has good application effect and popularization prospect.

Description

TECHNICAL FIELD

[0001] The present application relates to buried hill carbonate reservoir prediction research technical field, especially relates to a kind of fracture prediction method of fusion prestack and poststack seismic attribute. BACKGROUND

[0002] Fractured reservoir is one of the most important unconventional oil and gas reservoirs. Fractures can be source traps, migration pathways, and reservoirs during oil and gas accumulation. Therefore, the exploration and development of fractured reservoirs will become increasingly important, which is of great significance to the development of China's and even the world's oil industry.

[0003] Fractured carbonate reservoirs in Bohai Bay Basin generally experience multiple tectonic movements, with extremely complex internal structures, multiple types of reservoir spaces, including fractures, pores, and dissolution cavities, strong heterogeneity, deep burial, low resolution of seismic data, and great difficulty in structural interpretation. The lithology and physical properties of the reservoirs are complex and change dramatically, which adds many difficulties to the development of buried hill reservoirs.

[0004] Taking the eastern buried hill belt of Shengli offshore as an example, Chengbei 30 and Zhuanghai 10, among other complex buried hills, have been put into development, with 34 production wells and 45.14 million tons of reserves in production. Currently, only 6 oil wells can produce normally, with a recovery degree of only 5.9%. Carbonate reservoirs are mainly low-porosity and low-permeability reservoirs, with complex and diverse pore genesis and types. Porosity and permeability have a significant impact on seismic elastic modulus and velocity, which brings great difficulties to the prediction of lithology and oil and gas bearing properties of fractured and dissolution cavity reservoirs with strong heterogeneity.

[0005] Current status of fractured buried hill reservoir prediction technology:

[0006] (1) Deep burial, low main frequency of seismic data, and seismic response characterized by chaotic reflection and weak reflection, resulting in low accuracy of poststack prediction of fracture-cavity reservoirs.

[0007] (2) Prestack azimuthal anisotropy inversion fracture prediction method can better depict fracture development orientation and density, but it is difficult to effectively identify and depict dissolution cavity reservoirs with no obvious azimuthal anisotropy characteristics, and the prediction effect is not good.

[0008] (3) There are many prestack and poststack fracture prediction methods, each with its own advantages and disadvantages and limitations of application conditions. It is difficult to form a unified standard for reference in actual production on how to fully utilize the advantages of various methods and complement each other's shortcomings.

[0009] In the Chinese patent application with the application number CN201611177648.1, a method for fracture quantitative prediction based on seismic attributes is involved. The method comprises: obtaining three-dimensional seismic data and logging data of a target layer, obtaining a structural interpretation of the target layer based on the two data; performing denoising processing on the three-dimensional seismic data; performing coherent attribute calculation on the three-dimensional seismic data body after denoising processing to obtain a coherent attribute body; performing curvature attribute calculation on the three-dimensional seismic data body after denoising processing to obtain a curvature attribute body; performing value range correction on the coherent attribute body; based on the curvature attribute body and the corrected coherent attribute body, a fracture prediction attribute body is calculated; and the spatial distribution information of the underground fracture is predicted according to the fracture prediction attribute body.

[0010] In the Chinese patent application with the application number CN201510033131.4, a method for predicting a fracture density body by using a fracture density curve is involved. The method comprises the following steps: (1) preparing geological, logging and seismic data, and performing normalization processing, inversion and attribute data body optimization on the related data to obtain normalized fracture density curves of each well and an optimized data body; (2) establishing a plurality of fracture density calculation models, bringing the data on the normalized fracture density curves of each well and the curve data of the optimized data body into calculation, and selecting a fracture density calculation model whose calculation result is closest to the actual fracture density condition; and (3) bringing the optimized data body into the selected fracture density calculation model for calculation to obtain a normalized fracture density body, and then performing de-normalization processing to obtain a fracture density body in the time domain.

[0011] In the Chinese patent application with the application number CN201810965774.6, a method and device for quantitatively predicting carbonate rock fracture density are involved. The method comprises: calculating the fracture density at the well point of a corresponding drilled well by using imaging logging data; calculating the difference between the wave impedance of a formation and the wave impedance of a rock skeleton based on acoustic logging information, density logging information and a three-dimensional post-stack seismic data body; extracting a seismic attribute body from the three-dimensional post-stack seismic data body; sorting the extracted seismic attribute body by using a forward stepwise regression method based on the difference and the fracture density; dividing all drilled wells into training wells and verification wells as a training well data set and a verification well data set respectively; training a probability neural network model corresponding to different numbers of seismic attributes by using the training well data set; determining an optimal attribute set used for quantitative prediction by using the verification data set and the verification error of the verification wells; and calculating a quantitative data body of carbonate rock fracture density according to the optimal attribute set and the probability neural network model.

[0012] The above prior arts are quite different from the present application, and fail to solve the technical problems we want to solve. Therefore, we have invented a new fracture prediction method fusing pre-stack and post-stack seismic attributes. SUMMARY

[0013] The purpose of the present application is to provide a fracture prediction method fusing prestack and poststack seismic attributes, which can reasonably and accurately describe the distribution of buried hill internal reservoirs, and is of great significance for guiding the deployment of development well sites.

[0014] The purpose of the present application can be achieved by the following technical measures:

[0015] Step 1, prestack fracture development direction and fracture development density prediction is performed;

[0016] Step 2, the fracture system of the buried hill reservoir is finely characterized by extracting and analyzing the gradient structure tensor eigenvalue, the fracture zone detection attribute and the chaotic zone detection attribute through coherent scanning;

[0017] Step 3, the fracture development of the buried hill reservoir is effectively detected by using the poststack fracture and cavity detection and enhancement technology;

[0018] Step 4, the minimum distribution range and the maximum distribution range of the fracture are predicted by using two different fusion algorithms of "and" and "or", so as to realize the prediction of the fracture development area of the buried hill reservoir.

[0019] The purpose of the present application can also be achieved by the following technical measures:

[0020] In step 1, the amplitude change with azimuth is used to perform elliptical fitting by using the quantile angle gather data, and the prestack fracture development direction and fracture development density prediction is performed.

[0021] In step 1, the amplitude change with azimuth fracture detection is used to predict the development direction and density of the fracture based on the amplitude change with azimuth characteristics, the seismic wave shows the azimuth anisotropy characteristics through the fracture medium, and the azimuth anisotropy fracture detection technology is used to predict the azimuth and density of the fracture.

[0022] In step 1, when the reflected P wave passes through the fracture medium, for a fixed offset, the P wave reflection amplitude response is represented as:

[0023] R=A+Bcos2θ

[0024] Wherein, θ=φ-α is the included angle between the shot-receiver direction and the fracture trend, A is a bias factor related to the offset, B is a modulation factor related to the offset and the fracture characteristics, φ is the included angle between the fracture trend and the north direction, and α is the included angle between the shot-receiver direction and the north direction.

[0025] In the simple harmonic oscillation characteristic formula, A is regarded as the reflection intensity under the uniform medium, B is regarded as the amplitude modulation factor varying with the azimuth under the fixed offset, and B / A is a function of the fracture development density; the above relationship is approximately represented by an elliptical graph.

[0026] Let the fracture azimuth be θ, and N observation azimuths be α i , then the seismic reflection amplitude R i at azimuth α i is

[0027] R i = A + B cos 2(α i - θ) i = 1, 2, 3,..., N

[0028] For 3D wide-azimuth seismic data, divide the data into N > 3 azimuths, then the above equation becomes an over-determined equation, and the least square method is used to calculate θ, A and B.

[0029] In step 2, for 3D seismic data, the gradient structure tensor of the 3D neighborhood is expressed as follows:

[0030]

[0031]

[0032] where g is the gradient vector, u(x, y, z) is the 3D seismic data, g x , g y , g z are the gradients of the 3D seismic data in the x, y, and z directions, respectively, and T is the gradient structure tensor.

[0033] In step 2, for the eigenvalues based on the gradient structure tensor, because the tensor matrix is a positive semi-definite matrix, all its eigenvalues are greater than or equal to 0; the eigenvector corresponding to the largest eigenvalue among all the eigenvalues represents the normal direction of the local seismic horizon; sort all the eigenvalues from large to small, i.e. λ i > λ i+1 Then the eigenvalues based on the gradient structure tensor can be used to construct different attributes.

[0034] In step 2, for the fracture zone detection attribute based on the gradient structure tensor, according to the relationship between the eigenvalues and the construction, the detection attribute C fault of the fault fracture zone, referred to as the fault feature, can be constructed, which can not only identify fractures on slices, but also achieve good results in identifying fractures on profiles:

[0035]

[0036] where C fault is the detection attribute of the fault fracture zone, λ1, λ2, and λ3 are the first three largest eigenvalues, respectively.

[0037] In step 2, for the chaotic band detection attribute based on the gradient structure tensor, the chaotic attribute is based on the relative size of the eigenvalue of the local structure tensor feature, the determination of the combined parameters to reflect the boundary of the special attribute body, and the chaotic attribute can be used to detect the ordered reflection intermingled or the area without reflection;

[0038] In order to highlight the overall characteristics of the interlayer seismic reflection structure, the chaotic signal is calculated by using the eigenvector calculation method; in a given range, the gradient vector of each point is calculated, and the covariance matrix is established:

[0039]

[0040] In the formula: C ij is the cross-correlation of the i-th trace and the j-th trace;

[0041] The eigenvector corresponding to the maximum eigenvalue of the covariance matrix is calculated, that is, the gradient main direction of a certain point; if the signal-to-noise ratio of the seismic reflection wave in the stratum is high and the continuity is good, the maximum eigenvalue λ max of the gradient vector corresponding to the covariance matrix is much larger than the intermediate eigenvalue λ mid and the minimum eigenvalue λ min ; on the contrary, when the effective wave is not obvious, the maximum eigenvalue of the covariance matrix is not much different from the other two eigenvalues; therefore, the mutual relationship between the three eigenvalues can be used to distinguish the distribution rule and chaos of the amplitude value; the quantitative parameter of the chaos measurement is given:

[0042]

[0043] In the formula: J is the chaotic attribute, λ max is the maximum eigenvalue, λ mid is the intermediate value, and λ min is the minimum value;

[0044] For the structure tensor λ max , λ mid , λ min refers to the eigenvalues λ1, λ2, λ3 in the structure tensor, the larger the J value is, the more chaotic the corresponding amplitude value is, indicating that the irregular change area in the structure is larger; the smaller the J value is, the more regular the amplitude value is, indicating that the reflection rule in the structure and the stratum change is not large. Due to the mutual influence of various scattering and diffraction of stratum lithology heterogeneity, the reflection amplitude chaotic band is often shown near the boundary fault zone.

[0045] Step 3 includes:

[0046] (1) Based on the coherent body analysis algorithm, the stratum main survey line apparent dip angle and the contact line apparent dip angle scanning analysis are carried out to obtain the dip angle and azimuth information of each time point underground;

[0047] (2) Constructing the local coordinate system of each time point of the subsurface, based on the stratigraphic dip angle and azimuth angle control;

[0048] (3) Under the constraint of this local coordinate system, the spatial amplitude variation rate analysis is carried out;

[0049] (4) Filtering the obtained amplitude variation rate volume to obtain the final fracture-vug sculpture volume.

[0050] In step 4, two different fusion algorithms of "and" and "or" are used to fuse the fault and fracture-vug results detected by prestack and poststack methods to predict the minimum and maximum distribution ranges of the fractures, and to realize the prediction of the fracture development area of the buried hill reservoir.

[0051] In step 4, in the fusion algorithm, the minimum distribution range prediction is the intersection operation of multiple seismic attributes. As long as it is a possible area, whether it is a fracture or a fracture-vug, the intersection operation is carried out, and the attribute and fusion operation is carried out, which is used to predict the minimum distribution range of the fracture. The specific calculation formula is as follows:

[0052]

[0053] In the formula: A ij (k) represents the Kth seismic attribute value; i and j represent the line number and CDP number; W(k) represents the weight coefficient of the Kth attribute; B ij represents the fused attribute value.

[0054] In step 4, in the fusion algorithm, the maximum distribution range prediction is the union operation of multiple attributes. As long as it is a possible area, whether it is a fracture or a fracture-vug, the union operation is carried out, and the attribute or fusion operation is carried out, which is used to predict the maximum distribution range of the fracture. The specific calculation formula is as follows:

[0055]

[0056] In the formula: A ij (k) represents the Kth seismic attribute value; i and j represent the line number and CDP number; B ij represents the fused attribute value.

[0057] The fracture prediction method of the fusion prestack poststack seismic attribute in the application, by using the fault and fracture hole sensitive attributes detected by prestack poststack multiple methods, according to the two different fusion algorithms of 'and' and 'or', the minimum distribution range and the maximum distribution range of the fracture are predicted, the buried hill reservoir fracture development area prediction is realized, the reliability of the fracture development prediction is improved, and the fracture distribution characteristics are accurately described. The fracture prediction method of the 'and' and 'or' fusion prestack poststack seismic algorithm in the application effectively combines the prestack and poststack fracture detection methods, fully utilizes the fault and fracture hole sensitive attributes detected by prestack poststack multiple methods, and predicts the minimum distribution range and the maximum distribution range of the fracture through the two different fusion algorithms of 'and' and 'or', so that the fracture prediction research from different perspectives is realized, the reliability of the fracture development prediction is improved, the fracture distribution characteristics can be accurately described, the buried hill reservoir exploration target area is further determined, and it has important practical significance for accelerating the pace of oil and gas production capacity construction. BRIEF DESCRIPTION OF DRAWINGS

[0058] Figure 1 The flow chart of a specific embodiment of the fracture prediction method of the fusion prestack poststack seismic attribute in the application;

[0059] Figure 2 The schematic diagram of the change of the seismic reflection amplitude with the direction in a specific embodiment of the application;

[0060] Figure 3 The schematic diagram of the fracture development density and the fracture development direction of the top surface of the Paleozoic in a specific embodiment of the application;

[0061] Figure 4 The schematic diagram of the fault detection result of the top surface of the Paleozoic in a specific embodiment of the application;

[0062] Figure 5 The schematic diagram of the fracture hole detection result of the top surface of the Paleozoic in a specific embodiment of the application;

[0063] Figure 6 The schematic diagram of the fracture development prediction result of the top surface of the Paleozoic in a specific embodiment of the application. DETAILED DESCRIPTION

[0064] It should be noted that the following detailed description is exemplary and is intended to provide further explanation of the application. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the application belongs.

[0065] It is to be understood that the terms so far as the wordings are concerned in this disclosure are not limited to the specific embodiments described, but encompasses inherently known meanings to a person skilled in the art, pertaining to example embodiments in accordance with the present disclosure. Also, as used in this disclosure, the singular form "a", "an" and "the" include plural references unless the context clearly dictates otherwise. Further, it will be appreciated that the terms "comprises", "comprising", "includes", "including" and / or "contains", "containing" when used in this disclosure specify the presence of stated features, steps, operations, elements, and / or components but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or groups thereof.

[0066] As Figure 1 shown, Figure 1 is a flow chart of the fracture prediction method of the fusion of pre-stack and post-stack seismic attributes of the present disclosure. The fracture prediction method of the fusion of pre-stack and post-stack seismic attributes comprises the following steps:

[0067] Step 1, using the amplitude versus azimuth (AVA) for elliptical fitting by quantile angle gather data, pre-stack fracture development azimuth and fracture development density prediction is carried out;

[0068] Amplitude versus azimuth (AVA) fracture detection aims to predict the development azimuth and density of the fracture based on the amplitude versus azimuth characteristics, and the technical means is that the seismic wave shows azimuth anisotropy characteristics through the fracture medium, and the azimuth and density of the fracture are predicted by using the azimuth anisotropy fracture detection technology.

[0069] If the anisotropy in the rock medium is caused by a set of directional vertical fractures, according to the propagation theory of seismic wave, when P wave propagates in anisotropic medium parallel or perpendicular to the fracture direction, it has different travel speed, thereby causing the change of P wave amplitude response and the difference of travel time. The azimuth anisotropy characteristics shown by P wave through the fracture medium are the direct parameters for fracture detection. The P wave amplitude detection directional vertical fracture technology uses wide azimuth seismic observation data to study the periodic change of P wave amplitude with azimuth, and estimates the azimuth and density of the fracture.

[0070] When the reflected P wave passes through the fracture medium, for a fixed offset, the P wave reflection amplitude response can be expressed as:

[0071] R = A + Bcos2θ

[0072] Where: θ = φ - a is the included angle between the shot-receiver direction and the fracture strike, A is the offset factor related to the offset, B is the amplitude modulation factor related to the offset and the fracture characteristics, φ is the included angle between the fracture strike and the north direction, a is the included angle between the shot-receiver azimuth and the north direction.

[0073] By analogy with the formula of simple harmonic oscillation, A can be regarded as the reflection intensity under the uniform medium, B can be regarded as the amplitude modulation factor varying with the azimuth under the fixed offset, and B / A is a function of the fracture development density. The above relationship can be approximately represented by an ellipse, as shown in Figure 2 .

[0074] Let the fracture azimuth be θ, and N observation azimuths be α i (i = 1, 2,..., N) seismic trace sets, then the seismic reflection amplitude R i corresponding to the azimuth α i is

[0075] R i = A + Bcos2(α i - θ) i = 1, 2, 3,..., N

[0076] Generally, as long as there are 3 azimuths of seismic data, the equation can be solved accurately. For 3D wide-azimuth seismic data, it can be divided into multiple azimuths (N > 3) of seismic data, and the above equation becomes an over-determined equation, which can be fitted and calculated by the least square method. Therefore, the analysis method of P-wave amplitude variation with azimuth can quantitatively reveal the development azimuth and strength (development density) of the vertically aligned vertical fractures in the subsurface.

[0077] Step 2, extract and analyze the fault system of buried hill reservoir by coherent scanning, eigenvalue based on gradient structure tensor, broken zone detection attribute, and chaotic zone detection attribute;

[0078] Fracture development is closely related to fault distribution, so first of all, the fault system of the work area should be understood. Based on the lateral similarity characteristics of the post-stack 3D data body, the fault system of the buried hill reservoir is extracted and analyzed by coherent scanning, eigenvalue based on gradient structure tensor, broken zone detection attribute, and chaotic zone detection attribute.

[0079] The third generation of coherent scanning is based on the eigenvalue of the covariance matrix. Although the size of the eigenvalue can quantitatively describe the degree of change of the data body, it is not sufficient in information. Therefore, a coherent algorithm containing more information should be sought.

[0080] Coherent scanning, eigenvalue, fracture zone detection attribute, chaotic zone detection attribute structure tensor algorithm comes from image processing neighborhood. The common feature of this algorithm used in smoothing image or detecting image boundary filter is to use the characteristics of adjacent points of the center target point to separate homogeneous region. Seismic data can also be considered as gray scale image. Two-dimensional seismic data is obviously locally directional. One method to express local direction is to use angle value, which corresponds to the angle of rotation along the principal axis. Another method is to use vector to express. For three-dimensional seismic volume or three-dimensional data, the method of direction expression needs to meet several requirements, ①uniqueness, vector x and vector-x mapping value is the same value; ②average of angle extension; ③two-level separability, in the case of considering angle change, the mapping vector is the function of the initial vector. Mapping this vector with structure tensor can meet the above requirements. The mathematical property of structure tensor is that it can be extended to high bit image without manual intervention.

[0081] Structure tensor definition:

[0082]

[0083] In the above formula, - means the average value, and T represents the transpose of the vector. From the formula, it can be seen that the calculation of structure tensor is divided into two steps:

[0084] ① Use direction filter for each individual point in the image or window to calculate;

[0085] ② Output the tensor expression of the smoothed result.

[0086] Use ‖x‖ n The normalization of structure tensor makes it possible to obtain the corresponding direction by intensity contrast. Structure tensor T can be used to effectively analyze linear structure.

[0087] First, convolve the image with the first-order Gaussian difference function with a standard deviation of σ g , to obtain the gradient vector Assuming an image with a highest dimension of w, the gradient vector of the ith dimension image is:

[0088]

[0089] Then, do binary product on the gradient vector g i . Gaussian smoothing with a scale of σ T is done on the tensor, which is the structure tensor. The definition of the structure tensor of the gradient is as follows:

[0090]

[0091] The local average integral or spatial integral can be calculated by convolving the tensor element with a Gaussian kernel:

[0092]

[0093] For 3D seismic data, the structure gradient tensor of its 3D neighborhood is expressed as follows:

[0094]

[0095]

[0096] where g is the gradient vector, T is the gradient structure tensor, and gx, gy, gz are the gradients of the 3D seismic data in the x, y, and z directions, respectively.

[0097] In a plane, only one of the eigenvalues of the structure gradient tensor is nonzero, and the corresponding eigenvector is the plane orthogonal vector. When a tilt occurs in a planar stratum, the gradient energy level increases.

[0098] Based on the eigenvalues of the gradient structure tensor:

[0099] Because the tensor matrix is a positive semi-definite matrix, all its eigenvalues are greater than or equal to 0. The eigenvector corresponding to the largest eigenvalue among all the eigenvalues represents the normal direction of the local seismic horizon. All the eigenvalues are sorted in descending order, i.e., λ i >λ i+1 Then, different attributes can be constructed using the eigenvalues.

[0100] Fracture zone detection attribute based on the gradient structure tensor:

[0101] Bakker (2002) constructed a fracture zone detection attribute C fault , referred to as the fault feature, which can not only identify fractures on slices but also achieve good results in identifying fractures on profiles.

[0102]

[0103] where C fault is the fracture zone detection attribute, and λ1, λ2, and λ3 are the first three largest eigenvalues, respectively.

[0104] Chaotic zone detection attribute based on the gradient structure tensor:

[0105] The chaotic attribute reflects the boundaries of special attribute bodies based on the relative size and combination parameters of the eigenvalues of the local structure tensor. Using the chaotic attribute, we can detect areas with chaotic or no reflections between ordered reflections, which is a special lithology prediction method introduced recently.

[0106] Chaos can be used to carry out fault, dissolution zone imaging, and to classify seismic chaotic features, while chaos can reflect some special geological body boundary features, such as reflecting salt and gypsum intrusion, reef and beach structure, and river channel filling boundary features. Chaos in seismic data volume is an estimation method for measuring the lack of dip and azimuth organization structure, and is used to classify chaotic signal features.

[0107] In order to highlight the overall features of interlayer seismic reflection structure, a feature vector calculation method is used for calculation. In a given range, the gradient vector of each point is calculated, and a covariance matrix is established:

[0108]

[0109] In the formula: C ij is the cross-correlation of the i-th trace and the j-th trace.

[0110] The feature vector corresponding to the maximum eigenvalue of the covariance matrix is solved and calculated, that is, the gradient main direction of a point. If the signal-to-noise ratio of the seismic reflection wave in the stratum is high and the continuity is good, the gradient vector corresponds to the maximum eigenvalue λ max of the covariance matrix, which is much larger than the intermediate eigenvalue λ mid and the minimum eigenvalue λ min ; otherwise, when the effective wave is not obvious, the maximum eigenvalue of the covariance matrix is not much different from the other two eigenvalues. Therefore, the mutual relationship between the three eigenvalues can be used to distinguish the distribution rule and chaos of the amplitude value. Trygve Randen gives a quantitative parameter for measuring chaos:

[0111]

[0112] In the formula: J is the chaos attribute, λ max is the maximum eigenvalue, λ mid is the intermediate value, and λ min is the minimum value.

[0113] Obviously, for the structure tensor λ max , λ mid , λ min refers to the eigenvalues λ1, λ2, λ3 in the structure tensor, the larger the J value, the more chaotic the corresponding amplitude value, indicating that the irregular change area in the structure is larger; the smaller the J value, the more regular the amplitude value, indicating that the reflection in the structure is regular and the stratum change is small. Due to the mutual influence of various scattering and diffraction of stratum lithology heterogeneity, the reflection amplitude is often chaotic near the boundary fracture zone.

[0114] Step 3, using the post-stack fracture and cavity detection and enhancement technology to effectively detect the fracture development of buried hill reservoirs;

[0115] Fracture-vug enhancement technology is usually calculated on the time slice of seismic amplitude volume. It can achieve good results for relatively flat structural horizon, but when the structural horizon is complex, it will bring the influence of stratigraphic dip. In order to eliminate the influence of stratigraphic dip, the fracture-vug sculpture technology based on the local structural coordinate system constraint is developed.

[0116] The specific implementation steps are as follows:

[0117] (1) Based on the coherent body analysis algorithm, the stratigraphic main line apparent dip and the contact line apparent dip scanning analysis are carried out to obtain the dip and azimuth information of each time point underground.

[0118] (2) The local coordinate system based on stratigraphic dip and azimuth control is constructed for each time point underground.

[0119] (3) The spatial amplitude variation rate analysis is carried out under the constraint of this local coordinate system.

[0120] (4) The obtained amplitude variation rate volume is filtered to obtain the final fracture-vug sculpture volume.

[0121] Step 4, adopt two different fusion algorithms of "and" and "or" to fuse the faults and fracture-vugs detected by prestack and poststack methods to predict the minimum and maximum distribution ranges of fractures and realize the prediction of buried hill reservoir fracture development area.

[0122] The fault plane has poor lateral continuity and is mostly a fracture zone. The development of fracture-vugs is related to faults. The fusion result can reflect the development of faults and fracture-vugs. Moreover, the faults and fracture-vugs detected by prestack and poststack methods are fused together, which can accurately obtain the fracture-vug development area and its spatial distribution and realize the planar prediction of fractures.

[0123] "and" fusion algorithm

[0124] The minimum distribution range prediction is to perform intersection operation on multiple seismic attributes. No matter it is a fracture or a fracture-vug, as long as it is a possible area, the intersection is performed. Essentially, it is the fusion operation of attribute "and", which is used to predict the minimum distribution range of fractures. The specific calculation formula is as follows:

[0125]

[0126] In the formula, A ij (k) represents the Kth seismic attribute value; i, j represent line number and CDP number; W(k) represents the weight coefficient of the Kth attribute; B ij represents the fused attribute value.

[0127] "or" fusion algorithm

[0128] The maximum distribution range prediction is a union operation of multiple attributes, whether it is a fracture or a fracture-vug, as long as it is a possible area, and it is essentially an attribute "or" fusion operation, which is used to predict the maximum distribution range of the fracture.

[0129]

[0130] In the formula, A ij (k) represents the Kth seismic attribute value; i and j represent line number and CDP number; B ij represents the fused attribute value.

[0131] The fracture prediction method for fusing prestack and poststack seismic attributes provided by the application realizes the prediction of the fracture development direction and density by using the angle bin data to perform elliptical fitting through amplitude variation with azimuth (AVA); the spatial distribution of the fracture system is described by using the poststack fault scanning technology; the fracture development of the buried hill reservoir is effectively detected by using the poststack fracture-vug detection and enhancement technology; the fault and fracture-vug results detected by the prestack and poststack methods are comprehensively integrated, and the minimum distribution range and the maximum distribution range of the fracture are predicted according to two different fusion algorithms, so that the fracture development area of the buried hill reservoir is predicted, and the reliability of the fracture development prediction is improved. The method has good application effect and popularization prospect.

[0132] In a specific embodiment 1 of the application, the fracture prediction method for fusing prestack and poststack seismic attributes provided by the application includes the following steps:

[0133] Step 1, for the top surface of the Paleozoic in Chengbei 30 well area, the fracture development direction and the fracture development density (D) are predicted by using the angle bin data to perform elliptical fitting through amplitude variation with azimuth (AVA). Figure 3 The fracture development density predicted by the prestack method based on the angle bin set has high resolution, and the fracture-vug description is clear and reliable;

[0134] Step 2, the coherent scanning and the gradient structure tensor scanning are used to obtain the coherent attribute, the eigenvalue attribute, the broken zone detection attribute and the chaotic zone detection attribute (C) of the top surface of the Paleozoic. Figure 4 It can be seen that the fault plane of the top surface of the Paleozoic has poor lateral distribution continuity, and is mostly a broken zone.

[0135] Step 3, the poststack fracture-vug detection and enhancement technology is used to obtain the fracture-vug detection result (F) of the top surface of the Paleozoic. Figure 5 The fracture-vug development zone detected by the fracture-vug enhancement body is closely related to the distribution of the fault, and the fracture-vug scale is smaller than the fault scale.

[0136] Step 4, the fault and fracture results detected by prestack and poststack methods are fused by using two different fusion algorithms of "and" and "or", and the minimum distribution range and the maximum distribution range of the Paleozoic top surface fracture predicted Figure 6

[0137] Example 2:

[0138] The patent method is applied in the development of Zheng 4 buried hill, which is an old buried hill that has been developed. Through the development of fracture prediction, the fracture development advantage area is implemented, combined with the dynamic analysis of the open well, the remaining oil enrichment area is determined, and one horizontal well is deployed. After the implementation, it is self-flowing production, the initial daily oil is 32t, and the cumulative oil production is nearly 20,000 tons, which effectively promotes the secondary development of buried hill reservoir.

[0139] Example 3:

[0140] The patent is applied in the development of Shengli offshore pile oblique 473 buried hill, and the Paleozoic reservoir prediction research is carried out, and the well site is deployed. Five new wells are deployed, and the five large displacement inclined wells have been drilled, with a drilling success rate of 100%, and the reservoir thickness coincidence rate is 72%. The average single well production capacity is 55t / d in the initial stage, which meets the design requirements, and effectively supports the benefit development of Shengli buried hill reservoir.

[0141] In summary, the application researches a fracture prediction method of fusing prestack and poststack seismic attributes by "and" and "or", which realizes the prediction of prestack fracture development direction and density by using partial angle gather data through amplitude variation with azimuth (AVA) for elliptical fitting; the buried hill reservoir fracture system is finely described by extracting and analyzing through coherent scanning, eigenvalue based on gradient structure tensor, fracture zone detection attribute and chaotic zone detection attribute; the fracture development of buried hill reservoir is effectively detected by using poststack fracture detection and enhancement technology; the minimum distribution range and the maximum distribution range of the fracture are predicted according to the two different fusion algorithms of "and" and "or" by comprehensively fusing the fault and fracture results detected by prestack and poststack methods, the buried hill reservoir fracture development area is predicted, and the reliability of the fracture development prediction is improved.

[0142] Finally, it should be noted that: the above only describes the preferred embodiments of the application and is not used to limit the application, although the application has been described in detail with reference to the foregoing embodiments, and for those skilled in the art, the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the application shall be included in the protection scope of the application.

[0143] In addition to the technical features described in the specification, they are known to those skilled in the art.​

Claims

1. A method of fracture prediction by fusing prestack and poststack seismic attributes, characterized in that, The fracture prediction method of the fusion prestack poststack seismic attribute comprises: Step 1, prestack fracture development azimuth and fracture development density prediction is carried out; Step 2, through coherent scanning, the eigenvalue of gradient structure tensor, broken zone detection attribute and chaotic zone detection attribute are used to extract and analyze the fine description of buried hill reservoir fracture system; Step 3, the poststack fracture and hole detection and enhancement technology is used to effectively detect the fracture development of the buried hill reservoir; Step 4, the and or two different fusion algorithms are used to predict the minimum distribution range and the maximum distribution range of the fracture, and the fracture development area prediction of the buried hill reservoir is realized; In step 1, the amplitude change with azimuth fracture detection is based on the amplitude change with azimuth feature to predict the development azimuth and density of the fracture, the seismic wave shows the azimuth anisotropy feature through the fracture medium, and the azimuth and density of the fracture are predicted by using the azimuth anisotropy fracture detection technology; In step 1, when the reflected P wave passes through the fracture medium, for the fixed offset, the P wave reflection amplitude response is represented as: R=A+Bcos2θ Wherein θ=φ-α is the included angle between the shot direction and the fracture trend, A is the offset factor related to the offset, B is the amplitude modulation factor related to the offset and the fracture feature, φ is the included angle between the fracture trend and the north direction, and α is the azimuth of the shot direction and the north direction; In the simple harmonic oscillation characteristic formula, A is regarded as the reflection intensity under the uniform medium, B is regarded as the amplitude modulation factor which changes with the azimuth under the fixed offset, and B / A is the function of the fracture development density; the above relationship is approximately represented by an ellipse graph; Let θ be the azimuth of the fracture and let N be the number of observation azimuths α sorted in clockwise azimuth from north i The seismic trace gather is sorted into N azimuths α i The seismic reflection amplitude R i is R i = A + B cos 2(α i - θ) i = 1, 2, 3,..., N For the three-dimensional wide azimuth acquisition seismic data, the seismic data is divided into multiple azimuth angles N>3, at this time the above equation becomes an over-determined equation, and the least square method is used to fit and calculate θ and A and B values; In step 2, for the three-dimensional seismic data, the gradient structure tensor of the three-dimensional neighborhood is represented as follows: where: g is the gradient vector, u(x, y, z) is the three-dimensional seismic data, g x , g y , g z are the gradients of the three-dimensional seismic data in the x, y, z directions, respectively, and T is the gradient structure tensor; In step 2, for the fracture zone detection attribute based on the gradient structure tensor, according to the relationship between the eigenvalue and the construction, the detection attribute C describing the fault fracture zone is constructed fault In addition to being able to identify fractures on slices, the fault feature, referred to as the fault feature for short, also achieves good results in identifying fractures in profiles: In the formula: C fault λ1, λ2, λ3 are the first three eigenvalues respectively, and the detection attribute is a fault fracture zone. In step 4, the and or two different fusion algorithms are used to fuse the fault and fracture and hole results detected by the prestack poststack multiple methods, the minimum distribution range and the maximum distribution range of the fracture are predicted, and the fracture development area prediction of the buried hill reservoir is realized; In step 4, in the and fusion algorithm, the minimum distribution range prediction is the intersection operation of multiple seismic attributes, as long as the possible region is intersected, the attribute and fusion operation is carried out, and the minimum distribution range of the fracture is predicted; the specific calculation formula is as follows: In the formula: A ij (k) represents the Kth seismic attribute value; i, j represent line number, CDP number; W(k) represents the weight coefficient of the Kth attribute; B ij represents the fused attribute value; In step 4, in the or fusion algorithm, the maximum distribution range prediction is the union operation of multiple attributes, as long as the possible region is unioned, the attribute or fusion operation is carried out, and the maximum distribution range of the fracture is predicted; the specific calculation formula is as follows: wherein: A ij (k) denotes the Kth seismic attribute value; i, j denote line number, CDP number; B ij denotes the fused attribute value.

2. The method of claim 1, wherein the method further comprises: In step 1, the azimuthal bin data is used to carry out the ellipse fitting through the amplitude change with azimuth, and the prestack fracture development azimuth and fracture development density prediction is carried out.

3. The method of claim 1, wherein the method further comprises: In step 2, for the eigenvalues based on the gradient structure tensor, because the tensor matrix is a positive semi-definite matrix, all eigenvalues thereof are greater than or equal to 0; the eigenvector corresponding to the largest eigenvalue among all the eigenvalues represents the normal direction of the local seismic horizon; all the eigenvalues are sorted in descending order, i.e. λ i > λ i+1 Then, the eigenvalues based on the gradient structure tensor can be used to construct different attributes.

4. The method of claim 1, wherein the method further comprises: In step 2, for the chaotic zone detection attribute based on the gradient structure tensor, the chaotic attribute is based on the relative size of the local structure tensor eigenvalue, the determination of the combination parameter to reflect the boundary of the special attribute body, and the ordered reflection intermingled or non-reflection region can be detected by using the chaotic attribute; In order to highlight the overall characteristics of the interlayer seismic reflection structure, the eigenvector calculation method is used to calculate the chaotic signal. In a given range, the gradient vector of each point is calculated, and the covariance matrix is established: wherein: C ij is the cross-correlation of the ith track with the jth track; Solving the calculation of the covariance matrix corresponding to the largest eigenvalue of the eigenvector, that is, the gradient of a point; if the stratum is high in signal-to-noise ratio and good in continuity, the gradient vector corresponds to the largest eigenvalue λ max of the covariance matrix mid , the middle eigenvalue λ min is much larger; on the contrary, when the effective wave is not obvious, the largest eigenvalue of the covariance matrix is not much different from the other two eigenvalues; therefore, the mutual relationship between the three eigenvalues can be used to judge the distribution rule and chaos of the amplitude value; the quantitative parameter of chaos measurement is given: wherein: J is a chaotic attribute, λ max is a maximum eigenvalue, λ mid is an intermediate value, λ min is a minimum value; For structural tensor λ max ,λ mid ,λ min refers to the eigenvalue λ1, λ2, λ3 in the structural tensor, the greater the value of J, the more chaotic the corresponding amplitude value, indicating that the larger the irregular change area in the structure; the smaller the value of J, the more regular the amplitude value, indicating that the reflection is regular and the stratum changes little; due to the mutual influence of various scattering and diffraction of stratum lithology heterogeneity, it is often manifested as a chaotic reflection amplitude zone near the boundary fault zone.

Images