A geological three-dimensional exploration image reconstruction method

Abstract

This invention belongs to the field of image processing technology, specifically relating to a method for reconstructing geological 3D exploration images. The method includes: calculating the confidence weight of data points, adaptive adjustment parameters, and geological complexity to determine the contribution weight of each data point to a grid node; selecting neighboring points of the grid node and constructing an initial 3D grid based on geological attribute values; calculating the gradient vector of each node in the initial 3D grid, constructing a structure tensor matrix and decomposing it to obtain eigenvalues ​​and eigenvectors; dividing the structure into structural regions based on the eigenvalues; calculating the relative weights and normalization adjustment factors in each direction; combining a preset region type factor to obtain the smoothing adjustment coefficients of each grid node in each direction; constructing a diffusion matrix from the smoothing adjustment coefficients in each direction and the corresponding eigenvectors; iteratively optimizing the initial 3D grid using anisotropic diffusion technology to obtain the reconstructed geological 3D exploration image. This invention improves the accuracy of geological 3D exploration image reconstruction.

Description

Technical Field

[0001] This invention relates to the field of image processing technology. More specifically, this invention relates to a method for reconstructing three-dimensional geological exploration images. Background Technology

[0002] Geological 3D exploration image reconstruction is a core technology in geological resource exploration and is applied in fields such as underground structure analysis. As geological structures become more complex, traditional 2D exploration images are unable to fully reflect the spatial distribution characteristics of geological bodies, while 3D reconstruction technology can intuitively present the three-dimensional morphology and fault distribution of geological bodies, providing more comprehensive and accurate data support for various geological decisions.

[0003] Currently, traditional radial basis function interpolation algorithms are commonly used for the reconstruction of 3D geological exploration images. The specific operations include: acquiring discrete exploration data and using it as control points; selecting a preset radial basis function, such as a Gaussian function, as the kernel function, the value of which depends only on the distance between control points; determining the weight coefficient corresponding to each control point by solving a system of linear equations based on all control points; and using these weight coefficients and the kernel function to calculate the data value at any unknown location in the entire 3D space, thereby obtaining the reconstructed 3D geological exploration image.

[0004] Traditional radial basis function interpolation methods have inherent limitations: First, in the kernel function selection and parameter setting steps, the parameters mostly rely on manual experience and lack the ability to adaptively match geological structural features. This can easily lead to overly smooth interpolation results in complex structural regions, specifically manifested as blurred fault boundaries and loss of lithological stratification details. Second, in the overall process of constructing and solving the interpolation model, its mechanism relies solely on the pure mathematical mapping of discrete data points and fails to effectively incorporate prior knowledge of geological structures. Consequently, the reconstructed 3D geological exploration image deviates significantly from the actual geological structure. Summary of the Invention

[0005] To address the technical problems of traditional methods that rely on manual experience to set parameters and do not incorporate prior geological knowledge, resulting in overly smoothed reconstructed images in complex areas, blurred faults, and significant deviations from the actual structure, this invention provides solutions in the following aspects.

[0006] In a first aspect, the present invention provides a method for reconstructing a three-dimensional geological exploration image, comprising: preprocessing exploration data of an exploration area to obtain an original exploration dataset, and obtaining the confidence weight, adaptive adjustment parameter, and geological complexity of data points, used to calculate the contribution weight of each data point to the interpolation grid node; filtering a set of neighboring points based on the Euclidean distance between the interpolation grid node and each data point in the original exploration dataset; performing a weighted summation of the geological attribute values ​​of all data points in the neighboring point set based on the contribution weight to obtain the geological attribute values ​​of each interpolation grid node, used to construct an initial three-dimensional grid; and calculating the gradient of each grid node in the initial three-dimensional grid using a three-dimensional Sobel operator. The structure tensor matrix is ​​constructed using degree vectors and eigenvalue decomposition is performed to obtain the eigenvalues ​​and eigenvectors of each grid node. Based on the relative magnitude of the eigenvalues, layered structure regions, fractured structure regions, and homogeneous structure regions are divided. Based on the eigenvalues, the relative weights and normalization adjustment factors in each direction are calculated. Combined with the preset region type factor, the smoothing adjustment coefficients of each grid node in each direction are obtained. Based on the smoothing adjustment coefficients of each grid node in each direction and the corresponding eigenvectors, the diffusion matrix of each grid node is constructed. Anisotropic diffusion technology is used to iteratively optimize the initial 3D grid to obtain the final 3D grid. The final 3D grid is then mapped through 3D visualization to obtain the final reconstructed 3D geological exploration image.

[0007] This invention assigns differentiated contribution weights to data points with varying reliability through a reliability-weighted calculation. It dynamically adjusts the distance attenuation effect based on geological complexity and adaptive adjustment parameters, improving the accuracy of geological attribute value calculations for interpolated grid nodes. By pre-setting a search radius to filter neighboring point sets, it reduces computational data and improves the efficiency of initial 3D grid construction. Pre-set thresholds quantify and define layered, fractured, and homogeneous structural regions, enabling accurate identification and division of geological structures and providing a basis for differentiated grid optimization. Based on the relative weights calculated from eigenvalues ​​and normalization adjustment factors, it combines regional type factors to obtain adaptive smoothing adjustment coefficients. By fusing feature vectors to construct a diffusion matrix, it implements anisotropic diffusion iterative optimization, smoothing and reducing noise in homogeneous structural regions while preserving boundary features in layered and fractured structural regions to avoid structural distortion. Ultimately, it generates a geological 3D exploration image with clear layers and accurate attribute distribution, providing highly reliable visualization support for geological exploration analysis and resource development decisions.

[0008] Preferably, obtaining the confidence weight of the data points includes: In the formula, Data points Credibility weight; Data points The normalized signal-to-noise ratio has a weighting factor of 0.6; Data points The stability coefficient, with a weighting factor of 0.4; It is the index of data points in the original exploration dataset.

[0009] This method calculates the reliability weight of data points by weighting and fusing normalized signal-to-noise ratio (SNR) and stability coefficient. The weighting coefficient for SNR is set to 0.6, and the weighting coefficient for stability coefficient is set to 0.4. This approach highlights the dominant role of data signal quality, reflected by the normalized SNR, in the reliability of data points, while also incorporating the neighborhood data consistency characteristics reflected by the stability coefficient. This makes the calculation of reliability weights more comprehensive and the quantification results more accurate. This calculation method can effectively distinguish the reliability differences of different data points in the original exploration dataset, providing a scientific and reasonable basis for subsequently calculating the contribution weight of data points to the interpolation grid nodes. It avoids interference from low-reliability data points in the interpolation calculation, thereby improving the accuracy of the calculation of geological attribute values ​​of the interpolation grid nodes.

[0010] Preferably, the method for obtaining the stability coefficient includes: obtaining the geological attribute values ​​of all data points in the original exploration dataset based on the physical survey results; and using the data points... Set as center The three-dimensional neighborhood will be used to store data points. The normalized value of the ratio of a geological attribute value to the average of the geological attribute values ​​of all data points in its three-dimensional neighborhood is used as the data point's value. The stability coefficient.

[0011] Preferably, calculating the contribution weight of each data point to the interpolation grid node includes: In the formula, Data points Mesh nodes for interpolation Contribution weight; Data points With the grid nodes to be interpolated The Euclidean distance between them; It is the distance attenuation parameter; It is an adaptive adjustment parameter; It is the degree of geological complexity; Data points Credibility weight; It is a natural exponential function; It is the index of data points in the original exploration dataset.

[0012] This method constructs a contribution weight calculation model by integrating data point reliability weights and spatial location features. Reliability weights are used as the core influencing factor, while Euclidean distance is introduced in conjunction with a distance attenuation parameter to construct a distance attenuation term. The influence of this attenuation term is adjusted by adaptively adjusting parameters in conjunction with geological complexity. This ensures that the contribution weight of a data point to the interpolation grid node is directly related to the reliability of the data point itself, adapts to the spatial distance between the data point and the interpolation grid node, and dynamically optimizes the weight allocation rules based on the geological complexity of the exploration area. This calculation method effectively weakens the interference of low-reliability data points and distant data points on the interpolation calculation, strengthens the contribution ratio of high-reliability, nearby data points, and makes the allocation of contribution weights more scientific and reasonable, providing reliable support for the accurate calculation of the geological attribute values ​​of the subsequent interpolation grid nodes.

[0013] Preferably, the step of filtering the set of neighboring points based on Euclidean distance includes: calculating the grid nodes to be interpolated. The Euclidean distance between each data point in the original exploration dataset and the grid nodes to be interpolated. Centered on the target data point, all data points whose Euclidean distance is less than the preset search radius are used as grid nodes for interpolation. Obtain the grid nodes to be interpolated from the nearest neighboring points. The set of neighboring points.

[0014] Preferably, the step of dividing the layered structure region, fractured structure region, and uniform structure region according to the relative magnitude of feature values ​​includes: mesh nodes. The eigenvalues ​​are respectively ,in, Based on the actual needs of geological exploration, preset ratio thresholds, first preset difference thresholds, and second preset difference thresholds are set in advance; when and The ratio is greater than the preset ratio threshold, and and When the difference is less than a first preset difference threshold, the grid node is determined to be located in a layered structure region; when and The difference is less than the first preset difference threshold, and and When the ratio is greater than a preset ratio threshold, the grid node is determined to be located in a fractured structural region; when and The difference is less than the second preset difference threshold. and When the difference is less than the second preset difference threshold, the grid node is determined to be located in a uniform structure region.

[0015] This method employs a quantitative judgment approach that combines a preset ratio threshold and a difference threshold, based on feature values. Based on the inherent sorting rules, this method specifically designs the judgment logic for different geological structure regions, achieving precise definition of layered structure regions, fracture structure regions, and homogeneous structure regions. This method replaces the vague description of the relative magnitude of feature values, making the region division criteria clearer, quantifiable, and repeatable, effectively avoiding division errors caused by subjective judgment. At the same time, the threshold can be flexibly set based on the actual needs of geological exploration, adapting to the structural recognition needs of different exploration scenarios. This provides a clear and accurate structural basis for subsequent acquisition of differentiated smoothing adjustment coefficients based on region type and construction of diffusion matrices, ensuring the pertinence and effectiveness of subsequent anisotropic diffusion optimization strategies, thereby improving the structural reconstruction accuracy of the final geological 3D exploration image.

[0016] Preferably, the step of calculating the normalization adjustment factor for each direction of the eigenvector based on the eigenvalue includes: […]. Direction down mesh node The eigenvalues ​​divided by the mesh nodes in the three directions The sum of the eigenvalues ​​is the result of the calculation. The relative weight of directions, It is the index of the direction corresponding to the feature vector. The values ​​are 1, 2, and 3, representing the mesh nodes in the three directions. The eigenvalues ​​are respectively Subtract the first value from the value 1. The result obtained from the relative weights of the directions is subjected to maximum value normalization, and the final normalized value is the first value obtained. Directional normalization adjustment factor.

[0017] Preferably, the preset region type factor includes: when grid nodes When it belongs to a layered structure region, the region type factor When grid nodes When it belongs to a fractured structural region, the region type factor When grid nodes When it belongs to a homogeneous structure region, the region type factor .

[0018] This method configures differentiated region type factors for layered, fractured, and homogeneous geological structures, and sets specific factor values ​​based on the optimization needs of different geological structure regions. The factor values ​​for layered geological structures emphasize the smoothing effect of non-dominant directions, the factor values ​​for fractured geological structures weaken the overall smoothing adjustment amplitude, and the factor values ​​for homogeneous geological structures achieve a moderate global smoothing adjustment effect. This differentiated configuration allows the region type factors to accurately match the optimization needs of the corresponding geological structure regions, providing a targeted quantitative basis for the subsequent calculation of smoothing adjustment coefficients. It avoids problems such as distortion of layered structures, blurred fracture boundaries, or insufficient noise reduction of homogeneous structures caused by using a uniform adjustment strategy, thereby improving the adaptability and effectiveness of anisotropic diffusion iterative optimization.

[0019] Preferably, the step of constructing the diffusion matrix of each grid node based on the smoothing adjustment coefficient and eigenvector of each grid node in all directions includes: ;in, It is a grid node The diffusion matrix; It is the first Direction down mesh node The smoothing adjustment coefficient; It is the first Direction down mesh node eigenvectors; This indicates transpose.

[0020] This method constructs a diffusion matrix for grid nodes by fusing smoothing adjustment coefficients and eigenvectors in each direction. The product of the eigenvector and its transpose represents the structural characteristics in the corresponding direction. The weighting effect of the smoothing adjustment coefficients integrates multi-directional characteristics, and the constructed diffusion matrix can accurately correlate with the optimization requirements of the geological structure region to which the grid node belongs. This construction method enables the diffusion matrix to have both directional distinctiveness and the ability to adapt to the differentiated diffusion requirements of layered structure regions, fracture structure regions, and homogeneous structure regions through smoothing adjustment coefficients. In the subsequent anisotropic diffusion iteration optimization process, it can achieve directional weighted diffusion that matches the geological structure orientation, avoid the loss of structural details caused by indiscriminate diffusion operations, strengthen the consistency of bedding orientation in layered structure regions, weaken the excessive smoothing trend in fracture structure regions, and improve the noise reduction effect in homogeneous structure regions. This provides key matrix support for the final generation of geological 3D exploration images with clear structure and accurate attributes.

[0021] Preferably, obtaining the reconstructed three-dimensional geological exploration image includes: iteratively optimizing the initial three-dimensional mesh using anisotropic diffusion technology based on the diffusion matrix of all mesh nodes; setting an iteration stop condition: stopping the iteration when the number of iterations reaches 10, and finally obtaining the final three-dimensional mesh after optimization convergence; finally, inputting the geological attribute values ​​of each mesh node in the final three-dimensional mesh into a three-dimensional graphics rendering engine to map and generate the reconstructed three-dimensional geological exploration image.

[0022] The beneficial effects of this invention are as follows: By introducing a weighted calculation method that combines confidence weight with normalized signal-to-noise ratio and stability coefficient, differentiated contribution weights are assigned to data points with different reliability. Simultaneously, by dynamically adjusting the influence of distance attenuation on contribution weights based on geological complexity and adaptive adjustment parameters, the interference of low-confidence data points on interpolation results is effectively reduced, improving the calculation accuracy of geological attribute values ​​of the grid nodes to be interpolated. By pre-setting a search radius to filter the set of neighboring points, the number of data points involved in the calculation can be reduced while ensuring the effectiveness of interpolation, thus reducing the overall computational load and improving the construction efficiency of the initial 3D grid. By pre-setting a ratio threshold and a difference threshold to quantitatively define layered structure regions, fracture structure regions, and uniform structure regions, accurate identification and delineation of different types of geological structures are achieved. This provides a clear structural basis for subsequent differentiated grid optimization. By calculating relative weights and normalization adjustment factors based on eigenvalues, and combining them with preset region type factors, targeted smoothing adjustment coefficients are obtained. This allows the smoothing adjustment strategy to be adapted to the geological structure type of the grid nodes. Then, by implementing anisotropic diffusion iterative optimization through a diffusion matrix constructed by fusing the smoothing adjustment coefficients and eigenvectors, efficient smoothing and noise reduction can be achieved in homogeneous structure regions. In layered and fractured structure regions, the boundary features and detailed information of the geological structure are effectively preserved, avoiding structural distortion caused by excessive smoothing. The final generated 3D geological exploration image has both clear structural hierarchy and accurate geological attribute distribution, providing highly reliable visualization support for geological exploration analysis and resource development decisions. Attached Figure Description

[0023] Figure 1 This is a flowchart illustrating a geological three-dimensional exploration image reconstruction method according to the present invention;

[0024] Figure 2 This is a schematic illustration of a geological three-dimensional exploration image reconstruction method according to the present invention. Detailed Implementation

[0025] 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, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0026] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0027] This invention discloses a method for reconstructing three-dimensional geological exploration images, referring to... Figure 1 This includes steps S1-S3:

[0028] S1. Collect and preprocess exploration data to obtain raw exploration data; interpolate the raw exploration data to construct an initial three-dimensional mesh.

[0029] It should be noted that this step aims to transform discrete exploration data into a more regular three-dimensional mesh structure through an adaptive interpolation algorithm, thereby providing a foundation for subsequent detail enhancement and feature extraction. Preprocessing discrete exploration data can improve its quality. Based on the assumption of geological continuity, a weighted interpolation function is designed to balance the accuracy of data points in discrete areas and dense areas, ensuring that the initial three-dimensional mesh model can truly reflect the overall situation of the geological body.

[0030] Specifically, exploration data is obtained based on the scope of the exploration area and the exploration needs. For example, exploration data includes seismic wave data, borehole coordinate data, and lithological sample data. The acquisition methods are existing technologies and will not be described in detail here.

[0031] Furthermore, the acquired exploration data undergoes preprocessing, specifically including: data denoising, data alignment, outlier removal, and data standardization; for example, bandpass filtering is used to remove environmental noise from the exploration data; based on 3 σ The principle of statistical methods is used to identify and remove outlier data points in exploration data that significantly deviate from geological patterns; exploration data from different sources and coordinate systems are transformed into a unified Cartesian three-dimensional coordinate system to ensure that all data points are strictly aligned in space; and exploration data of different dimensions are mapped to a unified Cartesian three-dimensional coordinate system using the minimum-maximum normalization method. Intervals were defined to eliminate dimensional differences; all preprocessed exploration data were recorded as the original exploration dataset, where... It is the index of data points in the original exploration dataset.

[0032] It should be noted that the core functions of data denoising, data alignment, outlier removal, and data standardization in preprocessing operations are as follows: Data denoising can filter out invalid signals caused by environmental interference during exploration, improving data purity; Data alignment can unify exploration data from different sources to the same spatial reference, avoiding modeling misalignment caused by coordinate deviations; Outlier removal can remove extreme abnormal data caused by equipment failure and operational errors, reducing their interference with interpolation results; Data standardization improves the fusion and adaptability of different types of exploration data by unifying data dimensions, ultimately ensuring the quality of the original data and laying the foundation for subsequent high-precision interpolation.

[0033] Furthermore, based on the original exploration dataset, an initial 3D mesh covering the exploration area is constructed. Specific operations include: calculating the normalized signal-to-noise ratio (SNR) of all data points in the original exploration dataset using the Fourier transform method. The max-min normalization method can be used, for example, to normalize the data points... The normalized signal-to-noise ratio is denoted as The Fourier transform method is a well-known technique and will not be elaborated upon here. The three-dimensional coordinates of all data points are directly read from the preprocessed original exploration dataset. For example, the data points... The three-dimensional coordinates are as follows Based on the results of physical exploration, the geological attribute values ​​of all data points in the original exploration dataset are obtained.

[0034] It should be noted that for seismic wave data, properties such as porosity and rock density are obtained through wave impedance inversion; for borehole coordinate data, properties such as lithological hardness and mineral composition ratio are obtained through on-site sampling and testing; for lithological sample data, properties such as formation thickness and permeability measured in the laboratory are directly read, and these properties are combined as the geological attribute values ​​of the corresponding data points.

[0035] Specifically, using data points from the original exploration dataset Taking data points as an example for analysis Set as center The three-dimensional neighborhood is set to satisfy the layered distribution characteristics of a three-dimensional geological body; data points are... The normalized value of the ratio of a geological attribute value to the average of the geological attribute values ​​of all data points in its three-dimensional neighborhood is used as the data point's value. The stability coefficient is denoted as Normalization can be performed using the min-max normalization method; all data points in the original exploration dataset are traversed to obtain the stability coefficient of each data point, which is used to reflect the degree of geological consistency of each data point in the local area.

[0036] Furthermore, based on data points Normalized signal-to-noise ratio and stability coefficient Calculate and obtain data points The credibility weights are expressed by the following formula:

[0037] ;

[0038] In the formula, Data points Credibility weight; Data points The normalized signal-to-noise ratio has a weighting factor of 0.6; Data points The stability coefficient and weighting coefficient are 0.4; when data points When the normalized signal-to-noise ratio is higher and the stability coefficient is higher, the data points The greater the credibility weight, the higher the credibility weight, and vice versa.

[0039] It should be noted that setting data points The weight of the normalized signal-to-noise ratio is greater than that of the data points. The logic behind the weighting of the stability coefficient is as follows: the signal-to-noise ratio directly reflects the quality of the data acquisition substrate and is a decisive factor in data reliability; while the stability coefficient only reflects local geometric features. Relatively speaking, the quality of the substrate data contributes more to the overall reconstruction accuracy than local geometric consistency, so the signal-to-noise ratio is given a higher weight.

[0040] Specifically, the ratio of the total number of exploration data points in the exploration area to the total volume of the exploration area is used as the average data point density of the exploration area, and the reciprocal of this average data point density is used as the adaptive adjustment parameter. Its units are consistent with coordinates; based on the geological survey report of the exploration area, the geological complexity of the exploration area is set. Exemplary When the geological conditions are simple, it indicates that the strata in the exploration area are flat and there are no obvious faults or folds. When the time is right, it indicates that there are a few faults and folds in the exploration area, and the geological complexity is moderate. When the time is specified, it indicates that there are numerous faults and folds in the exploration area, and the geological complexity is classified as complex.

[0041] Furthermore, let the mesh nodes to be interpolated in the initial 3D mesh constructed be... The three-dimensional coordinates are ,in, It is the index of the grid node to be interpolated in the initial 3D mesh; the grid node to be interpolated is calculated. The Euclidean distance between each data point in the original exploration dataset and the grid nodes to be interpolated. Centered on a grid, a search radius is set, and all data points whose Euclidean distance is less than the search radius are used as grid nodes to be interpolated. Obtain the grid nodes to be interpolated from the nearest neighboring points. The set of neighboring points; where the search radius can be set according to actual needs, for example, the search radius is set to 5.

[0042] Furthermore, based on the geological continuity assumption, and utilizing the principle that data points that are closer in distance and have higher geological similarity have greater contribution weights, this approach is combined with the data points... Credibility weights, adaptively adjusted parameters and geological complexity Calculate and obtain data points Mesh nodes for interpolation The contribution weights are expressed as follows:

[0043] ;

[0044] In the formula, Data points Mesh nodes for interpolation Contribution weight; Data points With the grid nodes to be interpolated The Euclidean distance between them; This is a distance attenuation parameter, for example. ; It is an adaptive adjustment parameter; It is the degree of geological complexity; Data points Credibility weight; It is a complexity correction factor. The greater the geological complexity, the smaller the effective smoothing radius, and the more the interpolation focuses on the local area to protect the details. It is the index of data points in the original exploration dataset.

[0045] Among them, when data points The higher the credibility weight, the more data points... With the grid nodes to be interpolated The smaller the Euclidean distance between them, the larger the adaptive adjustment parameter; the smaller the geological complexity, the more data points... Mesh nodes for interpolation The greater the contribution weight, the higher the weight, and vice versa.

[0046] It should be noted that traditional interpolation typically assumes all control points have the same confidence level, while this relationship exists in the exponential term. Multiply by credibility weight This allows data points with high signal-to-noise ratio and good local geological stability to dominate the interpolation process, directly amplifying the influence of high-quality data while suppressing the interference of low-quality data on the reconstruction results, thus improving the overall robustness of the model. For the exponential term, a Gaussian kernel function is used as the main structure, utilizing the negative exponential decay characteristic of Euclidean distance and strictly adhering to the spatial autocorrelation assumption in geology. This ensures that the interpolation weights decrease smoothly with increasing spatial distance, reflecting the continuous spatial variation characteristics of geological bodies and avoiding erroneous geometric guidance of local structures from distant irrelevant points. The denominator... A dynamic and efficient smoothing radius was constructed, and the parameters were adaptively adjusted. With average data point density Related, with a relatively low average data point density, adaptive parameter adjustment. Larger values ​​mean a wider search range, ensuring regional connectivity; higher average data point density allows for adaptive parameter adjustment. Smaller, meaning a smaller search area, focusing on local details; when the geological complexity is relatively low. Increasing this value, approaching 1, allows for a larger smoothing radius to fill in blank areas, especially when the geological complexity is high. The radius of smoothing decreases and approaches zero, resulting in a sharp contraction of the effective smoothing radius. This effectively prevents excessive smoothing across faults and solves the core problem of blurred fault boundaries in traditional methods.

[0047] Specifically, iterate through all data points in the original exploration dataset to obtain the grid node for each data point to be interpolated. The contribution weights are based on the grid nodes to be interpolated. The set of neighboring points is used to perform a weighted summation of the geological attribute values ​​of all neighboring data points using contribution weights. The calculated weighted average is then used as the grid node to be interpolated. The geological attribute values ​​are obtained by traversing all the grid nodes to be interpolated and obtaining the geological attribute values ​​of all the grid nodes to be interpolated. Based on this, the initial three-dimensional grid is constructed.

[0048] S2. Based on the initial 3D mesh, obtain the gradient vector of each mesh node, construct the structure tensor matrix of each mesh node, perform eigenvalue decomposition on the structure tensor matrix to obtain the eigenvector and eigenvalue of each mesh node, and then divide the layered structure region, fracture structure region and uniform structure region; obtain the smoothing adjustment coefficient of each mesh node in all directions.

[0049] It should be noted that this step aims to overcome the inherent defects of traditional radial basis function interpolation algorithms, such as fixed kernel functions and lack of incorporation of prior geological knowledge. By extracting geological structure information from the initial 3D mesh and constructing a geological structure-guided tensor field, adaptive guidance for local interpolation strategies is achieved. While maintaining the overall continuity of the geological body, this method can highlight the boundary details of discontinuous structures such as faults, effectively reducing the over-smoothing phenomenon commonly found in traditional methods, thereby significantly improving the characterization accuracy of 3D exploration images in complex structural areas.

[0050] Specifically, based on the initial three-dimensional mesh constructed in step S1, the gradient vector of each mesh node on the geological attribute value is calculated. The specific operation includes: using the three-dimensional Sobel operator to calculate the partial derivative of the geological attribute value of each mesh node in the initial three-dimensional mesh to obtain the gradient vector of each mesh node. The magnitude of the gradient vector can reflect the intensity of the change of the geological attribute value of the mesh node in space. The larger the magnitude, the more significant the geological difference. The direction of the gradient vector is along the direction of the fastest change rate of the geological attribute value. In the geological scene, the direction of the gradient vector is usually perpendicular to the rock layer interface.

[0051] Furthermore, a three-dimensional Gaussian smoothing kernel is used to smooth the gradient vector of each grid node to obtain its structure tensor matrix. This structure tensor matrix is ​​a 3×3 symmetric positive semi-definite matrix, whose eigenvalues ​​and eigenvectors can quantitatively describe the geometric features of the geological structure within the neighborhood of each grid node, such as fault strike. The three-dimensional Gaussian smoothing kernel is existing technology and will not be elaborated upon here. Eigenvalue decomposition is performed on the structure tensor matrix of each grid node to obtain three eigenvalues. The structure tensor matrix of each grid node is then used to perform eigenvalue decomposition. Taking the example of mesh nodes, let's analyze them. The eigenvalues ​​are respectively ,in, Based on the actual needs of geological exploration, preset ratio thresholds, first preset difference thresholds, and second preset difference thresholds are set in advance; when and The ratio is greater than the preset ratio threshold, and and When the difference is less than a first preset difference threshold, the grid node is determined to be located in a layered structure region. At this time, the maximum eigenvalue... Corresponding feature vector Indicates the direction of the normal to the interface of rock strata, while The plane formed is parallel to the rock strata interface; when and The difference is less than the first preset difference threshold, and and When the ratio is greater than a preset ratio threshold, the mesh node is determined to be located in a fractured structural region. At this time, the minimum eigenvalue... Corresponding feature vector Indicates the strike of the fault zone, and The corresponding plane is perpendicular to the fracture surface; when and The difference is less than the second preset difference threshold. and When the difference is less than the second preset difference threshold, the grid node is determined to be located in a uniform structure region. At this time, the geological attribute values ​​change gently in all directions in space and there is no significant directional structure. Traverse all grid nodes to obtain grid nodes located in layered structure region, fracture structure region and uniform structure region.

[0052] It should be noted that the theoretical basis for the above division is that the eigenvalues ​​obtained from the structural tensor matrix directly reflect the spatial distribution of the local gradient field; in layered structural regions, geological attribute values ​​change drastically only in a single direction perpendicular to the bedding plane, causing the gradient vector group to be highly concentrated in this normal direction, thus appearing as the largest eigenvalue in the structural tensor, such as... In fractured regions, geological property values ​​may abruptly change within a two-dimensional plane perpendicular to the fracture surface, resulting in the gradient vector group being distributed within a single plane. This manifests as two similar and relatively large eigenvalues ​​in the structural tensor, such as... In regions with uniform structures, geological attribute values ​​change weakly and without directionality, and the gradient vector has a small magnitude and random direction, resulting in the three eigenvalues ​​being close and small in value. This partitioning method based on gradient statistical characteristics provides a solid mathematical and physical foundation for extracting quantitative and computable a priori knowledge of geological structures from numerical models.

[0053] Specifically, based on the feature vectors of grid nodes in layered, fractured, and homogeneous structural regions, structural correction coefficients are constructed for each region; using the grid nodes... Taking the example of analysis, let's take the first... Direction down mesh node The eigenvalues ​​divided by the mesh nodes in the three directions The sum of the eigenvalues ​​is the result of the calculation. The relative weight of the direction, which reflects the first Direction down mesh node The greater the relative weight of the geological attribute value, the more drastic the change; among them, the directional index... The values ​​are 1, 2, and 3, representing the mesh nodes in the three directions. The eigenvalues ​​are respectively Subtract the first value from the value 1. The result obtained from the relative weights of the directions is subjected to maximum value normalization, and the final normalized value is the first value obtained. The normalization adjustment factor for direction is denoted as ; It is the index of the direction corresponding to the feature vector.

[0054] Furthermore, set the region type factor. When grid nodes When it belongs to a layered structure region, This indicates suppression along the feature vector. Smoothing of direction, increasing along the feature vector Smoothing of direction; when grid nodes When it belongs to the fracture structure region, This indicates suppression along the feature vector. The smoothing of the direction protects and sharpens the boundaries of faults or cracks; when the mesh nodes... When it belongs to a homogeneous structure region, This indicates that the amount of time along the feature vector should be appropriately increased. Smoothing of direction; where 0.5 and 0.1 are exemplary values ​​and can be set according to actual needs; according to the first Directional normalization adjustment factor and region type factor Multiply to obtain the first... Direction down mesh node The smoothing adjustment coefficient is denoted as ; Traverse all grid nodes to obtain the first Smoothing adjustment coefficients for each grid node under the direction.

[0055] It should be noted that this step integrates prior geological knowledge with data-driven features to generate smoothing adjustment coefficients for different directions of each grid node. The preset region type factor sets a basic smoothing strategy for layered structure regions, fractured structure regions, and homogeneous structure regions. The adjustment factor calculated based on local gradient features then adaptively fine-tunes this basic strategy. Its core principle is to suppress smoothing in directions where geological attribute values ​​change drastically and to enhance smoothing in directions where geological attribute values ​​change gently. This effectively solves the technical problems of boundary ambiguity and loss of structural details in traditional radial basis function interpolation methods.

[0056] S3. Based on the smoothing adjustment coefficients of each grid node in all directions, the initial three-dimensional grid is iteratively optimized to obtain the final three-dimensional grid, and the reconstructed geological three-dimensional exploration image is output through visualization mapping.

[0057] It should be noted that the core of this step is to transform the smoothing adjustment coefficients and corresponding eigenvectors of each grid node obtained in S2 into precise grid optimization actions, thereby achieving a balance between noise removal and structural fidelity through directional diffusion. The diffusion matrix constructed from the smoothing adjustment coefficients and eigenvectors is essentially a customized optimization criterion for each grid node. Anisotropic diffusion technology can dynamically adjust the smoothing intensity in different directions based on this matrix, strengthening diffusion in homogeneous structural regions to filter noise, and suppressing diffusion at the boundaries of fractured and layered structural regions to preserve details. After a finite number of iterations and convergence, the quantified geological attribute values ​​of the grid nodes are transformed into a three-dimensional geological exploration image that intuitively reflects the underground geological structure through three-dimensional visualization mapping, providing precise support for geological analysis.

[0058] Specifically, according to the first The smoothing adjustment coefficients and corresponding eigenvectors of each grid node under the direction are used to construct the diffusion matrix of each grid node. Taking this as an example, the relational expression is constructed as follows:

[0059] ;

[0060] in, It is a grid node The diffusion matrix; It is the first Direction down mesh node The smoothing adjustment coefficient; It is the first Direction down mesh node The eigenvectors, i.e., the eigenvectors ; Indicates transpose; It is the index of the direction corresponding to the eigenvector; traverse all grid nodes to obtain the diffusion matrix of all grid nodes.

[0061] Furthermore, based on the diffusion matrix of all grid nodes, the initial 3D mesh is iteratively optimized using the well-known anisotropic diffusion technique. This technique dynamically controls the smoothing intensity and range in each direction during each iteration using the diffusion matrix of each grid node: in directions with a larger diffusion matrix, geological attribute values ​​diffuse more strongly, effectively smoothing and denoising areas of geological continuity; in directions with a smaller diffusion matrix, diffusion is strictly suppressed, ensuring that discontinuous boundaries such as faults remain sharp during iteration and preventing loss of detail. An iteration stop condition is set: iteration stops when the number of iterations reaches 10, where 10 is an example value and can be set according to actual needs. Finally, the optimized and converged final 3D mesh is obtained. The geological attribute values ​​of each grid node in the final 3D mesh are then input into the 3D graphics rendering engine to generate a reconstructed 3D geological exploration image. (See [link to relevant documentation]). Figure 2 , Figure 2 This is a schematic illustration of a geological three-dimensional exploration image reconstruction method according to the present invention.

[0062] Specifically, Figure 2 Different colors correspond to different geological attribute values. For example, dark areas may represent high-density rock layers, such as granite, while light-colored areas may represent low-density loose strata, such as soil. Gradual color transitions reflect continuous changes in geological attributes, and abrupt color boundaries correspond to discontinuities in geological structures, such as faults. Layered structure areas are represented by parallel, continuous color bands with clear boundaries and no blurred transitions, accurately restoring the distribution characteristics of rock layers, such as the bedding structure of sedimentary rocks, demonstrating the optimization effect of suppressing smoothing in the normal direction. Fault structure areas are represented by abrupt color bands that run through multiple color bands, with colors that are significantly different from the layered structures on both sides. The boundaries are sharp and without trailing, clearly showing the direction, depth, and cutting range of the fault zone, solving the problem of blurred fault boundaries in traditional methods. Uniform structure areas are represented by a gently distributed area of ​​a single color with no obvious color fluctuations, demonstrating the effect of uniform diffusion noise reduction, while avoiding the loss of details caused by excessive smoothing.

[0063] Among them, compared with the reconstruction results of the traditional radial basis function interpolation method, Figure 2 Its core advantage lies in avoiding the blurring of details in complex areas while removing noise in uniform areas, thus fully restoring the true spatial distribution characteristics of geological bodies and providing high-precision visualization support for geological decision-making such as underground resource exploration and fault risk assessment.

[0064] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for reconstructing three-dimensional geological exploration images, characterized in that, include: The confidence weight of a data point is obtained by weighted summation of its normalized signal-to-noise ratio and stability coefficient. The average data point density is calculated, and its reciprocal is used as an adaptive adjustment parameter. Based on the geological survey report, the geological complexity is set to calculate the contribution weight of each data point to the interpolation grid node. , ; Data points With the grid nodes to be interpolated The Euclidean distance between them; It is the distance attenuation parameter; It is an adaptive adjustment parameter; It is the degree of geological complexity; It is a natural exponential function; It is the index of data points in the original exploration dataset; based on Euclidean distance, the set of neighboring points is filtered; based on the contribution weight, the geological attribute values ​​of all data points in the set of neighboring points are weighted and summed to obtain the geological attribute values ​​of each grid node to be interpolated, which are used to construct the initial three-dimensional grid; The gradient vectors of each grid node in the initial three-dimensional mesh are obtained by calculating the three-dimensional Sobel operator. The structure tensor matrix is ​​constructed and eigenvalue decomposition is performed to obtain the eigenvalues ​​and eigenvectors of each grid node. The layered structure region, fractured structure region and uniform structure region are divided according to the relative magnitude of the eigenvalues. The normalization adjustment factor corresponding to each direction of the eigenvector is calculated based on the eigenvalue, and multiplied by the preset region type factor to obtain the smoothing adjustment coefficient of each grid node in each direction. Based on the smoothing adjustment coefficients and eigenvectors of each grid node in all directions, a diffusion matrix for each grid node is constructed. Anisotropic diffusion technology is then used to iteratively optimize the initial 3D grid to obtain the final 3D grid. Through 3D visualization mapping, a reconstructed 3D geological exploration image is obtained.

2. The method of claim 1, wherein, The obtained confidence weights for data points include: ； In the formula, Data points Credibility weight; Data points The normalized signal-to-noise ratio has a weighting factor of 0.6; Data points The stability coefficient, with a weighting factor of 0.4; It is the index of data points in the original exploration dataset.

3. The method of claim 1, wherein, The method for obtaining the stability coefficient includes: Based on the results of physical exploration, the geological attribute values ​​of all data points in the original exploration dataset are obtained; using the data points Set as center The three-dimensional neighborhood will be the data points The normalized value of the ratio of a geological attribute value to the average of the geological attribute values ​​of all data points in its three-dimensional neighborhood is used as the data point's value. The stability coefficient.

4. The method of claim 1, wherein, The filtering of the nearest point set based on Euclidean distance includes: Calculate the grid nodes to be interpolated The Euclidean distance between each data point in the original exploration dataset and the grid nodes to be interpolated. Centered on the target data point, all data points whose Euclidean distance is less than the preset search radius are used as grid nodes for interpolation. Obtain the grid nodes to be interpolated from the nearest neighboring points. The set of neighboring points.

5. The method of claim 1, wherein, The method of dividing the structure into layered, fractured, and uniform regions based on the relative magnitude of eigenvalues ​​includes: mesh nodes. The eigenvalues ​​are respectively ,in, Based on the actual needs of geological exploration, preset ratio thresholds, first preset difference thresholds, and second preset difference thresholds are set in advance; when and The ratio is greater than the preset ratio threshold, and and When the difference is less than a first preset difference threshold, the grid node is determined to be located in a layered structure region; when and The difference is less than the first preset difference threshold, and and When the ratio is greater than a preset ratio threshold, the grid node is determined to be located in a fractured structural region; when and The difference is less than the second preset difference threshold. and When the difference is less than the second preset difference threshold, the grid node is determined to be located in a uniform structure region.

6. The method of claim 1, wherein, The normalization adjustment factor for each direction of the eigenvector calculated based on eigenvalues ​​includes: The first Direction down mesh node The eigenvalues ​​divided by the mesh nodes in the three directions The sum of the eigenvalues ​​is the result of the calculation. The relative weight of directions, It is the index of the direction corresponding to the feature vector. The values ​​are 1, 2, and 3, representing the mesh nodes in the three directions. The eigenvalues ​​are respectively Subtract the first value from the value 1. The result obtained from the relative weights of the directions is subjected to maximum value normalization, and the final normalized value is the first value obtained. Directional normalization adjustment factor.

7. The geological three-dimensional exploration image reconstruction method according to claim 1, characterized in that, The preset region type factor includes: When grid nodes When it belongs to a layered structure region, the region type factor When grid nodes When it belongs to a fractured structural region, the region type factor When grid nodes When it belongs to a homogeneous structure region, the region type factor .

8. The method of claim 1, wherein, The process of constructing the diffusion matrix for each grid node based on the smoothing adjustment coefficients and eigenvectors of each grid node from all directions includes: ； in, It is a grid node The diffusion matrix; It is the first Direction down mesh node The smoothing adjustment coefficient; It is the first Direction down mesh node eigenvectors; This indicates transpose.

9. The geological three-dimensional exploration image reconstruction method according to claim 1, characterized in that, The obtained reconstructed three-dimensional geological exploration images include: Based on the diffusion matrix of all grid nodes, the initial 3D grid is iteratively optimized using anisotropic diffusion technology. An iteration stopping condition is set: the iteration stops when the number of iterations reaches 10, and the final 3D grid after optimization convergence is finally obtained. Finally, the geological attribute values ​​of each grid node in the final 3D grid are input into the 3D graphics rendering engine and mapped to generate the reconstructed 3D geological exploration image.

Images