Method for identifying shale fractures based on conventional well logging data

Abstract

The application provides a method for identifying mud shale fractures based on conventional logging data, comprising the following steps: step 1, obtaining acoustic travel time logging curves and natural gamma logging curves of a mud shale section; step 2, performing discrete wavelet transform and box dimension calculation on the acoustic travel time curves and the natural gamma curves to obtain wavelet coefficients and box dimensions; step 3, calculating the difference between the wavelet coefficients of the acoustic travel time curves and the wavelet coefficients of the natural gamma curves to obtain wavelet coefficient differences, and calculating the difference between the box dimensions of the acoustic travel time curves and the box dimensions of the natural gamma curves to obtain box dimension differences; step 4, averaging the wavelet coefficient differences and the fractal dimension differences to construct a fracture development response index; and step 5, identifying fractures in the mud shale section according to the relatively high value of the fracture development response index. The method for identifying mud shale fractures based on conventional logging data can accurately and reliably identify mud shale fracture development sections.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oilfield logging technology, and particularly relates to a method for identifying mud shale fractures based on conventional logging data. BACKGROUND

[0002] Mud shale reservoirs have the characteristics of low porosity and low permeability, and are generally good hydrocarbon source rocks and cap rocks. Fractures are discontinuous surfaces formed by the breaking of rocks under stress, and often develop in ultra-low permeability reservoirs. The reservoir space of mud shale reservoirs is usually composed of pores and fractures. Fractures are both migration channels and reservoir spaces, and are one of the main factors controlling the productivity of mud shale reservoirs.

[0003] Core and imaging logging data are important means for identifying fractures, and have the characteristics of intuitiveness and high identification accuracy. However, core and imaging logging data are limited in quantity and high in cost. Therefore, it is of great significance to identify fractures by using conventional logging data. The existing method for identifying mud shale fractures by using conventional logging data is to calibrate conventional logging data by using limited drilling core data and imaging logging data, such as logging response characteristic method, three porosity ratio method, and resistivity invasion correction difference ratio method. Conventional logging data are jointly affected by factors such as lithology, fluid, and fractures, and logging curves have a certain response to fractures, but the response characteristics are not obvious. At the same time, the response characteristics of conventional logging curves have multiple solutions, and fractures are not the only factor causing the response of logging curves. It is difficult to directly identify fractures by using conventional logging data.

[0004] When fractures appear in mud shale, transient characteristics appear on logging curves, which are reflected in the changes and details of logging signals. Therefore, we change our thinking and identify mud shale fractures by extracting and amplifying the logging response of fractures and eliminating the influence of non-fracture factors. We obtain the information of formation heterogeneity in logging curves, eliminate the vertical heterogeneity caused by the cyclical change of formation lithology, and highlight the vertical heterogeneity of the formation caused by fractures.

[0005] Based on the theoretical basis of wavelet analysis, we combine logging with multi-scale analysis theory to extract the high-frequency attributes of logging signals, amplify the oscillation characteristics of fracture information through high-frequency attributes, and identify mud shale fractures. Fractal dimension is a quantitative characterization parameter for describing the complexity of an object. When a reservoir section develops fractures, the complexity and abnormality of logging data are increased, which is reflected in the increase of fractal dimension.

[0006] In order to identify the fracture development section of the mud shale section by wavelet transform and fractal dimension analysis of conventional logging curves, it is necessary to combine core data and imaging logging data in the research area for calibration and correction to achieve the best effect.

[0007] In the Chinese patent application with application number CN202010010511.7, a method and device for identifying and characterizing fractures in shale reservoirs are disclosed. The method includes: based on the distribution differences of acoustic travel time and resistivity in fracture development segments and non-fracture development segments, respectively, a crossplot is obtained by weighing the influence of lithology, using a weighted algorithm to fuse natural gamma ray logging data and acoustic travel time logging data to obtain a lithology-property fusion parameter for fracture identification, and combining normalized resistivity to obtain a preliminary response chart for identifying fractures; using the preliminary response chart for identifying fractures to separate data points once, and separating data points in the overlapping part of the preliminary response chart for identifying fractures twice to obtain the final fracture identification result. The embodiments of the invention use conventional logging data to identify and characterize fractures, greatly improving the efficiency and reducing the cost of logging interpretation.

[0008] In the Chinese patent application with application number CN201710487285.X, a method for identifying fractures in tight limestone reservoirs is disclosed, which belongs to the field of well logging technology and solves the technical problems of low identification accuracy and easy misjudgment in the prior art. The invention uses acoustic porosity, natural gamma value and deep lateral resistivity relative value to improve the fracture comprehensive identification model and construct a new fracture comprehensive identification index, which is simple in form and high in identification accuracy, avoiding the drawbacks of easy misjudgment and low identification accuracy in traditional methods, and establishing a fracture grading identification standard according to the size of the fracture comprehensive identification index FRA to guide the exploration and development of such oil and gas reservoirs.

[0009] In the Chinese patent application with application number CN201911241396.8, a method for identifying fractures and quantitatively calculating porosity in tight reservoirs is disclosed, which includes: cutting the acoustic and density logging curves based on acoustic and density, cutting the extreme points of the curves, and taking the midpoint as the cutting point when a straight line segment is encountered; calculating the correlation coefficient of the cut curve segments according to the principle of correlation coefficient; calculating the AC-DEN correlation coefficient to identify the fracture development segment and the non-fracture development segment in the tight reservoir; evaluating the fractures, establishing an equation set using logging data to minimize the objective function, and obtaining the comprehensive skeleton of the reservoir acoustic travel time and density; calculating the shale content of the tight reservoir; and establishing a fracture porosity calculation model using acoustic travel time and density logging data. The fracture porosity results calculated by the invention have improved fracture identification accuracy compared with the results calculated by FMI, and the use of conventional logging data to identify and evaluate fractures greatly improves the efficiency of logging interpretation.

[0010] The above prior art has significant differences from the present invention and fails to solve the technical problems we want to solve, so we have invented a new method for identifying shale fractures based on conventional logging data. SUMMARY

[0011] The purpose of the present application is to provide a method for accurately and reliably identifying shale fractures based on conventional logging data.

[0012] The purpose of the present application can be achieved by the following technical measures: a method for identifying shale fractures based on conventional logging data, comprising:

[0013] Step 1: obtaining acoustic travel time logging curves and natural gamma logging curves of shale sections;

[0014] Step 2: performing discrete wavelet transform and box dimension calculation on the acoustic travel time curves and the natural gamma curves to obtain wavelet coefficients and box dimensions;

[0015] Step 3: calculating the difference between the wavelet coefficients of the acoustic travel time curves and the wavelet coefficients of the natural gamma curves to obtain the wavelet coefficient difference, and calculating the difference between the box dimensions of the acoustic travel time curves and the box dimensions of the natural gamma curves to obtain the box dimension difference;

[0016] Step 4: averaging the wavelet coefficient difference and the fractal dimension difference to construct a fracture development response index;

[0017] Step 5: identifying fractures in the shale section according to the relatively high value of the fracture development response index.

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

[0019] In step 1, analyze the logging curves of the shale section, and obtain the acoustic travel time logging curves AC and the natural gamma logging curves GR for shale fracture identification.

[0020] In step 1, according to the theory of well logging, the acoustic travel time AC curve can not only reflect the vertical heterogeneity caused by lithological cycle changes, but also reflect the local heterogeneity caused by secondary pores and fractures; the natural gamma GR curve mainly reflects the vertical heterogeneity caused by lithological cycle changes.

[0021] In step 2, the acoustic travel time curves and the natural gamma curves are subjected to discrete wavelet transform, and the logging curves are decomposed into approximate wavelet coefficients and detailed wavelet coefficients. The discrete wavelet transform is defined as the sum of the wavelet function ψ(t) after shifting and stretching of all original signals f(t) multiplied by the independent variable t. The calculation formula of the wavelet transform coefficient is as follows:

[0022]

[0023] In the formula, W f(a,b) are wavelet transform coefficients, f(t) is the original signal, ψ(t) is the wavelet function, and a and b are scale factor and shift factor respectively.

[0024] In step 2, the AC curve and the GR curve are subjected to three-layer wavelet decomposition by using a db 11 wavelet function, and the second layer of decomposition is analyzed in detail according to the wavelet coefficients cD2 in accordance with the fracture development condition obtained from the imaging logging of the shale section. For the convenience of subsequent analysis, the cD2 is normalized in the shale section, and the upper envelope of the normalized wavelet coefficients is obtained by Hilbert transform, and finally the wavelet coefficients of the acoustic traveltime curve and the natural gamma curve are obtained.

[0025] In step 2, the wavelet coefficients are normalized.

[0026]

[0027] In the formula, cD2N is the normalized wavelet coefficient, cD2 max is the maximum value of the wavelet coefficient of the shale section.

[0028] In step 2, the fractal dimension of the acoustic traveltime curve and the natural gamma curve is analyzed to obtain the box dimension of the acoustic traveltime curve and the natural gamma curve, and the box dimension is defined as:

[0029]

[0030] In the formula, r is the length of the small box, N(r) is the minimum number of boxes covering the logging curve, and D is the box dimension.

[0031] In step 2, the box dimension of the AC curve and the GR curve is obtained by using the sliding window method, the step is 1 logging data point, and the corresponding depth of the middle logging data point in the window is taken as the depth value of the box dimension. According to the fracture development condition obtained from the imaging logging of the shale section, the box dimension of the window length of 11 logging data points is mainly analyzed, and finally the box dimension of the acoustic traveltime curve and the natural gamma curve is obtained.

[0032] In step 2, the wavelet coefficients and the fractal dimension of the AC curve reflect the vertical heterogeneity caused by the local change of rock physical properties superimposed on the basis of lithology change; although the wavelet coefficients and the fractal dimension of the GR curve are not obvious in detecting the local physical property change of the stratum rock, they mainly reflect the vertical heterogeneity caused by the lithology cycle change of the stratum.

[0033] In step 3, the difference between the wavelet coefficient of the acoustic traveltime curve and the wavelet coefficient of the natural gamma curve is obtained, and the wavelet coefficient difference is obtained, and the formula is as follows:

[0034] cD2M = |c2AB-c2GB|

[0035] In the formula, cD2M is the wavelet coefficient difference, c2AB is the wavelet coefficient of the acoustic travel time curve, and c2GB is the wavelet coefficient of the natural gamma curve.

[0036] In step 3, the difference between the box dimension of the acoustic travel time curve and the box dimension of the natural gamma curve is obtained to obtain the box dimension difference, and the formula is as follows:

[0037] BFDM = |BFDA-BFDG|

[0038] In the formula, BFDM is the box dimension difference, BFDA is the box dimension of the acoustic travel time curve, and BFDG is the box dimension of the natural gamma curve.

[0039] In step 3, the difference between the wavelet coefficient of the acoustic travel time curve and the wavelet coefficient of the natural gamma curve, and the difference between the box dimension of the acoustic travel time curve and the box dimension of the natural gamma curve are used to eliminate the vertical heterogeneity caused by the lithological cycle change of the stratum and highlight the vertical heterogeneity of the stratum caused by the fracture.

[0040] In step 4, the response of the wavelet coefficient difference and the fractal dimension difference to the shale fracture is comprehensively considered, the wavelet coefficient difference and the fractal dimension difference are averaged, and a fracture development response index is constructed.

[0041]

[0042] In the formula, FDR is the fracture development response index, cD2M is the wavelet coefficient difference, and BFDM is the box dimension difference.

[0043] The method for identifying shale fractures based on conventional logging data in the application starts from extracting the logging response of enlarged fractures and eliminating the influence of non-fracture factors, obtains the wavelet coefficient and the box dimension of the acoustic travel time curve and the natural gamma curve to obtain the information of the stratum heterogeneity in the logging curve, obtains the wavelet coefficient difference and the box dimension difference of the acoustic travel time curve and the natural gamma curve to highlight the vertical heterogeneity of the stratum caused by the fracture, comprehensively considers the wavelet coefficient difference and the box dimension difference to construct the fracture development response index, and finally identifies the shale fracture according to the relatively high value of the fracture development response index, so that the shale fracture development can be accurately and reliably identified through the conventional logging data.

[0044] This invention, based on conventional well logging data, uses wavelet analysis and fractal theory to obtain formation heterogeneity information from conventional well logging curves. It then innovatively utilizes the difference in wavelet coefficients and box-counting dimension between the sonic transit time curve and the natural gamma curve to highlight the vertical heterogeneity of the formation caused by fractures, achieving the goal of fracture identification in shale with a novel approach and method. Core logging and imaging logging are important and accurate means of fracture identification, but only a few wells possess core and imaging logging data, making widespread application difficult. This invention, in the absence of core and imaging logging data, utilizes conventional well logging data available in most wells for fracture identification in shale, enabling widespread application. In shale reservoirs, fractures serve as both migration channels and reservoir spaces, and are one of the main factors controlling production capacity. Therefore, this invention is of great significance for the evaluation of shale reservoirs. Attached Figure Description

[0045] Figure 1 This is a flowchart of a specific embodiment of the method for identifying shale fractures based on conventional well logging data according to the present invention;

[0046] Figure 2 This is a diagram showing the logging curves and fracture identification results of the 3160-3255m mudstone and shale section of well A in this embodiment of the invention.

[0047] Figure 3 This is a diagram showing the logging curves and fracture identification results of the 3615-3675m mudstone and shale section of well B in this embodiment of the invention. Detailed Implementation

[0048] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0049] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0050] This invention relates to a method for identifying fractures in shale based on conventional well logging data. The following are several specific embodiments of the application of this invention.

[0051] Example 1

[0052] like Figure 1 As shown, Figure 1A flow chart of the method for identifying mud shale fractures based on conventional logging data of the present application. The method for identifying mud shale fractures based on conventional logging data comprises:

[0053] 1) analyzing the logging curves of the mud shale section, obtaining the acoustic travel time logging curve (AC) and the natural gamma logging curve (GR) for mud shale fracture identification;

[0054] 2) performing discrete wavelet transform on the acoustic travel time curve and the natural gamma curve, decomposing the logging curves into approximate wavelet coefficients and detailed wavelet coefficients. The discrete wavelet transform is defined as the sum of the wavelet function ψ(t) multiplied by all original signals f(t) after shifting and stretching the independent variable t, and the formula is as follows:

[0055]

[0056] In the formula, W f (a, b) is the wavelet transform coefficient, f(t) is the original signal, ψ(t) is the wavelet function, and a and b are the scale factor and the displacement factor, respectively;

[0057] The AC curve and the GR curve are decomposed into three layers using the db 11 wavelet function, and the detailed wavelet coefficients cD2 of the second layer decomposition are mainly analyzed. For the convenience of subsequent analysis, the cD2 is normalized in the mud shale section, and the upper envelope line of the normalized wavelet coefficient is obtained by Hilbert transform, and finally the wavelet coefficients of the acoustic travel time curve and the natural gamma curve are obtained.

[0058] Wavelet coefficient normalization processing:

[0059]

[0060] In the formula, cD2N is the normalized wavelet coefficient, cD2 max is the maximum value of the wavelet coefficient of the mud shale section.

[0061] The fractal dimension of the acoustic travel time curve and the natural gamma curve is analyzed to obtain the box dimension of the acoustic travel time curve and the natural gamma curve. The box dimension is defined as:

[0062]

[0063] In the formula, r is the length of the small box, N(r) is the minimum number of boxes covering the logging curve, and D is the box dimension;

[0064] The box dimension of the AC curve and the GR curve is obtained by using the sliding window method, and the step length is 1 logging data point. The corresponding depth of the middle logging data point in the window is taken as the depth value of the box dimension. The box dimension of the window length of 11 logging data points is mainly analyzed, and finally the box dimension of the acoustic travel time curve and the natural gamma curve is obtained.

[0065] 3) Calculate the difference between the wavelet coefficients of the acoustic time curve and the wavelet coefficients of the natural gamma curve. To avoid negative values, take the absolute value of the difference to obtain the wavelet coefficient difference, as follows:

[0066] cD2M = |c2AB-c2GB|

[0067] where cD2M is the wavelet coefficient difference, c2AB is the wavelet coefficient of the acoustic time curve, and c2GB is the wavelet coefficient of the natural gamma curve.

[0068] Calculate the difference between the box dimension of the acoustic time curve and the box dimension of the natural gamma curve. To avoid negative values, take the absolute value of the difference to obtain the box dimension difference, as follows:

[0069] BFDM = |BFDA-BFDG|

[0070] where BFDM is the box dimension difference, BFDA is the box dimension of the acoustic time curve, and BFDG is the box dimension of the natural gamma curve.

[0071] 4) Integrate the responses of the wavelet coefficient difference and the fractal dimension difference to mud shale fractures, average the wavelet coefficient difference and the fractal dimension difference, and construct a fracture development response index.

[0072]

[0073] where FDR is the fracture development response index.

[0074] 5) Identify mud shale fractures based on the relative high value of the fracture development response index in the mud shale section.

[0075] In step 1), according to well logging theory, the acoustic time curve (AC) can not only reflect the vertical heterogeneity caused by lithological cycle changes, but also reflect the local heterogeneity caused by secondary pores and fractures. The natural gamma (GR) curve mainly reflects the vertical heterogeneity caused by lithological cycle changes.

[0076] In step 2), it can be considered that the wavelet coefficients and fractal dimensions of the AC curve reflect the vertical heterogeneity caused by local changes in rock properties superimposed on the basis of lithological changes. Although the wavelet coefficients and fractal dimensions of the GR curve are not obviously advantageous in detecting local changes in the physical properties of the formation, they mainly reflect the vertical heterogeneity caused by lithological cycle changes.

[0077] In step 3), the difference between the wavelet coefficients of the acoustic transit time curve and the wavelet coefficients of the natural gamma curve, and the difference between the box dimension of the acoustic transit time curve and the box dimension of the natural gamma curve are used to eliminate the vertical heterogeneity caused by the cyclic changes in strata lithology and highlight the vertical heterogeneity of strata caused by fractures.

[0078] Example 2

[0079] This embodiment describes a method for identifying shale fractures based on conventional well logging data. Taking the shale reservoir of Well A in an oilfield as an example, the flowchart is shown below. Figure 1 This includes the following steps:

[0080] 1) Analyze the logging curves of the shale section and obtain the sonic transit time (AC) and natural gamma ray (GR) logging curves for shale fracture identification.

[0081] 2) The db11 wavelet function was used to perform three-level wavelet decomposition on the AC and GR curves. The detailed wavelet coefficients cD2 of the second level decomposition were mainly analyzed. To facilitate subsequent analysis, cD2 was normalized in the mudstone and shale section, and the upper envelope of the normalized wavelet coefficients was obtained by Hilbert transform. Finally, the wavelet coefficients of the acoustic time difference curve and the natural gamma curve were obtained.

[0082] The sliding window method was used to calculate the box dimension of AC and GR curves, with a step size of 1 logging data point. The depth corresponding to the logging data point in the middle of the window was taken as the depth value of the box dimension. The box dimension of the window with a length of 11 logging data points was mainly analyzed, and the box dimension of the sonic transit time curve and natural gamma curve were finally obtained.

[0083] 3) Calculate the difference between the wavelet coefficients of the acoustic time difference curve and the wavelet coefficients of the natural gamma curve, and calculate the difference between the box dimension of the acoustic time difference curve and the box dimension of the natural gamma curve.

[0084] 4) The response of shale fractures to the combined wavelet coefficient difference and fractal dimension difference is analyzed by averaging the wavelet coefficient difference and fractal dimension difference to construct a fracture development response index.

[0085] 5) Based on the relatively high value of the fracture development response index in the mudstone and shale section, fractures in the mudstone and shale are identified.

[0086] Fracture orientation and parameters were derived from imaging logging data to verify the effectiveness of this method for identifying fractures in shale. (See...) Figure 2 . Figure 2 The first channel is the lithology curve, the second channel is the depth channel, the third channel is the resistivity curve, the fourth channel is the porosity curve, the fifth channel is the fracture orientation obtained from imaging logging, the sixth channel is the fracture parameters obtained from imaging logging, the seventh and eighth channels are the wavelet coefficient difference and box dimension difference, respectively, and the ninth channel is the fracture development response index.

[0087] Well A is located in the lower part of the sag, the target formation has large and stable sedimentary thickness, shale is developed, and the abundance of organic matter is high. The 3160-3255m interval of Well A is a shale section. In this interval, GR is about 60gAPI, resistivity is about 10 ohm-m, density is about 2.55g / cm 3 , and a large number of fractures are developed with a high angle of about 60 degrees. Figure 2 Four fracture development intervals are obtained from the imaging logging, and the fracture development response index at the corresponding depth shows a relatively high value. The fracture development response index can be used to identify shale fractures.

[0088] Example 3

[0089] The method for identifying shale fractures based on conventional logging data in this example takes the shale reservoir of Well B in an oilfield as an example, and the flow chart is shown in Figure 1 , which includes the following steps:

[0090] 1) Analyze the logging curves of the shale section to obtain the acoustic travel time logging curve (AC) and the natural gamma logging curve (GR) for shale fracture identification.

[0091] 2) Use the db 11 wavelet function to perform three-layer wavelet decomposition on the AC curve and the GR curve, mainly analyze the detailed wavelet coefficients cD2 of the second layer decomposition, normalize cD2 in the shale section for subsequent analysis, and obtain the upper envelope of the normalized wavelet coefficients by Hilbert transform, and finally obtain the wavelet coefficients of the acoustic travel time curve and the natural gamma curve.

[0092] Use the sliding window method to calculate the box dimension of the AC curve and the GR curve, with a step of 1 logging data point, and take the corresponding depth of the middle logging data point in the window as the depth value of the box dimension, mainly analyze the box dimension with a window length of 11 logging data points, and finally obtain the box dimension of the acoustic travel time curve and the natural gamma curve.

[0093] 3) Calculate the difference between the wavelet coefficients of the acoustic travel time curve and the wavelet coefficients of the natural gamma curve, and calculate the difference between the box dimension of the acoustic travel time curve and the box dimension of the natural gamma curve.

[0094] 4) Synthesize the response of the wavelet coefficient difference and the fractal dimension difference to the shale fracture, average the wavelet coefficient difference and the fractal dimension difference, and construct a fracture development response index.

[0095] 5) Identify the shale fracture according to the relatively high value of the fracture development response index in the shale section.

[0096] The fracture occurrence and fracture parameters obtained from the imaging logging data are used to verify the shale fracture identification effect of the method, see Figure 3 . Figure 3 The first track is a lithology curve, the second track is a depth track, the third track is a resistivity curve, the fourth track is a porosity curve, the fifth track is a fracture occurrence obtained from imaging logging, the sixth track is a fracture parameter obtained from imaging logging, the seventh track and the eighth track are respectively a wavelet coefficient difference value and a box dimension difference value, and the ninth track is a fracture development response index.

[0097] Different from the well A, the well B is located on the east slope of the depression zone and develops a large set of mudstone and oil shale. The 3615-3675m well section of the well B is a mudstone section, in which the GR is about 75gAPI, the resistivity is about 6 ohm meters, the density is about 2.43g / cm 3 , and the fracture develops with a low fracture angle of about 30 degrees. Figure 3 In the well B, three fracture development sections are obtained from the imaging logging, and the fracture development response index at the corresponding depth appears a relatively high value, so that the shale fracture can be identified through the fracture development response index.

[0098] The present application is applied in different wells, and the fracture identification effect is relatively obvious although the specific reservoir parameters and fracture parameters are different, so that the method for identifying the shale fracture based on the conventional logging data is feasible and effective.

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

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

Claims

1. A method for identifying fractures in shale based on conventional well log data, characterized in that, The method for identifying mud shale fractures based on conventional logging data comprises the following steps: Step 1, obtaining acoustic travel time logging curves and natural gamma logging curves of the mud shale section; Step 2, performing discrete wavelet transform and box dimension calculation on the acoustic travel time curves and the natural gamma curves to obtain wavelet coefficients and box dimensions; Step 3, calculating the difference between the wavelet coefficients of the acoustic travel time curves and the wavelet coefficients of the natural gamma curves to obtain wavelet coefficient difference, and calculating the difference between the box dimensions of the acoustic travel time curves and the box dimensions of the natural gamma curves to obtain box dimension difference; Step 4, averaging the wavelet coefficient difference and the box dimension difference to construct a fracture development response index; Step 5, identifying fractures in the mud shale section according to the relative high value of the fracture development response index.

2. The method for identifying fractures in shale based on conventional well log data according to claim 1, wherein, In step 1, the logging curves of the mud shale section are analyzed to obtain acoustic travel time logging curves AC and natural gamma logging curves GR for mud shale fracture identification.

3. The method for identifying fractures in shale based on conventional well log data according to claim 2, wherein, In step 1, according to the theory of well logging, the acoustic travel time AC curve can not only reflect the vertical heterogeneity caused by the lithological cycle change of the formation, but also reflect the local heterogeneity caused by the secondary pore and fracture; the natural gamma GR curve mainly reflects the vertical heterogeneity caused by the lithological cycle change of the formation.

4. The method for identifying fractures in shale based on conventional well log data of claim 1, wherein, In step 2, the discrete wavelet transform is performed on the acoustic travel time curves and the natural gamma curves to decompose the logging curves into approximate wavelet coefficients and detailed wavelet coefficients, and the discrete wavelet transform is defined as the sum of the wavelet function ψ(t) after the original signal f(t) is shifted and stretched by the independent variable t, and the calculation formula of the wavelet transform coefficient is as follows: ； wherein W f (a,b) are the wavelet transform coefficients, f(t) is the original signal, ψ(t) is the wavelet function, and a and b are the scale and shift factors, respectively.

5. The method for identifying fractures in shale based on conventional well log data according to claim 4, wherein, In step 2, the AC curve and the GR curve are decomposed into three layers by using the db11 wavelet function, and the detailed wavelet coefficients cD2 of the second layer decomposition are mainly analyzed according to the coincidence with the fracture development obtained by the imaging logging of the mud shale section, for the convenience of subsequent analysis, the cD2 is normalized in the mud shale section, and the upper envelope line of the normalized wavelet coefficient is obtained by Hilbert transform, and finally the wavelet coefficients of the acoustic travel time curve and the natural gamma curve are obtained.

6. The method for identifying fractures in shale based on conventional well log data according to claim 5, wherein, In step 2, the wavelet coefficient normalization processing: ； In the formula, cD2N represents the normalized wavelet coefficients, and cD2 max This represents the maximum value of the wavelet coefficients for the mudstone and shale section.

7. The method for identifying fractures in shale based on conventional log data of claim 1, wherein, In step 2, the box dimension analysis is performed on the acoustic travel time curves and the natural gamma curves to obtain the box dimensions of the acoustic travel time curves and the natural gamma curves, and the box dimension is defined as: ； In step 2, the box dimension of the AC curve and the GR curve is calculated by using the sliding window method, the step is 1 logging data point, the corresponding depth of the middle logging data point in the window is taken as the depth value of the box dimension, and the box dimension with a window length of 11 logging data points is mainly analyzed according to the coincidence with the fracture development obtained by the imaging logging of the mud shale section, and finally the box dimensions of the acoustic travel time curve and the natural gamma curve are obtained.

8. The method for identifying fractures in shale based on conventional well log data according to claim 7, wherein, In step 2, the wavelet coefficients and the box dimensions of the AC curve reflect the vertical heterogeneity caused by the superposition of the local rock physical property change on the lithological change; although the wavelet coefficients and the box dimensions of the GR curve have no obvious advantage in detecting the local physical property change of the formation rock, they mainly reflect the vertical heterogeneity caused by the lithological cycle change of the formation.

9. The method for identifying fractures in shale based on conventional log data of claim 1, wherein, ​ 10. The method for identifying fractures in shale based on conventional well log data of claim 1, wherein, In step 3, the difference between the wavelet coefficients of the acoustic traveltime curve and the wavelet coefficients of the natural gamma curve is calculated to obtain the wavelet coefficient difference, and the formula is as follows: ； In the formula, cD2M is the wavelet coefficient difference, c2AB is the wavelet coefficient of the acoustic traveltime curve, and c2GB is the wavelet coefficient of the natural gamma curve.

11. The method for identifying fractures in shale based on conventional well log data according to claim 1, wherein, In step 3, the difference between the box dimension of the acoustic traveltime curve and the box dimension of the natural gamma curve is calculated to obtain the box dimension difference, and the formula is as follows: ； In the formula, BFDM is the box dimension difference, BFDA is the box dimension of the acoustic traveltime curve, and BFDG is the box dimension of the natural gamma curve.

12. The method for identifying fractures in shale based on conventional well log data of claim 1, wherein, In step 3, the difference between the wavelet coefficients of the acoustic traveltime curve and the wavelet coefficients of the natural gamma curve is calculated to obtain the wavelet coefficient difference, and the formula is as follows:

13. The method for identifying fractures in shale based on conventional well log data of claim 1, wherein, In step 4, the response of the wavelet coefficient difference and the box dimension difference to the shale fracture is integrated, the wavelet coefficient difference and the box dimension difference are averaged, and the fracture development response index is constructed; ； In the formula, FDR is the fracture development response index, cD2M is the wavelet coefficient difference, and BFDM is the box dimension difference.

Images