Natural crack identification method based on bivariate variable scale fractal technology

Abstract

The invention discloses a natural fracture identification method based on a double-variable scale-variable fractal technology, and the method comprises the steps: selecting a target well logging curve of a target layer series, obtaining positive sequence well logging data, and reversing the depth of the positive sequence well logging data to obtain reverse sequence well logging data; through R / S analysis, the range R and the standard deviation S of the two groups of data at different starting points m and sliding interval lengths n are calculated, and then the Hurst index H and the fractal dimension D are obtained. And converting the fractal dimension of the inverted-order data into a coordinate system which is the same as the coordinate system of the positive-order data, and drawing a two-dimensional fractal dimension D contour map. And according to the distribution of the fractal dimension D, selecting a top interface of the fracture development layer section in the positive sequence data area, and selecting a bottom interface in the reverse sequence data area. And the number, the top interface depth, the bottom interface depth and the thickness of the fracture development layer sections are quantitatively obtained through the top interfaces and the bottom interfaces of the selected fracture development layer sections. According to the invention, the recognition precision of the fracture development layer section can be improved.

Description

Technical Field

[0001] This application relates to the technical field of natural fracture identification, and particularly to a method for natural fracture identification. Background Art

[0002] The comprehensive characterization of natural fractures in a reservoir is the basis for predicting and modeling natural fractures in the reservoir, and is also an important reference for optimizing sweet spots in reservoir engineering and formulating exploration and development plans. At present, the methods for natural fracture identification and characterization mainly include geological and petrological methods, logging methods, drilling and mud logging, and production dynamic data analysis methods. Among them, modern mathematical means such as the R / S variable scale analysis technology based on the finite difference method of conventional logging data are widely used in the identification of natural fractures. This is because the acquisition cost of conventional logging data is relatively low among various types of logging data, and various types of logging curves have different degrees of response characteristics to the natural fracture development intervals, which can be used to quantitatively identify the number, bottom interface depth, and top interface depth of natural fracture development intervals. The essence of the R / S variable scale analysis technology is to simplify the identification of natural fractures into the identification of abnormal segments in logging curves, where R and S are the range and standard deviation respectively, as shown in the following formulas:

[0003]

[0004] Where Z is the logging curve that is more sensitive to the response of natural fractures; N is the number of data points of the logging curve; u is the data point number for calculating the range and standard deviation each time within the range of 0 to N. After calculating the range R(N) and standard deviation S(N), a curve graph of lg[R(N) / S(N)] versus lg(N) is plotted, and the curve slope (i.e., the Hurst exponent H) and fractal dimension D are calculated as follows:

[0005]

[0006] D = 2 - H

[0007] Moreover, in the prior art, the decreasing segment of the curve slope H (i.e., the increasing segment of the fractal dimension D) is considered to be the potential natural fracture development interval. However, there are several limitations in the current method:

[0008] Limitation 1: The identification of fracture development intervals requires artificial determination based on the curve shape, which is a heavy workload and highly subjective. Although some scholars have proposed an improved R / S-FD method to achieve automatic identification of fracture development intervals. However, this improved method directly identifies the decreasing segment of the curve slope as the concave segment of the curve, and the accuracy of this determination is not very high.

[0009] Limitation 2: In the prior art, the decreasing segment of the slope, i.e., the Hurst exponent H, is simultaneously determined as the natural fracture development interval, and this determination method is not accurate either.

[0010] Limitation 3: The traditional R / S variable scale analysis technology always calculates lg[R(N) / S(N)] starting from the top of the target interval. Therefore, the response of the abnormal data segment deep in the target interval is significantly weakened, and the recognition ability of the top interface of the fracture-developed interval in the target interval is stronger than that of the bottom interface of the fracture-developed interval, which limits the recognition accuracy of the thickness of the fracture-developed interval.

[0011] Therefore, the determination of the fracture-developed interval by the above method is not accurate, and its determination accuracy is still difficult to meet the actual research needs. Summary of the Invention

[0012] Based on this, in view of the above technical problems, a natural fracture recognition method based on the bivariate variable scale fractal technology is provided to solve the problem of low accuracy in determining the fracture-developed interval in the prior art.

[0013] In a first aspect, a natural fracture recognition method based on the bivariate variable scale fractal technology, the method includes:

[0014] Select the target logging curve of the target formation, and obtain the logging data corresponding to the target logging curve, denoted as the first logging data; reverse the first logging data according to the depth to obtain the second logging data;

[0015] By introducing the bivariate variable scale fractal technology, calculate the range R and standard deviation S of the first logging data and the second logging data at different starting points m and sliding interval lengths n based on the R / S analysis;

[0016] Calculate the Hurst exponent H of the first logging data and the second logging data, and thus obtain the fractal dimensions of the first logging data and the second logging data respectively; the Hurst exponent H is obtained by calculating the slope of the curve of lg[R(m,n) / S(m,n)] with respect to n; denote the fractal dimension obtained from the first logging data as the first fractal dimension;

[0017] Convert the fractal dimension obtained from the second logging data into the same coordinate system as the fractal dimension obtained from the first logging data through the coordinate conversion formula, and denote the fractal dimension corresponding to the second test data after coordinate conversion as the second fractal dimension;

[0018] Draw a two-dimensional fractal dimension D contour map according to the first fractal dimension and the second fractal dimension;

[0019] Select the top interface of the natural fracture development interval in the first fractal dimension region of the two-dimensional fractal dimension D contour map, and select the bottom interface in the second fractal dimension region; the rules for selecting the top and bottom interfaces of the natural fracture development interval include: near the starting point m*, when the sliding interval length n is within a specific interval (p, q), if the fractal dimension on the left side of the starting point m* is greater than 1 and the right side is less than 1; and when n is not within the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the top interface of the natural fracture development interval; and if when the sliding window length n is within the specific interval (p, q), if the fractal dimension on the left side of the starting point m* is less than 1 and the right side is greater than 1; and when n is not within the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the bottom interface of the natural fracture development interval;

[0020] Based on the selected top and bottom interfaces of the natural fracture development interval, obtain the number of natural fracture development intervals, the depth of the top interface, the depth of the bottom interface, and the thickness.

[0021] In the above solution, optionally, the target logging curves for the target formation series include: calculating the relationship between the number of natural fractures per well and the actual depth based on core data or imaging logging data;

[0022] Calculate the natural fracture line density based on the relationship between the number of natural fractures per well and the actual depth, conduct a sensitivity analysis of different types of logging curves on the natural fracture line density, and select the logging curve with the largest sensitivity coefficient as the target logging curve.

[0023] In the above solution, further optionally, the sensitivity analysis of different types of logging curves on the natural fracture line density includes: using linear goodness-of-fit analysis, Pearson correlation analysis, or Spearman rank correlation analysis to conduct the sensitivity analysis of different types of logging curves on the natural fracture line density.

[0024] In the above solution, further optionally, the natural fracture line density is calculated based on the relationship between the number of natural fractures per well and the actual depth through the following formula:

[0025]

[0026] Wherein, MD represents the depth at which the natural fracture line density is to be calculated; FD represents the relationship curve between the number of natural fractures per well and the depth; FD(MD) represents the value of the natural fracture line density to be measured; CL represents the cumulative number curve of natural fractures; and W is the measurement window length, usually taken as a constant greater than 1.

[0027] In the above solution, optionally, the initial value of the sliding interval length n and its increment are selected from 1 / 50 to 1 / 20 of the total thickness of the target formation series.

[0028] In the above solution, optionally, the range R and standard deviation S of the first logging data and the second logging data at different starting points m and sliding interval lengths n are calculated by the following formula:

[0029]

[0030] Where R and S are the range and standard deviation respectively, Z is the target logging curve; N is the total number of data points of the target logging curve; u is the data point number for calculating the range and standard deviation each time within the range of 0 to N.

[0031] In the above solution, optionally, the Hurst exponent H of the first logging data and the second logging data is calculated, and the fractal dimensions of the first logging data and the second logging data are obtained respectively through the following formula:

[0032]

[0033] D(m,n) = 2 - H

[0034] Where H is the Hurst exponent and D(m,n) is the fractal dimension.

[0035] In the above solution, optionally, the selection rules for the top and bottom interfaces of the natural fracture development interval further include: taking the vertical line of the starting point m* as the top interface of the natural fracture development interval, denoted as the initial top interface. If the numerical region with a fractal dimension less than 1 with the top interface as the left boundary intersects with the horizontal line represented by the sliding window length n equal to 2*n0, then the initial top interface is further determined as the final top interface of the natural fracture development interval; n0 is the initial value of the sliding interval length;

[0036] Taking the vertical line of the starting point m* as the bottom interface of the natural fracture development interval, denoted as the initial bottom interface. If the fractal dimension of the numerical region with the initial interface as the right boundary less than 1 is equal to 1 at the horizontal line represented by the sliding window length n equal to N - 2*n0, then the initial bottom interface is further determined as the final bottom interface of the natural fracture development interval.

[0037] In a second aspect, a computer device includes a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements the steps of the natural fracture identification method based on the bivariate variable-scale fractal technology described in the first aspect above.

[0038] In a third aspect, a computer-readable storage medium stores a computer program which, when executed by a processor, implements the steps of the natural fracture identification method based on the bivariate variable-scale fractal technology described in the first aspect above.

[0039] The present application has at least the following beneficial effects:

[0040] The present application uses the slope of the logarithm lg[R(m,n) / S(m,n)] of the ratio of the range to the standard deviation with respect to n to characterize the fractal dimension D(m,n), which more intuitively reflects the slope of lg[R(m,n) / S(m,n)] with respect to the depth of the target layer. As a result, the characterized fractal dimension D(m,n) is related to the depth, and can more intuitively display the characteristics of the natural fracture development section. At the same time, the invention uses two variables, the initial position m of the sampling point and the sampling interval length n, to carry out the R / S variable-scale fractal analysis, which overcomes the drawback in the prior art that the recognition accuracy of the natural fracture development layers at different depths varies greatly due to always starting the R / S variable-scale fractal analysis from the top interface of the target layer.

[0041] Moreover, the present invention uses a method that combines the forward sequence (from shallow to deep in terms of depth) and the reverse sequence (from deep to shallow in terms of depth) R / S variable-scale fractal analysis of the conventional well logs of the target layer, which overcomes the drawback in the prior art that the recognition accuracies of the top and bottom interfaces of natural fractures vary greatly.

[0042] After the bivariate R / S variable-scale fractal technology proposed by the present invention draws the fractal dimension plane, if there exists an interval (p,q) at a certain starting point mm*, such that the fractal dimension D(m,n) < 1 only holds within n ∈ (p,q), and D(m,n) > 1 outside the interval, and the fractal dimension on the left side of the starting point m* is greater than 1 in a specific interval and less than 1 on the right side, then the vertical line represented by the starting point m* serves as the top interface of the natural fracture development section. If the fractal dimension on the left side of the starting point m* is less than 1 in a specific interval and greater than 1 on the right side, then the vertical line of this point m* serves as the bottom interface of the natural fracture development section. Judging the boundary surface in this way can improve the accuracy and reliability of the selection of the boundary line. Description of the Drawings

[0043] Figure 1 It is a comparison chart of the analysis results of the R / S variable-scale fractal technology of the prior art for 4 sets of ideal data;

[0044] Figure 2 It is a schematic flowchart of a natural fracture identification method based on the bivariate variable-scale fractal technology provided by an embodiment of the present application;

[0045] Figure 3 It is a relationship diagram between the acoustic travel time (AC) well log data and the depth of the target layer series in the study area provided by an embodiment of the present application;

[0046] Figure 4 This is the analysis result of the fracture development interval of a natural fracture identification method using a dual-variable variable-scale fractal technology according to an embodiment of the present application. Detailed implementation manners

[0047] In order to make the objectives, technical solutions and advantages of the present application clearer, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0048] In the description of the present application: Unless otherwise specified, "a plurality of" means two or more. The terms "first", "second", "third", etc. in the present application are intended to distinguish the objects being referred to and do not have special meanings in terms of technical connotations (for example, it should not be understood as emphasizing the importance or order, etc.). Expressions such as "including", "comprising", "having", etc. also mean "not limited to" (certain units, components, materials, steps, etc.).

[0049] In the prior art, it is considered that the section where the slope, i.e., the Hurst exponent H, decreases is used to study the natural fracture development interval, but whether this determination can represent it remains doubtful. According to the attached Figure 1 As shown, 4 data groups are set up to carry out numerical experiments. In the data group Figure 1 shown in Attachment A, which represents the rock density, it is generally considered that the low-density area can represent the natural fracture development interval. Two data abnormal low-value sections that are far apart are set to represent the natural fracture development interval, namely ① and ②. In the attached Figure 1 data group shown in Attachment B, which represents the acoustic travel time, it is generally considered that the high-value area of the acoustic travel time can represent the natural fracture development interval. Two data abnormal high-value sections that are far apart are set to represent the natural fracture development interval, namely ① and ②. In the attached Figure 1 data group shown in Attachment C, which represents the rock density, two data abnormal low-value sections that are close to each other are set to represent the natural fracture development interval, namely ① and ②. In the attached Figure 1 data group shown in Attachment D, which represents the rock density, 6 sections with relatively small thicknesses (relative to the attached Figure 1For A and 1C, the abnormally low-value segments of the data represent the intervals with developed natural fractures. The density and acoustic travel time of the formation without developed natural fractures are set to 2.6 g / cm3 and 40 μs / ft respectively, while the density and acoustic travel time of the intervals with developed natural fractures are set to 1.3 g / cm3 and 80 μs / ft respectively. Perform R / S variable-scale fractal analysis on the above four data groups of the prior art and plot the curve of lg[R(N) / S(N)] against N (which can be converted to depth). Further, calculate the first derivative and second derivative of lg[R(N) / S(N)] with respect to N according to the finite difference method. H should be the slope of lg[R(N) / S(N)] with respect to lg(N), but lg(N) is not linearly related to depth, and the characteristics of the intervals with developed natural fractures cannot be intuitively displayed when corresponding the R / S variable-scale fractal analysis results with depth. Therefore, this application reflects the correlation between the R / S variable-scale fractal analysis results and the intervals with developed natural fractures by studying the slope of lg[R(N) / S(N)] with respect to N. The curve of lg[R(N) / S(N)] against N (which can be converted to depth), the first and second derivative curves of lg[R(N) / S(N)] with respect to N are respectively shown as in Figure 1 the green, purple and orange curves in A.

[0050] In Figure 1 A and 1B, the first derivative of lg[R(N) / S(N)] with respect to N at the natural fracture interval ① is significantly increased compared with the formation without developed natural fractures, while the first derivative is significantly decreased at the natural fracture interval ②. When the spacing of the intervals with developed natural fractures is small, such as in Figure 1 C, the first derivative of the natural fracture interval ① is also significantly increased, while the first derivative of the natural fracture interval ② is partly less than and partly greater than the first derivative of the formation without developed natural fractures. And in Figure 1 D where the thickness of the intervals with developed natural fractures is small and the spacing is small, except for the natural fracture interval ①, the first derivatives of all the intervals with developed natural fractures are less than the first derivative of the formation without developed natural fractures. Therefore, there is no obvious corresponding relationship between the decreasing segment of the Hurst exponent H and the intervals with developed natural fractures.

[0051] In the prior art, in the calculation formula of the R / S variable-scale analysis technology, when N takes values from 1 to the total number of data points, lg[R(N) / S(N)] is always calculated starting from the top of the target interval. Therefore, the response of the abnormal data segment in the deep part of the target layer is significantly weakened, and the recognition ability of the top interface of the fracture-developed interval in the target layer is stronger than that of the bottom interface of the fracture-developed interval, which limits the recognition accuracy of the thickness of the fracture-developed interval. Figure 1There are also significant differences in the recognition accuracy of natural fracture development intervals at different depths, which is particularly evident in formations with a large number of natural fracture development intervals. For example, in Appendix Figure 1 D, the law of the first derivative of lg[R(N) / S(N)] with respect to N for 6 natural fracture development intervals is as follows: the first derivative of interval ① is greater than that of the formation without natural fractures, and as the depth increases, the first derivative decreases rapidly; the first derivatives of intervals ②-⑥ are all less than that of the formation without natural fractures, and as the depth increases, the first derivative decreases slowly. The change amplitude of the first derivative at the top interface of all natural fracture development intervals is significantly greater than that at the bottom interface, indicating that in the existing technology, the calculation formula of the R / S variable-scale analysis technology cannot guarantee the recognition accuracy of both the top and bottom interfaces of natural fracture development intervals simultaneously.

[0052] In one embodiment, as Figure 2 shown, a natural fracture identification method based on the bivariate variable-scale fractal technology is provided, including the following steps:

[0053] Step S1: Select the target logging curve of the target formation, and obtain the logging data corresponding to the target logging curve, denoted as the first logging data; reverse the first logging data according to the depth to obtain the second logging data.

[0054] In step S1, select the target formation according to the actual needs of oil and gas field development, and clarify the depth distribution characteristics of the top interface and the bottom interface of the target formation in the study area. Mainly obtain the depth of the top interface and the bottom interface of the target formation from the drilling data of each well in the study area.

[0055] Various parameters such as the mechanical properties, quantity, and occurrence of fractures can be directly counted from core or imaging logging data, which is the most intuitive manifestation of the fracture development characteristics of the reservoir. However, due to its high cost, it is relatively scarce. In contrast, conventional logging data is usually more common. Therefore, count the number of wells drilled through the target formation in the study area, and divide the wells into two categories. The first category is the wells that have at least one of core and imaging logging data in the target formation and also have conventional logging data, and the second category is the wells that only have conventional logging data in the target formation.

[0056] For the first category of wells, if the logging of this category has imaging logging data, preferentially carry out natural fracture identification, characterization, and statistics based on the imaging logging data. If this category of wells does not have imaging logging, use the core data to carry out natural fracture identification, characterization, and statistics. Based on the imaging logging (core data), describe the natural fractures of the target formation in the study area for each well and count the number of fractures to obtain the relationship between the number of natural fractures per well and the actual depth.

[0057] The relationship between the number and depth of natural fractures obtained in a single well. The natural fracture line density is calculated according to the following formula (1):

[0058]

[0059] Wherein, MD represents the depth at which the natural fracture line density is to be calculated; FD represents the relationship curve between the number of natural fractures in a single well and the depth; FD(MD) represents the measured value of the natural fracture line density, in fractures / m; CL represents the cumulative number curve of natural fractures, that is, the cumulative number curve made according to the relationship between the number of natural fractures in a single well and the depth obtained in step S22; W is the measurement window length, usually a constant greater than 1.

[0060] For the first type of drilling wells, the natural fracture line density curve FD(MD) has been calculated. Since these drilling wells have both the natural fracture line density curve and the conventional logging curve, the sensitivity analysis of the conventional logging curve to the natural fracture line density can be carried out. Currently, the more commonly used methods are linear goodness-of-fit analysis, Pearson correlation analysis, and Spearman rank correlation analysis. All three methods can be used to study the correlation between the conventional logging curve and the natural fracture line density curve. The conventional logging curve with a high correlation with the natural fracture line density is the conventional logging curve that is more sensitive to the development of natural fractures. Taking the Spearman rank correlation coefficient as an example, its formula is:

[0061]

[0062] Where Rs(n) is the Spearman rank correlation coefficient; n is the number of data sets used for the quantitative relationship analysis of two variables; di is the sorting difference of the i-th variable set (X, Y). Each variable set (X, Y) is changed to a variable set (X’, Y’) by the following formulas (3) and (4) so that there are no variables with the same rank in X (or Y):

[0063] X' = X + Rand() / 10 8 (3)

[0064] Y' = Y + Rand() / 10 8 (4)

[0065] Wherein, Rand() is a random number between 0 and 1.

[0066] Taking the calculation results of the Spearman rank correlation coefficient as an example to illustrate how to select the conventional logging curves that are most sensitive to the degree of fracture development: The negative and positive values of the Spearman rank correlation coefficient Rs(n) calculated according to formulas (2) to (4) represent negative and positive correlations of the variable groups under study, respectively. The closer the absolute value is to 1, the stronger the negative or positive correlation between variables. Therefore, after calculating the Spearman rank correlation coefficients between various conventional logging data and the natural fracture line density, the conventional logging curve corresponding to the highest absolute value of the Spearman rank correlation coefficient can be selected as the optimal conventional logging curve for reflecting the natural fracture development degree of the target layer.

[0067] In one embodiment, the depth range of the target layer series is 0 - 500m. And through correlation analysis, it is considered that the correlation between the acoustic travel time logging data (AC) and the natural fracture line density is the highest (positive correlation). Therefore, the AC curve is selected as the optimal conventional logging curve. The relationship between the AC data and depth is as shown in the appendix Figure 3 In the figure, the five high - value segments numbered ① - ⑤ of the AC data are closely related to the development of natural fractures.

[0068] The calculation of the fractal dimension D(m,n) in the prior art is from the top interface of the target layer series to the bottom interface, with a strict order from shallow to deep. This may lead to differences in the recognition accuracy of the top and bottom interfaces of the natural fracture development intervals, and differences in the recognition accuracy of natural fractures in the shallow and deep parts of the target layer series. Therefore, in order to eliminate the influence of the current research order on the research results, it is necessary to conduct research on the conventional logging data from deep to shallow (i.e., reverse - order research). Reverse the conventional logging data of the target layer series in each well according to depth. All the original conventional logging data located in the shallow part are transformed into conventional logging data located in the deep part, and all the original conventional logging data located in the deep part are transformed into conventional logging data located in the shallow part. For example, in a well where the depth range of the target layer series is 500m to 1000m, after the conventional logging data is reversed, the original conventional logging data at a depth of 500m becomes the conventional logging data at a depth of 1000m, the original conventional logging data at a depth of 1000m becomes the conventional logging data at a depth of 500m, the original conventional logging data at a depth of 600m becomes the conventional logging data at a depth of 900m, the original conventional logging data at a depth of 800m becomes the conventional logging data at a depth of 700m, and so on.

[0069] Step S2: By introducing the variable - scale fractal technology of bivariate variables, calculate the range R and standard deviation S of the first logging data and the second logging data at different starting points m and sliding interval lengths n based on the R / S analysis.

[0070] In the traditional R / S variable scale analysis technique, the total number N of data points of the logging curve is the only independent variable, and the range of u in the formula is 0 to N, that is, the analysis always starts from the top interface of the target interval, and the length of the analysis interval is always equal to N. In the improved R / S variable scale analysis technique proposed by the present invention, the single variable in the traditional method is extended to two variables m and n, where m is the starting point of the R / S variable scale analysis, and n is the length of the sliding window, that is, the number of data points within the sampling interval. The variables m and n need to satisfy the following conditions:

[0071]

[0072] where N is the total number of conventional logging sampling points of the target interval.

[0073] The initial value and increment of the sliding interval length n should be determined comprehensively according to the thickness of the target interval and the requirements of experimental accuracy. The two are usually equal. The smaller the initial value and increment of the sliding interval length n, the higher the prediction accuracy, and the thinner natural fracture development intervals can be identified, but the computational workload also increases accordingly. The initial value and increment of the sliding interval length n are usually taken as 1 / 50 to 1 / 20 of the total thickness of the target interval.

[0074] In one embodiment, the total thickness of the target interval system is 500 m, and the initial value and increment of the sliding interval length n are set to 10 m.

[0075] After the single variable in the traditional method is extended to two variables m and n, the formulas of the traditional R / S variable scale analysis technique can be respectively rewritten as formulas (6) and (7) with two independent variables:

[0076]

[0077] Step S3: Calculate the Hurst exponent H of the first logging data and the second logging data, and respectively obtain the fractal dimensions of the first logging data and the second logging data; the Hurst exponent H is obtained by calculating the slope of the curve of lg[R(m,n) / S(m,n)] with respect to n; the fractal dimension obtained according to the first logging data is denoted as the first fractal dimension.

[0078] In step S3, the fractal dimension D presented in the prior art is directly related to the slope of the curve of lg[R(N) / S(N)] with respect to lg(N), but the relationship between lg(N) and depth is not linear, and the correlation between the fractal dimension D and depth cannot be intuitively reflected. In order to make the correlation between the fractal dimension D and depth more intuitive, the fractal dimension D is improved to formula (12), that is, the expression of the fractal dimension D(m,n) is:

[0079]

[0080] In formula (8), the Hurst exponent H is rewritten as the slope of the curve of lg[R(m,n) / S(m,n)] with respect to n, and n is in a strictly linear positive correlation with depth. Therefore, the corresponding relationship between the fractal dimension D(m,n) and depth is more intuitive. Each pair of starting point m and window n can correspond to a fractal dimension D(m,n).

[0081] Step S4: Convert the fractal dimension obtained from the second logging data into the same coordinate system as the fractal dimension obtained from the first logging data through the coordinate conversion formula, and denote the fractal dimension corresponding to the second test data after coordinate conversion as the second fractal dimension.

[0082] In step S4, perform an R / S variable scale analysis on the conventional logging data (second logging data) after deep and shallow inversion to obtain the relationship and distribution of the fractal dimension D(m’,n’) with respect to m’ and n’. Here, m’ and n’ are respectively the starting point and window when performing the R / S variable scale analysis on the inverted conventional logging data, and do not correspond to the starting point m and window n of the original conventional logging data. Therefore, after the following processing on m’ and n’, convert m’ and n’ into the same coordinate system as when studying in the positive order:

[0083]

[0084] Dr(m’,n’) is also thus converted to Dr(N + 1 - m, N - n), and is in the same coordinate system as D(m,n). Plot Dr(N + 1 - m, N - n) and D(m,n) as a plan view with respect to m and n in a plane rectangular coordinate system with m as the horizontal axis and n as the vertical axis, and the two-dimensional fractal dimension D contour map of “multiple starting points, multiple sliding windows, and positive and negative order superposition” can be plotted.

[0085] When the total thickness of the target formation is 500m and the initial value and increment of the sliding window length n are set to 10m, the two-dimensional fractal dimension D contour map calculated based on the AC data is as Figure 4 shown.

[0086] Step S5: Draw a two-dimensional fractal dimension D contour map according to the first fractal dimension and the second fractal dimension.

[0087] In step S5, first, draw a diagonal line to separate the fractal dimension D obtained from the analysis of the normal well logging data in the positive sequence from the fractal dimension D obtained from the analysis of the normal well logging data in the reverse sequence. Let the initial value of the sliding window length n be n0. Then, the abscissa of the upper left corner endpoint of the diagonal line is n0, and the ordinate is N - n0. The abscissa of the lower right corner of the diagonal line is N - n0, and the ordinate is n0. The triangular contour map below the diagonal line corresponds to the fractal dimension D obtained from the analysis of the normal well logging data in the positive sequence, and the triangular contour map above the diagonal line corresponds to the fractal dimension D obtained from the analysis of the normal well logging data in the reverse sequence. In the contour map of the two-dimensional fractal dimension D, the numerical region with a fractal dimension lower than 1 is filled with yellow to red, and the numerical region with a fractal dimension greater than or equal to 1 is filled with light blue to dark blue.

[0088] The abscissa of the upper left corner endpoint of the diagonal line of the two-dimensional fractal dimension D contour map is 10m, and the ordinate is 490m. The abscissa of the lower right corner of the diagonal line is 490m, and the ordinate is 10m. The numerical region with a fractal dimension lower than 1 is filled with yellow to red, and the numerical region with a fractal dimension greater than or equal to 1 is filled with light blue to dark blue, as Figure 4 shown.

[0089] Step S6: Select the top interface of the natural fracture development interval in the first fractal dimension region and the bottom interface in the second fractal dimension region in the two-dimensional fractal dimension D contour map; the rules for selecting the top interface and the bottom interface of the natural fracture development interval include: near the starting point m*, when the sliding interval length n is within a specific interval (p, q), if the fractal dimension on the left side of the starting point m* is greater than 1 and the right side is less than 1; and when n is not within the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the top interface of the natural fracture development interval; and if when the sliding window length n is within a specific interval (p, q), if the fractal dimension on the left side of the starting point m* is less than 1 and the right side is greater than 1; and when n is not within the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the bottom interface of the natural fracture development interval;

[0090] In step S6, for the region below the diagonal line in the two-dimensional fractal dimension D contour map, if for a specific starting point m of the R / S analysis, there exists a non-empty numerical interval (p, q) such that the following inequality group (10) holds when n ∈ (p, q), and when n ∉ (p, q), D(m, n) is greater than 1, then it is considered that the vertical line (m = m*) represented by the starting point m at this time can be used as an alternative effective starting boundary line for the abnormal region with a fractal dimension D lower than 1, that is, the alternative top interface of the natural fracture development interval.

[0091]

[0092] Where m- is the left limit of m (a value infinitely close to m but less than m), and m+ is the right limit of m (a value infinitely close to m but greater than m). If the fractal dimension low-value color block (red block) with this effective starting boundary line as the left boundary intersects with the horizontal line represented when the sliding window length n is equal to 2*n0, then the alternative effective starting boundary line represented by m = m* can be further selected as the effective starting boundary line (m = m**).

[0093] For the area above the diagonal of the two-dimensional fractal dimension D contour map, if for any m, there exists a non-empty numerical interval (p, q) such that the following inequality group (11) holds when n ∈ (p, q), and when at this time, D(m, n) is greater than 1, then it is considered that the vertical line (m = m*) represented by the starting point m at this time can be used as the alternative effective termination boundary line of the abnormal area where the fractal dimension D is less than 1, that is, the bottom interface of the alternative natural fracture development section.

[0094]

[0095] Where m- is the left limit of m (a value infinitely close to m but less than m), and m+ is the right limit of m (a value infinitely close to m but greater than m). If the fractal dimension at the intersection of the fractal dimension low-value color block (red block) with this effective starting boundary line as the right boundary and the horizontal line represented when the sliding window length n is equal to N - 2*n0 is equal to 1, then the alternative effective termination boundary line represented by m = m* can be further selected as the effective termination boundary line (m = m**). Thus, all effective starting boundary lines and effective termination boundary lines can be quantitatively identified, corresponding to the top interface and bottom interface of the natural fracture development section respectively.

[0096] As Figure 4 shown, for the vertical line L1-1, when n ∈ (p, q), the fractal dimension D on the left side of the vertical line L1-1 is greater than 1, and the fractal dimension D on the right side is less than 1, satisfying the inequality group (14). Therefore, the vertical line L1-1 can be used as the alternative top interface of the natural fracture development section, and the fractal dimension D low-value color block with the vertical line L1-1 as the left boundary intersects with the white short horizontal line n = 2*n0. Therefore, it can be determined as the top interface of the natural fracture development section. Further, the top interfaces and bottom interfaces of 5 natural fracture development sections can be respectively identified. The top interfaces are as shown by L1-1, L2-1, L3-1, L4-1, and L5-1 in the appendix Figure 4 and the bottom interfaces are as shown by L1-2, L2-2, L3-2, L4-2, and L5-2 in the appendix Figure 4 respectively.

[0097] Step S7: Based on the top and bottom interfaces of the selected natural fracture development intervals, determine the natural fracture development intervals, and obtain the number of natural fracture development intervals, the top interface depth, the bottom interface depth, and the thickness.

[0098] In step S7, the top and bottom interfaces of the quantitatively identified natural fracture development intervals can be used to quantitatively identify the natural fracture development intervals. The sampling starting point m is increased from 0 to N. When adjacent m1 and m2 are respectively the starting point and the ending point of a natural fracture development interval and m1 < m2, it is considered that m1 and m2 are respectively the starting point and the ending point of a certain natural fracture development interval, and the thickness of this natural fracture development interval is m2 - m1. If the top and bottom interface depths of the natural fracture development interval adjacent to this natural fracture development interval are m3 and m4 respectively, then the distance between these two adjacent natural fracture development intervals is m3 - m2. Repeat this operation until reaching the bottom boundary of the single-well conventional logging data of the target formation in the study area, and the number, top interface depth, bottom interface depth, thickness and other parameters of the natural fracture development intervals of the target formation at the corresponding well can be quantitatively identified.

[0099] In one embodiment, as Figure 4 shown, 5 natural fracture development intervals are identified. The top interface depths are the m values m(L1 - 1), m(L2 - 1), m(L3 - 1), m(L4 - 1) and m(L5 - 1) corresponding to the vertical lines L1 - 1, L2 - 1, L3 - 1, L4 - 1 and L5 - 1 respectively, the bottom interface depths are the m values m(L1 - 2), m(L2 - 2), m(L3 - 2), m(L4 - 2) and m(L5 - 2) corresponding to the vertical lines L1 - 2, L2 - 2, L3 - 2, L4 - 2 and L5 - 2 respectively, the thicknesses of the 5 natural fracture development intervals are m(L1 - 2) - m(L1 - 1), m(L2 - 2) - m(L2 - 1), m(L3 - 2) - m(L3 - 1), m(L4 - 2) - m(L4 - 1) and m(L5 - 2) - m(L5 - 1) respectively, and the distances are m(L2 - 1) - m(L1 - 2), m(L3 - 1) - m(L2 - 2), m(L4 - 1) - m(L3 - 2) and m(L5 - 1) - m(L4 - 2) respectively.

[0100] In the above natural fracture identification method based on the bivariate variable-scale fractal technology, the logarithm of the ratio of the range to the standard deviation, lg[R(m,n) / S(m,n)], is used to characterize the fractal dimension D(m,n) by the slope with respect to n, which more intuitively reflects the slope of lg[R(m,n) / S(m,n)] with respect to the depth of the target layer. Thus, the characterized fractal dimension D(m,n) is depth-related and can more intuitively display the characteristics of the natural fracture development section. At the same time, the invention uses two variables, the initial position m of the sampling points and the sampling interval length n, to carry out the R / S variable-scale fractal analysis, which overcomes the drawback in the prior art that the recognition accuracy of natural fracture development sections at different depths varies greatly due to always starting the R / S variable-scale fractal analysis from the top interface of the target layer.

[0101] Moreover, the invention combines the R / S variable-scale fractal analysis of the normal sequence (from shallow to deep in depth) and the reverse sequence (from deep to shallow in depth) of the target layer's conventional logging, which overcomes the drawback in the prior art that the recognition accuracy of the top and bottom interfaces of natural fractures varies greatly.

[0102] After the bivariate R / S variable-scale fractal technology proposed by the invention draws the fractal dimension plane, if there exists an interval (p,q) at a certain starting point mm*, such that the fractal dimension D(m,n)<1 only holds when n∈(p,q), and D(m,n)>1 outside the interval, and the fractal dimension on the left side of the starting point m* is greater than 1 in a specific interval and less than 1 on the right side, then the vertical line represented by the starting point m* serves as the top interface of the natural fracture development section. If the fractal dimension on the left side of the starting point m* is less than 1 in a specific interval and greater than 1 on the right side, then the vertical line of this point m* serves as the bottom interface of the natural fracture development section. Judging the boundary surface in this way can improve the accuracy and reliability of the selection of the boundary line.

[0103] In one embodiment, a computer device is provided. The computer device can be a server. The computer device includes a processor, a memory, and a network interface connected through a system bus. Among them, the processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system, a computer program, and a database. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The network interface of the computer device is used to communicate with an external terminal through a network connection. When the computer program is executed by the processor, it realizes the above natural fracture identification method based on the bivariate variable-scale fractal technology.

[0104] In one embodiment, a computer-readable storage medium is further provided, on which a computer program is stored. When the computer program is executed by the processor, it involves all or part of the processes in the method of the above embodiment.

[0105] Those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above various methods. Among them, any reference to a memory, storage, database, or other medium used in the various embodiments provided in the present application can include at least one of non-volatile and volatile memories. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical memory, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0106] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the various technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope described in this specification.

[0107] The above-described embodiments merely represent several implementation manners of the present application. The description is relatively specific and detailed, but it should not be construed as a limitation on the scope of the invention patent. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present application, several modifications and improvements can still be made, and these all belong to the protection scope of the present application. Therefore, the protection scope of the patent of the present application should be subject to the appended claims.

Claims

1. A natural fracture identification method based on dual variable variable scale fractal technology, characterized in that: The method comprises: Selecting a target logging curve of a target layer system, and obtaining logging data corresponding to the target logging curve, which is recorded as first logging data; reversing the first logging data according to depth to obtain second logging data; By introducing the variable scale fractal technology of two variables, the range R and standard deviation S of the first logging data and the second logging data at different starting points m and sliding interval lengths n are calculated based on R / S analysis; Calculating the Hurst index H of the first well logging data and the second well logging data, thereby obtaining the fractal dimensions of the first well logging data and the second well logging data respectively; the Hurst index H is obtained by calculating the slope of the curve of lg[R(m,n) / S(m,n)] with respect to n; recording the fractal dimension obtained according to the first well logging data as the first fractal dimension; The fractal dimension obtained from the second well logging data is converted into the same coordinate system as the fractal dimension obtained from the first well logging data by using a coordinate conversion formula, and the fractal dimension corresponding to the second well logging data after the coordinate conversion is recorded as the second fractal dimension; Draw a two-dimensional fractal dimension D contour map based on the first fractal dimension and the second fractal dimension; The top interface of the natural fracture development layer section is selected in the first fractal dimension area of ​​the two-dimensional fractal dimension D contour map, and the bottom interface is selected in the second fractal dimension area; the selection rules for the top interface and the bottom interface of the natural fracture development layer section include: if near the starting point m*, when the sliding interval length n is within the specific interval (p, q), if the left fractal dimension of the starting point m* is greater than 1, and the right side is less than 1; and when n is not in the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the top interface of the natural fracture development layer section; and if when the sliding window length n is within the specific interval (p, q), if the left fractal dimension of the starting point m* is less than 1, and the right side is greater than 1; and when n is not in the specific interval (p, q), the fractal dimension corresponding to m* is also greater than 1, then the vertical line of the starting point m* is used as the bottom interface of the natural fracture development layer section; According to the top interface and bottom interface of the selected natural fracture development layer segment, the natural fracture development layer segment is determined, and the number of the natural fracture development layer segments, the top interface depth, the bottom interface depth and the thickness are obtained.

2. The natural fracture identification method based on dual variable scale fractal technology according to claim 1 is characterized in that: The target logging curve of the target formation is selected by: calculating the relationship between the number of natural fractures in a single well and the actual depth based on core data or imaging logging data; The natural fracture line density is calculated according to the relationship between the number of natural fractures in the single well and the actual depth, and the sensitivity analysis of different types of logging curves to the natural fracture line density is performed, and the logging curve with the largest sensitivity coefficient is selected as the target logging curve.

3. The natural fracture identification method based on dual variable scale fractal technology according to claim 2 is characterized in that: The sensitivity analysis of different types of well logging curves to natural fracture line density includes: using linear fit goodness of fit analysis, Pearson correlation analysis or Spearman rank correlation analysis to perform sensitivity analysis of different types of well logging curves to natural fracture line density.

4. The natural fracture identification method based on dual variable scale fractal technology according to claim 2 is characterized in that: The natural fracture line density is calculated based on the relationship between the number of natural fractures in a single well and the actual depth using the following formula: Wherein, MD represents the depth of the natural fracture line density to be calculated; FD represents the relationship curve between the number of natural fractures in a single well and the depth; FD(MD) represents the value of the natural fracture line density to be measured; CL represents the cumulative number curve of natural fractures; W is the measurement window length, which is usually a constant greater than 1.

5. The natural fracture identification method based on dual variable variable scale fractal technology according to claim 1 is characterized in that: The initial value of the sliding interval length n and its increment are selected from 1 / 50 to 1 / 20 of the total thickness of the target layer system.

6. The natural fracture identification method based on dual variable variable scale fractal technology according to claim 1 is characterized in that: The range R and standard deviation S of the first logging data and the second logging data at different starting points m and sliding interval lengths n are calculated by the following formula: Among them, R and S are the range and standard deviation respectively, Z is the target logging curve; N is the total number of data points of the target logging curve; u is the data point number for each range and standard deviation calculation in the range of 0 to N.

7. The natural fracture identification method based on dual variable variable scale fractal technology according to claim 1 is characterized in that: The Hurst exponent H of the first well logging data and the second well logging data is calculated to obtain the fractal dimensions of the first well logging data and the second well logging data respectively through the following formula: D(m,n)=2-H Among them, H is the Hurst index and D(m,n) is the fractional dimension.

8. The natural fracture identification method based on dual variable variable scale fractal technology according to claim 1 is characterized in that: The selection rule of the top and bottom interfaces of the natural fracture development layer section also includes: taking the starting point m*vertical line as the top interface of the natural fracture development layer section, recorded as the initial top interface, if the numerical region of the fractal dimension of the top interface as the left boundary is less than 1 and the horizontal line represented when the sliding window length n is equal to 2*n0 has an intersection, then the initial top interface is further determined to be the top interface of the final natural fracture development layer section; the no is the initial value of the sliding interval length; The vertical line of the starting point m* is taken as the bottom interface of the natural fracture development layer section, which is recorded as the initial bottom interface. If the fractal dimension of the initial interface as the numerical area of ​​the right boundary is lower than 1 and the fractal dimension of the horizontal line represented when the sliding window length n is equal to N-2*n0 is equal to 1, then the initial bottom interface is further determined to be the final bottom interface of the natural fracture development layer section.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 8 are implemented.

10. A computer program product comprising a computer program / instructions, characterized in that When the computer program / instructions are executed by a processor, the steps of the method according to any one of claims 1 to 8 are implemented.

Images