A method and system for extracting fracture rods based on morphological erosion algorithm

Abstract

This invention relates to the field of petroleum geological exploration and geophysical data processing technology, and discloses a method and system for fault rod extraction based on morphological corrosion algorithm. By removing noise from low signal-to-noise ratio data and accurately segmenting fault fracture zones, corrosion calculations are performed based on directional anisotropic structural elements customized according to the tectonic stress field, ensuring the corrosion operation follows the dominant direction of the fault. Iterative corrosion and feedback correction methods constrained by the fault trend model are used to compress the fracture zone towards the fault core axis. Multi-attribute fusion verification and topology optimization methods are employed to remove prostheses on fault rods and fill in missing parts, giving the fault rods greater realism and spatial continuity. The generated data volume is input into 3D kriging interpolation for gap filling, and then transformed by a coordinate system, so that the output fault rod data volume directly meets the accuracy requirements for structural framework modeling. This invention proposes a fault rod extraction method suitable for low signal-to-noise ratio data.

Description

Technical Field

[0001] This invention relates to the field of petroleum geological exploration and geophysical data processing technology, and in particular to a method and system for extracting fault rods based on a morphological corrosion algorithm. Background Technology

[0002] In the field of seismic exploration data processing, coherence volume technology and ant volume technology are mainly used to detect faults. These two technologies use the continuous changes in the seismic reflection phase axis to represent the distribution characteristics of underground faults. However, due to the complex surface conditions, stratum absorption attenuation, and interference from acquisition noise, the signal-to-noise ratio of seismic data is low. Therefore, seismic data with low signal-to-noise ratio cannot accurately distinguish the fault trunk and the fracture zone around the fault. As a result, the output attribute volume is represented as a wide and continuous fault fracture zone, rather than the narrow strip fault response required.

[0003] Existing fault rod extraction techniques are mostly designed based on high signal-to-noise ratio data and generally employ a scheme of single threshold segmentation combined with simple morphological operations. These methods have two key drawbacks when dealing with wide fault fracture zones: first, they struggle to isolate the true fault core axis from the fracture zone, resulting in extraction results that either excessively include fracture zone noise or lose fault extension information due to excessively high thresholds; second, they lack constraints on geometric features such as fault strike and dip, leading to poor continuity of extracted fault rods and numerous spurious breaks, failing to directly meet the accuracy requirements for fault spatial morphology in lattice modeling.

[0004] To address the aforementioned technical challenges, there is an urgent need to develop a method for extracting tomographic bars from low signal-to-noise ratio data. Summary of the Invention

[0005] This invention provides a method and system for extracting fracture rods based on a morphological erosion algorithm, aiming to solve at least one of the above-mentioned technical problems.

[0006] To achieve the above objectives, this invention provides a method for extracting fracture rods based on a morphological erosion algorithm, the method comprising the following steps:

[0007] The original ant body is subjected to noise suppression, smoothing and enhancement optimization processing to generate an optimized ant body. The optimized ant body is then subjected to threshold segmentation to obtain a binary mask body.

[0008] Determine the fault axis and design different types of morphological structural elements; wherein, the morphological structural elements include principal direction structural elements, secondary direction structural elements and isotropic spherical structural elements;

[0009] The first round of erosion operation is performed on the binarized mask using the main direction structural element to generate a preliminary fault skeleton. With the fault trend model as a constraint, the auxiliary direction structural element and the isotropic spherical structural element are used to perform iterative erosion and feedback correction on the preliminary fault skeleton to generate an optimized fault skeleton.

[0010] The optimized fault skeleton is subjected to three-dimensional gradient enhancement and refinement to form the initial fault rod;

[0011] Multi-attribute fusion verification is performed on the initial fault bar to generate verified fault bars. Three-dimensional connectivity analysis and topology optimization, as well as spatial interpolation processing for lattice adaptation, are then performed on the verified fault bars to generate the final fault plane.

[0012] Optionally, the original ant body undergoes noise suppression, smoothing, and enhancement optimization processes to generate an optimized ant body, specifically including:

[0013] Obtain the global mean and global standard deviation of the original ant body. Based on the set outlier judgment threshold, traverse each voxel in the original ant body. When a voxel meets the judgment requirements, replace the voxel value with the median of all valid voxels in the preset area around the voxel to remove outliers.

[0014] For the ant body after outlier removal, the filter window size is determined according to the target layer burial depth, and the attribute body is smoothed point by point using the Gaussian filter formula.

[0015] Perform grayscale histogram statistics on the smoothed ant body and calculate the voxel cumulative distribution function for each grayscale level. Based on the histogram equalization operation, stretch the dynamic range of grayscale of the data volume to the standard range to output the optimized ant body.

[0016] Optionally, threshold segmentation is performed on the optimized ant body to obtain a binary mask, specifically including:

[0017] Based on the optimized gray-level histogram of ant bodies, the global initial segmentation threshold is calculated using the maximum inter-class variance method.

[0018] The optimized ant body is divided into multiple 10×10 local processing units. The local mean and local standard deviation of each local processing unit are calculated. Combined with the global initial segmentation threshold, the local segmentation threshold of the local processing unit is generated.

[0019] In the optimized ant body, each voxel is traversed and compared with the local segmentation threshold of its local unit. When the voxel value is greater than or equal to the local segmentation threshold, it is set to 1; otherwise, it is set to 0, so as to obtain the binary mask body.

[0020] Optionally, the fault axis is determined and different types of morphological structural elements are designed, including:

[0021] Based on the rose diagram statistical fault strike used to characterize the tectonic features of the study area, the dominant strike angle of the main fault in the study area is determined. Based on the dominant strike angle, different types of morphological structural elements are designed; wherein, the morphological structural elements include:

[0022] Main structural element: A long strip structure is adopted, with a length equal to 5-7 seismic sampling points and a width equal to 1 seismic sampling point. The extension direction of the structural element is oriented at the angle to the strike of the main fault. completely consistent;

[0023] Auxiliary directional structural elements: short strip structures are adopted, with a length of 3 seismic sampling points and a width of 1 seismic sampling point. The extension direction of the structural elements is oriented at the angle to the strike of the main fault. vertical;

[0024] Isotropic spherical structural element: Utilizing a spherical structure with a radius of 1 seismic sampling point.

[0025] Optionally, the first round of erosion operation is performed on the binarized mask using the main direction structuring element to generate a preliminary fault skeleton, specifically including:

[0026] The morphological first-round erosion operation is performed on the binarized mask using the main direction structuring element to obtain the first-round erosion result;

[0027] In the first round of corrosion calculation, multiple corrosion calculations are performed according to the set number of first round corrosion iterations;

[0028] In the first round of corrosion calculation, after each corrosion calculation, the area threshold method is used to determine the connected regions with fewer than 3 voxels as noise patches, and the voxel value corresponding to the noise patches is assigned to 0 in order to generate a preliminary fault skeleton.

[0029] Optionally, using the fault trend model as a constraint, and utilizing the auxiliary directional structural elements and the isotropic spherical structural elements, iterative erosion and feedback correction are performed on the preliminary fault skeleton to generate an optimized fault skeleton, specifically including:

[0030] A three-dimensional Hough transform was performed on the preliminary fault skeleton, and the peak value was determined in the three-dimensional parameter space to establish a regional fault trend model.

[0031] Using the fault trend model as a constraint, the auxiliary directional structural elements and spherical structural elements are alternately used to perform iterative corrosion calculations on the fault skeleton prototype after the first round of corrosion, so as to generate an optimized fault skeleton.

[0032] In the iterative erosion operation, the stopping criterion for iterative erosion is set as the rate of change of the number of skeleton voxels between two adjacent iterations being less than 5%;

[0033] In the iterative erosion operation, after each iteration of the erosion operation, the degree of fit between the current skeleton body and the fault trend model is calculated. When the degree of fit is less than 0.7, it is determined that the skeleton body in the region deviates from the fault trend, and morphological expansion correction of spherical structural elements is performed on the region.

[0034] Optionally, the optimized fault skeleton is subjected to three-dimensional gradient enhancement and refinement to form an initial fault rod, specifically including:

[0035] The gradient of the optimized fault skeleton is calculated using the three-dimensional Sobel gradient operator, and the solutions are obtained respectively. gradient components in three directions The gradient magnitude volume is obtained through vector synthesis;

[0036] The gradient magnitude volume is linearly normalized, and a parallel thinning algorithm is used to thin the normalized gradient magnitude volume to generate an initial fault bar with rod-like characteristics.

[0037] Optionally, multi-attribute fusion verification is performed on the initial fault rod to generate a verified fault rod, specifically including:

[0038] Three types of auxiliary attribute volumes—post-stack seismic amplitude volume, root mean square amplitude volume, and curvature volume—were selected for optimization processing including noise suppression, smoothing, and enhancement.

[0039] Spatial registration is performed between the initial fault rod and the three types of auxiliary attribute volumes to ensure consistency of three-dimensional coordinates, and the correlation coefficient is calculated for each voxel in the initial fault rod.

[0040] Based on the comparison results between the correlation coefficient and the preset correlation threshold, the fault segments are eliminated and completed to generate verified fault bars.

[0041] Optionally, the verified fault bars are subjected to three-dimensional connectivity analysis and topology optimization, and spatial interpolation processing for lattice adaptation is performed sequentially to generate the final fault plane, specifically including:

[0042] The connectivity of the verified fault rods is tracked using a three-dimensional seed filling algorithm. The spatial range of each connected branch is recorded using a 6-neighbor connectivity tracking method to obtain the geometric parameters of each connected branch. The geometric parameters include the length, orientation angle, and dip angle of each connected branch.

[0043] Short, invalid branches with fewer than 10 seismic sampling points are removed based on the connected component length threshold to obtain a set of candidate connected components;

[0044] The candidate connected branch set is subjected to topology simplification and optimization processing. At the intersection of fault bars, branches that match the fault trend model are retained. At the bifurcation point, small bifurcations with a width of less than 1 voxel are removed. At the twist point, the fault bar direction is corrected by smooth interpolation to generate topology-optimized fault bars.

[0045] Three-dimensional kriging interpolation is used to characterize the spatial correlation of voxel values ​​through a variogram function, so as to perform interpolation gap filling on the topology-optimized fault rods.

[0046] After completing the interpolation and gap filling, the grid parameter step size and origin coordinates of the structural lattice model are obtained, and the seismic sampling coordinates are converted into the node coordinates of the structural lattice model to generate the final fault plane.

[0047] Furthermore, to achieve the above objectives, the present invention also provides a tomographic rod extraction system based on a morphological erosion algorithm, comprising:

[0048] The segmentation module is used to perform noise suppression, smoothing and enhancement optimization processing on the original ant body to generate an optimized ant body. The optimized ant body is then segmented by threshold to obtain a binary mask body.

[0049] A determination module is used to determine the fault axis and design different types of morphological structural elements; wherein, the morphological structural elements include principal direction structural elements, secondary direction structural elements and isotropic spherical structural elements;

[0050] The generation module is used to perform the first round of erosion operation on the binary mask using the main direction structural element to generate a fault skeleton prototype. Using the fault trend model as a constraint, the auxiliary direction structural element and the isotropic spherical structural element are used to perform iterative erosion and feedback correction on the fault skeleton prototype to generate an optimized fault skeleton.

[0051] The module is used to perform three-dimensional gradient enhancement and refinement on the optimized fault skeleton to form the initial fault rod.

[0052] The verification module is used to perform multi-attribute fusion verification on the initial fault rod, generate the verified fault rod, and then perform three-dimensional connectivity analysis and topology optimization, as well as spatial interpolation processing for lattice adaptation on the verified fault rod to generate the final fault plane.

[0053] The beneficial effects of this invention are as follows:

[0054] (1) Solving the technical pain points under low signal-to-noise ratio: The combined noise reduction optimization method of multi-scale Gaussian filtering, standard deviation threshold denoising and histogram equalization effectively removes noise in low signal-to-noise ratio data; the method of combining global OTSU threshold segmentation and local dynamic correction segmentation is used to accurately segment fault fracture zones; the problem of traditional methods being difficult to apply effectively under low signal-to-noise ratio conditions is solved.

[0055] (2) Improve the accuracy and continuity of fault rod extraction: perform corrosion calculation based on the directional anisotropic structural elements customized by the tectonic stress field, so that the corrosion calculation is carried out along the dominant direction of the fault; use iterative corrosion and feedback correction methods constrained by the fault trend model to compress the fracture zone towards the core axis of the fault; use multi-attribute fusion verification and topology optimization methods to remove the prostheses on the fault rods and complete the missing parts on the fault rods, so that the fault rods have stronger realism and spatial continuity.

[0056] (3) Seamless integration of lattice modeling and construction: The generated data volume is input into the three-dimensional kriging interpolation and then transformed by the coordinate system. The final output fault rod data volume can directly meet the accuracy required for lattice modeling. No further format conversion or data modification is required, which greatly improves the work efficiency after the construction interpretation. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating the overall technical process of the present invention.

[0058] Figure 2 a represents the original ant body cross-section in this invention. Figure 2 b is the optimized and noise-suppressed ant body cross-section in this invention;

[0059] Figure 3 This is a cross-section of the binarized mask body in this invention;

[0060] Figure 4 This is the dominant strike rose diagram of the main fault in the study area of ​​this invention;

[0061] Figure 5 This is a cross-sectional view of the preliminary fault skeleton using the main directional structural elements to perform the first round of corrosion calculation in this invention.

[0062] Figure 6 This is the optimized fault skeleton cross-section display effect in this invention;

[0063] Figure 7 This is a 3D display effect diagram of the initial fracture bar after 3D gradient enhancement and parallel refinement processing in this invention;

[0064] Figure 8This is a three-dimensional display effect diagram of the topology-optimized fault rod in this invention;

[0065] Figure 9 A three-dimensional display effect diagram of the fault bar of the spatial interpolation surface adapted to the structural lattice in this invention;

[0066] Figure 10 This is a schematic diagram of the system structure of the present invention. Detailed Implementation

[0067] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0068] It should be noted that existing fault rod extraction techniques are mostly designed based on high signal-to-noise ratio data and generally adopt a scheme of single threshold segmentation combined with simple morphological operations. These methods have two key drawbacks when dealing with wide fault fracture zones: first, it is difficult to isolate the true fault core axis from the fracture zone, resulting in extraction results that either excessively include fracture zone noise or lose fault extension information due to excessively high thresholds; second, they lack constraints on geometric features such as fault strike and dip, resulting in poor continuity of extracted fault rods and numerous false breaks, failing to directly meet the accuracy requirements for fault spatial morphology in lattice modeling.

[0069] To address the aforementioned issues, this embodiment proposes a fault rod extraction method and system based on a morphological erosion algorithm. By removing noise from low signal-to-noise ratio (SNR) data and accurately segmenting fault fracture zones, erosion calculations are performed using directional anisotropic structural elements customized based on the tectonic stress field, ensuring the erosion operation follows the dominant fault direction. Iterative erosion and feedback correction methods constrained by a fault trend model are used to compress the fracture zone towards the fault core axis. Multi-attribute fusion verification and topology optimization methods are employed to remove prostheses from fault rods and fill in missing parts, enhancing the realism and spatial continuity of the fault rods. The generated data volume is input into a 3D kriging interpolation system for gap filling, and then transformed using a coordinate system, ensuring the output fault rod data volume directly meets the accuracy requirements for lattice modeling. This presents a fault rod extraction method suitable for low SNR data.

[0070] Specifically, such as Figure 1 As shown in this embodiment, a method for extracting fault rods based on a morphological erosion algorithm is proposed, which includes the following steps:

[0071] S1: Perform optimization processing on the ant body, such as noise suppression, smoothing, and enhancement.

[0072] In this embodiment of the invention, random noise and abnormal interference in the original ant body are suppressed by multiple means to enhance the property differences between the fault fracture zone and the surrounding rock, providing a high-quality data foundation for subsequent segmentation. The specific operations are as follows:

[0073] Outlier removal using standard deviation thresholding: Obtaining the initial ant body global mean Compared with global standard deviation Set the outlier threshold as follows: Iterate through each voxel in the data volume. If the voxel value satisfies... If a value is found to be an isolated noise point, its value is replaced with the median of all valid voxels within a 3×3×3 neighborhood around that voxel to avoid interference from a single extreme value in subsequent processing.

[0074] Multi-scale Gaussian filtering smoothing: For the ant-like data after outlier removal, the filter window size is determined based on the target layer depth—a 3×3×3 three-dimensional cubic filter window is used for shallow strata with a depth less than 500m; a 5×5×5 three-dimensional cubic filter window is used for deep strata with a depth greater than or equal to 500m. The attribute volume is then smoothed point-by-point using the Gaussian filtering formula:

[0075] ;

[0076] Histogram equalization attribute enhancement: Gray-level histogram statistics are performed on the smoothed ant volume, and the cumulative distribution function of voxels for each gray level is calculated. Histogram equalization stretches the gray-level dynamic range of the data volume to the standard interval [0, 1], so that the high attribute value area of ​​the fault fracture zone contrasts sharply with the low attribute value area of ​​the surrounding rock. The final output is the optimized ant volume. . Figure 2 a is a cross-section of the original ant body. Figure 2 b represents the optimized ant body cross-section;

[0077] S2: Threshold segmentation is performed on the optimized ant body to obtain a binary mask body.

[0078] In this embodiment of the invention, a precise segmentation of fault fracture zones and non-fault regions is achieved by combining global thresholding and local correction, thus solving the problem of threshold ambiguity under low signal-to-noise ratio. The specific operation is as follows:

[0079] OTSU method global threshold calculation: based on optimized ant body The grayscale histogram was used to calculate the global initial segmentation threshold using the Otsu's method (OTSU). The core principle of this method is to divide the grayscale values ​​of the attribute volume into a fault-bounded zone class (foreground) and a non-fault-bounded class (background). The optimal threshold is determined by maximizing the inter-class variance between the two classes. The formula for calculating the inter-class variance is:

[0080] ;

[0081] In the formula: The proportion of fault fracture zone-like voxels to the total number of voxels in the data volume; The proportion of non-tomographic voxels to the total number of voxels in the data volume, and satisfying the following conditions: ; The average gray value of the voxel-like structure in the fault fracture zone; The average gray value of the non-tomographic voxel; This represents the average grayscale value of all voxels in the data volume.

[0082] Dynamic correction of local window threshold: This optimizes the ant body... According to the seismic gather partitioning rules, it was divided into multiple 10×10 gather local processing units. For each local processing unit... Calculate its local mean With local standard deviation The local segmentation threshold of the unit is adjusted using a correction formula. The formula is:

[0083] ;

[0084] In the formula: This is the mean correction factor, with a value ranging from 0.1 to 0.3, used to correct local grayscale shifts. This is the standard deviation correction factor, with a value range of 0.05-0.15, used to match the degree of grayscale dispersion in local areas.

[0085] Binarization mask generation: in the attribute body For each voxel in the array, iterate through the voxels and match the voxels located in a local unit with the local threshold of that local unit. By comparison, if the voxel value is greater than or equal to the local threshold, the voxel is considered to belong to the fault fracture zone and its value is set to 1; otherwise, the voxel belongs to the non-fault region and its value is set to 0. Finally, the binarized mask volume is obtained. . Figure 3 This is a cross-section of the binary mask.

[0086] S3: Determine the fault orientation and design different types of morphological structural elements.

[0087] In this embodiment of the invention, morphological structural elements with anisotropic orientation are customized according to the structural features of the study area, so that the extension direction of the structural elements is highly matched with the fault strike, thereby improving the targeting of subsequent erosion operations. The specific operation is as follows:

[0088] Fault strike determination in the study area: Data such as structural geological survey reports and seismic profile interpretation results were collected. Fault strikes were statistically analyzed using rose diagrams to determine the dominant strike angle of the main faults in the study area. (With true north as 0°, clockwise rotation is the positive direction), from Figure 4 The dominant strike rose diagram of the main fault in the statistical study area shows that the fault strike ranges from approximately 25° to 65°.

[0089] Design of multiple structural elements: based on the strike angle of the main fault Three different types of morphological structural elements were designed to adapt to the corrosion requirements at different stages:

[0090] (1) Main direction structural elements It adopts a long strip structure, with a length of Take 5-7 seismic sampling points, width Take one seismic sampling point, and determine the angle between the extension direction of the structural element and the strike of the main fault. Completely consistent, used for precise contraction along the fault strike;

[0091] (2) Auxiliary directional structural elements It adopts a short strip structure, with a length of Take 3 seismic sampling points, width Take one seismic sampling point, and determine the angle between the extension direction of the structural element and the strike of the main fault. Vertical, used for fine trimming perpendicular to the fault strike;

[0092] (3) Isotropic spherical structural elements: Using a spherical structure with a radius of 1 seismic sampling point, the fault edges are smoothed to avoid sharp edges.

[0093] The orientation matching verification of the structural elements was carried out. The orientation matching verification of the three designed structural elements and the binary mask was performed. A typical fault profile was selected for trial calculation verification. The results showed that the extension direction of the structural elements was consistent with the direction of the fault fracture zone. The corrosion calculation direction determined by the result will be used for subsequent calculations.

[0094] S4: Based on the binary mask, the first round of erosion operation is carried out using the main direction structuring element to obtain the preliminary fault skeleton after the first round of erosion.

[0095] In this embodiment of the invention, the main direction structuring element is used to structurate the binary mask. The corrosion treatment is performed to achieve initial contraction of the fault fracture zone towards the fault core region. The specific operation is as follows:

[0096] Morphological erosion operation execution: using the main direction structuring element Binarized fault fracture zone mask Perform morphological erosion operation. The mathematical expression for the erosion operation is:

[0097] ;

[0098] In the formula: For morphological erosion operators; structural element Translate to position The set after; the physical meaning of this formula is: only when the structural element Translate to position At that time, all its voxel points are completely contained within the mask volume. Within the foreground region (value 1), the position Only then will it be preserved in the corrosion results middle.

[0099] First, the number of iterations for the first round of corrosion is set to 2 to 3, and multiple rounds of corrosion calculations are completed. After each round of corrosion calculations, connected region analysis is performed on the corrosion results, and isolated small patches are removed using the area threshold method. That is, connected regions with less than 3 voxels are identified as noise patches, and the voxel values ​​corresponding to these patches are assigned to 0 to avoid small patches interfering with the fault skeleton.

[0100] First-round corrosion results output: After completing multiple rounds of iterative corrosion and small patch removal, the preliminary fault skeleton after the first round of corrosion is output. At this point, the skeleton has initially achieved the contraction of the fracture zone and retained the core extension characteristics of the fault, but there are still some redundant areas that need further optimization. Figure 5 The effect of the cross-sectional view of the fault skeleton prototype for the first round of corrosion calculation using the main directional structural elements.

[0101] S5: Iterative erosion and feedback correction based on fault strike constraints yield the optimized fault skeleton.

[0102] In this embodiment of the invention, the first step is to build a fault trend model, which involves first removing the erosion and then going back to make feedback corrections. Based on this cyclical process, the fault skeleton is continuously refined.

[0103] The initial shape of the fault skeleton after the first round of corrosion. Perform a three-dimensional Hough transform on it in the three-dimensional parameter space (towards angle). Inclination angle ,inclination The corresponding peak value is found, and then the geometric parameters such as the strike angle, dip angle, and dip angle of each branch are obtained. Based on the above parameters and the tectonic background of the corresponding study area, a regional fault trend model is established. This clarifies the dominant distribution direction and dip range of faults in the study area.

[0104] Alternating Iterative Corrosion of Multiple Structural Elements: Using a Fault Trend Model To constrain the use of auxiliary directional structural elements, alternately use them. With spherical structural elements The initial shape of the fault skeleton after the first round of corrosion Perform iterative erosion calculations; the formula for iterative erosion is:

[0105] ;

[0106] In the formula: Let be the number of iterations, and be the initial value. ; For the first The skeleton after the next iteration; For the first The skeleton after the next iteration.

[0107] Setting and executing the iterative stopping criterion: The stopping criterion for iterative erosion is set as the rate of change of the number of skeleton voxels between two adjacent iterations. For values ​​less than 5%, the formula for calculating the rate of change in voxel count is:

[0108] ;

[0109] In the formula: For the first The total number of voxels in the skeleton after the next iteration; For the first The total number of voxels in the skeleton after each iteration. The rate of change is calculated sequentially during the iteration process. ,when At that point, it was determined that the morphology of the skeleton had become stable, and the iterative corrosion was stopped.

[0110] Trend fit feedback correction: After each iteration of erosion, the current skeleton body is calculated. With fault trend model The degree of matching The fit calculation uses a cosine similarity algorithm, which quantifies the degree of matching between the skeleton and the trend model by comparing the angle between their direction vectors. The system was determined to deviate from the fault trend in the region, and a morphological dilation correction was performed on the region. The dilation calculation used spherical structural elements. The expansion formula is:

[0111]

[0112] In the formula: This is a morphological dilation operator; through dilation correction, it appropriately expands the skeleton that deviates from the trend, bringing it back to the dominant fault distribution direction. After all iterations and corrections are completed, the optimized fault skeleton is output. . Figure 6 This is the optimized display effect of the fault skeleton profile.

[0113] S6: Perform three-dimensional gradient enhancement and refinement on the optimized fault skeleton to form the initial fault rod.

[0114] In this embodiment of the invention, the core of this step is to highlight the edge features of the fault skeleton through gradient enhancement, and then shrink the skeleton body to the width of a single element through a parallel thinning algorithm to form an initial fault rod. The specific operation is as follows:

[0115] 3D Sobel Gradient Operator Edge Enhancement: The 3D Sobel gradient operator is used to enhance the optimized fault skeleton. Perform gradient calculations and solve them separately. gradient components in three directions Then, the gradient magnitude volume is obtained through vector synthesis. The calculation formula is:

[0116] ;

[0117] Gradient magnitude volume It can accurately delineate the boundaries of the fault skeleton, with the gradient amplitude at the edges being much greater than that inside, which is beneficial for highlighting the fault edges.

[0118] Gradient magnitude volume normalization processing, for gradient magnitude volume Linear normalization transforms the range of gradient magnitude values ​​to the [0, 1] interval, thereby eliminating the adverse effects of excessive differences in gradient magnitude between regions, which is beneficial for subsequent thresholding and refinement operations.

[0119] The normalized gradient-enhanced volume skeleton is refined with the aid of the Zhang-Suen parallel refinement algorithm. The refinement process is completed in two parallel iterations, targeting voxels that meet the criteria for "deletable points" (based on connectivity and the distribution of surrounding voxels). This ensures that the fracture skeleton remains coherent and complete even as voxels are continuously deleted. After an infinite number of refinement iterations, the fracture skeleton degenerates into the width of a single voxel, and the final output value is an initial fracture bar with clearly defined rod-like features. . Figure 7This is a 3D rendering of the initial fault bars after 3D gradient enhancement and parallel refinement processing.

[0120] S7: Verify the authenticity of the initial fault rod by using multi-attribute fusion to obtain the verified fault rod.

[0121] In this embodiment of the invention, the core of this step is to introduce multiple types of auxiliary seismic attributes, verify the authenticity of the initial fault rod through correlation analysis, eliminate false fault segments, and supplement missing fault segments. The specific operations are as follows:

[0122] Auxiliary attribute volume selection and preprocessing: Selecting post-stack seismic amplitude volume Root mean square amplitude body curvature body Three types of auxiliary attribute volumes are closely related to fault development. Each type of auxiliary attribute volume undergoes the same noise reduction and optimization process as in step 1 to ensure that the quality of the auxiliary attribute volumes is consistent with that of the target attribute volume.

[0123] Correlation coefficient calculation and threshold determination: Initial fault rods Spatial registration was performed with three types of auxiliary attribute volumes to ensure that the three-dimensional coordinates of all attribute volumes were completely consistent. Each voxel in the initial fault rod was traversed, and the correlation coefficient between the voxel value and the corresponding auxiliary attribute voxel value was calculated using the following formula:

[0124] ;

[0125] In the formula: The first in the initial fault rod The values ​​of individual elements; Let be the value of the i-th voxel in the auxiliary attribute body; The mean value of all voxels in the initial fault rod; This is the mean value of the voxels at the corresponding positions of the auxiliary attribute body. A correlation threshold is set. If the average correlation coefficient between a fault segment and the three types of auxiliary attribute bodies is less than If so, it is determined to be a suspected fault segment.

[0126] Manual interactive verification and correction: Suspicious fault segments are superimposed, and manual interactive interpretation removes false fault segments caused by noise. Furthermore, high-value anomalies in auxiliary attribute volumes are used to locate fault segments lost due to excessive corrosion. Linear interpolation is used to fill in the voxel values ​​of missing segments, finally yielding the verified fault bars. .

[0127] S8: Perform three-dimensional connectivity analysis and topology optimization on the verified fault rods to obtain the topology-optimized fault rods.

[0128] In this embodiment of the invention, the fault rod obtained in step S7 is a result without any combination relationship. In reality, a fault rod is a cross-section composed of multiple faults. Therefore, it is necessary to perform connectivity tracing and topology simplification to improve the spatial continuity and geometric regularity of the fault rod. The specific operations are as follows:

[0129] Using a 3D seed filling algorithm to identify fault rod candidates Connectivity is traced by performing a seed-filling algorithm at each isolated high-value voxel selected as a candidate, and then using the filling start point as the center point for 6-neighborhood connectivity (up, down; left, right; front, back) tracing. The spatial extent of each traced connected component is recorded. After the tracing is complete, the length of each connected component is obtained. , directional angle Inclination angle Equal geometric parameters.

[0130] Invalid branch removal: Based on the geological characteristics of fault development in the study area, the threshold for the length of connected branches was determined: short, branch-like invalid branches with fewer than 10 seismic sampling points were removed; main branches with 10 or more seismic sampling points were retained, and short, noisy branches were removed, thereby reducing the impact of short, noisy branches on fault interpretation.

[0131] Topology simplification and optimization: Topology simplification methods are used to process topologically complex regions such as fault rod intersections, bifurcations, and twists. At fault rod intersections, branches with a high degree of matching to the fault trend model are retained. At bifurcation points, small bifurcations with a width less than one voxel are removed. At twist points, smooth interpolation is used to correct the orientation of the fault rods, causing them to extend in the dominant trend direction. Finally, the topology-optimized fault rods are output. . Figure 8 This is a 3D rendering of the fault rod after topology optimization.

[0132] S9: Perform spatial interpolation on the topology-optimized fault bars to construct a lattice for adaptation, and obtain the fault plane.

[0133] In this embodiment of the invention, spatial interpolation for lattice adaptation is performed on the topology-optimized fault bars to obtain the fault plane. The specific operation is as follows:

[0134] Three-dimensional kriging interpolation for gap filling: Three-dimensional kriging interpolation is used to fill gaps in the topology-optimized fault bars. Interpolation is performed to fill in the gaps in fault bars caused by discontinuities in seismic data. The core of Kriging interpolation is to characterize the spatial correlation of voxel values ​​using a variogram, the formula for which is:

[0135] ;

[0136] In the formula: The spatial distance between two sample voxels; For distance equal to The number of sample voxel pairs; For sample voxels The value at; To match the sample voxels Distance is The voxel values ​​are determined. Through variogram fitting and interpolation calculation, seamless filling of fault rod gaps is achieved, improving the spatial continuity of fault rods.

[0137] Coordinate system transformation: Obtain mesh parameters for subsequent lattice modeling, including the mesh's coordinate system transformation. Step size in three directions and the origin coordinates of the modeling area Seismic sampling coordinates of the fault rod Convert node coordinates to lattice modeling The conversion formula is:

[0138]

[0139] By transforming coordinates, the spatial position of the fault rods is made to perfectly match the mesh system used for modeling the structural lattice.

[0140] Final Fault Plane Output: After interpolation and coordinate transformation, the output can be directly used to construct the fault plane for lattice modeling. . Figure 9 A 3D view of the cross-sectional plane after interpolation to fit the lattice structure.

[0141] Reference Figure 10 , Figure 10 This is a schematic diagram of the tomographic rod extraction system based on the morphological corrosion algorithm according to an embodiment of the present invention.

[0142] like Figure 10 As shown, the tomographic rod extraction system based on morphological erosion algorithm proposed in this embodiment of the invention includes:

[0143] The segmentation module 10 is used to perform noise suppression, smoothing and enhancement optimization processing on the original ant body to generate an optimized ant body, and to perform threshold segmentation on the optimized ant body to obtain a binary mask body.

[0144] The determination module 20 is used to determine the fault axis and design different types of morphological structural elements; wherein, the morphological structural elements include main direction structural elements, auxiliary direction structural elements and isotropic spherical structural elements;

[0145] The generation module 30 is used to perform the first round of erosion operation on the binary mask using the main direction structural element to generate a fault skeleton prototype. Using the fault trend model as a constraint, the auxiliary direction structural element and the isotropic spherical structural element are used to perform iterative erosion and feedback correction on the fault skeleton prototype to generate an optimized fault skeleton.

[0146] Module 40 is used to perform three-dimensional gradient enhancement and refinement on the optimized fault skeleton to form an initial fault rod.

[0147] The verification module 50 is used to perform multi-attribute fusion verification on the initial fault rod, generate the verified fault rod, and sequentially perform three-dimensional connectivity analysis and topology optimization, and spatial interpolation processing for lattice adaptation on the verified fault rod to generate the final fault plane.

[0148] Other embodiments or specific implementations of the tomographic rod extraction system based on the morphological corrosion algorithm of the present invention can be referred to the above-described method embodiments, and will not be repeated here.

[0149] It is understood that in the description of this specification, references to terms such as "one embodiment," "another embodiment," "other embodiments," or "first embodiment to Nth embodiment," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0150] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or system that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or system. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or system that includes that element.

[0151] The above are merely preferred embodiments of the present invention and do not limit the scope of the patent. Any equivalent structural or procedural transformations made based on the description and drawings of the present invention, or direct or indirect applications in other related technical fields, are similarly included within the scope of patent protection of the present invention.

Claims

1. A method for extracting fracture rods based on a morphological erosion algorithm, characterized in that, The method includes the following steps: The original ant body is subjected to noise suppression, smoothing and enhancement optimization processing to generate an optimized ant body. The optimized ant body is then subjected to threshold segmentation to obtain a binary mask body. Determine the fault axis and design different types of morphological structural elements; wherein, the morphological structural elements include principal direction structural elements, secondary direction structural elements and isotropic spherical structural elements; The first round of erosion operation is performed on the binarized mask using the main direction structural element to generate a preliminary fault skeleton. With the fault trend model as a constraint, the auxiliary direction structural element and the isotropic spherical structural element are used to perform iterative erosion and feedback correction on the preliminary fault skeleton to generate an optimized fault skeleton. The optimized fault skeleton is subjected to three-dimensional gradient enhancement and refinement to form the initial fault rod; Multi-attribute fusion verification is performed on the initial fault bar to generate verified fault bars. Three-dimensional connectivity analysis and topology optimization, as well as spatial interpolation processing for lattice adaptation, are then performed on the verified fault bars to generate the final fault plane.

2. The method for extracting fracture rods based on morphological erosion algorithm as described in claim 1, characterized in that, The original ant body undergoes noise suppression, smoothing, and enhancement optimization processes to generate an optimized ant body, specifically including: Obtain the global mean and global standard deviation of the original ant body. Based on the set outlier judgment threshold, traverse each voxel in the original ant body. When a voxel meets the judgment requirements, replace the voxel value with the median of all valid voxels in the preset area around the voxel to remove outliers. For the ant body after outlier removal, the filter window size is determined according to the target layer burial depth, and the attribute body is smoothed point by point using the Gaussian filter formula. Perform grayscale histogram statistics on the smoothed ant body and calculate the voxel cumulative distribution function for each grayscale level. Based on the histogram equalization operation, stretch the dynamic range of grayscale of the data volume to the standard range to output the optimized ant body.

3. The method for extracting fault rods based on morphological erosion algorithm as described in claim 1, characterized in that, The optimized ant body is subjected to threshold segmentation to obtain a binary mask, specifically including: Based on the optimized gray-level histogram of ant bodies, the global initial segmentation threshold is calculated using the maximum inter-class variance method. The optimized ant body is divided into multiple 10×10 local processing units. The local mean and local standard deviation of each local processing unit are calculated. Combined with the global initial segmentation threshold, the local segmentation threshold of the local processing unit is generated. In the optimized ant body, each voxel is traversed and compared with the local segmentation threshold of its local unit. When the voxel value is greater than or equal to the local segmentation threshold, it is set to 1; otherwise, it is set to 0, so as to obtain the binary mask body.

4. The method for extracting fault rods based on morphological erosion algorithm as described in claim 1, characterized in that, Determine the fault axis and design different types of morphological structural elements, specifically including: Based on the rose diagram statistical fault strike used to characterize the tectonic features of the study area, the dominant strike angle of the main fault in the study area is determined. Based on the dominant strike angle, different types of morphological structural elements are designed; wherein, the morphological structural elements include: Main structural element: A long strip structure is adopted, with a length equal to 5-7 seismic sampling points and a width equal to 1 seismic sampling point. The extension direction of the structural element is oriented at the angle to the strike of the main fault. completely consistent; Auxiliary directional structural elements: short strip structures are adopted, with a length of 3 seismic sampling points and a width of 1 seismic sampling point. The extension direction of the structural elements is oriented at the angle to the strike of the main fault. vertical; Isotropic spherical structural element: Utilizing a spherical structure with a radius of 1 seismic sampling point.

5. The method for extracting fault rods based on morphological erosion algorithm as described in claim 4, characterized in that, The first round of erosion operation is performed on the binarized mask using the main directional structuring element to generate a preliminary fault skeleton, specifically including: The morphological first-round erosion operation is performed on the binarized mask using the main direction structuring element to obtain the first-round erosion result; In the first round of corrosion calculation, multiple corrosion calculations are performed according to the set number of first round corrosion iterations; In the first round of corrosion calculation, after each corrosion calculation, the area threshold method is used to determine the connected regions with fewer than 3 voxels as noise patches, and the voxel value corresponding to the noise patches is assigned to 0 in order to generate a preliminary fault skeleton.

6. The method for extracting fault rods based on morphological erosion algorithm as described in claim 5, characterized in that, Using the fault trend model as a constraint, and utilizing the auxiliary directional structural elements and the isotropic spherical structural elements, iterative erosion and feedback correction are performed on the preliminary fault skeleton to generate an optimized fault skeleton, specifically including: A three-dimensional Hough transform was performed on the preliminary fault skeleton, and the peak value was determined in the three-dimensional parameter space to establish a regional fault trend model. Using the fault trend model as a constraint, the auxiliary directional structural elements and spherical structural elements are alternately used to perform iterative corrosion calculations on the fault skeleton prototype after the first round of corrosion, so as to generate an optimized fault skeleton. In the iterative erosion operation, the stopping criterion for iterative erosion is set as the rate of change of the number of skeleton voxels between two adjacent iterations being less than 5%; In the iterative erosion operation, after each iteration of the erosion operation, the degree of fit between the current skeleton body and the fault trend model is calculated. When the degree of fit is less than 0.7, it is determined that the skeleton body in the region deviates from the fault trend, and morphological expansion correction of spherical structural elements is performed on the region.

7. The method for extracting fault rods based on morphological erosion algorithm as described in claim 1, characterized in that, The optimized fault skeleton is subjected to three-dimensional gradient enhancement and refinement to form initial fault rods, specifically including: The gradient of the optimized fault skeleton is calculated using the three-dimensional Sobel gradient operator, and the solutions are obtained respectively. gradient components in three directions The gradient magnitude volume is obtained through vector synthesis; The gradient magnitude volume is linearly normalized, and a parallel thinning algorithm is used to thin the normalized gradient magnitude volume to generate an initial fault bar with rod-like characteristics.

8. The method for extracting fault rods based on morphological erosion algorithm as described in claim 1, characterized in that, Perform multi-attribute fusion verification on the initial fault rod to generate a verified fault rod, specifically including: Three types of auxiliary attribute volumes—post-stack seismic amplitude volume, root mean square amplitude volume, and curvature volume—were selected for optimization processing including noise suppression, smoothing, and enhancement. Spatial registration is performed between the initial fault rod and the three types of auxiliary attribute volumes to ensure consistency of three-dimensional coordinates, and the correlation coefficient is calculated for each voxel in the initial fault rod. Based on the comparison results between the correlation coefficient and the preset correlation threshold, the fault segments are eliminated and completed to generate verified fault bars.

9. The method for extracting fault rods based on morphological erosion algorithm as described in claim 8, characterized in that, After verification, the fault bars are subjected to three-dimensional connectivity analysis and topology optimization, followed by spatial interpolation processing for lattice adaptation to generate the final fault plane. Specifically, this includes: The connectivity of the verified fault rods is tracked using a three-dimensional seed filling algorithm. The spatial range of each connected branch is recorded using a 6-neighbor connectivity tracking method to obtain the geometric parameters of each connected branch. The geometric parameters include the length, orientation angle, and dip angle of each connected branch. Short, invalid branches with fewer than 10 seismic sampling points are removed based on the connected component length threshold to obtain a set of candidate connected components; The candidate connected branch set is subjected to topology simplification and optimization processing. At the intersection of fault bars, branches that match the fault trend model are retained. At the bifurcation point, small bifurcations with a width of less than 1 voxel are removed. At the twist point, the fault bar direction is corrected by smooth interpolation to generate topology-optimized fault bars. Three-dimensional kriging interpolation is used to characterize the spatial correlation of voxel values ​​through a variogram function, so as to perform interpolation gap filling on the topology-optimized fault rods. After completing the interpolation and gap filling, the grid parameter step size and origin coordinates of the structural lattice model are obtained, and the seismic sampling coordinates are converted into the node coordinates of the structural lattice model to generate the final fault plane.

10. A tomographic rod extraction system based on a morphological erosion algorithm, characterized in that, The system includes: The segmentation module is used to perform noise suppression, smoothing and enhancement optimization processing on the original ant body to generate an optimized ant body. The optimized ant body is then segmented by threshold to obtain a binary mask body. A determination module is used to determine the fault axis and design different types of morphological structural elements; wherein, the morphological structural elements include principal direction structural elements, secondary direction structural elements and isotropic spherical structural elements; The generation module is used to perform the first round of erosion operation on the binary mask using the main direction structural element to generate a fault skeleton prototype. Using the fault trend model as a constraint, the auxiliary direction structural element and the isotropic spherical structural element are used to perform iterative erosion and feedback correction on the fault skeleton prototype to generate an optimized fault skeleton. The module is used to perform three-dimensional gradient enhancement and refinement on the optimized fault skeleton to form the initial fault rod. The verification module is used to perform multi-attribute fusion verification on the initial fault rod, generate the verified fault rod, and then perform three-dimensional connectivity analysis and topology optimization, as well as spatial interpolation processing for lattice adaptation on the verified fault rod to generate the final fault plane.

Images