Fracture network morphology data analysis method and system for oil and gas fields

Abstract

The present application provides a kind of oil and gas field fracture network form data analysis method and system, it is related to computer technology field, first acquisition target fracturing well area after fracturing operation downhole microseismic monitoring data set and interwell tracer monitoring data set;Wellbore microseismic monitoring data set is carried out spatial clustering processing, constructs initial fracture network spatial skeleton;According to initial fracture network spatial skeleton, interwell tracer monitoring data set is carried out migration path matching, generates fracture network connectivity topology;Based on fracture network connectivity topology, initial fracture network spatial skeleton is supplemented, and complete fracture network three-dimensional form is generated;Finally, fracture network parameter is extracted, and fracture network form data set is obtained.The present application can comprehensively and accurately analyze fracture network form, which helps to improve oil and gas recovery and reduce development cost.

Description

Technical Field

[0001] This invention relates to the field of computer technology, and more specifically, to a method and system for analyzing the morphology of fracture networks in oil and gas fields. Background Technology

[0002] In the field of oil and gas field development, fracturing is one of the key technologies for improving oil and gas recovery. Accurately understanding the morphology of fracturing fracture networks is crucial for evaluating fracturing effectiveness and optimizing subsequent development strategies. Currently, traditional methods for analyzing fracturing fracture network morphology have many limitations.

[0003] On the one hand, some methods rely solely on single monitoring data, such as downhole microseismic monitoring data. While these methods can obtain information such as the spatial location and timing of microseismic events, they cannot comprehensively reflect the actual fluid transport within the fracture network. This is because microseismic events primarily reflect the energy release caused by rock fracturing and cannot directly demonstrate the flow path and connectivity of tracers within the fractures.

[0004] On the other hand, some methods lack consideration for the integrity of the fracture network when constructing the fracture network morphology. The initial fracture skeleton constructed solely based on microseismic data may be incomplete, failing to accurately represent the complex relationship between the main fracture channels and secondary connected fracture channels. This results in an incomplete and inaccurate understanding of the fracture network morphology, which in turn affects the assessment of oil and gas field development effectiveness and the formulation of subsequent development plans. Summary of the Invention

[0005] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a method for analyzing fracture network morphology data in oil and gas fields, the method comprising: Acquire a set of downhole microseismic monitoring data and an inter-well tracer monitoring dataset collected after fracturing operations in the target fracturing well area. The downhole microseismic monitoring data set includes the spatial coordinates of multiple microseismic event points recorded in chronological order, as well as the earthquake time parameters and magnitude and energy parameters corresponding to each microseismic event point. The inter-well tracer monitoring dataset includes tracer breakthrough time records and tracer concentration change curves between the injection well and the production well. The downhole microseismic monitoring data set is subjected to spatial clustering of microseismic event points. Based on the spatial proximity between the spatial coordinates of the microseismic event points and the temporal order between the earthquake occurrence time parameters, an initial spatial skeleton of the hydraulic fracture network containing multiple spatial clusters of microseismic event points is constructed. The initial spatial skeleton of the hydraulic fracture network is used to characterize the orientation and extension trend of the main fracture channels formed by hydraulic fracturing in three-dimensional space and their branching and converging relationships. Based on the initial fracture network spatial skeleton, tracer migration paths are matched on the inter-well tracer monitoring dataset. The spatial location and connectivity direction of the dominant tracer migration channels in the main fracture channels are determined by combining the tracer breakthrough time record and the tracer concentration change curve, and a fracture network connectivity topology containing the connectivity relationships of the dominant migration channels is generated. Based on the fracture network connectivity topology, the initial pressure fracture network spatial skeleton is supplemented with fracture network integrity. The spatial location corresponding to the dominant transport channel is spatially associated and fused with the main fracture channel in the initial pressure fracture network spatial skeleton to generate a complete three-dimensional shape of the pressure fracture network containing the main fracture channel and the secondary connected fracture channel. The fracture network parameters are extracted from the complete three-dimensional morphology of the fracture network to obtain the fracture network morphology data set of the target fracture well area. The fracture network morphology data set includes the total surface area parameter of the fracture network, the total volume parameter of the fracture network, the branch density parameter of the fracture network, and the tortuosity parameter of the fracture network.

[0006] Furthermore, embodiments of the present invention also provide an oil and gas field fracture network morphology data analysis system, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to perform the above-described oil and gas field fracture network morphology data analysis method by executing the machine-executable instructions.

[0007] In another aspect, embodiments of the present invention also provide a computer program product, the computer program product including machine-executable instructions, the machine-executable instructions being stored in a computer-readable storage medium, the processor of the oil and gas field fracture network morphology data analysis system reading the machine-executable instructions from the computer-readable storage medium, the processor executing the machine-executable instructions, causing the oil and gas field fracture network morphology data analysis system to perform the above-described oil and gas field fracture network morphology data analysis method.

[0008] Based on the above, by performing spatial clustering of microseismic event points on the downhole microseismic monitoring data set, an initial spatial framework for the fracturing network is constructed. This framework can present the direction and extension trend of the main fracturing channels in three-dimensional space, as well as their branching and converging relationships. Secondly, based on the initial spatial framework, tracer migration paths are matched to the inter-well tracer monitoring dataset to generate a fracture network connectivity topology. This determines the spatial location and connectivity direction of the dominant tracer migration channels in the main fracture channels, compensating for the inability of microseismic data alone to reflect fluid migration. Then, based on the fracture network connectivity topology, the initial spatial framework for the fracturing network is supplemented with fracture network integrity data, generating a complete three-dimensional morphology of the fracturing network including main fracture channels and secondary connected fracture channels. Finally, fracture network parameters are extracted from the complete three-dimensional morphology of the fracturing network. The resulting fracturing network morphology data set can provide a valuable reference for evaluating oil and gas field development effectiveness and optimizing subsequent development plans, contributing to improved oil and gas recovery and reduced development costs. Attached Figure Description

[0009] Figure 1 This is a schematic diagram of the execution flow of the oil and gas field fracture network morphology data analysis method provided in the embodiments of the present invention.

[0010] Figure 2 This is a schematic diagram of exemplary hardware and software components of the oil and gas field fracture network morphology data analysis system provided in an embodiment of the present invention. Detailed Implementation

[0011] Figure 1 This is a flowchart illustrating a method for analyzing the morphology of fracture networks in oil and gas fields according to an embodiment of the present invention, which will be described in detail below.

[0012] Step S110: Obtain the downhole microseismic monitoring data set and the inter-well tracer monitoring data set collected after the fracturing operation in the target fracturing well area. The downhole microseismic monitoring data set includes the spatial coordinates of multiple microseismic event points recorded in chronological order, as well as the earthquake time parameters and magnitude and energy parameters corresponding to each microseismic event point. The inter-well tracer monitoring data set includes the tracer breakthrough time record and tracer concentration change curve between the injection well and the production well.

[0013] In this embodiment, for a target fracturing well area, multiple injection wells for fracturing operations and multiple production wells for extraction are deployed in the well area. After the hydraulic fracturing operation is completed, two types of key monitoring data need to be obtained to support the subsequent fracture network analysis.

[0014] The first type of data is a set of downhole microseismic monitoring data. Its generation process is as follows: Data acquisition trigger commands are sent to multiple downhole microseismic geophone arrays deployed in monitoring wells around the target fracturing well area. The raw microseismic waveform data, containing precise timestamps, is received from these geophone arrays. Then, the raw microseismic waveform data undergoes automatic microseismic event identification, detecting waveform anomaly intervals where the waveform amplitude exceeds a preset environmental noise limit. The start and end times of each anomaly interval are used as time windows for candidate microseismic events. Next, the waveform data within each candidate microseismic event time window is used to pick up the first arrival of the seismic phase. The energy ratio curve is calculated using the long-short time window average ratio method, and the time point corresponding to the first maximum value of the curve is taken as the first arrival time of the P-wave. The time point corresponding to the second maximum value is taken as the first arrival time of the shear wave. Then, based on the first arrival times of the p-wave and shear wave, and combined with the spatial coordinates of multiple geophone arrays in the geological coordinate system, a spatial positioning algorithm based on p-wave and shear wave velocity models is used to calculate the spatial coordinates of the candidate microseismic event. At the same time, the magnitude energy parameter is calculated based on the maximum amplitude of the original microseismic waveform data within the event's time window, and the first arrival time of the p-wave is used as the earthquake occurrence time parameter to generate a microseismic event point record containing the spatial coordinates of the microseismic event point, the earthquake occurrence time parameter, and the magnitude energy parameter. Finally, the microseismic event point records corresponding to all identified candidate microseismic events are sorted and combined according to the order of the earthquake occurrence time parameter to generate a downhole microseismic monitoring data set.

[0015] The second type of data is the inter-well tracer monitoring dataset, which is generated as follows: tracer solution is injected into the tracer injection wellhead in the target fracturing well area, and online tracer concentration monitoring instruments are deployed at multiple production wellheads to receive tracer concentration sampling data returned by these instruments at continuous time points; then, time series analysis is performed on the tracer concentration sampling data returned by each production wellhead to detect the sampling time point when the concentration value first exceeds the pre-configured background concentration threshold, and this sampling time point is taken as the tracer breakthrough time record between the production well and the tracer injection well; next, tracer concentration sampling data returned by each production wellhead at all time points after the tracer breakthrough time record are extracted, and the sampling time points and corresponding concentration values ​​are arranged in chronological order to generate the tracer concentration change curve of the production wellhead; finally, the tracer breakthrough time records between all injection wells and all production wells and the tracer concentration change curves of all production wellheads are correlated and combined to generate the inter-well tracer monitoring dataset.

[0016] Step S120: Perform spatial clustering processing on the downhole microseismic monitoring data set for microseismic event points. Based on the spatial proximity relationship between the spatial coordinates of the microseismic event points and the temporal order between the earthquake occurrence time parameters, construct an initial spatial skeleton of the hydraulic fracture network containing multiple spatial clusters of microseismic event points. The initial spatial skeleton of the hydraulic fracture network is used to characterize the orientation and extension trend and branching and converging relationship of the main fracture channels formed by hydraulic fracturing in three-dimensional space.

[0017] After obtaining the set of downhole microseismic monitoring data, the discrete event points need to be processed to construct a preliminary fracture network framework.

[0018] Step S121: Perform spatial coordinate analysis on each microseismic event point in the downhole microseismic monitoring data set, extract the three-dimensional spatial coordinate values ​​corresponding to each microseismic event point and the offset of the coordinate origin of the three-dimensional spatial coordinate values ​​in the geological coordinate system to which the target fracturing well area belongs, and use the extracted three-dimensional spatial coordinate values ​​of each microseismic event point as the spatial location identifier of the microseismic event point.

[0019] First, each event point record is traversed from the microseismic monitoring dataset. For each record, the values ​​of the three coordinate axes are parsed out, for example, denoted as the x-coordinate X_e, y-coordinate Y_e, and z-coordinate Z_e, respectively. Simultaneously, the origin offset of these coordinate values ​​is recorded. This origin offset describes the translation relationship of the geological coordinate system relative to a certain absolute geographic coordinate system. However, subsequent spatial calculations are all performed within the geological coordinate system based on this offset. Therefore, the three-dimensional spatial coordinate values ​​(X_e, Y_e, Z_e) of each event point are themselves its unique spatial identifier, and these coordinate values ​​will be used for all subsequent spatial distance calculations.

[0020] Step S122: Perform time parameter analysis processing on each microseismic event point in the downhole microseismic monitoring data set, extract the absolute time value and the relative time difference between the absolute time value and the start time of the fracturing operation from the earthquake time parameter corresponding to each microseismic event point, and use the extracted absolute time value of each microseismic event point as the time sequence identifier of the microseismic event point.

[0021] Next, for the same microseismic event record, the seismic occurrence time parameter T_p is analyzed. This seismic occurrence time parameter is usually an absolute time value accurate to the second or even millisecond, for example, in the format YYYY-MM-DDHH:MM:SS.sss. To facilitate the analysis of the temporal sequence of fracture propagation, the difference between this absolute time value T_p and the official start time of the fracturing operation T_frack_start is also calculated, resulting in a relative time difference Δt_rel, in seconds or minutes. The formula is Δt_rel equals T_p minus T_frack_start. This relative time difference can be used to determine the order in which ruptures occur, but as a time sequence identifier, the unique absolute time value T_p is still the standard.

[0022] Step S123: Calculate the spatial Euclidean distance between any two microseismic event points based on the spatial location identifier of each microseismic event point, and calculate the absolute value of the time interval between any two microseismic event points based on the time sequence identifier of each microseismic event point, thereby generating a microseismic event spatiotemporal correlation matrix containing a spatial Euclidean distance numerical matrix and a time interval absolute value matrix.

[0023] All microseismic event points in the dataset are paired. Assume there are N event points in the dataset. For any two event points, denoted as event point i and event point j, event point i has three-dimensional spatial coordinates (X_i, Y_i, Z_i) and an absolute time value T_i, and event point j has three-dimensional spatial coordinates (X_j, Y_j, Z_j) and an absolute time value T_j. First, calculate their spatial Euclidean distance D_ij, calculated as the square root of ((X_i minus X_j) squared plus (Y_i minus Y_j) squared plus (Z_i minus Z_j) squared). Second, calculate their absolute time interval Δt_ij, which is the absolute value of T_i minus T_j. Arrange all calculated D_ij values ​​in N rows and N columns according to the event point index i, forming an N x N matrix M_distance, where the diagonal element D_ii is zero. Similarly, all calculated Δt_ij are arranged into an N x N matrix M_time. These two matrices, M_distance and M_time, together constitute the spatiotemporal correlation matrix of microseismic events in this dataset.

[0024] Step S124: Call the pre-configured spatial clustering density threshold parameter to perform spatial proximity filtering on the spatiotemporal correlation matrix of the microseismic events, and mark the microseismic event point pairs whose spatial Euclidean distance values ​​are less than the spatial clustering density threshold parameter as spatial proximity event point pairs, thereby obtaining an initial spatial proximity relationship set containing all spatial proximity event point pairs.

[0025] A spatial clustering density threshold parameter, denoted as D_threshold, is preset. This threshold parameter is set based on the rock mechanical properties and fracturing operation parameters of the well area, representing the maximum distance between two event points that can be considered to belong to the same fracture channel in space. The spatial Euclidean distance matrix M_distance is traversed. For each pair of event points (i, j), if the spatial Euclidean distance D_ij between them is less than the spatial clustering density threshold parameter D_threshold, then the point pair (i, j) is marked as a spatially neighboring event point pair, and the identifiers i and j of these two event points are recorded. All spatially neighboring event point pairs that meet the condition form an initial spatial proximity relationship set Set_proximal. This initial spatial proximity relationship set can be stored as a list, where each element is an unordered pair (i, j).

[0026] Step S125: Perform time continuity verification on each spatial neighbor event point pair in the initial spatial neighbor relationship set, extract the absolute value of the time interval corresponding to each spatial neighbor event point pair, and retain spatial neighbor event point pairs with the absolute value of the time interval less than the pre-configured time continuity threshold parameter as valid neighbor event point pairs with spatiotemporal consistency, and generate a spatiotemporally consistent neighbor relationship set.

[0027] Spatially proximate event points should also be temporally continuous to belong to the same expanding crack. Therefore, for each spatially proximate event point pair (i, j) in the initial spatial proximate relationship set Set_proximal, their corresponding absolute time interval value Δt_ij is retrieved from the absolute time interval value matrix M_time. A preset temporal continuity threshold parameter, denoted as Δt_threshold, is used, representing the maximum allowable time interval between two adjacent rupture events during crack expansion. If the Δt_ij of the point pair is less than the temporal continuity threshold parameter Δt_threshold, then the two events are considered not only spatially proximate but also temporally continuous, and are considered valid. The point pair (i, j) is retained and added to the new set Set_temporal_consistent. Conversely, if Δt_ij is greater than or equal to Δt_threshold, then they are considered to possibly belong to ruptures of different phases or have no direct temporal correlation, and are removed from the set. After this filtering step, the set of spatiotemporal consistent proximity relationships, Set_temporal_consistent, is obtained.

[0028] Step S126: Construct an undirected graph structure based on the effective neighbor event point pairs in the spatiotemporally consistent neighbor relationship set. The undirected graph structure uses each microseismic event point as a graph node and the spatial neighbor relationship corresponding to the effective neighbor event point pair as the connecting edge between the graph nodes to generate an initial event point association graph representing the spatiotemporal association relationship between microseismic event points.

[0029] Each microseismic event point is abstracted as a graph node, with a total of N nodes. For each valid neighboring event point pair (i, j) in the spatiotemporal consistent proximity set `Set_temporal_consistent`, an undirected edge is established between their corresponding graph nodes i and j. This edge represents a direct, spatiotemporally consistent association between nodes i and j. All point pairs in the set are traversed, and edges are added between the corresponding nodes in turn. Finally, all N nodes and the edges between them together form an undirected graph structure, namely the initial event point association graph `Graph_events`. In this graph, connected nodes represent spatiotemporally consistent rupture event groups.

[0030] Step S127: Perform connectivity component detection processing on the initial event point association graph, identify all interconnected graph node groups in the initial event point association graph, and treat each interconnected graph node group as an independent microseismic event point spatial cluster, thereby obtaining an initial cluster set containing multiple microseismic event point spatial clusters.

[0031] A connected component detection algorithm from graph theory, such as depth-first search or breadth-first search, is used to analyze the initial event point graph Graph_events constructed in step S126. Starting from any unvisited node, the graph traverses along the edges, and all reachable nodes constitute a connected component. Once all nodes in a connected component have been visited, the search continues for the next unvisited node, repeating the process until all nodes in the graph are assigned to a connected component. Each identified connected component is denoted as Cluster_k, and the microseismic event points represented by all the graph nodes within it constitute a spatial cluster of microseismic event points. Thus, the entire initial event point graph Graph_events is divided into several independent clusters, and the set of these clusters is denoted as Set_clusters_initial, which is the initial cluster set.

[0032] Step S128: Based on the initial cluster set, perform cluster feature extraction and spatial connectivity analysis to construct the initial spatial skeleton of the hydraulic fracture network.

[0033] After obtaining the initial cluster set Set_clusters_initial, it is necessary to further analyze the characteristics of each cluster and the relationship between clusters in order to construct a continuous fracture network skeleton.

[0034] Step S1281: Extract the three-dimensional spatial coordinate values ​​of all microseismic event points in the spatial cluster of each microseismic event point, calculate the average value of the three-dimensional spatial coordinate values ​​of all microseismic event points in the spatial cluster of each microseismic event point as the spatial coordinates of the cluster center of the spatial cluster of the microseismic event point, and extract the earliest and latest earthquake occurrence times from the earthquake occurrence time parameters of all microseismic event points in the spatial cluster of each microseismic event point, and calculate the time span parameter of the spatial cluster of the microseismic event point.

[0035] For each cluster in the initial set of clusters, Set_clusters_initial, such as Cluster_k, all microseismic event points within the cluster are traversed. First, the x-coordinate values ​​of all event points are summed to obtain the x-coordinate sum Sum_X, the y-coordinate values ​​to obtain the y-coordinate sum Sum_Y, and the z-coordinate values ​​to obtain the z-coordinate sum Sum_Z. These are then divided by the total number of event points N_k within the cluster to obtain the spatial coordinates of the cluster center (Center_X_k, Center_Y_k, Center_Z_k). Second, the minimum value of the earthquake occurrence time parameter T_p among all event points within the cluster is identified as the earliest occurrence time T_earliest_k, and the maximum value is identified as the latest occurrence time T_latest_k. Subtracting T_earliest_k from T_latest_k yields the time span parameter Δt_span_k for the cluster, which reflects the duration of the crack channel formation.

[0036] Step S1282: Based on the spatial coordinates of the cluster center of each microseismic event point spatial cluster and the three-dimensional spatial coordinates of all microseismic event points in the spatial cluster, the principal component analysis algorithm is used to fit the spatial principal axis direction of each microseismic event point spatial cluster. The covariance matrix of the three-dimensional spatial coordinates of all microseismic event points in the spatial cluster is calculated. The covariance matrix is ​​decomposed into eigenvalues ​​to obtain the eigenvector corresponding to the largest eigenvalue. The spatial direction of the eigenvector corresponding to the largest eigenvalue is taken as the spatial orientation direction of the crack segment represented by the spatial cluster of the microseismic event point.

[0037] To determine the extension direction of the crack segment represented by each cluster, principal component analysis is performed on the coordinates of all event points within each cluster. First, an N_k x 3 matrix M_coords is constructed using the coordinate data of all event points within cluster_k. Each row represents an event point, and the three columns represent the x-coordinate, y-coordinate, and y-coordinate, respectively. The covariance matrix Cov_k of this matrix is ​​calculated. This covariance matrix is ​​a 3x3 matrix that describes the dispersion of the coordinate data in each direction. The covariance matrix is ​​calculated by first obtaining the mean value for each coordinate axis (obtained in step S1281), then calculating the deviation of each event point's coordinates from the mean, and finally calculating the covariance between these deviations. Then, the covariance matrix Cov_k is decomposed into eigenvalues ​​to obtain three eigenvalues ​​λ1_k, λ2_k, and λ3_k, where λ1_k is greater than or equal to λ2_k and greater than or equal to λ3_k, and the corresponding three eigenvectors v1_k, v2_k, and v3_k. The eigenvector v1_k corresponding to the largest eigenvalue λ1_k represents the direction of the greatest data dispersion, which is the direction in which the event points within the cluster are most distributed. The spatial orientation of this eigenvector v1_k, i.e., its directional components (v1x_k, v1y_k, v1z_k) on the three coordinate axes, is taken as the spatial orientation of the crack segment represented by the cluster Cluster_k.

[0038] Step S1283: Based on the spatial coordinates of the cluster center, spatial orientation direction, and time span parameters of each microseismic event point spatial cluster, and combined with the effective neighboring event point pairs that cross different microseismic event point spatial clusters in the spatiotemporally consistent proximity relationship set, perform spatial connection relationship inference processing on different microseismic event point spatial clusters, determine the sequential connection order between microseismic event point spatial clusters that have a temporal order and a continuous spatial orientation direction, and form a crack segment chain structure composed of multiple microseismic event point spatial clusters connected in sequence according to temporal order and spatial extension direction.

[0039] This step requires connecting the discrete clusters into a chain. First, review the set of spatiotemporal consistent proximity relationships, Set_temporal_consistent, obtained in step S125. This set contains not only event point pairs within the same cluster but also event point pairs belonging to different clusters. If the two event points i and j in a point pair (i, j) belong to clusters Cluster_A and Cluster_B respectively, and this point pair exists in Set_temporal_consistent, it means that clusters Cluster_A and Cluster_B are spatiotemporally adjacent. Then, combine the spatial coordinates of the cluster centers (Center_A, Center_B), the spatial direction (v1_A, v1_B), and their time span parameters (T_earliest_A, T_latest_A, T_earliest_B, T_latest_B) of clusters Cluster_A and Cluster_B. To determine the direction of the connection: Compare the latest oscillation time T_latest_A of cluster A with the earliest oscillation time T_earliest_B of cluster B. If T_latest_A is less than T_earliest_B, the possible flow direction is from cluster A to cluster B. Simultaneously, examine the vector pointing from the cluster center of cluster A to the cluster center of cluster B, denoted as Vec_AB. Calculate the dot product of Vec_AB and the spatial direction v1_A of cluster A. If the dot product is positive, it indicates that the direction of cluster A roughly points towards cluster B. Similarly, calculate the dot product of Vec_AB and the spatial direction v1_B of cluster B. If the dot product is negative, it indicates that the direction of cluster B roughly extends outward from the direction of cluster A. By comprehensively analyzing all spatiotemporal proximity pairs across clusters and using the aforementioned directional judgment, the order of connections between multiple clusters can be inferred, thus linking a series of discrete clusters together to form a chain structure Chain_m representing a continuous crack channel. Each chain structure Chain_m is an ordered list of clusters.

[0040] Step S1284: The spatial orientation direction and spatial coordinates of the cluster center of each microseismic event point represented by the spatial cluster of the crack segment in the crack segment chain structure are used as the spatial location and extension direction identifier of the crack segment. The spatial intersection points between different crack segment chain structures are used as the branch nodes of the crack network. Based on the serial relationship of the crack segment chain structures and the branch relationship of different crack segment chain structures at the spatial intersection points, an initial pressure crack network spatial skeleton containing multiple crack segment chain structures and branch nodes is constructed.

[0041] After forming multiple chain structures Chain_m, it is necessary to identify their intersection relationships. For example, if the spatial coordinates of the cluster centers of Chain_1 and Chain_2 are very close (less than a preset intersection distance threshold D_junction), and they are associated through point pairs in the set Set_temporal_consistent, or their spatial directions intersect in a certain spatial region, then this intersection region can be defined as a branch node Junction_p. This branch node Junction_p has a specific coordinate position in three-dimensional space, which can be the average of the coordinates of the cluster centers of the intersecting clusters. Finally, all the chain structures Chain_m, their branch nodes Junction_p, and their connection relationships are integrated. Each chain structure consists of a series of ordered clusters, each with its own spatial location and orientation. Different chain structures intersect at branch nodes. This constructs a preliminary network skeleton, namely the initial spatial skeleton of the hydraulic fracture network, Skeleton_initial, which describes the orientation, extension, and branching relationships of the main fracture channels. The data structure of this initial spatial skeleton of the hydraulic fracture network can be a graph, where nodes represent cluster centers or branch nodes, and edges represent the connections between clusters.

[0042] Step S130: Based on the initial fracture network spatial skeleton, perform tracer migration path matching on the inter-well tracer monitoring dataset, and combine the tracer breakthrough time record and the tracer concentration change curve to determine the spatial location and connection direction of the dominant tracer migration channel in the main fracture channel, and generate a fracture network connectivity topology that includes the connectivity relationship of the dominant migration channel.

[0043] After constructing the fracture skeleton based on microseismic data, tracer data is needed to verify and refine which channels are the actual dominant flow paths of the fluid.

[0044] Step S131: Perform time series analysis on the tracer breakthrough time records in the inter-well tracer monitoring dataset, extract the first tracer breakthrough time between each injection well and each production well, and calculate the tracer migration time difference between each injection well and each production well based on the injection start time of each injection well and the first tracer breakthrough time of each production well, generating a tracer migration time matrix containing the tracer migration time differences between all injection wells and production wells.

[0045] From the inter-well tracer monitoring dataset, the tracer breakthrough time record matrix Matrix_breakthrough is extracted. The element T_breakthrough_ij in this matrix represents the time point at which the tracer injected from injection well i is first detected in production well j. Simultaneously, the tracer injection start time T_inject_i for each injection well i is known. For each pair of injection wells i and production well j, the tracer migration time difference Δt_migration_ij is calculated, using the formula: Δt_migration_ij equals T_breakthrough_ij minus T_inject_i. All calculated Δt_migration_ij are arranged into a new matrix M_migration_time, with injection well i as the row and production well j as the column. This matrix is ​​the tracer migration time-effect matrix. This tracer migration time-effect matrix quantifies the fluid migration rate between different well pairs.

[0046] Step S132: Perform concentration peak feature extraction processing on the tracer concentration change curves in the well inter-well tracer monitoring dataset, identify the maximum concentration point and the corresponding occurrence time point in each tracer concentration change curve, and extract the maximum concentration value and the corresponding occurrence time point of each tracer concentration change curve as the peak concentration feature parameter of the tracer concentration change curve. The peak concentration feature parameter includes the peak concentration value and the peak arrival time point value.

[0047] For the tracer concentration variation curve of each production well j, which is a time series consisting of a series of (T_sample, C_sample) data pairs, we iterate through all data points, find the maximum value of the concentration C_sample, denoted as C_peak_j, and record the sampling time point corresponding to this maximum value, denoted as T_peak_j. C_peak_j and T_peak_j together constitute the peak concentration characteristic parameter of production well j. T_peak_j differs from T_breakthrough_ij; it represents the moment when the tracer concentration reaches its highest peak and contains richer transport information, such as the degree of channel dispersion.

[0048] Step S133: Based on the spatial orientation of the main fracture channels and the spatial coordinates of the branch nodes in the initial fracture network spatial skeleton, generate a three-dimensional fracture channel grid covering the spatial extension range of all main fracture channels in the geological coordinate system of the target fracture well area. The three-dimensional fracture channel grid is composed of multiple continuously arranged grid cells. Each grid cell corresponds to a spatial volume element and has the spatial coordinates of the grid cell center in the geological coordinate system.

[0049] Based on the distribution range of all major fracture channels in the initial fracture network spatial skeleton Skeleton_initial, a three-dimensional spatial region capable of covering all channels is determined. Within this region, a uniform three-dimensional mesh, denoted as Grid_fracture, is created. The mesh consists of numerous cubic mesh cells of the same size, with the side length of each mesh cell denoted as L_cell. Each mesh cell represents a spatial volume element, with a volume equal to the cube of L_cell. Each mesh cell has a unique index (i, j, k), and the spatial coordinates of its center point (X_center_ijk, Y_center_ijk, Z_center_ijk) can be calculated based on its index and the coordinates of the mesh origin.

[0050] Step S134: Map the tracer migration time difference between each injection well and each production well in the tracer migration time matrix to the fracture channel three-dimensional grid. Based on the injection wellhead spatial coordinates of each injection well in the geological coordinate system and the production wellhead spatial coordinates of each production well in the geological coordinate system, mark all possible paths from the injection wellhead spatial coordinates to the production wellhead spatial coordinates through consecutive adjacent grid cells in the fracture channel three-dimensional grid, and generate multiple candidate tracer migration path sets.

[0051] For each pair of injection wells i and production wells j, their wellhead spatial coordinates are known, denoted as P_inj_i and P_prod_j, respectively. On the 3D mesh Grid_fracture, a path search algorithm, such as the A* algorithm or Dijkstra's algorithm, is used to find all possible paths from the mesh cell containing P_inj_i to the mesh cell containing P_prod_j. The search space for paths can be limited to the mesh range defined in step S133. The path construction rule is: starting from the initial mesh cell, each step can move to adjacent mesh cells sharing a face or edge (i.e., six-neighborhood or twenty-six-neighborhood) until the target mesh cell is reached. Each path consists of a sequence of coordinates of the center points of consecutive adjacent mesh cells, denoted as Path_ij_k, where k is the path index. Thus, for each pair of wells, a candidate tracer migration path set Set_paths_ij can be obtained, containing multiple possible paths connecting injection well i and production well j.

[0052] Step S135: Perform path length calculation processing on each candidate tracer migration path set, count the total number of consecutive adjacent grid cells contained in each candidate tracer migration path set, and calculate the geometric path length value of each candidate tracer migration path based on the side length of the spatial volume element corresponding to each grid cell in the three-dimensional grid of the crack channel, and generate a candidate path length set containing the geometric path length values ​​of all candidate tracer migration paths.

[0053] For each candidate path Path_ij_k, it consists of a series of grid cell center point coordinates: P0, P1, P2, ..., Pm. The geometric path length L_path_ij_k of this candidate path is equal to the sum of the Euclidean distances between all adjacent points. Since the grid is uniform, the distance between the center points of adjacent grid cells is either L_cell (adjacent to shared faces), L_cell multiplied by √2 (adjacent to shared edges), or L_cell multiplied by √3 (adjacent to shared corners). Therefore, the total length can be easily calculated by counting the number of different types of movement steps in the path. Alternatively, the Euclidean distances of all adjacent point pairs can be directly calculated and summed. All calculated L_path_ij_k values ​​are associated with their corresponding path identifiers to generate a set of candidate path lengths, Set_path_lengths.

[0054] Step S136: Based on the tracer migration time difference corresponding to the tracer migration time matrix and the geometric path length of each candidate tracer migration path in the candidate path length set, calculate the average tracer migration velocity corresponding to each candidate tracer migration path. Compare the average tracer migration velocity with a pre-configured standard migration velocity range threshold for tracers in the crack channel, and select candidate tracer migration paths whose average tracer migration velocity falls within the standard migration velocity range threshold as valid tracer migration paths with velocity consistency, thus obtaining a set of valid tracer migration paths.

[0055] For each pair of injection wells i and production wells j, and their corresponding candidate path Path_ij_k, the geometric length L_path_ij_k of the path and the actual tracer migration time difference Δt_migration_ij for the well pair obtained from the tracer migration time matrix M_migration_time are known. The theoretical average tracer migration velocity v_ij_k corresponding to this path is calculated using the formula: v_ij_k equals L_path_ij_k divided by Δt_migration_ij. A standard migration velocity range threshold for the tracer in the fracture channel is preset, with a lower limit of V_min and an upper limit of V_max. This standard migration velocity range is determined based on parameters such as the rock physical properties and fracturing fluid viscosity of the region. The calculated v_ij_k is compared with [V_min, V_max]. If v_ij_k falls within this range, the path is considered reasonable from a velocity perspective and is marked as a valid tracer migration path. All paths that meet the criteria constitute the set of effective tracer transport paths, Set_effective_paths.

[0056] Step S137: Based on the set of effective tracer transport paths, perform spatial identification of advantageous transport channels and construction of connectivity relationships to generate a fracture network connectivity topology.

[0057] After obtaining the set of effective tracer transport paths (Set_effective_paths), it is necessary to identify which paths are the most advantageous flow channels for the fluid and construct their connectivity.

[0058] For example, step S1371: Perform grid tracer concentration allocation processing on each effective tracer migration path in the effective tracer migration path set. Based on the peak concentration value and peak arrival time value in the peak concentration characteristic parameters, and combined with the spatial distance ratio between the spatial coordinates of the center of the grid cell traversed by each effective tracer migration path and the spatial coordinates of the injection wellhead and the production wellhead, assign a tracer concentration contribution weight value to each grid cell traversed by each effective tracer migration path, and generate a path concentration distribution matrix containing the tracer concentration contribution weight value of each grid cell.

[0059] For each path Path_ij_k in the effective tracer migration path set Set_effective_paths, connecting injection well i and production well j, the peak concentration characteristic parameters (C_peak_j, T_peak_j) of production well j are obtained. For each grid cell traversed by path_ij_k, with its center point coordinates P_g, its weight in the entire path needs to be calculated. A distance-based linear decay model can be used. First, calculate the distance d_inj from P_g to injection well i and the distance d_prod to production well j. The total path length is L_path_ij_k. Then, the tracer concentration contribution weight w_g of this grid cell can be set as a function related to d_inj and d_prod, such as w_g equal to C_peak_j multiplied by (d_inj divided by L_path_ij_k) or a more complex model. Alternatively, the time factor can be considered; the tracer concentration contribution of grid cells closer to the injection well may also be related to the time T_peak_j when the peak occurs. Using the method described above, a weight value w_g is calculated for each grid cell on each path. The weight values ​​of all grid cells across all paths are then aggregated to form a mapping, Map_grid_weights, with the grid cell index as the key and the list of weight values ​​as the value.

[0060] Step S1372: The tracer concentration contribution weight values ​​of all effective tracer transport paths in the path concentration distribution matrix are superimposed and summed in the same grid cell. The comprehensive tracer concentration characterization value of each grid cell is calculated. The grid cells whose comprehensive tracer concentration characterization value is greater than the pre-configured concentration contribution significance threshold are marked as target contribution grid cells occupied by the dominant tracer transport channel, and a set of target contribution grid cells is generated.

[0061] For each grid cell G in Map_grid_weights, it may appear on multiple effective paths, thus having a list of weight values ​​[w_g1, w_g2, ...]. The comprehensive tracer concentration characterization value S_G for this grid cell is calculated by summing all the weight values ​​in the list; S_G equals w_g1 plus w_g2 plus .... A preset concentration contribution significance threshold S_threshold is used. All grid cells are iterated through, and those with S_G greater than S_threshold are marked as target contribution grid cells occupied by the dominant tracer transport channels and added to the set Set_dominant_cells. These grid cells represent the regions where multiple effective paths share the most concentrated fluid flow.

[0062] Step S1373: Perform spatial connectivity analysis on the target contribution grid cells in the target contribution grid cell set, identify all groups of target contribution grid cells in the target contribution grid cell set that are adjacent to each other through shared surfaces or shared edges, and take each group of adjacent target contribution grid cells as a tracer dominant transport channel segment to obtain multiple tracer dominant transport channel segments.

[0063] For the mesh cells in the set `Set_dominant_cells`, perform a 3D spatial connectivity analysis. Starting with any mesh cell, check if other mesh cells adjacent to it via shared faces (six-neighborhood) or shared edges (eighteen-neighborhood) are also in `Set_dominant_cells`. If so, group these adjacent cells into the same group and continue expanding the search from newly added cells until no new adjacent cells are found. This identifies a connected mesh cell group, denoted as `Segment_q`. Repeat this process until all mesh cells in `Set_dominant_cells` are assigned to a `Segment_q`. Each `Segment_q` represents a dominant tracer transport channel segment, which may be discontinuous, separated by non-dominant mesh cells.

[0064] Step S1374: Based on the spatial distribution of the multiple tracer dominant transport channel segments in the three-dimensional grid of the fracture channel, and combined with the spatial orientation of the main fracture channel in the initial fracture network spatial skeleton, determine the corresponding spatial interval of each tracer dominant transport channel segment on the main fracture channel, and connect all tracer dominant transport channel segments in the corresponding spatial interval according to their spatial adjacency to construct a continuous trajectory of dominant transport channels composed of continuous tracer dominant transport channel segments connected in series.

[0065] Project the dominant transport channel segments Segment_q obtained in step S1373 onto the initial fracture network spatial skeleton Skeleton_initial constructed in step S120. For each Segment_q, calculate the average of the center point coordinates of all its contained grid cells to obtain the center position of the segment. Find the fracture segment or branch node closest to this center position on Skeleton_initial to determine the corresponding spatial interval of Segment_q on the skeleton. Then, along the topology of Skeleton_initial, connect those Segments_q that are spatially adjacent (i.e., their projection intervals are continuous or close on the skeleton) end-to-end according to their order on the skeleton. During connection, it may be necessary to fill some grid cells on the skeleton path to enable connectivity between segments. Finally, a series of continuous dominant transport channel trajectories Trajectory_r, consisting of a sequence of grid cell coordinates, are formed. Each Trajectory_r corresponds to a main channel with dominant fluid flow.

[0066] Step S1375: Extract the coordinates of the starting and ending grid center points of each dominant migration channel continuous trajectory in the three-dimensional grid of the crack channel. Determine the migration start point and migration end point of the dominant migration channel continuous trajectory based on the coordinates of the starting and ending grid center points. Calculate the average direction vector of the line connecting the center points of all adjacent grid cells in the dominant migration channel continuous trajectory. Use the spatial direction of the average direction vector as the dominant migration direction of the dominant migration channel continuous trajectory, generating a set of dominant migration channels containing multiple dominant migration channel continuous trajectories and their corresponding dominant migration directions.

[0067] For each dominant migration channel's continuous trajectory *Trajectory_r*, it consists of a series of ordered grid cell center point coordinates: *P_start*, ..., *P_end*. *P_start* is the coordinate of the starting grid center point of the trajectory, typically near an injection well; *P_end* is the coordinate of the ending grid center point, typically near a production well. *P_start* and *P_end* are used as the migration start and end points of the channel, respectively. Then, the dominant migration direction of the trajectory is calculated. For all adjacent point pairs (*P_i*, *P_i+1*) on the trajectory, the direction vector *v_i* from *P_i* to *P_i+1* is calculated. All *v_i* vectors are vector-superimposed to obtain the total vector *V_total*. Then, *V_total* is normalized to obtain a unit vector *v_dominant*, whose spatial direction is the dominant migration direction of the continuous trajectory *Trajectory_r*. All Trajectory_r and their corresponding dominant transport directions v_dominant_r, along with the grid cell coordinate sequences they contain, are collectively formed into the dominant transport channel set Set_dominant_trajectories.

[0068] Step S1376: The coordinate sequence of the center point of the grid cell occupied by the continuous trajectory of each dominant migration channel in the set of dominant migration channels is used as the spatial location identifier of the dominant migration channel of the tracer. The connection relationship between different dominant migration channel segments connected by the continuous trajectory of the dominant migration channel is used as the channel connectivity identifier. Based on the channel connectivity identifier, a dominant migration channel connectivity graph is constructed with each continuous trajectory of the dominant migration channel as a node and the channel connectivity relationship as an edge, generating a crack network connectivity topology that includes the connectivity relationship of the dominant migration channels.

[0069] Each dominant transport path's continuous trajectory, `Trajectory_r`, is represented as a node `Node_r` in the graph. The node's attributes include its spatial location identifier (i.e., the coordinate sequence of its grid cell center points) and its dominant transport direction, `v_dominant_r`. If two trajectories, `Trajectory_r` and `Trajectory_s`, are spatially connected (e.g., the endpoint of one trajectory is adjacent to the starting point of another, or they share some grid cells), an edge `Edge_rs` is created between their corresponding nodes, `Node_r` and `Node_s`. The edge's attributes can include information such as the connection type (direct connection or indirect connection through branch nodes). The resulting directed or undirected graph, `Graph_connectivity`, is the fracture network connectivity topology. This fracture network connectivity topology reveals which tracer dominant transport paths are connected and the connections between them.

[0070] Step S140: Based on the fracture network connectivity topology, the initial pressure fracture network spatial skeleton is supplemented with fracture network integrity, and the spatial location corresponding to the dominant transport channel is spatially associated and fused with the main fracture channel in the initial pressure fracture network spatial skeleton to generate a complete three-dimensional shape of the pressure fracture network including the main fracture channel and the secondary connected fracture channel.

[0071] After obtaining the initial skeleton based on microseismic data and the connectivity topology based on tracer data, the two need to be merged to obtain a more complete fracture network.

[0072] Step S141: Analyze the dominant migration channel connectivity graph in the fracture network connectivity topology, extract the spatial coordinate sequence of the continuous trajectory of the dominant migration channel corresponding to each node in the dominant migration channel connectivity graph, and the channel connectivity relationship type parameter corresponding to the edge between each node. The channel connectivity relationship type parameter includes direct connectivity type and indirect connectivity type.

[0073] The crack network connectivity topology Graph_connectivity constructed in step S137 is parsed. All nodes Node_r in the graph are traversed, and the spatial coordinate sequence stored in their attributes is extracted, i.e., the coordinate sequence of the grid cell center points List_coords_r of Trajectory_r. All edges Edge_rs in the graph are traversed, and the channel connectivity type parameter Type_rs stored in their attributes is extracted. For example, Type_rs equal to "direct" indicates that two dominant channels are directly connected in space, and Type_rs equal to "indirect" indicates that they need to be indirectly connected through other channels.

[0074] Step S142: Perform spatial coordinate matching processing on the spatial coordinate sequence of the chain structure and branch nodes of each major fracture channel in the initial fracture network spatial skeleton and the spatial coordinate sequence of each node in the dominant migration channel connectivity diagram. Calculate the closest point pair between the spatial coordinate sequence of each major fracture channel and the spatial coordinate sequence of the continuous trajectory of each dominant migration channel to obtain the set of spatial proximity matching relationships between the major fracture channels and the tracer dominant migration channels.

[0075] The initial fracture network spatial skeleton, Skeleton_initial, consists of multiple fracture segment chain structures, Chain_m, and branch nodes, Junction_p. Each Chain_m consists of a series of cluster center coordinates. Similarly, the nodes Node_r in the dominant migration channel connectivity graph correspond to a series of grid cell center point coordinates, List_coords_r. For each Chain_m and each Node_r, the spatially closest point pair between them is calculated. This can be done by traversing all points C_i on Chain_m and all points G_j on List_coords_r, calculating the Euclidean distance of all point pairs (C_i, G_j), and finding the point pair with the smallest distance (C_min, G_min) and its distance value D_min. All point pairs that satisfy D_min being less than a certain large search radius are recorded, forming a spatial proximity matching set, Set_proximal_matches.

[0076] Step S143: Perform spatial overlap evaluation processing on each spatial proximity matching relationship in the spatial proximity matching relationship set. Based on the coordinate difference of the closest point pair between the main fracture channel and the tracer dominant transport channel and the side length of each grid cell in the three-dimensional grid of the fracture channel, calculate the overlap ratio parameter of the main fracture channel and the tracer dominant transport channel in the spatial dimension. Mark the spatial proximity matching relationships with the overlap ratio parameter greater than the pre-configured overlap threshold as spatial coincident matching relationships, and generate a set of spatial coincident matching relationships.

[0077] For each matching pair (Chain_m, Node_r) in Set_proximal_matches, its nearest point pair (C_min, G_min) and distance D_min are known. To assess whether the two are not merely neighbors but spatially overlapping, the overlap ratio needs to be calculated. A neighborhood can be considered centered on the grid cell containing G_min, for example, a cube region with side length k times L_cell centered on G_min. The number of points on Chain_m falling within this cube region is counted, denoted as Count_overlap. Simultaneously, the number of points on Node_r falling within this cube region is counted, denoted as Count_node_overlap. The overlap ratio parameter O_ratio can be defined as Count_overlap divided by Count_node_overlap, or a more complex comprehensive metric. A preset overlap threshold O_threshold is used. If O_ratio is greater than O_threshold, then Chain_m and Node_r are considered to be spatially coincident. The matching pair is marked as a spatially coincident matching relationship and added to the spatially coincident matching relationship set Set_coincident_matches.

[0078] Step S144: Based on each spatial overlap matching relationship in the set of spatial overlap matching relationships, the main fracture channel and the tracer dominant transport channel with spatial overlap matching relationship are merged. Using the spatial coordinate sequence of the main fracture channel as a reference, the coordinates of the center points of the grid cells in the spatial coordinate sequence of the continuous trajectory of the dominant transport channel that do not overlap with the spatial coordinate sequence of the main fracture channel are inserted into the corresponding positions of the spatial coordinate sequence of the main fracture channel to generate the main fracture fusion enhancement coordinate sequence.

[0079] For each pair (Chain_m, Node_r) in Set_coincident_matches, channel merging is performed. The spatial coordinate sequence List_chain (i.e., a series of cluster center coordinates) of Chain_m is used as a reference. The spatial coordinate sequence List_node of Node_r is a series of more densely packed grid cell center coordinates. Points in List_node that are not near List_chain (i.e., distance greater than a certain threshold) are identified; these points represent fracture segments that were not detected by microseismic data but are indicated by tracer data. These points are then inserted, according to their spatial location, between the two nearest neighbors in List_chain, forming a new, denser coordinate sequence List_enhanced_m. This coordinate sequence is the enhanced spatial coordinate sequence of the merged main fracture channel.

[0080] Step S145: Identify and process the tracer dominant transport channels in the dominant transport channel connectivity graph that have connectivity with the spatial overlap matching relationship set but do not have a spatial overlap matching relationship with any major fracture channel, and extract the tracer dominant transport channels as candidates for secondary fracture channels to be supplemented, thereby obtaining a secondary fracture channel candidate set.

[0081] Iterate through all nodes (Node_s) in the dominant migration channel connectivity graph (Graph_connectivity). If a node (Node_s) does not appear in any spatially overlapping match (Set_coincident_matches) (i.e., it does not coincide with any primary fracture channel), but it is connected to a node (Node_r) that appears in Set_coincident_matches via an edge (Edge_rs), then the dominant migration channel represented by node_s is likely a secondary fracture connecting to the primary channel. Identify these nodes (Node_s) as candidates for secondary fracture channels and add them to the secondary fracture channel candidate set (Set_candidate_secondary).

[0082] Step S146: Perform spatial extension direction analysis on each secondary crack channel candidate to be supplemented in the secondary crack channel candidate set, calculate the direction vector of the line connecting the center points of all adjacent grid cells in the spatial coordinate sequence of each secondary crack channel candidate, count the main direction distribution frequency of the direction vector, and take the direction vector with the highest frequency as the main extension direction of the secondary crack channel candidate.

[0083] For each candidate channel Node_c in Set_candidate_secondary, its spatial coordinate sequence is List_coords_c. Traverse all adjacent point pairs (P_i, P_i+1) in this sequence and calculate the direction vector v_i from P_i to P_i+1. Perform cluster analysis or statistics on all these v_i vectors, for example, projecting all vectors onto a unit sphere and counting which region has the most vectors. The direction with the highest frequency is the main extension direction v_main_c of this candidate channel. This main extension direction will be used to subsequently determine its compatibility with the main crack.

[0084] Step S147: Based on the main fracture fusion enhancement coordinate sequence, perform secondary fracture direction consistency verification and spatial fusion processing to generate a complete three-dimensional morphology of the pressure fracture network.

[0085] After obtaining the primary fracture fusion enhancement coordinate sequence and the secondary fracture channel candidate set, the final fusion is performed.

[0086] For example, in step S1471: the main extension direction of the secondary fracture channel candidate is compared with the spatial orientation direction of the adjacent main fracture channel in the initial pressure fracture network spatial skeleton. The cosine value of the direction angle between the main extension direction and the spatial orientation direction of the adjacent main fracture channel is calculated. The secondary fracture channel candidate with the direction angle cosine value greater than the pre-configured direction consistency threshold is marked as a valid secondary fracture channel with extension direction consistency, and a set of valid secondary fracture channels is generated.

[0087] For each candidate secondary fracture channel Node_c, based on its spatial location, find the nearest main fracture channel Chain_nearby in the initial fracture network spatial skeleton Skeleton_initial, and obtain the spatial orientation direction v_main_chain of this main fracture channel (obtainable from step S1282). Calculate the cosine value cosθ of the angle between v_main_c and v_main_chain, i.e., their dot product. A pre-set orientation consistency threshold cos_threshold is used, for example, 0.7 or 0.8. If cosθ is greater than cos_threshold, it indicates that the extension direction of the secondary fracture is basically consistent with the orientation of the adjacent main fracture, conforming to geomechanical laws. It is then marked as an effective secondary fracture channel and added to the effective secondary fracture channel set Set_effective_secondary.

[0088] Step S1472: Based on the spatial positional relationship between the spatial coordinate sequence of each effective secondary fracture channel in the effective secondary fracture channel set and the coordinate sequence of the main fracture fusion enhancement, determine the coordinates of the connection point between each effective secondary fracture channel and the main fracture channel in space. The coordinates of the connection point are the coordinates of the center point of the grid cell that is closest to the coordinate sequence of the main fracture fusion enhancement in the spatial coordinate sequence of the effective secondary fracture channel.

[0089] For each effective secondary fracture channel Node_eff in Set_effective_secondary, its spatial coordinate sequence is List_coords_eff. The primary fracture fusion enhancement coordinate sequence is List_enhanced (there may be multiple lists, corresponding to multiple primary fractures). Iterate through each point P_eff in List_coords_eff, calculate its Euclidean distance to all points in List_enhanced, find the point P_connect_on_main corresponding to the global minimum distance, and use P_eff itself as P_connect_on_secondary. Record this connection point pair (P_connect_on_secondary, P_connect_on_main). P_connect_on_secondary is the coordinate of the connection point between this secondary fracture channel and the primary fracture.

[0090] Step S1473: Using the coordinates of the connection point as the branch starting point, take the coordinates of all grid cell center points from the branch starting point to the end of the effective secondary crack channel in the spatial coordinate sequence of each effective secondary crack channel as the branch trajectory coordinate sequence of the effective secondary crack channel, and generate a set of secondary crack branches containing the branch trajectory coordinate sequences of all effective secondary crack channels and the corresponding branch starting point coordinates.

[0091] Starting from the connection point P_connect_on_secondary determined in step S1472, proceed along the spatial coordinate sequence List_coords_eff of the effective secondary fracture channels Node_eff, moving away from the main fracture (i.e., towards the end of the channel). Extract all coordinate points from P_connect_on_secondary to the last point of the sequence to form the branch trajectory coordinate sequence Branch_coords_eff of the channel. Simultaneously record the coordinates of the branch starting point P_connect_on_main (the corresponding point on the main fracture). Combine the Branch_coords_eff of all effective secondary fracture channels with their corresponding main fracture connection points P_connect_on_main to form the secondary fracture branch set Set_branches.

[0092] Step S1474: Spatially associate and fuse each branch trajectory coordinate sequence in the secondary fracture branch set with the main fracture fusion enhancement coordinate sequence according to its branch starting point coordinates. Mark the position corresponding to the starting point coordinate of each branch in the main fracture fusion enhancement coordinate sequence as a branch node, and attach the corresponding branch trajectory coordinate sequence as a new branch path to the fracture network backbone structure formed by the main fracture fusion enhancement coordinate sequence at the branch node, thereby forming an initial fused fracture network structure containing the main fracture channel enhancement coordinate sequence and the attached secondary fracture branch trajectory coordinate sequence.

[0093] For each branch in `Set_branches`, based on its branch starting point coordinates `P_connect_on_main`, find the corresponding position in the main fracture fusion enhancement coordinate sequence `List_enhanced` (this could be a point in the sequence or a point between two points), and mark it as a branch node `Node_branch`. Then, associate the branch trajectory coordinate sequence `Branch_coords_eff` as a new branch path with that branch node `Node_branch` on the main fracture. The final data structure is a tree or network structure containing the enhanced main fracture coordinate sequence and the secondary fracture branch coordinate sequences branching off from specific nodes on these main fractures. This structure is denoted as `Initial_fused_network`.

[0094] Step S1475: Perform spatial gridding resampling processing on all the main crack channel enhancement coordinate sequences and secondary crack branch trajectory coordinate sequences in the initial fused crack network structure. Discretize all coordinate sequences into grid cell center point coordinates aligned with the three-dimensional grid of the crack channel according to the grid cell size of the crack channel, eliminate the uneven distribution of spatial coordinate points caused by the difference in the original sampling accuracy between different coordinate sequences, and generate a unified grid coordinate sequence set with uniform spatial coordinate point distribution.

[0095] Since the primary crack coordinate sequence `List_enhanced` may originate from the cluster centers (relatively sparse), while the secondary crack coordinate sequence `Branch_coords_eff` originates from the grid cell centers (relatively dense), their spatial resolutions are inconsistent. To address this issue, all coordinate points in the entire `Initial_fused_network` are remapped onto the crack channel 3D grid `Grid_fracture` created in step S133. For each coordinate point, the center point coordinates of its respective grid cell are used as a replacement. In this way, all coordinate points are unified under the same grid coordinate system, forming a new set of sequences composed of grid cell center point coordinates. For example, a primary crack may be represented by a series of consecutive grid cell center point coordinates, and a secondary crack may also be represented by a series of consecutive grid cell center point coordinates. This new set is denoted as `Unified_grid_coords_set`.

[0096] Step S1476: Based on the position of the grid cell corresponding to each coordinate point in the unified grid coordinate sequence set in the three-dimensional grid of the fracture channel, and combined with the connectivity relationship between each tracer dominant migration channel in the dominant migration channel connectivity diagram, a spatial connection relationship is established between different branch trajectory coordinate sequences in the unified grid coordinate sequence set. Spatially adjacent grid cells belonging to different secondary fracture branch trajectory coordinate sequences are connected through shared surfaces or shared edges to form a complete fracture network connectivity path network, generating a complete three-dimensional morphology of the pressure fracture network containing the main fracture channel and the secondary connected fracture channels.

[0097] The `Unified_grid_coords_set` already contains the coordinates of all primary and secondary fractures' grid cells. Now, based on the spatial location of these grid cells within the 3D grid `Grid_fracture` and the connectivity relationships in the fracture network connectivity topology `Graph_connectivity` generated in step S137, the connections between all fracture segments need to be finalized. For example, if the terminal grid cells of two secondary fracture branches are spatially adjacent (sharing a face or edge), and according to `Graph_connectivity` they should be connected, then a connection is established between these two grid cells, indicating that they are connected. Finally, all fracture segments represented by grid cells are integrated into a fully connected path network composed of grid cells, based on spatial adjacency and the connectivity revealed by the tracer data. This path network, containing all primary and secondary fractures and their connections, constitutes the complete 3D morphology of the fracture network, `Final_network_3d`. Its data form can be a set of connected grid cell indices or a graph with grid cells as nodes and adjacency relationships as edges.

[0098] Step S150: Extract fracture network parameters from the complete three-dimensional morphology of the fracture network to obtain a set of fracture network morphology data for the target fracture well area. The set of fracture network morphology data includes the total surface area parameter of the fracture network, the total volume parameter of the fracture network, the branch density parameter of the fracture network, and the tortuosity parameter of the fracture network.

[0099] After obtaining the complete three-dimensional morphology of the fracture network, quantitative parameters need to be extracted from it to evaluate the fracturing effect.

[0100] Step S151: Perform coordinate point traversal processing on all main fracture enhancement coordinate sequences and all secondary fracture branch coordinate sequences in the complete three-dimensional morphology of the fracture network, and count the total number of spatial coordinates of the center points of all grid cells in the complete three-dimensional morphology of the fracture network as the parameter of the total number of grid cells occupied by the fracture network.

[0101] Traverse all mesh cells contained in the complete three-dimensional fracture network morphology Final_network_3d. Count the total number of these mesh cells, denoted as N_total_cells.

[0102] Step S152: Based on the actual physical size of the spatial volume element corresponding to each grid cell in the three-dimensional grid of the fracture channel, multiply the total number of grid cells occupied by the fracture network by the volume value of the spatial volume element corresponding to each grid cell to calculate the total volume parameter of the fracture network occupied by the complete three-dimensional morphology of the fracture network.

[0103] Given that the side length of each mesh cell is L_cell, the volume of each mesh cell, V_cell, is equal to L_cell multiplied by L_cell multiplied by L_cell. The total volume parameter of the fracture network, V_total_network, is equal to N_total_cells multiplied by V_cell. This parameter characterizes the total volume occupied by the fracture system formed by fracturing in three-dimensional space.

[0104] Step S153: Perform surface mesh element identification processing on all main fracture enhancement coordinate sequences and all secondary fracture branch coordinate sequences in the three-dimensional morphology of the complete fracture network. Mark the mesh elements located at the edge of the fracture network in each spatial coordinate sequence and having at least one adjacent face that is not shared with any other fracture network mesh element as fracture network surface mesh elements. Count the total number of fracture network surface mesh elements as the number of fracture network surface elements.

[0105] For each mesh cell in Final_network_3d, check if its six adjacent mesh cells (east, west, south, north, top, bottom) are also in Final_network_3d. If at least one of the adjacent mesh cells corresponding to a mesh cell's faces is not in the crack network, then the mesh cell is a surface mesh cell. Count the number of all surface mesh cells that meet the condition, denoted as N_surface_cells.

[0106] Step S154: Based on the actual physical dimensions of the spatial volume element corresponding to each grid cell in the three-dimensional mesh of the fracture channel, calculate the total area of ​​the exposed surface of each grid cell marked as a fracture network surface grid cell; sum up the total exposed surface areas of all fracture network surface grid cells to calculate the total surface area parameter of the fracture network exposed by the complete three-dimensional morphology of the fracture network.

[0107] For each mesh cell labeled as a surface mesh cell, calculate the total area of ​​its exposed surfaces. A mesh cell has six faces, each with an area equal to L_cell multiplied by L_cell. For this mesh cell, count how many faces have no adjacent fracture mesh cells, denoted as N_exposed_faces. Then, the exposed surface area A_cell of this cell is equal to N_exposed_faces multiplied by (L_cell multiplied by L_cell). Summing A_cell for all N_surface_cells surface mesh cells yields the total surface area parameter A_total_surface of the fracture network. This total surface area parameter reflects the contact area between the fracture and the rock matrix.

[0108] Step S155: Identify and statistically process all branch nodes in the complete three-dimensional morphology of the fracture network, extract the connection points of the coordinate sequences of the main fracture channels and secondary fracture branches in the complete three-dimensional morphology of the fracture network as fracture network branch nodes, and count the total number of fracture network branch nodes as the number of fracture network branch nodes.

[0109] In the complete three-dimensional network morphology Final_network_3d, mesh elements that simultaneously belong to both primary and secondary fractures, or connect multiple secondary fractures, are identified as branch nodes. Specifically, if a mesh element is in the fracture network and has three or more adjacent mesh elements (i.e., it extends from fractures in at least three directions), then that mesh element can be considered a branch node. The total number of all such branch nodes is counted and denoted as N_branch_nodes.

[0110] Step S156: Calculate the branch density and tortuosity parameters of the fracture network based on the complete three-dimensional morphology of the fracture network, and generate a set of fracture network morphology data.

[0111] For example, step S1561: Based on the connected path network formed by all main fracture enhancement coordinate sequences and all secondary fracture branch coordinate sequences in the three-dimensional morphology of the complete fracture network, identify independent path segments extending from any fracture network branch node to the end grid cell or another fracture network branch node, and count the total number of the independent path segments as the fracture network branch number parameter.

[0112] In the graph structure of Final_network_3d, branch nodes and terminal nodes (nodes connecting only one adjacent crack mesh cell) are used as key points. A path search algorithm from graph theory is employed to identify all simple paths (paths that do not repeat nodes) from one branch node to another, or from one branch node to one terminal node. The number of all these independent path segments is counted and denoted as N_branches. This branch count parameter represents the branch richness of the crack network.

[0113] Step S1562: Divide the number of branches of the crack network by the total volume of the crack network to calculate the branch density parameter of the crack network per unit volume.

[0114] The fracture network branch density parameter D_branch is equal to N_branches divided by V_total_network. This parameter characterizes the density of fracture branches within a unit volume of rock.

[0115] Step S1563: Perform path tortuosity calculation on each main fracture enhancement coordinate sequence in the three-dimensional morphology of the complete fracture network. For each main fracture enhancement coordinate sequence, calculate the straight-line Euclidean distance between the coordinates of the starting grid center point and the coordinates of the ending grid center point in the main fracture enhancement coordinate sequence as the straight-line distance parameter of the main fracture channel, and calculate the sum of the lengths of the lines connecting the center points of all adjacent grid units in the main fracture enhancement coordinate sequence as the actual path length parameter of the main fracture channel.

[0116] For each main crack channel (whose coordinate sequence is List_main_enhanced), its starting point is P_start_main and its ending point is P_end_main. Calculate the straight-line distance parameter L_straight_main, which is the Euclidean distance between P_start_main and P_end_main. Then, traverse all adjacent point pairs on List_main_enhanced, accumulate their distances, and obtain the actual path length parameter L_actual_main.

[0117] Step S1564: Divide the actual path length parameter of each main fracture channel by the straight-line distance parameter of the main fracture channel to obtain the channel tortuosity value of each main fracture channel.

[0118] For each main crack, its tortuosity τ_main is equal to L_actual_main divided by L_straight_main. τ_main is a value greater than or equal to 1; the closer it is to 1, the straighter the channel; the larger it is, the more tortuous the channel.

[0119] Step S1565: Calculate the arithmetic mean of the channel tortuosity values ​​of all main fracture channels to obtain the fracture network tortuosity parameters of the three-dimensional morphology of the complete pressure fracture network.

[0120] The fracture network tortuosity parameter Tau_network is equal to the sum of τ_main of all main fracture channels divided by the number of main fracture channels. This fracture network tortuosity parameter describes the overall degree of curvature of the fracture network.

[0121] Step S1566: Combine and encapsulate the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the fracture network to generate a set of pressure fracture network morphology data containing the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the fracture network.

[0122] Finally, the four key parameters V_total_network, A_total_surface, D_branch, and Tau_network, along with other potentially needed parameters such as N_branch_nodes, are combined into a structured dataset denoted as Fracture_network_parameters. This dataset is the fracture network morphology dataset, which quantitatively describes the geometric morphology and topological characteristics of the fracture network in the target fractured well area.

[0123] Step S160: Input the total volume parameters, total surface area parameters, branch density parameters, and tortuosity parameters of the fracture network in the fracture network morphology data set into the pre-training model of the fracturing effect. The pre-training model of the fracturing effect includes multiple sequentially connected convolutional layers, pooling layers, and fully connected layers.

[0124] In this embodiment, after completing step S150 and obtaining the fracture network morphology data set of the target fractured well area, the following operations are performed to comprehensively evaluate the fracturing operation effect of the well area. First, the four key parameters in the fracture network morphology data set Fracture_network_parameters generated in step S1566—namely, the total volume parameter of the fracture network (denoted as V_total_network), the total surface area parameter of the fracture network (denoted as A_total_surface), the branch density parameter of the fracture network (denoted as D_branch), and the tortuosity parameter of the fracture network (denoted as Tau_network)—are organized and formatted to meet the preset input format requirements. These four parameters are combined into a one-dimensional input feature vector, denoted as Input_vector, with a shape of [1, 4], i.e., a vector containing four elements, with the element arrangement fixed as [V_total_network, A_total_surface, D_branch, Tau_network]. This input feature vector Input_vector is then fed into a pre-trained fracturing effect pre-training model, denoted as Model_evaluation. The architecture of this model consists of multiple interconnected neural network layers, specifically including a first convolutional layer, a first pooling layer, a second convolutional layer, a second pooling layer, and a last fully connected layer.

[0125] Step S161: In the pre-trained model of the fracturing effect, the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the fracture network are extracted by one-dimensional convolutional feature extraction through the first convolutional layer to generate a first convolutional feature map containing local correlation features of the four input parameters.

[0126] The input feature vector, Input_vector, first enters the first convolutional layer, denoted as Conv1. Conv1 is a one-dimensional convolutional layer configured with multiple convolutional kernels, each with a size of [1, K1], where K1 is a small integer, such as 2 or 3, representing that each kernel covers K1 adjacent parameters in the input vector. Since the input vector length is 4, the convolutional kernels slide along the parameter dimension (i.e., the second dimension). Each convolutional kernel performs a convolution operation with Input_vector, that is, the weights of the convolutional kernel are multiplied by the corresponding input element, summed, and a bias term is added. Then, it is passed through a non-linear activation function, such as a linear rectified function, to obtain a value in the output feature map. Through the parallel computation of multiple such convolutional kernels (denoted as F1), the Conv1 layer generates an output feature map, denoted as Feature_map_1. Feature_map_1 has a shape of [1, L1, F1], where L1 is the length dimension obtained after convolution and padding. Due to the appropriate padding method, L1 can be the same as or slightly different from the input length. F1 is the number of convolution kernels, representing the number of extracted feature channels. Each feature channel in Feature_map_1 contains some combination features of the original input parameters in the local neighborhood, such as the local association between V_total_network and A_total_surface, or the local association between D_branch and Tau_network, etc.

[0127] Step S162: Max pooling downsampling is performed on the first convolutional feature map through the first pooling layer to reduce the feature dimension of the first convolutional feature map and retain the main feature responses, thereby generating the first pooling feature map.

[0128] The Feature_map_1 output from the first convolutional layer is then fed into the first pooling layer, denoted as Pool1. Pool1 employs max pooling with a pooling window size of [1, P1], where P1 can be 2, and the stride is also set to P1. The pooling window slides along the length dimension of Feature_map_1, retaining only the maximum value within each window's coverage area and discarding other values. This compresses the length dimension of Feature_map_1; for example, if L1 is 4, P1 is 2, and the stride is 2, the output length L2 becomes 2. Simultaneously, the number of feature channels F1 remains unchanged. Max pooling effectively reduces the feature map size while preserving the most salient feature responses in each local region, reducing computational complexity in subsequent layers and enhancing the model's translation invariance to some extent. The output of the Pool1 layer is denoted as Feature_map_pool1, with a shape of [1, L2, F1].

[0129] Step S163: Perform deep one-dimensional convolutional feature extraction on the first pooling feature map through the second convolutional layer to capture the high-order interaction relationship between the total volume parameter of the crack network, the total surface area parameter of the crack network, the branch density parameter of the crack network, and the tortuosity parameter of the crack network, and generate the second convolutional feature map.

[0130] The output of the first pooling layer, Feature_map_pool1, is used as the input to the second convolutional layer, denoted as Conv2. Conv2 is also a one-dimensional convolutional layer, and its configuration can have more convolutional kernels, for example, F2, where F2 is greater than F1. The kernel size can be set to [1, K2]. Conv2 performs convolution operations on the length dimension and feature channel dimension of Feature_map_pool1, which is actually performing cross-channel, deeper-level feature combination and abstraction. Each convolutional kernel performs a weighted combination of all feature channels of the input, thereby learning the higher-order interaction relationships between the original four parameters, such as more complex non-linear combination patterns between the four parameters V_total_network, A_total_surface, D_branch, and Tau_network. After the output of Conv2 is processed by the activation function, the second convolutional feature map is obtained, denoted as Feature_map_2, with a shape of [1, L3, F2], where L3 is the length dimension after the second convolution.

[0131] Step S164: Perform global average pooling on the second convolutional feature map through the second pooling layer, and average all feature values ​​of each feature channel in the second convolutional feature map to generate a one-dimensional global feature vector.

[0132] The output of the second convolutional layer, Feature_map_2, is then fed into the second pooling layer, denoted as Pool2. Pool2 employs global average pooling. Unlike the previous local max pooling, global average pooling operates along the entire length dimension of the feature map. Specifically, for each feature channel (i.e., the third dimension) of Feature_map_2, the arithmetic mean of all L3 feature values ​​in that channel is calculated. After global average pooling, each feature channel receives a single value, and the shape of Feature_map_2 changes from [1, L3, F2] to [1, 1, F2]. Then, the first two dimensions are typically compressed into a one-dimensional vector, denoted as Global_feature_vector, with a length of F2. This vector encapsulates the most essential and global feature information extracted from the input parameters by the entire convolutional neural network, preparing it for the final regression or classification task.

[0133] Step S165: Input the one-dimensional global feature vector into the fully connected layer of the fracturing effect pre-training model, and perform linear transformation and nonlinear activation processing on the one-dimensional global feature vector through the fully connected layer to calculate the comprehensive scoring parameters of the fracturing effect of the target fracturing well area.

[0134] The one-dimensional global feature vector, Global_feature_vector, obtained in step S164, is fed into the fully connected layer of the model, denoted as FullyConnected. This fully connected layer consists of multiple neurons, each fully connected to all elements in Global_feature_vector. Specifically, for a neuron in the fully connected layer, the calculation process is as follows: multiply each element in Global_feature_vector by its corresponding weight coefficient, then sum all the products, add a bias term, and finally pass it through a non-linear activation function, such as a linear rectified function or a sigmoid function, to obtain the neuron's output value. If the fully connected layer has only one output neuron, the output value of that neuron is a scalar, representing the comprehensive score parameter for fracturing effect, denoted as Score_final. If the fully connected layer has multiple neurons, another output layer can be connected to obtain the final scalar score. Score_final is a dimensionless numerical value, the magnitude of which reflects the quality of fracturing effect; a higher value indicates a better fracturing effect.

[0135] Step S166: Based on the comparison between the comprehensive scoring parameters of the fracturing effect and the preset fracturing effect level threshold, determine the fracturing effect level of the target fracturing well area, and generate a fracturing effect evaluation result that includes the fracturing effect level and the comprehensive scoring parameters of the fracturing effect.

[0136] A set of threshold values ​​for classifying fracturing effectiveness is preset, for example, three thresholds can be set: Threshold_excellent, Threshold_good, and Threshold_fair, classifying fracturing effectiveness into four levels: "Excellent," "Good," "Pass," and "Fail." The comprehensive fracturing effectiveness score parameter Score_final calculated in step S165 is compared with these threshold values ​​sequentially. If Score_final is greater than or equal to Threshold_excellent, the fracturing effectiveness level is determined to be "Excellent"; if Score_final is between Threshold_good and Threshold_excellent, the level is "Good"; if Score_final is between Threshold_fair and Threshold_good, the level is "Pass"; if Score_final is less than Threshold_fair, the level is "Fail." The determined fracturing effectiveness level is recorded as Level_evaluation, and combined with Score_final to form the fracturing effectiveness evaluation result Result_evaluation, which can be stored, for example, as a data structure containing two fields.

[0137] Step S167: Associate and store the fracturing effect evaluation results with the fracturing fracture network morphology data set, and generate a visualization report containing the fracturing effect evaluation results. The visualization report includes a three-dimensional rendering of the complete three-dimensional fracturing fracture network morphology and a labeling of the fracturing effect level.

[0138] Finally, the fracturing effect evaluation result Result_evaluation generated in step S166 is associated with the fracture network morphology data set Fracture_network_parameters obtained in step S150, and the complete three-dimensional fracture network morphology Final_network_3d generated in step S1476, and stored together in the database for subsequent querying and analysis. Simultaneously, a three-dimensional visualization engine is invoked to generate a three-dimensional rendering of the fracture network for the target fracturing well area based on the grid cell coordinate data of Final_network_3d. In the generated visualization report page, this three-dimensional rendering is displayed along with the fracturing effect level Level_evaluation, the comprehensive fracturing effect score parameter Score_final, and four key fracture network morphology parameters (V_total_network, A_total_surface, D_branch, Tau_network), forming an intuitive and comprehensive fracturing effect analysis and evaluation report.

[0139] Step S170: Perform connectivity quality assessment on the fracture network connectivity topology, extract the grid center point coordinate sequence of the continuous trajectory of each dominant transport channel in the fracture network connectivity topology, and calculate the sum of the Euclidean distances between the center points of all adjacent grid cells in the continuous trajectory of each dominant transport channel as the actual transport path length of the tracer dominant transport channel.

[0140] After generating the crack network connectivity topology in step S130, the following operations are performed to evaluate the quality of this connectivity. First, from the crack network connectivity topology Graph_connectivity generated in step S1376, each node Node_r is traversed, and each node corresponds to a dominant transport channel continuous trajectory Trajectory_r. The spatial coordinate sequence List_coords_r of Trajectory_r is extracted, which consists of a series of grid cell center point coordinates. The actual transport path length L_actual_tracer_r of this trajectory is calculated. Specifically, it is calculated by traversing all adjacent coordinate point pairs in List_coords_r, calculating the Euclidean distance between each pair of points, and summing all these distance values. The sum is L_actual_tracer_r.

[0141] Step S171: Obtain the straight-line Euclidean distance between the spatial coordinates of the injection wellhead and the spatial coordinates of the production wellhead corresponding to the continuous trajectory of each dominant migration channel as the straight-line connection distance of the dominant migration channel of the tracer.

[0142] For the same dominant migration channel continuous trajectory *Trajectory_r*, based on the injection well *i* and production well *j* associated with it during construction in step S137, the wellhead spatial coordinates *P_inj_i* of injection well *i* and *P_prod_j* of production well *j* are obtained. The straight-line Euclidean distance *L_straight_ij* between these two points is calculated, which is the spatial distance between *P_inj_i* and *P_prod_j*. This distance represents the ideal straight-line connection distance between the injection well and the production well, without considering fracture curvature.

[0143] Step S172: Divide the actual migration path length of each tracer dominant migration channel by the straight-line connectivity distance of that tracer dominant migration channel to calculate the channel tortuosity connectivity index of each tracer dominant migration channel.

[0144] For each dominant transport channel's continuous trajectory *Trajectory_r*, the tortuosity connectivity index *τ_tracer_r* is calculated using the actual transport path length *L_actual_tracer_r* calculated in step S170 and the straight-line connectivity distance *L_straight_ij* calculated in step S171. The formula is: *τ_tracer_r* = *L_actual_tracer_r* divided by *L_straight_ij*. This index reflects the curvature of the tracer's actual transport path relative to an ideal straight line. A larger *τ_tracer_r* indicates a more tortuous channel, potentially implying a more complex flow path or higher flow resistance.

[0145] Step S173: Based on the first breakthrough time of the tracer in the tracer breakthrough time record corresponding to each tracer dominant transport channel, calculate the average tracer transport velocity corresponding to the tracer dominant transport channel, divide the average tracer transport velocity by the pre-configured standard tracer transport velocity, and calculate the transport efficiency connectivity index of each tracer dominant transport channel.

[0146] For each dominant migration channel's continuous trajectory *Trajectory_r*, its associated injection well *i* and production well *j* are retrieved again. From the tracer migration time matrix *M_migration_time* generated in step S131, the tracer migration time difference Δt_migration_ij for that well pair is extracted. Using the actual migration path length *L_actual_tracer_r* obtained in step S170, the actual average tracer migration velocity *v_actual_r* of the channel is calculated, with the formula: *v_actual_r* equals *L_actual_tracer_r* divided by Δt_migration_ij. A theoretical standard tracer migration velocity *v_standard* is preset based on the rock physics and fluid properties of the region. The migration efficiency connectivity index η_r of the channel is calculated, with the formula: η_r equals *v_actual_r* divided by *v_standard*. η_r greater than 1 indicates a migration velocity higher than the standard, suggesting a potentially very smooth channel; η_r less than 1 indicates a migration velocity lower than the standard, suggesting possible flow obstructions or a tortuous channel.

[0147] Step S174: Calculate the arithmetic mean of the channel tortuosity connectivity index of all tracer dominant transport channels in the fracture network connectivity topology to generate a comprehensive value of fracture network tortuosity.

[0148] Traverse all nodes (i.e., all dominant migration channels) in the fracture network connectivity topology Graph_connectivity, and sum the tortuosity connectivity index τ_tracer_r for each channel calculated in step S172 to obtain a sum Sum_τ. Count the total number of nodes N_nodes. Calculate the fracture network tortuosity composite value Tau_connectivity, which is equal to Sum_τ ​​divided by N_nodes. This composite value comprehensively evaluates the tortuosity of the entire fracture network along connected paths.

[0149] Step S175: Calculate the weighted average of the transport efficiency connectivity index of all tracer dominant transport channels in the fracture network connectivity topology, and use the peak concentration value in the tracer concentration change curve of each tracer dominant transport channel as the weighting coefficient to generate the fracture network transport efficiency value.

[0150] For each node (dominant migration channel) in Graph_connectivity, the peak concentration value C_peak_j of its corresponding production well j is obtained from step S132. C_peak_j is used as the weight of this channel. The migration efficiency connectivity index η_r of each channel calculated in step S173 is multiplied by its corresponding weight C_peak_j to obtain a weighted value. The weighted values ​​of all nodes are summed to obtain the weighted sum Sum_weighted_η. At the same time, the weights C_peak_j of all nodes are summed to obtain the total weight Sum_weights. The fracture network migration efficiency value Eta_connectivity is calculated as Eta_connectivity equal to Sum_weighted_η divided by Sum_weights. This weighted average value better reflects the migration efficiency of dominant channels that contribute significantly to fluid production.

[0151] Step S176: Combine and encapsulate the comprehensive value of the tortuosity of the crack network and the value of the transport efficiency of the crack network to generate a crack connectivity quality assessment report.

[0152] Finally, the combined fracture network tortuosity value Tau_connectivity generated in step S174 and the fracture network transport efficiency value Eta_connectivity generated in step S175 are combined and encapsulated. A data structure Quality_report can be created, containing two fields: Tau_connectivity and Eta_connectivity. This report is the fracture network connectivity quality assessment report, which quantitatively describes the tortuosity and flow efficiency of the entire fracture network in terms of connectivity.

[0153] Step S180: Perform spatial distribution analysis of the fracture network in the three-dimensional morphology of the complete fracture network. Based on the coordinate sequences of all main fracture reinforcements and all secondary fracture branches in the three-dimensional morphology of the complete fracture network, calculate the maximum and minimum values ​​of all coordinate sequences in the three-dimensional coordinate system along the three coordinate axes, and determine the minimum circumscribed cuboid range of the three-dimensional morphology of the complete fracture network in space.

[0154] After generating the complete three-dimensional morphology of the fracture network, Final_network_3d, in step S1476, the following operations are performed to further analyze the spatial distribution characteristics of the fracture network. First, obtain the center point coordinates of all mesh elements from Final_network_3d. Traverse all these coordinate points and find their minimum value along the horizontal axis, denoted as X_min, and maximum value as X_max; the minimum value along the vertical axis, denoted as Y_min, and maximum value as Y_max; and the minimum value along the vertical axis, denoted as Z_min, and maximum value as Z_max. The cuboid determined by points (X_min, Y_min, Z_min) and (X_max, Y_max, Z_max) is the smallest bounding cuboid that can completely contain the entire fracture network.

[0155] Step S181: Multiply the side lengths of the minimum circumscribed cuboid along the three coordinate axes to obtain the fracture network envelope volume of the complete three-dimensional fracture network.

[0156] Based on the minimum circumscribed cuboid range determined in step S180, calculate its side lengths along the three coordinate axes. The side length Lx along the horizontal axis is equal to X_max minus X_min; the side length Ly along the vertical axis is equal to Y_max minus Y_min; and the side length Lz along the vertical axis is equal to Z_max minus Z_min. Calculate the fracture network envelope volume V_envelope, which is calculated as V_envelope equal to Lx multiplied by Ly multiplied by Lz. This volume parameter characterizes the extent of the fracture system's impact in space and is an important indicator for measuring the scale of fracturing.

[0157] Step S182: Perform spatial symmetry analysis on the three-dimensional morphology of the complete fracture network. Based on the spatial distribution of the coordinate sequences of all main fracture reinforcement and all secondary fracture branches, calculate the average value of the spatial coordinates of all coordinate points in the three coordinate axes as the coordinates of the center point of the fracture network.

[0158] Iterate through the center point coordinates of all grid cells in Final_network_3d again. Sum the x-coordinates of all points to obtain a total of Sum_X_all; sum the y-coordinates of all points to obtain a total of Sum_Y_all; sum the z-coordinates of all points to obtain a total of Sum_Z_all. Let the total number of grid cells in the crack network be N_total_cells (same as step S151). Calculate the coordinates of the crack network center point P_center, where its x-coordinate Center_X equals Sum_X_all divided by N_total_cells, its y-coordinate Center_Y equals Sum_Y_all divided by N_total_cells, and its z-coordinate Center_Z equals Sum_Z_all divided by N_total_cells.

[0159] Step S183: Using the coordinates of the center point of the fracture network as the center of symmetry, calculate the second moment of all coordinate points in the three-dimensional shape of the complete fracture network relative to the coordinates of the center point of the fracture network in the three coordinate axis directions, and generate the anisotropy matrix of the fracture network.

[0160] Using the crack network center point P_center obtained in step S182 as the reference point, the coordinates of each grid cell center point P_i in Final_network_3d are (X_i, Y_i, Z_i). The deviation vector (dX_i, dY_i, dZ_i) of P_i relative to P_center is calculated, where dX_i equals X_i minus Center_X, dY_i equals Y_i minus Center_Y, and dZ_i equals Z_i minus Center_Z. Then, a 3x3 crack network anisotropy matrix M_anisotropy is constructed. The elements of the matrix are calculated as follows: First, calculate the product of dX_i and dY_i for all points and sum them to get Sum_dXdX; calculate the product of dX_i and dY_i for all points and sum them to get Sum_dXdY; and so on, calculate Sum_dXdZ, Sum_dYdY, Sum_dYdZ, and Sum_dZdZ respectively. Then, the six independent elements of the matrix M_anisotropy (only six need to be calculated due to symmetry) are: M_11 equals Sum_dXdX divided by N_total_cells, M_12 equals Sum_dXdY divided by N_total_cells, M_13 equals Sum_dXdZ divided by N_total_cells, M_22 equals Sum_dYdY divided by N_total_cells, M_23 equals Sum_dYdZ divided by N_total_cells, and M_33 equals Sum_dZdZ divided by N_total_cells. This matrix describes the discreteness and anisotropy characteristics of the crack points distributed around their center in three-dimensional space.

[0161] Step S184: Perform eigenvalue decomposition on the anisotropy matrix of the crack network to obtain three eigenvalues ​​and corresponding eigenvectors. Use the ratio of the maximum value to the minimum value of the three eigenvalues ​​as the anisotropy coefficient of the crack network.

[0162] The anisotropy matrix M_anisotropy of the crack network generated in step S183 is decomposed using eigenvalues. The result of eigenvalue decomposition is three eigenvalues, denoted as λ1, λ2, and λ3, satisfying that λ1 is greater than or equal to λ2 and greater than or equal to λ3, along with three corresponding eigenvectors. These three eigenvalues ​​represent the dispersion of the crack network along three mutually perpendicular principal axes. The largest eigenvalue, λ1, corresponds to the direction in which the crack network extends the longest, and the smallest eigenvalue, λ3, corresponds to the direction in which the crack network distribution is most concentrated. The anisotropy coefficient A_ratio of the crack network is calculated using the formula A_ratio equals λ1 divided by λ3. A larger A_ratio indicates stronger spatial anisotropy of the crack network, meaning it mainly extends along a specific direction (the direction corresponding to λ1), while its distribution range is smaller in other directions perpendicular to that direction. A_ratio closer to 1 indicates that the spatial distribution of the crack network tends to be more isotropic, meaning it is uniformly distributed in all directions.

[0163] Step S185: Combine the crack network envelope volume and the crack network anisotropy coefficient to generate a crack network distribution descriptor.

[0164] Finally, the fracture network envelope volume V_envelope calculated in step S181 and the fracture network anisotropy coefficient A_ratio calculated in step S184 are combined and encapsulated. A data structure Spatial_extent_descriptor can be created, containing two fields: envelope volume V_envelope and anisotropy coefficient A_ratio. This descriptor is the spatial distribution descriptor of the fracture network, which quantitatively characterizes the size of the sweeping range of the pressure fracture network in space and the directional characteristics of its morphology.

[0165] For example, in step S190: the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the fracture network in the fracture network morphology data set are respectively correlated with the pre-configured target oil and gas reservoir properties, which include the reservoir effective porosity parameter and reservoir permeability parameter.

[0166] In this embodiment, after obtaining the fracture network morphology data set in step S150, the following operations are performed to investigate the influence of reservoir properties on the final fracture network morphology. First, the reservoir properties of the oil and gas reservoir to which the target fractured well area belongs are obtained. These parameters are typically obtained through core analysis, well logging interpretation, etc. For this well area, the effective porosity parameter, denoted as φ_reservoir, and the reservoir permeability parameter, denoted as K_reservoir, are obtained. These two parameters are used as independent variables. Simultaneously, four parameters from the fracture network morphology data set Fracture_network_parameters generated in step S150—the total fracture network volume parameter V_total_network, the total fracture network surface area parameter A_total_surface, the fracture network branch density parameter D_branch, and the fracture network tortuosity parameter Tau_network—are used as dependent variables. The goal is to analyze the degree of influence of the independent variables φ_reservoir and K_reservoir on each dependent variable.

[0167] Step S191: Construct a fracture network parameter prediction model using a multiple linear regression algorithm, with the reservoir effective porosity parameter and the reservoir permeability parameter as independent variables, and the fracture network total volume parameter, the fracture network total surface area parameter, the fracture network branch density parameter, and the fracture network tortuosity parameter as dependent variables.

[0168] For each dependent variable, a multiple linear regression model is constructed. Taking the total volume parameter V_total_network of the fracture network as an example, the constructed model is V_total_network equal to β0_V plus β_φ_V multiplied by φ_reservoir plus β_K_V multiplied by K_reservoir. Here, β0_V is the intercept term, β_φ_V is the weighting coefficient of porosity on volume, and β_K_V is the weighting coefficient of permeability on volume. The model construction process is based on historical fractured well data, i.e., there are multiple known sets of (φ_reservoir, K_reservoir) and corresponding V_total_network observations. Algorithms such as least squares are used to estimate the optimal β coefficients, minimizing the sum of squared errors between the model predictions and the actual observations. Similarly, multiple linear regression models are constructed to predict A_total_surface, D_branch, and Tau_network, respectively, obtaining their corresponding coefficients.

[0169] Step S192: Calculate the influence weighting coefficients of the reservoir effective porosity parameter and the reservoir permeability parameter on the total volume parameter, the total surface area parameter, the branch density parameter, and the tortuosity parameter of the fracture network according to the fracture network parameter prediction model.

[0170] For each model constructed in step S191, its regression coefficients are the weighting coefficients of the influence of the corresponding independent variable on the dependent variable. For example, for the V_total_network model, β_φ_V is the weighting coefficient of the influence of the reservoir effective porosity parameter φ_reservoir on the total fracture network volume V_total_network, and β_K_V is the weighting coefficient of the influence of the reservoir permeability parameter K_reservoir on V_total_network. Extracting these coefficients from all four models yields a set of original influence weighting coefficients.

[0171] Step S193: Normalize the influence weight coefficients to generate the influence weight value of each reservoir physical property parameter on each fracture network morphology parameter.

[0172] Since the coefficients in different models may have different dimensions and significantly different numerical ranges, the influence weight coefficients obtained in step S192 need to be normalized for easier comparison. For each dependent variable, such as V_total_network, the absolute values ​​of its two corresponding coefficients β_φ_V and β_K_V are taken, and then divided by the sum of these two absolute values ​​to obtain the normalized influence weight values. For example, the normalized influence weight W_φ_V of porosity on V_total_network is equal to |β_φ_V| divided by (|β_φ_V| plus |β_K_V|), and the normalized influence weight W_K_V of permeability on V_total_network is equal to |β_K_V| divided by (|β_φ_V| plus |β_K_V|). The sum of W_φ_V and W_K_V after this processing is 1, and their magnitudes directly reflect the relative importance of porosity and permeability in influencing V_total_network. The same normalization process was performed on the other three dependent variables (A_total_surface, D_branch, Tau_network) to obtain their respective normalized influence weight values.

[0173] Step S194: Sort the influence weight values ​​in descending order, and identify the reservoir physical property parameters that have the greatest impact on the total volume parameter of the fracture network, the reservoir physical property parameters that have the greatest impact on the total surface area parameter of the fracture network, the reservoir physical property parameters that have the greatest impact on the branch density parameter of the fracture network, and the reservoir physical property parameters that have the greatest impact on the tortuosity parameter of the fracture network.

[0174] For the total volume parameter V_total_network of the fracture network, compare the magnitudes of W_φ_V and W_K_V obtained in step S193. If W_φ_V is greater than W_K_V, then the reservoir property parameter with the greatest impact on V_total_network is considered to be the reservoir effective porosity parameter φ_reservoir; otherwise, it is the reservoir permeability parameter K_reservoir. Similarly, for A_total_surface, compare W_φ_A and W_K_A; for D_branch, compare W_φ_D and W_K_D; for Tau_network, compare W_φ_T and W_K_T. Identify the reservoir property parameter with the greatest impact on each fracture morphology parameter.

[0175] Step S195: Based on the reservoir physical property parameters that have the greatest impact on each fracture network morphology parameter, generate a list of reservoir physical property sensitive parameters to guide subsequent fracturing construction design, and associate and store the list of reservoir physical property sensitive parameters with the fracture network morphology data set.

[0176] Finally, the identification results from step S194 are summarized. For example, if φ_reservoir has the greatest impact on V_total_network, K_reservoir has the greatest impact on A_total_surface, φ_reservoir has the greatest impact on D_branch, and K_reservoir has the greatest impact on Tau_network, then the generated list of reservoir property sensitive parameters, List_sensitive_params, will contain this information: for fracturing targets aiming for larger fracture volumes, reservoir porosity should be the focus; for targets aiming for larger fracture surface areas, reservoir permeability should be the focus, etc. This list, List_sensitive_params, is then stored in association with the original fracture network morphology data set, Fracture_network_parameters. For example, fracturing fluid viscosity, sand ratio, injection rate, and other construction parameters can be adjusted specifically based on the reservoir property sensitive parameters.

[0177] In one exemplary embodiment, an oil and gas field fracture network morphology data analysis system is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2As shown, the system includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores the operating system and computer programs. The internal memory provides the environment for the operation of the operating system and computer programs in the non-volatile storage medium. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a method for analyzing fracture network morphology data in oil and gas fields. The display unit is used to generate a visually visible image and can be a display screen, projection device, or virtual reality imaging device. It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A method for analyzing fracture network morphology data in oil and gas fields, characterized in that, The method includes: Acquire a set of downhole microseismic monitoring data and an inter-well tracer monitoring dataset collected after fracturing operations in the target fracturing well area. The downhole microseismic monitoring data set includes the spatial coordinates of multiple microseismic event points recorded in chronological order, as well as the earthquake time parameters and magnitude and energy parameters corresponding to each microseismic event point. The inter-well tracer monitoring dataset includes tracer breakthrough time records and tracer concentration change curves between the injection well and the production well. The downhole microseismic monitoring data set is subjected to spatial clustering of microseismic event points. Based on the spatial proximity between the spatial coordinates of the microseismic event points and the temporal order between the earthquake occurrence time parameters, an initial spatial skeleton of the hydraulic fracture network containing multiple spatial clusters of microseismic event points is constructed. The initial spatial skeleton of the hydraulic fracture network is used to characterize the orientation and extension trend of the main fracture channels formed by hydraulic fracturing in three-dimensional space and their branching and converging relationships. Based on the initial fracture network spatial skeleton, tracer migration paths are matched on the inter-well tracer monitoring dataset. The spatial location and connectivity direction of the dominant tracer migration channels in the main fracture channels are determined by combining the tracer breakthrough time record and the tracer concentration change curve, and a fracture network connectivity topology containing the connectivity relationships of the dominant migration channels is generated. Based on the fracture network connectivity topology, the initial pressure fracture network spatial skeleton is supplemented with fracture network integrity. The spatial location corresponding to the dominant transport channel is spatially associated and fused with the main fracture channel in the initial pressure fracture network spatial skeleton to generate a complete three-dimensional shape of the pressure fracture network containing the main fracture channel and the secondary connected fracture channel. Fracture network parameters are extracted from the complete three-dimensional morphology of the fracture network to obtain a set of fracture network morphology data for the target fracture well area.

2. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The process of performing spatial clustering of microseismic event points on the downhole microseismic monitoring data set, and constructing an initial spatial framework of the fracture network containing multiple microseismic event point spatial clusters based on the spatial proximity between the spatial coordinates of the microseismic event points and the temporal order between the earthquake occurrence time parameters, includes: Spatial coordinate analysis is performed on each microseismic event point in the downhole microseismic monitoring data set to extract the three-dimensional spatial coordinate values ​​corresponding to each microseismic event point and the coordinate origin offset of the three-dimensional spatial coordinate values ​​in the geological coordinate system to which the target fracturing well area belongs. The extracted three-dimensional spatial coordinate values ​​of each microseismic event point are used as the spatial location identifier of the microseismic event point. For each microseismic event point in the downhole microseismic monitoring data set, time parameter analysis processing is performed to extract the absolute time value of the earthquake occurrence time parameter corresponding to each microseismic event point and the relative time difference of the absolute time value relative to the fracturing operation start time. The extracted absolute time value of each microseismic event point is used as the time sequence identifier of the microseismic event point. The spatial Euclidean distance between any two microseismic event points is calculated based on the spatial location identifier of each microseismic event point, and the absolute value of the time interval between any two microseismic event points is calculated based on the temporal sequence identifier of each microseismic event point, generating a microseismic event spatiotemporal correlation matrix containing a spatial Euclidean distance numerical matrix and a time interval absolute value matrix. The spatial proximity relationship of the microseismic event spatiotemporal correlation matrix is ​​filtered by calling the pre-configured spatial clustering density threshold parameter. The microseismic event point pairs whose spatial Euclidean distance value is less than the spatial clustering density threshold parameter are marked as spatially nearby event point pairs, thus obtaining an initial spatial proximity relationship set containing all spatially nearby event point pairs. For each spatial neighbor event point pair in the initial spatial neighbor relationship set, perform time continuity verification processing, extract the absolute value of the time interval corresponding to each spatial neighbor event point pair, and retain spatial neighbor event point pairs with the absolute value of the time interval less than the pre-configured time continuity threshold parameter as valid neighbor event point pairs with spatiotemporal consistency, and generate a spatiotemporally consistent neighbor relationship set. An undirected graph structure is constructed based on the effective neighbor event point pairs in the spatiotemporally consistent neighbor relationship set. The undirected graph structure uses each microseismic event point as a graph node and the spatial neighbor relationship corresponding to the effective neighbor event point pair as the connecting edge between the graph nodes to generate an initial event point association graph representing the spatiotemporal association between microseismic event points. The initial event point association graph is subjected to connectivity component detection processing to identify all interconnected graph node groups in the initial event point association graph. Each interconnected graph node group is taken as an independent microseismic event point spatial cluster, resulting in an initial cluster set containing multiple microseismic event point spatial clusters. Based on the initial cluster set, cluster feature extraction and spatial connectivity analysis are performed to construct the initial spatial skeleton of the hydraulic fracture network.

3. The method for analyzing fracture network morphology data in oil and gas fields according to claim 2, characterized in that, The step of extracting cluster features and analyzing spatial connectivity based on the initial cluster set to construct the initial spatial framework of the hydraulic fracture network includes: Extract the three-dimensional spatial coordinates of all microseismic event points in the spatial cluster of each microseismic event point, calculate the average value of the three-dimensional spatial coordinates of all microseismic event points in the spatial cluster of each microseismic event point as the spatial coordinates of the cluster center of the spatial cluster of the microseismic event point, and extract the earliest and latest earthquake occurrence times from the earthquake occurrence time parameters of all microseismic event points in the spatial cluster of each microseismic event point to calculate the time span parameter of the spatial cluster of the microseismic event point. Based on the spatial coordinates of the cluster center of each microseismic event point spatial cluster and the three-dimensional spatial coordinates of all microseismic event points in that spatial cluster, principal component analysis is used to fit the spatial principal axis direction of each microseismic event point spatial cluster. The covariance matrix of the three-dimensional spatial coordinates of all microseismic event points in that spatial cluster is calculated. The covariance matrix is ​​then decomposed into eigenvalues ​​to obtain the eigenvector corresponding to the largest eigenvalue. The spatial direction of the eigenvector corresponding to the largest eigenvalue is taken as the spatial orientation direction of the crack segment represented by the microseismic event point spatial cluster. Based on the spatial coordinates of the cluster center, spatial orientation direction, and time span parameters of each microseismic event point spatial cluster, and combined with the effective neighboring event point pairs that cross different microseismic event point spatial clusters in the spatiotemporally consistent proximity relationship set, spatial connection relationship inference processing is performed on different microseismic event point spatial clusters to determine the sequential connection order between microseismic event point spatial clusters with temporal order and spatial orientation continuity, forming a crack segment chain structure composed of multiple microseismic event point spatial clusters connected in chronological order and spatial extension direction; The spatial orientation direction and spatial coordinates of the cluster center of each microseismic event point in the chain structure of the fracture segment are used as the spatial location and extension direction identifier of the fracture segment. The spatial intersection points between different chain structures of the fracture segment are used as the branch nodes of the fracture network. Based on the serial relationship of the chain structures of the fracture segment and the branch relationship of different chain structures of the fracture segment at the spatial intersection points, an initial spatial skeleton of the pressure fracture network containing multiple chain structures of fracture segments and branch nodes is constructed.

4. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The step involves matching tracer migration paths in the inter-well tracer monitoring dataset based on the initial fracture network spatial skeleton, determining the spatial location and connectivity direction of dominant tracer migration channels in the main fracture channels by combining the tracer breakthrough time records and the tracer concentration change curves, and generating a fracture network connectivity topology that includes the connectivity relationships of dominant migration channels. The tracer breakthrough time records in the well-to-well tracer monitoring dataset are subjected to time series analysis processing to extract the first tracer breakthrough time between each injection well and each production well. The tracer migration time difference between each injection well and each production well is calculated based on the injection start time of each injection well and the first tracer breakthrough time of each production well. A tracer migration time matrix containing the tracer migration time difference between all injection wells and production wells is generated. The concentration peak feature extraction process is performed on the tracer concentration change curves in the well-to-well tracer monitoring dataset. The maximum concentration point and the corresponding occurrence time point in each tracer concentration change curve are identified. The maximum concentration value and the corresponding occurrence time point of each tracer concentration change curve are extracted as the peak concentration feature parameters of the tracer concentration change curve. The peak concentration feature parameters include the peak concentration value and the peak arrival time point value. Based on the spatial orientation of the main fracture channels and the spatial coordinates of the branch nodes in the initial fracture network spatial skeleton, a three-dimensional fracture channel grid covering the spatial extension range of all main fracture channels is generated in the geological coordinate system of the target fracture well area. The three-dimensional fracture channel grid is composed of multiple continuously arranged grid units, each grid unit corresponds to a spatial volume element and has the spatial coordinates of the grid unit center in the geological coordinate system. The tracer migration time difference between each injection well and each production well in the tracer migration time matrix is ​​mapped to the fracture channel three-dimensional mesh. Based on the injection wellhead spatial coordinates of each injection well in the geological coordinate system and the production wellhead spatial coordinates of each production well in the geological coordinate system, all possible paths that can reach the production wellhead spatial coordinates from the injection wellhead spatial coordinates through consecutive adjacent mesh cells are marked in the fracture channel three-dimensional mesh, generating multiple candidate tracer migration path sets. For each candidate tracer migration path set, the path length is calculated. The total number of consecutive adjacent grid cells contained in each candidate tracer migration path set is counted. The geometric path length of each candidate tracer migration path is calculated based on the side length of the spatial volume element corresponding to each grid cell in the three-dimensional grid of the crack channel. A candidate path length set containing the geometric path length values ​​of all candidate tracer migration paths is generated. Based on the tracer migration time difference corresponding to the tracer migration time matrix and the geometric path length of each candidate tracer migration path in the candidate path length set, the average tracer migration velocity corresponding to each candidate tracer migration path is calculated. The average tracer migration velocity is compared with a pre-configured standard migration velocity range threshold for tracers in the crack channel. Candidate tracer migration paths whose average tracer migration velocity falls within the standard migration velocity range threshold are selected as effective tracer migration paths with velocity consistency, thus obtaining a set of effective tracer migration paths. Based on the set of effective tracer transport paths, spatial identification of dominant transport channels and construction of connectivity relationships are performed to generate a fracture network connectivity topology.

5. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The method involves supplementing the initial pressure fracture network spatial skeleton based on the fracture network connectivity topology, spatially associating and fusing the spatial locations corresponding to the dominant migration channels with the main fracture channels in the initial pressure fracture network spatial skeleton, and generating a complete three-dimensional morphology of the pressure fracture network including main fracture channels and secondary connected fracture channels, including: The dominant migration channel connectivity graph in the fracture network connectivity topology is analyzed, and the spatial coordinate sequence of the continuous trajectory of the dominant migration channel corresponding to each node in the dominant migration channel connectivity graph and the channel connectivity relationship type parameter corresponding to the edge between each node are extracted. The channel connectivity relationship type parameter includes direct connectivity type and indirect connectivity type. The spatial coordinate sequence of the chain structure and branch nodes of each major fracture channel in the initial fracture network spatial skeleton is matched with the spatial coordinate sequence of each node in the dominant migration channel connection diagram. The closest point pair between the spatial coordinate sequence of each major fracture channel and the spatial coordinate sequence of the continuous trajectory of each dominant migration channel is calculated to obtain the set of spatial proximity matching relationships between the major fracture channels and the tracer dominant migration channels. For each spatial proximity matching relationship in the set of spatial proximity matching relationships, spatial overlap evaluation processing is performed. Based on the coordinate difference of the closest point pair between the main fracture channel and the tracer dominant transport channel and the side length of each grid cell in the three-dimensional mesh of the fracture channel, the overlap ratio parameter of the main fracture channel and the tracer dominant transport channel in the spatial dimension is calculated. Spatial proximity matching relationships with the overlap ratio parameter greater than the pre-configured overlap threshold are marked as spatial coincident matching relationships, and a set of spatial coincident matching relationships is generated. According to each spatial coincidence matching relationship in the set of spatial coincidence matching relationships, the main fracture channel with the tracer dominant transport channel with the spatial coincidence matching relationship is merged. Based on the spatial coordinate sequence of the main fracture channel, the coordinates of the center point of the grid cell in the spatial coordinate sequence of the continuous trajectory of the dominant transport channel that does not coincide with the spatial coordinate sequence of the main fracture channel are inserted into the corresponding position of the spatial coordinate sequence of the main fracture channel to generate the main fracture fusion enhancement coordinate sequence. The tracer dominant transport channels in the dominant transport channel connectivity graph that have connectivity with the spatially overlapping matching relationship set but do not have a spatially overlapping matching relationship with any major fracture channel are identified and processed. The tracer dominant transport channels are extracted as candidates for secondary fracture channels to be supplemented, thus obtaining a candidate set of secondary fracture channels. For each secondary crack channel candidate to be added in the secondary crack channel candidate set, spatial extension direction analysis is performed. The direction vector of the line connecting the center points of all adjacent grid cells in the spatial coordinate sequence of each secondary crack channel candidate is calculated. The main direction distribution frequency of the direction vector is counted, and the direction vector with the highest frequency is taken as the main extension direction of the secondary crack channel candidate. Based on the enhanced coordinate sequence of the main fracture fusion, the consistency of the secondary fracture direction is verified and spatial fusion is performed to generate a complete three-dimensional morphology of the pressure fracture network.

6. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The step of extracting fracture network parameters from the complete three-dimensional morphology of the fracture network to obtain a set of fracture network morphology data for the target fractured well area includes: The coordinate point traversal process is performed on all the main fracture enhancement coordinate sequences and all the secondary fracture branch coordinate sequences in the complete three-dimensional morphology of the fracture network. The total number of spatial coordinates of the center points of all grid cells in the complete three-dimensional morphology of the fracture network is counted and used as the parameter of the total number of grid cells occupied by the fracture network. Based on the actual physical size of the spatial volume element corresponding to each grid cell in the three-dimensional mesh of the fracture channel, the total number of grid cells occupied by the fracture network is multiplied by the volume value of the spatial volume element corresponding to each grid cell to calculate the total volume parameter of the fracture network occupied by the complete three-dimensional morphology of the fracture network. Surface mesh element identification processing is performed on all main fracture enhancement coordinate sequences and all secondary fracture branch coordinate sequences in the complete three-dimensional morphology of the fracture network. Mesh elements located at the edge of the fracture network and having at least one adjacent face not shared with any other fracture network mesh element in each spatial coordinate sequence are marked as fracture network surface mesh elements. The total number of fracture network surface mesh elements is counted as the number of fracture network surface elements. Based on the actual physical dimensions of the spatial volume element corresponding to each grid cell in the three-dimensional mesh of the fracture channel, the total area of ​​the exposed surface of each grid cell marked as the surface of the fracture network is calculated; the total exposed surface areas of all grid cells of the fracture network are summed to obtain the total surface area parameter of the fracture network exposed by the complete three-dimensional morphology of the fracture network. All branch nodes in the complete three-dimensional morphology of the fracture network are identified and statistically processed. The connection points of the coordinate sequences of the main fracture channels and secondary fracture branches in the complete three-dimensional morphology of the fracture network are extracted as fracture network branch nodes. The total number of fracture network branch nodes is counted as the number of fracture network branch nodes. Based on the complete three-dimensional morphology of the fracture network, the branch density and tortuosity parameters of the fracture network are calculated to generate a set of fracture network morphology data.

7. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, After extracting fracture network parameters from the complete three-dimensional morphology of the fracture network to obtain the fracture network morphology data set of the target fractured well area, the method further includes: The total volume parameters, total surface area parameters, branch density parameters, and tortuosity parameters of the fracture network in the fracture network morphology data set are input into the pre-training model of the fracturing effect. The pre-training model of the fracturing effect includes multiple sequentially connected convolutional layers, pooling layers, and fully connected layers. In the pre-trained model of the fracturing effect, the first convolutional layer performs one-dimensional convolutional feature extraction on the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the fracture network, generating a first convolutional feature map containing local correlation features of the four input parameters. The first convolutional feature map is subjected to max pooling downsampling through the first pooling layer to reduce the feature dimension of the first convolutional feature map while retaining the main feature responses, thereby generating the first pooling feature map. The first pooling feature map is extracted by deep one-dimensional convolution through the second convolutional layer to capture the high-order interaction relationship between the total volume parameter, total surface area parameter, branch density parameter, and tortuosity parameter of the crack network, and to generate the second convolutional feature map. The second convolutional feature map is subjected to global average pooling through the second pooling layer, which averages all feature values ​​of each feature channel in the second convolutional feature map to generate a one-dimensional global feature vector. The one-dimensional global feature vector is input into the fully connected layer of the fracturing effect pre-training model. The one-dimensional global feature vector is then subjected to linear transformation and nonlinear activation processing through the fully connected layer to calculate the comprehensive fracturing effect score parameters of the target fracturing well area. The fracturing effect level of the target fracturing well area is determined by comparing the comprehensive fracturing effect scoring parameters with the preset fracturing effect level threshold, and a fracturing effect evaluation result containing the fracturing effect level and the comprehensive fracturing effect scoring parameters is generated. The fracturing effect evaluation results are associated and stored with the fracturing fracture network morphology data set, and a visualization report containing the fracturing effect evaluation results is generated. The visualization report includes a three-dimensional rendering of the complete three-dimensional fracturing fracture network morphology and a labeling of the fracturing effect level.

8. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The step of matching tracer migration paths in the inter-well tracer monitoring dataset based on the initial fracture network spatial skeleton, determining the spatial location and connectivity direction of the dominant tracer migration channels in the main fracture channels by combining the tracer breakthrough time records and the tracer concentration change curves, and generating a fracture network connectivity topology containing the connectivity relationships of the dominant migration channels, further includes: The connectivity quality of the fracture network topology is assessed, and the coordinate sequence of the grid center point of the continuous trajectory of each dominant transport channel in the fracture network connectivity topology is extracted. The sum of the Euclidean distances between the center points of all adjacent grid cells in the continuous trajectory of each dominant transport channel is calculated as the actual transport path length of the tracer's dominant transport channel. The straight-line Euclidean distance between the spatial coordinates of the injection wellhead and the spatial coordinates of the production wellhead corresponding to the continuous trajectory of each dominant migration channel is taken as the straight-line connection distance of the dominant migration channel of the tracer. The channel tortuosity connectivity index of each tracer dominant transport channel is calculated by dividing the actual transport path length of each tracer dominant transport channel by the straight-line connectivity distance of that tracer dominant transport channel. Based on the first breakthrough time of the tracer in the breakthrough time record corresponding to each tracer dominant transport channel, the average transport velocity of the tracer corresponding to the dominant transport channel is calculated. The average transport velocity of the tracer is divided by the pre-configured standard transport velocity of the tracer to obtain the transport efficiency connectivity index of each tracer dominant transport channel. The channel tortuosity connectivity index of all tracer dominant transport channels in the fracture network connectivity topology is calculated by arithmetic mean to generate a comprehensive value of fracture network tortuosity. The transport efficiency connectivity index of all tracer dominant transport channels in the fracture network connectivity topology is calculated by weighted average. The peak concentration value in the tracer concentration change curve of each tracer dominant transport channel is used as the weighting coefficient to generate the fracture network transport efficiency value. The combined value of the fracture network tortuosity and the fracture network migration efficiency value are combined and packaged to generate a fracture connectivity quality assessment report.

9. The method for analyzing fracture network morphology data in oil and gas fields according to claim 1, characterized in that, The step of supplementing the initial pressure fracture network spatial skeleton with fracture network integrity based on the fracture network connectivity topology, spatially associating and fusing the spatial location corresponding to the dominant migration channel with the main fracture channel in the initial pressure fracture network spatial skeleton, and generating a complete three-dimensional morphology of the pressure fracture network including the main fracture channel and secondary connected fracture channels, further includes: A spatial distribution analysis of the fracture network is performed on the three-dimensional morphology of the complete fracture network. Based on the coordinate sequences of all main fracture reinforcements and all secondary fracture branches in the three-dimensional morphology of the complete fracture network, the maximum and minimum values ​​of all coordinate sequences in the three-dimensional coordinate system are calculated in the three-dimensional coordinate system, and the minimum circumscribed cuboid range of the three-dimensional morphology of the complete fracture network in space is determined. The product of the side lengths of the minimum circumscribed cuboid along the three coordinate axes is used to calculate the fracture network envelope volume of the complete three-dimensional fracture network. A spatial symmetry analysis of the complete three-dimensional fracture network was performed. Based on the spatial distribution of the coordinate sequences of all main fracture reinforcement and all secondary fracture branches, the average value of the spatial coordinates of all coordinate points in the three coordinate axes was calculated as the coordinates of the center point of the fracture network. Using the coordinates of the center point of the fracture network as the center of symmetry, calculate the second moments of all coordinate points in the three-dimensional shape of the complete fracture network with respect to the coordinates of the center point of the fracture network in the three coordinate axes, and generate the anisotropy matrix of the fracture network. The anisotropy matrix of the crack network is subjected to eigenvalue decomposition to obtain three eigenvalues ​​and corresponding eigenvectors. The ratio of the maximum value to the minimum value of the three eigenvalues ​​is used as the anisotropy coefficient of the crack network. The crack network envelope volume and the crack network anisotropy coefficient are combined to generate a crack network distribution descriptor.

10. A data analysis system for fracture network morphology in oil and gas fields, characterized in that, include: processor; A machine-readable storage medium for storing machine-executable instructions of the processor; The processor is configured to execute the oil and gas field fracture network morphology data analysis method according to any one of claims 1 to 9 by executing the machine-executable instructions.

Images