Oil reservoir fracture identification method based on deep learning

Abstract

The invention discloses an oil reservoir fracture identification method based on deep learning, and the method comprises the following steps: S1, receiving earthquake or underground imaging data, and carrying out the filtering, time-depth conversion and multi-scale wavelet denoising to generate a multi-scale feature body; s2, based on the feature body, constructing a crack characterization tensor containing a spatial gradient channel, an anisotropic channel, a trend sensitive channel and a similarity channel; s3, inputting the representation tensor into a cross-layer attention depth network to fuse high and low layer features, and outputting a crack probability graph and a trend tensor; s4, constructing a path graph based on the probability and the trend, and generating a crack master control path by using a graph neural network; s5, performing disturbance on the main control path to generate an anti-fact path; and S6, outputting a structured crack model containing crack trend, length, connectivity and density through comparative learning. According to the method, automatic identification and interpretable modeling of a complex oil reservoir fracture system are realized, and the fracture quantitative characterization precision is improved.

Description

Technical Field

[0001] This invention relates to the fields of oil and gas exploration geophysics and deep learning, and in particular to a deep learning-based method for identifying reservoir fractures. Background Technology

[0002] In oil and gas exploration and development, fractures are a crucial factor affecting reservoir permeability, reservoir performance, and oil and gas production capacity. The accuracy of fracture identification directly determines the reliability of geological modeling and development planning. Currently, fracture identification methods mainly rely on manual interpretation or geometric attribute volume analysis, commonly using seismic attributes such as seismic coherence volume, strike consistency volume, curvature volume, and quasi-fracture volume. Geological experts manually delineate the fracture direction and distribution on profiles or volume views. However, manual methods are not only inefficient and highly subjective, but also prone to false detections and missed detections in areas with large variations in fracture direction, dense intersections, and high noise levels. Furthermore, it is difficult to quantify key parameters such as fracture direction, length, connectivity, and density, often limiting fracture modeling results to qualitative interpretation and failing to meet the needs of refined reservoir description.

[0003] To improve automation, some methods employ threshold segmentation based on geometric attributes such as strike consistency, coherence, and curvature for crack detection in seismic volumes. While these methods perform well for simple cracks, they exhibit poor connectivity, strike discontinuities, and numerous false cracks in complex crack zones. Furthermore, they rely on a single-scale attribute volume, failing to simultaneously characterize cracks of different sizes. In addition, when the signal-to-noise ratio of seismic data is low or crack strikes change abruptly, traditional attribute-based methods are prone to misidentifying noise anomalies as cracks or missing actual cracks, making it difficult to guarantee the stability and reliability of the identification results.

[0004] In recent years, deep learning-based crack detection methods have gradually emerged, outputting crack probability maps by inputting seismic profiles into convolutional neural networks such as U-Net and FCN. Compared with traditional methods, deep learning has shown better crack enhancement capabilities in complex background noise. However, most studies are limited to shallow or single-scale networks, failing to simultaneously learn the global morphology and local texture features of cracks, especially at abrupt changes in crack direction or branches where directional fractures are prone to occur. Furthermore, existing methods only output probabilistic volumes or binary segmentation maps, lacking the ability to model the main crack control path, branch structure, and connectivity, making it difficult to automatically generate structured crack models containing quantified information on direction, length, connectivity, and density.

[0005] In fracture path modeling, existing technologies often employ simple connected component algorithms to merge high-probability regions within the probability volume, lacking path planning strategies based on joint optimization of fracture orientation consistency and node probabilities. They also fail to incorporate graph neural networks to update fracture node features and plan the dominant fracture path. Regarding the issue of false fractures in fracture identification, current methods have not introduced counterfactual sample generation mechanisms to compare and learn positive and negative samples, thus limiting the model's ability to identify noise and false fractures. Furthermore, most current fracture detection methods lack interpretable output of fracture identification results, lacking orientation distribution, length distribution, and connectivity indicators, failing to meet the needs of geological modeling, reservoir evaluation, and production deployment for quantitative fracture models.

[0006] Therefore, how to provide a deep learning-based method for reservoir fracture identification is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0007] One objective of this invention is to propose a deep learning-based method for reservoir fracture identification. This invention acquires and receives seismic or downhole imaging data, processes it to obtain multi-scale feature volumes, and then constructs a fracture representation tensor based on these feature volumes, incorporating spatial gradients, anisotropy, strike sensitivity, and similarity channels. This representation tensor is then input into a cross-layer attention deep network to fuse high- and low-layer features, outputting a fracture probability map and strike tensor. A graph neural network is used to generate the dominant fracture path, and perturbations are applied to this path to generate counterfactual paths. Finally, through comparative learning, a structured fracture model containing fracture strike, length, connectivity, and density is output. This invention achieves automatic identification and interpretable modeling of complex reservoir fracture systems, improving the accuracy of quantitative fracture characterization.

[0008] A method for path planning in a discrete patch traversal survey under the influence of multiple factors, according to an embodiment of the present invention, includes the following steps: S1. Receive seismic data or downhole imaging data, perform adaptive predictive filtering on the received data to suppress strip noise, remove isolated noise points through morphological filtering, perform time-depth conversion based on velocity model on the filtered data to convert time domain data into depth domain data, and use symmetric orthogonal wavelet basis to perform three-level multi-scale wavelet denoising processing on the depth domain data to generate multi-scale feature bodies. S2. Construct a crack representation tensor based on multi-scale feature volumes, including calculating the gradients of each voxel in the three spatial directions to form a spatial gradient tensor, using an elliptic anisotropy operator to calculate the anisotropic properties of voxels to generate an anisotropic tensor, using multiple directional orientation-sensitive convolution kernels covering the range of 0° to 180° to calculate the response values ​​of voxels in each direction to generate orientation-sensitive tensors, calculating the local similarity value of each voxel based on the gray-level similarity of the voxel neighborhood to generate a similarity tensor, and concatenating the tensors in the channel dimension to obtain the crack representation tensor. S3. Input the crack representation tensor into the cross-layer attention depth network. The cross-layer attention depth network includes a high-level crack morphology perception branch and a low-level crack texture perception branch. The high-level and low-level features are fused through bidirectional cross-layer attention matrix interaction to output the crack probability map and crack orientation tensor. S4. Construct a crack path planning graph based on the crack probability graph and crack orientation tensor. Use voxels with crack probability values ​​higher than a preset threshold as graph nodes. Calculate edge weights based on the orientation consistency between nodes, the mean crack probability, and the spatial geometric distance. Construct edges between nodes. Use a graph neural network to update node features multiple times. Use the shortest path search algorithm based on edge weights to generate the main crack control path. S5. Perform structural disturbance operations on the main crack control path to generate a set of counterfactual crack paths, including perturbation of the order of path nodes, perturbation of the direction of path nodes, perturbation of the path length, perturbation of the path nodes by interpolation, and filtering out counterfactual paths with morphological similarity higher than a preset threshold. S6. Input the main control path and the counterfactual crack path into the comparative embedding encoder, and embed the path node sequence, crack direction tensor and path connectivity vector into the encoder respectively. Concatenate them to form a path-level embedding vector. Maximize the embedding distance between the main control path and the counterfactual path through comparative learning, and output a structured crack distribution model containing crack direction, length, connectivity and crack density information.

[0009] Optionally, S1 specifically includes: S11. Perform adaptive prediction filtering on the received seismic data or downhole imaging data, divide the input data into blocks according to the gather direction, adaptively adjust the filtering parameters based on the average frequency drift within each block, apply a prediction error filter to the data of each block, and generate preliminary denoised data. S12. Perform morphological filtering on the data generated in S11. Use a 3×3 structuring element to perform grayscale expansion and erosion operations respectively to eliminate isolated high-frequency noise and smooth small-scale random noise. S13. Based on the regional velocity model, perform time-depth conversion on the data generated in S12, and convert the time-domain samples of each profile... Convert to depth domain samples The conversion formula is: in, Indicates time Longitudinal wave velocity at that location For depth location; S14. Perform three-level multi-scale wavelet decomposition on the depth domain data generated in S13, using a symmetric orthogonal wavelet basis. Decomposed into low-frequency approximate components and high-frequency detail components The number of decomposition layers The multi-scale wavelet decomposition formula is: in, and These are the low-pass and high-pass filter coefficients of the wavelet basis, respectively; S15. Based on the decomposition results of S14, a threshold denoising method is used to denoise each high-frequency component. Soft thresholding is performed, and the denoised low-frequency feature components are reconstructed. ; S16. The reconstruction result of S15 is weighted and fused with the original depth domain data to form a multi-scale feature volume containing multi-resolution information, which serves as the input for the subsequent crack characterization tensor construction.

[0010] Optionally, S2 specifically includes: S21. Calculate the gradient value of each voxel in the three-dimensional space direction for the multi-scale feature volume generated in S1, and obtain the gradient components in the X, Y, and Z directions using the Sobel operator. Calculate the spatial gradient magnitude The calculation formula is as follows: in, These represent the gradient components of the voxel in the X, Y, and Z directions, respectively. Indicates the magnitude of the spatial gradient; S22. Based on the gradient information from S21, the anisotropy of voxels is calculated using the elliptic anisotropy operator. The Based on the ratio of the major and minor radii of the principal axes of voxels The definition and calculation formula are as follows: in, , These are the largest and second largest eigenvalues ​​of the voxel gray-level covariance matrix; S23. Apply orientation-sensitive convolution kernels covering the range of 0° to 180° to the multi-scale feature volume generated in S1, with K directions. Calculate the voxel response value for each direction. This forms a direction-sensitive tensor, where each direction is 15° when K=12; S24. The mean square error of gray level in the 3×3×3 neighborhood of each voxel in the multi-scale feature volume generated based on S1. Calculate local similarity of voxels And normalize it. The formula is defined as follows: in, Represents the similarity value. The gray-level variance within the neighborhood; S25. The spatial gradient tensor generated in S21 The anisotropic tensor generated by S22 S23 generated direction-sensitive tensor Similarity tensor generated by S24 The cracks are spliced ​​along the channel dimension to form a crack representation tensor containing multiple crack features, which serves as the input for subsequent S3 processing.

[0011] Optionally, S3 specifically includes: S31. Input the crack representation tensor generated in S2 into the high-level crack morphology perception branch of the cross-layer attention deep network, and extract the global crack morphology features through five three-dimensional convolutional layers in sequence. The kernel size of each convolutional layer is 3×3×3, the stride is 1, and the receptive field covers a crack range of 50m. S32. Simultaneously input the crack representation tensor generated in S2 into the low-level crack texture perception branch of the cross-layer attention deep network. Extract local crack texture details through three 3D convolutional layers and multi-scale pyramid pooling. Pyramid pooling generates feature tensors at different scales. ; S33, the feature tensors output from the higher layer of S31 respectively and the feature tensor output by the lower layer of S32 Calculate the cross-layer bidirectional attention matrix Note the matrix The calculation formula is: in, Represents the high-level feature tensor. This represents the low-level feature tensor. For learnable weight matrix, Represents a bidirectional cross-layer attention matrix; S34, will Acting on and The system uses matrix multiplication to perform interactive updates of high-level and low-level features, and outputs the updated high-level features. and low-level features ; S35, will and The data is concatenated along the channel dimension and then integrated using a 3D convolutional fusion network to output a crack probability map. and crack orientation tensor ,in This represents the probability value of each voxel belonging to a crack. This indicates the angle of the crack direction for each voxel.

[0012] Optionally, S4 specifically includes: S41. Crack probability map based on S3 output and crack orientation tensor Select crack probability value Greater than the set threshold voxels as candidate nodes for cracks ,in The probability threshold for crack voxels; S42, For each pair of adjacent nodes Based on the consistency of the direction between nodes Mean probability of nodal cracks Euclidean distance between nodes Calculate edge weights The formula for calculating the edge weight is: in, This indicates that nodes Nᵢ and Nⱼ follow the same direction. This represents the mean probability of cracks at two nodes. For Euclidean distance, , , These are weighting coefficients; S43, based on and Constructing a crack voxel map ,in For a set of nodes, Let it be the set of edges; S44, in the crack voxel diagram above The layered graph neural network updates the node state, and the feature update formula for each layer is as follows: in, Represents a node In the Layer feature representation, Pay attention to weights between nodes. It is a learnable transformation matrix; S45, in the updated diagram Based on the shortest path search algorithm, using a weighted sum of path consistency and node probability as the cost function, the algorithm searches for the optimal master crack path connecting the start and end nodes, and outputs the master crack path sequence. .

[0013] Optionally, S5 specifically includes: S51, The main control path of the crack generated in S4 Perform node order misalignment perturbation, proportionally adjusting the path node sequence. Random misalignment, This is a misaligned proportional parameter; S52, to Perform a path offset perturbation to change the path of each node in the path. By maximum offset angle exist Random disturbances within the range This is the offset angle threshold; S53, to Execute path length truncation perturbation to trim the path by a certain percentage from either the beginning or end. The node, This is the truncation ratio parameter; S54, to Perform node interpolation perturbation in the node sequence with average node spacing. Insert a new node to make the smoothness of the disturbance path approximate that of the master path, where for Average Euclidean distance between mid-nodes; S55. Calculate each disturbance path and morphological similarity The formula for morphological similarity is: in, for and The distance to Hausdorff; S56, will Greater than the set threshold The perturbation paths are preserved, forming a set of counterfactual crack paths with reasonable morphology. This is used for subsequent comparative learning processes.

[0014] Optionally, S6 specifically includes: S61, Convert the crack master control path generated in S4 The set of counterfactual paths generated by S5 Each counterfactual path in Input the path feature encoder separately to the path node sequence Crack orientation sequence and path connectivity vector Encode the nodes to generate node embedding vectors. Towards Embedded Vectors and connectivity embedding vectors ; S62, will , , Path-level embedding vectors are formed by concatenating the data along the channel dimension. ,in S63. Employing a contrastive learning network to... Training is performed to construct triplet samples. ,in Embedded in the master control path Positive samples for the main control path. For negative samples of the counterfactual path, a triplet loss function is used during training, and the loss formula is: in This represents the Euclidean distance between embedding vectors. This refers to the boundary margin hyperparameter; S64. Based on the trained contrastive learning network, the master path... and The model is input via a mid-path, and the output includes the crack orientation. Crack length Crack connectivity and crack density Structured crack distribution model It is used for the quantitative interpretation and visualization of crack systems.

[0015] The beneficial effects of this invention are: (1) It can automatically achieve accurate identification of complex reservoir fracture systems, significantly improving the quantitative characterization accuracy of geometric features such as fracture orientation, length, connectivity and density; (2) By using multi-scale feature extraction and cross-layer attention network, high and low layer information are integrated, taking into account both the global morphology of cracks and local texture details, thus improving the adaptability to cracks of different scales. (3) By using the crack path planning method based on graph neural network, the continuity of crack direction and the coherence of path can be effectively guaranteed, and a reliable crack main control path can be generated. (4) The introduction of counterfactual perturbation samples and contrastive learning mechanism enhances the model’s ability to distinguish noise and false cracks, and improves the robustness and generalization ability of crack identification results. (5) It can output structured fracture models containing quantitative indicators such as fracture orientation, length, connectivity, and density, providing accurate and interpretable fracture information for geological modeling, reservoir evaluation, and development design; (6) It can reduce the reliance on manual interpretation and improve the automation level and overall work efficiency of crack identification. Attached Figure Description

[0016] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is an overall flowchart of a deep learning-based reservoir fracture identification method proposed in this invention. Figure 2 This diagram illustrates the comparative learning of the master control path and counterfactual fracture paths in a deep learning-based reservoir fracture identification method. It includes the master control path and multiple counterfactual paths generated through node sequence misalignment, strike offset, length truncation, and node interpolation. It also demonstrates the distance calculation process between the master control path and the counterfactual paths by the comparative learning network. Detailed Implementation

[0017] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0018] refer to Figure 1 and Figure 2 A deep learning-based method for reservoir fracture identification includes the following steps: S1. Receive seismic data or downhole imaging data, perform adaptive predictive filtering on the received data to suppress strip noise, remove isolated noise points through morphological filtering, perform time-depth conversion based on velocity model on the filtered data to convert time domain data into depth domain data, and use symmetric orthogonal wavelet basis to perform three-level multi-scale wavelet denoising processing on the depth domain data to generate multi-scale feature bodies. S2. Construct a crack representation tensor based on multi-scale feature volumes, including calculating the gradients of each voxel in the three spatial directions to form a spatial gradient tensor, using an elliptic anisotropy operator to calculate the anisotropic properties of voxels to generate an anisotropic tensor, using multiple directional orientation-sensitive convolution kernels covering the range of 0° to 180° to calculate the response values ​​of voxels in each direction to generate orientation-sensitive tensors, calculating the local similarity value of each voxel based on the gray-level similarity of the voxel neighborhood to generate a similarity tensor, and concatenating the tensors in the channel dimension to obtain the crack representation tensor. S3. Input the crack representation tensor into the cross-layer attention depth network. The cross-layer attention depth network includes a high-level crack morphology perception branch and a low-level crack texture perception branch. The high-level and low-level features are fused through bidirectional cross-layer attention matrix interaction to output the crack probability map and crack orientation tensor. S4. Construct a crack path planning graph based on the crack probability graph and crack orientation tensor. Use voxels with crack probability values ​​higher than a preset threshold as graph nodes. Calculate edge weights based on the orientation consistency between nodes, the mean crack probability, and the spatial geometric distance. Construct edges between nodes. Use a graph neural network to update node features multiple times. Use the shortest path search algorithm based on edge weights to generate the main crack control path. S5. Perform structural disturbance operations on the main crack control path to generate a set of counterfactual crack paths, including perturbation of the order of path nodes, perturbation of the direction of path nodes, perturbation of the path length, perturbation of the path nodes by interpolation, and filtering out counterfactual paths with morphological similarity higher than a preset threshold. S6. Input the main control path and the counterfactual crack path into the comparative embedding encoder, and embed the path node sequence, crack direction tensor and path connectivity vector into the encoder respectively. Concatenate them to form a path-level embedding vector. Maximize the embedding distance between the main control path and the counterfactual path through comparative learning, and output a structured crack distribution model containing crack direction, length, connectivity and crack density information.

[0019] This invention acquires and receives seismic or downhole imaging data, processes it to obtain multi-scale feature volumes, and then constructs a fracture characterization tensor based on these feature volumes, incorporating spatial gradients, anisotropy, strike sensitivity, and similarity channels. This characterization tensor is then input into a cross-layer attention depth network to fuse high- and low-layer features, outputting a fracture probability map and strike tensor. A graph neural network is used to generate the dominant fracture path, and perturbations are applied to this path to generate counterfactual paths. Finally, through comparative learning, a structured fracture model containing fracture strike, length, connectivity, and density is output. This invention achieves automatic identification and interpretable modeling of complex reservoir fracture systems, improving the accuracy of quantitative fracture characterization. In this embodiment, S1 specifically includes: S11. Perform adaptive prediction filtering on the received seismic data or downhole imaging data, divide the input data into blocks according to the gather direction, adaptively adjust the filtering parameters based on the average frequency drift within each block, apply a prediction error filter to the data of each block, and generate preliminary denoised data. S12. Perform morphological filtering on the data generated in S11. Use a 3×3 structuring element to perform grayscale expansion and erosion operations respectively to eliminate isolated high-frequency noise and smooth small-scale random noise. S13. Based on the regional velocity model, perform time-depth conversion on the data generated in S12, and convert the time-domain samples of each profile... Convert to depth domain samples The conversion formula is: in, Indicates time Longitudinal wave velocity at that location For depth location; S14. Perform three-level multi-scale wavelet decomposition on the depth domain data generated in S13, using a symmetric orthogonal wavelet basis. Decomposed into low-frequency approximate components and high-frequency detail components The number of decomposition layers The multi-scale wavelet decomposition formula is: in, and These are the low-pass and high-pass filter coefficients of the wavelet basis, respectively; S15. Based on the decomposition results of S14, a threshold denoising method is used to denoise each high-frequency component. Soft thresholding is performed, and the denoised low-frequency feature components are reconstructed. ; S16. The reconstruction result of S15 is weighted and fused with the original depth domain data to form a multi-scale feature volume containing multi-resolution information, which serves as the input for the subsequent crack characterization tensor construction.

[0020] In this embodiment, S2 specifically includes: S21. Calculate the gradient value of each voxel in the three-dimensional space direction for the multi-scale feature volume generated in S1, and obtain the gradient components in the X, Y, and Z directions using the Sobel operator. Calculate the spatial gradient magnitude The calculation formula is as follows: in, These represent the gradient components of the voxel in the X, Y, and Z directions, respectively. Indicates the magnitude of the spatial gradient; S22. Based on the gradient information from S21, the anisotropy of voxels is calculated using the elliptic anisotropy operator. The Based on the ratio of the major and minor radii of the principal axes of voxels The definition and calculation formula are as follows: in, , These are the largest and second largest eigenvalues ​​of the voxel gray-level covariance matrix; S23. Apply orientation-sensitive convolution kernels covering the range of 0° to 180° to the multi-scale feature volume generated in S1, with K directions. Calculate the voxel response value for each direction. This forms a direction-sensitive tensor, where each direction is 15° when K=12; S24. The mean square error of gray level in the 3×3×3 neighborhood of each voxel in the multi-scale feature volume generated based on S1. Calculate local similarity of voxels And normalize it. The formula is defined as follows: in, Represents the similarity value. The gray-level variance within the neighborhood; S25. The spatial gradient tensor generated in S21 The anisotropic tensor generated by S22 S23 generated direction-sensitive tensor Similarity tensor generated by S24 The cracks are spliced ​​along the channel dimension to form a crack representation tensor containing multiple crack features, which serves as the input for subsequent S3 processing.

[0021] In this embodiment, S3 specifically includes: S31. Input the crack representation tensor generated in S2 into the high-level crack morphology perception branch of the cross-layer attention deep network, and extract the global crack morphology features through five three-dimensional convolutional layers in sequence. The kernel size of each convolutional layer is 3×3×3, the stride is 1, and the receptive field covers a crack range of 50m. S32. Simultaneously input the crack representation tensor generated in S2 into the low-level crack texture perception branch of the cross-layer attention deep network. Extract local crack texture details through three 3D convolutional layers and multi-scale pyramid pooling. Pyramid pooling generates feature tensors at different scales. ; S33, the feature tensors output from the higher layer of S31 respectively and the feature tensor output by the lower layer of S32 Calculate the cross-layer bidirectional attention matrix Note the matrix The calculation formula is: in, Represents the high-level feature tensor. This represents the low-level feature tensor. For learnable weight matrix, Represents a bidirectional cross-layer attention matrix; S34, will Acting on and The system uses matrix multiplication to perform interactive updates of high-level and low-level features, and outputs the updated high-level features. and low-level features ; S35, will and The data is concatenated along the channel dimension and then integrated using a 3D convolutional fusion network to output a crack probability map. and crack orientation tensor ,in This represents the probability value of each voxel belonging to a crack. This indicates the angle of the crack direction for each voxel.

[0022] In this embodiment, S4 specifically includes: S41. Crack probability map based on S3 output and crack orientation tensor Select crack probability value Greater than the set threshold voxels as candidate nodes for cracks ,in The probability threshold for crack voxels; S42, For each pair of adjacent nodes Based on the consistency of the direction between nodes Mean probability of nodal cracks Euclidean distance between nodes Calculate edge weights The formula for calculating the edge weight is: in, This indicates that nodes Nᵢ and Nⱼ follow the same direction. This represents the mean probability of cracks at two nodes. For Euclidean distance, , , These are weighting coefficients; S43, based on and Constructing a crack voxel map ,in For a set of nodes, Let it be the set of edges; S44, in the crack voxel diagram above The layered graph neural network updates the node state, and the feature update formula for each layer is as follows: in, Represents a node In the Layer feature representation, Pay attention to weights between nodes. It is a learnable transformation matrix; S45, in the updated diagram Based on the shortest path search algorithm, using a weighted sum of path consistency and node probability as the cost function, the algorithm searches for the optimal master crack path connecting the start and end nodes, and outputs the master crack path sequence. .

[0023] In this embodiment, S5 specifically includes: S51, The main control path of the crack generated in S4 Perform node order misalignment perturbation, proportionally adjusting the path node sequence. Random misalignment, This is a misaligned proportional parameter; S52, to Perform a path offset perturbation to change the path of each node in the path. By maximum offset angle exist Random disturbances within the range This is the offset angle threshold; S53, to Execute path length truncation perturbation to trim the path by a certain percentage from either the beginning or end. The node, This is the truncation ratio parameter; S54, to Perform node interpolation perturbation in the node sequence with average node spacing. Insert a new node to make the smoothness of the disturbance path approximate that of the master path, where for Average Euclidean distance between mid-nodes; S55. Calculate each disturbance path and morphological similarity The formula for morphological similarity is: in, for and The distance to Hausdorff; S56, will Greater than the set threshold The perturbation paths are preserved, forming a set of counterfactual crack paths with reasonable morphology. This is used for subsequent comparative learning processes.

[0024] In this embodiment, S6 specifically includes: S61, Convert the crack master control path generated in S4 The set of counterfactual paths generated by S5 Each counterfactual path in Input the path feature encoder separately to the path node sequence Crack orientation sequence and path connectivity vector Encode the nodes to generate node embedding vectors. Towards Embedded Vectors and connectivity embedding vectors ; S62, will , , Path-level embedding vectors are formed by concatenating the data along the channel dimension. ,in S63. Employing a contrastive learning network to... Training is performed to construct triplet samples. ,in Embedded in the master control path Positive samples for the main control path. For negative samples of the counterfactual path, a triplet loss function is used during training, and the loss formula is: in This represents the Euclidean distance between embedding vectors. This refers to the boundary margin hyperparameter; S64. Based on the trained contrastive learning network, the master path... and The model is input via a mid-path, and the output includes the crack orientation. Crack length Crack connectivity and crack density Structured crack distribution model It is used for the quantitative interpretation and visualization of crack systems.

[0025] By receiving seismic or downhole imaging data and performing filtering, time-depth conversion, and multi-scale wavelet denoising, a fracture feature volume that takes into account multi-resolution characteristics is generated. Based on this, a fracture representation tensor integrating spatial gradient, anisotropy, orientation sensitivity, and similarity information is constructed. This representation tensor is then input into a cross-layer attention depth network to extract global fracture morphology and local texture features, outputting a fracture probability map and orientation tensor. A path graph is then constructed based on the probability and orientation information, and a graph neural network is used to plan the dominant fracture path. Perturbations such as node misalignment, orientation shift, length truncation, and node interpolation are performed to generate counterfactual paths. Finally, contrastive learning is used to maximize the discriminative power between the dominant path and the counterfactual path, outputting a structured fracture model containing fracture orientation, length, connectivity, and density. This method achieves automatic modeling of fracture continuity and orientation consistency, significantly improving the accuracy and stability of fracture identification results. It can effectively distinguish between noise and pseudo-fractures, while automatically outputting quantified fracture information, providing reliable support for geological modeling, reservoir evaluation, and well placement, reducing reliance on manual interpretation, and improving work efficiency.

[0026] Example 1: To verify the feasibility of this invention in practice, it was applied to a carbonate reservoir in the Tarim Basin of China. This block is located in the Tarim North Uplift area, with reservoir depths ranging from 4000 to 4800 meters. The area is characterized by well-developed and complex fracture morphology, exhibiting features such as intersecting strikes, high fracture density, and association with faults. Simultaneously, the seismic data has a low signal-to-noise ratio, with a trace spacing of 12.5 meters, 1200 traces, a 3D seismic volume coverage area of ​​50 square kilometers, over 400 million sample points, and a maximum sampling window of 3500 milliseconds. In this block, manual fracture identification methods based on coherence and curvature volumes generally suffer from discontinuous fracture strike descriptions, inaccurate fracture connectivity assessments, and high sensitivity to noise. This has led to repeated misjudgments of fracture connectivity during well placement, resulting in unsatisfactory subsequent production measures. Therefore, this block became an ideal scenario for verifying the effectiveness of the method of this invention.

[0027] In the application, the block-based 3D seismic volume is first used as input data. Adaptive predictive filtering removes strip noise, morphological filtering removes isolated high-frequency noise points, and time-depth conversion is performed based on the velocity model to transform the time-domain data into the depth domain, with a maximum depth of 4700 meters. Subsequently, the depth-domain data is subjected to three-level multi-scale wavelet decomposition and denoising using the symmetric orthogonal wavelet basis db4 to generate multi-scale feature volumes. Based on this, the spatial gradient magnitude of each voxel in the X, Y, and Z directions is calculated, and the anisotropy attribute is calculated using the elliptic anisotropy operator. Orientation-sensitive convolution kernel operations are performed on 12 directions within the range of 0° to 180° to extract the voxel response in each direction. Combined with the gray-level similarity of the 3×3×3 voxel neighborhood, a crack characterization tensor containing spatial gradient, anisotropy, orientation sensitivity, and similarity channels is formed.

[0028] The crack representation tensor is input into a cross-layer attention deep network. Higher-level branches extract global crack morphology covering a 50-meter range, while lower-level branches extract local crack texture features using three convolutional layers combined with pyramid pooling. These high- and low-level features are then fused through a bidirectional cross-layer attention matrix to output a crack probability map and orientation tensor. By setting a crack probability threshold of 0.6, crack voxel nodes are extracted. A crack voxel map is constructed based on orientation consistency, mean crack probability, and Euclidean distance between nodes. After updating node features using a three-layer graph neural network, the master crack path is searched and generated using the weighted sum of orientation consistency and node probabilities as the cost function.

[0029] To enhance the model's ability to distinguish pseudo-cracks, various perturbations were applied to the main crack path to generate counterfactual paths. These included misaligning the node sequence order by 20%, randomly shifting the direction within ±10°, pruning the path length by 10%, and smoothing the inter-node interpolation. Paths with a Hausdorff distance and morphological similarity to the main path greater than 0.7 were used as valid counterfactual samples input to the contrastive learning model. Triple samples were constructed during the contrastive learning process, and a triple loss function with a boundary margin of 0.5 was used during training. The model converged after 200,000 iterations.

[0030] Actual test results show that, verified by well logging data from 8 drilling wells, the model of this invention has an average error of only 6.8° in identifying fracture orientation, which is 68.4% lower than the average error (21.5°) of the traditional coherence-based method. The accuracy of fracture connectivity prediction has increased from 61% to 93% of the traditional method. The relative error between fracture density estimation and core statistics has decreased from 35% of the traditional method to 8%. Furthermore, for noise areas that are easily misidentified as fractures in manual interpretation, the fracture probability value output by this method is less than 0.2, effectively eliminating the influence of noise. In terms of processing efficiency, the entire process of processing a 50-square-kilometer 3D seismic body in this block took 17 hours, while manual interpretation of fractures in a block of this size took more than a week and showed poor consistency. This verifies the automated, high-precision, and quantitative modeling capabilities of this invention in complex fracture zones.

[0031] The experiment demonstrates that the method of this invention effectively solves existing problems such as poor continuity in identifying complex-oriented fractures, difficulty in quantifying connectivity, and susceptibility to false detection due to noise. It can not only output a structured fracture model that includes orientation, length, connectivity, and density, but also significantly reduce reliance on manual intervention and interpretation time, providing a reliable technical means for fine reservoir description and development design.

[0032] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for reservoir fracture identification based on deep learning, characterized in that, Includes the following steps: S1. Receive seismic data or downhole imaging data, perform adaptive predictive filtering on the received data to suppress strip noise, remove isolated noise points through morphological filtering, perform time-depth conversion based on velocity model on the filtered data to convert time domain data into depth domain data, and use symmetric orthogonal wavelet basis to perform three-level multi-scale wavelet denoising processing on the depth domain data to generate multi-scale feature bodies. S2. Construct a crack representation tensor based on multi-scale feature volumes, including calculating the gradients of each voxel in the three spatial directions to form a spatial gradient tensor, using an elliptic anisotropy operator to calculate the anisotropic properties of voxels to generate an anisotropic tensor, using multiple directional orientation-sensitive convolution kernels covering the range of 0° to 180° to calculate the response values ​​of voxels in each direction to generate orientation-sensitive tensors, calculating the local similarity value of each voxel based on the gray-level similarity of the voxel neighborhood to generate a similarity tensor, and concatenating the tensors in the channel dimension to obtain the crack representation tensor. S3. Input the crack representation tensor into the cross-layer attention depth network. The cross-layer attention depth network includes a high-level crack morphology perception branch and a low-level crack texture perception branch. The high-level and low-level features are fused through bidirectional cross-layer attention matrix interaction to output the crack probability map and crack orientation tensor. S4. Construct a crack path planning graph based on the crack probability graph and crack orientation tensor. Use voxels with crack probability values ​​higher than a preset threshold as graph nodes. Calculate edge weights based on the orientation consistency between nodes, the mean crack probability, and the spatial geometric distance. Construct edges between nodes. Use a graph neural network to update node features multiple times. Use the shortest path search algorithm based on edge weights to generate the main crack control path. S5. Perform structural disturbance operations on the main crack control path to generate a set of counterfactual crack paths, including perturbation of the order of path nodes, perturbation of the direction of path nodes, perturbation of the path length, perturbation of the path nodes by interpolation, and filtering out counterfactual paths with morphological similarity higher than a preset threshold. S6. Input the main control path and the counterfactual crack path into the comparison embedding encoder, and embed the path node sequence, crack direction tensor and path connectivity vector into the encoder respectively. Concatenate them to form a path-level embedding vector. Maximize the embedding distance between the main control path and the counterfactual path through comparison learning, and output a structured crack distribution model containing crack direction, length, connectivity and crack density information.

2. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S1 specifically includes: S11. Adaptive prediction filtering is performed on the received seismic data or downhole imaging data. The input data is divided into blocks according to the gather direction. The filtering parameters are adaptively adjusted based on the average frequency drift within each block. A prediction error filter is applied to the data of each block to generate preliminary denoised data. S12. Perform morphological filtering on the data generated in S11. Use a 3×3 structuring element to perform grayscale expansion and erosion operations respectively to eliminate isolated high-frequency noise and smooth small-scale random noise. S13. Based on the regional velocity model, perform time-depth conversion on the data generated in S12, and convert the time-domain samples of each profile... Convert to depth domain samples The conversion formula is: in, Indicates time Longitudinal wave velocity at that location For depth location; S14. Perform three-level multi-scale wavelet decomposition on the depth domain data generated in S13, using a symmetric orthogonal wavelet basis. Decomposed into low-frequency approximate components and high-frequency detail components The number of decomposition layers The multi-scale wavelet decomposition formula is: in, and These are the low-pass and high-pass filter coefficients of the wavelet basis, respectively; S15. Based on the decomposition results of S14, a threshold denoising method is used to denoise each high-frequency component. Soft thresholding is performed, and the denoised low-frequency feature components are reconstructed. ; S16. The reconstruction result of S15 is weighted and fused with the original depth domain data to form a multi-scale feature volume containing multi-resolution information, which serves as the input for the subsequent crack characterization tensor construction.

3. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S2 specifically includes: S21. Calculate the gradient value of each voxel in the three-dimensional space direction for the multi-scale feature volume generated in S1, and obtain the gradient components in the X, Y, and Z directions using the Sobel operator. Calculate the spatial gradient magnitude The calculation formula is as follows: in, These represent the gradient components of the voxel in the X, Y, and Z directions, respectively. Indicates the magnitude of the spatial gradient; S22. Based on the gradient information from S21, the anisotropy of voxels is calculated using the elliptic anisotropy operator. The Based on the ratio of the major and minor radii of the principal axes of voxels The definition and calculation formula are as follows: in, , These are the largest and second largest eigenvalues ​​of the voxel gray-level covariance matrix; S23. Apply orientation-sensitive convolution kernels covering the range of 0° to 180° to the multi-scale feature volume generated in S1, with K directions. Calculate the voxel response value for each direction. This forms a direction-sensitive tensor, where each direction is 15° when K=12; S24. The mean square error of gray level in the 3×3×3 neighborhood of each voxel in the multi-scale feature volume generated based on S1. Calculate local similarity of voxels And normalize it. The formula is defined as follows: in, Represents the similarity value. The gray-level variance within the neighborhood; S25. The spatial gradient tensor generated in S21 The anisotropic tensor generated by S22 S23 generated direction-sensitive tensor Similarity tensor generated by S24 The cracks are spliced ​​along the channel dimension to form a crack representation tensor containing multiple crack features, which serves as the input for subsequent S3 processing.

4. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S3 specifically includes: S31. Input the crack representation tensor generated in S2 into the high-level crack morphology perception branch of the cross-layer attention deep network, and extract the global crack morphology features through five three-dimensional convolutional layers in sequence. The kernel size of each convolutional layer is 3×3×3, the stride is 1, and the receptive field covers a crack range of 50m. S32. Simultaneously input the crack representation tensor generated in S2 into the low-level crack texture perception branch of the cross-layer attention deep network. Extract local crack texture details through three 3D convolutional layers and multi-scale pyramid pooling. Pyramid pooling generates feature tensors at different scales. ; S33, the feature tensors output from the higher layers of S31 respectively and the feature tensor output by the lower layer of S32 Calculate the cross-layer bidirectional attention matrix Note the matrix The calculation formula is: in, Represents the high-level feature tensor. Represents the low-level feature tensor. For learnable weight matrix, Represents a bidirectional cross-layer attention matrix; S34, will Acting on and The system uses matrix multiplication to perform interactive updates of high-level and low-level features, and outputs the updated high-level features. and low-level features ; S35, will and The features are concatenated along the channel dimension and integrated using a 3D convolutional fusion network to output a crack probability map. and crack orientation tensor ,in This represents the probability value of each voxel belonging to a crack. This indicates the angle of the crack direction for each voxel.

5. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S4 specifically includes: S41. Crack probability map based on S3 output and crack orientation tensor Select crack probability value Greater than the set threshold voxels as candidate nodes for cracks ,in The probability threshold for crack voxels; S42, For each pair of adjacent nodes Based on the consistency of the direction between nodes Mean probability of nodal cracks Euclidean distance between nodes Calculate edge weights The formula for calculating the edge weight is: in, This indicates that nodes Nᵢ and Nⱼ follow the same direction. This represents the mean probability of cracks at two nodes. For Euclidean distance, , , These are weighting coefficients; S43, based on and Constructing a crack voxel map ,in For a set of nodes, Let it be the set of edges; S44, in the crack voxel diagram above The layered graph neural network updates the node state, and the feature update formula for each layer is as follows: in, Represents a node In the Layer feature representation, Pay attention to weights between nodes. It is a learnable transformation matrix; S45, in the updated diagram Based on the shortest path search algorithm, using a weighted sum of path consistency and node probability as the cost function, the algorithm searches for the optimal master crack path connecting the start and end nodes, and outputs the master crack path sequence. .

6. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S5 specifically includes: S51, The main control path of the crack generated in S4 Perform node order misalignment perturbation, proportionally adjusting the path node sequence. Random misalignment, This is a misaligned proportional parameter; S52, to Perform a path offset perturbation to change the path of each node in the path. By maximum offset angle exist Random disturbances within the range This is the offset angle threshold; S53, to Perform path length truncation perturbation, trimming the path proportionally from either the beginning or end. The node, This is the truncation ratio parameter; S54, to Perform node interpolation perturbation in the node sequence with average node spacing. Insert a new node to make the smoothness of the disturbance path approximate that of the master path, where for Average Euclidean distance between mid-nodes; S55. Calculate each disturbance path and morphological similarity The formula for morphological similarity is: in, for and The distance to Hausdorff; S56, will Greater than the set threshold The perturbation paths are preserved, forming a set of counterfactual crack paths with reasonable morphology. This is used for subsequent comparative learning processes.

7. The method for reservoir fracture identification based on deep learning according to claim 1, characterized in that, S6 specifically includes: S61, Convert the crack master path generated in S4 The set of counterfactual paths generated by S5 Each counterfactual path in Input the path feature encoder separately to the path node sequence Crack orientation sequence and path connectivity vector Encode the nodes to generate node embedding vectors. Towards Embedded Vectors and connectivity embedding vectors ; S62, will , , Path-level embedding vectors are formed by concatenating the data along the channel dimension. ,in S63. Employing a contrastive learning network to... Training is performed to construct triplet samples. ,in Embedded in the master control path Positive samples for the main control path. For negative samples of the counterfactual path, a triplet loss function is used during training, and the loss formula is: in This represents the Euclidean distance between embedding vectors. This refers to the boundary margin hyperparameter; S64. Based on the trained contrastive learning network, the master path... and The model is input via a mid-path, and the output includes the crack orientation. Crack length Crack connectivity and crack density Structured crack distribution model It is used for the quantitative interpretation and visualization of crack systems.

Images