A natural fracture identification method based on dual-variable variable-scale fractal technology

Abstract

The present application discloses a natural fracture identification method based on dual-variable variable scale fractal technology, comprising: selecting a target logging curve of the target formation, obtaining forward logging data, and reversing its depth to obtain reverse logging data. Through R / S analysis, the range R and standard deviation S of the two sets of data at different starting points m and sliding interval lengths n are calculated, and then the Hurst index H and fractal dimension D are obtained. The fractal dimension of the reverse data is converted to the same coordinate system as the forward data, and a two-dimensional fractal dimension D contour map is drawn. According to the distribution of fractal dimension D, the top interface of the fracture development layer segment is selected in the forward data area, and the bottom interface is selected in the reverse data area. Then, through the top and bottom interfaces of the selected fracture development layer segment, the number of fracture development layer segments, the top interface depth, the bottom interface depth and the thickness are quantitatively obtained. Through the present application, the recognition accuracy of fracture development layer segments can be improved.

Description

Technical Field

[0001] The present application relates to the technical field of natural fracture identification, and in particular to a natural fracture identification method. Background Art

[0002] The comprehensive characterization of natural fractures in reservoirs is the basis for prediction and modeling of natural fractures in reservoirs, and is also an important reference for the optimization of reservoir engineering sweet spots and the formulation of exploration and development plans. At present, the methods for identifying and characterizing natural fractures mainly include geological and petrological methods, well logging methods, drilling and logging, and production dynamic data analysis methods. Among them, modern mathematical methods 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 conventional logging data has a relatively low acquisition cost among all types of logging data, and each type of logging curve has different degrees of response characteristics to the natural fracture development layer, which can be used to quantitatively identify the number of natural fracture development layers, the bottom interface depth, and the top interface depth. The essence of the R / S variable scale analysis technology is to simplify the identification of natural fractures to the identification of abnormal segments in the logging curve, where R and S are the range and standard deviation, respectively, as shown in the following formulas:

[0003]

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

[0005]

[0006] D=2-H

[0007] Moreover, the existing technology believes that the decreasing section of the curve slope H (i.e. the increasing section of the fractal dimension D) is the potential natural fracture development section. However, the current method has several limitations:

[0008] Limitation 1: Identifying fracture-developed intervals requires manual judgment based on curve morphology, which is labor-intensive and highly subjective. Although researchers have proposed an improved R / S-FD method for automated fracture identification, this method directly identifies segments with decreasing curve slope as concave, resulting in limited accuracy.

[0009] Limitation 2: In the existing technology, the natural fracture development layer is determined by the decreasing slope, that is, the Hurst index H. This determination method is not accurate.

[0010] Limitation 3: Traditional R / S variable scaling analysis technology always calculates lg[R(N) / S(N)] starting from the top of the target layer. Therefore, the response of the abnormal data segment located deep in the target layer is significantly weakened. In addition, the recognition ability of the top interface of the fracture-developed layer in the target layer is stronger than that of the bottom interface of the fracture-developed layer. This limits the thickness recognition accuracy of the fracture-developed layer.

[0011] Therefore, the above method for determining fracture development intervals is not accurate, and its determination accuracy is still difficult to meet the needs of actual research. Summary of the Invention

[0012] Based on this, in order to solve the above technical problems, a natural fracture identification method based on dual-variable variable scale fractal technology is provided to solve the problem of low accuracy in determining fracture development intervals in the existing technology.

[0013] In a first aspect, a natural fracture identification method based on a bivariate variable scale fractal technique comprises:

[0014] Selecting a target well logging curve of a target layer, and obtaining well logging data corresponding to the target well logging curve, which is recorded as first well logging data; inverting the first well logging data according to depth to obtain second well logging data;

[0015] By introducing a bivariate variable scale fractal technique, 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;

[0016] Calculating the Hurst exponent H of the first well logging data and the second well logging data to obtain fractal dimensions of the first well logging data and the second well logging data, respectively; the Hurst exponent H is obtained by calculating the slope of a curve of lg[R(m,n) / S(m,n)] with respect to n; and recording the fractal dimension obtained based on the first well logging data as a first fractal dimension;

[0017] 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 test data after the coordinate conversion is recorded as the second fractal dimension;

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

[0019] The top interface of the natural fracture development layer segment 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 and bottom interfaces of the natural fracture development layer segment include: if near the starting point m*, when the sliding interval length n is within the 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 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 segment; 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 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 segment;

[0020] 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 natural fracture development layer segments, the top interface depth, the bottom interface depth and the thickness are obtained.

[0021] In the above solution, optionally, selecting the target logging curve of the target formation includes: 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;

[0022] The natural fracture line density is calculated based on the relationship between the number of natural fractures in the single well and the actual depth, and 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.

[0023] In the above scheme, further optionally, the sensitivity analysis of different types of 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 logging curves to natural fracture line density.

[0024] In the above solution, further optionally, the natural fracture linear density is calculated according to the relationship between the number of natural fractures in the single well and the actual depth by the following formula:

[0025]

[0026] Where MD represents the depth of the natural fracture linear density to be calculated; FD represents the relationship between the number of natural fractures in a single well and the depth; FD(MD) represents the value of the natural fracture linear density to be measured; CL represents the cumulative number of natural fractures; and W is the measurement window length, which is usually 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 layer system.

[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; and u is the data point number within the range of 0 to N for each calculation of the range and standard deviation.

[0031] In the above solution, optionally, 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 by the following formula:

[0032]

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

[0034] Among them, H is the Hurst index and D(m,n) is the fractal dimension.

[0035] In the above scheme, optionally, the selection rule of the top and bottom interfaces of the natural fracture development layer section further 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 with a fractal dimension less than 1 of the top interface as the left boundary intersects with the horizontal line represented when the sliding window length n is equal to 2*n0, then the initial top interface is further determined to be the top interface of the final natural fracture development layer section; and no is the initial value of the sliding interval length;

[0036] The vertical line of the starting point m* is used as the bottom interface of the natural fracture development layer segment, which is recorded as the initial bottom interface. If the fractal dimension of the initial interface as the right boundary is lower than 1 in the numerical area 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 segment.

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

[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 dual-variable variable-scale fractal technology described in the first aspect.

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

[0040] The present application uses the slope of the logarithm of the ratio of the range and the standard deviation lg[R(m,n) / S(m,n)] 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, thereby making the characterized fractal dimension D(m,n) related to the depth, and more intuitively showing 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 length n of the sampling interval, to carry out R / S variable scale fractal analysis, which overcomes the disadvantage of the existing technology that the R / S variable scale fractal analysis always starts from the top interface of the target layer, resulting in large differences in the identification accuracy of natural fracture development sections at different depths.

[0041] In addition, the present invention adopts a method that combines the positive sequence (according to depth from shallow to deep) and the reverse sequence (according to depth from deep to shallow) R / S variable scale fractal analysis of conventional well logging of the target layer, which overcomes the disadvantage of large differences in the identification accuracy of the top and bottom interfaces of natural fractures in the existing technology.

[0042] After plotting the fractal dimension plane using the dual-variable R / S variable-scaling fractal technique proposed in this invention, if there exists an interval (p,q) at a certain starting point mm* such that the fractal dimension D(m,n) < 1 only within n∈(p,q) and D(m,n) > 1 outside the interval, and the fractal dimension to the left of the starting point m* is greater than 1 within a specific interval and less than 1 to the right, then the vertical line represented by the starting point m* is used as the top interface of the natural fracture development layer segment. If the fractal dimension to the left of the starting point m* is less than 1 within a specific interval and greater than 1 to the right, then the vertical line represented by the starting point m* is used as the bottom interface of the natural fracture development layer segment. This method of determining the boundary surface can improve the accuracy and reliability of boundary line selection. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 A comparison chart of the analysis results of four ideal data sets using the existing R / S variable scale fractal technology;

[0044] Figure 2 A flowchart of a natural fracture identification method based on a dual-variable variable-scale fractal technology is provided for one embodiment of the present application;

[0045] Figure 3 A graph showing the relationship between acoustic transit time (AC) logging data and depth of a target formation in a study area provided in one embodiment of the present application;

[0046] Figure 4 This is an analysis result of a fracture development layer using a natural fracture identification method based on a dual-variable variable scale fractal technology in one embodiment of the present application. DETAILED DESCRIPTION

[0047] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0048] In the description of this application: unless otherwise specified, the meaning of "plurality" is two or more. The terms "first", "second", "third", etc. in this application are intended to distinguish the objects referred to and do not have any special meaning in terms of technical connotation (for example, they 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 believed that the section where the slope, i.e., the Hurst index H, decreases is the section where natural fractures develop, but whether this identification is representative is still in doubt. Figure 1 As shown in the figure, four data sets are set to carry out numerical experiments. Figure 1 In the data set representing rock density shown in A, it is generally believed that the low-density area can represent the natural fracture development layer. Two data abnormal low-value segments far apart are set to represent the natural fracture development layer, namely ① and ②. Figure 1 In the data set representing the acoustic time difference shown in B, it is generally believed that the high value area of ​​the acoustic time difference can represent the natural fracture development layer. Two data abnormal high value segments far apart are set to represent the natural fracture development layer, namely ① and ②. Figure 1 In the data set representing rock density shown in C, two sections with abnormally low values ​​of data that are close to each other are set to represent the natural fracture development layer sections, namely ① and ②. Figure 1 In the data set representing rock density shown in D, six smaller thickness (relative to the adjacent Figure 1The abnormally low values ​​in the data (for A and 1C) represent layers with well-developed natural fractures. The density and acoustic transit time for layers without natural fractures were set to 2.6 g / cm³ and 40 μs / ft, respectively, while those for layers with well-developed natural fractures were set to 1.3 g / cm³ and 80 μs / ft, respectively. Existing techniques for R / S variable-scaling fractal analysis were performed on these four data sets, plotting the graph of lg[R(N) / S(N)] with respect to N (which can be converted to depth). The first and second derivatives of lg[R(N) / S(N)] with respect to N were further calculated using 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 correlated with depth. Therefore, when the R / S variable-scaling fractal analysis results are mapped to depth, the characteristics of layers with well-developed natural fractures cannot be intuitively displayed. Therefore, this application reflects the correlation between the R / S variable scale fractal analysis results and the natural fracture development interval by studying the slope of lg[R(N) / S(N)] with respect to N. The curve of lg[R(N) / S(N)] with respect to N (which can be converted into depth), the first and second order derivative curves of lg[R(N) / S(N)] with respect to N are shown in the attached figure. Figure 1 The green, purple, and orange curves in A.

[0050] In the attached Figure 1 The first-order derivative of lg[R(N) / S(N)] with respect to N at the natural fracture interval ① in A and 1B is significantly larger than that in the formation without natural fractures, while the first-order derivative is significantly reduced at the natural fracture interval ②. When the spacing between the natural fracture intervals is small, as shown in the attached figure, Figure 1 C. The first-order derivative of the natural fracture developed section ① also increases significantly, while the first-order derivative of the natural fracture developed section ② is partially smaller than the first-order derivative of the stratum without natural fractures, and partially larger than the first-order derivative of the stratum without natural fractures. Figure 1 In D, except for the natural fracture developed section ①, the first-order derivatives of all natural fracture developed sections are smaller than the first-order derivatives of the strata without natural fractures. Therefore, there is no obvious correspondence between the sections with decreasing Hurst exponent H and the sections with natural fractures.

[0051] In the existing R / S variable scaling analysis technology, when N ranges from 1 to the total number of data points, the calculation of lg[R(N) / S(N)] always starts from the top of the target layer. As a result, the response of abnormal data segments located deep in the target layer is significantly weakened. In addition, the recognition ability of the top interface of the fracture-developed layer in the target layer is stronger than that of the bottom interface of the fracture-developed layer. This limits the thickness recognition accuracy of the fracture-developed layer. Figure 1The results also show that the recognition accuracy of natural fracture development intervals at different depths is significantly different, which is particularly evident in the formations with more natural fracture development intervals. Figure 1 In D, the first-order derivatives of lg[R(N) / S(N)] with respect to N for the six naturally fractured sections show the following pattern: the first-order derivative of section ① is greater than that of the stratum without natural fractures, and decreases rapidly with increasing depth; the first-order derivatives of sections ② through ⑥ are all smaller than those of the stratum without natural fractures, and decrease slowly with increasing depth. The magnitude of the first-order derivative at the top interface of all naturally fractured sections is significantly greater than that at the bottom interface, indicating that the existing R / S variable scaling analysis formula cannot accurately identify both the top and bottom interfaces of naturally fractured sections.

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

[0053] Step S1: Select a target logging curve of a target formation, and obtain logging data corresponding to the target logging curve, which is recorded as first logging data; invert the first logging data according to depth to obtain second logging data.

[0054] In step S1, a target stratum is selected based on the actual needs of oil and gas field development, and the top and bottom depth distribution characteristics of the target stratum within the study area are determined. The top and bottom depths of the target stratum are primarily obtained from drilling data within the study area.

[0055] Core or imaging logging data can directly quantify various parameters such as fracture mechanical properties, number, and occurrence, providing the most intuitive representation of reservoir fracture development. However, due to its high cost, it is relatively scarce. In contrast, conventional logging data is generally more widely available. Therefore, we need to count the number of wells that encounter the target formation in the study area and divide the wells into two categories: the first category is wells that have at least one of core and imaging logging data, as well as conventional logging data, in the target formation; the second category is wells that only have conventional logging data in the target formation.

[0056] For the first type of wells, if imaging logging data is available, natural fracture identification, characterization, and statistics are prioritized based on these data. If imaging logging is not available, core data are used for this purpose. Based on imaging logging (core data), natural fractures are described and counted in individual wells within the target formations of the study area, revealing the relationship between the number of natural fractures in each well and their actual depth.

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

[0058]

[0059] Wherein, MD represents the depth of the natural fracture linear 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 linear density to be measured, in bars / m; CL represents the cumulative number curve of natural fractures, i.e., the cumulative number curve prepared based on 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, which is usually a constant greater than 1.

[0060] For the first type of wells, the natural fracture line density curve FD (MD) has been calculated. Since these wells have both natural fracture line density curves and conventional logging curves, a sensitivity analysis of the conventional logging curve to the natural fracture line density can be carried out. Currently, the most commonly used methods are linear fit goodness of fit analysis, Pearson correlation analysis, and Spearman rank correlation analysis. All three methods can be used to study the correlation between conventional logging curves and natural fracture line density curves. Conventional logging curves with a high correlation with natural fracture line density are conventional logging curves that are 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 to analyze the quantitative relationship between two variables; and di is the rank difference of the i-th variable set (X, Y). Each variable set (X, Y) is transformed into a variable set (X', Y') using the following formulas (3) and (4), so that no variables in X (or Y) have the same rank:

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

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

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

[0066] The calculation results of the Spearman rank correlation coefficient are used as an example to illustrate how to select the conventional logging curve that is most sensitive to the degree of fracture development: According to formulas (2) to 4), the negative and positive values ​​of the Spearman rank correlation coefficient Rs(n) respectively represent the negative and positive correlations between the studied variable groups. The closer its absolute value is to 1, the stronger the negative or positive correlation between the variables. Therefore, after calculating the Spearman rank correlation coefficient between each 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 reflecting the degree of natural fracture development in the target layer.

[0067] In one embodiment, the depth range of the target formation is 0-500m, and after correlation analysis, it is believed that the correlation between the acoustic transit time logging data (AC) and the natural fracture line density is the highest (positive correlation), so the AC curve is selected as the optimal conventional logging curve. The relationship between AC data and depth is shown in the attached figure. Figure 3 As shown in the figure, the five segments numbered ① to ⑤ are high-value segments of AC data, which are closely related to the development of natural fractures.

[0068] The fractal dimension D(m,n) in existing technologies is calculated from the top interface of the target formation to the bottom interface of the target formation, 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 layer, and differences in the recognition accuracy of natural fractures located in the shallow and deep parts of the target formation. Therefore, in order to eliminate the influence of the current research order on the research results, it is necessary to study the conventional logging data from deep to shallow (i.e., reverse the order). The conventional logging data of the target formation in each drilling well is reversed according to the depth. All conventional logging data originally located in the shallow part is converted to conventional logging data located in the deep part, and all conventional logging data originally located in the deep part is converted to conventional logging data located in the shallow part. For example, if the depth range of the target formation in a certain drilling is 500m to 1000m, then after the conventional logging data is inverted, the conventional logging data originally at a depth of 500m becomes the conventional logging data at a depth of 1000m, the conventional logging data originally at a depth of 1000m becomes the conventional logging data at a depth of 500m, the conventional logging data originally at a depth of 600m becomes the conventional logging data at a depth of 900m, the conventional logging data originally at a depth of 800m becomes the conventional logging data at a depth of 700m, and so on.

[0069] Step S2: By introducing a bivariate variable scale fractal technique, 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.

[0070] In traditional R / S variable scaling analysis, the total number of well logging data points, N, is the only independent variable, and u in the formula ranges from 0 to N. This means that the analysis always begins at 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 scaling analysis proposed in this paper, the single variable in the traditional method is upgraded to a pair of variables, m and n, where m is the starting point of the R / S variable scaling analysis, and n is the sliding window length, i.e., the number of data points within the sampling interval. The variables m and n must meet the following conditions:

[0071]

[0072] Where N is the total number of conventional logging sampling points in the target layer.

[0073] The initial value of the sliding interval length n and its increment should be determined based on the thickness of the target layer and the experimental accuracy requirements. The two are usually equal. Smaller initial values ​​of the sliding interval length n and its increment increase the prediction accuracy and the ability to identify thinner natural fracture zones, but the computational workload also increases accordingly. The initial value of the sliding interval length n and its increment are typically 1 / 50 to 1 / 20 of the total thickness of the target layer.

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

[0075] After upgrading the single variable in the traditional method to two variables m and n, the traditional R / S variable scaling analysis formula can be rewritten as formulas (6) and (7) with two independent variables:

[0076]

[0077] Step S3: Calculate the Hurst exponent H of the first well logging data and the second well logging data to obtain the fractal dimension of the first well logging data and the second well 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; the fractal dimension obtained based on the first well logging data is recorded 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). However, the relationship between lg(N) and depth is not linear, and the correlation between fractal dimension D and depth cannot be intuitively reflected. In order to make the correlation between fractal dimension D and depth more intuitive, fractal dimension D is improved to formula (12), that is, the expression of fractal dimension D(m,n) is:

[0079]

[0080] In formula (8), the Hurst exponent H is rewritten as the slope of the curve lg[R(m,n) / S(m,n)] with respect to n. Since n and depth are strictly linearly positively correlated, the corresponding relationship between the fractal dimension D(m,n) and depth is more intuitive. Each pair of starting point m and window n corresponds to a fractal dimension D(m,n).

[0081] Step S4: converting the fractal dimension obtained from the second well logging data into the same coordinate system as the fractal dimension obtained from the first well logging data using a coordinate conversion formula, and recording the fractal dimension corresponding to the second test data after coordinate conversion as the second fractal dimension.

[0082] In step S4, an R / S variable scaling analysis is performed on the conventional logging data (second logging data) that has undergone deep and shallow inversion, and the relationship and distribution of the fractal dimension D(m', n') with respect to m' and n' is obtained. Here, m' and n' are the starting point and window for performing the R / S variable scaling analysis on the inverted conventional logging data, respectively. They do not correspond to the starting point m and window n of the original conventional logging data. Therefore, after performing the following processing on m' and n', m' and n' are converted into the same coordinate system as in the forward sequence study:

[0083]

[0084] Dr(m', n') is thus converted to Dr(N+1-m, Nn), which is in the same coordinate system as D(m, n). Plotting Dr(N+1-m, Nn) and D(m, n) as a plane graph with respect to m and n in a rectangular coordinate system with m as the horizontal axis and n as the vertical axis creates a two-dimensional fractal dimension D contour map with multiple starting points, multiple sliding windows, and forward and reverse order superposition.

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

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

[0087] In step S5, a diagonal line is first drawn to distinguish the fractal dimension D obtained from the analysis of the forward-sequence conventional logging data from the fractal dimension D obtained from the analysis of the reverse-sequence conventional logging data. Assuming the initial value of the sliding window length n to be n0, the abscissa of the upper left corner 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 line below the diagonal line corresponds to the fractal dimension D obtained from the analysis of the forward-sequence conventional logging data, while the triangular contour line above the diagonal line corresponds to the fractal dimension D obtained from the analysis of the reverse-sequence conventional logging data. In the two-dimensional fractal dimension D contour line, the region with a fractal dimension less than 1 is filled with yellow to red, and the region with a fractal dimension greater than or equal to 1 is filled with light blue to dark blue.

[0088] The horizontal coordinate of the upper left corner of the diagonal line in the two-dimensional fractal dimension D contour map is 10m, and the vertical coordinate is 490m. The horizontal coordinate of the lower right corner of the diagonal line is 490m, and the vertical coordinate is 10m. The value area of ​​the fractal dimension less than 1 is filled with yellow to red, and the value area of ​​the fractal dimension greater than or equal to 1 is filled with light blue to dark blue, such as Figure 4 shown.

[0089] Step S6: selecting the top interface of the natural fracture development layer segment in the first fractal dimension area and the bottom interface in the second fractal dimension area in the two-dimensional fractal dimension D contour map; the selection rules for the top interface and the bottom interface of the natural fracture development layer segment include: if near the starting point m*, when the sliding interval length n is within the 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 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 segment; 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 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 segment;

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

[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 low-value fractal dimension block (red block) with the effective starting boundary line as the left boundary intersects 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 line of the contour map of the two-dimensional fractal dimension D, if for any m, there exists a non-empty numerical interval (p,q) such that when n∈(p,q) the following inequality group (11) holds, and when When D(m,n) is greater than 1, the vertical line represented by the starting point m (m=m*) is considered to be an alternative effective termination boundary line of the abnormal area with a fractal dimension D lower than 1, that is, the bottom interface of the alternative natural fracture development layer.

[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 of the intersection of the low-valued fractal dimension block (red block) with the effective starting boundary 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 ending boundary represented by m=m* can be further selected as the effective ending boundary (m=m**). In this way, all effective starting and ending boundaries can be quantitatively identified, corresponding to the top and bottom interfaces of the natural fracture development layer, respectively.

[0096] like Figure 4 As 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 an alternative top interface of the natural fracture development layer segment, and the low-value color block of the fractal dimension D with the vertical line L1-1 as the left boundary has an intersection with the white short horizontal line n=2*n0, so it can be determined as the top interface of the natural fracture development layer segment. Furthermore, the top interface and bottom interface of the five natural fracture development layer segments can be identified respectively. The top interface is shown in the attached figure. Figure 4 As shown in L1-1, L2-1, L3-1, L4-1 and L5-1, the bottom interface is as shown in the attached Figure 4 As shown in L1-2, L2-2, L3-2, L4-2 and L5-2.

[0097] Step S7: Based on the top and bottom interfaces of the selected natural fracture development intervals, 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 and ending points of a natural fracture development interval and m1 < m2, m1 and m2 are considered to be the starting and ending points 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, 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(L(2-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-mentioned natural fracture identification method based on dual-variable variable scale fractal technology, the slope of the logarithm of the ratio of the range and the standard deviation lg[R(m,n) / S(m,n)] with respect to n is used 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, thereby making the characterized fractal dimension D(m,n) related to the depth, and more intuitively showing 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 length n of the sampling interval, to carry out R / S variable scale fractal analysis, which overcomes the disadvantage of the existing technology that the R / S variable scale fractal analysis always starts from the top interface of the target layer, resulting in large differences in the identification accuracy of natural fracture development sections at different depths.

[0101] In addition, the present invention adopts a method that combines the positive sequence (according to depth from shallow to deep) and the reverse sequence (according to depth from deep to shallow) R / S variable scale fractal analysis of conventional well logging of the target layer, which overcomes the disadvantage of large differences in the identification accuracy of the top and bottom interfaces of natural fractures in the existing technology.

[0102] After plotting the fractal dimension plane using the dual-variable R / S variable-scaling fractal technique proposed in this invention, if there exists an interval (p,q) at a certain starting point mm* such that the fractal dimension D(m,n) < 1 only within n∈(p,q) and D(m,n) > 1 outside the interval, and the fractal dimension to the left of the starting point m* is greater than 1 within a specific interval and less than 1 to the right, then the vertical line represented by the starting point m* is used as the top interface of the natural fracture development layer segment. If the fractal dimension to the left of the starting point m* is less than 1 within a specific interval and greater than 1 to the right, then the vertical line represented by the starting point m* is used as the bottom interface of the natural fracture development layer segment. This method of determining the boundary surface can improve the accuracy and reliability of boundary line selection.

[0103] In one embodiment, a computer device is provided, which may be a server. The computer device includes a processor, memory, and a network interface connected via a system bus. The processor of the computer device is configured 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 computer program in the non-volatile storage medium. The network interface of the computer device is configured to communicate with an external terminal via a network connection. When executed by the processor, the computer program implements the aforementioned natural fracture identification method based on dual-variable 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 a processor, all or part of the processes in the above-mentioned embodiment method are involved.

[0105] Those skilled in the art will appreciate that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing the relevant hardware through a computer program, and 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-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may include read-only memory (ROM), magnetic tape, floppy disk, flash memory or optical memory, etc. Volatile memory may include random access memory (RAM) or external cache memory. As an illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM).

[0106] The technical features of the above embodiments can be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0107] The above-described embodiments merely represent several implementation methods of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that a person skilled in the art could make various modifications and improvements without departing from the spirit of the present application, all of which fall within the scope of protection of the present application. Therefore, the scope of protection of the present patent application shall be determined by the appended claims.

Claims

1. A natural fracture identification method based on bivariate variable scale fractal technology, characterized in that: The method comprises: Selecting a target well logging curve of a target layer, and obtaining well logging data corresponding to the target well logging curve, which is recorded as first well logging data; inverting the first well logging data according to depth to obtain second well logging data; By introducing the variable scale fractal technology of two variables, based on Analyze and calculate the first logging data and the second logging data at different starting points and the sliding interval length The extreme difference and standard deviation ; Calculate the Hurst index of the first logging data and the second logging data , thereby obtaining the fractal dimensions of the first logging data and the second logging data respectively; the Hurst index It is calculated by about The slope of the curve is obtained; the fractal dimension obtained according to the first logging data is recorded 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 two-dimensional fractal dimension based on first fractal dimension and second fractal dimension contour plots; In the two-dimensional fractal dimension The top interface of the natural fracture development layer section is selected in the first fractal dimension area of ​​the contour map, and the bottom interface is selected in the second fractal dimension area; the selection rules of the top interface and the bottom interface of the natural fracture development layer section include: Nearby, when the sliding interval length In a specific range If the starting point The left side fractal dimension is greater than 1, and the right side is less than 1; and when Not in a specific range hour, The corresponding fractal dimension is also greater than 1, then the starting point The vertical line is the top interface of the natural fracture development layer; if the sliding interval length In a specific range If the starting point The left fractal dimension is less than 1, and the right fractal dimension is greater than 1; and when Not in a specific range hour, The corresponding fractal dimension is also greater than 1, then the starting point The vertical line is taken as the bottom interface of the natural fracture development interval; 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, top interface depth, bottom interface depth and thickness of the natural fracture development layer segment are obtained.

2. The natural fracture identification method based on bivariate 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 based on the relationship between the number of natural fractures in the single well and the actual depth, and 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 bivariate 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: performing sensitivity analysis of different types of well logging curves to natural fracture line density using linear fit goodness of fit analysis, Pearson correlation analysis or Spearman rank correlation analysis.

4. The natural fracture identification method based on bivariate 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: ， in, Indicates the depth of the natural fracture linear density to be calculated; A curve showing the relationship between the number of natural fractures and depth in a single well; Indicates the linear density value of the natural fracture to be measured; represents the cumulative number curve of natural fractures; To measure the window length, a constant greater than 1 is usually taken.

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

6. The natural fracture identification method based on bivariate variable scale fractal technology according to claim 1 is characterized in that: The first logging data and the second logging data at different starting points are calculated by the following formula and the sliding interval length The extreme difference and standard deviation : ， ， in, and are the range and standard deviation, is the target logging curve; is the total number of data points of the target logging curve; From 0 to The data point number for each range and standard deviation calculation.

7. The natural fracture identification method based on bivariate variable scale fractal technology according to claim 1 is characterized in that: The calculation of the Hurst exponent of the first well logging data and the second well logging data , and the fractal dimensions of the first logging data and the second logging data are obtained by the following formula: ， in, is the Hurst exponent, is the fractal dimension.

8. The natural fracture identification method based on bivariate variable scale fractal technology according to claim 6 is characterized in that: The selection rules for the top and bottom interfaces of the natural fracture development layer section also include: The vertical line is taken as the top interface of the natural fracture development layer, which is recorded as the initial top interface. If the fractal dimension of the top interface as the left boundary is less than 1, the value of the sliding interval length is equal When the horizontal lines represented by have an intersection, the initial top interface is further determined to be the top interface of the final natural fracture development layer; is the initial value of the sliding interval length; Starting from the above The vertical line is taken as the bottom interface of the natural fracture development layer, which is recorded as the initial bottom interface. If the fractal dimension of the bottom interface as the right boundary is lower than the numerical area of ​​1 and the sliding interval length equal When the fractal dimension at the horizontal line represented by is equal to 1, it is further determined that the initial bottom interface is the bottom interface of the final natural fracture development layer.

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