A gravity three-dimensional imaging method and system based on improved radon transform and regularization inversion

Abstract

The application discloses a kind of based on improved Radon transform and regularized inversion's gravity three-dimensional imaging method and system, it is related to geophysical gravity exploration and potential field imaging technical field, this method mainly includes: using improved Radon transform to establish the stable mapping between space domain and Radon domain, in Radon domain, traditional three-dimensional gravity imaging problem is converted into a group of structure consistent two-dimensional regularized inversion sub-problems, and through different depth merging and inverse Radon transform restores space domain three-dimensional density image. Based on the gravity three-dimensional imaging method and system of improved Radon transform and regularized inversion provided in the application, large area large data gravity imaging can be used by parallel computing, and the amount of calculation and memory occupation are reduced.

Description

Technical Field

[0001] This invention relates to the field of geophysical gravity exploration and potential field imaging technology, and more specifically, to a gravity three-dimensional imaging method and system based on improved Radon transform and regularized inversion. Background Technology

[0002] Gravity exploration is an important passive source geophysical exploration method, characterized by its large detection depth, wide coverage, and relatively low observation cost. It has wide applications in regional mineral resource surveys, oil and gas and salt dome exploration, basement undulation and basin structure research, fault structure identification, geothermal and groundwater surveys, engineering geological evaluation, urban underground space exploration, and marine and planetary geophysical research. With the development of high-precision gravity instruments, airborne gravity measurements, and satellite gravity and gravity gradient observation technologies, the scale of gravity data has increased rapidly, and three-dimensional gravity inversion and three-dimensional imaging have gradually become important technical directions in current potential field research and applications.

[0003] Existing gravity 3D inversion methods typically establish a mapping relationship between surface observation data and subsurface 3D density models directly in the spatial domain, forming a massive 3D forward kernel matrix. This is then combined with methods such as Tikhonov regularization, Lp norm constraints, depth weighting, prior information constraints, machine learning, or neural networks to solve for the subsurface density structure. While these methods can obtain relatively complete 3D results, they generally suffer from the following bottlenecks: First, the 3D kernel matrix is ​​extremely large, requiring significant memory, making it difficult for ordinary computers to store directly under large-area, high-resolution grid conditions. Second, the high coupling between different directions and layers in the spatial domain 3D inversion leads to low solution efficiency, making it difficult to fully utilize the advantages of parallel computing. Third, in inversions involving large areas, multiple survey lines, and large sample data, the kernel matrix is ​​frequently called repeatedly, and the computational load increases dramatically with the number of observation points and model units. Fourth, deep anomalies are easily affected by shallow anomalies and regularization smoothing effects in conventional inversions, resulting in insufficient boundary recovery capabilities.

[0004] Existing research has shown that the Radon transform possesses the characteristics of line integration and dimensionality reduction representation, enabling the mapping of two-dimensional spatial domain data to the Radon domain composed of angle and distance parameters. This facilitates orientation analysis, feature separation, and dimensionality reduction solutions to certain inverse problems. Early literature has proposed using the Radon transform to analyze the relationship between three-dimensional and two-dimensional anomalies and explored the idea of ​​using the Radon transform to reduce the dimensionality of three-dimensional problems. However, most existing publicly available solutions remain at the level of theoretical relationships, idea verification, or local processing, lacking a stable, parallelizable, and directly applicable large-scale gravity three-dimensional imaging technology route suitable for engineering implementation.

[0005] On the other hand, the classical Radon transform is essentially an integral transform, and the integration process amplifies low-frequency components, which can easily suppress high-frequency boundary information in numerical processing, leading to anomalous boundary diffusion, local peak distortion, and even the introduction of artifacts. Furthermore, directly performing the inverse transform through the Fourier domain still faces issues of computational time and numerical stability. For large-area and large-data gravity inversion tasks, relying solely on the classical Radon transform is insufficient to simultaneously achieve computational efficiency, boundary resolution, noise robustness, and engineering feasibility.

[0006] Therefore, there is an urgent need to propose a new gravity three-dimensional imaging method that can not only take advantage of the dimensionality reduction and decoupling of Radon transform, but also overcome the problems of low-frequency enhancement and unstable reconstruction of classical Radon transform. Furthermore, it is necessary to develop a system scheme that is suitable for angle task decomposition, unified two-dimensional kernel function reuse, parallel scheduling and large-area big data computing, so as to meet the actual needs of modern gravity exploration and deep imaging. Summary of the Invention

[0007] The purpose of this invention is to provide a gravity three-dimensional imaging method and system based on improved Radon transform and regularized inversion, which can utilize parallel computing to perform gravity imaging of large areas and large datasets, reducing computational load and memory consumption.

[0008] This invention provides a gravity three-dimensional imaging method based on improved Radon transform and regularized inversion, comprising the following steps: S1: Based on the gravity anomaly observation data of the target area, perform gridding preprocessing to obtain preprocessed gravity anomaly observation data, and construct the corresponding observation grid and underground three-dimensional model grid. S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the preprocessed gravity anomaly observation data is subjected to an improved Radon transform. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. S3: Based on the observation grid, the underground 3D model grid, and the data corresponding to each discrete angle in the Radon domain representation, and using the forward modeling relationship of the gravity field in the Radon domain, establish two-dimensional regularized inversion subproblems based on a unified two-dimensional forward modeling kernel function, perform two-dimensional regularized inversion calculations, solve for the minimum value of the two-dimensional inversion objective function containing model constraint terms, and obtain the Radon domain density distribution corresponding to each discrete angle. S4: Merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer. Perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer. Reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer.

[0009] Furthermore, the gridding preprocessing includes at least one of detrending, coordinate unification, regular sampling interpolation, outlier removal, and noise suppression.

[0010] Furthermore, the expression for the classical two-dimensional Radon transform is:

[0011] in, For gravity observation data, For the Dirac function, The horizontal coordinate of the spatial domain, The vertical coordinate of the spatial domain, For Radon domain scale parameters, For Radon domain angle parameters, for The classical Radon transform results; The expression for the improved Radon field definition is:

[0012] in, This represents the derivative convolution kernel of the Hilbert transform. This represents the convolution operator. To improve the results of the Radon transform.

[0013] Furthermore, the expression for the system of linear equations is:

[0014]

[0015] in, For the observed data vector, Let be the Radon field parameter vector to be determined. For Radon domain discretization angle subdivision unit, For discrete mapping operators defined by the parameter functions of the improved Radon transform or inverse transform, the superscript... This indicates transpose.

[0016] Furthermore, the solution method for the linear equation system is the conjugate gradient algorithm or the least squares iterative algorithm.

[0017] Furthermore, the expression for the forward modeling relation of the gravitational field in the Radon domain is:

[0018] in, The Radon domain representation of gravity anomalies, Radon domain representation for different depth layer densities. Here is the forward modeling kernel function for the two-dimensional surface density at depth z. The depth of the three-dimensional model space.

[0019] Furthermore, the expression for the two-dimensional inversion objective function is:

[0020] in, The objective function is a two-dimensional inversion function. For a unified two-dimensional forward kernel function, For the corresponding two-dimensional profile density parameters, For the first k Radon domain data from various angles, For regularization parameters, These are model constraint terms.

[0021] Furthermore, the model constraints include at least one of smoothing constraints, sparsity constraints, depth-weighted constraints, and total variation constraints.

[0022] Furthermore, the expression for the inverse Radon transform is:

[0023] in, To improve the results of the Radon inverse transform, To improve the results of the Radon transform.

[0024] This invention also provides a gravity three-dimensional imaging system based on improved Radon transform and regularized inversion, the system comprising the following modules: The gravity data entry and storage module is used for S1: based on the gravity anomaly observation data of the target area, to perform gridded preprocessing to obtain preprocessed gravity anomaly observation data, and to construct the corresponding observation grid and underground three-dimensional model grid; The Radon forward and inverse transformation module is used in S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the improved Radon transform is applied to the preprocessed gravity anomaly observation data. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. The two-dimensional parallel forward modeling and regularized inversion module is used in S3: Based on the observation grid, the underground three-dimensional model grid and the data corresponding to each discrete angle in the Radon domain representation, the two-dimensional regularized inversion subproblems are established based on the forward modeling relationship of the gravity field in the Radon domain and a unified two-dimensional forward modeling kernel function. The two-dimensional regularized inversion calculation is performed to solve the minimum value of the two-dimensional inversion objective function containing the model constraint term and obtain the Radon domain density distribution corresponding to each discrete angle. The three-dimensional density imaging module is used in S4 to: merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer; perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer; and reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer.

[0025] The gravity three-dimensional imaging method and system based on improved Radon transform and regularized inversion provided by this invention have the following beneficial effects: This invention addresses the problems of existing spatial domain 3D gravity inversion methods, such as large kernel matrix size, high memory consumption, large computational load, strong coupling between different directions and depths, difficulty in achieving efficient parallel computing, unsuitability for large-area and large-data inversion, and boundary blurring and local anomaly distortion caused by low-frequency enhancement of the classical Radon integral transform. It utilizes an improved Radon transform to establish a stable mapping between the spatial domain and the Radon domain, transforming the traditional 3D gravity imaging problem into a set of structurally consistent 2D regularized inversion subproblems within the Radon domain, and recovering the spatial domain 3D density image through merging at different depths and inverse Radon transform. The method includes: acquiring gravity anomaly observation data of the study area and establishing a subsurface three-dimensional model mesh; constructing an improved Radon forward and inverse transform and establishing a discrete linear mapping relationship from the spatial domain to the Radon domain, and obtaining the Radon domain representation of the observation data by solving a system of linear equations; performing two-dimensional gravity regularization inversion calculations on the data of each discrete angle in the Radon domain, and using a unified two-dimensional kernel function to obtain the two-dimensional density distribution of the Radon domain corresponding to each discrete angle; merging the Radon domain density results of each depth layer and performing an inverse Radon transform to obtain the spatial domain surface density distribution of each depth layer; and reconstructing the three-dimensional density imaging results of the target area based on the surface density distribution of all depth layers.

[0026] This application fully considers the practical requirements of computational efficiency and storage capacity for gravity 3D imaging under conditions of large-area, large-data observation and deep exploration. By transforming the spatial domain 3D problem into multiple 2D subproblems within the Radon domain, the size of the kernel matrix is ​​reduced and redundant computations are minimized. The inherent independence and parallel decomposition of the 2D inversion subproblems from different perspectives are utilized. Through parallel task scheduling, unified 2D kernel function reuse, and hierarchical result aggregation, the overall computational speed is improved and the engineering feasibility for large-scale applications is enhanced.

[0027] This application fully considers the differentiated response characteristics of gravity anomalies in different directions and depths. By introducing angle and distance parameters to characterize the anomaly distribution, information from different directions can be separated and focused in the Radon domain, providing a foundation for directional solutions and result reconstruction. By constructing an improved Radon transform and inverse transform, the application solves the numerical defects caused by the classical Radon transform integration process, which leads to low-frequency enhancement and high-frequency boundary suppression. This enhances the ability to preserve anomaly boundary and local mutation information, thereby improving reconstruction stability and imaging resolution.

[0028] This application fully considers the practical interpretation needs of gravity 3D imaging results in depth layer reconstruction, spatial domain recovery, and anomaly boundary characterization. By performing inverse transformation and 3D combination of Radon domain density results at each depth layer, spatial density imaging results that are easy to interpret geologically are obtained. This method is not only applicable to model experiments but also to large-area gravity 3D imaging tasks under complex geological conditions.

[0029] Compared with existing technologies, this invention improves the preservation of anomalous boundaries, high-frequency information, and local peaks by constructing an improved Radon transform, thus alleviating the low-frequency enhancement problem caused by classical Radon integrals. By decoupling the 3D gravity imaging problem into multiple 2D regularized inversion subproblems in the Radon domain according to angles, the difficulty of solving 3D coupled problems is significantly reduced. Through the independence between tasks at different angles, it is naturally suitable for multi-threaded, multi-process, GPU, or distributed parallel computing, solving the problem that traditional methods are difficult to achieve efficient parallel computing. By unifying the 2D kernel function construction and caching mechanism, the number of times the kernel matrix is ​​repeatedly constructed is greatly reduced, lowering memory usage and total computational load, making it particularly suitable for large-area, large-grid, and large-data inversion. Through depth-by-depth layer merging and inverse Radon transform, a clear 3D density image in the spatial domain can be obtained, enhancing the ability to identify deep anomalies.

[0030] This invention can transform the traditional spatial domain three-dimensional gravity imaging problem into multiple two-dimensional regularized inversion problems with consistent structure. It has the advantages of easy parallel computing, low computational cost, small memory footprint, and suitability for large-area and large-data gravity imaging. Attached Figure Description

[0031] The present invention will be further described below with reference to the accompanying drawings and embodiments. In the accompanying drawings: Figure 1 This is a flowchart of the gravity three-dimensional imaging method based on improved Radon transform and regularized inversion provided by the present invention; Figure 2 This is a schematic diagram of the spatial domain model provided by the present invention; Figure 3 This is a plane diagram of gravity forward modeling data corresponding to the spatial domain model provided by this invention; Figure 4 This invention provides Figure 2 The Radon transform result of the gravity data is shown in the figure. Figure 5 This is the perspective provided by the present invention. The Radon domain two-dimensional gravity regularization inversion result at 0°; Figure 6 This is the perspective provided by the present invention. The Radon domain two-dimensional gravity regularization inversion result at 150°; Figure 7 This is the first layer depth result image after merging the Radon domain two-dimensional inversion results from all angles provided by this invention; Figure 8 This is the 20th layer depth result image after merging the Radon domain two-dimensional inversion results from all angles provided by this invention; Figure 9 This invention provides a method for... Figure 7 The spatial domain density distribution map of the first layer after performing an inverse Radon transform on the first layer merging results is shown. Figure 10 This invention provides a method for... Figure 8 The spatial domain density distribution map of the 20th layer obtained after performing an inverse Radon transform on the merged results of the 20th layer is shown. Figure 11 This is a diagram showing the three-dimensional inversion density distribution results within the spatial domain provided by this invention; Figure 12 This is a block diagram of the gravity three-dimensional imaging system based on improved Radon transform and regularized inversion provided by the present invention. Detailed Implementation

[0032] To provide a clearer understanding of the technical features, objectives, and effects of the present invention, specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0033] Figure 1A schematic diagram of the gravity three-dimensional imaging method based on improved Radon transform and regularized inversion of this embodiment is shown. In this embodiment, the gravity three-dimensional imaging method based on improved Radon transform and regularized inversion includes the following steps: S1: Based on the gravity anomaly observation data of the target area, perform gridding preprocessing to obtain preprocessed gravity anomaly observation data, and construct the corresponding observation grid and underground three-dimensional model grid. In one exemplary embodiment, the gridding preprocessing includes at least one of detrending, coordinate unification, regular sampling interpolation, outlier removal, and noise suppression; S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the preprocessed gravity anomaly observation data is subjected to an improved Radon transform. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. In one exemplary embodiment, the expression for the classical two-dimensional Radon transform is:

[0034] in, For gravity observation data, For the Dirac function, The horizontal coordinate of the spatial domain, The vertical coordinate of the spatial domain, For Radon domain scale parameters, For Radon domain angle parameters, for The classical Radon transform results; The expression for the improved Radon field definition is:

[0035] in, This represents the Hilbert transform derivative convolution kernel. This represents the convolution operator. To improve the results of the Radon transform; In one exemplary embodiment, the expression for the system of linear equations is:

[0036]

[0037] in, For the observed data vector, Let be the Radon field parameter vector to be determined. For Radon domain discretization angle subdivision unit, For discrete mapping operators defined by the parameter functions of the improved Radon transform or inverse transform, the superscript... This indicates transpose.

[0038] In one exemplary embodiment, the solution method for the linear equation system is the conjugate gradient algorithm or the least squares iterative algorithm; S3: Based on the observation grid, the underground 3D model grid, and the data corresponding to each discrete angle in the Radon domain representation, and using the forward modeling relationship of the gravity field in the Radon domain, establish two-dimensional regularized inversion subproblems based on a unified two-dimensional forward modeling kernel function, perform two-dimensional regularized inversion calculations, solve for the minimum value of the two-dimensional inversion objective function containing model constraint terms, and obtain the Radon domain density distribution corresponding to each discrete angle. In one exemplary embodiment, the expression for the forward modeling of the gravitational field in the Radon domain is:

[0039] in, The Radon domain representation of gravity anomalies. Radon domain representation for different depth layer densities. Here is the forward modeling kernel function for the two-dimensional surface density at depth z. The depth of the three-dimensional model space; In one exemplary embodiment, the expression for the two-dimensional inversion objective function is:

[0040] in, The objective function is a two-dimensional inversion function. For a unified two-dimensional forward kernel function, For the corresponding two-dimensional density profile parameters, For the first k Radon domain data from various angles, For regularization parameters, These are model constraint terms; In one exemplary embodiment, the model constraints include at least one of smoothing constraints, sparsity constraints, depth-weighted constraints, and total variation constraints. As an exemplary embodiment, the two-dimensional inversion tasks of different discrete angles in step S3 are independent of each other and are respectively assigned to different computing nodes, central processing unit threads, graphics processing unit threads or heterogeneous computing units for execution.

[0041] As an exemplary embodiment, step S3 also includes a unified two-dimensional kernel function cache, angle task queue construction, node load balancing, angle result recycling, and inter-layer merging to reduce the overhead of repeated kernel construction, data migration, and parallel scheduling.

[0042] S4: Merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer. Perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer. Reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer.

[0043] In one exemplary embodiment, the expression for the inverse Radon transform is:

[0044] in, To improve the results of the Radon inverse transform, To improve the results of the Radon transform.

[0045] This embodiment provides a gravity three-dimensional imaging system based on improved Radon transform and regularized inversion. The system includes the following modules: The gravity data entry and storage module is used for S1: based on the gravity anomaly observation data of the target area, to perform gridded preprocessing to obtain preprocessed gravity anomaly observation data, and to construct the corresponding observation grid and underground three-dimensional model grid; The Radon forward and inverse transformation module is used in S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the improved Radon transform is applied to the preprocessed gravity anomaly observation data. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. The two-dimensional parallel forward modeling and regularized inversion module is used in S3: Based on the observation grid, the underground three-dimensional model grid and the data corresponding to each discrete angle in the Radon domain representation, the two-dimensional regularized inversion subproblems are established based on the forward modeling relationship of the gravity field in the Radon domain and a unified two-dimensional forward modeling kernel function. The two-dimensional regularized inversion calculation is performed to solve the minimum value of the two-dimensional inversion objective function containing the model constraint term and obtain the Radon domain density distribution corresponding to each discrete angle. The three-dimensional density imaging module is used in S4 to: merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer; perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer; and reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer.

[0046] In some embodiments, the above-described gravity three-dimensional imaging system based on improved Radon transform and regularized inversion can also be implemented in the following ways.

[0047] In this embodiment, the gravity three-dimensional imaging system based on improved Radon transform and regularized inversion includes: The gravity data entry and storage module is used to receive and store gravity anomaly observation data, as well as observation grid and model grid information in the study area; The Radon forward and inverse transform module is used to establish a discrete mapping relationship between the spatial domain and the Radon domain based on the improved Radon transform with Hilbert transform derivative convolution processing, and to obtain the Radon domain representation of the observation data. At the same time, it performs the inverse Radon transform on the merged results of each depth layer. The two-dimensional parallel forward modeling and regularized inversion module is used to establish a two-dimensional regularized inversion subproblem based on a unified two-dimensional gravity forward modeling kernel function for Radon domain data from various discrete angles, and to perform inversion calculations through angle task splitting and parallel scheduling. The 3D density imaging module is used to merge the Radon domain density distributions of all discrete angles by depth layer, and reconstruct the 3D density imaging results based on the spatial domain surface density distribution after inverse Radon transform.

[0048] In one exemplary embodiment, the two-dimensional parallel forward modeling and regularized inversion module includes a unified two-dimensional kernel function construction and caching unit, an angle task splitting unit, a parallel scheduling unit, a two-dimensional regularized inversion unit, and a result collection unit.

[0049] In one exemplary embodiment, the parallel scheduling unit employs multi-threaded parallelism, multi-process parallelism, GPU parallelism, MPI distributed parallelism, or any combination thereof.

[0050] In some embodiments, the above-described gravity three-dimensional imaging method based on improved Radon transform and regularized inversion can also be implemented in the following ways.

[0051] In this embodiment, the gravity three-dimensional imaging method based on improved Radon transform and regularized inversion includes the following steps: Step 1: Acquire gravity anomaly observation data of the target area. Perform one or more of the following preprocessing techniques on the data: detrending, coordinate unification, regular sampling interpolation, outlier removal, and noise suppression. Establish an observation grid and a subsurface 3D model grid based on the observation point spacing, survey area range, and target imaging depth. The gravity anomaly observation data can be Bouguer gravity anomaly, residual gravity anomaly, or target anomaly data after regional field separation.

[0052] Step 2: Based on the classical two-dimensional Radon transform and the definition of the improved Radon domain, perform the improved Radon transform on the preprocessed gravity anomaly observation data, and establish the connection between the spatial domain and the Radon domain through discrete mapping relationship.

[0053] The classical two-dimensional Radon transform satisfies:

[0054] The improved Radon field is defined as follows:

[0055] The improved Radon inverse transform satisfies:

[0056] The discrete inverse transform is a linear transform, therefore it can be written in the form of a system of linear equations:

[0057] The improved Radon transform can be formed into an implicit functional form:

[0058] in: For gravity observation data, This is a vector of gravity observation data. For the Dirac function, for The classical Radon transform results, This represents the Hilbert transform derivative convolution kernel. This represents the convolution operator. The result of this improved Radon transform, For Radon domain discretization angle subdivision unit, For discrete mapping operators defined by the parameter functions of the improved Radon transform or inverse transform. The operator exists implicitly and the kernel matrix cannot be obtained directly. However, the improved Radon transform can be obtained by iteratively solving the linear equation system using the conjugate gradient algorithm, the least squares iterative algorithm, or other iterative solution algorithms. The value of . By introducing the improved Radon transform in the form of Hilbert transform derivative convolution in step 2, the high-frequency components can be enhanced and the low-frequency enhancement effect caused by the classical line integral process can be offset.

[0059] Step 3: For the data corresponding to each discrete angle in the Radon domain representation, establish two-dimensional regularized inverse subproblems based on a unified two-dimensional forward kernel function.

[0060] The three-dimensional gravitational potential satisfies the Poisson equation:

[0061] in, It is the Laplace operator. For a point in three-dimensional space Gravity potential, For a point in three-dimensional space density, is the gravitational constant.

[0062] Along any direction Integral, vertical The plane of direction is set with two perpendicular coordinates, respectively and We can obtain the two-dimensional equivalence relation:

[0063] in, It is a second derivative operator. For vertical A point on the plane of direction Gravity potential, For a certain point areal density, is the gravitational constant.

[0064] The above relationships show that the three-dimensional gravitational potential and its density distribution, after line integration along any direction, satisfy an equivalent Poisson equation relationship with the two-dimensional case, thus providing a theoretical basis for reducing the three-dimensional imaging problem to a two-dimensional inversion in the Radon domain. The forward modeling relationship of the gravitational field in the Radon domain satisfies:

[0065] in: The Radon domain representation of gravity anomalies. Radon domain representation for different depth layer densities. This is the corresponding two-dimensional surface density forward modeling kernel function.

[0066] Establish a two-dimensional inversion objective function for the k-th angle:

[0067] in: For the first k Radon domain data from various angles, For the corresponding two-dimensional density profile parameters, For a unified two-dimensional forward kernel function, For regularization parameters, These are model constraint terms; the model constraint terms may include at least one of smoothness constraint, sparsity constraint, depth-weighted constraint, and total variation constraint.

[0068] Step 4: Merge the Radon domain density distributions of all discrete angles according to depth layers to form the Radon domain merged density map of layers 1 to L, and perform inverse Radon transform on the merged results of each depth layer to obtain the spatial domain surface density distribution of each depth layer.

[0069] Step 5: Reconstruct the 3D density imaging results of the target region based on the spatial domain surface density distribution of each depth layer, and output slice maps, isosurface maps, or volume rendering results. The 3D density structure of the spatial domain can be recovered through a technical chain of "improved Radon transform—2D regularized inversion—layer merging—inverse Radon transform reconstruction".

[0070] In some embodiments, the above-described gravity three-dimensional imaging system based on improved Radon transform and regularized inversion can also be implemented in the following ways.

[0071] In this embodiment, the gravity 3D imaging system based on improved Radon transform and regularized inversion includes a gravity data input and storage module, a Radon forward and inverse transform module, a 2D parallel forward modeling and regularized inversion module, and a 3D density imaging module. The gravity data input and storage module receives, stores, and manages gravity anomaly observation data of the target area, as well as observation grid and model grid information. The Radon forward and inverse transform module performs the improved Radon transform, establishes a discrete mapping relationship between the spatial domain and the Radon domain, and performs an inverse Radon transform on the depth layer aggregation results. The 2D parallel forward modeling and regularized inversion module establishes a 2D regularized inversion subproblem for Radon domain data at each discrete angle and performs parallel solution; it may also include a unified 2D kernel function construction and caching unit, an angle task splitting unit, a parallel scheduling unit, a 2D regularized inversion unit, and a result retrieval unit. The 3D density imaging module aggregates the results at each angle according to the depth layer and reconstructs the 3D density imaging results based on the spatial domain surface density distribution after the inverse Radon transform.

[0072] In some embodiments, the above-described gravity three-dimensional imaging method based on improved Radon transform and regularized inversion can also be implemented in the following ways.

[0073] In this embodiment, the gravity 3D imaging method based on improved Radon transform and regularized inversion includes: establishing a model, calculating spatial domain gravity forward modeling data, implementing the improved Radon transform, establishing 2D regularized inversion sub-problems by angle, merging all angle inversion results by depth, implementing the inverse Radon transform, and recovering the spatial domain 3D density distribution. Please refer to... Figure 1 , Figure 1 This is a flowchart illustrating the method of the present invention.

[0074] Specifically, the following steps are included: Step 1: Obtain gravity anomaly observation data of the study area, and perform one or more gridding preprocessing on the gravity anomaly observation data, including detrending, coordinate unification, regular sampling interpolation, outlier removal and noise suppression. Establish corresponding observation grids and underground 3D model grids according to the distance between observation points, the range of the survey area and the target imaging depth. It should be noted that the gravity anomaly observation data in step one can be forward modeling data from model experiments or ground or airborne gravity anomaly observation data collected in the actual study area.

[0075] Gravity anomaly observation data can be Bouguer gravity anomaly, residual gravity anomaly, or target anomaly data separated by regional field.

[0076] like Figure 2 As shown, a three-dimensional subsurface model mesh is established, and the anomaly model is designed to include a plate-shaped distributed density volume and a locally concentrated density volume, which are spatially adjacent to each other. The density range is given by a color scale. This combined model contains both laterally continuous structures and locally concentrated high-value anomalies, and can be used to verify the imaging capability of this invention when simultaneously processing laterally continuous structures and locally concentrated anomalies.

[0077] Based on the established underground 3D model mesh and the designed anomaly morphology and density values, ground gravity observation data are obtained using the gravity forward modeling analytical formula of the upright prism. Figure 3 Given Figure 2 The model's gravity forward modeling data on the observation surface shows that the observed anomalies appear as two coupled peak regions on the plane, exhibiting obvious lateral diffusion and superposition characteristics. This indicates that it is difficult to directly recover the true three-dimensional density morphology relying solely on two-dimensional observation data in the spatial domain.

[0078] Step 2: Construct an improved Radon transform based on the classical two-dimensional Radon transform. This is achieved by introducing Hilbert transform derivative convolution processing on the Radon domain parameter functions to suppress the low-frequency enhancement effect during integration. Based on the improved Radon transform, a discrete mapping relationship between the spatial domain and the Radon domain is established. The Radon domain representation of the gravity anomaly observation data is obtained by solving a system of linear equations.

[0079] The classical two-dimensional Radon transform satisfies:

[0080] The improved Radon field is defined as follows:

[0081] The improved Radon inverse transform satisfies:

[0082] The discrete inverse transform is a linear transform, therefore it can be written in the form of a system of linear equations:

[0083] The improved Radon transform can be formed into an implicit functional form:

[0084] in: For gravity observation data, This is a vector of gravity observation data. For the Dirac function, for The classical Radon transform results, This represents the Hilbert transform derivative convolution kernel. This represents the convolution operator. The result of this improved Radon transform, For Radon domain discretization angle subdivision unit, For discrete mapping operators defined by the parameter functions of the improved Radon transform or inverse transform. The operator exists implicitly and the kernel matrix cannot be obtained directly. However, the improved Radon transform can be obtained by iteratively solving the linear equation system using the conjugate gradient algorithm, the least squares iterative algorithm, or other iterative solution algorithms. The value of . By introducing the improved Radon transform in the form of Hilbert transform derivative convolution in step S2, the high-frequency components can be enhanced and the low-frequency enhancement effect caused by the classical line integral process can be offset.

[0085] Figure 4 The following is given Figure 3The results of transforming the gravity data to the Radon domain are shown. It can be seen that the gravity anomaly forms a continuous energy distribution band at different angles α and distance parameters τ. The components of different directions are re-encoded in the Radon domain, and two-dimensional regularized inversion subproblems can be established at each discrete angle.

[0086] Step 3: For the data corresponding to each discrete angle in the Radon domain representation, establish two-dimensional regularized inversion subproblems based on a unified two-dimensional forward kernel function, perform two-dimensional regularized inversion calculations, and obtain the Radon domain density distribution corresponding to each discrete angle.

[0087] The forward relation of the gravitational field in the Radon domain satisfies:

[0088] in: The Radon domain representation of gravity anomalies. Radon domain representation for different depth layer densities. This is the forward modeling kernel function for the two-dimensional surface density at depth z.

[0089] Establish a two-dimensional inversion objective function for the k-th angle:

[0090] in: For the first k Radon domain data from various angles, For the corresponding two-dimensional density profile parameters, For a unified two-dimensional forward kernel function, For regularization parameters, These are model constraint terms; the model constraint terms may include at least one of smoothness constraint, sparsity constraint, depth-weighted constraint, and total variation constraint.

[0091] The conjugate gradient algorithm is used to solve for the minimum value of the two-dimensional inversion objective function, and the Radon domain two-dimensional inversion results corresponding to each angle are obtained.

[0092] The inherent independence and parallel decomposition of computational characteristics between different two-dimensional inversion subproblems, through parallel task scheduling and unified two-dimensional kernel function reuse, will significantly improve the inversion calculation speed and greatly save computer memory consumption.

[0093] Figure 5 The results of the two-dimensional gravity-regularized inversion in the Radon domain under the condition of angle α=0° are presented. The upper curve represents the fitting situation at this angle, and the lower image represents the density distribution of the inverted two-dimensional Radon domain. It can be seen that the main anomaly is concentrated around τ≈0 and gradually diffuses with increasing depth, reflecting the projection information of the dominant structure of the model at this angle. Figure 6 Two-dimensional gravity regularized inversion results in the Radon domain under the condition of angle α=150° are presented. Figure 5 In contrast, the anomaly center shifts at this angle, and the morphology of the two-dimensional boundaries changes significantly, indicating that the two-dimensional inversions at different angles respectively characterize the projection features of the model in the corresponding directions. The red line represents the pseudo-Radon domain profile data of the original gravity data, and the blue line represents the fitted data of the pseudo-Radon domain two-dimensional inversion.

[0094] Step 4: Merge the Radon domain density distributions of all discrete angles according to depth layers to form the Radon domain merged density map of layers 1 to L. Perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer.

[0095] Figure 7 The results of the first layer, obtained by merging the two-dimensional inversion results of all discrete angles according to depth layers, are presented. The energy distribution of this layer is relatively complex, preserving both the shallow main anomaly and the peripheral response formed by the superposition of multiple angles. It is an important input for the subsequent inverse Radon reconstruction of the shallow spatial domain.

[0096] Figure 8 The results for layer 20, after merging the angle inversion results, are presented. Figure 7 In comparison, this layer mainly retains the stable anomaly response in the middle and deep layers, and the high-frequency interference in the periphery is significantly reduced, indicating that the deep layer results have more stable structural indicative significance after multi-angle aggregation.

[0097] Figure 9 The image shows the spatial density distribution of the first layer after applying the inverse Radon transform to the merged Radon domain of the first layer. It can be seen that the shallow image already shows two anomalous centers in the spatial domain, but still exhibits some tailing and boundary diffusion, consistent with the characteristic of shallow aggregation results containing a significant amount of angular mixing information.

[0098] Figure 10 The spatial domain density distribution of layer 20 is obtained after applying the inverse Radon transform to the merged results of layer 20. The result shows a relatively regular elliptical high-value region, indicating that the overall outline of the anomalous body is more stable and clear after multi-angle aggregation and inverse Radon reconstruction at deeper layers.

[0099] Step 5: Reconstruct the 3D density imaging results of the study area based on the spatial domain surface density distribution at each depth layer. Output slice maps, isosurface maps, or volume rendering results. Through the technical link of "improved Radon transform—2D regularized inversion—layer merging—inverse Radon transform reconstruction", the 3D density structure of the spatial domain can be recovered.

[0100] Figure 11This is a diagram showing the three-dimensional inversion density distribution results in the spatial domain. From an overall morphological perspective, the obtained three-dimensional imaging results recover the spatial distribution and distribution trend of the main body of the model quite well, verifying that the present invention can recover the three-dimensional density structure in the spatial domain through the technical link of "improved Radon transform - angle decoupling two-dimensional inversion - inter-layer merging - inverse Radon reconstruction".

[0101] Figure 12 The diagram illustrates the modular composition of the system of this invention. As shown in the figure, the system adopts a workflow organization of "gravity anomaly data - Radon forward and inverse transform - Radon domain two-dimensional parallel forward and inverse calculation - three-dimensional density imaging". Among them, the tasks of each angle in the two-dimensional parallel forward and regularized inverse modules can be executed independently and concurrently, which is an important foundation for the efficient parallel computing of this invention.

[0102] In the system implementation, the gravity data input and storage module is used to receive measured gravity anomaly data, perform coordinate unification, outlier removal, regular interpolation and noise suppression on the raw data, and establish observation grid and model grid.

[0103] The Radon forward and inverse transform module establishes discrete mapping relationships based on the improved Radon transform, and expresses the relationship between the observed data and the Radon domain parameters as a system of linear equations. The module solves the Radon domain representation at each discrete angle using the conjugate gradient algorithm. It also performs the inverse Radon transform on the aggregation results of each depth layer to recover the spatial domain surface density distribution.

[0104] The two-dimensional parallel forward modeling and regularized inversion modules establish a two-dimensional surface density forward modeling kernel function based on the model mesh scale and depth layer settings, and cache it for reuse in two-dimensional inversion tasks at various angles; simultaneously, according to the discrete angle set { , ,…, A task queue is generated, and the two-dimensional inversion subproblem corresponding to each angle is allocated to different computing resources. Since each angle's subproblem is only related to the Radon domain data and the unified two-dimensional kernel function of the corresponding angle, there is no strong coupling across angles, making it suitable for multi-threaded, multi-process, GPU, and MPI distributed parallelism.

[0105] The 3D density imaging module sorts the 2D results returned by all computation nodes by angle and reorganizes them by depth layer to form Radon domain aggregated density maps from layer 1 to layer L. Then, it performs an inverse Radon transform on the aggregated results of each depth layer to restore the spatial domain surface density distribution. Finally, it superimposes the results of each depth layer to generate a 3D density volume and outputs slice maps, isosurface maps, or volume plots.

[0106] Without departing from the basic idea of ​​this invention, the model constraint terms in two-dimensional regularized inversion Smoothing constraints, sparsity constraints, depth-weighted constraints, total variation constraints, or combinations thereof can be flexibly selected; this invention can also be extended to three-dimensional imaging of magnetic anomalies, combined gravity and magnetic imaging, combined gravity and gravity gradient imaging, and imaging scenarios constrained by multi-source prior information.

[0107] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the present invention without departing from the spirit and scope of the claims. All of these forms are within the protection scope of the present invention.

Claims

1. A gravity three-dimensional imaging method based on improved Radon transform and regularized inversion, characterized in that, Includes the following steps: S1: Based on the gravity anomaly observation data of the target area, perform gridding preprocessing to obtain preprocessed gravity anomaly observation data, and construct the corresponding observation grid and underground three-dimensional model grid. S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the preprocessed gravity anomaly observation data is subjected to an improved Radon transform. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. S3: Based on the observation grid, the underground 3D model grid, and the data corresponding to each discrete angle in the Radon domain representation, and utilizing the forward modeling relationship of the gravity field in the Radon domain, establish two-dimensional regularized inversion subproblems based on a unified two-dimensional forward modeling kernel function, perform two-dimensional regularized inversion calculations, solve for the minimum value of the two-dimensional inversion objective function containing model constraint terms, obtain the Radon domain density distribution corresponding to each discrete angle, and execute the tasks for each angle concurrently. S4: Merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer, perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer, and reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer. The expression for the classical two-dimensional Radon transform is: in, For gravity observation data, For the Dirac function, The horizontal coordinate of the spatial domain, The vertical coordinate of the spatial domain, For Radon domain scale parameters, For Radon domain angle parameters, for The classical Radon transform results; The expression for the improved Radon field definition is: in, This represents the Hilbert transform derivative convolution kernel. This represents the convolution operator. To improve the results of the Radon transform.

2. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The gridding preprocessing includes at least one of the following: detrending, coordinate unification, regular sampling interpolation, outlier removal, and noise suppression.

3. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The expression for the system of linear equations is: in, For the observed data vector, Let be the Radon field parameter vector to be determined. For Radon domain discretization angle subdivision unit, For discrete mapping operators defined by the parameter functions of the improved Radon transform or inverse transform, the superscript... This indicates transpose.

4. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The solution method for the linear equation system is either the conjugate gradient algorithm or the least squares iterative algorithm.

5. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The expression for the forward modeling relation of the gravitational field in the Radon domain is as follows: in, The Radon domain representation of gravity anomalies. Radon domain representation for different depth layer densities. Here is the forward modeling kernel function for the two-dimensional surface density at depth z. The depth of the three-dimensional model space.

6. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The expression for the two-dimensional inversion objective function is: in, The objective function is a two-dimensional inversion function. For a unified two-dimensional forward kernel function, For the corresponding two-dimensional profile density parameters, k For serial number, For the first k Radon domain data from various angles, For regularization parameters, These are model constraint terms.

7. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The model constraints include at least one of the following: smoothing constraints, sparsity constraints, depth-weighted constraints, and total variation constraints.

8. The gravity three-dimensional imaging method based on improved Radon transform and regularized inversion according to claim 1, characterized in that, The expression for the inverse Radon transform is: in, To improve the results of the Radon inverse transform, To improve the results of the Radon transform.

9. A gravity three-dimensional imaging system based on improved Radon transform and regularized inversion, characterized in that, The system includes the following modules: The gravity data entry and storage module is used for S1: based on the gravity anomaly observation data of the target area, to perform gridded preprocessing to obtain preprocessed gravity anomaly observation data, and to construct the corresponding observation grid and underground three-dimensional model grid; The Radon forward and inverse transformation module is used in S2: Based on the classical two-dimensional Radon transform and the improved Radon domain definition, the improved Radon transform is applied to the preprocessed gravity anomaly observation data. Using discrete mapping relationships, a set of linear equations is established to establish the discrete mapping relationship between the gravity anomaly observation data in the spatial domain and the Radon domain parameters. Solving the set of linear equations yields the Radon domain representation of the gravity anomaly observation data. The two-dimensional parallel forward modeling and regularized inversion module is used in S3: Based on the observation grid, the underground three-dimensional model grid and the data corresponding to each discrete angle in the Radon domain representation, the two-dimensional forward modeling relationship of the gravity field in the Radon domain is used to establish two-dimensional regularized inversion subproblems based on a unified two-dimensional forward modeling kernel function, perform two-dimensional regularized inversion calculations, solve for the minimum value of the two-dimensional inversion objective function containing model constraint terms, obtain the Radon domain density distribution corresponding to each discrete angle, and execute the tasks of each angle concurrently; The three-dimensional density imaging module is used in S4 to: merge the Radon domain density distributions corresponding to each discrete angle by depth layer to obtain the aggregation results of each depth layer; perform inverse Radon transform on the aggregation results of each depth layer to obtain the spatial domain surface density distribution of each depth layer; and reconstruct the three-dimensional density imaging results of the target area based on the spatial domain surface density distribution of each depth layer. The expression for the classical two-dimensional Radon transform is: in, For gravity observation data, For the Dirac function, The horizontal coordinate of the spatial domain, The vertical coordinate of the spatial domain, For Radon domain scale parameters, For Radon domain angle parameters, for The classical Radon transform results; The expression for the improved Radon field definition is: in, This represents the Hilbert transform derivative convolution kernel. This represents the convolution operator. To improve the results of the Radon transform.

Images