Gravity exploration ultrahigh-precision inversion method based on multi-stage simulated annealing

Abstract

The invention belongs to the technical field of geophysical exploration, and particularly provides a gravity exploration ultrahigh-precision inversion method based on multi-stage simulated annealing, and the method specifically comprises the steps: 1, carrying out the parameterization of an exploration region, and building an initial model; 2, performing gravity forward modeling simulation calculation; 3, performing multi-stage simulated annealing inversion optimization; and step 4, performing inversion result evaluation and quality analysis. Compared with the prior art, the method has the advantages that the problem of gravity inversion is solved, and the method has serious multiplicity of solutions; the abnormal body position identification precision is seriously insufficient; an optimization algorithm is easy to fall into a local optimal solution; and a contradiction which is difficult to reconcile exists between the calculation efficiency and the inversion precision. The calculation efficiency is improved by more than 30%; and the average number of iterations is reduced by 25%. The global search capability is enhanced, and the probability of finding a global optimal solution reaches 95% or above. A complex neural network is not needed, the method is simplified, and high precision is still kept.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration technology, and specifically provides an ultra-high precision inversion method for gravity exploration based on multi-stage simulated annealing. Background Technology

[0002] Gravity exploration, as an important branch of geophysical exploration, infers the density distribution differences of subsurface media by measuring the spatial variations of the Earth's gravity field, thereby identifying subsurface structures and mineral occurrence states. Traditional gravity inversion methods mainly include linear inversion methods, least squares inversion methods, stochastic search methods, and gradient optimization-based inversion methods.

[0003] Existing technologies have many technical defects and limitations in practical engineering applications: 1. Gravity inversion problems suffer from severe ambiguity. Due to the equivalent source characteristics of the gravity field, different density distribution models may produce extremely similar surface gravity responses, leading to significant uncertainties in the inversion results. Especially under complex geological conditions, traditional methods struggle to accurately distinguish the spatial locations and physical properties of multiple anomalies.

[0004] 2. Insufficient accuracy in identifying the location of anomalies. Conventional gravity inversion methods have limited ability to identify the boundaries of anomalies, with location errors typically reaching tens of meters or even greater, failing to meet the demands of modern high-precision exploration. In mineral resource exploration, such location errors can lead to drilling positioning deviations, resulting in significant economic losses.

[0005] 3. Optimization algorithms are prone to getting trapped in local optima. Traditional optimization methods based on gradient descent are highly susceptible to getting trapped in local minima in complex nonlinear inversion problems, failing to guarantee a globally optimal solution. Although simulated annealing theoretically possesses global search capabilities, in practical gravity inversion applications, due to unreasonable parameter settings and imperfect convergence mechanisms, it is still difficult to effectively avoid local optima problems.

[0006] 4. There is an irreconcilable contradiction between computational efficiency and inversion accuracy. High-precision inversion typically requires significant computational resources and time, while fast algorithms often sacrifice accuracy. This contradiction is particularly pronounced in large-scale exploration projects, severely limiting the practical application of gravity exploration technology.

[0007] Therefore, there is an urgent need to develop a gravity exploration inversion method that can significantly improve the accuracy of anomaly location identification while ensuring computational efficiency. Summary of the Invention

[0008] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows: The ultra-high precision inversion method for gravity exploration based on multi-stage simulated annealing includes the following steps: Step 1: Parameter initialization, parameterization of the exploration area and establishment of the initial model; Discretize the exploration area into a grid and establish a parameterized density model; Step 2: Forward gravity simulation calculation; Based on the established density model, calculate the gravity anomaly response at each measuring point on the surface; Step 3: Multi-stage simulated annealing inversion optimization; Set optimization strategies for multiple stages, and adopt different perturbation mechanisms and convergence criteria for each stage; The inversion process is as follows: Stage 1: Define the stage transition mechanism using a piecewise function; Stage 2: Define the model parameter vector; Stage 3: Adaptive temperature update; Stage 4: Design of the comprehensive objective function; Stage 5: Acceptance probability criterion; Stage 6: Optimization mechanism for the position of the anomaly body; Stage 7: Optimization update of the density parameter; Step 4: Evaluate the inversion results and analyze the quality, and establish a complete quality evaluation system.

[0009] Furthermore, in Stage 1 of Step 3, the stage transition mechanism is defined using the piecewise function S(k);​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​​k stage This represents the stage-related model update amount in the k-th iteration; The stage-dependent perturbation operator is defined as: ; In the formula, T is a temperature parameter; relatively speaking, T k The temperature at the k-th iteration; Θ represents the set of control parameters; P is the update policy function; m k This represents the current density model at the k-th iteration.

[0011] Furthermore, in stage three of step three, the temperature update employs a two-factor control mechanism based on stage and location error: ; The cooling coefficient function is specifically defined as follows: ; In the formula, σ k pos The position error at the k-th iteration reflects the degree of deviation between the current model and the actual anomaly position.

[0012] Furthermore, in stage five of step three, a modified Metropolis criterion is adopted: ; In the formula, P accept For the probability of acceptance; ΔE represents the energy change; T represents the current temperature; ; ζ(m k ) is the adaptive adjustment factor; When the objective function is improved (ΔE<0), the new solution is always accepted; When the objective function deteriorates (ΔE≥0), the probability is exp(-ΔE / (T·ζ(m)). k Accepting the new interpretation; The acceptance probability is dynamically adjusted based on the current model quality.

[0013] Furthermore, in step three, the inversion process terminates and outputs the final result when any of the following conditions are met: the total number of iterations exceeds the limit; stage four is fully converged; there is no improvement for a long time; or the temperature is too low.

[0014] The beneficial effects of using this invention are: Computational efficiency is improved by more than 30%; The average number of iterations was reduced by 25%; The global search capability is enhanced, and the probability of finding the global optimal solution reaches over 95%. It eliminates the need for complex neural networks, simplifying the method while maintaining high accuracy. Attached Figure Description

[0015] Figure 1 This is a flowchart of the present invention; Figure 2 This is the result of the actual model; Figure 3 This is the inversion result of the present invention; Figure 4 This is a schematic diagram of the absolute density error of the present invention; Figure 5 The diagram shows the convergence curve and error analysis of this invention. Detailed Implementation

[0016] The present invention will now be described in detail with reference to the accompanying drawings.

[0017] Reference Figure 1 A high-precision inversion method for gravity exploration based on multi-stage simulated annealing. The specific steps include: Step 1: Parameter initialization, parameterization of the exploration area and establishment of the initial model; The exploration area is discretized into a grid.

[0018] Establish a parameterized density model: ; The grid size parameter is defined as follows: ; ; G is the set of matrices for the density model, used to describe the discretized grid of the exploration area; n is a vector representing the number of grid cells; d is the size of the mesh cell; Δx is the horizontal dimension of each grid cell; Δz is the vertical dimension of each grid cell.

[0019] Define the density fluctuation range ρ in the physical property parameter space. range : ; In the formula, ρ0 is the background density; Specifically, the horizontal dimension Δx of the grid cell should be less than or equal to the distance between the field measurement points; The vertical dimension Δz of the grid cell matches the detection depth.

[0020] Step 2: Gravity forward simulation calculation; Based on the established density model, calculate the gravity anomaly responses at each measurement point on the surface.

[0021] The gravity forward calculation uses the following formula: ; In the formula, g is the gravity anomaly response value; ρ ij represents the density value at the grid (i,j); The physical constants G and the conversion factor K conv are respectively: ; The preprocessing of the observation data includes noise filtering and trend term elimination: ; In the formula, g obs is the observed gravity data; g raw is the original gravity data; g trend is the trend term gravity data; ε is the noise or error term; Step 3: Multi-stage simulated annealing inversion optimization; Set the optimization strategies for multiple stages, and each stage adopts different perturbation mechanisms and convergence criteria.

[0022] Stage 1: Define the stage transition mechanism using the piecewise function S(k); ; In the formula, k is the iteration number; N1, N2, and N3 are respectively the boundary values of each stage; Determine the current optimization stage S(k) according to the k value, When k ≤ N1, it is in the first stage; When N1 < k ≤ N2, it is in the second stage; When N2 < k ≤ N3, it is in the third stage; When k > N3, it is in the fourth stage.

[0023] Among them, N1, N2, and N3 are jointly determined by the iteration number and the convergence state: ; In the formula, N totalThis represents the total number of iterations.

[0024] Phase 2: m is the model parameter vector. In each iteration, the update formula for m is: ; In the formula, Δm k stage This represents the stage-related model update amount in the k-th iteration; The stage-dependent perturbation operator is defined as: ; In the formula, T is a temperature parameter; relatively speaking, T k The temperature at the k-th iteration; Θ represents the set of control parameters; P is the update policy function; m k This represents the current density model at the k-th iteration. in, P coarse It is a coarse disturbance used in the first stage, characterized by large disturbance amplitude and strong exploratory nature; Parameter configuration Θ1: large step size (0.3), strong temperature coupling (1.5), high random perturbation (0.4); Behavioral rules: uniform random perturbation, high acceptance rate of inferior solutions, no strict convergence conditions; In the early stages of the inversion process, a global "coarse search" is performed to quickly find potential regions that may contain the global optimum, thus avoiding getting trapped in local optima too early. P fine This is a fine perturbation used in the second stage; the perturbation amplitude is less than P. coarse ; Parameter configuration Θ2: medium step size (0.15), medium temperature coupling (1.0), medium disturbance (0.2); Behavioral rules: Gaussian distribution perturbation, gradient-assisted direction, and relaxed convergence conditions; Within the roughly located potential area, a more refined search is conducted, and convergence begins. P micro This is a micro-perturbation, used in the third stage, with an amplitude less than P. fine ; Parameter configuration Θ3: small step size (0.05), weak temperature coupling (0.5), small perturbation (0.05); Behavioral rules: directional perturbation, conjugate direction, strict convergence condition; In regions that are close to the optimal solution, "fine-tuning" is performed to precisely optimize the model parameters; P ultraThis is an ultra-micro perturbation; used in the fourth stage, the perturbation amplitude is less than P. micro ; Parameter configuration Θ4: microstep size (0.01), extremely weak temperature coupling (0.1), minimal perturbation (0.01); Behavioral rules: Newtonian direction, extremely strict convergence conditions, and acceptance of almost only improved solutions; In the final stage of the inversion, a final refinement is performed to obtain the most accurate inversion results.

[0025] This is the switch to implement a multi-stage strategy. Each stage automatically switches based on the temperature threshold, number of iterations, and convergence state. As the iteration progresses, the algorithm switches from one stage to the next, forming a complete optimization pipeline.

[0026] This design enables a smooth transition from global search to local fine-grained optimization and is a core component of multi-stage optimization strategies.

[0027] Phase 3: Adaptive Temperature Update; Temperature updates employ a two-factor control mechanism based on stage and location errors: ; The cooling coefficient function is specifically defined as follows: ; In the formula, σ k pos The position error at the k-th iteration reflects the degree of deviation between the current model and the actual anomaly position.

[0028] Phase Four: Design of the Comprehensive Objective Function; The objective function of the inversion problem is in a multi-constraint weighted form: ; In the formula, E is the weighting coefficient; w is the weighting coefficient, which determines the importance or influence of this term in the overall objective function.

[0029] The specific calculation formulas for each component are as follows: Data fitting terms: ; Position error term: ; Penalty for the number of abnormal entities: ; Model smoothness constraint terms: ; Prior information constraints: ; Weighting coefficients are adaptively adjusted: ; Phase 5: Acceptance Probability Criterion; Adopting the improved Metropolis criterion: ; In the formula, P accept For the probability of acceptance; ΔE represents the energy change; T represents the current temperature; ; ζ(m k ) is the adaptive adjustment factor; When the objective function is improved (ΔE<0), the new solution is always accepted; When the objective function deteriorates (ΔE≥0), the probability is exp(-ΔE / (T·ζ(m)). k Accepting the new interpretation; The acceptance probability is dynamically adjusted based on the current model quality. ; Phase Six: Anomaly Location Optimization Mechanism; For each detected anomaly, position optimization employs a stage-adaptive step size r: ; The displacement is specifically defined as: ; The position perturbation amount is determined based on the current optimization stage S(k). The perturbation amplitude is controlled by a different random integer range for each stage. From stage one to stage four, the perturbation amplitude gradually decreases, realizing the transition from coarse adjustment to fine adjustment. Phase 7: Density parameter optimization and update; The update of the anomalous volume density parameter employs a temperature-dependent adaptive mechanism: ; The formula for calculating the change in density is: ; The magnitude of the disturbance related to the stage: ; The inversion process terminates and outputs the final result when any of the following conditions are met: 1. Total number of iterations exceeded: k≥N max ; Normally, N is taken. max=N total This ensures that the algorithm fully executes the four preset stages; N can also be set max >N total Provide additional convergence buffer for Phase 4; 2. Stage 4 complete convergence (must satisfy simultaneously): The objective function value changes less than ε in M ​​consecutive iterations. obj ; All parameter changes are less than ε. param ; Numerical gradient norm less than ε grad ; At least N have been executed in Phase 4. min The next iteration; 3. No improvement over a long period: P consecutive iterations failed to improve the optimal objective function value; 4. Temperature too low: System temperature T k <T min This indicates that the ability to continue optimization has been lost.

[0030] Step 4: Evaluation and quality analysis of the inversion results, and establishment of a complete quality evaluation system; Reference Figure 4 and Figure 5 Data fit quality assessment: ; The fitting correlation coefficient R: ; In the formula, g obs The actual observed gravity anomaly value at the measuring point; g inv The gravity anomaly value at the measuring point is predicted (calculated by forward modeling) by the density model obtained from the inversion. Location recognition accuracy assessment: ; Model complexity assessment: ; Uncertainty analysis: ; Reference Figure 2 and Figure 3 Taking a 500×300 exploration case as an example, that is ; The comparison between the two figures shows that the results after inversion processing highly overlap with the actual results, thus verifying that the inversion results have high accuracy.

[0031] The above content is only a preferred embodiment of the present invention. For those skilled in the art, many changes can be made in the specific implementation and application scope based on the concept of the present invention. As long as these changes do not depart from the concept of the present invention, they all fall within the protection scope of the present invention.

Claims

1. A high-precision inversion method for gravity exploration based on multi-stage simulated annealing, characterized in that: The specific steps are as follows: Step 1: Parameter initialization, parameterization of the exploration area, and establishment of an initial model; Perform discretized grid division on the exploration area and establish a parameterized density model; Step 2: Forward gravity simulation calculation; Based on the established density model, calculate the gravity anomaly responses of each measuring point on the surface; Step 3: Multi-stage simulated annealing inversion optimization; Set optimization strategies for multiple stages, and adopt different perturbation mechanisms and convergence criteria for each stage; The inversion process is as follows: Stage 1: Define the stage transition mechanism using a piecewise function; Stage 2: Define the model parameter vector; Stage 3: Adaptive temperature update; Stage 4: Design of the comprehensive objective function; Stage 5: Acceptance probability criterion; Stage 6: Optimization mechanism for the position of the anomaly body; Stage 7: Optimization update of density parameters; Step 4: Evaluate the inversion results and analyze the quality, and establish a complete quality evaluation system.

2. The ultra-high-precision inversion method for gravity exploration based on multi-stage simulated annealing according to claim 1, wherein: In stage 1 of step 3, the stage transition mechanism is defined using the piecewise function S(k); ； In the formula, k is the iteration number; N total This represents the total number of iterations. N1, N2, and N3 are the boundary values of each stage respectively, and N1, N2, and N3 are both determined by the iteration number and the convergence state; Determine the current optimization stage S(k) according to the value of k, When k ≤ N1, it is in the first stage; When N1 < k ≤ N2, it is in the second stage; When N2 < k ≤ N3, it is in the third stage; When k > N3, it is in the fourth stage.

3. The ultra-high-precision inversion method for gravity exploration based on multi-stage simulated annealing according to claim 2, wherein: 。 4. The ultra-high-precision inversion method for gravity exploration based on multi-stage simulated annealing according to claim 2, wherein: [[ID= ​ ； ​ Δm k stage This represents the stage-related model update amount in the k-th iteration; ​ ； ​ T is a temperature parameter; relatively speaking, T k The temperature at the k-th iteration; ​ ​ m k This represents the current density model at the k-th iteration. ​ ​ ； ​ ； ​ σ k pos The position error at the k-th iteration reflects the degree of deviation between the current model and the actual anomaly position. ​ ​ ； ​ P accept For the probability of acceptance; ​ ​ ； ζ(m k ) is the adaptive adjustment factor; ​ When the objective function deteriorates (ΔE≥0), the probability is exp(-ΔE / (T·ζ(m)). k Accepting the new interpretation; ​ ​ ​

Images