Fractal feature prediction method for buried hill fracture system

Abstract

The invention relates to a fractal feature prediction method for a buried hill fracture system. The method comprises the following steps: S1, acquiring original seismic data of a buried hill region; s2, according to the original seismic data, calculating a curvature body attribute, a maximum likelihood body attribute and an ant body attribute of the buried hill reservoir, and extracting horizon slices of the three attributes; s3, according to the horizon slices of the three attributes of the buried hill reservoir, fusing by using a wavelet transform method to obtain seismic multi-scale fracture prediction plane attributes; s4, according to seismic multi-scale fracture prediction plane attributes, utilizing a shear wave transformation method to extract fracture zone vectors in the seismic multi-scale fracture prediction plane attributes; s5, calculating a crack fractal dimension corresponding to each plane grid according to a grid coverage method and the extracted crack zone vector; and S6, performing two-dimensional interpolation and smoothing on the fracture fractal dimension corresponding to each plane grid to obtain fracture system fractal dimension plane attributes.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration and development, and in particular to a method for predicting fractal characteristics of a buried hill fracture system. Background Art

[0002] The Archean large-scale metamorphic buried-hill oil and gas reservoirs in the Bohai Sea have become an important area for increasing offshore reserves and production. The multi-stage tectonic movement of the buried-hill has created a complex fracture system, and the development of fractures is the key to the effective development of the fractured reservoirs in the buried-hill. Due to the complex geological conditions and data quality constraints, the effective evaluation of the fracture system faces huge technical challenges.

[0003] In the existing technology, fractal methods are widely used in the field of earth sciences, in which the complex morphology and distribution of fractures have obvious fractal characteristics. The larger the fractal dimension, the more complex the underground fracture system, and the potential reservoir thickness, fracture network connectivity, production capacity and other indicator parameters may be better.

[0004] However, the existing calculation of the fractal dimension of the fracture system mainly uses the interpreted planar fault information or the logging fracture imaging interpretation under dense well network conditions. Due to the insufficient accuracy of the data basis interpretation and the large well spacing, the calculated fractal dimension has limited guiding significance for well deployment and production capacity design. Summary of the invention

[0005] In view of the above problems, the present invention provides a method for predicting fractal characteristics of a buried-hill fracture system, which uses seismic multi-scale prediction of fracture attributes to calculate fractal dimensions and achieve effective evaluation of the complexity of the buried-hill fracture system.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] In a first aspect, the present application provides a method for predicting fractal characteristics of a buried hill fracture system, the method comprising:

[0008] S1, obtain the original seismic data of the buried hill area;

[0009] S2, based on the original seismic data, calculate the curvature volume attributes, maximum likelihood volume attributes, and ant volume attributes of the buried hill reservoir, and extract the slices along the layer of the above three attributes;

[0010] S3, based on the slices along the layer of the three attributes of the buried hill reservoir, the wavelet transform method is used to fuse them to obtain the plane attributes of seismic multi-scale fracture prediction;

[0011] S4, based on the plane attributes of earthquake multi-scale fracture prediction, the fracture zone vector is extracted using the shear wave transform method;

[0012] S5, calculating the fracture fractal dimension corresponding to each plane grid according to the grid covering method and the extracted fracture zone vector;

[0013] S6, perform two-dimensional interpolation and smoothing on the fracture fractal dimension corresponding to each plane grid to obtain the fractal dimension plane attribute of the fracture system.

[0014] In one implementation, in S3, the fusion processing includes:

[0015] The curvature volume attribute, maximum likelihood volume attribute, and ant volume attribute are decomposed in multiple levels to obtain low-frequency sub-bands and high-frequency sub-bands of the three attributes;

[0016] According to the preset fusion rule, the low-frequency sub-band and the high-frequency sub-band are superimposed according to the designed fusion weight.

[0017] In one implementation, the process of determining the fusion weight of the low frequency sub-band includes:

[0018] A low-pass filter is used to perform low-pass filtering on the curvature body attribute slices, the maximum likelihood body attribute slices, and the ant body attribute slices, and the root mean square values ​​of the curvature body attribute slices, the maximum likelihood body attribute slices, and the ant body attribute slices after low-pass filtering are calculated;

[0019] R1, R2, and R3 are the root mean square values ​​of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice after low-pass filtering, respectively. The low-frequency fusion weights of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice are: 1 / R1: 1 / R2: 1 / R3.

[0020] In one implementation, the process of determining the fusion weight of the high frequency sub-band includes:

[0021] Use a bandpass filter to perform bandpass filtering on the curvature body attribute slices, maximum likelihood body attribute slices, and ant body attribute slices, and calculate the standard deviation values ​​of the curvature body attribute slices, maximum likelihood body attribute slices, and ant body attribute slices after bandpass filtering;

[0022] σ1, σ2, and σ3 are the standard deviation values ​​of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice after bandpass filtering, respectively. The high-frequency fusion weights of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice are: 1 / σ1: 1 / σ2: 1 / σ3.

[0023] In one implementation, the S5 includes: designing a plane grid size and a sliding step size;

[0024] Specifically: Calculate the distribution range histogram of the plane extension length of the interpreted fault in the study area, and calculate L mean , L Δ are the mean and standard deviation of the plane extension length of the fault, respectively, and the plane grid size is determined as L mean-2·L Δ , to ensure that all faults with a non-small probability of development length can be covered by the grid, and the sliding step length is determined to be L Δ / 4; in S6, the crack fractal dimension corresponding to each plane grid is placed at the exact center of the plane grid to form two-dimensional scattered points of the fractal dimension. The total number of points of the two-dimensional scattered points of the fractal dimension is recorded as m·n, where m is the number of fractal dimension points in the horizontal direction and n is the number of fractal dimension points in the vertical direction. The two-dimensional cubic spline interpolation method is used to interpolate the m·n two-dimensional scattered points of the fractal dimension into an M·N two-dimensional grid, where M is the number of horizontal sample points of the fractal dimension plane attribute and N is the number of vertical sample points of the fractal dimension plane attribute. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 A flow chart of a method for calculating the fractal dimension of a buried hill fracture system based on earthquake multi-scale fracture prediction attributes provided for the implementation of the present invention;

[0026] Figure 2 The slice diagram of maximum likelihood volume attributes, curvature volume attributes, and ant volume attributes calculated for the buried hill target layer of Bohai A oilfield;

[0027] Figure 3 This is the plane attribute map of earthquake multi-scale fracture prediction in Bohai A oilfield;

[0028] Figure 4 Extract fracture zone vector diagram from seismic multi-scale fracture prediction plane attributes for Bohai A Oilfield;

[0029] Figure 5 This is a schematic diagram of the grid covering method;

[0030] Figure 6 Schematic diagram for calculating fractal dimension;

[0031] Figure 7 This is the fractal dimension plane attribute diagram of the fracture system of Bohai A oilfield. DETAILED DESCRIPTION

[0032] In order to make the purpose, technical solution and advantages of the embodiment of the present invention clearer, the technical solution of the embodiment of the present invention will be clearly and completely described below in conjunction with the drawings of the embodiment of the present invention. Obviously, the described embodiment is a part of the embodiment of the present invention, not all of the embodiments. Based on the described embodiment of the present invention, all other embodiments obtained by ordinary technicians in this field belong to the scope of protection of the present invention.

[0033] In view of the defects and problems of the prior art, the present application provides a method for predicting fractal characteristics of a buried hill fracture system, the method comprising:

[0034] S1, obtain the original seismic data of the buried hill area;

[0035] S2, based on the original seismic data, calculate the curvature volume attributes, maximum likelihood volume attributes, and ant volume attributes of the buried hill reservoir, and extract the slices along the layer of the above three attributes;

[0036] S3, based on the slices along the layer of the three attributes of the buried hill reservoir, the wavelet transform method is used to fuse them to obtain the plane attributes of seismic multi-scale fracture prediction;

[0037] S4, based on the plane attributes of earthquake multi-scale fracture prediction, the fracture zone vector is extracted using the shear wave transform method;

[0038] S5, calculating the fracture fractal dimension corresponding to each plane grid according to the grid covering method and the extracted fracture zone vector;

[0039] S6, perform two-dimensional interpolation and smoothing on the fracture fractal dimension corresponding to each plane grid to obtain the fractal dimension plane attribute of the fracture system.

[0040] The above method flow is described in more detail in the following embodiments in conjunction with the accompanying drawings of the present invention.

[0041] like Figure 1 As shown in the flow chart, taking the Bohai A oilfield as an example, the method of the present invention is used to realize the fractal dimension prediction of the fracture reservoir of the oilfield. The specific steps are:

[0042] (1) The curvature volume attributes, maximum likelihood volume attributes, and ant volume attributes of the buried hill reservoir in oil field A are calculated. The seismic interpretation horizon of the Archean buried hill top is used, and the vertical time window range from the Archean buried hill top to 120 ms below the top surface is used to extract the slices along the layer of the above three attributes, as shown in the following figure: Figure 2 As shown;

[0043] (2) The three attribute slices are fused by wavelet transform method to obtain the fusion results of the three attribute slices as the plane attribute map of earthquake multi-scale fracture prediction, such as Figure 3 As shown;

[0044] The wavelet transform method refers to: using Haar wavelet transform to decompose the curvature volume attribute, maximum likelihood volume attribute, and ant volume attribute slices in multiple levels to obtain low-frequency subband coefficients of the three attributes and high-frequency subband coefficients in three directions (horizontal direction, vertical direction, and diagonal direction).

[0045] The fusion rules of low-frequency sub-band coefficients and high-frequency sub-band coefficients are designed, and the low-frequency and high-frequency sub-band coefficients of the decomposed curvature volume attribute, maximum likelihood volume attribute, and ant volume attribute slices are fused according to the low-frequency fusion rules and high-frequency fusion rules respectively, to obtain the low-frequency sub-band and high-frequency sub-band fusion results in the Haar wavelet domain. Finally, the Haar wavelet inverse transform is used to obtain the final fusion result slices of the three attributes in the time-space domain.

[0046] (3) The fracture zone vector is extracted from the final fusion result slice of the three attributes in the time-space domain using shearlet transform, such as Figure 4 As shown, the vector image is a 0, 1 binary image. The pixels covered by black lines represent that there is a crack vector in the pixel (value is 1), and the pixels not covered by black lines are the background (value is 0);

[0047] The shearlet transform method refers to: using the classic shearlet transform algorithm, mapping the original image to different directions through the shear matrix, and establishing an adaptive multi-directional directional filtering operator, so as to extract the crack zone vector from the final fusion result of the three attributes in the time-space domain.

[0048] (4) Design the plane grid size and sliding step size, use the grid covering method to calculate the fractal dimension of the extracted fracture zone vector, and obtain the fracture fractal dimension corresponding to each plane grid, as follows: Figure 5 , Figure 6 As shown;

[0049] The grid covering method for calculating fractal dimension means: the side length of the plane grid of the crack zone vector diagram is L pixels, the initial radius of the sub-grid is designed to be r pixels, the number of covering times is n, the growth step is (Lr) / (n-1) pixels, and the sub-grid radius is R k =r+(k-1)·(Lr) / (n-1), k = 1 to n, for each crack zone vector plane grid, start from the lower left corner and start with R k Calculate R for the subgrid radius k ·R k The cumulative sum of pixel data of the vector graph in the square area. The specific calculation process is: sum all pixel data of the vector graph in the area, where the black covered pixel data is 1 and the white covered pixel data is 0. k The accumulated sum of pixel data in the corresponding square area is S k , calculate R k , S k The natural logarithm of LR k =ln(R k ),LS k =ln(S k ), where ln represents the natural logarithm operation. k Data is the horizontal axis, with LS k Data is the vertical axis, plotting LR k For LS k Scatter plot, perform linear regression on the scatter plot, and get the slope K of the straight line, which is the fractal dimension of the crack corresponding to the current plane grid.

[0050] Figure 5The process of subdividing a plane grid into subgrids when calculating the fractal dimension is shown. The initial radius of the grid is 1, the growth step is 1, the number of coverages is n=13, and the square areas corresponding to 1, 2, 3, ... 13 are the subgrids. Figure 6 Shows the Figure 5 The process of calculating the fractal dimension of a plane grid, the horizontal axis is Figure 5 The logarithm of the radius of the 13 sub-grids is taken, and the vertical axis is the logarithm of the cumulative sum of pixel data. The linear regression obtains that the slope of the straight line is K = 1.89, which is Figure 5 Fractal dimension of a midplane grid.

[0051] According to the sliding step length, the crack fractal dimension corresponding to each plane grid is calculated in turn.

[0052] (5) Perform two-dimensional interpolation and smoothing on the fractal dimensions of the cracks corresponding to each plane grid to obtain the fractal dimension plane properties of the crack system. Figure 7 The normalized fractal dimension plane attribute map of the fracture system of Bohai A Oilfield and the fractal dimension contour lines are displayed. The fractal dimension of the fractures in the darker areas is larger, indicating that the fracture network is more complex and the fractures are more developed, indicating that it is a high-quality reserve area.

[0053] In the several embodiments provided by the present invention, it should be understood that the disclosed method can be implemented in other ways. For example, the device embodiments described above are only schematic. For example, the division of the above units is only a logical function division. There may be other division methods in actual implementation. For example, multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0054] The above-mentioned integrated unit implemented in the form of a software functional unit can be stored in a computer-readable storage medium. The above-mentioned software functional unit is stored in a storage medium, including a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) or a processor (Processor) to perform some steps of the above-mentioned method of each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (Read-Only Memory, ROM), random access memory (Random Access Memory, RAM), disk or optical disk and other media that can store program codes.

[0055] The above are only preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for predicting fractal characteristics of a buried hill fracture system, characterized in that: The method comprises: S1, obtain the original seismic data of the buried hill area; S2, based on the original seismic data, calculate the curvature volume attributes, maximum likelihood volume attributes, and ant volume attributes of the buried hill reservoir, and extract the layer-by-layer slices of the above three attributes; S3, based on the slices along the layer of the three attributes of the buried hill reservoir, the wavelet transform method is used to fuse them to obtain the plane attributes of seismic multi-scale fracture prediction; S4, based on the plane attributes of earthquake multi-scale fracture prediction, the fracture zone vector is extracted using the shear wave transform method; S5, calculating the fracture fractal dimension corresponding to each plane grid according to the grid covering method and the extracted fracture zone vector; S6, perform two-dimensional interpolation and smoothing on the fracture fractal dimension corresponding to each plane grid to obtain the fractal dimension plane attribute of the fracture system.

2. The method for predicting fractal characteristics of buried hill fracture systems according to claim 1, characterized in that: In S3, the fusion processing includes: The curvature volume attribute, maximum likelihood volume attribute, and ant volume attribute are decomposed in multiple levels to obtain low-frequency sub-bands and high-frequency sub-bands of the three attributes; According to the preset fusion rule, the low-frequency sub-band and the high-frequency sub-band are superimposed according to the designed fusion weight.

3. The method for predicting fractal characteristics of buried hill fracture systems according to claim 2, characterized in that: The process of determining the fusion weight of the low frequency sub-band includes: A low-pass filter is used to perform low-pass filtering on the curvature body attribute slices, the maximum likelihood body attribute slices, and the ant body attribute slices, and the root mean square values ​​of the curvature body attribute slices, the maximum likelihood body attribute slices, and the ant body attribute slices after low-pass filtering are calculated; R1, R2, and R3 are the root mean square values ​​of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice after low-pass filtering, respectively. The low-frequency fusion weights of the curvature volume attribute slice, maximum likelihood volume attribute slice, and ant volume attribute slice are: 1 / R1: 1 / R2: 1 / R3.

4. The method for predicting fractal characteristics of buried hill fracture systems according to claim 2, characterized in that: The process of determining the fusion weight of the high frequency sub-band includes: Use a bandpass filter to perform bandpass filtering on the curvature body attribute slices, maximum likelihood body attribute slices, and ant body attribute slices, and calculate the standard deviation values ​​of the curvature body attribute slices, maximum likelihood body attribute slices, and ant body attribute slices after bandpass filtering; σ1, σ2, and σ3 are the standard deviations of the curvature body attribute slices, maximum likelihood body attribute slices, and ant body attribute slices after bandpass filtering. High frequency fusion weight It is: 1 / σ1: 1 / σ2: 1 / σ3.

5. The method for predicting fractal characteristics of buried hill fracture systems according to claim 1, characterized in that: The S5 includes: designing the plane grid size and the sliding step size; Specifically: Calculate the distribution range histogram of the plane extension length of the interpreted fault in the study area, and calculate L mean , L Δ are the mean and standard deviation of the plane extension length of the fault, respectively, and the plane grid size is determined as L mean -2·L Δ , to ensure that all faults with a non-small probability of development length can be covered by the grid, and the sliding step length is determined to be L Δ / 4; in S6, the crack fractal dimension corresponding to each plane grid is placed at the exact center of the plane grid to form two-dimensional scattered points of the fractal dimension. The total number of points of the two-dimensional scattered points of the fractal dimension is recorded as m·n, where m is the number of fractal dimension points in the horizontal direction and n is the number of fractal dimension points in the vertical direction. The two-dimensional cubic spline interpolation method is used to interpolate the m·n two-dimensional scattered points of the fractal dimension into an M·N two-dimensional grid, where M is the number of horizontal sample points of the fractal dimension plane attribute and N is the number of vertical sample points of the fractal dimension plane attribute.

Images