Unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation

Abstract

The invention discloses an unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation. The unsupervised three-dimensional acoustic logging reflector imaging interpretation method sequentially comprises the steps that underground acoustic wave full wave train data collected by a multi-pole sub-array acoustic logging instrument is obtained; wave field separation and reflected wave extraction; multidirectional migration imaging is carried out; adaptive image enhancement driven by evolutionary computation; constructing a three-dimensional point cloud; adaptive spatial clustering driven by evolutionary computation and geometric fitting and interpretation of a three-dimensional structure are carried out; according to the scheme, the processing parameters are adaptively optimized through an evolutionary computation algorithm, so that unsupervised three-dimensional imaging of the geological reflector around the well, adaptive optimization of the automatic segmentation and processing parameters and true three-dimensional quantitative interpretation of the geological structure are realized; the technical problems that parameter adjustment depends on artificial experience and cannot adapt to a dynamic noise environment, geometric distortion exists in two-dimensional interpretation, and supervised learning is difficult to use due to lack of labeled data in an existing sound wave remote detection technology are solved.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and well logging data processing technology, and in particular to an unsupervised three-dimensional acoustic well logging reflector imaging interpretation method based on evolutionary computation. Background Technology

[0002] In oil, gas, and mineral resource exploration, the detection of geological structures around the wellbore is crucial. Dipole shear wave remote logging technology radiates shear waves from a dipole acoustic source in the wellbore to the outside of the wellbore and receives the reflected echo signals from the surrounding geological bodies. It is one of the core technologies for detecting geological structures around wells, and its detection accuracy and effectiveness are widely recognized in the industry.

[0003] Existing processing workflows typically include: data acquisition, wavefield separation (time-slowness-correlation method or FK filtering), migration imaging (generating 2D profiles), image enhancement (fixed threshold / gain), and manual interpretation. However, the following drawbacks remain: (1) Image enhancement and segmentation heavily rely on manual parameter tuning, which is highly subjective and inefficient: In traditional methods, filter thresholds, gain coefficients, etc., need to be set manually. Due to the drastic changes in downhole noise with depth (e.g., high noise at the enlarged diameter), fixed parameters often have unintended consequences: high thresholds may accidentally delete micro-fractures, while low thresholds may retain artifacts. Manual parameter tuning requires repeated trial and error, and processing a single well (3000m deep) can take 8-12 hours. Moreover, due to differences in experience, the interpretation error is generally large, making it difficult to meet the needs of efficient exploration.

[0004] (2) Lack of automated quantitative interpretation methods in three-dimensional space: Existing commercial software processes well perimeter data into multiple two-dimensional variable density maps. Interpreters need to "draw lines" on the two-dimensional map to infer the three-dimensional structure, which results in geometric distortion and makes it difficult to intuitively describe the volume of non-planar structures (such as caves), leading to misjudgment of reservoir connectivity.

[0005] (3) Difficulty in implementing supervised learning methods: Existing deep learning methods rely on a large number of labeled samples, while obtaining real geological labels in deep underground exploration is extremely costly and samples are scarce, resulting in poor model generalization ability and difficulty in promotion. Summary of the Invention

[0006] To address the problems in existing acoustic remote sensing technologies, such as parameter adjustment relying on human experience, inability to adapt to dynamic noise environments, geometric distortion in two-dimensional interpretation, and difficulty in using supervised learning due to lack of labeled data, this invention provides an unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation, which achieves adaptive optimization of processing parameters and true three-dimensional quantitative interpretation of geological structures.

[0007] This invention is achieved using the following technical solution: An unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation includes the following steps: Step S1: Use an orthogonal dipole acoustic logging tool to log within the wellbore depth range. Acquire sonic logging data internally; sonic logging data is represented as a four-dimensional tensor. ,in The unit is , indicating the signal sampling time; , indicating the receiver azimuth angle; The unit is meters (m), which represents the source distance between the sound source and the receiver. Step S2: Wave field separation and reflected wave extraction. Using dispersion correction and multi-scale filtering techniques, direct waves and mode waves are suppressed, and the reflected wave field is separated. Step S3: Multi-azimuth migration imaging, performing azimuth-gated Kirchhoff integral migration imaging on the separated reflected wave data to generate a two-dimensional reflection image set covering the entire well perimeter; Step S4: Construct an evolutionary optimization model with image information entropy and local contrast as the joint fitness function, and use the optimization algorithm to adaptively search for the optimal filtering threshold and nonlinear gain parameter to enhance the two-dimensional reflection image; Step S5: Extract effective reflection points from the enhanced image based on the energy threshold, map them from the well perimeter cylindrical coordinate system to the three-dimensional Cartesian coordinate system, and construct a three-dimensional sparse point cloud; Step S6: Construct an evolutionary optimization model with clustering effectiveness as the fitness function, and use the particle swarm optimization algorithm to adaptively search for the optimal neighborhood parameters of the density clustering algorithm. Based on the core point threshold, the scattered point cloud is automatically segmented into independent geological body clusters; Step S7: Perform principal component analysis on each segmented point cloud cluster, automatically identify the structure type, and calculate its orientation, dip, tilt angle, extension length, and volumetric geometric parameters respectively.

[0008] Specifically, step S2 includes the following sub-steps: Step S21: Extract formation shear wave slowness from the raw full wave train data using STC slowness-time correlation filtering. Slowness of shear waves in formation Defined as transverse wave velocity The reciprocal of, that is The unit is It is used to characterize the ease or difficulty of propagation of shear waves in different lithological strata; Step S22: Design a sector cutoff filter in the FK frequency-wavenumber domain Preserve the apparent speed And satisfy , of which For the shear wave velocity of the formation, For signal apparent speed; Step S23: The target output after wavefield separation is a reflected wave signal in the FK frequency-wavenumber domain. Only the reflected wave component reflecting the geological structure around the well is retained for subsequent migration imaging, expressed as: ; in, This is the original full-wave train data of the acoustic wave in the frequency domain, derived from the original data in the time domain. Obtained through a two-dimensional Fourier transform, it contains all wavefield components of direct waves, mode waves, and reflected waves; The reflected wave signal is in the frequency-number domain (FK domain) and is the target output after wavefield separation; Step S24: Inverse transform back to the time domain to obtain the reflected wave data.

[0009] Specifically, step S3 includes the following sub-steps: Step S31: Employ azimuth-focused Kirchhoff integral migration technique. Divide the well perimeter into... Each directional sector; Step S32: For the first Azimuth Imaging points Amplitude at point Defined as: ; in, Indicates belonging to the first A set of receivers in each direction; Indicating two-way travel, ; Indicates the amplitude correction factor. This indicates the depth position of the receiver, where z is the depth of the imaging point and r is the radial distance. Indicates the first Each receiver at depth ,time The raw acoustic wave data collected at the location is the input data source for offset imaging; This represents the speed of sound propagation in the formation medium, and its value needs to be determined based on the geological characteristics of the formation lithology and porosity. Step S33: Output representation Two-dimensional cross-sectional images are used to form an image set. .

[0010] Specifically, step S4 is adaptive image enhancement based on evolutionary computation, which includes the following steps: Step S41: Define the optimization parameters. The parameter vector to be optimized is defined as follows: Indicates a low threshold. Indicates the exponent of the power-law transformation. This represents the Laplacian sharpening factor; the processing formula is expressed as: ; Step S42: Define the no-reference image quality evaluation function The objective is to maximize this function, expressed as: ; This represents the fitness function value, used to evaluate the image enhancement effect corresponding to the parameter vector X. The larger the value, the better the image enhancement effect; X represents the parameter vector to be optimized, and ,in For low threshold, The power-law transformation exponent. This is the Laplacian sharpening factor; These represent the normalized weight coefficients, which correspond to the information entropy. Local contrast and noise estimation The degree of contribution to the fitness function, and satisfying ; Indicates the enhanced image Information entropy ; This represents the enhanced image; Indicates the enhanced image Local contrast; Indicates the enhanced image The noise estimate is used to measure the noise level in the image. The smaller the noise estimate, the better the noise suppression effect during image enhancement. Indicates the enhanced image medium gray level is The probability of a pixel appearing; the enhanced image The total number of gray levels is expressed as ; Step S43: Evolutionary optimization, initial generation A population of parameter vectors is established; a differential mutation operation is performed on the population to generate experimental vectors; if the experimental vectors produce an image... If the value is higher, replace the original vector; iterate until the maximum algebraic number is reached, and then output the globally optimal parameters. And apply this parameter to process all Zhang cross-section.

[0011] Specifically, step S5 includes the following sub-steps: Step S51: Convert the enhanced 2D image set into 3D point cloud data; Step S52: For the first Sections in each direction any pixel in If its amplitude Then it is mapped to a point in three-dimensional space. : ; in, The third in three-dimensional Cartesian coordinate system The spatial coordinates of the reflection point and The coordinates are in the horizontal direction. Depth coordinates Radial distance from the reflection point to the well axis No. The azimuth angle of each directional sector, with a value range of [value missing]. The depth position of the reflection point is directly used as the three-dimensional coordinate system. value Azimuth The cosine and sine values ​​are used to represent the radial distance. Decomposed into horizontal directions and Quantity; Step S53: Finally obtain the 3D point cloud set Each point carries an attribute value. , which represents the intensity of reflection.

[0012] Specifically, step S6 is adaptive spatial clustering based on evolutionary computation, which includes the following steps: Step S61: Define optimization parameters, the particle position vector is defined as follows: (Epsilon) represents the neighborhood search radius. Indicates the core point threshold; Step S62: Use the silhouette coefficient SC as a clustering effectiveness index, defined as: ; Among them, SC is used to evaluate the effectiveness of clustering results, and its value range is [value range missing]. The closer the SC value is to 1, the better the clustering effect; the closer it is to -1, the worse the clustering effect. K represents the number of clusters after point cloud clustering. Indicates sample The average distance to other points in the same cluster; Indicates sample The average distance to all points in the nearest heterogeneous cluster; the objective is to maximize Value, range [-1, 1]; Step S63: Evolutionary optimization, particles are optimized based on their own history. and global optimal Update speed and location representation: ; ; in, This represents the particle's current velocity vector; This represents the inertia weight, and its value typically ranges from 0.4 to 0.9. The learning factor represents the particle's path to its historical best position. and the global optimal position of the population The weight of learning; Represents a random number within the interval [0,1]. This represents the vector of the particle's own historical best position. Represents the global optimal position vector of the population; This represents the particle's current position vector; This represents the updated velocity vector of the particle; This represents the updated position vector of the particle; in each generation, the current position vector is used. Point cloud Run the DBSCAN algorithm; calculate the clustering results. The value is used as the fitness; the final output is the optimal clustering result, which yields... Cluster .

[0013] Specifically, step S7 applies to each point cloud cluster. Geometric morphology analysis is performed, specifically including: Morphological discrimination: Calculate the covariance matrix of the clusters and its eigenvalues ;like It is determined to be a planar structure; if It is determined to be a body-like structure; if It was determined to be a linear structure; Calculation of attitude of planar structures: Constructing plane equations The spatial morphology of the fitted planar structure, in which These are the components of the plane normal vector. , For constant terms; inclination angle The angle between the planar structure and the well axis reflects the inclination of the planar structure, expressed through the plane normal vector. The angle with the Z-axis is calculated using the following formula: The range of values ​​is ; Volume calculation of volumetric structures: Calculating the volume of volumetric structures. The volumetric shape of the point cloud is constructed using the AlphaShapes algorithm or the convex hull algorithm, and then the volume of the polyhedron is calculated.

[0014] Specifically, the covariance matrix of the clusters in the morphological discrimination The calculations include: For point cloud clusters all points First, calculate the centroid of the point cloud, expressed as: ; in, ; The number of points within the cluster; Construct the covariance matrix, represented as: ; in, , They are respectively Any coordinate component in the matrix; Define the criteria for judging similar eigenvalues: if Then it is determined to be .

[0015] The beneficial effects of this invention are as follows: (1) Full-process automation: Through the evolutionary calculation global optimization capability of step S4 (DE algorithm to optimize image enhancement parameters) and step S6 (PSO algorithm to optimize clustering parameters), the traditional manual trial and error parameter tuning is replaced, human experience bias is eliminated, the standard deviation of single well processing results is reduced, and the processing results are standardized. (2) Strong environmental adaptability: The algorithm estimates the noise in the fitness function in step S4. Real-time sensing of data signal-to-noise ratio changes; automatic adjustment in high-noise well sections such as wellbore enlargement. Tighten the control to suppress artifacts from the maximum amplitude; automatically adjust in weak signal well sections. It enhances weak reflection signals and increases the effective reflection point extraction rate, which is significantly higher than that of the fixed parameter method. (3) True 3D quantitative interpretation: Geometric fitting is performed directly in the 3D point cloud space, which solves the distortion problem of 2D interpretation and can accurately quantify the volume of the cave. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the structures shown in these drawings without creative effort.

[0017] Figure 1 This is a flowchart of the unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation in an embodiment of the present invention; Figure 2 This is a schematic diagram of the evolutionary calculation iteration process in an embodiment of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0019] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0020] The following is in conjunction with the appendix Figure 1-2 The following describes some embodiments of the present invention in detail. Unless otherwise specified, the following embodiments and features can be combined with each other.

[0021] This invention proposes an unsupervised three-dimensional acoustic logging reflector imaging interpretation method based on evolutionary computation. In a preferred embodiment, the method is as follows: Figure 1 As shown, it includes: Step S1: Acquire the full waveform data of downhole acoustic waves collected by the multipole array acoustic logging instrument. The data format is defined as a four-dimensional tensor. ,in: (Unit: m), indicating the range of wellbore depth; (unit: ), indicating the signal sampling time; , indicating the receiver azimuth angle; (Unit: m), representing the source distance between the sound source and the receiver.

[0022] Step S2: Wavefield separation and reflected wave extraction. Using dispersion correction and multi-scale filtering techniques (such as time-slowness-correlation method combined with FK filtering), direct waves and mode waves are suppressed, and the reflected wavefield is separated.

[0023] Step S3: Multi-azimuth migration imaging. Azimuth-gated Kirchhoff integral migration imaging is performed on the separated reflected wave data to generate a two-dimensional reflected wave image set covering the entire well perimeter. Number of azimuth sectors. 16-72 (preferably 36, per sector) ) Step S4: Evolutionary Computation-Driven Adaptive Image Enhancement. An evolutionary optimization model is constructed with image information entropy and local contrast as the joint fitness function. An optimization algorithm adaptively searches for the optimal filtering threshold and nonlinear gain parameters to enhance the two-dimensional reflective image. This step replaces manual parameter tuning and automatically adjusts parameters based on the signal-to-noise ratio. The preferred optimization algorithm is Differential Evolutionary Algorithm (DE), but a Genetic Algorithm (GA) can also be used; the vector of parameters to be optimized is defined as follows: The value range and physical meaning of each parameter are as follows: Low threshold Used to suppress background noise; power-law transform index Used for nonlinear enhancement of weakly reflected signals; Laplace sharpening coefficient Used to enhance the edges of geological structures along the same phase axis. The fitness function comprehensively considers image information entropy, local contrast, and noise level.

[0024] Step S5: 3D point cloud construction. Effective reflection points are extracted from the enhanced image based on energy thresholding, and mapped from the wellbore cylindrical coordinate system to a 3D Cartesian coordinate system to construct a 3D sparse point cloud.

[0025] Step S6: Evolutionary Computation-Driven Adaptive Spatial Clustering. An evolutionary optimization model is constructed with a clustering effectiveness metric (such as silhouette coefficient) as the fitness function. The particle swarm optimization algorithm is used to adaptively search for the optimal neighborhood parameters of the density clustering algorithm (DBSCAN). The system uses a core point threshold and other parameters to automatically segment scattered point clouds into independent geological clusters. Optimization algorithms include Particle Swarm Optimization (PSO) or Ant Colony Optimization (ACO); the parameter vector to be optimized is used. Where: neighborhood search radius : 0.1-2.0 meters (preferably 0.3-0.8 meters), core point threshold : 3-20 points (preferably 5-10 points), fitness function: silhouette coefficient.

[0026] Step S7: 3D structural geometric fitting and interpretation. Principal component analysis is performed on each segmented point cloud cluster to automatically determine the structure type (area, volume, line) and calculate its geometric parameters such as orientation, dip, tilt angle, extension length, and volume.

[0027] In one specific embodiment, the steps are as follows: Step S1: Acquiring acoustic logging data Using an orthogonal dipole acoustic logging tool, in the depth range Data is collected internally. The data is represented as a four-dimensional tensor. ,in For depth, Sampling time, For the receiver azimuth angle, The source distance.

[0028] Step S2: Data Preprocessing (Wavefield Separation) In order to extract the reflected wave from the full wave train It is necessary to remove the sliding waves (direct P-waves and S-waves) and Stoneley waves propagating along the wellbore. A combination of slow-time correlation (STC) filtering and frequency-wavenumber (FK) filtering is used. 1. Extracting formation shear wave slowness using STC Formation shear wave slowness is defined as shear wave velocity. The reciprocal of ( ), the unit is usually The parameter is used to characterize the ease of propagation of shear waves in different lithological strata and is a key parameter for designing FK domain filters. It is extracted from the original full wave train data through STC analysis.

[0029] 2. Design a sector cutoff filter in the FK domain Preserve the apparent speed satisfy The signal components (reflected waves typically appear as hyperbolic in-phase axes within a specific apparent velocity range) include... For the shear wave velocity of the formation, The apparent speed of the signal. ; The reflected wave signal in the frequency-wavenumber domain (FK domain) is the target output after wavefield separation, retaining only the reflected wave components reflecting the surrounding geological structure for subsequent migration imaging. This is derived from the raw acoustic full-wave train data in the frequency-wavenumber domain. With sector cut-off filter It is obtained through multiplication. And the parameter... The corresponding frequency is the result of a Fourier transform in the time dimension, primarily used to describe the oscillation characteristics of a signal in the time dimension; parameters The corresponding wavenumber is the result of a Fourier transform of the spatial dimension, used to characterize the wave characteristics of a signal in the spatial dimension. It is related to the wavelength. Between The relationship. : The original full-wavelength acoustic data in the frequency domain, derived from the original data in the time domain. It is obtained by two-dimensional Fourier transform and contains all wave field components, including direct waves, mode waves and reflected waves.

[0030] 3. The reflected wave data is obtained by inverse transformation back to the time domain.

[0031] Step S3: Multi-azimuth migration imaging A Kirchhoff integral migration technique with azimuth focusing is employed. The well perimeter is divided into... Each directional sector (e.g.) Each sector For the first... Azimuth Imaging points Amplitude at point for: ; in : Belongs to the A set of receivers in each direction. When traveling in both directions, . Amplitude correction factor, relative to receiver depth The imaging point depth z is related to the radial distance r, which is used to compensate for the amplitude loss of sound waves during propagation due to factors such as distance attenuation and medium absorption, so that the amplitude of imaging points at different locations is comparable. : No. Each receiver at depth ,time The raw acoustic data collected at the site is the input data source for migration imaging, containing raw signal information of the downhole geological structure. The depth position of the receiver is an important parameter for locating the receiver's spatial coordinates and is used to determine the spatial source of the signal in data acquisition and subsequent imaging calculations. The speed of sound propagation in the formation medium is used to calculate the two-way travel time. The core parameters, whose values ​​need to be determined based on geological characteristics such as strata lithology and porosity, directly affect the accuracy of calculating the depth and radial distance of the imaging point. Output: Two-dimensional cross-sectional images are used to form an image set. .

[0032] Step S4: Adaptive Image Enhancement Based on Evolutionary Computation like Figure 2 This step aims to automatically improve the signal-to-noise ratio of the image set without relying on manual intervention. This embodiment uses the Differential Evolution (DE) algorithm.

[0033] 1. Optimization parameter definition (gene): The parameter vector to be optimized. . : Low threshold, used to set background noise to zero. : Power-law transformation (GammaCorrection) exponent, used for nonlinear enhancement of weakly reflected signals. Laplacian sharpening factor, used to enhance the edges of the in-phase axis. Processing formula: .

[0034] 2. Fitness Function Design: To evaluate the enhancement effect, a no-reference image quality evaluation function is defined. The goal is to maximize this function: ; The fitness function value is used to evaluate the image enhancement effect corresponding to the parameter vector X. It is a core indicator for selecting the optimal parameters in evolutionary algorithms. A larger value indicates a better image enhancement effect. X: The vector of parameters to be optimized, specifically... ,in For low threshold, For the power-law transformation exponent, The Laplacian sharpening coefficient is the object of optimization in evolutionary algorithms, directly determining the processing method and effect of image enhancement. Normalized weight coefficients, corresponding to information entropy respectively. Local contrast and noise estimation The contribution of the signal to the fitness function needs to be adjusted according to the requirements of image information content, signal saliency, and noise suppression in practical applications, and must satisfy the following: (Normalization conditions). Enhanced image Information entropy is used to measure the amount of information contained in an image (texture richness). The higher the information entropy, the richer the detailed information related to geological structures in the image. The calculation method is... ). The enhanced image is the output of image enhancement processing. Its quality (information content, contrast, noise level) is determined by a fitness function. To conduct an evaluation. Enhanced image The local contrast is obtained by calculating the variance of local regions of the image. It is used to evaluate the significance of effective quality signals. The greater the local contrast, the more obvious the distinction between geological structures and background, which facilitates the subsequent extraction of reflection points. Enhanced image The noise estimate is obtained by calculating the energy proportion of the high-frequency components in the image. It is used to measure the noise level in the image. The smaller the noise estimate, the better the noise suppression effect during the image enhancement process. Enhanced image medium gray level is The probability of a pixel appearing is used to calculate information entropy. The basic parameters are obtained by statistically analyzing the proportion of pixels at each gray level in the image. Enhanced image The total number of gray levels, that is, the number of different gray values ​​that a pixel in an image can have, is used to calculate information entropy. The upper limit of the summation range depends on the number of bits quantized in the image (e.g., L=256 for an 8-bit image).

[0035] 3. Evolutionary optimization process: Initialization: Generation A population of parameter vectors. Mutation and crossover: Differential mutation is performed on the population to generate experimental vectors. Selection: If the experimental vectors produce an image... If the value is higher, replace the original vector. Convergence: Iterate to the maximum algebra and output the globally optimal parameters. And apply this parameter to process all Zhang cross-section.

[0036] Step S5: 3D Point Cloud Construction The enhanced 2D image set is converted into 3D point cloud data. For the first... Sections in each direction any pixel in If its amplitude Then it is mapped to a point in three-dimensional space. : ; The third in three-dimensional Cartesian coordinate system The spatial coordinates of the reflection points, where and The coordinates are in the horizontal direction. Depth-direction coordinates (relative to wellbore depth) Consistency is a core parameter for constructing 3D point clouds, used to characterize the spatial location of reflection points. The radial distance from the reflection point to the well axis is a key parameter for the transformation from a two-dimensional profile (cylindrical coordinate system) to a three-dimensional Cartesian coordinate system, reflecting the horizontal distance between the reflector and the wellbore. : No. The azimuth angle of each directional sector, with a value range of [value missing]. The position of the reflection point on the horizontal plane is used to determine the orientation of the reflection point and is an important parameter for stitching together multi-directional two-dimensional profiles into a three-dimensional point cloud. The depth position of the reflection point is consistent with the depth coordinates in the two-dimensional profile, and is directly used as the depth coordinate in the three-dimensional coordinate system. This value ensures the consistency and accuracy of the depth information of the reflection point. azimuth The cosine and sine values ​​are in cylindrical coordinates (radial distance). Azimuth ,depth To the three-dimensional Cartesian coordinate system The trigonometric function factor of the transformation is used to convert the radial distance. Decomposed into horizontal directions and Components. The final result is a 3D point cloud set. Each point carries an attribute value. (Reflection intensity).

[0037] Step S6: Adaptive Spatial Clustering Based on Evolutionary Computation This step aims to transform the discrete point cloud Automatic clustering is performed on the DBSCAN clusterer as independent geological bodies. In this embodiment, the Particle Swarm Optimization (PSO) algorithm is used to optimize the DBSCAN clusterer.

[0038] 1. Optimization parameter definition (particle position): Particle position vector . (Epsilon): Neighborhood search radius, which physically represents the maximum distance between points within the same geological body. Core point threshold, physically meaning the minimum number of reflection points required to form an effective structure (noise reduction parameter).

[0039] 2. Fitness function design: The Silhouette Coefficient (SC) is used as the clustering effectiveness index. SC: Used to evaluate the validity of clustering results, with a value range of [value missing]. The closer the SC value is to 1, the better the clustering effect (closer samples within clusters and more separated samples between clusters); the closer it is to -1, the worse the clustering effect. K: The number of clusters after point cloud clustering, that is, the number of independent geological bodies segmented by the DBSCAN algorithm into the discrete point cloud. It is the upper limit of the summation range in the contour coefficient calculation and reflects the complexity of the downhole geological structure. :sample The average distance to other points in the same cluster (the closer the cluster, the higher the average distance to other points in the cluster). The smaller). :sample The average distance to all points in the nearest heterogeneous cluster (the more separated the clusters, the better). (The larger). Goal: Maximize Value (range [-1, 1]).

[0040] 3. Evolutionary optimization process: Particles optimize based on their own historical best. and global optimal Update speed and location: . v: The particle's current velocity vector, corresponding to the parameter to be optimized. The adjustment step size (and MinPts), whose dimension is consistent with the parameter vector Y, is used to control the search direction and amplitude of the particle in the parameter space. Inertial weight, typically ranging from 0.4 to 0.9, is used to balance the global exploration capability of particles (large). ) and local development capabilities (small In the initial stage, a larger value is used to expand the search range, and in the later stage of iteration, a smaller value is used to finely converge to the optimal solution. The learning factor (acceleration factor), usually set to 2.0, represents the particle's trajectory towards its historical best position. and the global optimal position of the population The weight of learning determines the degree to which a particle is influenced by its own experience and the experience of the group. The random number in the interval [0,1] is used to increase the randomness of particle search and avoid the algorithm from getting stuck in local optima. It is generated randomly in each iteration to ensure the diversity of the search process. The particle's own historical best position vector records the parameter combination corresponding to the highest fitness (maximum silhouette coefficient SC) obtained by the particle in all previous iterations, guiding the particle to update towards its own optimal direction. : The global optimal position vector of the population, recording the parameter combination that achieves the highest fitness across all iterations, guiding all particles to update towards the optimal direction for the population. Y: The current position vector of the particle, i.e., the DBSCAN clustering parameters to be optimized ( The validity of clustering results is directly determined by the clustering parameters (MinPts). The updated velocity vector of the particle, obtained by combining inertia, self-learning, and group learning, provides the step size for the next round of position update. The updated position vector of the particle, i.e., the new DBSCAN clustering parameter combination, is used for the next round of clustering calculation and fitness evaluation. In each generation, the current position vector is used. Point cloud Run the DBSCAN algorithm. Calculate the clustering results. The value is used as the fitness. The final output is the optimal clustering result. Cluster .

[0041] Step S7: 3D structural interpretation and geometric parameter calculation For each point cloud cluster Perform geometric morphological analysis: 1. Morphological discrimination (based on PCA): Calculate the covariance matrix of the clusters. and its eigenvalues .like (Extending significantly in two directions and less so in one direction) is identified as a planar structure (fracture, fault). If (With similar extensions in three directions), it is determined to be a body-like structure (cavity). If (Extending significantly in one direction), identified as a linear structure (drill string scratches or narrow channels around the wellbore). Cluster covariance matrix. Calculation method: For point cloud clusters all points First, calculate the centroid of the point cloud. ( , (where the number of points is within the cluster), then construct the covariance matrix. ,in ( They are respectively (Any coordinate component in the equation). Additionally, a criterion for determining "similar" eigenvalues ​​is provided: if... Then it is determined to be This is to avoid misjudging the construction type due to vague judgment criteria.

[0042] 2. Calculation of attitude of planar structures: Plane equations : Used to fit the spatial morphology of planar structures, where The components of the plane normal vector ( ), The constant term, obtained by fitting the point cloud of the planar structure using the least squares method, forms the basis for calculating the dip angle, dip direction, and other attitude parameters of the planar structure. Dip angle The angle between the planar structure and the well axis (Z-axis) reflects the degree of inclination of the planar structure, expressed through the plane normal vector. The angle with the Z-axis is calculated using the formula: The range of values ​​is .

[0043] 3. Volume calculation of body-shaped structures: Volume of body-shaped structures (such as caves) The volumetric structure is obtained by constructing the envelope polyhedron of the point cloud of the volumetric structure using the AlphaShapes algorithm or the ConvexHull algorithm, and then calculating the volume of the polyhedron. It is a key parameter for quantitatively describing the scale of the volumetric structure and can be directly used for reservoir evaluation and resource estimation.

[0044] For the foregoing embodiments, in order to simplify the description, they are all described as a series of actions. However, those skilled in the art should understand that this application is not limited to the described order of actions, because according to this application, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to this application.

[0045] The above embodiments describe the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Modifications and variations made by those skilled in the art without departing from the spirit and scope of the invention should be within the protection scope of the appended claims.

Claims

1. An unsupervised 3D sonic well-logging reflector imaging interpretation method based on evolutionary computation, characterized in that, The method comprises the following steps: Step S1: Use an orthogonal dipole acoustic logging tool to log within the wellbore depth range. Acquire sonic logging data internally; sonic logging data is represented as a four-dimensional tensor. ,in The unit is , indicating the signal sampling time; , indicating the receiver azimuth angle; The unit is meters (m), which represents the source distance between the sound source and the receiver. Step S2: wave field separation and reflected wave extraction, dispersion correction and multi-scale filtering technology are used to suppress direct wave and mode wave, and separate out reflected wave field; Step S3: multi-azimuth migration imaging, azimuth gating Kirchhoff integral migration imaging is performed on the separated reflected wave data to generate a two-dimensional reflection image set covering the full azimuth around the well; Step S4: an evolutionary optimization model is constructed with image information entropy and local contrast as a joint fitness function, and an optimization algorithm is used to adaptively search for optimal filtering threshold and nonlinear gain parameters to enhance the two-dimensional reflection image; Step S5: based on the energy threshold, effective reflection points in the enhanced image are extracted, which are mapped from the well column coordinate system to the three-dimensional Cartesian coordinate system to construct a three-dimensional sparse point cloud; Step S6: Constructing an evolutionary optimization model with the cluster validity index as the fitness function, and using the particle swarm algorithm to adaptively search for the optimal neighborhood parameters of the density clustering algorithm and the core point threshold value, the scattered point cloud is automatically segmented into independent geological body clusters. Step S7: principal component analysis is performed on each segmented point cloud cluster to automatically identify the structure type, and the trend, tendency, dip angle, extension length and volume geometric parameters are calculated respectively.

2. The evolutionary computation based unsupervised 3D sonic well- log reflectivity body imaging interpretation method of claim 1, wherein, The step S2 specifically comprises the following sub-steps: Step S21: Extract formation shear wave slowness from the raw full wave train data using STC slowness-time correlation filtering. Slowness of shear waves in formation Defined as transverse wave velocity The reciprocal of, that is The unit is It is used to characterize the ease or difficulty of propagation of shear waves in different lithological strata; Step S22: design a sector cut filter in F-K frequency-wavenumber domain , retaining the signal components with apparent velocity and satisfying where is the S-wave velocity of the formation, is the apparent velocity of the signal; Step S23: the target output after wave field separation is the reflected wave signal in the F-K frequency-wavenumber domain, only the reflected wave component reflecting the geological structure around the well is retained, which is used for subsequent migration imaging, and is expressed as: ； wherein, is the original acoustic full wavefield data in the frequency-wavenumber domain, which is derived from the time-domain original data is obtained by two-dimensional Fourier transform, and contains all wavefield components of direct wave, mode wave and reflected wave; is the reflected wave signal in the frequency-wavenumber domain F-K domain, which is the target output after wavefield separation; Step S24: inverse transform back to the time domain to obtain reflected wave data.

3. The evolutionary computation based unsupervised 3D sonic well- log reflectivity body imaging interpretation method of claim 2, wherein, The step S3 specifically comprises the following sub-steps: Step S31: Kirchhoff integral migration with azimuth focusing. The well is divided into azimuth sectors. Step S32: For the first Azimuth Imaging points Amplitude at point Defined as: ； wherein, represents a set of receivers belonging to the th azimuth; represents the two-way travel time, ; represents the amplitude correction factor, represents the depth position of the receiver, z is the depth of the imaging point, and r is the radial distance; represents the original acoustic wave data collected by the th receiver at depth , time , which is the input data source of the migration imaging; represents the acoustic wave propagation speed in the formation medium, which needs to be determined according to the geological characteristics of the formation lithology and porosity. Step S33: outputting the representation a two-dimensional profile view, forming an image set .

4. The evolutionary computation based unsupervised 3D sonic well- log reflectivity body imaging interpretation method of claim 3, wherein, The step S4 is adaptive image enhancement based on evolutionary calculation, specifically comprising the following steps: Step S41: defining optimization parameters, the parameter vector to be optimized is defined as denotes a low threshold, denotes a power law transform exponent, denotes a Laplacian sharpening coefficient; the processing formula is represented as: ； Step S42: defining the no-reference image quality assessment function The objective is to maximize this function, denoted as: ； represents fitness function value, used to evaluate the image enhancement effect corresponding to parameter vector X, The greater the value, the better the image enhancement effect; X represents the parameter vector to be optimized, and , wherein is a low threshold, is a power law transformation index, is a Laplace sharpening coefficient; represents a normalized weight coefficient, corresponding to the information entropy , local contrast and noise estimation in the fitness function respectively, and satisfies ; represents the information entropy of the enhanced image ; ; represents the enhanced image; represents the local contrast of the enhanced image ; represents the noise estimation value of the enhanced image , used to measure the noise level in the image, the smaller the noise estimation value, the better the noise suppression effect in the image enhancement process; represents the probability of the pixel with gray level in the enhanced image ; the total number of gray levels of the enhanced image is represented as ; Step S43: evolution optimization, initialization generation of parameter vector population; differential mutation operation is performed on the population to generate test vector; if the image value generated by the test vector is higher, the original vector is replaced; iteration is performed to the maximum number of iterations, and the global optimal parameters are output and the parameters are applied to process all​​ 5. The evolutionary computation based unsupervised 3D sonic well- log reflectivity body imaging interpretation method of claim 4, wherein, The step S5 specifically comprises the following sub-steps: Step S51: convert the enhanced two-dimensional image set into three-dimensional point cloud data; Step S52: For the first Sections in each direction any pixel in If its amplitude Then it is mapped to a point in three-dimensional space. : ； in, The third in three-dimensional Cartesian coordinate system The spatial coordinates of the reflection point and The coordinates are in the horizontal direction. Depth coordinates Radial distance from the reflection point to the well axis No. The azimuth angle of each directional sector, with a value range of [value missing]. The depth position of the reflection point is directly used as the three-dimensional coordinate system. value Azimuth The cosine and sine values ​​are used to represent the radial distance. Decomposed into horizontal directions and Quantity; Step S53: final set of three-dimensional point clouds is obtained each point carrying a value of an attribute representing the intensity of reflection.

6. The evolutionary computation based unsupervised 3D sonic well- log reflectivity body imaging interpretation method of claim 5, wherein, The step S6 is adaptive spatial clustering based on evolutionary calculation, specifically comprising the following steps: Step S61: define optimization parameters, the particle position vector is defined as (Epsilon) represents the neighborhood search radius, represents the core point threshold value; Step S62: the contour coefficient SC is used as a clustering effectiveness index, which is defined as: ； Among them, SC is used to evaluate the effectiveness of clustering results, and its value range is [value range missing]. The closer the SC value is to 1, the better the clustering effect; the closer it is to -1, the worse the clustering effect. K represents the number of clusters after point cloud clustering. Indicates sample The average distance to other points in the same cluster; Indicates sample The average distance to all points in the nearest heterogeneous cluster; the objective is to maximize Value, range [-1, 1]; Step S63: Evolutionary optimization, particles according to their own history optimal and global optimal Update speed and position representation: ； ； where, represents the current velocity vector of the particle; represents the inertia weight, usually in the range of 0.4-0.9; represents the learning factor, respectively represents the particle to its own historical optimal position and the global optimal position of the population weight of learning; represents a random number in the interval [0, 1]; represents the particle's own historical optimal position vector; represents the global optimal position vector of the population; represents the current position vector of the particle; represents the updated velocity vector of the particle; represents the updated position vector of the particle; in each generation, the current point cloud runs the DBSCAN algorithm; calculates the value of the clustering result as the fitness; finally outputs the optimal clustering result, obtains clusters .

7. The evolutionary-computing-based, unsupervised, 3D sonic well- log reflector imaging interpretation method of claim 6, wherein, The step S7 is performed for each point cloud cluster carrying out geometric shape analysis, specifically including: Shape discrimination: Compute the covariance matrix of the cluster and its eigenvalues If , decide on planar structure; if , decide on volumetric structure; if , decide on linear structure. Facetted structure occurrence calculation: construction of plane equation , fitting the spatial form of the facetted structure, wherein are components of the plane normal vector, , is a constant term; dip angle is the angle between the facetted structure and the well axis, reflecting the inclination degree of the facetted structure, calculated by the angle between the plane normal vector and the Z axis, the formula is: , the value range is ; Body shaped construct volume calculation: calculating the volume of a body shaped construct The volume of the convex hull is calculated by the AlphaShapes algorithm or the convex hull algorithm.

8. The evolutionary-computing-based, unsupervised, 3D sonic well- log reflector imaging interpretation method of claim 7, wherein, The calculation of the covariance matrix of the clusters in the shape discrimination includes: For all points in a point cloud cluster , the centroid of the point cloud is first computed and denoted as:​ ； wherein ; is the number of points within the cluster; The covariance matrix is constructed, which is expressed as: ； wherein , are any of the coordinate components in , respectively; The determination criterion of the similar eigenvalue is defined as follows: if , it is determined that .

Images