A comprehensive monitoring and identification method for the propagation range of hydraulic fracturing fractures

Abstract

This application relates to the field of oil and gas field development engineering technology, and discloses a comprehensive monitoring and identification method for the propagation range of hydraulic fracturing fractures, including the following steps: First, a unified three-dimensional coordinate system is established, and microseismic and transient electromagnetic monitoring data are spatially aligned and sliced ​​along the depth direction; then, for each horizontal slice layer, fluid domain polygons based on transient electromagnetic inversion and microseismic event domain polygons based on microseismic location are extracted; subsequently, geometric intersection operations are performed on the fluid domain and microseismic event domain at the same level to identify the overlapping area, which is defined as the effective stimulation range; finally, the effective stimulation volume, fracturing fluid loss rate, and far-field disturbance rate are calculated based on the overlapping area of ​​each layer. This invention achieves accurate quantification and evaluation of the effective stimulation volume of hydraulic fracturing by eliminating ineffective dry fractures and fluid loss zones through geometric topological fusion of multiphysics data.

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development engineering technology, specifically to a method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures. Background Technology

[0002] Hydraulic fracturing technology is a core method for developing unconventional oil and gas resources such as shale gas and tight oil. Its core objective is to construct a complex and conductive fracture network system within the formation to improve reservoir permeability. Accurately monitoring the spatial propagation range of fracturing fractures and scientifically assessing the effective stimulation volume based on this information is of crucial engineering significance for optimizing fracturing operation parameters, predicting oil and gas well productivity, and ensuring operational safety.

[0003] Among existing monitoring technologies, microseismic monitoring is the mainstream method for evaluating fracturing effectiveness. This technology involves deploying geophone arrays on the surface or in wells to receive elastic wave signals generated by rock fractures, and then inverting the source location to infer fracture morphology. However, microseismic monitoring essentially reflects the shear fracture or slip movement of the formation rock. The detected microseismic event points only represent stress release at that location and do not directly indicate the formation of a conductive fracture with a certain width and effectively filled by fracturing fluid and proppant. In actual geomechanical environments, some microseismic events are isolated fractures not connected to the main fracture, or disturbances of far-field dry fractures caused by changes in geostress. Therefore, the stimulation volume calculated solely based on microseismic data often includes a large number of invalid areas, leading to assessment results that are generally overestimated in terms of the actual production contribution area.

[0004] Meanwhile, electromagnetic monitoring technologies such as transient electromagnetic methods utilize the electrical differences between fracturing fluid (usually a low-resistivity fluid) and the surrounding rock medium to track the migration trajectory of underground fluids by observing the secondary field response. Although this method is highly sensitive to fluid distribution and can reflect the affected area of ​​the fracturing fluid, it lacks the ability to directly identify the specific geometric structure and mechanical opening state of fractures. Under underground high pressure, fracturing fluid is prone to microscopic filtration along rock pores or natural weak surfaces, forming a large-scale low-resistivity halo. This filtration range is not equivalent to an effective supporting fracture network. If electromagnetic monitoring results are relied upon alone, ineffective fluid filtration zones can easily be misjudged as effective modification areas, making it impossible to accurately define the physical boundaries of fractures.

[0005] In current engineering practice, although some projects attempt to apply microseismic and electromagnetic monitoring technologies simultaneously, the data processing for both is often independent, or merely involves simple overlaying of results and qualitative comparative analysis. Existing technologies lack a mechanism for rigorous geometric topological fusion and quantitative calculation of heterogeneous monitoring data within a unified three-dimensional coordinate system. This makes it difficult to effectively eliminate interference signals from dry fractures in microseismic monitoring and fluid loss artifacts in electromagnetic monitoring. Consequently, it is impossible to accurately identify the true modified areas where both physical rock fracturing and effective fluid filling have occurred, limiting the accuracy and reliability of fracturing effect evaluation. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides a comprehensive monitoring and identification method for the propagation range of hydraulic fracturing fractures. This method solves the problem that existing single monitoring methods are unable to distinguish between interference from microseismic dry fractures and the artifacts of fluid loss in electromagnetic methods, which leads to the inability to accurately define the true modified volume that has both physical rupture and effective fluid filling.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures, comprising the following steps:

[0008] Step S1: Establish a microseismic monitoring subsystem and a transient electromagnetic monitoring subsystem in the fracturing construction area, establish a unified three-dimensional Cartesian coordinate system, and map the microseismic monitoring data and transient electromagnetic monitoring data into the three-dimensional Cartesian coordinate system to achieve spatial reference alignment of heterogeneous data;

[0009] Step S2: Discretize the monitoring target area along the depth direction into multiple horizontal slice layers, and decompose the irregular crack network in three-dimensional space into geometric features on a two-dimensional plane for processing;

[0010] Step S3: For each horizontal slice layer, based on the apparent resistivity distribution obtained by inversion from transient electromagnetic monitoring data, extract the fluid domain polygon representing the distribution range of fracturing fluid. This polygon characterizes the sweep range of fracturing fluid.

[0011] Step S4: For each horizontal slice layer, based on the source location information in the microseismic monitoring data, extract the microseismic event domain polygon representing the range of rock strata rupture. This polygon characterizes the physical boundary of rock mechanical rupture.

[0012] Step S5: Perform geometric intersection operation on the fluid domain polygon and the microseismic event domain polygon on the same horizontal slice layer to identify the overlapping area of ​​the two, and define the overlapping area as the effective modification range of the horizontal slice layer, thereby eliminating unsupported dry crack areas and ineffective fluid diffusion areas.

[0013] Step S6: Calculate the area of ​​the effective modification range of each horizontal slice layer, and in combination with the thickness of the horizontal slice layer, calculate the effective modification volume of the fracturing operation by integration or summation.

[0014] Step S7: Based on the area of ​​the fluid domain polygon, the area of ​​the microseismic event domain polygon, and the area of ​​the effective modification range, calculate the fracturing fluid loss rate and far-field disturbance rate, and generate comprehensive monitoring and evaluation results.

[0015] Furthermore, in step S2, the system sets the starting and ending depths of the monitoring target area to determine the total thickness, and sets the thickness of each horizontal slice layer. The system calculates the center depth coordinates of each horizontal slice layer on the depth axis, and defines the depth range of the horizontal slice layer based on these center depth coordinates. Subsequent data processing is performed independently on the two-dimensional plane corresponding to the depth range. This layered processing strategy effectively reduces the computational complexity of modeling complex three-dimensional geological bodies.

[0016] Further, the process of extracting the fluid domain polygon in step S3 includes: obtaining the background apparent resistivity matrix of the corresponding horizontal slice layer before fracturing and the monitored apparent resistivity matrix after fracturing, and calculating the resistivity change rate of both. The system uses an adaptive threshold algorithm to determine the apparent resistivity change threshold, selects regions whose resistivity change rate meets the threshold condition and exhibits low resistivity characteristics as low resistivity anomaly regions, and extracts the edge contour of these regions to construct a fluid domain polygon with vertices arranged in a counterclockwise order.

[0017] Furthermore, in step S4, noise reduction processing is required before extracting the polygons. The system filters out microseismic event points whose source depth coordinates fall within the depth range of the current horizontal slice layer, and projects them to form an initial set of microseismic events. Subsequently, a density-based spatial clustering algorithm is used to count the number of neighboring points of each event point within a preset neighborhood radius. If the number of neighboring points is less than the preset minimum number of included points, and the point is not in the neighborhood of any core point, it is marked as a noise point and removed, thereby obtaining an effective set of microseismic points that can truly reflect the main fracture zone.

[0018] Furthermore, in step S4, when constructing the microseismic event domain polygon, the Alpha-Shape algorithm is used to construct a non-convex hull for the effective microseismic point set. This process involves performing Delaunay triangulation on the point set to generate a set of triangles, calculating the circumcircle radius of each triangle, and filtering using a preset rolling sphere radius parameter to remove triangles with excessively large circumcircle radii. Finally, the boundary edges of the remaining triangle set are extracted and connected to form a closed polygon. This method can accurately depict the concave regions between crack branches, avoiding the problem of traditional convex hull algorithms exaggerating the crack extent.

[0019] Furthermore, the geometric intersection operation in step S5 employs the vector cross product algorithm. The system establishes the vector parametric equations of the polygonal edges in the fluid domain, uses the two-dimensional vector cross product operation to solve for the line segment parameter factors, and determines whether the factors are within the valid interval to identify the geometric intersection points. After all valid intersection points are inserted into the original vertex sequence in topological order, the closed path enclosed by the two sets of polygonal boundary segments is extracted using the bidirectional boundary tracing method to accurately define the effective modification range.

[0020] Furthermore, the area calculation in step S6 utilizes the shoelace formula principle, which involves summing the differences in the cross products of the horizontal and vertical coordinates of adjacent vertices in the overlapping plane vertex set and taking the absolute value to obtain the accurate geometric area.

[0021] Furthermore, specific quantitative evaluation indicators are introduced in step S7: the fracturing fluid loss rate is obtained by calculating the ratio of the difference between the total area of ​​the fluid domain polygon and the effective modification area, which is used to evaluate the effective utilization of the fracturing fluid; the far-field disturbance rate is obtained by calculating the ratio of the difference between the total area of ​​the microseismic event domain polygon and the effective modification area, which is used to evaluate the proportion of indirect hydraulic action in the microseismic signal.

[0022] Furthermore, the present invention also includes a three-dimensional reconstruction step, which involves spatially stacking and voxelizing the effective modification range of all horizontal slice layers in a three-dimensional Cartesian coordinate system, assigning different confidence attribute values ​​according to the location of the voxels, and generating a three-dimensional visualization cloud map to intuitively display the spatial distribution characteristics of the cracks.

[0023] Furthermore, this invention supports dynamic evolution monitoring. By dividing the entire fracturing process into continuous time windows, the above steps are repeatedly executed to calculate the effective modified volume curve that changes over time, and the volume growth rate is obtained by differentiating the curve, thereby determining the rate and state of fracture propagation.

[0024] This invention provides a comprehensive monitoring and identification method for the propagation range of hydraulic fracturing fractures. It has the following beneficial effects:

[0025] 1. This invention establishes a unified three-dimensional Cartesian coordinate system and performs geometric intersection operations on the fluid domain polygons obtained from transient electromagnetic inversion and the event domain polygons obtained from microseismic monitoring. This effectively overcomes the limitations of single microseismic monitoring, which cannot distinguish whether a crack is supported, and single electromagnetic monitoring, which is difficult to characterize the fine morphology of cracks. This technical feature can accurately eliminate dry crack areas without fluid filling and filtration areas that have not formed effective support, ensuring that the final determined effective modification area is a substantial production-enhancing area that has both physically fractured and effectively filled with fluid.

[0026] 2. This invention employs a deep layered slicing strategy, combined with density-based clustering noise reduction and the Alpha-Shape non-convex hull algorithm, to independently construct fracture boundaries for each slice. This technique effectively solves the problem that traditional convex hull algorithms easily misjudge blank areas between fracture branches as fracturing regions. It can accurately capture the concave details and non-connected structures in the fracture network, making the 3D reconstruction results more consistent with the objective law of irregular expansion of hydraulic fracturing fractures in geological space.

[0027] 3. Based on the topological differences of polygon sets, this invention innovatively defines and calculates the fracturing fluid loss rate and far-field perturbation rate. By quantifying the area differences between the fluid domain and the effective domain, and between the event domain and the effective domain, the system can intuitively reflect the degree of ineffective diffusion of fracturing fluid and the proportion of indirect hydraulic action in microseismic signals. This helps engineers evaluate the construction effect from the perspective of geometric space utilization and provides a scientific basis for adjusting subsequent fracturing pumping procedures and process parameters. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the system structure of the present invention;

[0029] Figure 2 This is the main flowchart of the method of the present invention;

[0030] Figure 3 This is a schematic diagram illustrating the principle of data fusion and effective modification range identification of the present invention. Detailed Implementation

[0031] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] See attached document Figure 1 This invention provides a comprehensive monitoring and identification system for the propagation range of hydraulic fracturing fractures. The system includes a microseismic monitoring subsystem, a transient electromagnetic monitoring subsystem, and a data processing center. The microseismic monitoring subsystem is communicatively connected to the data processing center and is used to acquire microseismic signals generated by rock fractures. The transient electromagnetic monitoring subsystem is also communicatively connected to the data processing center and is used to acquire electromagnetic signals reflecting changes in the conductivity of the formation.

[0033] The microseismic monitoring subsystem comprises multiple microseismic detectors positioned at predetermined locations on the ground above the fracturing area. These detectors are arranged according to a pre-defined observation system geometry, such as a star-shaped or grid-like arrangement, to cover the target fracturing area. The microseismic detectors are configured to receive P-wave and S-wave signals from the subsurface rock layers and convert analog signals into digital signals for transmission to the data processing center.

[0034] The transient electromagnetic monitoring subsystem includes a transmitting loop and a receiving coil positioned on the ground above the fracturing area. The transmitting loop is configured to transmit a primary pulse magnetic field into the ground, and the receiving coil is configured to observe the induced electromotive force generated by the secondary eddy current field during the gap between the primary pulse magnetic field cutoff. The transient electromagnetic monitoring subsystem obtains formation apparent resistivity information by measuring the decay curve of the induced electromotive force over time.

[0035] The data processing center includes a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, cause the processor to perform the steps of the integrated monitoring and identification method for the propagation range of hydraulic fracturing fractures. The processor first establishes a unified three-dimensional spatial monitoring environment and performs spatial registration and synchronization processing on the data from the microseismic monitoring subsystem and the transient electromagnetic monitoring subsystem.

[0036] Specifically, the processor constructs a unified three-dimensional Cartesian coordinate system. In this coordinate system, The origin is usually set as the projection point of the fracturing wellhead on the ground or the center point of the observation area; shaft and The axes are located in the horizontal plane and point due east and due north, respectively; The axis is perpendicular to the horizontal plane and points towards the depth of the subsurface. The source location coordinates acquired by the microseismic monitoring subsystem... Volume coordinates of resistivity data obtained by inversion from the transient electromagnetic monitoring subsystem All are mapped to this unified three-dimensional Cartesian coordinate system. In order to achieve spatial alignment of heterogeneous data.

[0037] The specific calculation process for coordinate mapping involves the conversion from geographic coordinates to projected coordinates. First, the latitude and longitude coordinates of the fracturing wellhead and each monitoring station are obtained. as high as the earth The Gauss-Kruger projection algorithm is used to convert spherical latitude and longitude into Cartesian coordinates. The projection point of the fracturing wellhead onto the ground is set as the local origin. Then any monitoring point unified coordinates The calculation formula is:

[0038] ;

[0039] in, This is the coordinate rotation angle, used to rotate the coordinate system. Align the axis to the direction of maximum horizontal principal stress or the direction of main crack extension in the construction design; This is the elevation of the origin point; To monitor the vertical depth of the target point relative to the Earth's surface. For microseismic data, the source coordinates... It is obtained through travel-time positioning algorithms (such as the Geiger method) based on velocity model inversion, and its depth reference surface needs to be corrected to... For transient electromagnetic data, the inversion depth is usually calculated downwards from surface measurement points and needs to be corrected accordingly. Using absolute depth coordinates as the reference, thus ensuring that both are in sync. The zero points on the axis are consistent.

[0040] To achieve precise identification of the three-dimensional propagation morphology of fracturing fractures, the processor is configured to monitor the target area in the depth direction. Discretization and layering are performed along the axial direction. The starting depth of the monitoring target area is set as... The termination depth is Total thickness The processor will monitor the target area along... The axis is divided into There are 1 horizontal slice layer, and the thickness of each horizontal slice layer is marked as . .

[0041] No. Depth position of each horizontal slice layer Defined as this layer in The center depth coordinates on the axis, where For layer sequence numbering, For any _th Layer, the corresponding depth range is Subsequent data processing steps, including fluid domain extraction, microseismic event domain extraction, and effective modification range identification, are all performed at each discrete stage. Execute independently on the horizontal slice plane of the layer.

[0042] During the data acquisition phase, the microseismic monitoring subsystem and the transient electromagnetic monitoring subsystem maintain time synchronization under the control of the processor. The specific time synchronization mechanism is as follows: both the microseismic monitoring subsystem and the transient electromagnetic monitoring subsystem are equipped with a high-precision time synchronization module (e.g., a GPS or BeiDou satellite time synchronization module) to receive PPS (pixel per second) signals and UTC time information transmitted by satellites. The processor interacts with each subsystem via the Network Time Protocol (PTP) to periodically correct local clock frequency drift. During acquisition, the sampling rate for microseismic data is typically set to 2kHz to 10kHz, and the sampling rate for transient electromagnetic data is set to 10kHz to 100kHz. To achieve data alignment, the processor operates with microsecond-level precision. Record the timestamp of each data packet. In subsequent data processing, the processor sets a uniform time window based on absolute time, and extracts the micro-seismic event sequence and electromagnetic response sequence within the time window to eliminate timing deviations caused by inconsistent device startup times.

[0043] The processor acquires microseismic event data streams and transient electromagnetic response data streams during the fracturing operation period. For the microseismic data, the processor extracts the three-dimensional coordinates of each microseismic event. For transient electromagnetic data, the processor calculates the apparent resistivity distribution matrix of the subsurface three-dimensional space using an inversion algorithm. These raw data are stored in memory as the basis for subsequent determination of the propagation range of hydraulic fracturing fractures. Through the above system architecture and environment construction, integrated monitoring of microseismic fields and electromagnetic fields with the same spatial and temporal reference is realized.

[0044] After establishing a unified three-dimensional spatial monitoring environment and completing the first... After the layered definition of the horizontal slices, the processor further performs transient electromagnetic data inversion and fluid domain extraction operations to determine the coverage of fracturing fluid underground.

[0045] The processor first retrieves the raw transient electromagnetic response data stored in memory, which includes background field data before fracturing and monitoring field data during fracturing. The processor then uses transient electromagnetic inversion algorithms, such as one-dimensional Occam inversion, two-dimensional pseudo-seismic inversion, or three-dimensional full-space inversion, to convert the time-domain induced electromotive force decay curve into a volume of apparent resistivity distribution data in the three-dimensional subsurface space.

[0046] Before executing the inversion algorithm, the processor first processes the raw induced electromotive force data. Preprocessing is performed. The processor applies digital filters (such as a 50Hz power frequency notch filter) to filter out power line interference and uses a five-point moving average algorithm to smooth the attenuation curve and remove high-frequency random noise. During the inversion process, the processor constructs the objective function. The minimum value of the objective function is found using the iterative least squares method. The objective function is defined as:

[0047] ;

[0048] in, The observed data vector (i.e., the preprocessed induced electromotive force); The vector of model parameters to be solved (i.e., the logarithm of the resistivity of each grid cell in the subsurface); This is the forward response function; The data weighting matrix is ​​typically the reciprocal of the observation error; λ is the model roughness weighting matrix; λ is the regularization factor, used to balance data fit poorness and model smoothness. The processor refines the model parameter vector through multiple iterations. until the objective function The convergence is achieved to a preset accuracy, thereby obtaining a stable three-dimensional apparent resistivity distribution matrix. .

[0049] For each of the first Each horizontal slice layer allows the processor to extract the corresponding planar resistivity matrix from the three-dimensional apparent resistivity distribution data volume. Each element in this matrix represents the depth. Location, corresponding to plane coordinates Apparent resistivity of the rock strata at the location.

[0050] To eliminate the interference of formation background resistivity inhomogeneities on fluid identification, the processor performs background differential processing. The processor calculates the apparent resistivity at the same location after fracturing. Background apparent resistivity before fracturing The difference between them. Based on the physical premise that fracturing fluids typically exhibit low resistivity, the processor sets a threshold for apparent resistivity change. The threshold The determination is based on statistical analysis of the background noise level, and is usually set as a multiple of the standard deviation of the background resistivity. The processor traverses the... Filter all grid points in the layer plane to find those that meet the requirements. The region where resistivity shows a decreasing trend is marked as a low-resistivity anomaly region. This low-resistivity anomaly region physically corresponds to the area swept and filled by fracturing fluid, i.e., the fluid domain.

[0051] After determining the grid range of the low-resistivity anomaly region, the processor uses a contour tracing algorithm or an edge detection algorithm (such as the MarchingSquares algorithm) to extract the outer edge contour of the low-resistivity anomaly region. The processor then vectorizes the extracted closed edge contour and defines it as the first closed edge contour. Fluid domain polygons on layer slices .

[0052] To perform subsequent geometric topology calculations, the processor will process the fluid domain polygons. Indicated as by The set of vertices is an ordered set. According to an embodiment of the invention, this vertex set... The mathematical expression is

[0053] ;

[0054] in, Represents polygons The first on the boundary There are vertices, whose coordinates are represented as follows: The coordinates here Corresponding to the unified Cartesian coordinate system The processor stores the planar coordinates of the polygons. To ensure consistent orientation of the polygons and facilitate subsequent vector cross product calculations, the processor stores the vertex set... The vertices in the set are arranged and indexed in counter-clockwise order. The last vertex in the set... With the first vertex There are connections between them, forming a geometrically closed ring structure that precisely defines the distribution boundary of the fracturing fluid at that depth. If multiple discontinuous low-resistivity anomalous clumps exist on the same slice, the processor defines them as independent sets of sub-polygons. And store them separately.

[0055] See attached document Figure 2 Simultaneously or subsequently, the processor processes the data acquired by the microseismic monitoring subsystem to define the effective areas where physical fracturing of the rock strata occurs. Since hydraulic fracturing fracture networks typically exhibit high irregularity and non-convexity, the processor employs a strategy combining density-based clustering with the Alpha-Shape algorithm to reconstruct the fracture geometry from discrete source points.

[0056] For each of the first In each horizontal slice layer, the processor first performs data mapping and extraction operations. The processor then iterates through all monitored microseismic events, filtering out the focal depth coordinates. Falling within the depth range of this slice layer The processor projects the selected event points onto the plane of the slice layer, forming the initial set of microseismic events corresponding to that layer. According to an embodiment of the present invention, the initial set is composed of... It consists of discrete points, and its mathematical expression is:

[0057] ;

[0058] in, Representing the The microseismic event point is at the [number]th [location]. The projected position on the layer plane, its coordinates are represented as follows: .

[0059] To eliminate interfering points caused by positioning errors, equipment background noise, or isolated stress release events far from the main fracture zone, the processor performs a process on the initial set. Perform density-based spatial clustering noise reduction. The processor pre-sets two key parameters: neighborhood radius. and minimum number of contained points The processor computes any two points in the set. and The Euclidean distance between them. For each point in the set. The processor counts its... The number of other data points contained within a circular neighborhood of radius .

[0060] The processor classifies and distinguishes points based on statistical results: if point of The number of points in the neighborhood is greater than or equal to Then Marked as a core point, indicating the presence of dense fracturing activity at that location; if If it is not a core point and its neighborhood does not contain any core points, then it will be... The data points marked as noise are then processed by the processor from the initial set. Remove the core points and retain the remaining density-reachable boundary points to obtain the purified effective microseismic point set. .

[0061] In obtaining an effective microseismic point set Subsequently, the processor further determines the geometric edge envelope of the point set, i.e., the monitoring range of the microseismic fracture. Given that hydraulic fracturing fractures often exhibit complex branching structures, conventional convex hull algorithms are insufficient to describe the concave regions (blank unruptured areas) between fractures. Therefore, the processor is configured to execute the Alpha-Shape algorithm to construct a non-convex hull.

[0062] Specifically, the processor first processes the effective microseismic point set. Perform Delaunay triangulation on all points to generate a set of triangles covering all points. The processor then calculates each triangle in the set. circumradius .

[0063] The processor introduces the rolling ball radius parameter. This parameter serves as the scale benchmark for spatial filtering. The maximum radius of the allowed blank region is defined. The processor traverses the triangle set and executes... Filtering operation: If a certain triangle The circumcircle radius satisfies If the triangle crosses an excessively large blank area (not a crack area), the processor will remove the triangle from the set. Delete it.

[0064] After filtering, the processor retrieves all edges from the remaining set of triangles and extracts edges shared only by one triangle (i.e., edges not adjacent to other triangles). The processor then connects these edges end-to-end according to geometric adjacency relationships, forming one or more closed polygonal contours, defined as the first polygonal contour. Microseismic event domain polygons on layer slices The polygon It accurately depicts the specific distribution range of rock strata fractures, including details of depressions between fracture branches.

[0065] It is worth noting that in generating polygons During the process, the Alpha-Shape algorithm generates pores within the point set, which are regions located within fracture zones but where microseismic events are relatively sparse. Physically, if these pores are surrounded by a high density of microseismic events, it is generally assumed that the rock mass as a whole has also suffered damage or deformation. Therefore, the processor is configured to perform the pore-filling operation.

[0066] Processor analyzes polygons The processor identifies all inner ring boundaries based on the topological structure of the crack. For inner ring cavities with an area smaller than a preset threshold (e.g., 10% of the total crack area), the processor treats them as connected parts of the crack network and fills them, i.e., it deletes the boundaries constituting the inner ring and merges them into the main polygon region. This step avoids the effective range being artificially fragmented due to local missing monitoring data, and is more consistent with the macroscopic fracture characteristics of rock mechanics.

[0067] The processor will divide the microseismic event domain polygon It is stored as an ordered set of vertices. According to an embodiment of the invention, this vertex set... The mathematical expression is:

[0068] ;

[0069] in, Represents polygons The first on the boundary There are vertices, whose coordinates are represented as follows: The coordinate values ​​are based on the unified Cartesian coordinate system. This set of vertices This will be used as one of the input data for subsequent region overlap calculations.

[0070] In response to the Each horizontal slice layer yielded a fluid domain polygon representing the fluid coverage area. and the microseismic event domain polygon representing the extent of rock strata fracture. The processor then performs a geometric topology fusion calculation to identify the spatially overlapping region between the two. This overlapping region represents the effective modification area at that depth level where both physical fracturing and effective fracturing fluid filling have occurred.

[0071] The processor is configured to perform polygon Boolean intersection operations, i.e., to calculate... To handle fluid domain polygons and microseismic event domain polygon For complex geometric features such as non-convex, multi-branch, and multi-connected domains, the processor uses a general polygon clipping algorithm based on vector algebra (such as an improved version of the Weiler-Atherton algorithm) to solve for accurate results.

[0072] In actual engineering data, polygons and These intersections are often extremely complex, containing self-intersections, overlapping edges, or tiny spikes. To ensure the stability of Boolean operations, the processor performs polygon regularization on the two polygons before performing the intersection calculation. This includes: merging closely spaced adjacent vertices (removing redundant points), eliminating extremely small, elongated, acute triangles (removing spikes), and correcting the clockwise direction of the vertices (unifying them to counterclockwise).

[0073] Furthermore, the processor incorporates anomaly handling mechanisms for degradation scenarios. For example, when the boundaries of the fluid domain and the microseismic domain completely coincide locally (collinear) or only have single-point contact, conventional cross-product algorithms may encounter division-by-zero errors or logical failures. The processor addresses this by introducing small random perturbations (e.g., adding a small amount of perturbation to the coordinate system). The offset (by a factor of magnitude) transforms the degenerate position into a normal position, thus ensuring that the intersection calculation and path tracing algorithms can converge smoothly. This robust design at the numerical computation level guarantees that the system will not crash or enter an infinite loop when processing massive amounts of fracturing monitoring data of various forms.

[0074] Specifically, the processor first performs the intersection point solution step. The processor then traverses the fluid domain polygon. Each edge and microseismic event domain polygon For each edge, calculate all potential geometric intersections. Let the fluid domain polygon... The current edge being processed is composed of vertices. point to The line segment is represented as a vector. Let the microseismic event domain be a polygon. The current edge being processed is composed of vertices. point to The line segment is represented as a vector. Here, and These are the indices of the vertices in their respective sets. The processor constructs the intersection points using vector parametric equations. Computational model:

[0075] ;

[0076] in, For line segments The normalized position parameters on. To solve for the parameters... The processor applies the two-dimensional vector cross product rule.

[0077] Calculate using the following formula:

[0078] ;

[0079] In this embodiment, the symbol × represents the cross product operation of two-dimensional vectors. For any two two-dimensional vectors... and Its cross product is defined as an algebraic sum: .

[0080] The processor calculates the parameters And for line segments Corresponding parameters Perform a validity check. A result is made if and only if the following conditions are met. and At that time, the processor determines the intersection point. It lies on the solid portion of the two line segments and is marked as a valid geometric intersection.

[0081] After acquiring all valid geometric intersections, the processor performs a vertex sequence reassembly step. The processor inserts all valid intersections into the fluid domain polygon according to their geometric positions on the line segments. The original set of vertices and microseismic event domain polygon The original set of vertices In this process, an expanded vertex list containing intersections is formed. The processor further marks each intersection as an in point (entering from the outside to the inside) or an out point (exiting from the inside to the outside) based on the in-and-out relationships of the vertices.

[0082] Subsequently, the processor performs a bidirectional boundary tracing step to construct the overlapping region. The processor then extracts data from the fluid domain polygon. A polygon located in the microseismic event domain Starting from an internal vertex or an inlet point, along the polygon... The processor traces along the boundary direction until it encounters an exit point. At the exit point, the processor switches the tracing path to the microseismic event domain polygon. The boundary, and along The process continues tracking along the boundary direction until it encounters the next ingress point. At the ingress point, the processor switches back to polygons. The processor repeats the switching and tracing process described above until the path closes and returns to the starting point, thus forming a closed sub-region.

[0083] The processor repeats the above tracing process until all unvisited intersections have been processed. All generated closed paths constitute the overlapping plane. The overlapping plane It consists of a single connected region, but also contains multiple discrete sub-regions that are not connected to each other. The processor will overlap the plane. All vertices are stored sequentially to form the final set of vertices within the effective transformation range.

[0084] ;

[0085] in, For the first overlapping region boundary Each vertex has a coordinate system that is consistent with the aforementioned unified Cartesian coordinate system. Maintain consistency.

[0086] As another optional implementation of this embodiment, when the microseismic event domain polygon... When a region or the entire system is determined to be a convex polygon, the processor can be configured to employ the Sutherland-Hodqman clipping algorithm to improve computational efficiency. In this case, the processor will treat the fluid domain polygon as a convex polygon. As the object to be clipped, polygons are used sequentially. The straight line containing each edge is used as the clipping boundary. Perform iterative cutting. For each cutting operation, the processor calculates the intersection points using the aforementioned vector parametric equations and retains those located inside the cutting line (i.e., close to the polygon). The sequence of vertices (on one side of the interior). After processing for polygons... After all edges are clipped, the remaining sequence of vertices forms the overlapping plane. Through the above steps, the system reaches the first... The fusion calculation from heterogeneous monitoring data to precise geometric range was completed on the layer slice.

[0087] See attached document Figure 2 and attached Figure 3 In response to each of the first Each horizontal slice layer completed the effective modification of the overlapping plane. After geometric identification, the processor is configured to perform further numerical calculations to quantify the geometric parameters of the fracturing transformation and reconstruct the three-dimensional spatial entity, thereby outputting a comprehensive evaluation index.

[0088] The processor first targets each individual... Calculate the effective area of ​​a horizontally sliced ​​layer within its planar region. This is done for the set of overlapping plane vertices corresponding to that layer. The processor employs a polygon area calculation algorithm, specifically the shoelace formula, to solve for the geometric area of ​​the closed region. According to an embodiment of the present invention, the first... Effective modification area of ​​layer slices The calculation formula is:

[0089] ;

[0090] In this formula, Represents overlapping planes The total number of vertices; Representing the The x and y coordinates of each vertex, based on the aforementioned unified three-dimensional Cartesian coordinate system. To ensure the closure of the summation calculation, the processor sets boundary conditions, specifying... This associates the last vertex with the first vertex. (Symbol) This represents the absolute value operation to ensure that the calculated area value is always positive regardless of the vertex arrangement direction (clockwise or counterclockwise). If the layer has multiple disconnected overlapping sub-regions, the processor calculates the area of ​​each sub-region separately and sums them to obtain the total effective transformation area of ​​the layer. .

[0091] After acquiring all data within the monitoring target area Effective modification area of ​​each horizontal slice layer Then, the processor performs 3D spatial reconstruction and volume calculation. The processor uses integration or discrete accumulation methods to calculate the area of ​​the 2D slice along the depth axis. Stacking is performed to calculate the effective modified volume (SRV) of the entire fracturing zone. According to an embodiment of the invention, the effective modified volume... The calculation model is as follows:

[0092] ;

[0093] in, The thickness of a single horizontal slice layer, i.e., the preset discretization step size; This represents the total number of slice layers. This volume... Physically, it represents the total volume of rock mass where the rock strata have both fractured and been effectively supported by fracturing fluid, and is a key quantitative indicator for evaluating the potential for increased fracturing production. The processor also interpolates and connects the vertex data of the overlapping planes of each layer in three-dimensional space to construct a three-dimensional geological model of the effective modification range, and then visualizes and renders it through a display device.

[0094] In terms of 3D visualization, the processor goes beyond simple geometric stacking; instead, it employs voxelization technology to construct high-precision attribute models. The processor divides the underground space into tiny cubic units (voxels), each of which is assigned a transformation confidence attribute.

[0095] For those located in overlapping regions Voxels within the domain are assigned high confidence values ​​(e.g., 1.0); for those located only in the fluid domain... Or only located in the microseismic zone The voxels are assigned medium confidence values ​​(e.g., 0.5), while external regions are assigned zero values. The processor then applies a color map, rendering high-confidence areas in warm tones (e.g., red) and low-confidence areas in cool tones (e.g., blue), and sets transparency gradients. The resulting 3D cloud map visually demonstrates the spatial hierarchy between the core area of ​​the fracturing fracture and the affected area of ​​secondary fractures, helping engineers identify key engineering risks such as asymmetric fracture propagation, abnormal breakthroughs in overlying strata, or ineffective extension into the floor.

[0096] To further evaluate the process compatibility of fracturing operations and the composition characteristics of monitoring data, the processor constructs dimensionless evaluation indicators based on geometric set relationships, including fracturing fluid loss rate. and far-field perturbation rate .

[0097] The processor calculates the fracturing fluid filtration rate. This characterizes the proportion of fracturing fluid injected into the formation that fails to form effective supporting fractures, instead filtering out through rock pores or ineffectively diffusing along natural weak surfaces. The processor acquires the fluid domain polygon. Total area And calculate according to the following formula:

[0098] ;

[0099] in, This refers to the effective renovation area at the corresponding level or overall. If... A high value indicates a significant loss of fracturing fluid, requiring adjustment of the fracturing fluid viscosity or the discharge rate.

[0100] At the same time, the processor calculates the far-field perturbation rate. This is used to characterize the proportion of far-field stress release or noise signals that are not directly caused by hydraulic forces in rupture events detected by microseismic monitoring. The processor acquires the microseismic event domain polygon. Total area And calculate according to the following formula:

[0101] ;

[0102] in, This represents the area of ​​the rupture zone delineated solely through microseismic monitoring. This indicator quantifies the error redundancy of a single microseismic monitoring method under the current geological environment, verifying the corrective effect of the integrated monitoring method of this invention.

[0103] Finally, the processor will calculate the effective area of ​​a single layer. Total effective volume of renovation Fracturing fluid filtration loss rate and far-field perturbation rate The system compiles and generates a comprehensive evaluation report on fracturing performance, which is then stored in a memory or transmitted to an external terminal via a communication interface, providing accurate data support for subsequent optimization of fracturing process parameters and production capacity prediction.

[0104] Furthermore, this embodiment also supports process evaluation based on time-volume curves. The processor applies steps S1 to S6 in each time window. The curve of effective modification volume changing over time was calculated above. The processor differentiates the curve to obtain the volume growth rate. .

[0105] This dynamic indicator has significant engineering guidance value: in the early stages of fracturing, the growth rate is usually high, indicating that the main fracture is forming rapidly; as fracturing progresses, if the growth rate suddenly drops to near zero, it indicates that fracture propagation is hindered or has reached its limit; if the growth rate is high but the filtration rate is low... A synchronized surge in fracturing fluid indicates that the fracturing fluid has connected to a natural fault or a highly permeable channel. The system correlates these dynamic characteristics with real-time pumping procedures (displacement rate, sand ratio) to generate real-time pump shutdown or fracturing rerouting suggestions, thus achieving a technological leap from post-event assessment to real-time optimization. This is not only about monitoring fracturing effectiveness but also about proactively controlling and assisting fracturing operations.

Claims

1. A method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures, characterized in that, Includes the following steps: S1. Establish a microseismic monitoring subsystem and a transient electromagnetic monitoring subsystem in the fracturing construction area, establish a unified three-dimensional Cartesian coordinate system, and map the microseismic monitoring data and transient electromagnetic monitoring data into the three-dimensional Cartesian coordinate system; S2. Discretize the monitoring target area into multiple horizontal slice layers along the depth direction; S3. For each horizontal slice layer, extract the fluid domain polygon representing the distribution range of fracturing fluid based on the apparent resistivity distribution obtained from the transient electromagnetic monitoring data. S4. For each horizontal slice layer, based on the source location information in the microseismic monitoring data, extract the microseismic event domain polygon representing the extent of rock rupture. S5. Perform geometric intersection operation on the fluid domain polygon and the microseismic event domain polygon on the same horizontal slice layer, identify the overlapping area of ​​the fluid domain polygon and the microseismic event domain polygon, and define the overlapping area as the effective modification range of the horizontal slice layer; S6. Calculate the area of ​​the effective modification range of each horizontal slice layer, and calculate the effective modification volume of the fracturing operation in combination with the thickness of the horizontal slice layer. S7. Based on the area of ​​the fluid domain polygon, the area of ​​the microseismic event domain polygon, and the area of ​​the effective modification range, calculate the fracturing fluid loss rate and far-field disturbance rate, and generate comprehensive monitoring and evaluation results, specifically including: Obtain the total area of ​​the fluid domain polygon, the total area of ​​the microseismic event domain polygon, and the effective modification area of ​​the effective modification range; Calculate the difference between the total area of ​​the fluid domain polygon and the effective modified area, and divide the difference by the total area of ​​the fluid domain polygon to obtain the fracturing fluid loss rate; The difference between the total area of ​​the microseismic event domain polygon and the effective modified area is calculated, and the difference is divided by the total area of ​​the microseismic event domain polygon to obtain the far-field perturbation rate.

2. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, Step S2, which discretizes the monitoring target area into multiple horizontal slices along the depth direction, specifically includes: Set the starting and ending depths of the monitoring target area and determine the total thickness of the monitoring target area; Set the thickness of a single horizontal slice layer; Calculate the center depth coordinates of each horizontal slice layer on the depth axis, and define the depth range of the horizontal slice layer based on the center depth coordinates. The subsequent processing steps S3 to S5 are all performed independently on the two-dimensional plane corresponding to the depth range.

3. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, Step S3 specifically includes: Obtain the background apparent resistivity matrix of the corresponding horizontal slice layer before fracturing and the monitored apparent resistivity matrix after fracturing; Calculate the rate of change of resistivity of the monitored apparent resistivity matrix relative to the background apparent resistivity matrix; The apparent resistivity change threshold is determined using an adaptive threshold algorithm; Regions whose resistivity change rate meets the apparent resistivity change threshold and exhibit low resistivity characteristics are selected as low resistivity anomaly regions. Extract the edge contour of the low-resistivity anomaly region and construct the fluid domain polygon with vertices arranged in a counterclockwise order.

4. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, In step S4, before extracting the microseismic event domain polygon representing the extent of rock strata fracture, a step of denoising the microseismic monitoring data is also included: Microseismic event points whose focal depth coordinates fall within the current horizontal slice layer depth range are selected, and these microseismic event points are projected onto the horizontal slice layer plane to form an initial set of microseismic events. The initial set of microseismic events was processed using a density-based spatial clustering algorithm. Count the number of neighboring points within a preset neighborhood radius for each microseismic event point; If the number of adjacent points is less than the preset minimum number of included points, and the microseismic event point is not in the neighborhood of any core point, then the microseismic event point is marked as a noise point and removed from the initial microseismic event set to obtain an effective microseismic point set.

5. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 4, characterized in that, Step S4, specifically the step of extracting the microseismic event domain polygon representing the extent of rock strata fracture, includes: The Alpha-Shape algorithm is used to construct a non-convex envelope for the effective microseismic point set; Perform Delaunay triangulation on the effective microseismic point set to generate a set of triangles; Calculate the circumcircle radius of each triangle in the triangle set; The set of triangles is filtered using a preset rolling ball radius parameter to remove triangles whose circumscribed circle radius is greater than the rolling ball radius parameter; Extract the boundary edges of the remaining triangle set and connect them to form a closed polygon representing the microseismic event domain.

6. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, Step S5, specifically the steps for performing the geometric intersection operation, include: The geometric intersection points of the edges of the fluid domain polygon and the edges of the microseismic event domain polygon are solved using the vector cross product algorithm. Establish the vector parametric equations of the edges of the fluid domain polygon, and solve for the line segment parameter factors using the two-dimensional vector cross product operation; Determine whether the line segment parameter factor is within the valid range. If it is within the valid range, then it is determined that there is a valid intersection point. All valid intersection points are inserted into the original vertex sequence in topological order, and the closed path formed by the polygonal boundary segments of the fluid domain and the polygonal boundary segments of the microseismic event domain is extracted using the bidirectional boundary tracing method. The closed path is the effective modification range.

7. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, Step S6, specifically calculating the area of ​​the effective modification range of each horizontal slice layer, includes: Obtain the set of overlapping plane vertices within the effective modification range; Using the shoelace formula, the difference between the cross products of the horizontal and vertical coordinates of adjacent vertices in the overlapping plane vertex set is summed and the absolute value is taken to obtain the effective modification area of ​​the horizontal slice layer.

8. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, It also includes a three-dimensional reconstruction step of the effective modification range: Based on a three-dimensional Cartesian coordinate system, the effective modification range of all horizontal slice layers is spatially stacked. The stacked spatial data is voxelized, and different confidence attribute values ​​are assigned to voxels located within the effective modification range, voxels located only within the fluid domain polygon, and voxels located only within the microseismic event domain polygon. A 3D visualization cloud map is generated based on the confidence attribute value.

9. The method for comprehensive monitoring and identification of the propagation range of hydraulic fracturing fractures according to claim 1, characterized in that, It also includes dynamic evolution monitoring steps: The entire fracturing process is divided into multiple consecutive time windows; For each time window, acquire the microseismic monitoring data and transient electromagnetic monitoring data within that time window; Repeat steps S3 to S6 to calculate the effective modified volume curve as a function of time. Differentiating the effective modified volume curve yields the volume growth rate, which is used to determine the crack propagation state.

Images