Discrete fracture network automatic extraction method based on geometric model fitting algorithm

Abstract

The invention discloses a discrete fracture network automatic extraction method based on a geometric model fitting algorithm, and aims at a microseism event three-dimensional space position point cloud and direction information (trend, inclination angle and the like). Through the steps of initial model generation, direction consistency check, multi-model energy minimization iterative optimization, direction clustering and abnormal value elimination, plane updating and merging, convergence criterion termination and the like, automatic and accurate extraction of multiple crack planes is realized. According to the method, a target function is jointly constructed by using normal vectors and position errors of event points in the discrete fracture network, and abnormal values are eliminated or merged by fusing direction information in an iteration process, so that the method has relatively high anti-noise performance and convergence stability. According to the method, the automatic extraction efficiency and precision of the discrete fracture network are effectively improved, and the method is suitable for fracture analysis and modeling of a complex fracture structure or an unconventional reservoir.

Description

Technical Field

[0001] The present invention relates to the technical field of microseismic monitoring and oil and gas reservoir geological engineering, and in particular to a method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm. Background Art

[0002] In the development of oil and gas fields and geological structure research, discrete fracture networks (DFN) are of great significance to reservoir transformation and fluid migration. Traditional methods usually rely on manual interpretation or core analysis, which are highly subjective and time-consuming. With the development of microseismic monitoring technology, higher-density three-dimensional spatial event data and their directional information (strike, dip, etc.) can be obtained. However, in an environment with high noise and uncertainty, it is still challenging to automatically and accurately extract multiple superimposed or intersecting fracture planes from a three-dimensional point cloud. To this end, the present invention proposes to integrate position and direction information into the same iterative optimization framework, and automatically output multiple fracture planes through multi-model energy minimization iteration, direction clustering and outlier removal, thereby improving the accuracy and efficiency of discrete fracture network analysis. Summary of the Invention

[0003] The purpose of the present invention is to solve the above problems. A method for automatic extraction of discrete fracture networks based on a geometric model fitting algorithm is designed. The coordinate and direction information of microseismic events is utilized to reduce threshold dependence while improving the accuracy of multi-fracture plane analysis through adaptive model merging and outlier removal. The method is widely applicable to oil and gas reservoir development and geological structure analysis.

[0004] To achieve the above objectives, this application provides the following technical solutions:

[0005] A method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm comprises the following steps:

[0006] S1. Obtain a three-dimensional point cloud of microseismic events from microseismic monitoring data, wherein each microseismic event point includes at least one direction parameter in addition to coordinate information, and the direction parameter includes strike and dip;

[0007] S2. Use a random sampling method to perform plane fitting on the microseismic event points in step S1 to screen out possible fracture planes. The possible fracture planes and the points they cover are incorporated into the initial plane model set. For each possible fracture plane, its normal vector is calculated and compared with the average direction of the event points assigned to the plane. If the direction difference exceeds a set angle threshold, the possible fracture plane is eliminated.

[0008] S3. Assign the initial plane in step S2 as a label to each microseismic event point in the 3D point cloud, and iteratively update the objective function based on multi-model energy minimization. Each iteration re-estimates the plane equation parameters and minimizes the objective function based on the distance from the point to the plane and the consistency of the direction.

[0009] S4. During the iteration of step S3, cluster the directional information of the points under the same plane label and mark the points with significant differences with outlier labels to avoid interfering with plane fitting;

[0010] S5. Perform merge re-estimation for plane models whose normal vectors are approximately parallel and close in space; regenerate new plane models for those with a large number of points and a high directional concentration in the outlier labels in step S4;

[0011] S6. When the change in the objective function in step S3 remains small or converges to a stable state within a preset number of iterations, the iteration is stopped, and the final multi-crack plane model and the event information contained in each plane are output.

[0012] Furthermore, the objective function E(ξ) of multi-model energy minimization in step S3 is:

[0013]

[0014] Among them, p i Represents the three-dimensional coordinates (x i ,y i ,z i ); L(p i ) indicates p i The plane label of the plane to which it belongs; ||p i -L(p i )|| represents the distance between the event point and the assigned plane; N represents the set of adjacent points, (p i ,p j )∈N means (p i ,p j ) are adjacent to each other in the geometric space; δ(·) represents the indicator function; when p i With p j When corresponding to different plane labels, L(p i )≠L(p j ) takes 1, when p i With p j When corresponding to the same plane label, L(p i )≠L(p j ) is 0; λ is the weight of the smoothing term; P Φ Represents a set of outliers that cannot be properly fitted by any plane; L ΦRepresents the outlier label, which does not participate in the plane update and will be penalized if assigned to an outlier; is the outlier penalty coefficient,

[0015] Furthermore, in step S4, a density-based spatial clustering algorithm is used to cluster the directional information. The angle threshold ε and the minimum number of points are selected, and the event points assigned to the same plane label are classified according to the differences in strike and dip. Only the main cluster with the most consistent direction with the plane normal vector is retained, and the remaining event points are classified as outlier labels.

[0016] Furthermore, in step S5, the plane merging conditions include that the angle between the normal vectors of the two planes is less than 5° and the closest distance between the two planes is less than 2 meters, and the merging is performed; after the event point sets of the two are merged, the least squares method is used to re-estimate the new plane.

[0017] Furthermore, in step S5, after several iterations, the event points in the outlier label will try to generate a new plane model again based on their quantity scale and direction range. If the matching conditions are still not met, they will eventually be retained in the outlier label and marked as outliers when output.

[0018] Furthermore, in step S6, the iterative convergence judgment condition includes that the change of the objective function is less than a set threshold, or the number of iterations reaches an upper limit and the event point allocation label tends to be stable.

[0019] In summary, the technical effects and advantages of the present invention are as follows: the present invention integrates "event position-fracture direction" into the same iterative optimization framework, and realizes the simultaneous extraction of multiple fracture planes by minimizing the objective function containing the geometric model and directional error; this method uses the normal vector and position error of the event point in the discrete fracture network to jointly construct the objective function, and integrates the directional information during the iterative process to eliminate or merge outliers, with strong anti-noise performance and convergence stability; the present invention effectively improves the efficiency and accuracy of automatic extraction of discrete fracture networks, and is suitable for fracture analysis and modeling of complex fault structures or unconventional reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present application. For those skilled in the art, other drawings can be obtained based on these drawings without paying any creative work.

[0021] Figure 1 is a flow chart of a method according to an embodiment of the present invention;

[0022] Figure 21 is a schematic diagram of plane fitting of multiple cracks according to an embodiment of the present invention. DETAILED DESCRIPTION

[0023] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0024] In addition, the technical features involved in different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0025] The present invention is based on the position point cloud of microseismic events in three-dimensional space, combined with the directional data (including strike and dip) obtained based on moment tensor or focal mechanism analysis. The core idea of the present invention is to integrate "event position-crack direction" into the same iterative optimization framework, and realize the simultaneous extraction of multiple fracture planes by minimizing the objective function that includes geometric model and directional errors.

[0026] This embodiment proposes a method for automatically extracting discrete crack networks based on a geometric model fitting algorithm. Figure 1-Figure 2 As shown, the following steps are included:

[0027] S1. Obtaining a three-dimensional point cloud of microseismic events from microseismic monitoring data, wherein each microseismic event point includes at least one direction parameter in addition to coordinate information;

[0028] Obtain 3D point cloud data {p1,p2,...,p n}, where p i Contains (x i ,y i ,z i ) coordinate information, and obtain the corresponding direction information for each point (for example, the direction of θ i and inclination angle Φ i ).

[0029] The microseismic monitoring data used in this example comes from an unconventional reservoir block. During the hydraulic fracturing process, three-component geophones collected raw data for location analysis and extraction of significant microseismic events. Each event contains three-dimensional spatial coordinates (x, y, z). The strike and dip information (which can also be used to characterize the main fracture plane) is obtained using focal mechanism solutions or moment tensor inversion.

[0030] Before officially entering the automatic extraction process, the raw data needs to be preprocessed, including detecting and eliminating events with complete errors or positioning anomalies due to low signal-to-noise ratio. Events with excessively large positioning errors (greater than 10 meters) or no valid direction information will not be included in this analysis.

[0031] S2. Use a random sampling method to perform plane fitting on the microseismic event points in step S1 to screen out possible fracture planes. The possible fracture planes and the points they cover are incorporated into the initial plane model set. For each possible fracture plane, its normal vector is calculated and compared with the average direction of the event points assigned to the plane. If the direction difference exceeds a set angle threshold, the possible fracture plane is eliminated.

[0032] In this embodiment, valid microseismic event points are used as input data. Three microseismic event points are randomly selected to fit the plane, and the distance between the fitting plane and the threshold d is counted. ini If the number of points in the fitting plane exceeds a certain minimum threshold (for example, 10 points), the fitting plane is considered as a possible crack plane, and the possible crack plane and the points covered by it are included in the initial plane model set (initial label set L0). The above sampling and detection operations are repeated to obtain the initial plane model. ini There is no specific limit on the value. In actual application, the fracture development situation at different fracturing scales is extracted, and the threshold d ini The value will change.

[0033] For each possible crack plane, calculate its normal vector n and compare it with the average direction of the event points assigned to that plane. If the angle between the two exceeds 30°, the plane is considered to be severely deviated from the direction data and is not suitable for further iteration. It must be discarded to avoid serious discrepancies with the direction data. Here, the normal vector n can be calculated using Equation 6; the average direction of the event points assigned to the plane is calculated using Equation 3.

[0034] S3. Assign the initial plane in step S2 as a label to each microseismic event point in the 3D point cloud, and iteratively update the objective function based on multi-model energy minimization. Each iteration re-estimates the plane equation parameters and minimizes the objective function based on the distance from the point to the plane and the consistency of the direction.

[0035] The initial plane is regarded as a label (such as L1, L2, L3, L4...) and assigned to each event point in the 3D point cloud. Each event point is temporarily assigned a best matching label based on the geometric distance and direction difference with each plane model.

[0036] After the initial plane model is generated, it is necessary to fine-tune the plane model and its corresponding point allocation through iteration. In order to measure the matching degree between the plane and the point, the present invention designs the following objective function: iand its label L(p i ) is optimized and solved, and the objective function is minimized in each iteration. The objective function is shown in Formula 1:

[0037]

[0038] Among them, p i Represents the three-dimensional coordinates (x i ,y i ,z i ); L(p i ) indicates p i The geometric model of the plane (or plane label) can be represented by the normal vector and position of the plane equation (e.g., Ax+By+Cz+D=0), or simply expressed as the geometric solution of the label; || p i -L(p i )|| represents the distance between the event point and the assigned plane (Euclidean distance or perpendicular distance from the point to the plane). The smaller the distance, the better the p i The higher the matching degree with the selected plane, the better. N represents the set of adjacent points. (p i ,p j )∈N means (p i ,p j ) are "adjacent" to each other in geometric space; δ(L(p i )≠L(p j ) is a discrete indicator function; when p i With p j When corresponding to different plane labels, L(p i )≠L(p j ) takes 1, when p i With p j When corresponding to the same plane label, L(p i )≠L(p j ) takes 0; P Φ Represents a set of outliers that cannot be properly fitted by any plane; L Φ Represents the outlier label, which does not participate in the plane update and will be penalized if assigned to an outlier; is the outlier penalty coefficient when A smaller value allows more points to be labeled as outliers, while a larger value forces the algorithm to try to find a crack plane for these points. This ensures that points that are too outliers or conflict with all plane directions are managed uniformly.

[0039] P Φ With L ΦThe difference between is that the former represents a set of outliers that cannot be properly fitted by any plane. These are the three-dimensional coordinate points identified as "outliers" during the geometric model fitting process; the latter represents the outlier label, which is a label rather than a point set; this label does not participate in plane updates. If a point is assigned to an outlier label, the point will be penalized (controlled by the penalty coefficient φ).

[0040] δ(L(p i )≠L(p j )) is a simple discrete indicator function, which only focuses on the situation where adjacent points are assigned to different planes; is a composite term. ||p i -L Φ || also represents distance. This "distance" can be understood as a penalty metric used to quantify the distance between points p and p. i The cost of being classified as an outlier.

[0041] During each iteration, the algorithm will try to update which crack plane label each event point should be assigned to; the shape parameters of each crack plane (mainly normal vector and position).

[0042] Label attribution of event points: Based on comprehensive indicators such as the vertical distance from the point to the plane, the degree of fit between the event direction and the plane normal vector (angle), if the new plane model can provide a smaller objective function value, the event point is transferred to the new plane label; if none of them meet the threshold d ini , then its label is set as an outlier.

[0043] Plane model parameters: For each event point assigned to a plane label, a new plane equation is estimated using least squares or weighted least squares (updating the normal vector n and intercept D). The orientation error can be used as a weighting factor to impose a higher penalty on data points that deviate significantly from the normal vector, making the plane more inclined to the main cluster with the most consistent orientation.

[0044] S4. During the iteration of step S3, cluster the directional information of points under the same plane label (e.g., based on angular density clustering), and mark points with significant differences with outlier labels to avoid interfering with plane fitting;

[0045] After each iteration of expansion or re-estimation, the event point assigned to a plane label has its direction (towards θ i , inclination angle Φ i) If the difference is significant, it means that cracks with other dip angles may be mixed in. The density-based spatial clustering algorithm is used to cluster the directions, and the angle threshold ε and the minimum number of points are selected. The angle threshold refers to the difference threshold between strike and dip, and the minimum number of points refers to the minimum number of points that a valid cluster needs to contain. This parameter is used to filter out small-scale point sets that may be noise or outliers, ensuring that only direction sets containing enough points are considered to be valid plane direction categories. The event points assigned to the same plane label are classified according to the difference between strike and dip. The specific process is to first calculate the normal vector direction (strike θ) of each point. i and inclination angle Φ i )

[0046] Mapping direction information into angle space;

[0047] Identify high-density regions in angular space, which represent the main orientation distributions;

[0048] Use the angle threshold ε to divide the dense area and the sparse area;

[0049] Filter noise points according to the minimum number of points m to ensure that each cluster contains enough points;

[0050] Finally, points with similar directions are grouped into the same set to form a directional clustering result. Only the main cluster containing the most event points and the most consistent direction with the current plane normal vector is retained; the remaining event points are transferred to the outlier label P Φ Through this process, points with inconsistent directions are eliminated, reducing the interference of incorrect plane fitting on the overall results.

[0051] To further improve the accuracy of fracture direction estimation, the present invention introduces data cluster analysis of event direction (strike, dip) during the iterative process. Specifically, it includes:

[0052]

[0053] Among them, θp i It represents the direction value of the i-th event assigned to the current crack plane p, C, S are the sine and cosine values, and R is the length of the resultant vector.

[0054] The average direction is calculated as shown in Formula 3:

[0055]

[0056] in, Indicates the average angle of the event points on the plane. S is the average sine value of all angles of the event points on the plane, and C is the average cosine value of all angles of the event points on the plane. It can handle angles in four quadrants, that is, -180° to +180°, and can correctly return the quadrant in which the angle is located.

[0057] Specifically, the calculations of S and C are shown in Formula 4 and Formula 5:

[0058]

[0059] The advantage of this calculation method is that the arithmetic mean of 359° and 1° is 180°, but the actual average direction should be 0°.

[0060] Joint optimization of orientation and geometry

[0061] When the direction of the plane is defined by the normal vector n, and the average event direction is When , the normal vector can be evaluated by The decision of whether to discard the current fitting result of the plane is based on whether the angle between the event point and the plane is too large; or the direction constraint is incorporated into the additional term of the geometric objective function by minimizing the direction deviation (such as 1-cos(Δθ)). Δθ represents the direction of the event point and the average direction of the plane. This angular difference is used to measure the degree of consistency between the direction of the point and the overall direction of the plane. In the objective function, it is converted into a constraint term using the form 1-cos(Δθ). When Δθ is 0 (completely consistent directions), the minimum value is 0; the greater the direction difference, the closer the value is to 2, thus imposing a greater penalty on inconsistent directions.

[0062] The method of incorporating the direction constraint into the additional term of the geometric objective function is to minimize the direction deviation 1-cos(Δθ). Specifically:

[0063] Calculate the direction of each point and the average direction of the plane The angle difference Δθ between them;

[0064] Use the 1-cos(Δθ) function to convert the angle difference into a constraint term;

[0065] Add this constraint term to the geometric objective function as an additional penalty term;

[0066] This can encourage the algorithm to prioritize plane fitting results with higher directional consistency.

[0067] S5. Perform merge re-estimation on plane models whose normal vectors are approximately parallel and spatially close. For those plane models with a large number of points and a high directional concentration in the outlier labels in step S4, regenerate a new plane model. After several iterations, the event points in the outlier labels will be re-generated based on their number, scale, and directional range. For example, if a new plane model is generated using the least squares method, if the matching conditions are still not met, it will be retained in the outlier label and marked as an outlier in the output.

[0068] When a large-scale redistribution of data points occurs, the plane equation (normal vector n and intercept D) needs to be updated with the new set of points.

[0069] The plane can be expressed in the form of Formula 5:

[0070] Ax+By+Cz+D=0 (6)

[0071] Where n = (A, B, C).

[0072] Given a set of points assigned to the plane {p i}, can be obtained by minimizing ||p i -L(p i )|| (i.e. the error from the point to the plane) to solve n and D;

[0073] If the angle between the normal vectors n1 and n2 of two planes (labels) is less than a set threshold (for example, less than 5°), and the closest distance is less than a set threshold (for example, less than 2 meters), they can be considered as the same crack and merged:

[0074] The event point sets of the two planes are combined, and the least squares method is used to re-estimate the new plane equation, and local fine-tuning is performed in subsequent iterations;

[0075] If there are only a few points left in the plane label (for example, less than 5), it may mean that the plane does not exist or the fit is unreliable. This plane label can be discarded, and these points can be transferred to the outlier label or assigned to other planes;

[0076] For the outlier label P Φ If the number of points in the image is, for example, 3 and the direction consistency is high, a new plane label can be generated again to avoid missing cracks.

[0077] S6. If the change in the objective function in step S3 remains small or converges to a stable state within a preset number of iterations, the iterations are terminated and the final multi-fracture plane model and the event information contained in each plane are output. Iterative convergence criteria include the change in the objective function being less than a set threshold, the number of iterations reaching an upper limit, and the event point labels being stable.

[0078] In this embodiment, when the number of iterations reaches 15, it is found that the change in the objective function E is less than 0.5% and the label attribution of the event points is basically stable (on average, only 1-2 points have label changes in each iteration), and the algorithm is determined to have converged.

[0079] Determine the convergence of the objective function E: If the change in E |ΔE| is small enough (less than the convergence threshold) or the change is very weak within the specified maximum number of iterations (for example, 15 times), stop the iteration.

[0080] Output:

[0081] The equation parameters (A, B, C, D) of each fracture plane can be converted into strike and dip (or a description similar to the normal vector n);

[0082] How many event points are included in each plane; how many points are finally included in the outlier label;

[0083] If necessary, the connectivity between planes (if there are merges or intersections) can be provided.

[0084] Determine label changes: If the distribution of event points is stable (label changes are minimal) after several consecutive iterations, it can also be considered convergence.

[0085] Finally complete all processes of the patent for this invention.

[0086] like Figure 2 Shown are multiple crack planes fitted using this method.

[0087] In summary, the technical effects and advantages of the present invention are as follows: the present invention integrates "event position-fracture direction" into the same iterative optimization framework, and realizes the simultaneous extraction of multiple fracture planes by minimizing the objective function containing the geometric model and directional error; this method uses the normal vector and position error of the event point in the discrete fracture network to jointly construct the objective function, and integrates the directional information during the iterative process to eliminate or merge outliers, with strong anti-noise performance and convergence stability; the present invention effectively improves the efficiency and accuracy of automatic extraction of discrete fracture networks, and is suitable for fracture analysis and modeling of complex fault structures or unconventional reservoirs.

[0088] Although the preferred embodiments of the present invention have been described above in conjunction with the accompanying drawings, the present invention is not limited to the above-mentioned specific embodiments. The above-mentioned specific embodiments are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms of specific changes without departing from the scope of protection of the invention and the claims. These all fall within the scope of protection of the present invention.

Claims

1. A method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm, characterized in that: The following steps are involved: S1. Obtain a three-dimensional point cloud of microseismic events from microseismic monitoring data, wherein each microseismic event point includes at least one direction parameter in addition to coordinate information, and the direction parameter includes strike and dip; S2. Use a random sampling method to perform plane fitting on the microseismic event points in step S1 to screen out possible fracture planes. The possible fracture planes and the points they cover are incorporated into the initial plane model set. For each possible fracture plane, its normal vector is calculated and compared with the average direction of the event points assigned to the plane. If the direction difference exceeds a set angle threshold, the possible fracture plane is eliminated. S3. Assign the initial plane in step S2 as a label to each microseismic event point in the 3D point cloud, and iteratively update the objective function based on multi-model energy minimization. Each iteration re-estimates the plane equation parameters and minimizes the objective function based on the distance from the point to the plane and the consistency of the direction. S4. During the iteration of step S3, cluster the directional information of the points under the same plane label and mark the points with significant differences with outlier labels to avoid interfering with plane fitting; S5. Perform merge re-estimation for plane models whose normal vectors are approximately parallel and close in space; regenerate new plane models for those with a large number of points and a high directional concentration in the outlier labels in step S4; S6. When the change in the objective function in step S3 remains small or converges to a stable state within a preset number of iterations, the iteration is stopped, and the final multi-crack plane model and the event information contained in each plane are output.

2. The method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm according to claim 1, characterized in that: The objective function E(ξ) of multi-model energy minimization in step S3 is: Among them, p i Represents the three-dimensional coordinates (x i ,y i ,z i ); L(p i ) indicates p i The plane label of the plane to which it belongs; ||p i -L(p i )|| represents the distance between the event point and the assigned plane; N represents the set of adjacent points, (p i ,p j )∈N means (p i ,p j ) are adjacent to each other in the geometric space; δ(·) represents the indicator function; when p i With p j When corresponding to different plane labels, L(p i )≠L(p j ) takes 1, when p i With p j When corresponding to the same plane label, L(p i )≠L(p j ) is 0; λ is the weight of the smoothing term; P Φ Represents a set of outliers that cannot be properly fitted by any plane; L Φ Represents the outlier label, which does not participate in the plane update and will be penalized if assigned to an outlier; is the outlier penalty coefficient, 3. The method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm according to claim 1, characterized in that: In step S4, the direction information is clustered using a density-based spatial clustering algorithm. The angle threshold ε and the minimum number of points are selected, and the event points assigned to the same plane label are classified according to the differences in strike and dip angles. Only the main cluster with the most consistent direction with the plane normal vector is retained, and the remaining event points are classified as outlier labels.

4. The method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm according to claim 1, characterized in that: In step S5, the plane merging conditions include that the angle between the normal vectors of the two planes is less than 5° and the closest distance between the two planes is less than 2 meters. After merging the event point sets of the two planes, the least squares method is used to re-estimate the new plane.

5. The method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm according to claim 1, characterized in that: In step S5, after several iterations, the event points in the outlier label will try to generate a new plane model again based on their quantity scale and direction range. If the matching conditions are still not met, they will eventually be retained in the outlier label and marked as outliers when output.

6. The method for automatically extracting discrete fracture networks based on a geometric model fitting algorithm according to claim 1, characterized in that: In step S6, the iterative convergence judgment conditions include that the change in the objective function is less than a set threshold, or the number of iterations reaches an upper limit and the event point allocation label tends to be stable.

Images