CN112698397B - Method for describing reservoir stratum with sliding fracture cavity in basin area - Google Patents

Abstract

The invention discloses a method for describing a reservoir with a sliding fracture cavity in a basin area. Aiming at obvious control characteristics of the carbonate reservoir in the deep buried region, the invention creatively provides two parts of treatment of decomposing a fracture cavity into a boundary line and an internal structure, finely inscribing the boundary line outline outside the fracture zone cavity controlled by sliding fracture, and providing targeted different identification inscription methods for the reservoir with the internal hole structure and the reservoir with the fracture structure, thereby finally realizing quantitative inscription of the beneficial regions of the fracture zone cavity and the fracture reservoir with the internal hole. In addition, logging data is not necessary in the calculation flow of the invention, only the existing drilling is used for threshold value correction and fitness verification when the prediction is finished, and the invention provides constructive opinion for the exploration deployment of the less wells and the no well areas.

Description

Method for describing reservoir stratum with sliding fracture cavity in basin area

Technical Field

The invention relates to the technical field of petroleum exploration, in particular to a method for describing a reservoir with a sliding fracture cavity in a basin area.

Background

The carbonate reservoir prediction technology mainly adopts geostatistical inversion at present, takes constrained sparse pulse inversion results as input data, takes logging data as main, uses geostatistical rules and seismic data constraint for interwell changes, generates attribute simulation results of a plurality of equiprobable time domains, preferably takes the most coincident results with logging and geological analysis as research objects, and the prediction resolution can reach 2-6m.

The main disadvantage of geostatistical inversion is that the method is mainly aimed at obvious beaded reflective reservoirs formed by the dissolution of atmospheric fresh water after the formation lifting stratum is degraded, the beaded reservoirs can be obviously identified on the seismic section, and the method has stronger energy differentiation, but has poor reflection on the reservoir with sliding fracture control and the micro beaded with the aperture less than 10 m. One of the most prominent types of carbonate reservoirs is the cavity formed by the control of the slip fracture of the trunk, no obvious beading reflection, the reservoir seismic response is a combination of a series of strong peaks or valleys distributed vertically along the slip fracture, and the exploration area lacks a large amount of well logging data, mostly few or no well data, and geostatistical inversion has obvious disadvantages for this situation.

Disclosure of Invention

The invention aims to provide a method for describing a sliding fracture cavity reservoir in a basin area, which aims to solve the problems of identification and description of the fracture cavity reservoir controlled by sliding fracture in the current production, has high accuracy, accords with geological recognition, and provides constructive opinion for exploration and deployment in a few-well and no-well area.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

The invention provides a method for describing a reservoir with a sliding fracture cavity in a basin area, which comprises the following steps:

s100, acquiring drilling data, logging data and seismic data of a basin area sliding fracture cavity reservoir;

s200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing bead space distribution by adopting a wavelet decomposition-based frequency-reduction speed inversion method;

s300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify a fracture-cavity abnormal body;

s400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir;

s500, calculating the discontinuous attribute of the three-dimensional earthquake to describe a crack structure in a fracture cavity reservoir;

s600, carrying out geostatistical trend inversion on the basis of wave impedance inversion of the hole reservoir, so as to improve the identification and characterization capability of the small fracture-hole reservoir, and combining with S500, completing the characterization of the inside of the fracture cavity reservoir.

Each step is described in detail below.

S200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing the space distribution of the beads by adopting a wavelet decomposition-based frequency-reduction speed inversion method.

S300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify the abnormal body of the fracture-cavity.

The configuration attribute combines a new attribute analysis method-gradient configuration analysis method introduced into seismic interpretation from the image processing field in recent years, a structural gradient vector is constructed through a gradient vector of post-stack seismic amplitude data, the characteristic value of the vector is decomposed to obtain a characteristic value and a characteristic vector, and the sedimentary landform units of different texture features are characterized according to the image structural meaning of the characteristic value.

The configuration attribute utilizes the seismic tensor field to judge the trend of the seismic event, and the core is to judge the occurrence of the seismic event by referring to the method for judging the fluid potential energy to achieve the aim of identifying the abnormal body of the fracture and tunnel.

Preferably, the process of calculating the gradient configuration attribute comprises the following steps of:

s301, establishing gradient configuration attributes:

calculating a gradient vector (directional derivative) for each point of the three-dimensional seismic data volume:

wherein x, y and z are the line, channel position and double-journey travel time direction of the three-dimensional earthquake voxels respectively; g is the amplitude gradient field; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; u is the derivative matrix in the corresponding direction;

The gradient is smoothed using a three-dimensional gaussian window:

wherein G is a Gaussian kernel function; g is the amplitude gradient field; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; sigma (sigma) _g Is the length corresponding to the i direction;

s302, calculating gradient configuration attribute T:

s303, smoothing each component of the gradient configuration attribute by using a Gaussian window:

the meaning of the average gradient configuration attribute is the following two aspects: (1) characterizing texture features within a region with average gradient topography properties; (2) for suppressing abrupt changes in the profile properties due to noise.

S304, calculating characteristic values and characteristic vectors of gradient configuration attributes:

|Tv-λv|＝0

wherein λ and ν are a eigenvalue and eigenvector, respectively; because the matrix T is a real symmetric matrix, lambda is satisfied ₁ ≥λ ₂ ≥λ ₃ >0. From the nature of the real symmetric matrix, three eigenvectors v ₁ 、ν ₂ 、ν ₃ Two pairs of the two-way cross-linking device are orthogonal, 1, 2 and 3 refer to x, y and z directions respectively.

1)λ ₁ ＝λ ₂ ＝λ ₃ ＞0

In this case, the system is isotropic, and since the seismic data is layered, it is rare that even in the case of very good seismic data it is rare that it is present.

2)λ ₁ ＞λ ₂ ＝λ ₃ ＝0

One feature value is large and the other 2 are very small close to 0, which is usually the case in the case of good seismic data quality and relatively flat seismic event, where lambda ₁ The corresponding vector direction approaches the Z direction.

3)λ ₁ ＞λ ₂ ＞λ ₃ ＝0

This situation occurs when the quality of the seismic data is poor, and the filtering difficulty is high, so that the filtering vector needs to be accurately calculated according to the eigenvalue and the eigenvector, and meanwhile, the fracture and the like need to be effectively stored.

The tensor attribute is well reflected on the tensor attribute through calculating the change gradient of the data in three directions of space, namely the bead string or the fracture, namely the severe change of amplitude or the phase axis dislocation, but the transverse resolution of the tensor body is too low, and the tensor attribute is well carved only as a fracture zone shape. The feature value 1 and the feature value 2 are good in appearance on the broken belt profile according to the calculated tensor body, and the signal to noise ratio of the feature value 3 is too low to be recommended.

S400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir.

Based on the method for describing the reservoir stratum of the sliding fracture cavity of the basin area, preferably, the transverse fracture range of fracture is calibrated by utilizing the drilling time curve of the sidetracking horizontal well in S400 so as to determine the threshold value range of the configuration attribute;

the drilling curve changes to show a linear change characteristic, the broken erosion zone shows a low drilling condition and is accompanied by a blow-out and leakage phenomenon, the position is used for calibrating a reservoir threshold value, namely the reservoir threshold value is used as a reservoir boundary point, the reservoir is arranged inside, and the matrix is arranged outside.

The drilling curve of the bedrock shows a smooth change rule, when the drilling curve enters the range of the broken zone (namely, when the drilling curve is at the boundary of a broken cavity reservoir), the drilling curve changes to show a linear change characteristic, the broken corrosion zone shows a low drilling state and is accompanied by a blow-out loss phenomenon, the position is used for calibrating a reservoir threshold value, the reservoir is arranged inside, and the matrix is arranged outside.

S500, calculating the discontinuous attribute of the three-dimensional earthquake to describe a crack structure in a fracture cavity reservoir;

the present invention uses the discontinuity properties as a fracture structure inside the fracture cavity reservoir. The discontinuity property analysis technology is developed from the similarity of coherence based on cross-correlation calculation of the first generation and inclination superposition calculation based on similarity of the second generation to the structure analysis technology based on the discontinuous detection of the intrinsic structure of the data covariance matrix of the third generation. The data of one analysis channel consisting of N sample point elements in a time window and J channels in the window taking the channel as a unit are arranged together to form a data matrix:

wherein u is _ij Representing the amplitude value of the ith seismic wave on the jth trace; n rows of matrix U areAnalyzing the nth amplitude value on each seismic trace in the body, and then the variance matrix of the nth sample is:

So that covariance matrix formed by covariance of each element in the time window of coherent analysis is:

let the covariance matrix M have J eigenvalues lambda ordered from big to small ₁ ，λ ₂ ，…，λ _J The corresponding feature vector is v ₁ ，ν ₂ ，…，ν _J . Then feature-based coherent volume computationThe formula is:

wherein lambda is _j And calculating eigenvalues of the covariance matrix of the seismic data in the time window. In theory, the third generation algorithm C3 based on the eigenvalue is superior to the second generation algorithm C2 based on the multi-channel cross-correlation and the C1 based on the two-channel cross-correlation, mainly because the noise which is higher than the additional noise level is automatically removed when the signal is higher than the additional noise level, and the transverse resolution is greatly improved. But the C3 calculation amount is large.

The fine coherence technique utilizes multi-channel similarity to convert a three-dimensional data body into a correlation coefficient data body through calculation, and compares the similarity of local seismic waveforms (the point with a lower coherence value is better related to the discontinuity of the reflected wave waveform) through fine coherence processing of the three-dimensional data body, so that various geological phenomena can be intuitively distinguished.

S600, inversion of geology statistics trend is carried out so as to improve the identification and characterization capability of the small fracture-cavity reservoir, and the characterization of the inside of the fracture cavity reservoir is completed by combining S500.

In order to further improve the identification and characterization capability of the seismic inversion result on the small fracture-cavity reservoir, the invention adopts geology statistical trend inversion. The trend analysis method is a method for simulating (or fitting) the spatial distribution of the geographic data and the regional variation trend thereof by using a mathematical method to decompose trend values and residual values in the actual data of the spatial geographic element distribution to construct a mathematical model. The specific principle of the method can be divided into two steps:

s601, establishing a trend model

Let the actual observation data of a certain space geographic element be Z _i (x _i |y _i ) (i=1, 2, …, n) with a trend fit of Z _i (x _i |y _i ) Then there is Z _i (x _i |y _i )＝Z^ _i (x _i |y _i )+∈ _i Wherein: epsilon _i I.e. the residual value (residual value). When (x _i |y _i ) When spatially varying, the above-mentioned method is just carvedAnd drawing the interactive relation among the actual distribution curved surface, the trend surface and the residual surface of the geographic elements.

The trend surface (or body) is calculated from the actual observation value, and a regression analysis method is generally adopted to minimize the sum of squares of residual errors, namely

According to the observed value x _i ,y _i ,z _i (i=1, 2, …, n) determining the coefficients a of the polynomial ₀ ，a ₁ ，…，a _p The sum of squares of the residuals is minimized.

S602, moderately testing the trend surface model, wherein the moderately testing comprises the degree of fit coefficient R2 testing and F testing of the trend surface and the actual surface.

The fitting degree coefficient R2 of the trend surface and the actual surface is an important index for measuring the fitting goodness of the regression model. The goodness of fit of the regression model is generally expressed as the sum of squares of the total dispersion of the variable z and the proportion of the sum of squares of the regression. The larger R2, the higher the fitting degree of the trend surface.

F, testing the significance of the whole trend surface regression model. And determining whether the regression relationship between the variable z and the independent variables x and y is obvious or not by using the ratio of the regression square sum to the residual square sum in the total dispersion square sum of the variable z.

Solving the difference between the regression square sum of the higher order polynomial equation and the regression square sum of the lower order polynomial equation; dividing the difference by the difference in degrees of freedom of the sum of squares of the regression to obtain a regression mean square error due to the increase of the polynomial degree; dividing the mean square error by the residual mean square error of the higher order polynomial to obtain a moderately comparative test value F of the trend surface models of the two successive orders.

The trend analysis technology is introduced into the formation uniformity research and mainly based on the heterogeneity of carbonate rock, the background value of the whole impedance is given, the opposite low-resistance area is searched in the high-low main impedance area to indicate the favorable development area of the small fracture hole, the interference of the background impedance value on the information of the small fracture hole is eliminated, and the favorable area can be indicated more accurately.

Based on the basin area sliding fracture cavity reservoir layer describing method, preferably, the trend body generating method in the geostatistical trend inversion comprises the following steps:

And carrying out large grid filtering on the absolute impedance attribute body to form a gradual change trend impedance trend body, subtracting the original absolute impedance attribute body from the gradual change trend body to obtain a residual impedance body, and representing the residual group of antibodies in the form of conventional seismic attributes.

According to the principle of the trend surface method, the absolute impedance is subtracted from the trend body to obtain a relatively low geology statistics trend wave impedance inversion body.

Based on the method for describing the bench basin area sliding fracture cavity reservoir, preferably, the method for describing the bench basin area sliding fracture cavity reservoir comprises the following steps of:

s100, acquiring drilling data, logging data and seismic data of a basin area sliding fracture cavity reservoir;

s200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing bead space distribution by adopting a wavelet decomposition-based frequency-reduction speed inversion method; the method specifically comprises the following steps:

s201, establishing an initial wave impedance low-frequency model

Establishing a space wave impedance model by combining logging data and seismic data; obtaining a layer speed by utilizing the superposition speed in seismic processing and through speed conversion, and then converting the density by utilizing a Gardner formula to finally obtain a wave impedance model;

S202, inversion processing

The method comprises the steps of restraining sparse pulse inversion, wherein each channel of initial wave impedance calculated according to an objective function is adjusted, and the initial wave impedance comprises the adjustment of reflection coefficients, so that a full-area relative wave impedance data body is obtained;

the high frequency band comes from a constraint sparse pulse inversion algorithm, and the low frequency band is combined by a wave impedance model established in S201;

s300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify a fracture-cavity abnormal body;

wherein the process of calculating gradient profile attributes comprises:

s301, establishing gradient configuration attributes:

calculating gradient vectors of each point of the three-dimensional seismic data volume:

wherein x, y and z are the line, channel position and double-journey travel time direction of the three-dimensional earthquake voxels respectively; g is the amplitude gradient field; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; u is the derivative matrix in the corresponding direction;

the gradient is smoothed using a three-dimensional gaussian window:

wherein G is a Gaussian kernel function; g is the amplitude gradient field; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; 1. 2 and 3 respectively refer to x, y and z directions; sigma (sigma) _g Is the length corresponding to the i direction;

s302, calculating gradient configuration attribute T:

s303, smoothing each component of the gradient configuration attribute by using a Gaussian window:

S304, calculating characteristic values and characteristic vectors of gradient configuration attributes:

|Tv-λv|＝0

wherein λ and ν are a eigenvalue and eigenvector, respectively; lambda (lambda) ₁ ≥λ ₂ ≥λ ₃ >0, three eigenvectors v ₁ 、ν ₂ 、ν ₃ Two pairs of the two orthogonal directions, 1, 2 and 3 respectively refer to x, y and z directions; the method comprises the following steps:

1)λ ₁ ＝λ ₂ ＝λ ₃ > 0, in this case isotropic;

2)λ ₁ ＞λ ₂ ＝λ ₃ =0, at this time λ ₁ The corresponding vector direction approaches the Z direction;

3)λ ₁ ＞λ ₂ ＞λ ₃ =0, in which case the filter vector is accurately calculated from the eigenvalues and eigenvectors, while preserving valid breaks;

s400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir;

the method comprises the steps of calibrating a fracture transverse crushing range by using a drilling curve of a sidetracking horizontal well so as to determine a configuration attribute threshold value range; the process for determining the configuration attribute threshold value range comprises the following steps:

the drilling curve changes and shows a linear change characteristic, the broken corrosion zone shows a low drilling state and is accompanied by a blow-out leakage phenomenon, the position is used for calibrating a reservoir threshold value, namely the reservoir threshold value is used as a reservoir boundary point, the reservoir is arranged inside, and the matrix is arranged outside;

s500, calculating the discontinuous attribute of the three-dimensional earthquake to describe a crack structure in a fracture cavity reservoir;

Wherein the step of calculating the seismic discontinuity attributes comprises:

the data of one analysis channel consisting of N sample point elements in a time window and J channels in the window taking the channel as a unit are arranged together to form a data matrix:

wherein u is _ij Representing the amplitude value of the ith seismic wave on the jth trace; n rows of matrix U areAnalyzing the nth amplitude value on each seismic trace in the body, and then the variance matrix of the nth sample is:

so that covariance matrix formed by covariance of each element in the time window of coherent analysis is:

let the covariance matrix M have J eigenvalues lambda ordered from big to small ₁ ，λ ₂ ，…，λ _J The corresponding feature vector is v ₁ ，ν ₂ ，…，ν _J The coherent volume calculation formula based on the feature structure is:

wherein lambda is _j Calculating eigenvalues of a covariance matrix of the seismic data in the time window;

s600, carrying out geostatistical trend inversion on the basis of wave impedance inversion of the hole reservoir so as to improve the identification and characterization capability of the small fracture hole reservoir, and completing the characterization of the inside of the fracture cavity reservoir by combining with S500;

wherein, the process of geostatistical trend inversion comprises:

s601, establishing a trend model;

let the actual observation data of a certain spatial geographic element be: z is Z _i (x _i |y _i ) (i=1, 2, …, n) with a trend fit of Z _i (x _i |y _i ) Then there is Z _i (x _i |y _i )＝Z^ _i (x _i |y _i )+∈ _i Wherein: epsilon _i Namely, a residual value, namely, a residual value;

calculating a trend surface or a trend body from an actual observed value, and adopting a regression analysis method to enable the square sum of residual errors to be minimum, namely:

according to the observed value x _i ,y _i ,z _i (i=1, 2, …, n) determining the coefficients a of the polynomial ₀ ，a ₁ ，…，a _p Minimizing the sum of squares of the residuals;

s602, moderately testing the trend surface model, wherein the moderately testing comprises the degree of fit coefficient R2 testing and F testing of the trend surface and the actual surface.

Aiming at the situation that the fracture control characteristics of the carbonate reservoir in the deep buried region are obvious, the invention creatively provides two parts of treatment for decomposing the fracture cavity into a boundary line and an internal structure, finely describes the boundary line outline outside the fracture zone cavity controlled by sliding fracture, provides different targeted identification and description methods for the reservoir with the internal pore structure and the reservoir with the fracture structure, and finally realizes quantitative description of the fracture zone cavity and the favorable region of the fracture reservoir with the internal pore.

According to the invention, the geological concept is fully integrated into the earthquake reservoir prediction technology, the boundary and the internal structure of the fracture cavity are characterized by using different earthquake attributes, the coincidence degree with the original earthquake data is extremely high, logging data in a calculation flow is unnecessary, and only the existing drilling is used for threshold value correction and coincidence degree verification when the prediction is finished.

Drawings

FIG. 1 is a flow chart of the inversion of the down-converted wave impedance based on wavelet decomposition.

FIG. 2a is a plan view of a fracture cavity reservoir boundary determined from the northbound 5 well configuration property in an example.

FIG. 2b is a cross-sectional view of an example northbound 5-well region configuration property determination fracture cavity reservoir boundary property.

FIG. 2c is a graph of variation in drilling of different configurations of a reservoir with a fracture during drilling of a northbound 5-well in an example.

FIG. 2d is a graph of different reservoir types for the fracture cavity according to the profile properties for the north-south 5 well region in the example.

FIG. 3a is a seismic section of a full depth 1 well in an example application.

Figure 3b is a profile cross-section of a full depth 1 well in an example application.

FIG. 3c is an external profile of a full deep block reservoir determined using configuration property thresholds in an embodiment.

Fig. 4 is a three-dimensional representation of a full depth block configuration attribute fracture cavity reservoir quantitative characterization in an application example.

FIG. 5a is a three-dimensional representation of a seismic section of a full depth 1 well in an application example.

Fig. 5b is a three-dimensional representation of a full depth 1 well fracture cavity reservoir boundary line in an application example.

FIG. 5c is a second three-dimensional representation of a reservoir boundary line of a full depth 1 well fracture cavity in an application example.

Fig. 6a is a full depth seismic section of an example application.

Fig. 6b is a cross-sectional view of a full depth discontinuity property in an application example.

FIG. 6c is a cross-sectional view of a full depth discontinuity attribute and seismic overlap in an embodiment.

FIG. 7a is a full depth segmented seismic section view of an application example.

FIG. 7b is a conventional wave impedance inversion profile of a full depth segment in an application example.

Fig. 7c is a cross-sectional view of the residual wave impedance of a full depth segment in an application example.

Fig. 8a is a three-dimensional representation of the engraving (topographical properties) of the broken cavity boundary line of the obliquely-split overlapping extrusion section 1 in the application example.

Fig. 8b is a three-dimensional representation of the internal structure (inversion of geostatistical trend + discontinuity) of a fracture cavity with dual parameters characterizing the fracture cavity of the oblique fracture overlapping extrusion section 1 in an application example.

Fig. 9a is a three-dimensional representation of the engraving (form attribute) of the broken cavity boundary line of the pigtail section 2 in the application example.

Fig. 9b is a three-dimensional representation of the dual-parameter characterization of the fracture cavity internal structure (residual wave impedance + discontinuity) of the fracture cavity of the pigtail segment 2 in the application example.

Detailed Description

In order to more clearly illustrate the present invention, the present invention will be further described with reference to preferred embodiments. It is to be understood by persons skilled in the art that the following detailed description is illustrative and not restrictive, and that this invention is not limited to the details given herein.

The embodiment of the invention provides an optimization scheme, namely a method for describing a reservoir with a sliding fracture cavity in a basin area, which comprises the following steps of:

s100, drilling data, logging data and seismic data of the basin area sliding fracture cavity reservoir are obtained.

S200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing string bead space distribution by adopting a wavelet decomposition-based frequency-reduction speed inversion method:

because of the well-free inversion, the space distribution of the beads is characterized by adopting a down-conversion speed inversion method based on wavelet decomposition. The main flow is shown in figure 1.

(1) Establishment of initial wave impedance low-frequency model

Establishing a wave impedance model as close as possible to the actual formation conditions is a fundamental way to reduce the polynomials of the inversion end result. The logging data reveals details of wave impedance change of the rock stratum in longitudinal direction, the seismic data continuously records depth change of a wave impedance interface, and the combination of the two provides necessary conditions for accurately establishing a space wave impedance model. The process of establishing the wave impedance model is actually a process of correctly combining the information of the earthquake interface and the wave impedance of the well logging, namely correctly explaining the wave impedance interface with control function for the earthquake, and endowing proper wave impedance information for the stratum among the wave impedance interfaces for the well logging.

The method comprises the steps of obtaining a layer speed by utilizing the superposition speed in seismic processing and through speed conversion, obtaining a layer speed by utilizing the Gardner formula to convert the density, and finally obtaining the wave impedance model.

(2) Inversion processing

And (3) carrying out constraint sparse pulse inversion on each path of calculated initial wave impedance according to an objective function, wherein the adjustment comprises the adjustment of reflection coefficients, and obtaining a full-area relative wave impedance data body. The wave impedance data is not full band, it lacks the low and high frequency information that seismic data lacks. The high frequency band comes from the constraint sparse pulse inversion algorithm, and the low frequency band is combined by the wave impedance model established in the front to obtain the wave impedance antibody with more abundant low frequency information.

In the inversion process, the same method is applied, sparse pulse wave impedance inversion processing is also carried out on the basis of original seismic data, the longitudinal and transverse resolution of inversion results are very low, and the inversion results are distributed in a sheet-shaped continuous mode and do not accord with the geological rule of a carbonate fracture-cavity reservoir in the area, so that the wave impedance inversion based on wavelet decomposition is proved to be an effective inversion method suitable for the area.

S300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify the abnormal body of the fracture-cavity.

The configuration attribute combines a new attribute analysis method-gradient configuration analysis method introduced into seismic interpretation from the image processing field in recent years, a structural gradient vector is constructed through a gradient vector of post-stack seismic amplitude data, the characteristic value of the vector is decomposed to obtain a characteristic value and a characteristic vector, and the sedimentary landform units of different texture features are characterized according to the image structural meaning of the characteristic value.

The configuration attribute utilizes the seismic tensor field to judge the trend of the seismic event, and the core is to judge the occurrence of the seismic event by referring to the method for judging the fluid potential energy to achieve the aim of identifying the abnormal body of the fracture and tunnel.

The specific calculation process is as follows:

s301, establishment of gradient configuration attribute

Calculating a gradient vector (directional derivative) for each point of the three-dimensional seismic data volume:

wherein x, y and z are the line, channel position and double-journey travel time direction of the three-dimensional earthquake voxels respectively; g is the amplitude gradient field; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; u is the derivative matrix in the corresponding direction;

the gradient is smoothed using a three-dimensional gaussian window:

Wherein G is a Gaussian kernel function; g is the amplitude gradient field; g1, g2, g3 are directional derivatives in the x, y, z directions, respectively; 1. 2 and 3 respectively refer to x, y and z directions; sigma (sigma) _g Is the length corresponding to the i direction.

S302, calculating gradient configuration attribute T:

s303, smoothing each component of the gradient configuration attribute by using a Gaussian window:

the meaning of the average gradient configuration attribute is the following two aspects: (1) characterizing texture features within a region with average gradient topography properties; (2) for suppressing abrupt changes in the profile properties due to noise.

S304, calculating characteristic values and characteristic vectors of gradient configuration attributes:

|Tv-λv|＝0

where λ and ν are eigenvalues and eigenvectors, respectively. Because the matrix T is a real symmetric matrix, lambda is satisfied ₁ ≥λ ₂ ≥λ ₃ >0. From the nature of the real symmetric matrix, three eigenvectors v ₁ ,ν ₂ ,ν ₃ Two by two.

1)λ ₁ ＝λ ₂ ＝λ ₃ ＞0

In this case, the system is isotropic, and since the seismic data is layered, it is rare that even in the case of very good seismic data it is rare that it is present.

2)λ ₁ ＞λ ₂ ＝λ ₃ ＝0

One feature value is large and the other 2 are very small close to 0, which is usually the case in the case of good seismic data quality and relatively flat seismic event, where lambda ₁ The corresponding vector direction approaches the Z direction.

3)λ ₁ ＞λ ₂ ＞λ ₃ ＝0

This situation occurs when the quality of the seismic data is poor, and the filtering difficulty is high, so that the filtering vector needs to be accurately calculated according to the eigenvalue and the eigenvector, and meanwhile, the fracture and the like need to be effectively stored.

The tensor attribute is well reflected on the tensor attribute through calculating the change gradient of the data in three directions of space, whether the data is beaded or broken, whether the amplitude is changed severely or the phase axis is misplaced, but the transverse resolution of the tensor body is too low, and the tensor attribute is well carved only as a broken belt shape. From the calculated tensor, the eigenvalue 1 (λ ₁ ) And a characteristic value of 2 (lambda) ₂ ) The profile of the fracture zone is better represented, and the characteristic value is 3 (lambda ₃ ) The signal to noise ratio is too low to be recommended.

S400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir.

The transverse breaking range of the fracture is marked by utilizing the curve of the sidetracking horizontal well during drilling, most sidetracking horizontal wells which are just broken are empty and lost, no logging curve is generally available, and the width of the broken core and the breaking zone range can be judged by utilizing the change rule of the curve during drilling. The drilling curve of the bedrock generally shows a stable change rule, when the drilling curve enters the range of the broken belt, the drilling curve changes, according to different wells, the drilling curve shows different change trends, but the whole of the drilling curve shows linear change characteristics, and the broken corrosion belt shows low drilling and is accompanied by the phenomenon of emptying leakage.

The example shown in fig. 2a depicts the planar boundary of a fracture-cavity reservoir using the profile attribute for the northbound 5 well region, which is the top surface of an otto-room set, and the profile attribute threshold used is the as-drilled curve for the example (as shown in fig. 2 c). The example in fig. 2b is the use of topographical properties to delineate the cross-sectional boundaries of a fractured-vuggy reservoir, which is shown to include an atrial set of the otto system and a hawk-mountain set of the barrage set. In the application example of the invention, the drilling leakage point is used for calibrating the threshold value of the gradient configuration attribute, and the plane accords with the overall geological knowledge. The transverse breaking range of the fracture is marked by utilizing the curve of the sidetracking horizontal well during drilling, most sidetracking horizontal wells which are just broken are empty and lost, no logging curve is generally available, and the width of the broken core and the breaking zone range can be judged by utilizing the change rule of the curve during drilling. The drilling curve of the bedrock generally shows a stable change rule, when the drilling curve enters the range of the broken belt, the drilling curve changes, according to different wells, the drilling curve shows different change trends, but the whole of the drilling curve shows linear change characteristics, and the broken corrosion belt shows low drilling and is accompanied by the phenomenon of emptying leakage. Calibrating the drilling time and tensor data body, wherein different attribute values represent the fracture transverse difference, combining experience of drilling time and threshold value determination in the North China 1 three-dimensional drilling time, and drilling a well in the North China 5 well region because of SHB5 and SHB5-2 directly at the core part of the fracture zone can not calibrate the threshold value, carrying out threshold value calibration by adopting the SHB5-4H well, and finally determining that the tensor attribute threshold value is 47. FIG. 2d shows a method for distinguishing different reservoir types of a fracture cavity according to the configuration attribute, wherein the reservoir types of the fracture cavity are divided into three types according to the size of the configuration attribute and the geological knowledge, the configuration attribute is a type of reservoir with the configuration attribute being >20, and the reservoir types belong to bead reflection with stronger energy in seismic data; 20> the structural attribute >10 is a type II reservoir, which is a slightly weak energy lamellar reflection in seismic data; three types of reservoirs are the 10> topography attribute, clutter in the seismic data.

S500, calculating the discontinuous attribute of the three-dimensional earthquake to describe the crack structure in the fracture cavity reservoir.

The present invention uses the discontinuity properties as a fracture structure inside the fracture cavity reservoir.

The discontinuity property analysis technology is developed from the similarity of coherence based on cross-correlation calculation of the first generation and inclination superposition calculation based on similarity of the second generation to the structure analysis technology based on the discontinuous detection of the intrinsic structure of the data covariance matrix of the third generation. The data of one analysis channel consisting of N sample point elements in a time window and J channels in the window taking the channel as a unit are arranged together to form a data matrix:

wherein u is _ij Representing the amplitude value of the ith seismic wave on the jth trace. N rows of matrix U areAnalyzing the nth amplitude value on each seismic trace in the body, and then the variance matrix of the nth sample is:

so that covariance matrix formed by covariance of each element in the time window of coherent analysis is:

let the covariance matrix M have J eigenvalues ordered from big to smallλ ₁ ，λ ₂ ，…，λ _J The corresponding feature vector is v ₁ ，ν ₂ ，…，ν _J . The coherent volume calculation formula based on the feature structure is:

wherein lambda is _j And calculating eigenvalues of the covariance matrix of the seismic data in the time window. In theory, the third generation algorithm C3 based on the characteristic value is superior to the second generation algorithm C2 based on the multi-channel cross-correlation and the C1 based on the two-channel cross-correlation, mainly because the noise which is higher than the additional noise level is automatically removed when the signal is higher than the additional noise level, and the transverse resolution is greatly improved; but the C3 calculation amount is large.

The fine coherence technique utilizes multi-channel similarity to convert a three-dimensional data body into a correlation coefficient data body through calculation, and compares the similarity of local seismic waveforms (the point with a lower coherence value is better related to the discontinuity of the reflected wave waveform) through fine coherence processing of the three-dimensional data body, so that various geological phenomena can be intuitively distinguished.

The eigenvalue coherence technology is used as the latest coherence algorithm, can form covariance matrixes by using multi-channel seismic data, and obtains the correlation among the multi-channel data by using multi-channel feature decomposition technology.

S600, carrying out geostatistical trend inversion on the basis of wave impedance inversion of the hole reservoir, so as to improve the identification and characterization capability of the small fracture-hole reservoir, and combining with S500, namely completing the characterization of the inside of the fracture cavity reservoir.

In order to further improve the identification and characterization capability of the seismic inversion result on the small fracture-cavity reservoir, the invention adopts geology statistical trend inversion. The trend analysis method is a method for simulating (or fitting) the spatial distribution of the geographic data and the regional variation trend thereof by using a mathematical method to decompose trend values and residual values in the actual data of the spatial geographic element distribution to construct a mathematical model. The specific principle of the method can be divided into two steps:

S601, establishing a trend model

Let the actual observation data of a certain spatial geographic element be: z is Z _i (x _i |y _i ) (i=1, 2, …, n) with a trend fit of Z _i (x _i |y _i ) Then there is Z _i (x _i |y _i )＝Z^ _i (x _i |y _i )+∈ _i Wherein: epsilon _i I.e. the residual value (residual value). When (x _i |y _i ) When the geographic elements are spatially changed, the interactive relation among the actual distribution curved surface, the trend surface and the residual surface of the geographic elements is described.

The trend surface (or body) is calculated from the actual observed value, and a regression analysis method is generally adopted, so that the sum of squares of residual errors tends to be minimum, namely:

according to the observed value x _i ,y _i ,z _i (i=1, 2, …, n) determining the coefficients a of the polynomial ₀ ，a ₁ ，…，a _p The sum of squares of the residuals is minimized.

S602, moderately checking trend surface model

The fitting degree coefficient R2 of the trend surface and the actual surface is an important index for measuring the fitting goodness of the regression model. The goodness of fit of the regression model is generally expressed as the sum of squares of the total dispersion of the variable z and the proportion of the sum of squares of the regression. The larger R2, the higher the fitting degree of the trend surface.

F, testing the significance of the whole trend surface regression model. And determining whether the regression relationship between the variable z and the independent variables x and y is obvious or not by using the ratio of the regression square sum to the residual square sum in the total dispersion square sum of the variable z.

Solving the difference between the regression square sum of the higher order polynomial equation and the regression square sum of the lower order polynomial equation; dividing the difference by the difference in degrees of freedom of the sum of squares of the regression to obtain a regression mean square error due to the increase of the polynomial degree; dividing the mean square error by the residual mean square error of the higher order polynomial to obtain a moderately comparative test value F of the trend surface models of the two successive orders. If the resulting F value is significant, the higher order polynomial makes a new contribution to regression; if the F value is not significant, the higher order polynomials do not make a new contribution to the regression.

The trend analysis technology is introduced into the formation uniformity research and mainly based on the heterogeneity of carbonate rock, the background value of the whole impedance is given, the opposite low-resistance area is searched in the high-low main impedance area to indicate the favorable development area of the small fracture hole, the interference of the background impedance value on the information of the small fracture hole is eliminated, and the favorable area can be indicated more accurately.

The trend body generation method in the geology statistics trend inversion comprises the following steps: and carrying out large grid filtering on the absolute impedance attribute body to form a gradual change trend impedance trend body, and then subtracting the original absolute impedance attribute body from the gradual change trend body to obtain a residual impedance body, wherein the residual impedance body represents the residual group of antibodies in the form of conventional seismic attributes.

According to the principle of the trend surface method, the absolute impedance is subtracted from the trend body to obtain a relatively low geology statistics trend wave impedance inversion body.

Application example 1

The preferable scheme is adopted to describe the fracture cavity reservoir boundary line by carrying out configuration attribute on the full-depth 1 well.

The seismic profile for a full depth 1 well is shown in FIG. 3a, the profile is shown in FIG. 3b, and the profile threshold is ultimately used to determine the profile of the full depth zone reservoir is shown in FIG. 3 c.

The transverse breaking range of the fracture is marked by using the curve of the sidetracking horizontal well during drilling, most sidetracking horizontal wells which are just broken are empty and lost, no logging curve is usually used, and the width of the broken core and the breaking zone range can be judged by using the change rule of the curve during drilling. The drilling curve of the bedrock generally shows a stable change rule, when the drilling curve enters the range of the broken belt, the drilling curve changes, according to different wells, the drilling curve shows different change trends, but the whole of the drilling curve shows linear change characteristics, and the broken corrosion belt shows low drilling and is accompanied by the phenomenon of emptying leakage.

And calibrating a threshold value of gradient configuration attribute by using a full-depth 1 well drilling lost circulation point, wherein the plane accords with the overall geological recognition. In the drilling process of the full-depth 1 well, three times of emptying leakage occur, and through geological analysis, the primary leakage point is positioned at the fracture of the trunk, the secondary leakage point is positioned in a miniature hole in the fracture cavity reservoir, and the three times of leakage are positioned at the bottom of the well.

Fig. 4 is a quantitative characterization three-dimensional representation of a full depth block configuration attribute fracture cavity reservoir, and fig. 5 a-5 c are graphical representations of full depth 1 well fracture cavity reservoir boundary lines. The gradient topography is used to quantitatively engrave the spatial morphology of the full deep-block fracture cavity reservoir and predict volume and resource quantity. According to the embodiment, VVA software is used for quantitative depiction, after seismic data of gradient configuration attributes are input, data meeting a threshold range are displayed by using a small grid of 25m multiplied by 5ms according to the threshold value of the gradient configuration attributes, and data not meeting the threshold value range are removed completely, so that a three-dimensional reservoir model is formed.

Application example 2

And (5) describing the internal structure of the full-depth fracture cavity reservoir by using the discontinuity attribute and the geology statistical trend inversion double parameters.

Fracture cavity internal crack structure: the discontinuity property may be used to detect a boundary, such as a discontinuity in the phase axis amplitude in the horizontal direction. If a short calculation window is applied, the discontinuity properties will have a good reflection for some formation deposit characteristics, such as reef phase, river course, breach fan, etc. The cracking degree of the crack structure in the cracking zone can be better reacted in carbonate rock. The fine coherence technique utilizes multi-channel similarity to convert a three-dimensional data body into a correlation coefficient data body through calculation, compares the similarity of local seismic waveforms (the point with a lower coherence value is better related to the discontinuity of a reflected wave waveform) through fine coherence processing of the three-dimensional data body, and can intuitively distinguish geological phenomena related to fracture, crack, sedimentary facies, lithology changes, even fluid changes and the like. The discontinuous coherence technology is used as the latest coherence algorithm, can form covariance matrixes by using multiple channels of seismic data, and obtains the correlation among the multiple channels of data by using a multiple channel feature decomposition technology.

FIG. 6a is a full depth seismic section, FIG. 6b is a full depth discontinuity profile, and FIG. 6c is a full depth discontinuity and seismic overlap profile; from the full depth seismic section shown in the figure, the full depth trunk has longer fracture extension length, more secondary fracture development is associated, the extension length is shorter, and the continuity is poor. The advantage of the discontinuity property is that the method is very advantageous for detecting the fracture with relatively large scale, the fracture grid is very clear, and the method is one of the common technical methods with better effect for fracture identification at present.

Inversion of geologic statistics trend: the outstanding is the comprehensive response of various high-quality reservoirs formed by broken bands, corrosion, hydrothermal corrosion transformation and the like caused by breaking activities, and the reflection of the unobvious beads formed by the holes in the broken cavity is clearer.

FIG. 7a is a full depth segmented seismic profile, FIG. 7b is a conventional wave impedance inversion profile, FIG. 7c is a geostatistical trend inversion profile; as can be seen from the figure, the reservoir in the area has various types and complex causes, and mainly comprises later-stage corrosion hole seams and structural cracks, and also comprises a few primary pores. The development degree of the reservoir is closely related to factors such as a layer system, lithology, a sedimentary facies belt, structural curvature, fracture and the like. Inversion results show that the reservoir develops in one room group and the hawk mountain group, and the reservoir is mainly provided with one room group and has strong non-homogeneous transverse characteristic. The reservoir types are various, both cave and hole and fracture reservoirs are developed, and inversion results indicate that the reservoirs are generally distributed 0-300 ms below the top surface of a compartment set.

The cave type reservoir is generally in a bead shape or a sheet shape on a post-earthquake stack section and an inversion section, mainly develops in a room and an eagle mountain group in the vertical direction, and has the characteristics of multi-stage and multi-layer development. The pore type reservoir layer widely develops in the area, and is reflected in three seismic phases of a bead phase, a lamellar phase and a clutter phase on a seismic section, and the lamellar phase and the clutter phase are the main; the inversion section is in a strip shape or a layered spread.

The whole full-depth block is positioned in the surface karst superposition fracture transformation area, and the fracture control and storage characteristics are obvious. On the plane, the fracture-cavity type reservoir layer is developed along the sliding fracture zone, is far away from the fracture zone and is underdeveloped, and is mainly controlled by the northeast sliding fracture. The inversion result has good consistency with the root mean square amplitude, the amplitude change rate attribute and the fracture prediction result, and the inversion result is reliable.

Fig. 8a and 8b are three-dimensional representations of the fracture cavity delineation of the obliquely-split overlapping extrusion section 1, wherein fig. 8a is a fracture cavity boundary line engraving (topographical attribute) and fig. 8b is a two-parameter characterization fracture cavity internal structure (geostatistical trend inversion+discontinuity attribute). Fig. 9a and 9b are three-dimensional representations of a broken cavity depiction of a pigtail pull segment 2, wherein fig. 9a is a broken cavity boundary line engraving (topographical properties) and fig. 9b is a dual parameter characterization broken cavity internal structure (residual wave impedance + discontinuity properties). From the figure, the advantage of using the dual parameter pore reservoir and the discontinuous property to characterize the internal structure of the reservoir is that: (1) the structure of different types of reservoirs can be clearly described; (2) the lower longitudinal limit of the reservoir can be predicted, the possible hydrocarbon height can be predicted, and the full depth 1 predicts the hydrocarbon height 840m.

According to the analysis and research of the internal structure of the fracture cavity reservoir, the main cave reservoir is distributed near the sliding fracture, the fracture reservoir has good fracture correlation and is distributed in the inclined overlapping extrusion section and the plait section; the inclined-column overlapped extrusion section and the plait-shaped pulling section are divided into high-quality reservoir distribution areas, and the linear section reservoir is relatively poor.

It should be understood that the foregoing examples of the present invention are provided merely for clearly illustrating the present invention and are not intended to limit the embodiments of the present invention, and that various other changes and modifications may be made therein by one skilled in the art without departing from the spirit and scope of the present invention as defined by the appended claims.

Claims (5)

1. The method for describing the reservoir stratum of the sliding fracture cavity in the basin area comprises the following steps of:

s100, acquiring drilling data, logging data and seismic data of a basin area sliding fracture cavity reservoir;

s200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing bead space distribution by adopting a wavelet decomposition-based frequency-reduction speed inversion method;

S300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify a fracture-cavity abnormal body;

s400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir;

s500, calculating the discontinuous attribute of the three-dimensional earthquake to describe a crack structure in a fracture cavity reservoir;

s600, carrying out geostatistical trend inversion on the basis of wave impedance inversion of the hole reservoir so as to improve the identification and characterization capability of the small fracture hole reservoir, and completing the characterization of the inside of the fracture cavity reservoir by combining with S500;

the process of calculating gradient profile attributes in S300 includes:

s301, establishing gradient configuration attributes:

calculating gradient vectors of each point of the three-dimensional seismic data volume:

wherein x, y and z are the line, channel position and double-journey travel time direction of the three-dimensional earthquake voxels respectively; g is a gradient vector; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; u is the derivative matrix in the corresponding direction;

smoothing the three-dimensional data volume gradient with a three-dimensional gaussian window:

wherein G is a Gaussian kernel function; g is the amplitude gradient field calculated for the data volume; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; 1. 2 and 3 respectively refer to x, y and z directions; sigma (sigma) _g Is the length corresponding to the i direction;

s302, calculating gradient configuration attribute T:

s303, smoothing each component of the gradient configuration attribute by using a Gaussian window:

s304, calculating characteristic values and characteristic vectors of gradient configuration attributes:

wherein λ and v are a eigenvalue and eigenvector, respectively; lambda (lambda) ₁ ≥λ ₂ ≥λ ₃ > 0, three eigenvectors v ₁ 、v ₂ 、v ₃ Two pairs of the two orthogonal directions, 1, 2 and 3 respectively refer to x, y and z directions; the method comprises the following steps:

1)λ ₁ ＝λ ₂ ＝λ ₃ > 0, in this case isotropic;

2)λ ₁ ＞λ ₂ ＝λ ₃ =0, at this time λ ₁ The corresponding vector direction approaches the Z direction;

3)λ ₁ ＞λ ₂ ＞λ ₃ =0, in which case the filter vector is accurately calculated from the eigenvalues and eigenvectors, while preserving valid breaks;

the step of calculating the seismic discontinuity attributes in S500 includes:

the data of one analysis channel consisting of N sample point elements in a time window and J channels in the window taking the channel as a unit are arranged together to form a data matrix:

wherein u is _ij Representing the amplitude value of the ith seismic wave on the jth trace; n rows of matrix U areAnalyzing the nth amplitude value on each seismic trace in the body, and then the variance matrix of the nth sample is:

so that covariance matrix formed by covariance of each element in the time window of coherent analysis is:

Let the covariance matrix M have J eigenvalues lambda ordered from big to small ₁ ，λ ₂ ，...，λ _J The corresponding feature vector is v ₁ ，v ₂ ，...，v _J The coherent volume calculation formula based on the feature structure is:

wherein lambda is _j Calculating eigenvalues of a covariance matrix of the seismic data in the time window;

the process of inversion of the geostatistical trend in S600 includes:

s601, establishing a trend model;

let the actual observation data of a certain spatial geographic element be: z is Z _i (x _i |y _i ) (i=1, 2,., n), the trend fit value isThere is->Wherein: epsilon _i Namely, a residual value, namely, a residual value;

calculating a trend surface or a trend body from an actual observed value, and adopting a regression analysis method to enable the square sum of residual errors to be minimum, namely:

according to the observed value x _i ，y _i ，z _i (i=1, 2,) n) determining coefficients a of the polynomial ₀ ，a ₁ ，...，a _p Minimizing the sum of squares of the residuals;

s602, moderately inspecting a trend surface model, wherein the moderately inspecting comprises an R2 and F coefficient of fitting degree of a trend surface and an actual surface;

the trend body generation method in the geostatistical trend inversion comprises the following steps:

performing large grid filtering on the absolute impedance attribute body to form a gradual change trend impedance trend body, and subtracting the original absolute impedance attribute body from the gradual change trend body to obtain a residual impedance body, wherein the residual group of antibodies are expressed in the form of conventional seismic attributes;

According to the principle of the trend surface method, the absolute impedance is subtracted from the trend body to obtain a relatively low geology statistics trend wave impedance inversion body.

2. The method for characterizing a basin-zone slip fracture cavity reservoir according to claim 1, wherein S200 specifically comprises:

s201, establishing an initial wave impedance low-frequency model

Establishing a space wave impedance model by combining logging data and seismic data; obtaining a layer speed by utilizing the superposition speed in seismic processing and through speed conversion, and then converting the density by utilizing a Gardner formula to finally obtain a wave impedance model;

s202, inversion processing

The method comprises the steps of restraining sparse pulse inversion, wherein each channel of initial wave impedance calculated according to an objective function is adjusted, and the initial wave impedance comprises the adjustment of reflection coefficients, so that a full-area relative wave impedance data body is obtained;

the high frequency band is from a constraint sparse pulse inversion algorithm, and the low frequency band is combined by the wave impedance model established in the step S201.

3. The bench basin area slip fracture cavity reservoir characterization method of claim 1, wherein the fracture lateral fracture extent is calibrated using a time-of-drilling curve of the sidetrack horizontal well in S400 to determine a profile attribute threshold range.

4. A counter basin area slip fracture cavity reservoir characterization method according to claim 3 wherein said determining a profile attribute threshold range comprises:

The drilling curve changes and shows a linear change characteristic, the broken corrosion zone shows a low drilling state and is accompanied by a blow-out leakage phenomenon, the position is used for calibrating a reservoir threshold value, namely the reservoir threshold value is used as a reservoir boundary point, the reservoir is arranged inside, and the matrix is arranged outside.

5. The counter basin area walk-slip fracture cavity reservoir characterization method of claim 1, wherein the counter basin area walk-slip fracture cavity reservoir characterization method comprises the steps of:

s100, acquiring drilling data, logging data and seismic data of a basin area sliding fracture cavity reservoir;

s200, carrying out wave impedance inversion on a hole type reservoir in a fracture cavity reservoir, and describing bead space distribution by adopting a wavelet decomposition-based frequency-reduction speed inversion method; the method specifically comprises the following steps:

s201, establishing an initial wave impedance low-frequency model

Establishing a space wave impedance model by combining logging data and seismic data; obtaining a layer speed by utilizing the superposition speed in seismic processing and through speed conversion, and then converting the density by utilizing a Gardner formula to finally obtain a wave impedance model;

s202, inversion processing

The method comprises the steps of restraining sparse pulse inversion, wherein each channel of initial wave impedance calculated according to an objective function is adjusted, and the initial wave impedance comprises the adjustment of reflection coefficients, so that a full-area relative wave impedance data body is obtained;

The high frequency band comes from a constraint sparse pulse inversion algorithm, and the low frequency band is combined by a wave impedance model established in S201;

s300, calculating gradient configuration attributes of the three-dimensional seismic data body to identify a fracture-cavity abnormal body;

wherein the process of calculating gradient profile attributes comprises:

s301, establishing gradient configuration attributes:

calculating gradient vectors of each point of the three-dimensional seismic data volume:

wherein x, y and z are the line, channel position and double-journey travel time direction of the three-dimensional earthquake voxels respectively; g is a gradient vector; g ₁ 、g ₂ 、g ₃ Respectively along the x, y and z directionsDirectional derivative; u is the derivative matrix in the corresponding direction;

smoothing the three-dimensional data volume gradient with a three-dimensional gaussian window:

wherein G is a Gaussian kernel function; g is the amplitude gradient field calculated for the data volume; g ₁ 、g ₂ 、g ₃ Directional derivatives along the x, y, z directions, respectively; 1. 2 and 3 respectively refer to x, y and z directions; sigma (sigma) _g Is the length corresponding to the i direction;

s302, calculating gradient configuration attribute T:

s303, smoothing each component of the gradient configuration attribute by using a Gaussian window:

s304, calculating characteristic values and characteristic vectors of gradient configuration attributes:

wherein λ and v are a eigenvalue and eigenvector, respectively; lambda (lambda) ₁ ≥λ ₂ ≥λ ₃ > 0, three eigenvectors v ₁ 、v ₂ 、v ₃ Two pairs of the two orthogonal directions, 1, 2 and 3 respectively refer to x, y and z directions; the method comprises the following steps:

1)λ ₁ ＝λ ₂ ＝λ ₃ > 0, in this case isotropic;

2)λ ₁ ＞λ ₂ ＝λ ₃ =0, at this time λ ₁ The corresponding vector direction approaches the Z direction;

3)λ ₁ ＞λ ₂ ＞λ ₃ =0, in which case the filter vector is accurately calculated from the eigenvalues and eigenvectors, while preserving valid breaks;

s400, using the threshold value range of the well-drilled calibration configuration attribute to delineate the fracture zone boundary of the fracture cavity reservoir;

the method comprises the steps of calibrating a fracture transverse crushing range by using a drilling curve of a sidetracking horizontal well so as to determine a configuration attribute threshold value range; the process for determining the configuration attribute threshold value range comprises the following steps:

the drilling curve changes and shows a linear change characteristic, the broken corrosion zone shows a low drilling state and is accompanied by a blow-out leakage phenomenon, the position is used for calibrating a reservoir threshold value, namely the reservoir threshold value is used as a reservoir boundary point, the reservoir is arranged inside, and the matrix is arranged outside;

s500, calculating the discontinuous attribute of the three-dimensional earthquake to describe a crack structure in a fracture cavity reservoir;

wherein the step of calculating the seismic discontinuity attributes comprises:

the data of one analysis channel consisting of N sample point elements in a time window and J channels in the window taking the channel as a unit are arranged together to form a data matrix:

Wherein u is _ij Representing the amplitude value of the ith seismic wave on the jth trace; n rows of matrix U areAnalyzing the nth amplitude value on each seismic trace in the body, and then the variance matrix of the nth sample is:

so that covariance matrix formed by covariance of each element in the time window of coherent analysis is:

let the covariance matrix M have J eigenvalues lambda ordered from big to small ₁ ，λ ₂ ，...，λ _J The corresponding feature vector is v ₁ ，v ₂ ，...，v _J The coherent volume calculation formula based on the feature structure is:

wherein lambda is _j Calculating eigenvalues of a covariance matrix of the seismic data in the time window;

s600, carrying out geostatistical trend inversion on the basis of wave impedance inversion of the hole reservoir so as to improve the identification and characterization capability of the small fracture hole reservoir, and completing the characterization of the inside of the fracture cavity reservoir by combining with S500;

wherein, the process of geostatistical trend inversion comprises:

s601, establishing a trend model;

let the actual observation data of a certain spatial geographic element be: z is Z _i (x _i |y _i ) (i=1, 2,., n), the trend fit value isThere is->Wherein: epsilon _i Namely, a residual value, namely, a residual value;

calculating a trend surface or a trend body from an actual observed value, and adopting a regression analysis method to enable the square sum of residual errors to be minimum, namely:

according to the observed value x _i ，y _i ，z _i (i=1, 2,) n) determining coefficients a of the polynomial ₀ ，a ₁ ，...，a _p Minimizing the sum of squares of the residuals;

s602, moderately testing the trend surface model, wherein the moderately testing comprises the degree of fit coefficient R2 testing and F testing of the trend surface and the actual surface.