Method for iterative terrain compensation of middle and far zone gravity anomalies based on dynamic density model

Abstract

The present application relates to the technical field of gravity measurement, in particular to a method for long-mid zone terrain iterative compensation of gravity anomaly based on a dynamic apparent density model, comprising: performing multi-scale decomposition on initial Bouguer gravity anomaly and DEM data respectively to obtain gravity anomaly components and terrain components corresponding to different wavelengths; screening at least one type of target wavelength according to the correlation coefficient between the gravity anomaly components and terrain components under the same wavelength; reconstructing the gravity anomaly components under all target wavelengths to obtain a pure terrain gravity anomaly field; and performing multi-round iterative compensation correction based on the pure terrain gravity anomaly field until the termination condition is met, and outputting the final updated Bouguer gravity anomaly. The present application realizes high-precision variable-density terrain compensation by constructing a complete technical link of adaptive anomaly separation, multi-constraint density inversion, partition coupling compensation and iterative optimization closed loop.

Description

Technical Field

[0001] This invention relates to the field of gravity measurement technology, and more specifically to an iterative compensation method for distant terrain in gravity anomalies based on a dynamic apparent density model. Background Technology

[0002] Gravity exploration is one of the core technologies of geophysical exploration, widely used in mineral resource exploration, geological structure surveys, oil and gas reservoir detection, and regional geological mapping. With the rapid improvement of gravity instrument observation accuracy and surveying technology, gravity exploration has entered a high-precision era. However, technological breakthroughs in terrain correction errors have lagged behind other aspects. Among them, terrain correction in the mid-to-far range (2km to 166.7km) has become a core bottleneck restricting the reliability of high-precision gravity exploration results due to its wide coverage and complex influencing factors. Dynamic apparent density models, as a technical means that can accurately characterize the spatiotemporal distribution of equivalent density of surface and near-surface materials, have a natural advantage in solving the problem of variable density terrain compensation. However, current technologies have not yet effectively combined them with layered iterative error elimination mechanisms.

[0003] Existing problems: Traditional mid-to-long-range terrain correction techniques generally use 2.67 g / cm³. 3 The calculation method of the unified stratigraphic density assumption (grams per cubic centimeter) combined with the low-resolution digital elevation grid has led to existing improvement technologies focusing on marginal optimization directions such as earth curvature correction and DEM (Digital Elevation Model) data benchmark conversion. Although some three-dimensional apparent density inversion technologies can reconstruct the subsurface density field, they do not introduce a dynamic apparent density model to perform layered iterative processing of topographic errors in the mid-to-far area. Instead, they uniformly process the complex sources of topographic errors, resulting in mutual cancellation and overlap during error elimination. The correction accuracy is difficult to meet the engineering application requirements of current high-precision gravity exploration. Summary of the Invention

[0004] This invention provides an iterative compensation method for distant terrain in gravity anomalies based on a dynamic apparent density model, in order to solve existing problems.

[0005] The gravity anomaly distant region terrain iterative compensation method based on dynamic apparent density model of the present invention adopts the following technical solution: One embodiment of the present invention provides an iterative compensation method for distant terrain in gravity anomalies based on a dynamic apparent density model. The method includes the following steps: Obtain the initial Bouguer gravity anomaly and DEM data corresponding to the target area; Multi-scale decomposition is performed on the initial Bouguer gravity anomaly and DEM data to obtain gravity anomaly components and terrain components corresponding to different wavelengths. At least one type of target wavelength is selected based on the cross-correlation coefficient between the anomalous gravity component and the terrain component at the same wavelength. The anomalous gravity components at all target wavelengths are reconstructed to obtain a pure terrain gravity anomaly field. Based on the pure terrain gravity anomaly field, the terrain compensation value of the first iteration in the first round is calculated to compensate and correct the initial Bouguer gravity anomaly, the Bouguer gravity anomaly of the first iteration in the first round is determined, and the first round of iterative compensation and correction is carried out until the convergence condition is met, and the updated Bouguer gravity anomaly of the first round is output. Based on the pure terrain gravity anomaly field after the first round of updated Bouguer gravity anomaly decomposition and reconstruction, multiple rounds of iterative compensation and correction are performed until the termination condition is met, and the final updated Bouguer gravity anomaly is output.

[0006] Furthermore, the specific steps for screening out at least one type of target wavelength are as follows: Calculate the cross-correlation coefficient between the anomalous gravity component and the topographic component at the same wavelength, and denot it as the gravity-topographic cross-correlation coefficient at each wavelength. If the gravity-terrain cross-correlation coefficient at each wavelength is greater than or equal to the preset cross-correlation coefficient threshold, then each wavelength is recorded as the target wavelength.

[0007] Furthermore, the specific steps involved in determining the Bouguer gravity anomaly in the first round of the first iteration are as follows: Inputting the pure topographic gravity anomaly field into the improved The model is inverted to obtain the overall average apparent density value corresponding to the target area. Combined with DEM data, the terrain compensation value of the first iteration of the first round is calculated using the terrain correction formula that takes into account the curvature of the earth. The first iteration of the first round of Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first round of the first iteration onto all elements in the initial Bouguer gravity anomaly.

[0008] Furthermore, the specific steps involved in the first round of updating the Bouguer gravity anomaly are as follows: If the relative deviation between the terrain compensation value of the first iteration of the first round and the preset initial density compensation value is less than or equal to the preset first round convergence deviation, then the Bouguer gravity anomaly of the first iteration of the first round is recorded as the first round updated Bouguer gravity anomaly. If the relative deviation between the terrain compensation value in the first iteration of the first round and the preset initial density compensation value is greater than the preset first round convergence deviation, then this process is repeated iteratively to determine the updated Bouguer gravity anomaly in the first round.

[0009] Furthermore, if the relative deviation between the terrain compensation value in the first iteration of the first round and the preset initial density compensation value is greater than the preset first round convergence deviation, then the calculation is repeated iteratively to determine the updated Bouguer gravity anomaly in the first round. The specific steps include the following: The terrain compensation value and the corresponding Bouguer gravity anomaly are solved successively in the first iteration until the first iteration... The terrain compensation value in the next iteration is the same as that in the first round. The iteration stops when the relative deviation of the terrain compensation value in the next iteration is less than or equal to the preset first-round convergence deviation, and the first-round iteration is stopped after the convergence condition is met. The Bouguer gravity anomaly in the next iteration is denoted as the first round updated Bouguer gravity anomaly. If the preset convergence deviation requirement for the first round is still not met after the first round of consecutive iterations reaches the preset convergence number, the iteration is forcibly terminated, and the Bouguer gravity anomaly of the last iteration of the first round is recorded as the first round updated Bouguer gravity anomaly.

[0010] Furthermore, the output ultimately updates the Bouguer gravity anomaly, including the following specific steps: The pure terrain gravity anomaly field, obtained from the first round of Bouguer gravity anomaly decomposition and reconstruction, is input into the improved... The model is inverted to obtain a two-dimensional grid-by-grid apparent density field with a preset grid size resolution. Combined with DEM data, a three-dimensional apparent density field is constructed. Based on the three-dimensional apparent density field, a terrain correction formula considering the curvature of the earth is adopted. The terrain compensation value of the first iteration of the second round is calculated by integrating grid-by-grid unit. The updated Bouguer gravity anomaly of the first round is compensated and corrected. The Bouguer gravity anomaly of the first iteration of the second round is determined. After the second round of iteration is completed and converged, the updated Bouguer gravity anomaly of the second round is output. Based on the pure terrain gravity anomaly field after the decomposition and reconstruction of the second round of updated Bouguer gravity anomaly, the terrain compensation value of the first iteration of the third round is calculated to compensate and correct the second round of updated Bouguer gravity anomaly, and the first iteration of the third round of Bouguer gravity anomaly is determined. After the third round of iteration is completed and converged, the final updated Bouguer gravity anomaly is output.

[0011] Furthermore, the specific steps involved in determining the Bouguer gravity anomaly in the first iteration of the second round are as follows: The first iteration of the second round of Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first iteration of the second round on all elements of the first round of updated Bouguer gravity anomaly.

[0012] Furthermore, the specific steps involved in outputting the second round of updates to the Bouguer gravity anomaly are as follows: The terrain compensation value corresponding to the first round of updating the Bougu gravity anomaly will be obtained and recorded as the first round terrain compensation value. If the relative deviation between the terrain compensation value of the first iteration of the second round and the terrain compensation value of the first round is less than or equal to the preset convergence deviation of the second round, then the Bouguer gravity anomaly of the first iteration of the second round is recorded as the updated Bouguer gravity anomaly of the second round. If the relative deviation between the terrain compensation value in the first iteration of the second round and the terrain compensation value in the first round is greater than the preset convergence deviation for the second round, then this process is repeated iteratively to solve for the terrain compensation value and the corresponding Bouguer gravity anomaly in each iteration of the second round, until the second round... The terrain compensation value of the second iteration and the second round The iteration stops when the relative deviation of the terrain compensation value in the second iteration is less than or equal to the preset second-round convergence deviation, and the second-round iteration is stopped after the convergence condition is met. The Bouguer gravity anomaly in the next iteration is denoted as the second-round updated Bouguer gravity anomaly. If the second round of consecutive iterations reaches the preset convergence number but still fails to meet the preset second round convergence deviation requirement, the iteration is forcibly terminated, and the Bouguer gravity anomaly of the last iteration of the second round is recorded as the second round updated Bouguer gravity anomaly.

[0013] Furthermore, the specific steps involved in determining the Bouguer gravity anomaly in the first iteration of the third round are as follows: The pure terrain gravity anomaly field, after the second round of updated Bouguer gravity anomaly decomposition and reconstruction, is input into the improved... The model, based on the two-dimensional grid-by-grid apparent density field with a preset grid size resolution obtained by inversion, is combined with DEM data to correct and expand to obtain the final high-precision three-dimensional apparent density field; based on the final high-precision three-dimensional apparent density field, the terrain correction formula considering the curvature of the earth is adopted, and the terrain compensation value of the first iteration of the third round is calculated by grid-by-grid integration. The third round of first iteration Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first iteration on all elements in the second round of updated Bouguer gravity anomaly.

[0014] Furthermore, the output ultimately updates the Bouguer gravity anomaly, including the following specific steps: The terrain compensation value corresponding to the second round of update of the Bougu gravity anomaly will be obtained and recorded as the second round terrain compensation value. If the relative deviation between the terrain compensation value of the first iteration of the third round and the terrain compensation value of the second round is less than or equal to the preset convergence deviation of the third round, then the Bouguer gravity anomaly of the first iteration of the third round is recorded as the final updated Bouguer gravity anomaly. If the relative deviation between the terrain compensation value in the first iteration of the third round and the terrain compensation value in the second round is greater than the preset convergence deviation for the third round, then this process is repeated iteratively to solve for the terrain compensation value and the corresponding Bouguer gravity anomaly in each iteration of the third round, until the third round... The terrain compensation value of the second iteration and the third round The iteration stops when the relative deviation of the terrain compensation value in the third iteration is less than or equal to the preset convergence deviation in the third round, and the convergence condition is met. The next iteration of the Bouguer gravity anomaly is denoted as the final updated Bouguer gravity anomaly. If the third round of consecutive iterations reaches the preset convergence number but still fails to meet the preset convergence deviation requirement, the iteration is forcibly terminated, and the Bouguer gravity anomaly in the last iteration of the third round is recorded as the final updated Bouguer gravity anomaly.

[0015] The beneficial effects of the technical solution of the present invention are: In this embodiment of the invention, the initial Bouguer gravity anomaly and DEM data are decomposed into multiple scales to obtain gravity anomaly components and terrain components corresponding to different wavelengths. Based on the magnitude of the cross-correlation coefficient between the anomalous gravity components and terrain components at the same wavelength, at least one type of target wavelength is selected. The anomalous gravity components at all target wavelengths are reconstructed to obtain a pure terrain gravity anomaly field. Based on the pure terrain gravity anomaly field, the first iteration terrain compensation value of the first round is calculated to compensate and correct the initial Bouguer gravity anomaly. The first iteration Bouguer gravity anomaly of the first round is determined, and the first round of iterative compensation and correction is carried out until the convergence condition is met. The first round of updated Bouguer gravity anomaly is output, thereby eliminating the largest uniform density assumption error. Based on the pure terrain gravity anomaly field reconstructed from the first round of updated Bouguer gravity anomaly decomposition, multiple rounds of iterative compensation and correction are performed until the termination condition is met, outputting the final updated Bouguer gravity anomaly. After eliminating the uniform density error, the density spatial heterogeneity error becomes apparent. The second round of iteration eliminates this density spatial heterogeneity error, minimizing the residual anomaly separation error. This residual error can only be effectively removed after the first two types of errors are eliminated. The third round of iteration eliminates the residual anomaly separation error. Thus, this invention achieves high-precision variable-density terrain compensation by constructing a complete technical chain of adaptive anomaly separation, multi-constraint density inversion, partitioned coupled compensation, and iterative optimization closed loop. Attached Figure Description

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

[0017] Figure 1 This is a flowchart illustrating the steps of the gravity anomaly distant area terrain iterative compensation method based on a dynamic apparent density model according to the present invention. Figure 2 This is a schematic diagram of a hierarchical, progressive, multi-round terrain compensation iterative architecture. Detailed Implementation

[0018] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the following, in conjunction with the accompanying drawings and preferred embodiments, details the specific implementation, structure, features, and effects of the gravity anomaly distant region terrain iterative compensation method based on a dynamic apparent density model proposed in this invention. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0019] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0020] The following description, in conjunction with the accompanying drawings, details the specific scheme of the gravity anomaly distant terrain iterative compensation method based on a dynamic apparent density model provided by this invention.

[0021] Please see Figure 1 The diagram illustrates a flowchart of an iterative compensation method for distant terrain based on a dynamic apparent density model for gravity anomalies, according to an embodiment of the present invention. The method includes the following steps: Step S001: Obtain the initial Bouguer gravity anomaly and DEM data corresponding to the target area.

[0022] To address the systematic errors caused by the uniform density assumption, the loss of detail due to low-resolution terrain characterization, and the inability to provide feedback optimization in a single calculation in existing mid-to-far region terrain correction technologies, this embodiment provides an iterative compensation method for gravity anomalies in mid-to-far region terrain based on a dynamic apparent density model. By constructing a complete technical chain of adaptive anomaly separation, multi-constraint density inversion, partitioned coupling compensation, and iterative optimization closed loop, high-precision variable density terrain compensation is achieved.

[0023] First, gravity data was collected using a data acquisition unit (high-precision relative gravimeter and dual-frequency GPS receiver) combined with national gravity big data, covering the target area and a surrounding range of 166.7 km, with a gravity measurement point density of ≥1 / 10 km. 2 The observation accuracy is ≤20 microgal. Only the three-dimensional coordinates (longitude, latitude, and elevation) and measured gravity values ​​of the gravity measurement points are retained, and duplicate and abnormal measurement points (gravity values ​​exceeding the mean gravity value of the area by ±3 times the standard deviation) are removed.

[0024] Then, the Global Digital Elevation Model (GDEM) V3 30m resolution digital elevation model (DEM) data was used, with coverage completely consistent with the gravity data and elevation accuracy ≤10m. This resolution can accurately characterize terrain undulations with wavelengths ≥60m, fully meeting the accuracy requirements for terrain compensation in the mid-to-far range, while avoiding the computational explosion caused by excessively high resolution.

[0025] It should be noted that in this embodiment, only borehole density data with a depth ≤ 500m within the target area are collected (density measurement accuracy ≤ 0.02g / cm³). 3 ) and 1:200,000 geological maps. Gravity data preprocessing: Wavelet threshold denoising was used, with the db4 wavelet function selected for 3-level decomposition. The threshold was set to 3 times the noise standard deviation to remove random noise. The signal-to-noise ratio of the denoised data was ≥30dB. Normal gravity values ​​were calculated using the international normal gravity formula with a calculation accuracy ≤0.1 microgal, eliminating gravity variations caused by Earth's rotation and ellipsoidal shape. Initial elevation correction was performed using the Bouguer correction formula, with preset initial density compensation values. It is 2.67 g / cm 3 The calculation accuracy is ≤0.1 microgal, and the initial Bouguer gravity anomaly is obtained. This refers to the standardized gravity data after multiple corrections. DEM data preprocessing: The coordinate system of the DEM data is converted to the WGS84 geodetic coordinate system with a conversion error ≤0.5m, ensuring consistency with the coordinates of the gravity measurement points; Kriging interpolation is used to fill in missing values ​​in the DEM data with an interpolation error ≤5m, ensuring data continuity. Wavelet thresholding denoising, the international normal gravity formula, the Bouguer correction formula, the WGS84 geodetic coordinate system, and Kriging interpolation are all well-known techniques, and their specific methods will not be described here.

[0026] This completes the data collection and standardized preprocessing.

[0027] Step S002: Perform multi-scale decomposition on the initial Bouguer gravity anomaly and DEM data respectively to obtain gravity anomaly components and terrain components corresponding to different wavelengths; based on the magnitude of the cross-correlation coefficient between the anomalous gravity components and terrain components under the same wavelength, select at least one type of target wavelength; reconstruct the anomalous gravity components under all target wavelengths to obtain the pure terrain gravity anomaly field; based on the pure terrain gravity anomaly field, calculate the terrain compensation value of the first iteration of the first round, compensate and correct the initial Bouguer gravity anomaly, determine the Bouguer gravity anomaly of the first iteration of the first round, carry out the first round of iterative compensation and correction until the convergence condition is met, and output the updated Bouguer gravity anomaly of the first round.

[0028] It should be noted that: furthermore, a closed-loop iterative process for eliminating orientation errors in the mid-to-long-range region is implemented. The three types of errors in mid-to-long-range terrain correction have clear hierarchical order and causal relationship: the uniform density assumption error has the largest magnitude and is the primary target for elimination; the density spatial heterogeneity error is the next largest and becomes apparent after eliminating the uniform density error; the anomaly separation residual error has the smallest magnitude and can only be effectively removed after the first two types of errors are eliminated. Traditional single-calculation treats the three types of errors together, which cannot achieve accurate elimination. Therefore, in this embodiment, a hierarchical iterative approach is adopted to gradually construct dynamic apparent density models at different accuracy levels to reflect the true terrain gravity effect.

[0029] First iteration: Eliminate errors in the uniform density assumption. Specifically: It should be noted that using a globally uniform density assumption can introduce significant systematic errors during terrain compensation in the mid-to-far regions. Therefore, we first obtain the regional average apparent density through inversion and construct a first-order dynamic apparent density model to replace the traditional globally uniform density, thus eliminating the maximum systematic error.

[0030] First, perform anomalous gravity decomposition: The db4 wavelet function was used to analyze the initial Bouguer gravity anomaly corresponding to the target region. Three-level wavelet decomposition was performed to obtain the anomalous gravity components at different wavelengths, specifically: three detail components and one approximate component.

[0031] It should be noted that different wavelengths correspond to different gravity anomalies. In this embodiment, the first wavelength is set to less than or equal to 60m (first layer detail, shallow interference), the second wavelength is greater than 60m and less than or equal to 120m (second layer detail), the third wavelength is greater than 120m and less than or equal to 240m (third layer detail), and the fourth wavelength is greater than 240m (approximate component, deep structural anomaly).

[0032] Following the above method, the db4 wavelet function is then used to perform three-level wavelet decomposition on the DEM data corresponding to the target area to obtain the terrain components at different wavelengths.

[0033] Then, the anomalous gravity components related to the terrain are filtered out and reconstructed: It should be noted that there are many causes of gravity anomalies. In order to extract gravity anomalies related to topographic relief and perform topographic compensation, the spatial cross-correlation coefficient between different components and topographic data is calculated. The spatial cross-correlation coefficient reflects the correlation between gravity anomalies and topography.

[0034] Calculate the cross-correlation coefficient between the anomalous gravity component and the topographic component at the same wavelength, and denot it as the gravity-topographic cross-correlation coefficient at each wavelength.

[0035] It should be noted that calculating the cross-correlation coefficient between the anomalous gravity component and the topographic component at the same wavelength is a classic and well-known routine calculation in the fields of geoscience and geophysics. The larger the cross-correlation coefficient at this wavelength, the more likely the gravity anomaly at this wavelength is caused by topographic relief.

[0036] In this embodiment, the preset cross-correlation coefficient threshold is 0.7, which is obtained through data statistics. When the cross-correlation coefficient is greater than or equal to 0.7, it indicates that the correlation between abnormal gravity and terrain is significant.

[0037] If the gravity-terrain cross-correlation coefficient at each wavelength is greater than or equal to the preset cross-correlation coefficient threshold, then each wavelength is recorded as the target wavelength.

[0038] If the gravity-topography cross-correlation coefficient at each wavelength is less than the preset cross-correlation coefficient threshold, then each wavelength is recorded as an invalid wavelength.

[0039] It should be noted that gravity anomalies at the target wavelength are mainly caused by topographic relief and should be retained. Non-topographic anomalies at invalid wavelengths should be discarded. In this embodiment, at least one target wavelength is selected, which is ensured by adjusting the preset cross-correlation coefficient threshold. This ensures the success of subsequent wavelet reconstruction.

[0040] By performing inverse wavelet transform on the anomalous gravity components at all target wavelengths, the pure topographic gravity anomaly field corresponding to the target area can be reconstructed.

[0041] Finally, a first-order dynamic apparent density model is constructed and the gravity field is updated using the reconstructed pure terrain gravity anomaly field: Input the pure topographic gravity anomaly field corresponding to the target area into the improved The model is used to invert the overall average apparent density value corresponding to the target area, construct a first-order dynamic apparent density model (regional average apparent density model), and combine it with the DEM data corresponding to the target area to calculate the first iteration of the first round of terrain compensation value using a terrain correction formula that considers the curvature of the earth. .

[0042] The formula for topographic correction considering the curvature of the Earth (spherical topographic correction) is a well-known technique, and the specific method will not be introduced here.

[0043] The initial Bouguer gravity anomaly corresponding to the target area All elements are superimposed with the terrain compensation value from the first iteration of the first round. The first iteration of the first round of Bouguer gravity anomaly was obtained. .

[0044] The preset convergence deviation for the first round is: The default convergence count is 50, and this will be used as an example for explanation.

[0045] Calculate the terrain compensation value in the first iteration of the first round. Compared with the preset initial density compensation value relative deviation .

[0046] in, It is an absolute value function. It is 2.67 g / cm 3 The denominator is not 0.

[0047] like The first iteration of the first round will then show the Bouguer gravity anomaly. This is recorded as the first update of the Bouguer gravity anomaly. The terrain compensation value of the first iteration of the first round This is recorded as the first round of terrain compensation value. .

[0048] like Then continue to solve for the terrain compensation value in the second iteration of the first round. And the Bouguer gravity anomaly in the first round of the second iteration .

[0049] It should be noted that, following the above method, the Bouguer gravity anomaly in the first iteration of the first round... The system is decomposed, the anomalous gravity components at the target wavelength are screened out, and then reconstructed. Based on the reconstructed data... Calculate the terrain compensation value for the first round of the second iteration. , used to Compensation was performed to obtain the Bouguer gravity anomaly in the first round of the second iteration. .

[0050] Calculate the terrain compensation value in the first round of the second iteration. Compared with the terrain compensation value in the first iteration of the first round relative deviation .

[0051] like Then the Bouguer gravity anomaly in the first round of the second iteration will be... This is recorded as the first update of the Bouguer gravity anomaly. The terrain compensation value of the first round of the second iteration This is recorded as the first round of terrain compensation value. .

[0052] like Then, this process is repeated iteratively to calculate the terrain compensation value and the corresponding Bouguer gravity anomaly for each iteration, until the first iteration... Sub-iteration terrain compensation value Compared with the first round Sub-iteration terrain compensation value relative deviation ,in Once the convergence condition is met, the iteration stops, and the first round of iterations is completed. The next iteration of the Bouguer gravity anomaly This is recorded as the first update of the Bouguer gravity anomaly. The first round Sub-iteration terrain compensation value This is recorded as the first round of terrain compensation value. .

[0053] If the preset convergence deviation requirement for the first round is not met after 50 consecutive iterations in the first round, the iteration will be forcibly terminated, and the Bouguer gravity anomaly in the 50th iteration of the first round will be recorded. This is recorded as the first update of the Bouguer gravity anomaly. The terrain compensation value of the 50th iteration in the first round This is recorded as the first round of terrain compensation value. .

[0054] It should be noted that the first round of updates to the Bouguer gravity anomaly... Once the uniform density error is eliminated, proceed to the second iteration.

[0055] Step S003: Based on the pure terrain gravity anomaly field after the decomposition and reconstruction of the Bouguer gravity anomaly in the first round, perform multiple rounds of iterative compensation and correction until the termination condition is met, and output the final updated Bouguer gravity anomaly.

[0056] Second iteration: Eliminate density spatial heterogeneity error. Specifically: It should be noted that, due to the extremely large coverage area of ​​the mid-to-far zone, the aforementioned regional average apparent density cannot reflect the differences between different lithological units in different regions. Therefore, it is necessary to further replace the overall average apparent density value corresponding to the target area in the first round with a grid-by-grid apparent density field with a resolution of 30m, construct a second-order dynamic apparent density model, capture local lithological changes within the same geological unit, and eliminate the spatial heterogeneity error of density in the topographic correction of the mid-to-far zone.

[0057] The first round of updates will include the Bouguer gravity anomaly. The corresponding pure terrain gravity anomaly field is denoted as the second round initial pure terrain gravity anomaly field.

[0058] It should be noted that, following the above method, the first round of updates to the Bouguer gravity anomaly... The system is decomposed, the anomalous gravity components at the target wavelength are selected, and then reconstructed to obtain... The corresponding pure topographic gravity anomaly field.

[0059] The default grid size is 30m×30m, and this will be used as an example for explanation.

[0060] The second round initial pure terrain gravity anomaly field is input into the improved... The model was inverted to obtain a two-dimensional grid-by-grid apparent density field with a resolution of 30m×30m. A threshold constraint was applied to the inverted density output, limiting the apparent density value to between 2.3 and 2.8 g / cm³. 3 Within a reasonable range, a second-order dynamic apparent density model (two-dimensional grid-by-grid apparent density model) is constructed. Combined with the DEM data corresponding to the target area, the elevation information corresponding to each grid is spatially coupled with the two-dimensional grid-by-grid apparent density field to expand and construct a three-dimensional apparent density field corresponding to the target area, and earth curvature correction processing is completed. Each 30m×30m DEM grid cell corresponds to one true apparent density value.

[0061] Based on the three-dimensional apparent density field corresponding to the target area, the terrain compensation value for the first iteration of the second round is calculated by integrating grid-by-grid cell using a terrain correction formula that takes into account the curvature of the Earth. .

[0062] The first round of updates will include the Bouguer gravity anomaly. All elements are superimposed with the terrain compensation value from the first iteration of the second round. The second round of the first iteration yielded the Bouguer gravity anomaly. .

[0063] The preset convergence bias for the second round is: Let's take this as an example.

[0064] Calculate the terrain compensation value in the first iteration of the second round. Compared with the first round of terrain compensation values relative deviation .

[0065] like Then the Bouguer gravity anomaly in the first iteration of the second round will be... This is recorded as the second update of the Bouguer gravity anomaly. The terrain compensation value of the first iteration of the second round This is recorded as the second round of terrain compensation value. .

[0066] like Then continue to solve for the terrain compensation value in the second iteration of the second round. And the second round, second iteration of Bouguer gravity anomaly .

[0067] It should be noted that, following the above method, the Bouguer gravity anomaly in the first iteration of the second round was analyzed. The system is decomposed, the anomalous gravity components at the target wavelength are screened out, and then reconstructed. Based on the reconstructed data... Calculate the terrain compensation value for the second round of the second iteration. , used to Compensation was performed to obtain the second round, second iteration of Bouguer gravity anomaly. .

[0068] Calculate the terrain compensation value in the second round of the second iteration. Compared with the terrain compensation value of the first iteration of the second round relative deviation .

[0069] like Then the second round, second iteration of the Bouguer gravity anomaly will be... This is recorded as the second update of the Bouguer gravity anomaly. The second round of terrain compensation values ​​in the second iteration This is recorded as the second round of terrain compensation value. .

[0070] like Then, this process is repeated iteratively to calculate the terrain compensation value and the corresponding Bouguer gravity anomaly for each iteration of the second round, until the second round... Sub-iteration terrain compensation value Compared with the second round Sub-iteration terrain compensation value relative deviation ,in Once the convergence condition is met, the iteration stops, and the second round of iteration... The next iteration of the Bouguer gravity anomaly This is recorded as the second update of the Bouguer gravity anomaly. The second round Sub-iteration terrain compensation value This is recorded as the second round of terrain compensation value. .

[0071] If the second round of iterations reaches the preset convergence count of 50 and still fails to meet the preset convergence deviation requirement, the iteration will be forcibly terminated, and the Bouguer gravity anomaly in the 50th iteration of the second round will be recorded. This is recorded as the second update of the Bouguer gravity anomaly. The terrain compensation value of the 50th iteration in the second round This is recorded as the second round of terrain compensation value. .

[0072] It should be noted that the second round of updates to the Bouguer gravity anomaly... The spatial heterogeneity error of medium density is eliminated, and the third iteration begins.

[0073] Third iteration: Eliminate residual errors from abnormal separation. Specifically: It should be noted that the scales of mid-to-far-range topography and deep structures overlap (the wavelengths of some deep structures are similar to those of large-scale topography), causing deep structure signals to deeply couple with topography signals. Traditional anomaly separation methods, due to interference from the aforementioned two density differences, cannot distinguish these coupled signals. Only after all density-related errors (uniform density error, density spatial heterogeneity error) are eliminated will the residual anomaly separation error become the main error source and be detected and corrected. Therefore, it is further necessary to refine the second-order dynamic apparent density model by stripping away the residual deep structure signals deeply coupled with the topography signal, constructing the final third-order dynamic apparent density model (high-precision three-dimensional apparent density model), eliminating the residual anomaly separation error in mid-to-far-range topography correction, and achieving final accuracy convergence.

[0074] The second round of updates to the Bouguer gravity anomaly will be performed. The corresponding pure topographic gravity anomaly field is denoted as the third round initial pure topographic gravity anomaly field.

[0075] It should be noted that, following the above method, the second round of updates to the Bouguer gravity anomaly... The system is decomposed, the anomalous gravity components at the target wavelength are selected, and then reconstructed to obtain... The corresponding pure topographic gravity anomaly field contains almost no deep structural signals in the third round of initial pure topographic gravity anomaly field.

[0076] The third round initial pure terrain gravity anomaly field is input into the improved... The model is based on the grid-by-grid two-dimensional apparent density field obtained in the second round of inversion. Combined with the elevation information provided by the DEM data corresponding to the target area, the two-dimensional apparent density field is refined grid-by-grid and optimized in three dimensions through the model’s feature learning and accurate fitting capabilities, resulting in the final high-precision three-dimensional apparent density field (high-precision third-order dynamic apparent density model).

[0077] Based on the final high-precision three-dimensional apparent density field, the terrain compensation value for the first iteration of the third round is calculated by integrating grid-by-grid cell using a terrain correction formula that takes into account the curvature of the Earth. .

[0078] The second round of updates to the Bouguer gravity anomaly will be performed. All elements are superimposed with the terrain compensation value from the first iteration of the third round. The third round of first iteration yielded the Bouguer gravity anomaly. .

[0079] The preset convergence bias for the third round is: Let's take this as an example.

[0080] Calculate the terrain compensation value in the first iteration of the third round. Compared with the second round of terrain compensation values relative deviation .

[0081] like Then the third round of the first iteration of the Bouguer gravity anomaly will be... This is recorded as the final update of the Bouguer gravity anomaly. .

[0082] like Then continue to solve for the terrain compensation value in the second iteration of the third round. And the Bouguer gravity anomaly in the third round, second iteration .

[0083] It should be noted that, following the above method, the Bouguer gravity anomaly in the first iteration of the third round was analyzed. The system is decomposed, the anomalous gravity components at the target wavelength are screened out, and then reconstructed. Based on the reconstructed data... Calculate the terrain compensation value for the second iteration of the third round. , used to Compensation was performed to obtain the Bouguer gravity anomaly in the third round, second iteration. .

[0084] Calculate the terrain compensation value in the second iteration of the third round. Compared with the terrain compensation value of the first iteration of the third round relative deviation .

[0085] like Then the third round, second iteration of the Bouguer gravity anomaly will be... This is recorded as the final update of the Bouguer gravity anomaly. .

[0086] like Then, this process is repeated iteratively to calculate the terrain compensation value and the corresponding Bouguer gravity anomaly for each iteration in the third round, until the third round... Sub-iteration terrain compensation value Compared with the third round Sub-iteration terrain compensation value relative deviation ,in Once the convergence condition is met, the iteration stops, and the third round of iteration... The next iteration of the Bouguer gravity anomaly This is recorded as the final update of the Bouguer gravity anomaly. .

[0087] If the preset convergence deviation requirement for the third round is not met after 50 consecutive iterations in the third round, the iteration will be forcibly terminated, and the Bouguer gravity anomaly in the 50th iteration of the third round will be recorded. This is recorded as the final update of the Bouguer gravity anomaly. .

[0088] It should be noted that during a single round of continuous iterative calculations, if the terrain compensation value obtained in a certain iteration is 0, it indicates that there is no new terrain correction in that iteration, and the Bouguer gravity anomaly in this iteration is completely consistent with the Bouguer gravity anomaly in the previous iteration. At this point, the residual terrain influence has been completely corrected, and there is no need to continue iterating and correcting. The entire iteration process is then terminated, and the Bouguer gravity anomaly of this iteration is designated as the final updated Bouguer gravity anomaly. This judgment rule can effectively avoid the problem of zero denominator and calculation failure in subsequent adjacent iterations of relative deviation calculation due to the previous compensation value being 0, thus ensuring the stability of iterative calculation. In this embodiment, a schematic diagram of the hierarchical progressive multi-round terrain compensation iteration architecture is shown below. Figure 2 As shown.

[0089] It should be further noted that this embodiment employs an improved method. The model constructs a three-dimensional apparent density field inversion model under multiple constraints. Traditional three-dimensional apparent density inversion aims to interpret deep underground structures, lacking spatial continuity constraints and geological hard constraints for terrain bodies. Inversion results often exhibit unreasonable density jumps and mixed physical meanings. Therefore, this embodiment specifically optimizes the model structure and loss function for terrain density inversion, strictly limiting the inversion space to above sea level to ensure that the output of the dynamic apparent density model uniquely corresponds to the true equivalent apparent density of the terrain body, thereby adapting to the aforementioned three-layer iterative method and improving terrain correction accuracy. To improve the adaptability of the dynamic apparent density model to mid-to-far-field terrain, adaptation modifications are needed to the model structure and loss function. Specifically: (1) Adaptation Model: In the basic A spatial attention module is added between the encoder and decoder of the model. This module adaptively amplifies features in areas with dramatic terrain undulations and suppresses interference from irrelevant features by calculating attention weights for each spatial location. The activation function used is... negative half-axis slope coefficient Add after each convolutional layer Layers are used to accelerate model convergence and prevent overfitting. The model input is a 32×32×1 pure topographic gravity anomaly grid block, and the output is a 32×32×1 average apparent density grid block, corresponding to the average apparent density of the strata above sea level.

[0090] (2) Design a multi-loss function that integrates mean square error, spatial gradient regularization, and borehole hard constraint: To ensure the prediction accuracy, spatial continuity, and borehole hard constraint of the dynamic apparent density model, a multi-loss function combining mean squared error, spatial gradient regularization, and borehole hard constraint is used. In this embodiment, the label apparent density is used as the target ground truth, and the mean squared error loss between the predicted apparent density and the label apparent density is calculated. To constrain the spatial stationarity of the apparent density field and suppress local density abrupt changes, spatial gradient regularization is introduced to construct a continuity loss function. Since the measured data points have the highest reliability, the boundary rigid constraints are constructed using the prior information from the measured points, and then the rigid constraint loss is obtained. The basic calculation methods for mean squared error loss and continuity loss, as well as the basic design principles of rigid constraints, are common knowledge in this field, and the specific methods will not be introduced here.

[0091] Final loss function for:

[0092] In the formula, , To preset the first weighting coefficient, , Using a pre-defined second weighting coefficient as an example, we will describe the mean squared error loss. The larger the value, the greater the difference between the overall predicted visual density and the labeled visual density. (Continuity loss function) The larger the value, the more severe the density jumps between different regions. (Rigid constraint loss) The larger the value, the greater the deviation between the predicted visual density and the measured visual density.

[0093] (3) Dataset construction and model training: The pure topographic gravity anomaly field was divided into 32×32 grid blocks, with adjacent grid blocks overlapping by 16 points to avoid boundary effects. Label data consisted of borehole apparent density data and apparent density distribution inferred from geological maps. Labels for borehole locations were the measured apparent density values, while labels for non-borehole locations ranged from 2.3 to 2.8 g / cm³, based on the corresponding lithology on the geological map. 3 The dataset is assigned values ​​within a specified range and divided into training and validation sets in a 7:3 ratio. The Adam optimizer is used with an initial learning rate of 0.001. A cosine annealing learning rate update strategy is employed (period Tmax = 100), with a batch size of 16 and 1000 iterations. Training stops when the validation set loss function value is ≤0.001 and shows no significant decrease after 50 consecutive iterations, resulting in the trained dynamic apparent density inversion model.

[0094] (4) Generation of three-dimensional apparent density field: The pure terrain gravity anomaly field is input into the trained model, and a two-dimensional average apparent density field with a resolution of 30m×30m is output. Combined with the elevation information of the DEM data, the two-dimensional average apparent density field is expanded into a three-dimensional apparent density field, with each 30m×30m×1m voxel corresponding to an apparent density value.

[0095] This invention is now complete.

[0096] In summary, in this embodiment of the invention, the initial Bouguer gravity anomaly and DEM data are decomposed at multiple scales to obtain gravity anomaly components and terrain components corresponding to different wavelengths. Based on the cross-correlation coefficient between the anomalous gravity components and terrain components at the same wavelength, at least one type of target wavelength is selected. The anomalous gravity components at all target wavelengths are reconstructed to obtain a pure terrain gravity anomaly field. Based on the pure terrain gravity anomaly field, the first iteration terrain compensation value of the first round is calculated to compensate and correct the initial Bouguer gravity anomaly. The first iteration Bouguer gravity anomaly is determined, and the first round of iterative compensation and correction is carried out until the convergence condition is met. The first round updated Bouguer gravity anomaly is output. Based on the pure terrain gravity anomaly field after decomposition and reconstruction of the first round updated Bouguer gravity anomaly, multiple rounds of iterative compensation and correction are carried out until the termination condition is met. The final updated Bouguer gravity anomaly is output. This invention achieves high-precision variable density terrain compensation by constructing a complete technical link of adaptive anomaly separation, multi-constraint density inversion, partitioned coupling compensation, and iterative optimization closed loop.

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

Claims

1. An iterative compensation method for distant terrain in gravity anomalies based on a dynamic apparent density model, characterized in that, The method includes the following steps: Obtain the initial Bouguer gravity anomaly and DEM data corresponding to the target area; Multi-scale decomposition is performed on the initial Bouguer gravity anomaly and DEM data to obtain gravity anomaly components and terrain components corresponding to different wavelengths. At least one type of target wavelength is selected based on the cross-correlation coefficient between the anomalous gravity component and the terrain component at the same wavelength. The anomalous gravity components at all target wavelengths are reconstructed to obtain a pure terrain gravity anomaly field. Based on the pure terrain gravity anomaly field, the terrain compensation value of the first iteration in the first round is calculated to compensate and correct the initial Bouguer gravity anomaly, the Bouguer gravity anomaly of the first iteration in the first round is determined, and the first round of iterative compensation and correction is carried out until the convergence condition is met, and the updated Bouguer gravity anomaly of the first round is output. Based on the pure terrain gravity anomaly field after the first round of updated Bouguer gravity anomaly decomposition and reconstruction, multiple rounds of iterative compensation and correction are performed until the termination condition is met, and the final updated Bouguer gravity anomaly is output.

2. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 1, characterized in that, The specific steps involved in selecting at least one type of target wavelength are as follows: Calculate the cross-correlation coefficient between the anomalous gravity component and the topographic component at the same wavelength, and denot it as the gravity-topographic cross-correlation coefficient at each wavelength. If the gravity-terrain cross-correlation coefficient at each wavelength is greater than or equal to the preset cross-correlation coefficient threshold, then each wavelength is recorded as the target wavelength.

3. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 1, characterized in that, The specific steps involved in determining the Bouguer gravity anomaly in the first iteration of the first round are as follows: Inputting the pure topographic gravity anomaly field into the improved The model is inverted to obtain the overall average apparent density value corresponding to the target area. Combined with DEM data, the terrain compensation value of the first iteration of the first round is calculated using the terrain correction formula that takes into account the curvature of the earth. The first iteration of the first round of Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first round of the first iteration onto all elements in the initial Bouguer gravity anomaly.

4. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 1, characterized in that, The specific steps involved in the first round of updating the Bouguer gravity anomaly are as follows: If the relative deviation between the terrain compensation value of the first iteration of the first round and the preset initial density compensation value is less than or equal to the preset first round convergence deviation, then the Bouguer gravity anomaly of the first iteration of the first round is recorded as the first round updated Bouguer gravity anomaly. If the relative deviation between the terrain compensation value in the first iteration of the first round and the preset initial density compensation value is greater than the preset first round convergence deviation, then this process is repeated iteratively to determine the updated Bouguer gravity anomaly in the first round.

5. The method for iterative compensation of gravity anomaly in distant regions based on a dynamic apparent density model according to claim 4, characterized in that, If the relative deviation between the terrain compensation value in the first iteration of the first round and the preset initial density compensation value is greater than the preset first round convergence deviation, then this process is repeated iteratively to determine the updated Bouguer gravity anomaly in the first round. The specific steps include the following: The terrain compensation value and the corresponding Bouguer gravity anomaly are solved successively in the first iteration until the first iteration... The terrain compensation value in the next iteration is the same as that in the first round. The iteration stops when the relative deviation of the terrain compensation value in the next iteration is less than or equal to the preset first-round convergence deviation, and the first-round iteration is stopped after the convergence condition is met. The Bouguer gravity anomaly in the next iteration is denoted as the first round updated Bouguer gravity anomaly. If the preset convergence deviation requirement for the first round is still not met after the first round of consecutive iterations reaches the preset convergence number, the iteration is forcibly terminated, and the Bouguer gravity anomaly of the last iteration of the first round is recorded as the first round updated Bouguer gravity anomaly.

6. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 1, characterized in that, The output ultimately updates the Bouguer gravity anomaly, and the specific steps involved are as follows: The pure terrain gravity anomaly field, obtained from the first round of Bouguer gravity anomaly decomposition and reconstruction, is input into the improved... The model is inverted to obtain a two-dimensional grid-by-grid apparent density field with a preset grid size resolution. Combined with DEM data, a three-dimensional apparent density field is constructed. Based on the three-dimensional apparent density field, a terrain correction formula considering the curvature of the earth is adopted. The terrain compensation value of the first iteration of the second round is calculated by integrating grid-by-grid unit. The updated Bouguer gravity anomaly of the first round is compensated and corrected. The Bouguer gravity anomaly of the first iteration of the second round is determined. After the second round of iteration is completed and converged, the updated Bouguer gravity anomaly of the second round is output. Based on the pure terrain gravity anomaly field after the decomposition and reconstruction of the second round of updated Bouguer gravity anomaly, the terrain compensation value of the first iteration of the third round is calculated to compensate and correct the second round of updated Bouguer gravity anomaly, and the first iteration of the third round of Bouguer gravity anomaly is determined. After the third round of iteration is completed and converged, the final updated Bouguer gravity anomaly is output.

7. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 6, characterized in that, The specific steps involved in determining the Bouguer gravity anomaly in the first iteration of the second round are as follows: The first iteration of the second round of Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first iteration of the second round on all elements of the first round of updated Bouguer gravity anomaly.

8. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 6, characterized in that, The specific steps involved in the second round of updating the Bouguer gravity anomaly are as follows: The terrain compensation value corresponding to the first round of updating the Bougu gravity anomaly will be obtained and recorded as the first round terrain compensation value. If the relative deviation between the terrain compensation value of the first iteration of the second round and the terrain compensation value of the first round is less than or equal to the preset convergence deviation of the second round, then the Bouguer gravity anomaly of the first iteration of the second round is recorded as the updated Bouguer gravity anomaly of the second round. If the relative deviation between the terrain compensation value in the first iteration of the second round and the terrain compensation value in the first round is greater than the preset convergence deviation for the second round, then this process is repeated iteratively to solve for the terrain compensation value and the corresponding Bouguer gravity anomaly in each iteration of the second round, until the second round... The terrain compensation value of the second iteration and the second round The iteration stops when the relative deviation of the terrain compensation value in the second iteration is less than or equal to the preset second-round convergence deviation, and the second-round iteration is stopped after the convergence condition is met. The Bouguer gravity anomaly in the next iteration is denoted as the second-round updated Bouguer gravity anomaly. If the second round of consecutive iterations reaches the preset convergence number but still fails to meet the preset second round convergence deviation requirement, the iteration is forcibly terminated, and the Bouguer gravity anomaly of the last iteration of the second round is recorded as the second round updated Bouguer gravity anomaly.

9. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 6, characterized in that, The specific steps involved in determining the Bouguer gravity anomaly in the first iteration of the third round are as follows: The pure terrain gravity anomaly field, after the second round of updated Bouguer gravity anomaly decomposition and reconstruction, is input into the improved... The model, based on the two-dimensional grid-by-grid apparent density field with a preset grid size resolution obtained by inversion, is combined with DEM data to correct and expand to obtain the final high-precision three-dimensional apparent density field; based on the final high-precision three-dimensional apparent density field, the terrain correction formula considering the curvature of the earth is adopted, and the terrain compensation value of the first iteration of the third round is calculated by grid-by-grid integration. The third round of first iteration Bouguer gravity anomaly is obtained by superimposing the terrain compensation value of the first iteration on all elements in the second round of updated Bouguer gravity anomaly.

10. The method for iterative compensation of distant terrain in gravity anomalies based on a dynamic apparent density model according to claim 6, characterized in that, The output ultimately updates the Bouguer gravity anomaly, and the specific steps involved are as follows: The terrain compensation value corresponding to the second round of update of the Bougu gravity anomaly will be obtained and recorded as the second round terrain compensation value. If the relative deviation between the terrain compensation value of the first iteration of the third round and the terrain compensation value of the second round is less than or equal to the preset convergence deviation of the third round, then the Bouguer gravity anomaly of the first iteration of the third round is recorded as the final updated Bouguer gravity anomaly. If the relative deviation between the terrain compensation value in the first iteration of the third round and the terrain compensation value in the second round is greater than the preset convergence deviation for the third round, then this process is repeated iteratively to solve for the terrain compensation value and the corresponding Bouguer gravity anomaly in each iteration of the third round, until the third round... The terrain compensation value of the second iteration and the third round The iteration stops when the relative deviation of the terrain compensation value in the third iteration is less than or equal to the preset convergence deviation in the third round, and the convergence condition is met. The next iteration of the Bouguer gravity anomaly is denoted as the final updated Bouguer gravity anomaly. If the third round of consecutive iterations reaches the preset convergence number but still fails to meet the preset convergence deviation requirement, the iteration is forcibly terminated, and the Bouguer gravity anomaly in the last iteration of the third round is recorded as the final updated Bouguer gravity anomaly.

Images