Method for calculating fracture density data volume

Abstract

The invention provides a method and equipment for calculating a fracture density data volume, a medium and a program product. The method comprises the following steps: acquiring normalized fracture density curves of high and low angles in a well and a preferable attribute data volume; a threshold value is set for the fracture density curve, a fracture phase is established, optimization is performed based on the fracture relative optimization attribute data volume to obtain a first attribute data volume, and a first fracture density calculation model is established to calculate the first fracture density data volume; and establishing a second fracture density calculation model for the first fracture density data volume, the first attribute data volume and the fracture density curve, calculating a second fracture density data volume of each fracture phase, and carrying out reverse normalization and weighted fusion processing on the second fracture density data volume to obtain a third fracture density data volume. The method can achieve the calculation of the fracture density values of different fracture development scales in the target interval, improves the fracture prediction precision and accuracy of the target interval of a related exploration area, reduces the drilling risk, and improves the economic benefits of oil-gas exploration.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil and gas exploration, and in particular to a method for calculating a fracture density data volume. BACKGROUND

[0002] Fractures are important pathways for oil and gas accumulation and migration in the subsurface, and fracture prediction refers to predicting the development strength of fractures or parameters related to fracture analysis. In most cases, fractures are mainly tectonic fractures, which are due to or associated with local tectonic events, including fault-related fracture systems, uplift-related fracture systems, and fold-related fracture systems. The high-quality shale of the marine Wufeng-Longmaxi Formation and the tight sandstone reservoirs in the southern Sichuan Basin are basically developed with micro-fractures, which play a role in the conventional and unconventional reservoirs. Therefore, finding fracture-type reservoirs is one of the important targets of marine and terrestrial oil and gas exploration. In general, fractures are a relatively important indicator for evaluating reservoir performance, and fracture prediction in reservoirs is also an important problem faced by geophysical exploration.

[0003] At present, most conventional fracture prediction techniques use pre-stack or post-stack seismic data to predict fractures, and some use finite element analysis, tectonic stress field analysis and other geological experience analysis techniques to predict fractures. In fracture prediction, multi-attribute fusion is also used to calculate fracture density. The inversion or attribute sensitive to fractures mainly includes coherence, curvature, P-wave anisotropy strength, amplitude and frequency attribute, etc., each of which has advantages and disadvantages, and the accuracy of micro-fracture prediction is relatively poor. In some fracture prediction, statistical methods are mainly used, and the main operation method is to cross the relevant elastic parameter curve with the fracture density, optimize the data points with relatively concentrated distribution, i.e. the elastic parameters with good correlation, and perform fitting to obtain the relevant relationship, and then convert the obtained elastic parameter data volume and relationship to obtain the fracture density data volume. Some patents, such as CN1737607A, entitled "Advantageous frequency band coherence processing method for fine fault interpretation", disclose a method for predicting small faults and fractures by using coherence processing. After the seismic profile is processed by the advantageous frequency band coherence technology, new missed small faults are displayed, the breakpoint position of the small faults is more accurate, and the extension direction is continuous. CN101907725A, entitled "Fracture prediction method and device", discloses a method for predicting the direction and density of fractures by picking up the reflection amplitude of each seismic trace in the target layer window, and performing elliptical fitting based on the obtained azimuth and reflection amplitude.

[0004] (1) The conventional fracture prediction technology is difficult to effectively describe the changes of micro-fractures in the vertical and horizontal directions, and the accuracy is often not ideal;

[0005] (2) The results of the conventional fracture prediction technology for micro-fracture prediction often do not match the measured results in the well, and the prediction accuracy is relatively low.

[0006] (3) The fracture density is calculated through a related linear fitting formula, which cannot eliminate the influence of some non-fracture factors, resulting in low accuracy of the prediction results. SUMMARY

[0007] The present application aims to provide a method, device, medium and program product for calculating fracture density data volume to solve the problems existing in the conventional fracture prediction technology.

[0008] In a first aspect, the present application provides a method for calculating fracture density data volume, comprising the following steps:

[0009] Step 100, establishing high and low angle fracture density curves in a known well, calculating related seismic attribute data volume and optimizing the seismic attribute data volume to obtain an optimized attribute data volume, and then normalizing the high and low angle fracture density curves in the well and the optimized attribute data volume to obtain normalized high and low angle fracture density curves in the well and the optimized attribute data volume;

[0010] Step 200, setting a related threshold value for the normalized high and low angle fracture density curves in the well respectively and establishing a fracture phase, optimizing the optimized attribute data volume based on the fracture phase to obtain a first attribute data volume, and establishing a related first fracture density calculation model to calculate the first fracture density data volume of each fracture phase of high and low angle fractures;

[0011] Step 300, establishing a second fracture density calculation model for the first fracture density data volume, the first attribute data volume and the fracture density curve and performing fracture density calculation to obtain a second fracture density data volume of each fracture phase, and then performing reverse normalization and weighted fusion processing on the second fracture density data volume to obtain a third fracture density data volume.

[0012] In some embodiments, step 100 comprises the following sub-steps:

[0013] Step 101, preparing geological data and logging data, and forming high and low angle fracture density curves according to the geological data and logging data respectively;

[0014] Step 102, preparing seismic data, and performing inversion and attribute extraction on the seismic data to obtain seismic inversion and attribute data volumes, and then extracting seismic inversion and attribute curve data of the target layer at each well point from the seismic inversion and attribute data volumes, and performing time-depth conversion and resampling calculation on the seismic inversion and attribute curve data to form a fracture density-seismic attribute correspondence between the seismic inversion and attribute curve data in the depth domain and the data on the fracture density curves corresponding to the well points;

[0015] Step 103, normalizing and calculating the correlation coefficients of the seismic inversion and attribute curve data of the target layer at each well point and the data on the high-angle and low-angle fracture density curves of the target layer at the corresponding well points, and selecting M and L seismic inversion and attribute data volumes with higher correlation coefficients as preferred attribute data volumes; wherein the M seismic inversion and attribute data volumes correspond to the data on the low-angle fracture density curves, and the L seismic inversion and attribute data volumes correspond to the data on the high-angle fracture density curves.

[0016] Step 104, obtaining respective normalized function formulas by normalizing the data on the high-angle and low-angle fracture density curves of the target layer at each well point and the M and L seismic inversion and attribute data volumes.

[0017] Step 105, performing normalization processing on the M and L seismic inversion and attribute data volumes according to the normalized function formulas to obtain normalized M and L seismic inversion and attribute data volumes.

[0018] In some embodiments, in step 103, the M and L seismic inversion and attribute data volumes with higher correlation coefficients refer to: first performing a first correlation coefficient calculation, then performing a second correlation coefficient calculation between the M and L seismic inversion and attribute data volumes with higher correlation coefficients obtained from the first correlation coefficient calculation, and again selecting M and L seismic inversion and attribute data volumes with lower correlation coefficients obtained from the second correlation coefficient calculation.

[0019] In some embodiments, in step 200, establishing a fracture facies includes: extracting the data on the high-angle and low-angle fracture density curves in the well after normalization, respectively, to establish high-angle and low-angle fracture density data sets of the target layer; setting a series of related threshold values for the high-angle and low-angle fracture density data sets of the target layer and dividing different fracture density ranges, so as to divide into different fracture density ranges; the fracture density data in different fracture density ranges belong to different fracture facies, so as to obtain each fracture facies with respect to high-angle and low-angle.

[0020] In some embodiments, in step 200, the first attribute data volume is obtained by selecting the preferred attribute data volume based on the fracture facies, including:

[0021] According to the fracture density data, a fracture density sample value of each fracture is determined, and a preferred attribute data volume corresponding to the fracture density sample value of each fracture facies is determined: a time-depth relationship of each well is determined according to acoustic wave and density logging curves in the well and post-stack seismic data; the curve data of the preferred attribute data volume is converted from the time domain to the depth domain by using the time-depth relationship, and the curve data of the preferred attribute data volume in the depth domain is resampled to make it consistent with the sampling interval of the fracture density curve, and a one-to-one correspondence between the data values is formed, and through the correspondence, the preferred attribute data volume corresponding to the fracture density sample value of the fracture facies is determined;

[0022] Through cross-plot analysis, the attribute curve of the preferred attribute data volume with the highest correlation coefficient of the fracture density curve in each well of the fracture facies is taken as the first attribute curve; then the first attribute curve and other preferred attribute curves are sequentially taken as a curve pair, and the best matching function determined by the least square method is used to calculate the fitting fracture density curve, and the correlation coefficient between the fitting fracture density curve and the fracture density curve of the fracture facies is calculated, and the attribute curve with the highest correlation coefficient is the second attribute curve; then the third attribute curve, the fourth attribute curve, …, and the Nth attribute curve are continuously searched for in the same way, and each attribute curve is the first attribute data volume.

[0023] In some embodiments, in the cross-plot analysis, if the data points of the preferred attribute data volume corresponding to the fracture density sample value of the fracture facies do not meet the correlation requirements, the fracture density sample value of the fracture facies is expanded to multiple fracture density values within a certain range.

[0024] In some embodiments, in step 200, establishing the first fracture density calculation model includes: determining the number of preferred attribute data volumes by using the method of blind well verification, that is, using each well as a verification well, using the number of attribute curves and the best matching function determined by other wells to calculate the fitting fracture density curve of the verification well, and calculating the error between the fitting fracture density curve and the fracture density curve of the fracture facies; after obtaining the error when each well is used as a verification well, the attribute curve combination and the best matching function with the lowest average verification error are the first fracture density calculation model of the fracture facies.

[0025] In some embodiments, in step 300, establishing the second fracture density calculation model includes: using the first fracture density data volume of each fracture facies, the second attribute data volume, and the corresponding fracture density curve to establish a corresponding attribute parameter cross-plot, and establishing a function formula related to fracture density calculation according to the result of the attribute parameter cross-plot, so as to establish the second fracture density calculation model; wherein the second attribute data volume is obtained by combining and calculating the first attribute data volume.

[0026] In some embodiments, the function is a probability density distribution function.

[0027] In some embodiments, the first attribute data volume is mathematically calculated in relation to two or three combinations to obtain a second attribute data volume.

[0028] In summary, due to the adoption of the technical solutions described above, the present application has the following beneficial effects:

[0029] 1. The present application can calculate the fracture density values of different fracture development scales in the target interval in a well, improve the fracture prediction accuracy and precision of the target interval in the relevant exploration area, and thus reduce the drilling risk and improve the economic benefits of oil and gas exploration.

[0030] 2. The present application has universality and can be applied to fracture prediction in shale gas, carbonate rock exploration and other aspects, and can also achieve good geological effects. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 The flowchart of the method for calculating the fracture density data volume in the embodiments of the present application. DETAILED DESCRIPTION

[0032] To make the objectives, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described in detail below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, but not all the embodiments of the present application. The components of the embodiments of the present application described and shown in the drawings herein can be arranged and designed in various different configurations.

[0033] Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present application without creative labor are within the scope of protection of the present application.

[0034] As shown in the drawings, Figure 1 The embodiments of the present application propose a method for calculating a fracture density data volume, which comprises the following steps:

[0035] Step 100: Establishing the high and low angle fracture density curves in the known well, calculating the related seismic attribute data volume, and optimizing the seismic attribute data volume to obtain an optimized attribute data volume, and then normalizing the high and low angle fracture density curves in the well and the optimized attribute data volume to obtain the normalized high and low angle fracture density curves in the well and the optimized attribute data volume;

[0036] Step 200, setting a correlation threshold value for the normalized high-angle and low-angle fracture density curves respectively and establishing a fracture phase, performing optimization on the preferred attribute data volume based on the fracture phase and the preferred attribute data volume, and establishing a first fracture density calculation model to calculate a first fracture density data volume for each fracture phase of the high-angle and low-angle fractures;

[0037] Step 300, establishing a second fracture density calculation model for the first fracture density data volume, the first attribute data volume and the fracture density curve and performing fracture density calculation to obtain a second fracture density data volume for each fracture phase, and then performing inverse normalization and weighted fusion processing on the second fracture density data volume to obtain a third fracture density data volume. The third fracture density data volume obtained by the above steps can be used to relatively accurately predict the fracture development in the research area by using seismic technology.

[0038] In some embodiments, step 100 includes the following sub-steps:

[0039] Step 101, preparing geological data and logging data, and forming high-angle and low-angle fracture density curves according to the geological data and the logging data.

[0040] The geological data includes core logging data, and the logging data includes sonic and density logging curves, FMI formation micro-resistivity scanning imaging data, etc. The geological data and the logging data can be obtained through geophysical exploration, core data and logging data.

[0041] The fracture density curve can be obtained through core observation, fracture data collection and post-processing of logging data. The fracture density curve is calculated according to a certain sampling interval. The number of fractures in the sampling interval is divided by the sampling interval to obtain a fracture density value, which is assigned to the midpoint depth of the sampling interval to obtain a depth value-fracture density value data pair. The fracture density curve is obtained by repeating the calculation. When the dip angle of the fracture surface is greater than or equal to 45°, the fracture is classified as a high-angle fracture and participates in the formation of a high-angle fracture density curve. Conversely, when the dip angle of the fracture surface is less than 45°, the fracture is classified as a low-angle fracture and participates in the formation of a low-angle fracture density curve.

[0042] Step 102, preparing seismic data, performing inversion and attribute extraction on the seismic data to obtain seismic inversion and attribute data volume, extracting seismic inversion and attribute curve data of the target layer at each well point from the seismic inversion and attribute data volume, and performing time-depth conversion and resampling calculation on the seismic inversion and attribute curve data to form a fracture density-seismic attribute correspondence between the data in the depth domain and the corresponding fracture density curve.

[0043] The seismic data is conventional three-dimensional pre-stack gather or post-stack seismic data volume.

[0044] The seismic inversion and attribute extraction can be realized by related geophysical exploration commercial software, such as the PAL module of the Landmark company, which can extract amplitude, frequency, and instantaneous seismic attributes from three-dimensional post-stack seismic data; the jason software can calculate wave impedance inversion data; the VVA software can extract curvature, coherence volume, and frequency volume data; and the FRS software can use three-dimensional pre-stack gather data to calculate P-wave anisotropy strength and extract P-wave anisotropy strength data.

[0045] In step 103, the seismic inversion and attribute curve data of the target layer at each well point are normalized (normalized to a certain value range) and correlation coefficient calculation is performed on the data on the high-angle and low-angle fracture density curves of the target layer at the corresponding well point, and M and L seismic inversion and attribute data volumes with higher correlation coefficients are selected as the preferred attribute data volumes for the next step, wherein the M seismic inversion and attribute data volumes correspond to the data on the low-angle fracture density curve, and the L seismic inversion and attribute data volumes correspond to the data on the high-angle fracture density curve.

[0046] In step 104, the data on the high-angle and low-angle fracture density curves of the target layer at the well point are normalized to obtain respective normalized function formulas.

[0047] In step 105, the M and L seismic inversion and attribute data volumes are normalized according to the normalized function formulas to obtain normalized M and L seismic inversion and attribute data volumes. In principle, the normalized data should be positively correlated with the M and L seismic inversion and attribute data volumes.

[0048] Preferably, in step 103, the M and L seismic inversion and attribute data volumes with higher correlation coefficients are obtained by first performing first correlation coefficient calculation, then performing second correlation coefficient calculation between the M and L seismic inversion and attribute data volumes with higher correlation coefficients (for example, greater than or equal to 0.5) obtained by the first correlation coefficient calculation, and then preferably selecting the M and L seismic inversion and attribute data volumes with lower correlation coefficients (for example, less than or equal to 0.3) obtained by the second correlation coefficient calculation.

[0049] Preferably, the normalization processing is a dimensionless processing method that changes the absolute value of the physical system data to a certain relative value relationship, and the normalization processing is a calculation using addition, subtraction, multiplication, division, or a combination thereof; in this embodiment, the normalization processing is specifically setting the sample data as x p(p = 1, 2,..., P), define x max = max{x p}, x min = min{x p}, the normalization processing calculation converts the data into the data in the interval of n~m, and the normalization processing calculation formula is as follows:

[0050]

[0051] In formula (1), x p is the data before normalization processing, x pi is the data after normalization processing, n and m are positive integers, and m>n≥0.

[0052] Preferably, the calculation formula of the correlation coefficient is:

[0053]

[0054] In formula (2), X i and Y i are the i-th data values of two kinds of data for correlation coefficient calculation, and are the average values of the rank order of the two kinds of data, respectively, and the value range of r is 0 to 1.

[0055] Preferably, the target interval refers to the interval where the reservoir is located, and the top and bottom positions of the target layer are obtained after the top and bottom positions are calibrated by the relevant well-seismic synthetic records, and the seismic data interpretation and related horizon data interpolation processing are performed.

[0056] In some embodiments, step 200 includes the following sub-steps:

[0057] Step 201, different threshold values are set for high-angle and low-angle fractures respectively, and each fracture facies is established according to the threshold values for the normalized high-angle and low-angle fracture density curves in the well.

[0058] In some embodiments, the main operation of step 201 is to extract the data on the normalized high-angle and low-angle fracture density curves in the well respectively, establish the high-angle and low-angle fracture density data sets of the target interval, set a series of related threshold values for the high-angle and low-angle fracture density data sets of the target interval respectively, and divide the different fracture density ranges, so as to divide into different fracture density ranges; the fracture density data in different fracture density ranges belongs to different fracture facies, so as to obtain each fracture facies of high-angle and low-angle.

[0059] Wherein, the setting of each threshold value of the high and low angle fracture density data set of the target interval should be determined according to the actual situation, expert experience and fracture density prediction accuracy, etc. For example, for low angle fractures, K fracture density ranges are divided for the fracture density data set and the establishment of relevant fracture facies is performed, thereby obtaining K fracture facies. In general, the number of known wells and data in the study area should be relatively sufficient to meet the requirements of the number of sampling points of the fracture density curve in different fracture facies of high and low angle fractures. In principle, the number of wells in the study area is generally required to be sixteen or more, the corresponding attribute data points of each fracture facies are required to be at least 20 or more, the number of fracture facies should be 2-3 times the number of sedimentary facies in the study area, and should be more than three or more.

[0060] Preferably, the number of threshold values should be greater than or equal to one, and the specific number and range of threshold values should be determined according to seismic data, fracture density development on well, expert experience and fracture density calculation accuracy, etc. For example, in the setting of a certain fracture facies, K i and K i+1 are threshold values and K i is less than K i+1 , then the fracture density value in the fracture facies should belong to [K i , K i+1 ]. In addition, the determination of relevant threshold values and numbers can also be implemented according to the histogram statistical results of the fracture density data values of each well.

[0061] Step 202, based on the fracture facies, the preferred attribute data body is obtained, and a relevant first fracture density calculation model is established to calculate the first fracture density data body of each fracture facies of high and low angle fractures. The specific operation is as follows:

[0062] First, the fracture density sample value of each fracture is determined according to the fracture density data, and the preferred attribute data body corresponding to the fracture density sample value of each fracture facies is determined. The implementation method is as follows: the time-depth relationship of each well is determined by calibrating the well-seismic according to the acoustic and density logging curves in the well and the post-stack seismic data; the curve data of the preferred attribute data body is converted from the time domain to the depth domain by using the time-depth relationship, and the curve data of the preferred attribute data body in the depth domain is resampled to make it consistent with the sampling interval of the fracture density curve, and a one-to-one correspondence between the data values is formed, and the preferred attribute data body corresponding to the fracture density sample value of the fracture facies is determined through the correspondence.

[0063] Wherein, each fracture facies corresponds to a fracture density sample value. The following methods are provided in this embodiment to determine the fracture density sample value:

[0064] ① The average of each fracture density data in the fracture phase is calculated as the fracture density sample value of the fracture phase.

[0065] ② The histogram of each fracture density data in the fracture phase is displayed and analyzed, and the fracture density data corresponding to the maximum peak value is selected as the fracture density sample value of the fracture phase.

[0066] ③ The fracture density sample value of each fracture phase is determined according to the actual fracture prediction and expert experience, etc.

[0067] In this way, the fracture density sample values of each fracture phase related to high and low angle fractures are determined through the above three methods.

[0068] Then, the first attribute data body is obtained by selecting the fracture phase based on the preferred attribute data body. This step mainly includes the following steps: the attribute curve of the preferred attribute data body with the highest correlation coefficient of the fracture density curve of each well in the fracture phase is selected as the first attribute curve through cross-plot analysis; then the first attribute curve and other preferred attribute curves are sequentially combined to form a curve pair, and the best matching function determined by the least square method is used to calculate the fitting fracture density curve and the correlation coefficient between the fitting fracture density curve and the fracture density curve of the fracture phase, and the attribute curve with the highest correlation coefficient is the second attribute curve; then the third attribute curve, the fourth attribute curve, …, the Nth attribute curve are sequentially searched according to the same method, and each attribute curve is the first attribute data body. It should be noted that if the data points of the preferred attribute data body corresponding to the fracture density sample value of the fracture phase do not meet the correlation requirements in the cross-plot analysis, the fracture density sample value of the fracture phase can be expanded to multiple fracture density values within a certain range to form a fracture density sample range, and the preferred attribute data body corresponding to each fracture density sample value in the fracture phase is used for attribute cross-plot in the two-dimensional coordinate system and as the first calculation fracture density model of the fracture phase. The calculation formula of the fracture density sample range of the fracture phase is as follows:

[0069]

[0070] In formula (3) and formula (4), is the maximum fracture density value of the fracture density sample range of the i th fracture phase, is the fracture density sample value of the i th fracture phase, and ΔF i is the fracture density increment of the fracture density sample range of the i th fracture phase, is the minimum fracture density value of the fracture density sample range of the i th fracture phase. The minimum fracture density value and the maximum fracture density value of the fracture density sample range of the related fracture phase should be included in the fracture density value range of the fracture phase.

[0071] Then, a first crack density calculation model is established. In actual operation, in order to avoid the problem of overfitting of the crack density curve, the number of preferred attribute data bodies can be determined by using the method of blind well verification, that is, each well is used as a verification well in turn, the number of attribute curves determined by other wells and the best matching function are used to calculate the fitting crack density curve of the verification well, and the error between the fitting crack density curve and the crack density curve of the crack phase is calculated. After obtaining the error when each well is used as a verification well, the attribute curve combination with the lowest average verification error and the best matching function are the first crack density calculation model of the crack phase, which is represented as:

[0072] F j = f1(A1, A2, …, A n ) (5)

[0073] In formula (5), F j is the first crack density data body of high and low angle cracks, A i is the i-th first attribute data body of the j-th crack phase, and f1 is the functional relationship between the crack phase and the corresponding first attribute data body curve (i.e. the best matching function determined by the least square method). Generally, the functional relationship includes linear relationship, polynomial relationship, exponential relationship and logarithmic relationship, etc. In different regions, the number of preferred attribute curves and the function relationship with the highest crack density correlation coefficient can be determined by using the least square method according to actual data.

[0074] Finally, the first crack density calculation model is used to calculate the first attribute data body to obtain the first crack density data body of each crack phase of high and low angle cracks.

[0075] It should be noted that the first crack density calculation model can also be a calculation model based on deep learning or neural network, or a combination of the two. The neural network can be a hybrid network that combines convolutional neural network and recurrent neural network. The related calculation model can be implemented by commercial software, and therefore will not be described in detail here.

[0076] In some embodiments, step 300 includes the following sub-steps:

[0077] Step 301, establishing a second crack density calculation model.

[0078] This step mainly uses the first crack density data body of each crack phase, the second attribute data body and the corresponding crack density curve to establish the corresponding attribute parameter crossplot, and establishes a function formula related to crack density calculation according to the result of the attribute parameter crossplot, thereby establishing a second crack density calculation model. The second attribute data body is obtained by combining and calculating the first attribute data body.

[0079] Generally, the second fracture density calculation model can use the corresponding relationship between the fracture density sample value of each fracture phase on the well and the corresponding first fracture density data body and second attribute data body to establish an attribute parameter crossplot of the corresponding two-dimensional coordinate; and determine the probability density distribution function corresponding to each fracture phase on the attribute parameter crossplot, which can be used as the second fracture density calculation model.

[0080] In some embodiments, the step of combining the first attribute data body to obtain the second attribute data body is specifically: performing two or three combination mathematical calculation on the first attribute data body to obtain the second attribute data body, so as to achieve the optimization purpose of the first attribute data body. Wherein, the two or three combination mathematical calculation on the first attribute data body is performed by using the mathematical function relationship, and the second attribute data body obtained by the calculation is an attribute data with better correlation with the fracture density curve on the well. In this embodiment, the mathematical function relationship is as follows:

[0081] B i = f2(E1, E2, …, En) (6) w

[0082] In formula (6), B i is the i-th attribute data in the second attribute data body of a certain fracture phase, E w is the first attribute data body of the w-th fracture phase (w < n, n is the number of the first attribute data body), and f2 is the mathematical function relationship (the correlation coefficient value with the fracture density curve on the well is larger after the calculation by the mathematical function relationship). Generally, the mathematical function relationship includes linear relationship, polynomial relationship, exponential relationship and logarithmic relationship, etc.

[0083] In step 302, the first fracture density data body and the second attribute data body of each fracture phase (high-angle fracture and low-angle fracture) are input into the related second fracture density calculation model for calculation to obtain the second fracture density data body related to each fracture phase.

[0084] Generally, the main operation is to convert and calculate the related first fracture density data body and second attribute data body by using the second fracture density calculation model (i.e. the probability density distribution function) corresponding to each fracture phase, so as to obtain the second fracture density data body corresponding to each fracture phase, which is expressed as:

[0085] P j = f3(B1, B2, …, Bn) (7) y

[0086] In formula (7), P j ​​B is the second fracture density data volume of the jth fracture phase with respect to high and low angles y B1 is the first fracture density data volume of the yth fracture phase, and f3 is a function relationship between the first attribute data volume and the second attribute data volume of the fracture phase (consistent with f1, that is, the best matching function determined by the least square method). In actual operation, the correlation between the statistical fracture density of the relevant fracture phase and the first fracture density data volume and the second attribute data volume on the probability density distribution function can be used as a constraint condition for Markov chain Monte Carlo + Bayesian probability simulation, so as to obtain the second fracture density data volume corresponding to the fracture phase through calculation.

[0087] In step 303, the second fracture density data volume is subjected to inverse normalization and weighted fusion processing to obtain a third fracture density data volume.

[0088] In some embodiments, for inverse normalization, the formula for inverse normalization used is mainly an inverse operation formula of the normalization processing of the fracture density on the well, and the calculation formula is as follows:

[0089]

[0090] In formula (8), y (x,y,t) is the fracture density data of a sampling point on the second fracture density data volume; y ( ′ x,y,t) is the fracture density data in the time domain corresponding to the inverse normalization x max = max{x p}, x min = min{x p}, n and m are positive integers, n and m are positive integers, and m > n ≥ 0, which is the same as formula (1).

[0091] In some embodiments, for weighted fusion processing, the calculation formula is as follows:

[0092]

[0093] In formula (9), Y is the third fracture density data volume after weighted fusion calculation processing, C i is the second fracture density data volume of the ith fracture phase; k i is a weighted coefficient of the second fracture density data volume, and k i ≥ 0, and n is the total number of the entire fracture phase (including high and low angle fractures) for weighting. Generally, k i = 1 / n.

[0094] Preferably, the weighting coefficients of each fracture phase can be designed by using the high and low angle fracture density curves of the relevant well. The main operation is to set the test well, design a series of weighting coefficients of each fracture phase respectively, and calculate the fracture density of the well; when the sum of the error with the fracture density curve of the well is the smallest, the corresponding weighting coefficient is the best weighting coefficient participating in the weighted fusion processing. In addition, when the number of test wells is generally greater than or equal to the number of fracture phases, the relevant equation set can be designed to calculate the best weighting coefficient of each fracture phase, and the best weighting coefficient of each fracture phase obtained by calculation is participated in the weighted fusion processing.

[0095] One example:

[0096] According to the technical process of the present application as Figure 1 shown, the working steps are formulated, and an example is to calculate the fracture density data volume of the tight sandstone section of a certain three-dimensional work area in the Sichuan Basin, so as to overall evaluate the fracture development in the channel sand body reservoir, and to optimize the relevant favorable area to implement the horizontal well design and drilling.

[0097] In step 100, the fracture density values in each well in the exploration area are calculated by well logging resampling to obtain fracture density curves of each well, and the fracture density curves are divided into high and low angle fracture density curves according to the occurrence of the fractures. In actual operation, the sampling interval of the fracture density curve is artificially defined as 2m; the well-seismic calibration is performed by using the logging data of each well to determine the position of the sandstone target layer in the seismic data, and the top and bottom interface horizons are interpreted to obtain the relevant time-depth relationship table; for fracture prediction, the P-wave anisotropy data volume, instantaneous amplitude, curvature volume, and coherence volume, maximum likelihood attribute data volume and chaos attribute data volume, instantaneous frequency volume, coherent data volume calculated based on 35hz frequency division data volume, curvature data volume calculated based on 35hz frequency division data volume, and ant attribute data volume based on related fracture prediction data are calculated and extracted by using conventional commercial software-VVA software and FRS software, totaling 18 seismic inversion and attribute data volumes. The interpreted horizons are used to extract the inversion and attribute data of the target layer section of each seismic inversion and attribute data volume, and the high and low angle fracture density values in the target layer section of each well are normalized and calculated to obtain the relevant normalization function calculation formula; then the corresponding normalization function calculation formula is used to normalize the fracture density curve and the seismic inversion and attribute data volume, and all kinds of data are normalized to the value range of (5, 100), wherein the normalization function calculation formula is y=4.2x+8, and x is the measured fracture density data; y is the normalized fracture density data. According to the fracture density data of the target layer section of the known well point after normalization and the attribute data on the seismic inversion and attribute data volume well point, the correlation coefficient is calculated, the data volume with higher correlation coefficient is optimized, the correlation coefficient between these data volumes is calculated, and the seismic inversion and attribute data volume with lower correlation coefficient is optimized for the next step. In the example, according to the operation procedure of step 1, 12 seismic inversion and attribute data volumes including P-wave anisotropy inversion data volume, curvature volume, and ant volume, maximum likelihood attribute volume, etc. are optimized, and the normalized fracture density data on the well are calculated for the next step.

[0098] In step 200, according to the histogram statistics results of the high and low angle fracture density data of the sandstone target interval of each well, both are set to three threshold values, thereby obtaining four fracture phases of high and low angle fractures respectively, a total of eight fracture phases, and determining the related first attribute data body and the first fracture density calculation model, and calculating eight first fracture density data bodies respectively. The specific operation is as follows: first, extract 12 attribute curves such as P-wave anisotropy data body and curvature body attribute at the well point, convert the corresponding normalized P-wave anisotropy data body, curvature body and 12 attribute curves from time domain to depth domain according to the time-depth relationship of each well, and resample the attribute curves in the depth domain according to the sampling interval of the fracture density curve to obtain the attribute curves after depth domain conversion and resampling processing. In the example, the sampling interval of various attribute curves in the depth domain is 2m, which is consistent with the sampling interval of the set fracture density curve; then, according to the related first attribute optimization and first fracture density calculation model determination method in step 202, the first attribute data body and the first fracture density calculation model of the eight fracture phases are respectively obtained from the 12 attribute data bodies; after the related first attribute data body and the first fracture density calculation model are used for fracture density calculation, the first fracture density data body of the eight fracture phases is obtained. In actual operation, for high angle fractures, six seismic attributes of the 12 optimized attributes are selected for calculation in the first fracture phase; four seismic attributes are selected for calculation in the second fracture phase, seven seismic attributes are selected for calculation in the third fracture phase, and five seismic attributes are selected for calculation in the fourth fracture phase.

[0099] In step 300, according to the first fracture density data body of each fracture phase of the high and low angle fractures, the second attribute data body and the fracture density data on the well, the second fracture density calculation model is established and the second fracture density data body is calculated, and then the second fracture density data body is inversely normalized and weighted and fused to obtain a third fracture density data body related to the fracture density calculation. In actual operation, the measured fracture density data of the related fracture phase on the well and the corresponding second attribute data value optimized on the basis of the first attribute data body, the first fracture density data body are used to establish a two-dimensional coordinate system, the probability density distribution function corresponding to each fracture phase is calculated by using the Bayesian classification method, so that the second fracture density calculation model is obtained; and the second fracture density data body of each fracture phase is calculated by using the second fracture density calculation model and the second attribute data body and the first fracture density data body, so that eight second fracture density data bodies are obtained; the inverse calculation (inverse normalization processing) of the fracture density normalization calculation is performed on the eight second fracture density data bodies, so that eight inverse normalized second fracture density data bodies are obtained; the weighted fusion processing is performed on the eight second fracture density data bodies by setting the weighting coefficients, so that a third fracture density data body related to the fracture density calculation is obtained. In the example, the weighting coefficient data values of the eight fracture density data bodies are all set to 0.125. For the optimization of the second attribute data body, in actual operation, after the optimization combination calculation of the six first attribute data bodies selected for the first fracture phase of the high angle fracture, three optimized second attribute data bodies are obtained, which participate in the establishment of the second fracture density model and are used as the input of the model calculation together with the first fracture density data body.

[0100] The fracture density calculation results obtained by using the present application are verified by the sandstone cores and FMI logging data of subsequent drilling, and the coincidence rate can reach more than 84.5%, which is better than the results obtained by using the conventional fracture prediction technology (such as curvature technology, coherence technology and P-wave anisotropy technology, etc.). It is proved that the present application is effective for the fracture density calculation of different scales, and is worth using in the fracture prediction and evaluation of some tight sandstone, marine shale and carbonate rock, etc.

[0101] The above only describes the preferred embodiments of the present application and is not used to limit the present application. For those skilled in the art, the present application can have various modifications and changes. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application shall be included in the protection scope of the present application.

Claims

1. A method of computing a fracture density data volume, characterized by, The method comprises the following steps: Step 100, establishing high and low angle fracture density curves in a known well, calculating related seismic attribute data bodies and optimizing the seismic attribute data bodies to obtain optimized attribute data bodies, and then normalizing the high and low angle fracture density curves in the well and the optimized attribute data bodies to obtain normalized high and low angle fracture density curves in the well and the optimized attribute data bodies; Step 200, setting related threshold values for the normalized high and low angle fracture density curves in the well respectively and establishing fracture facies, optimizing the optimized attribute data bodies based on the fracture facies to obtain first attribute data bodies, and establishing related first fracture density calculation models to calculate first fracture density data bodies about each fracture facies of high and low angle fractures; Step 300, establishing second fracture density calculation models for the first fracture density data bodies, the first attribute data bodies and the fracture density curves and performing fracture density calculation to obtain second fracture density data bodies of each fracture facies, and then performing reverse normalization and weighted fusion processing on the second fracture density data bodies to obtain third fracture density data bodies.

2. The method of computing a fracture density data volume of claim 1, wherein, Step 100 comprises the following sub-steps: Step 101, preparing geological data and logging data, and forming high and low angle fracture density curves according to the geological data and the logging data respectively; Step 102, preparing seismic data, performing inversion and attribute extraction on the seismic data to obtain seismic inversion and attribute data bodies, extracting seismic inversion and attribute curve data of a target layer at each well point from the seismic inversion and attribute data bodies, and performing time-depth conversion and resampling calculation on the seismic inversion and attribute curve data to form fracture density-seismic attribute corresponding relationship between the seismic inversion and attribute curve data and data on the corresponding fracture density curves in the depth domain; Step 103, performing normalization and correlation coefficient calculation on the seismic inversion and attribute curve data of the target layer at each well point and data on the high and low angle fracture density curves of the target layer at the corresponding well point respectively, and taking M and L seismic inversion and attribute data bodies with higher correlation coefficients as optimized attribute data bodies; wherein, the M seismic inversion and attribute data bodies correspond to data on the low angle fracture density curve, and the L seismic inversion and attribute data bodies correspond to data on the high angle fracture density curve; Step 104, performing normalization calculation on the data on the high and low angle fracture density curves of the target layer at the well point and the M and L seismic inversion and attribute data bodies to obtain respective normalized function formulas; Step 105, performing normalization processing calculation on the M and L seismic inversion and attribute data bodies according to the normalized function formulas to obtain normalized M and L seismic inversion and attribute data bodies.

3. The method of computing a fracture density data volume of claim 2, wherein, In step 103, the M and L seismic inversion and attribute data bodies with higher correlation coefficients refer to: firstly, performing first correlation coefficient calculation, then performing second correlation coefficient calculation between the M and L seismic inversion and attribute data bodies with higher correlation coefficients obtained by the first correlation coefficient calculation, and thirdly, again optimizing the M and L seismic inversion and attribute data bodies with lower correlation coefficients obtained by the second correlation coefficient calculation.

4. The method of computing a fracture density data volume of claim 1, wherein, In step 200, the fracture facies are established by: extracting data on the normalized high-angle and low-angle fracture density curves respectively, and establishing high-angle and low-angle fracture density data sets of the target interval; setting a series of related threshold values for the high-angle and low-angle fracture density data sets of the target interval respectively, and dividing different fracture density ranges, so as to divide into different fracture density ranges; the fracture density data in different fracture density ranges belong to different fracture facies, so as to obtain each fracture facies with respect to high-angle and low-angle.

5. The method of computing a fracture density data volume of claim 1, wherein, In step 200, the first attribute data body is obtained by optimizing the preferred attribute data body based on the fracture facies, including: According to the fracture density data, the fracture density sample value of each fracture is determined, and the preferred attribute data body corresponding to the fracture density sample value of each fracture facies is determined: the well-seismic calibration is performed according to the acoustic and density logging curves and the post-stack seismic data, and the time-depth relationship of each well is determined; the curve data of the preferred attribute data body is converted from the time domain to the depth domain by using the time-depth relationship, and the curve data of the preferred attribute data body in the depth domain is resampled and calculated to make it consistent with the sampling interval of the fracture density curve, and a one-to-one correspondence between data values is formed, and through the correspondence, the preferred attribute data body corresponding to the fracture density sample value of the fracture facies is determined; Through the cross-plot analysis, the attribute curve of the preferred attribute data body with the highest correlation coefficient with the fracture density curve of each well in the fracture facies is taken as the first attribute curve; then the first attribute curve and other preferred attribute curves are sequentially taken as curve pairs, and the best matching function determined by the least square method is used to calculate the fitting fracture density curve, and the correlation coefficient between the fitting fracture density curve and the fracture density curve of the fracture facies is calculated, and the highest correlation coefficient is the second attribute curve; then the third attribute curve, the fourth attribute curve, …, and the Nth attribute curve are continuously searched for in the same way, and each attribute curve is the first attribute data body.

6. The method of computing a fracture density data volume of claim 5, wherein, In the cross-plot analysis, if the data points of the preferred attribute data body corresponding to the fracture density sample value of the fracture facies do not meet the correlation requirements, the fracture density sample value of the fracture facies is expanded to multiple fracture density values within a certain range.

7. The method of computing a fracture density data volume of claim 5, wherein, In step 200, the first fracture density calculation model is established by determining the number of preferred attribute data bodies by using the method of blind well verification, i.e. using each well as a verification well in turn, calculating the fitting fracture density curve of the verification well by using the number of attribute curves determined by other wells and the best matching function, and calculating the error between the fitting fracture density curve and the fracture density curve of the fracture facies; after obtaining the error when each well is used as a verification well, the attribute curve combination and the best matching function with the lowest average verification error are the first fracture density calculation model of the fracture facies.

8. The method of computing a fracture density data volume of claim 1, wherein, In step 300, the second fracture density calculation model is established by using the first fracture density data body of each fracture facies, the second attribute data body and the corresponding fracture density curve, establishing the corresponding attribute parameter crossplot, and establishing the function formula related to the fracture density calculation according to the result of the attribute parameter crossplot, so as to establish the second fracture density calculation model; wherein the second attribute data body is obtained by combination calculation of the first attribute data body.

9. The method of computing a fracture density data volume of claim 8, wherein, The function formula is a probability density distribution function.

10. The method of computing a fracture density data volume of claim 8, wherein, The second attribute data body is obtained by mathematical calculation of the first attribute data body in combination of two or three. The second attribute data body is obtained by mathematical calculation of the first attribute data body in combination of two or three.

Images