A density anomaly body probability identification method based on a Gaussian mixture model

Abstract

The application discloses a density anomaly body probability identification method based on a Gaussian mixture model, and belongs to the technical field of geophysical exploration and data interpretation. The method comprises the following steps: obtaining gravity anomaly observation data and discretizing underground space, and constructing a forward relationship model; a Gaussian mixture model is constructed to describe the distribution of density parameters, and a model parameter sample set is obtained by solving; the statistical characteristics of the density parameters of each grid unit are calculated, at least including mean value and standard deviation; based on the statistical characteristics and the model parameters, the probability values of each grid unit belonging to different categories are calculated through first-order approximation, and an identification result is generated; in the first-order approximation, the obtained posterior distribution and prior distribution are applied in the first-order approximation based on moments, and the transmission of model solving uncertainty to the probability calculation result is realized. The quantitative expression ability of the density anomaly body identification result and the geological interpretation reliability are improved, and the problems that the existing technology cannot directly obtain the probability distribution of the anomaly body based on the model statistical characteristics and cannot realize the transmission of uncertainty are solved.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and data analysis technology, and in particular to a method for probabilistic identification of density anomalies based on a Gaussian mixture model. Background Technology

[0002] Gravity exploration is a crucial method for obtaining information on subsurface density distribution and is widely used in mineral resource exploration and geological structure research. Its basic principle is to measure gravity anomaly data at the Earth's surface or in the air, invert the density distribution characteristics of subsurface space, and thus identify geological bodies with different density properties. However, the gravity inversion problem is essentially a geophysical inverse problem, exhibiting severe ill-posedness. This means that even small perturbations in the observational data can lead to significant changes in the inversion results, and the same set of observational data may correspond to multiple different subsurface density distribution models.

[0003] To address the aforementioned issues, existing technologies primarily employ two types of processing methods. One type is the deterministic inversion method, which establishes a mathematical mapping relationship between observed data and subsurface density parameters, and introduces constraints such as Tikhonov regularization and total variational regularization to solve for stable density distribution results. This type of method has high computational efficiency and is widely used in engineering practice. However, the introduction of Tikhonov regularization constraints often leads to overly smoothed inversion results, making it difficult to accurately characterize the boundary morphology and internal structural features of subsurface anomalies, especially for geological bodies with complex geometries or multiple superimposed phases, where its resolution is significantly insufficient. Furthermore, deterministic inversion typically provides only a single optimal estimate, failing to quantify the uncertainty of the inversion results and thus failing to meet the reliability assessment requirements of refined geological interpretation.

[0004] Another category is statistical modeling-based methods, which treat subsurface density parameters as random variables and use probability distribution models to characterize their spatial distribution features. Among these, Gaussian mixture models have gained increasing attention due to their multimodal expressive capabilities and ability to describe the density distribution differences corresponding to different types of geological bodies. Some methods further utilize techniques such as Markov chain Monte Carlo sampling to obtain posterior statistical information of the model parameters. However, existing statistical modeling methods often apply Gaussian mixture models as prior constraints to the inversion process, or only perform clustering and classification on the results after inversion. The integration of model construction and parameter solution is not tight enough, failing to fully utilize the statistical information generated during the solution process. More importantly, existing methods cannot directly calculate the probability value of each spatial location belonging to a specific geological anomaly based on the statistical characteristics of the model parameters. This results in a lack of quantitative uncertainty expression in the identification results, failing to answer the core question in geological interpretation: "What is the probability that a certain location belongs to an anomaly?" Furthermore, in high-dimensional parameter spaces, traditional random sampling methods suffer from slow convergence speeds and strong correlations in sample sequences, leading to high computational costs and making it difficult to balance identification accuracy and computational efficiency.

[0005] In summary, current technologies for identifying underground density anomalies are unable to accurately characterize the boundaries of complex geological bodies while simultaneously achieving probabilistic identification based on model statistical features and maintaining computational efficiency. There is an urgent need for a method that can solve these problems within a unified framework. Summary of the Invention

[0006] The purpose of this invention is to provide a density anomaly probabilistic identification method based on Gaussian mixture models. By constructing a Gaussian mixture model to perform multimodal modeling of density parameters, and by introducing the uncertainty of the solution into the probability calculation in a first-order approximation based on the statistical characteristics of the model parameter samples, the quantitative expression ability and geological interpretation reliability of the density anomaly identification results are improved. This solves the problem that existing technologies cannot directly obtain the probability distribution of anomalies based on the statistical characteristics of the model and cannot realize the transmission of uncertainty.

[0007] To address the aforementioned technical problems, a first aspect of this invention provides a method for probabilistic identification of density anomalies based on a Gaussian mixture model, comprising the following steps: S1: Obtain gravity anomaly observation data of the target area, and discretize the underground space of the target area into multiple grid cells to construct a forward model of the relationship between the gravity anomaly observation data and the density parameters of each grid cell; S2: Construct a Gaussian mixture model to describe the distribution of density parameters of each grid cell, and solve the Gaussian mixture model based on the forward model and the gravity anomaly observation data to obtain a model parameter sample set consisting of multiple sets of model parameter values; wherein, the Gaussian mixture model is composed of a weighted combination of multiple Gaussian components, and the model parameter values ​​include the mean, variance and weight of each Gaussian component; Step S2 involves solving the Gaussian mixture model based on the forward model and the gravity anomaly observation data to obtain a model parameter sample set consisting of multiple sets of model parameter values, including: Step S21: Construct an objective function that includes data fitting terms and model constraint terms, wherein the data fitting terms are used to measure the degree of agreement between the model forward response and the gravity anomaly observation data; Step S22: Using a gradient-based sampling method, iterative sampling is performed under the guidance of the gradient of the objective function to generate the model parameter sample set; S3: Based on the model parameter sample set, calculate the density parameter statistical characteristics corresponding to each grid cell, wherein the statistical characteristics include at least the mean and standard deviation; S4: Based on the mean and standard deviation of the density parameters corresponding to each grid cell, as well as the mean, variance and weight of each Gaussian component, the probability value of the grid cell belonging to the category corresponding to different Gaussian components is calculated by first-order approximation to generate the density anomaly identification result. The first-order approximation refers to applying a moment-based first-order approximation to the posterior and prior distributions to realize the transfer of uncertainty in model solution to probability calculation results.

[0008] Furthermore, the gradient-based sampling method described in step S22 is the Hamiltonian Monte Carlo sampling method; Step S22 employs a gradient-based sampling method, performing iterative sampling guided by the gradient of the objective function to generate the model parameter sample set, including: Step S221a: Introduce an auxiliary momentum variable to construct a Hamiltonian dynamical system with the objective function as the potential energy function; Step S222a: Simulate the motion trajectory of the Hamiltonian dynamic system using a preset numerical solution method to generate candidate samples; Step S223a: Based on the preset acceptance criteria, the candidate samples are screened, and the samples that meet the conditions are retained to form the model parameter sample set.

[0009] Further, the construction of the Gaussian mixture model used to describe the density parameter distribution of each grid cell in step S2 includes: Step S201: Obtain the geological prior information of the target area, wherein the geological prior information includes at least the background medium type and the potential anomaly type; Step S202: Determine the number of Gaussian components in the Gaussian mixture model based on the geological prior information, and configure an initial mean and initial variance for each Gaussian component to characterize the density distribution characteristics of the corresponding geological body type; Step S203: Based on the configured Gaussian components, construct the complete Gaussian mixture model.

[0010] Further, in step S4, based on the mean and standard deviation of the density parameter corresponding to each grid cell, and the mean, variance, and weight of each Gaussian component, the probability value of the grid cell belonging to the category corresponding to different Gaussian components is calculated by first-order approximation, including: Step S411: Based on the mean, variance and weight of each Gaussian component, calculate the initial probability density of the grid cell under each category; Step S412: Introduce the standard deviation into the variance parameter of the initial probability density, and obtain the corrected probability density by superimposing the variances; Step S413: Normalize the corrected probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

[0011] Further, step S412, which involves incorporating the standard deviation into the variance parameter of the initial probability density and obtaining the corrected probability density through variance superposition, includes: Step S4121: Obtain the standard deviation of the density parameter corresponding to the grid cell and the variance of the Gaussian component corresponding to the current category; Step S4122: Perform a sum of squares operation on the standard deviation and the standard deviation of the Gaussian component, and use the result of the sum of squares operation as the corrected variance parameter to construct the corrected probability density.

[0012] Further, the normalization process for the corrected probability density corresponding to each category in step S413 includes: Step S4131: Calculate the sum of the corrected probability densities of the grid cell under all categories, and use it as the normalized denominator; Step S4132: Divide the corrected probability density corresponding to each category by the normalized denominator to obtain the probability value of the grid cell belonging to that category.

[0013] Furthermore, the statistical features also include the covariance matrix; Step S4, which involves calculating the probability value of the mesh cell belonging to the category corresponding to different Gaussian components using a first-order approximation, includes: Step S4141: Obtain the density parameter covariance matrix corresponding to each grid cell, and the variance of the Gaussian component corresponding to the current category; Step S4142: Fuse the covariance matrix with the variance of the Gaussian components to construct a corrected covariance matrix; Step S4143: Based on the corrected covariance matrix, calculate the corrected probability density of the grid cell under the current category; Step S4144: Normalize the modified probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

[0014] Further, the generation of density anomaly identification results in step S4 includes: Step S421: For each grid cell, obtain its probability value belonging to different categories, and compare it with the preset probability threshold corresponding to each category respectively; Step S422: Simultaneously label the grid cell as all categories that meet the probability threshold conditions to form a multi-category joint identification result, which characterizes the grid cell's probability of belonging to different geological body types.

[0015] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described density anomaly probabilistic identification method based on Gaussian mixture model.

[0016] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described method for probabilistic identification of density anomalies based on a Gaussian mixture model.

[0017] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By employing a gradient-based sampling method to efficiently solve the model parameters of a Gaussian mixture model, and using the gradient of the objective function to guide the sampling path, a rapid exploration of the model parameter distribution is achieved in the high-dimensional parameter space of gravity inversion. This significantly improves the sampling convergence speed and reduces the correlation of sample sequences. Thus, while ensuring the representativeness of the model parameter sample set, it effectively solves the problem of low computational efficiency of traditional random sampling methods in solving high-dimensional density parameters, providing a reliable data foundation for subsequent uncertainty analysis and probability identification. 2. By explicitly introducing the statistical characteristics, especially the standard deviation, obtained from solving the model parameters during the probability calculation process into the probability density function of the Gaussian component in a first-order approximation manner, the uncertainty generated in the model solving stage is quantitatively transferred to the anomaly identification stage. This results in the final output no longer being a single deterministic classification result, but rather a specific probability value for each spatial grid cell to belong to different geological body categories. This solves the problem that existing technologies cannot directly obtain the probability distribution of anomalies based on model statistical characteristics and that the identification results lack quantitative uncertainty expression, significantly improving the reliability of the geological interpretation of density anomaly identification results. 3. By establishing a correspondence between each Gaussian component in the Gaussian mixture model and the geological prior information of the target area, and determining the number of Gaussian components and configuring their initial parameters according to the background medium type and potential anomaly type, multimodal density parameter modeling is no longer a purely mathematical distribution fitting, but a physical modeling process with clear geological meaning. The probability values ​​of each grid cell belonging to the corresponding category of different Gaussian components can be directly interpreted as the probability of belonging to different geological body types. This solves the problem of the disconnect between existing statistical modeling methods and geological cognition, and improves the correspondence and interpretability between the identification results and geological reality. Attached Figure Description

[0018] Figure 1 This is a flowchart of the density anomaly probability identification method based on Gaussian mixture model provided in an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0020] Please refer to Figure 1 The first aspect of this invention provides a method for probabilistic identification of density anomalies based on a Gaussian mixture model, comprising the following steps: S1: Obtain gravity anomaly observation data of the target area and discretize the underground space of the target area into multiple grid cells to construct a forward model of the relationship between gravity anomaly observation data and the density parameters of each grid cell.

[0021] Within the target area for gravity exploration, gravity anomaly observation data is acquired through field measurements or data acquisition systems. This data reflects the changes in the gravity field at different underground locations due to differences in density distribution, serving as fundamental input information for inferring the underground density structure. To establish a computable mathematical relationship, the underground space of the target area is discretized into multiple grid cells in three dimensions. Each grid cell represents a specific spatial location, and its internal density parameters are assumed to be uniform or represent the equivalent value of that cell. Based on gravity forward modeling theory, a forward modeling relationship between gravity anomaly observation data and the density parameters of each grid cell is constructed, using the density parameters of each grid cell as unknowns. This forward modeling relationship describes the physical process of "given the density distribution of each underground grid cell, what kind of gravity anomaly response occurs on the surface," generally expressed as a discrete form of a system of linear equations or linear integral equations. Each element in the coefficient matrix represents the contribution of unit density at the corresponding grid cell to the gravity response at a specific observation point. This provides a physical framework for the connection between the parameter space and the observation data for subsequent steps. In practice, the size of the grid cells can be set according to the geological complexity of the target area, the expected size of the anomaly, and the spatial sampling density of gravity observation data. For example, a larger grid can be used in sedimentary basins with relatively uniform surface cover, while a smaller grid is required in areas such as metal mining areas to improve resolution.

[0022] S2: Construct a Gaussian mixture model to describe the distribution of density parameters in each grid cell, and solve the Gaussian mixture model based on the forward model and gravity anomaly observation data to obtain a sample set of model parameters consisting of multiple sets of model parameter values. The Gaussian mixture model is composed of a weighted combination of multiple Gaussian components, and the model parameter values ​​include the mean, variance, and weight of each Gaussian component.

[0023] Considering that subsurface media are typically composed of various types of lithology or geological bodies, the density parameters of these different geological bodies often exhibit differentiated statistical distribution characteristics. Step S2 constructs a Gaussian mixture model to describe the density parameter distribution of all grid cells. A Gaussian mixture model is a probabilistic model composed of multiple Gaussian components linearly combined with certain weights. Each Gaussian component is defined by three model parameters: its mean, variance, and weight within the entire model. Specifically, the mean of each Gaussian component represents the central tendency of the density parameters of that geological body type, the variance represents the dispersion of density parameters within that geological body type, and the weight represents the volume proportion of that geological body type within the entire study area. This multi-Gaussian component combination modeling method can characterize the multimodal distribution characteristics of density distributions of different types of geological bodies in the subsurface media, such as the density peak corresponding to the background surrounding rock, the high density peak corresponding to the target ore body, and the low density peak corresponding to the alteration zone.

[0024] After constructing the Gaussian mixture model, the model parameters need to be solved under the constraints of gravity anomaly observation data and the forward model. The solution process involves searching the model parameter space for parameter combinations that best match the model's forward response to the gravity anomaly observation data, resulting in a sample set of model parameters consisting of multiple sets of parameter values. This sample set does not point to a single optimal estimate but contains a certain number of parameter combinations that can be statistically interpreted reasonably from the observation data, thus reflecting the variable range and statistical distribution characteristics of the model parameters. In the high-dimensional parameter space of gravity inversion, obtaining the sample set of model parameters provides a direct data foundation for subsequent analysis of the uncertainties in the model parameters.

[0025] S3: Based on the model parameter sample set, calculate the statistical characteristics of the density parameters corresponding to each grid cell. The statistical characteristics include at least the mean and standard deviation.

[0026] Based on the model parameter sample set obtained in step S2, statistical analysis is performed on the density parameters corresponding to each grid cell. Since the model parameter sample set contains different model parameter values ​​obtained from multiple solutions or sampling, each grid cell corresponds to multiple possible density parameter values. Calculating the statistical characteristics of these values ​​quantifies the estimation tendency and uncertainty range of the density parameter at that location. As basic statistical characteristics, the mean and standard deviation are calculated: the mean represents the expected estimate of the density parameter of that grid cell, reflecting the central tendency of multiple solutions; the standard deviation measures the dispersion of the density parameter values ​​obtained from each solution around the mean, i.e., the uncertainty level of the model parameter at that location.

[0027] In real-world geological scenarios, the uncertainty of model parameters can stem from multiple sources, including noise levels in gravity anomaly observation data, approximation errors introduced by grid discretization, and the ill-posedness of the inversion problem itself. By using the mean and standard deviation, these multiple sources of uncertainty are condensed into two quantifiable numerical indicators on each grid cell, preparing the ground for further propagating model uncertainty to the anomaly identification stage in subsequent steps. In exploration projects with high requirements for interpretability, areas with larger standard deviations often indicate lower reliability of the inferred density parameters at those locations.

[0028] S4: Based on the mean and standard deviation of the density parameter corresponding to each grid cell, and the mean, variance, and weight of each Gaussian component, the probability value of the grid cell belonging to the corresponding category of different Gaussian components is calculated through a first-order approximation, generating the density anomaly identification result. The uncertainty information obtained in step S3 is quantitatively transferred to the anomaly identification result, realizing the probabilistic discrimination of anomaly categories at each spatial location. Specifically, for each grid cell, based on the mean and standard deviation of its density parameter, combined with the original mean, variance, and weight of each Gaussian component in the Gaussian mixture model, the probability value of the grid cell belonging to the corresponding category of different Gaussian components is calculated using a moment-based first-order approximation method.

[0029] The first-order approximation refers to applying a moment-based first-order approximation to the posterior and prior distributions, thus transferring the uncertainty of the model solution to the probability calculation results. The moment-based first-order approximation substitutes the mean as the first-order moment estimate of the grid cell density parameter into the location parameter of the Gaussian component to characterize the current best estimate of the grid cell. Simultaneously, when constructing the relationship between the posterior and prior distributions, instead of performing a complete analytical solution for the posterior distribution or costly numerical integration, it utilizes the mean and standard deviation—two statistical moments sufficient to characterize the main features of the distribution—to perform a first-order approximation of the posterior distribution. This approach ensures that the superimposed probability density function reflects two parts of uncertainty simultaneously: first, the inherent dispersion of the Gaussian component's density distribution across different geological bodies; and second, the model uncertainty remaining from the estimation of the grid cell density parameter during step S2. Based on this, the probability density of the grid cell under each category is calculated, and the probability densities corresponding to all categories are normalized to obtain the probability values ​​of the grid cell belonging to each category. The higher the calculated probability value, the higher the confidence that the grid cell belongs to that category.

[0030] Based on the probability values ​​of each grid cell belonging to different categories, density anomaly identification results can be generated. These results not only indicate the spatial location of the anomaly but also quantitatively express the confidence level of each location belonging to a specific anomaly category in the form of a probability distribution, thus adding an uncertainty dimension to traditional classification results. For example, at a certain grid cell, if the calculated probability of belonging to the high-density ore body category is high, while the probability of belonging to the background surrounding rock category is low, the location can be identified as a mineralized anomaly with high confidence. Conversely, if the probabilities of each category are close, it indicates that the classification results have significant uncertainty, and subsequent interpretation should be based on a comprehensive judgment combined with other geological data.

[0031] Through the above steps, this method uses a Gaussian mixture model to perform multimodal modeling of subsurface density parameters. In the probability calculation stage, the statistical uncertainty generated during the model solution process is explicitly introduced into the identification results through a first-order approximation based on moments. This allows the anomaly category of each spatial location to be presented in the form of a probability distribution. At the same time, it takes into account the sample acquisition efficiency in the high-dimensional parameter space of the solution stage, and improves the uncertainty expression ability of the density anomaly identification results and the reliability of geological interpretation.

[0032] Specifically, in step S2, the Gaussian mixture model is solved based on the forward modeling model and gravity anomaly observation data to obtain a model parameter sample set consisting of multiple sets of model parameter values, including: Step S21: Construct an objective function that includes data fitting terms and model constraint terms. The data fitting terms are used to measure the degree of agreement between the model's forward response and the gravity anomaly observation data.

[0033] In gravity inversion problems, directly solving for the density parameters of each grid cell often faces difficulties such as underdetermined or ill-conditioned equations, meaning that a large number of different parameter combinations can mathematically fit the observed data. To impose reasonable constraints on the solution process and select physically meaningful parameter combinations, step S21 constructs an objective function consisting of two components: a data fitting term and a model constraint term. The data fitting term quantifies the degree of difference between the model's forward response and the observed gravity anomaly data; that is, the deviation between the theoretical gravity response calculated based on the current model parameters and the field-measured gravity anomaly value. The smaller the deviation, the stronger the interpretability of the current model parameters for the observed data. In practice, the data fitting term can be constructed as the sum of squares of the residuals between the observed data and the forward response, or a more refined form such as weighted least squares can be used based on the noise statistical characteristics of the observed data.

[0034] The role of model constraints is to quantify prior knowledge or expected properties regarding the distribution of subsurface density parameters into mathematical constraints that are introduced into the solution process, limiting the search range of invalid or non-physical regions in the model parameter space. Within the Gaussian mixture model framework constructed in step S2, model constraints can be designed based on the distribution characteristics of Gaussian components. For example, constraints can be imposed on the separation between the means of each Gaussian component to avoid excessive overlap of different components in the parameter space, thus losing the ability to distinguish different geological body types. Appropriate constraints can also be imposed on the spatial continuity or smoothness of the density parameters of each grid cell to suppress spurious high-frequency oscillations introduced by grid discretization. By combining the data fitting terms and model constraint terms with appropriate weighting coefficients into a unified objective function, the solution process establishes a balance between fitting the observed data and adhering to prior constraints, ensuring that each set of parameter values ​​in the model parameter sample set can both form a reasonable mathematical fit to the observed data and conform to a geologically reasonable distribution.

[0035] Step S22: Using a gradient-based sampling method, iterative sampling is performed under the guidance of the gradient of the objective function to generate a set of model parameter samples.

[0036] Once the objective function is constructed, it is necessary to find parameter combinations that minimize the objective function value within the high-dimensional model parameter space, which is composed of the density parameters of each grid cell and the parameters of each Gaussian component, and generate a statistically representative set of model parameter samples. Since the dimension of the model parameter space is typically proportional to the number of grid cells and Gaussian components, in embodiments with large regions or fine mesh divisions, the parameter space dimension can reach thousands or even higher. Traditional random sampling methods suffer from slow convergence speed and strong statistical correlation between adjacent samples in such high-dimensional spaces, making it difficult to obtain a sample set that stably reflects the parameter distribution characteristics within an acceptable computational time.

[0037] To address the aforementioned issues, step S22 employs a gradient-based sampling method, iteratively sampling under the guidance of the gradient direction of the objective function on each model parameter. The gradient information of the objective function indicates which direction each model parameter should be adjusted in at the current parameter position to cause the objective function value to decrease most steeply, i.e., pointing to the parameter region with higher probability density. By using gradient information to guide the sampling direction and step size in each iteration, the sampling path can be more effectively concentrated in the parameter space, targeting regions that have a high degree of agreement with the observed data and satisfy the model constraints, avoiding ineffective wandering in the empty space outside the high-probability region. In a typical iteration, candidate update quantities are calculated based on the current model parameter value and the gradient information at that location, generating a set of candidate parameter values. Then, a preset probability criterion is used to determine whether to accept the candidate parameter value as a new member of the sample set. After multiple iterations, multiple sets of parameter values ​​sufficient to characterize the distribution features of the model parameters can be collected, forming a model parameter sample set. Each set of parameter values ​​in this sample set statistically follows a probability distribution defined by the objective function, thus reflecting the possible range and distribution pattern of each model parameter under given gravity anomaly observation data and model constraints.

[0038] By constructing an objective function that includes data fitting terms and model constraint terms, and obtaining a set of model parameter samples through iterative sampling guided by gradient information, the above solution process balances sampling efficiency and sample representativeness in a high-dimensional parameter space. This enables the acquisition of multiple sets of parameter values ​​that effectively reflect the statistical distribution characteristics of model parameters at a reasonable computational cost in practical applications of gravity inversion, providing reliable data support for subsequent steps of statistical feature calculation and uncertainty propagation.

[0039] In an optional embodiment of the present invention, the sampling method based on gradient information in step S22 is the Hamilton Monte Carlo sampling method.

[0040] Accordingly, step S22 employs a gradient-based sampling method, performing iterative sampling guided by the gradient of the objective function to generate a model parameter sample set, including: Step S221a: Introduce an auxiliary momentum variable to construct a Hamiltonian dynamical system with the objective function as the potential energy function.

[0041] Within the framework of the Hamiltonian Monte Carlo sampling method, the sampling process of model parameters is modeled as a dynamic system. Specifically, a set of auxiliary momentum variables is introduced for each model parameter to be solved (i.e., the mean, variance, and weight of each Gaussian component in the Gaussian mixture model, and the density parameters of each grid cell). These momentum variables only serve the numerical calculation of the sampling process and have no actual physical meaning; their role is to endow the parameter space with dynamic evolution capabilities. Based on this, the objective function constructed in step S21 is regarded as a potential energy function, which, together with the momentum variables, constitutes a Hamiltonian dynamic system. In this system, the potential energy function is determined by the objective function and represents the "height" of different positions in the parameter space. Regions with smaller objective function values ​​correspond to lower potential energy, which are more likely to be reasonable values ​​for model parameters; the momentum variables contribute the kinetic energy of the system. Through this construction, the sampling problem in the model parameter space is transformed into an evolution problem of a mechanical system in phase space, providing a mathematical framework for subsequent efficient trajectory exploration using gradient information. In high-dimensional embodiments with a large number of grids, such as metal ore exploration, although the phase space dimension of the Hamiltonian system is further multiplied, its dynamic structure allows the exploration of the parameter space to proceed along a physical path determined by the gradient, rather than a random blind walk.

[0042] Step S222a: Simulate the motion trajectory of the Hamiltonian dynamic system using a preset numerical integration method to generate candidate samples.

[0043] After constructing the Hamiltonian dynamics system, a pre-defined numerical method is used to numerically simulate the system's trajectory. The choice of numerical method must balance computational accuracy and efficiency; in practice, the leapfrog method or a higher-order symplectic scheme can be used. The numerical simulation starts from the existing model parameters and corresponding momentum variable values, proceeding along the dynamic path of the Hamiltonian system according to a pre-defined time step and number of integration steps. The model parameter and momentum variable values ​​at the end of the trajectory are calculated, and these end-point model parameter values ​​are used as candidate samples. Since the dynamic evolution of the Hamiltonian system is essentially along the gradient direction of the objective function, the simulated trajectory can traverse regions with low probability density in the parameter space, maintaining a large displacement interval between candidate samples from consecutive sampling, thus effectively reducing the statistical autocorrelation of the sample sequence. In practical applications, the time step and number of integration steps can be adjusted according to the parameter space size determined by the number of grid cells and Gaussian components in the target region to balance the computational cost of trajectory simulation with the sufficiency of sample exploration.

[0044] Step S223a: Based on the preset acceptance criteria, the candidate samples are screened, and the samples that meet the conditions are retained to form a model parameter sample set.

[0045] Due to the discretization error inherent in numerical methods, not all candidate samples simulated along the Hamiltonian dynamic trajectory precisely satisfy the probability distribution defined by the objective function. Therefore, step S223a filters candidate samples based on a preset acceptance criterion. In practice, the Metropolis acceptance criterion, using an equal-probability testing mechanism, determines whether to include a current candidate sample in the model parameter sample set with a certain probability. The acceptance criterion's determination process introduces a probabilistic test, making candidate samples generated by trajectory segments with more accurate numerical integration more likely to be retained, while candidate samples with larger integration errors have a corresponding probability of being rejected. After screening and discrimination, the retained samples statistically constitute a sample set that follows the posterior distribution of the parameters defined by the objective function. As the number of iterations increases, this sample set gradually accumulates a sufficient number of model parameter values, covering the main probabilistic mass regions of density parameters in high-dimensional space, ultimately forming a model parameter sample set that reflects the distribution characteristics of density volume parameters for each category, providing a sufficient data foundation for the calculation of statistical features in step S3.

[0046] By constructing the solution process of model parameters as a Hamiltonian dynamic system, effectively exploring the parameter space using gradient-guided trajectory simulation, and ensuring the statistical rationality of the samples with the help of a preset acceptance criterion, this invention achieves an effective balance between sampling efficiency and sample quality in the high-dimensional parameter space of gravity inversion. The generated model parameter sample set has low autocorrelation and sufficient spatial coverage, which helps to improve the accuracy and stability of subsequent uncertainty statistical inference.

[0047] As another optional implementation of this invention, the gradient-based sampling method in step S22 is the Langevin dynamics sampling method. This method samples the model parameter space by performing a noisy random walk along the gradient direction of the objective function. Specifically, based on the current model parameter values ​​and their corresponding objective function gradient information, a Langevin diffusion equation with random noise terms is constructed to describe the dynamic process of the random walk of model parameters in the parameter space. In this dynamic process, the update direction of the model parameters is determined by two parts: the first part is the negative gradient direction of the objective function at the current parameter position, which guides the sampling towards the parameter region with lower objective function values ​​and higher probability density; the second part is the Gaussian random noise term, which ensures the ergodicity of the sampling process in the parameter space, enabling the sampling path to cover the main probability mass region of the target distribution and avoiding premature entry into local extrema. Subsequently, the Langevin diffusion equation is discretized and solved using a preset discretized numerical format to generate candidate samples. In specific implementations, common solution methods for stochastic differential equations, such as the Euler-Maruyama method, can be used. Starting from the current model parameter values, the Langevin diffusion equation is numerically integrated in one or more steps according to a preset time step to obtain new model parameter values ​​as candidate samples. Since numerical discretization may introduce approximation errors, the distribution of candidate samples may deviate from the target distribution. Therefore, candidate samples are screened based on a preset acceptance criterion. The Metropolis acceptance criterion, an equal probability test mechanism, determines whether to include the current candidate sample in the model parameter sample set. After multiple iterations, a sufficient number of retained model parameter values ​​are collected to form the model parameter sample set. Compared with the Hamiltonian Monte Carlo sampling method, the Langevin dynamic sampling method does not require the introduction of auxiliary momentum variables, and the calculation steps are simpler. It has a computational efficiency advantage in application scenarios with extremely high parameter space dimensionality and relatively mild objective functions.

[0048] As another optional implementation of this invention, the sampling method based on gradient information in step S22 is the Stein variational gradient descent method. This method approximates the target posterior distribution by maintaining a set of model parameter particles and using the gradient information of the objective function to drive the collaborative updating of each particle. Specifically, firstly, a set of particles containing multiple sets of model parameter values ​​is initialized. Each set of model parameter values ​​includes the mean, variance, and weight of each Gaussian component, as well as the density parameter corresponding to each grid cell. The initial value of each particle can be set based on the preliminary analysis results of geological prior information and gravity anomaly observation data of the target area, or by random initialization within a reasonable range of the parameter space. Subsequently, the particle set is iteratively updated. In each iteration, the gradient direction corresponding to each particle is calculated based on the gradient information of the objective function for the current set of model parameter values, and the gradient direction is corrected using the Stein operator to obtain the update direction of each particle. The correction ensures that the update directions of different particles continue to approach the high probability density region, while introducing mutual repulsion between particles through the kernel function to maintain the overall diversity of the particle set and prevent particles from collapsing into a single mode. Each particle updates synchronously according to the update direction. After each iteration, the objective function value of each particle is evaluated. The iteration terminates when the entire particle set satisfies the preset convergence condition or reaches the preset upper limit of the number of iterations, and the particle set is used as the model parameter sample set. Unlike the aforementioned sampling method based on the Markov chain Monte Carlo framework, the Stein variational gradient descent method generates multiple sets of model parameter values ​​simultaneously through particle optimization. Each particle cooperates with each other during the iteration process, combining the statistical completeness of the sampling method with the convergence efficiency of the optimization method. In actual exploration scenarios with high parameter dimensionality and strict requirements for computation time, it can obtain a model parameter sample set that reasonably reflects the probability distribution characteristics of the target within a fewer iterations.

[0049] It should be noted that the aforementioned gradient-based sampling methods share a common technical feature: they utilize the first-order gradient information of the objective function on the model parameters to guide the sampling process, causing the sampling path to favor regions with higher probability density in the parameter space, thereby achieving efficient sampling in high-dimensional gravity inversion problems. All of the above methods can achieve the technical effects of improved sampling efficiency and enhanced sample representativeness claimed in this invention; these alternative implementations are all extensions of the core concept of this invention.

[0050] Further, step S2 involves constructing a Gaussian mixture model to describe the distribution of density parameters for each grid cell, including: Step S201: Obtain the geological prior information of the target area. The geological prior information includes at least the background medium type and the type of potential anomalies.

[0051] Before constructing a Gaussian mixture model, prior geological information about the target area must first be obtained. This prior geological information comes from existing regional geological survey data, borehole logging results, laboratory data from core samples, or geophysical interpretation profiles from adjacent known mining areas. It includes at least an understanding of the background medium type in the subsurface space of the target area and the identification of potential anomaly types. The background medium type typically refers to the widely distributed surrounding rocks or stratigraphic units within the target area, whose density parameters constitute the baseline density level for the region. Potential anomaly types refer to geological bodies related to mineralization, tectonic activity, or alteration, such as high-density massive sulfide ore bodies, intrusive intrusive bodies, or low-density fault fracture zones, salt dome structures, etc., whose density parameters deviate significantly from the background medium. The purpose of obtaining this prior geological information is to provide physically grounded structural guidance for constructing the Gaussian mixture model, ensuring that the model is not merely a mathematical distribution fitting tool, but incorporates the geophysical exploration worker's targeted understanding of the geological background of the study area.

[0052] Step S202: Determine the number of Gaussian components in the Gaussian mixture model based on the geological prior information, and configure the initial mean and initial variance for each Gaussian component to characterize the density distribution characteristics of the corresponding geological body type.

[0053] Based on the geological prior information obtained in step S201, the key structural parameters of the Gaussian mixture model are configured in a targeted manner. First, the number of Gaussian components in the Gaussian mixture model is determined according to the total number of background medium types and potential anomaly types in the geological prior information. For example, if one main surrounding rock and two different types of mineralization anomalies are known to exist in the study area, the number of Gaussian components can be set to three, so that each Gaussian component physically corresponds to a known or inferred geological body type. Second, an initial mean and initial variance characterizing the density distribution of the corresponding geological body type are configured for each Gaussian component. The initial mean can be set with reference to the typical density values ​​of the corresponding geological body type in existing borehole or laboratory measurements; for example, the mean for background sedimentary rocks is set to approximately 2,600 kg / m³, and the mean for massive sulfide ore bodies is set to over 4,000 kg / m³. The initial variance setting reflects a priori estimation of the density variation of this type of geological body. In embodiments lacking detailed variation information, it can be adjusted according to the degree of heterogeneity of the geological body type. For example, a smaller initial variance is assigned to relatively homogeneous sedimentary layers, while a larger initial variance is assigned to complex ore bodies that have undergone multiple phases of alteration after mineralization, to accommodate a greater range of density variations. Through this configuration process, each Gaussian component obtains a clear geological attribute label, establishing a direct correspondence between its mathematical parameters and the physical characteristics of a specific geological body type.

[0054] Step S203: Construct a complete Gaussian mixture model based on the configured Gaussian components.

[0055] After the quantity, initial mean, and initial variance of each Gaussian component have been configured based on prior geological information, all Gaussian components are weighted and combined according to their respective weight parameters to construct a complete Gaussian mixture model. The weight parameters of each Gaussian component can also be initially assigned based on the estimated volume proportion of different types of geological bodies in the study area from the prior geological information. Thus, the constructed Gaussian mixture model is mathematically a probability density function composed of multiple Gaussian distributions superimposed with weights; geologically, it is a parameterized expression of different density types of geological bodies in the subsurface space of the target area. When this model enters the subsequent solution stage in step S2, the mean, variance, and weight of each Gaussian component will be used as model parameters to be solved, further updated and refined under the constraints of gravity anomaly observation data and forward modeling to reflect the density distribution information contained in the observation data.

[0056] By acquiring prior geological information and determining the number of Gaussian components and initial parameters accordingly, the constructed Gaussian mixture model establishes a semantic correspondence between each Gaussian component and a specific geological body type during the modeling stage. This allows the probability value of each grid cell belonging to the corresponding category of different Gaussian components calculated in subsequent steps to be directly interpreted as the probability of belonging to the corresponding geological body type, thereby enhancing the physical interpretability of the modeling process and the correspondence between the final identification results and geological reality.

[0057] Specifically, in step S4, based on the mean and standard deviation of the density parameter corresponding to each grid cell, and the mean, variance, and weight of each Gaussian component, the probability value of the grid cell belonging to the corresponding category of different Gaussian components is calculated using a first-order approximation, including: Step S411: Calculate the initial probability density of the grid cell under each category based on the mean, variance and weight of each Gaussian component.

[0058] After obtaining the statistical characteristics of the density parameters corresponding to each grid cell, the category probability calculation for that grid cell is initiated. Based on the mean, variance, and weight parameters of each Gaussian component in the Gaussian mixture model constructed and solved in step S2, the initial probability density of the current grid cell under the category corresponding to each Gaussian component is calculated. Each Gaussian component defines a probability density function centered on its own mean and with its own variance as the width. The calculation of the initial probability density essentially measures the likelihood of the density parameter of the grid cell appearing under the current category density distribution model. In specific implementation, if a Gaussian component corresponds to the background surrounding rock type, its mean is within the common rock density range, and its variance is small, then grid cells with density parameters close to this mean will obtain a higher initial probability density under this component; conversely, if the density parameter of the grid cell deviates significantly from the mean of this component, the initial probability density will be lower. Similarly, the Gaussian component corresponding to a high-density ore body will assign a higher initial probability density to the high-density grid cell. The weights of each Gaussian component are used at this stage to adjust the basic contribution ratio of different categories in the overall model. Components with larger weights provide higher basic probability density weights for their corresponding categories.

[0059] Step S412: Introduce the standard deviation into the variance parameter of the initial probability density, and obtain the corrected probability density by summing the variances.

[0060] The initial probability density calculated in step S411 only reflects the distribution characteristics of each component of the Gaussian mixture model and does not yet include the uncertainty remaining in the density parameter estimation from step S2. This uncertainty is introduced into the probability calculation through variance superposition. Specifically, the standard deviation of the density parameter of the grid cell is introduced into the variance parameter of the Gaussian component used in the initial probability density, so that the superimposed variance parameter contains two parts of information: one is the inherent dispersion of the density of this type of geological body represented by the variance of the Gaussian component under investigation, and the other is the uncertainty of the density parameter estimation of the grid cell represented by the square of the standard deviation calculated in step S3. The essence of variance superposition is to merge two types of uncertainties with different sources and different properties under the same probability density calculation framework. The superimposed variance parameter is larger than the original Gaussian component variance, and the corresponding probability density function spread range is correspondingly widened, reflecting a mitigating effect on the lack of confidence in the density parameter estimation of this grid cell in terms of morphology. In geological exploration practice, when a grid cell is located in an area with sparse borehole control or low signal-to-noise ratio of observation data, the standard deviation of the density parameter of the grid cell is often large. After variance superposition correction, the probability density distribution of the grid cell under each category will show a more gradual shape, and the probability difference between each category will tend to be mild, which intuitively reflects the increased classification uncertainty caused by insufficient data constraints.

[0061] Step S413: Normalize the corrected probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

[0062] After step S412 corrects the variance of the probability densities corresponding to each Gaussian component through variance superposition, the resulting corrected probability density values ​​cannot be directly used as the probability basis for category discrimination because they have not been normalized, and the corrected probability densities corresponding to different categories are not comparable in magnitude. The corrected probability densities calculated for the current grid cell under all categories are normalized, converting the corrected probability density values ​​corresponding to each category into a probability value sequence whose sum is a fixed constant. The normalized value is the probability value of the grid cell belonging to the category corresponding to different Gaussian components, ranging from zero to one, and the sum of the probabilities of all categories is one. In actual geological interpretation, for a given grid cell, if the probability value of belonging to the high-density ore body category is close to 0.9, and the probability value of belonging to the background surrounding rock category is 0.1, the grid cell can be identified as a mineralized anomaly with relatively clear confidence. If the two probability values ​​are approximately 0.45 and 0.55 respectively, it indicates that there is significant uncertainty in the category at this location, and the identification conclusion should be adopted with caution. By using the normalized probability distribution, different grid cells can be compared horizontally under a unified probability metric, providing comparable and quantifiable statistical basis for subsequent spatial localization and reliability assessment of anomalies.

[0063] Through the three interconnected steps of initial probability density calculation, variance superposition correction, and normalization, step S4 integrates the statistical characteristics of Gaussian mixture model parameters, grid cell density parameters, and uncertainties in the model solving process within a unified probabilistic framework. This ensures that the final probability value not only reflects the prior distribution characteristics of various geological bodies but also explicitly incorporates the uncertainties in the model solving process into the discrimination criteria. This achieves the quantitative transmission of uncertainty from data inversion to identification results, enabling the density anomaly identification results to have a quantifiable probabilistic interpretation capability.

[0064] Further, in step S412, the standard deviation is introduced into the variance parameter of the initial probability density, and the corrected probability density is obtained by summing the variances, including: Step S4121: Obtain the standard deviation of the density parameter corresponding to the mesh cell and the variance of the Gaussian component corresponding to the current category.

[0065] Before performing the specific calculation of variance superposition, step S4121 first prepares and acquires two sets of input data. The first set of data is the standard deviation of the density parameter corresponding to the current grid cell. This standard deviation comes from the statistical characteristics of the grid cell calculated based on the model parameter sample set in step S3. It reflects the degree of dispersion of the estimated density parameter value of the grid cell around its mean during the multiple solutions in step S2. The larger the standard deviation, the more dispersed the solution results of the density parameter at that location, that is, there is a large uncertainty in the model's inference of the density value at that location. The second set of data is the variance of the Gaussian component corresponding to the currently examined category. This variance is one of the inherent parameters of the Gaussian mixture model in step S2. It characterizes the natural variation range of the density parameter within this type of geological body. It belongs to the model's prior estimate of the dispersion of the density distribution of the geological body or the estimate after data correction. In the actual context of geophysical exploration, this variance may come from regional rock property statistics, statistical analysis of borehole core measured data, or analogical experience from nearby known deposits. This step explicitly juxtaposes two discreteness metrics with different physical origins, preparing the data layer for subsequent uncertainty fusion within the same computational framework.

[0066] Step S4122: Perform a sum of squares operation on the standard deviation and the standard deviation of the Gaussian component, and use the result of the sum of squares operation as the corrected variance parameter to construct the corrected probability density.

[0067] The standard deviation of the current grid cell density parameter obtained in step S4121 is summed with the variance of the Gaussian component of the current category. That is, the square of the standard deviation plus the variance of the Gaussian component is used as the corrected variance parameter. From a mathematical and statistical perspective, the square of the standard deviation is the estimated variance of the grid cell density parameter. The summation of the squares of this standard deviation and the variance of the Gaussian component itself is equivalent to superimposing the variances of two independent sources of dispersion. The result of this superposition is numerically equal to the sum of the two variances, fully encompassing information from both the estimation uncertainty introduced during the model parameter solution process and the density distribution dispersion of the geological body type itself. In practical applications, if a grid cell is located in an area with dense observation points and high signal-to-noise ratio, its density parameter standard deviation is usually small. The result of the sum of squares operation is close to the original Gaussian component variance, and the correction magnitude is limited, indicating that the uncertainty in solving the problem has a weak impact on the classification conclusion of the grid cell. Conversely, if the grid cell is located in a sparsely observed deep region, its standard deviation is large. After the sum of squares operation, the corrected variance will be significantly greater than the original Gaussian component variance. In subsequent probability density calculations, this manifests as a significant broadening of the distribution shape, which in turn reduces the probability difference between different categories for that grid cell. This reflects the underlying logic that insufficient data information leads to reduced classification discriminative power. Based on the corrected variance parameter obtained from the above sum of squares operation, the corrected probability density of the grid cell under the current category is reconstructed, completing the quantitative transfer of uncertainty information to the probability calculation stage.

[0068] Further, step S413 normalizes the corrected probability density corresponding to each category, including: Step S4131: Calculate the sum of the corrected probability densities of the grid cell under all categories, and use it as the normalized denominator.

[0069] After calculating the corrected probability density for all categories, the corrected probability densities of the current grid cell under each category are aggregated and summed. Specifically, the corrected probability density values ​​of the grid cell obtained in step S412 corresponding to the background rock category, various potential anomaly categories, and other possible geological body categories are summed one by one to obtain the sum of these corrected probability densities, and this sum is used as the denominator for subsequent normalization operations. This summation operation summarizes the comprehensive probability density level of the grid cell under all Gaussian components of the current Gaussian mixture model, and its magnitude reflects the degree to which the density parameter of the grid cell is explained by the current overall model. In real-world geological scenarios, if a grid cell is located in the transitional boundary region between anomalies and surrounding rocks, its density parameter may be explained to a moderate degree by multiple Gaussian components. In this case, the denominator of the summation will contain moderate contributions from multiple components. If the density parameter of the grid cell is extremely low in probability density under all Gaussian components due to extreme numerical anomalies, the denominator of the summation will approach a small value. This situation suggests that the density characteristics of the grid cell exceed the range of geological body types covered by the current Gaussian mixture model, and may correspond to unknown geological types or data anomalies.

[0070] Step S4132: Divide the corrected probability density corresponding to each category by the normalized denominator to obtain the probability value of the grid cell belonging to that category.

[0071] After obtaining the normalized denominator, the final calculation of the probability values ​​for each category is performed one by one. For each category of the current grid cell, the modified probability density value corresponding to that category is divided by the normalized denominator calculated in step S4131. The quotient is the probability value of the grid cell belonging to that category. After this division operation, the modified probability densities of each category, which were originally not directly comparable at the absolute value level, are transformed into a standardized probability sequence that takes values ​​in the interval between zero and one, and the sum of the probability values ​​of all categories is one. Each probability value directly represents the degree of confidence that the grid cell belongs to the corresponding geological body type under the current Gaussian mixture model framework, after taking into account the uncertainty of the model solution. In practical geological interpretation, the obtained probability values ​​can provide a clear basis for subsequent spatial location and reliability assessment of anomalies. For example, at a certain grid cell, if the probability of belonging to the high-density ore body category is 0.85, the probability of belonging to the background surrounding rock category is 0.10, and the probability of belonging to the low-density alteration zone category is 0.05, then this probability distribution structure provides quantifiable confidence support for inferring that the location is a mineralized anomaly. If the probability values ​​of each category are relatively close, it indicates that the classification result is uncertain, and the identification conclusion of this location should be treated cautiously based on a comprehensive review of more geological data. The normalization operation also provides a unified probability measurement benchmark for horizontal comparison between different grid cells, enabling the certainty and uncertainty of the classification results of each spatial location within the entire study area to be evaluated and displayed using a consistent standard.

[0072] Specifically, probability value The calculation formula is: ; in, For density parameters, Standard deviation, For category numbering, As weight, For category The corresponding mean, For category The corresponding standard deviation Let j be the mean of category j. Let be the standard deviation corresponding to category j, i be the voxel index, and j be the weight index. Let be the probability density function of a normal distribution. M This represents the total number of Gaussian components in the Gaussian mixture model.

[0073] In actual calculations, the determination of the above probability values ​​is essentially based on the result of correcting and normalizing the probability density functions of each component of the Gaussian mixture model using the mean and standard deviation obtained in step S3. For each grid cell, the mean density parameter represents the best estimate of the density value at that location, replacing the traditional method of directly using the parameters to be inverted; its standard deviation fully preserves the residual uncertainty in the solution process of step S2. By introducing the standard deviation into the variance parameter of the Gaussian components in a variance superposition manner, the probability density function corresponding to each component is appropriately increased in terms of spread width. This preserves more sufficient inter-class transition space for grid cells with weaker data constraints or greater solution discrepancies in the discrimination stage, avoiding premature and overly certain classification decisions when information is insufficient. This approach ensures the physical rationality and geological interpretability of the final probability values ​​from the underlying mechanism.

[0074] It should be noted that the statistical features calculated here are not limited to the mean and standard deviation. Depending on the needs of the actual geological problem, the covariance matrix of the density parameters of each grid cell can be further calculated based on the model parameter sample set to characterize the correlation structure of density parameter estimation errors between different spatial locations; or higher-order joint statistical features can be calculated among multiple grid cells to more comprehensively describe the uncertainty coupling relationship of the model solution process in multiple spatial locations. The introduction of these extended statistical features in the subsequent step S4 is similar to that of the standard deviation, both of which quantitatively transmit the uncertainty information of the corresponding dimension to the anomaly probability identification stage through variance superposition or covariance structure fusion. The various statistical feature quantities covered in the terminology explanation section can all be regarded as scalable configurations of this method under different accuracy requirements.

[0075] Accordingly, step S4, which calculates the probability value of the mesh cell belonging to the category corresponding to different Gaussian components using a first-order approximation, includes: Step S4141: Obtain the density parameter covariance matrix corresponding to each grid cell, and the variance of the Gaussian component corresponding to the current category.

[0076] In the implementation of extending statistical features to the covariance matrix, two types of input data with different structural levels are prepared for subsequent uncertainty fusion calculations. The first type of data is the density parameter covariance matrix corresponding to the current grid cell. This covariance matrix is ​​derived from a higher-order statistical feature obtained in step S3 based on the statistical analysis of the density parameters of the current grid cell using the model parameter sample set. Its main diagonal elements represent the variance of the density parameter of the grid cell itself, while the off-diagonal elements describe the covariance relationship between the density parameters of the grid cell and other grid cells, i.e., the degree of correlation between density estimation errors at different spatial locations. In the geophysical inversion process, due to the volume effect of the gravitational field, the density parameters of adjacent grid cells often have a compensatory relationship of mutual increase and decrease when solving. This spatial correlation is quantified and expressed through the off-diagonal elements of the covariance matrix. The second type of data is the variance of the Gaussian component corresponding to the currently examined category. This variance is an inherent parameter in the Gaussian mixture model that describes the natural dispersion of density within this type of geological body. By obtaining the covariance matrix and Gaussian component variance data structures side by side, a data foundation is provided for uncertainty fusion at a more complete statistical structure level in step S4142.

[0077] Step S4142: Fuse the covariance matrix with the variance of the Gaussian components to construct the corrected covariance matrix.

[0078] The covariance matrix of the current grid cell density parameters obtained in step S4141 is fused with the variance of the Gaussian component of the current category to construct a corrected covariance matrix. In specific implementation, this fusion operation can superimpose the variance of the Gaussian component into the structure of the covariance matrix in an appropriate manner. For example, the variance of the Gaussian component can be fused with the corresponding main diagonal elements in the covariance matrix to form the corrected main diagonal terms, while retaining the off-diagonal covariance structure in the original covariance matrix, thus obtaining the complete corrected covariance matrix. This corrected covariance matrix simultaneously accommodates two types of statistical information from different sources: the inherent dispersion of the density distribution of this type of geological body represented by the Gaussian component variance, and the uncertainty of the estimation of density parameters at various spatial locations during the model solution process and their interrelationship patterns represented by the model parameter sample covariance matrix. Compared with the original covariance matrix, the modified covariance matrix not only improves the overall discreteness level of the main diagonal direction, but also more completely inherits the off-diagonal information describing the spatial correlation structure, enabling subsequent probability density calculations to be performed within a more adequate framework of uncertainty description.

[0079] Step S4143: Based on the corrected covariance matrix, calculate the corrected probability density of the grid cell under the current category.

[0080] Using the corrected covariance matrix constructed in step S4142 as the variance parameter, the corrected probability density of the grid cell under the current category is calculated. Due to the introduction of the covariance matrix, the probability density calculation here is no longer limited to the independent uncertainty represented by a single standard deviation. Instead, within the unified framework of the covariance structure, the estimation uncertainty of the current grid cell's density parameter and the error correlation between this location and other spatial locations are incorporated into the construction of the probability density function. In actual geological scenarios, if a grid cell is located deep within the profile and the observation data has weak constraints on this location, the main diagonal terms of its density parameter covariance matrix are often large. Simultaneously, since multiple adjacent grid cells in the deep region often exhibit unidirectional or opposite-directional linkage characteristics during gravity inversion, their off-diagonal terms may also have significant values. This information is incorporated into the probability density calculation through the corrected covariance matrix, allowing the corrected probability density distribution of the grid cell under the current category to more realistically reflect the uncertainty structure remaining after model solving. For regions with strong spatial correlation, using the covariance matrix can more precisely characterize the joint uncertainty caused by the inversion compensation effect between different grid cells compared to using only the standard deviation.

[0081] Step S4144: Normalize the modified probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

[0082] After calculating the corrected probability density for each category based on the corrected covariance matrix, a normalization process is performed to convert the corrected probability density of the current grid cell under all categories into a probability value sequence with a sum of one. This normalization process is operationally consistent with the process described in step S413, using the sum of the corrected probability densities of each category as the normalization denominator to calculate the normalized probability value for each category. The probability value sequence obtained after normalization still represents the degree of confidence that the grid cell belongs to different geological body types within the current Gaussian mixture model framework, after considering the uncertainty of the model solution. However, because the covariance matrix was incorporated into the probability density calculation process in the previous steps, the probability values ​​produced here have integrated the spatial correlation uncertainty structure into the discrimination conclusion. This makes the final probability identification result not only reflect the independent classification tendency of each grid cell in a statistical sense, but also reflect the linkage characteristics of classification uncertainty at spatially adjacent locations. This provides a richer information dimension for subsequent geological interpretation, especially near the boundaries of anomalies, where the introduction of covariance information can provide additional quantitative reference for the reliability of boundary location inference.

[0083] By fusing the covariance matrix of the model parameter samples with the Gaussian component variance and performing probability density calculation and normalization based on the corrected covariance matrix, this implementation extends the uncertainty propagation mechanism from the standard deviation dimension of a single grid cell to the covariance dimension that includes spatial correlation structures. This achieves joint quantitative propagation of the uncertainty in the estimation of density parameters at different spatial locations and their interrelationship patterns during the model solving process. As a result, the final output category probability value can reflect more complete statistical uncertainty structure information, further improving the comprehensiveness and reliability of the geological interpretation of density anomaly identification results.

[0084] Simultaneously, step S4 generates density anomaly identification results, including: Step S421: For each grid cell, obtain the probability value of its belonging to different categories, and compare it with the preset probability threshold corresponding to each category.

[0085] After obtaining the probability values ​​of each grid cell belonging to different Gaussian component categories through the aforementioned steps, the conversion process from probability values ​​to identification results is initiated. For the current grid cell, the probability values ​​of its belonging to each category calculated in step S413 or step S4144 are obtained, and these probability values ​​are compared one by one with the preset probability threshold corresponding to that category. The preset probability threshold is a discrimination threshold set for each category, representing the minimum confidence level required to classify the grid cell into a certain category. In specific geological exploration implementation, the preset probability thresholds for different categories can be set differently according to the risk preference and fault tolerance requirements of the exploration target. For example, for high-value target ore body categories, a relatively low threshold can be set to improve the capture rate and reduce the risk of missed detection; for background surrounding rock categories, a higher threshold can be set to maintain the purity of the identification results and avoid too many non-abnormal areas being mislabeled. The specific value of the threshold can be selected as needed between 0.5 and 0.9. By comparing the probability value with the threshold class by class, the discrimination result of the grid cell in each category is no longer a vague probability number, but is transformed into a clear binary judgment conclusion that either the discrimination condition is met or not.

[0086] Step S422: Simultaneously label the grid cell as all categories that meet the probability threshold conditions to form a multi-category joint identification result. The multi-category joint identification result characterizes the grid cell's probability of belonging to different geological body types.

[0087] After comparing the probability values ​​of each category with their corresponding thresholds, a multi-category joint labeling operation is performed. For the current grid cell, all categories that meet their respective preset probability threshold conditions are simultaneously assigned to that grid cell as its identification label, forming a multi-category joint identification result. This multi-category joint labeling method differs substantially from the conventional practice of single hard classification: single hard classification usually only assigns grid cells to the category with the highest probability value, even if that highest probability value itself is not significantly higher than other categories; while multi-category joint identification allows a grid cell to belong to multiple categories simultaneously under certain conditions. For example, a grid cell in a certain boundary area may be simultaneously labeled as both background rock and high-density ore body. This joint labeling result objectively reflects the actual situation that the density parameter at this location is in the transition range between the two types of geological bodies, and both geological interpretations have a certain degree of credibility. In practical geological interpretation, multi-category joint identification results provide geophysical exploration workers with richer information dimensions. They can identify high-confidence areas with clearly defined categories, marked by only a single category, as well as low-confidence mixed areas where conditions for multiple categories are met. Users can thus form a more comprehensive understanding of the spatial distribution of anomaly boundaries. For grid cells simultaneously marked by multiple categories, they can be prioritized for verification in subsequent borehole drilling or further exploration deployments, allowing for more refined discrimination by integrating other geophysical or geological data.

[0088] It should be noted that the above embodiments are only some specific examples of the technical solution of the present invention and are not intended to limit the scope of protection of the present invention. In practical applications, in addition to Gaussian mixture models, other probability distribution models with multimodal characteristics can also be used to describe the density distribution of underground media, as long as the model can characterize the multimodal distribution features of different types of geological bodies. Regarding model parameter solving, in addition to the Hamiltonian Monte Carlo sampling method and Langevin diffusion sampling method mentioned above, any sampling method that can efficiently explore the high-dimensional parameter space based on the gradient information of the objective function can be used, as long as it can obtain a sample set that reflects the statistical distribution characteristics of the model parameters. In the probability calculation process, different forms of uncertainty fusion methods can also be used, as long as they can achieve the quantitative transfer of model solution uncertainty to anomaly identification results. Regarding the generation of anomaly identification results, in addition to the maximum probability criterion, probability threshold segmentation, and multi-class joint discrimination, other criteria can be used to extract results according to the actual needs of the geological task. These alternative methods are all extensions of the core inventive concept of the present invention.

[0089] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor and a memory connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, cause the at least one processor to perform the aforementioned density anomaly probabilistic identification method based on a Gaussian mixture model.

[0090] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described method for probabilistic identification of density anomalies based on a Gaussian mixture model.

[0091] The embodiments of this invention aim to protect a probability identification method for density anomalies based on a Gaussian mixture model, which has the following effects: 1. By employing a gradient-based sampling method to efficiently solve the model parameters of a Gaussian mixture model, and using the gradient of the objective function to guide the sampling path, a rapid exploration of the model parameter distribution is achieved in the high-dimensional parameter space of gravity inversion. This significantly improves the sampling convergence speed and reduces the correlation of sample sequences. Thus, while ensuring the representativeness of the model parameter sample set, it effectively solves the problem of low computational efficiency of traditional random sampling methods in solving high-dimensional density parameters, providing a reliable data foundation for subsequent uncertainty analysis and probability identification. 2. By explicitly introducing the statistical characteristics, especially the standard deviation, obtained from solving the model parameters during the probability calculation process into the probability density function of the Gaussian component in a first-order approximation manner, the uncertainty generated in the model solving stage is quantitatively transferred to the anomaly identification stage. This results in the final output no longer being a single deterministic classification result, but rather a specific probability value for each spatial grid cell to belong to different geological body categories. This solves the problem that existing technologies cannot directly obtain the probability distribution of anomalies based on model statistical characteristics and that the identification results lack quantitative uncertainty expression, significantly improving the reliability of the geological interpretation of density anomaly identification results. 3. By establishing a correspondence between each Gaussian component in the Gaussian mixture model and the geological prior information of the target area, and determining the number of Gaussian components and configuring their initial parameters according to the background medium type and potential anomaly type, multimodal density parameter modeling is no longer a purely mathematical distribution fitting, but a physical modeling process with clear geological meaning. The probability values ​​of each grid cell belonging to the corresponding category of different Gaussian components can be directly interpreted as the probability of belonging to different geological body types. This solves the problem of the disconnect between existing statistical modeling methods and geological cognition, and improves the correspondence and interpretability between the identification results and geological reality.

[0092] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0093] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0094] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0095] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0096] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A method for probabilistic identification of density anomalies based on Gaussian mixture models, characterized in that, Includes the following steps: S1: Obtain gravity anomaly observation data of the target area, and discretize the underground space of the target area into multiple grid cells to construct a forward model of the relationship between the gravity anomaly observation data and the density parameters of each grid cell; S2: Construct a Gaussian mixture model to describe the distribution of density parameters of each grid cell, and solve the Gaussian mixture model based on the forward model and the gravity anomaly observation data to obtain a model parameter sample set consisting of multiple sets of model parameter values; wherein, the Gaussian mixture model is composed of a weighted combination of multiple Gaussian components, and the model parameter values ​​include the mean, variance and weight of each Gaussian component; Step S2 involves solving the Gaussian mixture model based on the forward model and the gravity anomaly observation data to obtain a model parameter sample set consisting of multiple sets of model parameter values, including: Step S21: Construct an objective function that includes data fitting terms and model constraint terms, wherein the data fitting terms are used to measure the degree of agreement between the model forward response and the gravity anomaly observation data; Step S22: Using a gradient-based sampling method, iterative sampling is performed under the guidance of the gradient of the objective function to generate the model parameter sample set; S3: Based on the model parameter sample set, calculate the density parameter statistical characteristics corresponding to each grid cell, wherein the statistical characteristics include at least the mean and standard deviation; S4: Based on the mean and standard deviation of the density parameters corresponding to each grid cell, as well as the mean, variance and weight of each Gaussian component, the probability value of the grid cell belonging to the category corresponding to different Gaussian components is calculated by first-order approximation to generate the density anomaly identification result. The first-order approximation refers to applying a moment-based first-order approximation to the posterior and prior distributions to realize the transfer of uncertainty in model solution to probability calculation results; Step S4, which involves calculating the probability value of a grid cell belonging to a category corresponding to different Gaussian components using a first-order approximation based on the mean and standard deviation of the density parameters corresponding to each grid cell, as well as the mean, variance, and weight of each Gaussian component, includes: Step S411: Based on the mean, variance and weight of each Gaussian component, calculate the initial probability density of the grid cell under each category; Step S412: Introduce the standard deviation into the variance parameter of the initial probability density, and obtain the corrected probability density by superimposing the variances; Step S413: Normalize the corrected probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

2. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 1, characterized in that, The gradient-based sampling method mentioned in step S22 is the Hamiltonian Monte Carlo sampling method; Step S22 employs a gradient-based sampling method, performing iterative sampling guided by the gradient of the objective function to generate the model parameter sample set, including: Step S221a: Introduce an auxiliary momentum variable to construct a Hamiltonian dynamical system with the objective function as the potential energy function; Step S222a: Simulate the motion trajectory of the Hamiltonian dynamic system using a preset numerical solution method to generate candidate samples; Step S223a: Based on the preset acceptance criteria, the candidate samples are screened, and the samples that meet the conditions are retained to form the model parameter sample set.

3. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 1, characterized in that, Step S2, which involves constructing a Gaussian mixture model to describe the distribution of density parameters for each grid cell, includes: Step S201: Obtain the geological prior information of the target area, wherein the geological prior information includes at least the background medium type and the potential anomaly type; Step S202: Determine the number of Gaussian components in the Gaussian mixture model based on the geological prior information, and configure an initial mean and initial variance for each Gaussian component to characterize the density distribution characteristics of the corresponding geological body type; Step S203: Based on the configured Gaussian components, construct the complete Gaussian mixture model.

4. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 1, characterized in that, Step S412, which involves incorporating the standard deviation into the variance parameter of the initial probability density and obtaining the corrected probability density through variance superposition, includes: Step S4121: Obtain the standard deviation of the density parameter corresponding to the grid cell and the variance of the Gaussian component corresponding to the current category; Step S4122: Perform a sum of squares operation on the standard deviation and the standard deviation of the Gaussian component, and use the result of the sum of squares operation as the corrected variance parameter to construct the corrected probability density.

5. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 1, characterized in that, The normalization process for the corrected probability density corresponding to each category in step S413 includes: Step S4131: Calculate the sum of the corrected probability densities of the grid cell under all categories, and use it as the normalized denominator; Step S4132: Divide the corrected probability density corresponding to each category by the normalized denominator to obtain the probability value of the grid cell belonging to that category.

6. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 4, characterized in that, The statistical features also include the covariance matrix; Step S4, which involves calculating the probability value of the mesh cell belonging to the category corresponding to different Gaussian components using a first-order approximation, includes: Step S4141: Obtain the density parameter covariance matrix corresponding to each grid cell, and the variance of the Gaussian component corresponding to the current category; Step S4142: Fuse the covariance matrix with the variance of the Gaussian components to construct a corrected covariance matrix; Step S4143: Based on the corrected covariance matrix, calculate the corrected probability density of the grid cell under the current category; Step S4144: Normalize the modified probability density corresponding to each category to obtain the probability value of the grid cell belonging to different categories.

7. The density anomaly probabilistic identification method based on Gaussian mixture model according to claim 1, characterized in that, The generation of density anomaly identification results in step S4 includes: Step S421: For each grid cell, obtain its probability value belonging to different categories, and compare it with the preset probability threshold corresponding to each category respectively; Step S422: Simultaneously label the grid cell as all categories that meet the probability threshold conditions to form a multi-category joint identification result, which characterizes the grid cell's probability of belonging to different geological body types.

8. An electronic device, characterized in that, include: At least one processor; And a memory connected to at least one processor; wherein the memory stores instructions executable by at least one processor, the instructions being executed by at least one processor to cause at least one processor to perform the density anomaly probabilistic identification method based on a Gaussian mixture model as described in any one of claims 1-7.

9. A computer-readable storage medium, characterized in that, It stores computer instructions that, when executed by a processor, implement the density anomaly probabilistic identification method based on a Gaussian mixture model as described in any one of claims 1-7.

Images