A method for obtaining gravity anomaly blank area data based on potential field characteristics

Abstract

The present application belongs to the field of gravity exploration, and particularly relates to a method for obtaining data of gravity anomaly blank area based on characteristics of potential field, which comprises the following steps: obtaining gravity and terrain data of an exploration area and its neighborhood, generating standard grid data and identifying blank area after pretreatment and gridding; selecting optimal model to obtain reliable Bouguer gravity anomaly value of the blank area by parallel computing and cross-validation quantitative evaluation of multiple interpolation algorithms; eliminating noise by mathematical analysis and filtering and smoothing processing of the interpolation result to obtain Bouguer gravity anomaly value meeting physical continuity and smoothness of the gravity field; and superimposing the processed Bouguer gravity anomaly and terrain gravity effect to obtain complete spatial gravity anomaly field. The present application successfully realizes filling and balancing of gravity anomaly value in the area with complex terrain and sparse data by system processing of existing Bouguer gravity anomaly and terrain data and using interpolation, analysis, smoothing and reduction techniques.

Description

Technical Field

[0001] This invention belongs to the field of gravity exploration, specifically relating to a method for obtaining gravity anomaly blank areas based on potential field characteristics. Background Technology

[0002] In geophysical exploration and crustal structure research, gravity anomalies are crucial physical information revealing subsurface mass distribution, tectonic morphology, and geological evolution. Spatial gravity anomalies reflect the relationship between topographic relief and gravity variations, while Bouguer anomalies primarily reveal differences in subsurface media density and crustal tectonic characteristics, serving as fundamental data for studying isostatics, surface structure, and crustal thickness. However, in areas with complex terrain or harsh environments, limitations in transportation, climate, and equipment often result in significant gaps or uneven distribution of measured gravity data. This not only affects the integrity of the local gravity field but also introduces systematic errors, reducing the accuracy of regional gravity field models and the reliability of interpretation.

[0003] Currently, existing technical solutions to the problem of missing or sparse gravity anomaly data in areas with complex terrain can be summarized into three categories:

[0004] The first is spatial interpolation, which uses the gravity anomaly values, geographic coordinates and auxiliary variables of the measured points to estimate the gravity values ​​of blank points. Common algorithms include inverse distance weighting, kriging interpolation and interpolation based on equivalent sources. These methods are suitable for areas with relatively uniform data distribution, but in complex areas such as areas with drastic terrain undulations, they often lead to interpolation distortion, abnormal ambiguity or abrupt boundary changes, making it difficult to reflect the real terrain gravity coupling relationship and the calculation results lack physical consistency.

[0005] The second method is numerical extrapolation. Based on the physical property that the gravity potential field satisfies the Laplace equation, it uses mathematical extrapolation to convert gravity data between different elevation planes. Common methods include upward extrapolation, downward extrapolation, and equivalent layer extrapolation. This method has good physical consistency, but it is extremely sensitive to observation errors in actual calculations. It is prone to noise amplification, artifact superposition, and local anomalous distortion, making it difficult to achieve stable and accurate gravity anomaly compensation in data gaps.

[0006] Thirdly, numerical simulation and inversion methods, based on the forward modeling equations of the gravity field, construct underground density models and use numerical optimization or error minimization methods to invert underground structures. Common methods include two-dimensional or three-dimensional grid inversion and layered model inversion. However, the inversion process involves large computational loads and complex algorithms, and is significantly dependent on prior models, regularization parameters and constraints. When the distribution of observation data is sparse or there are gaps in the survey area, the inversion results are prone to multiple solutions and instability, and the characterization of the influence of terrain is not sufficient. When applied to areas with variable terrain, the stability and applicability are limited.

[0007] In summary, existing methods generally have three limitations: First, they rely too heavily on the density of observation data, resulting in a significant decrease in interpolation accuracy in areas with missing data; second, they fail to fully consider the physical coupling relationship between topography and gravity anomalies, leading to insufficient physical consistency of the results; and third, they are computationally complex or have poor stability, making it difficult to balance accuracy and practicality. Summary of the Invention

[0008] This invention overcomes the above-mentioned defects and provides a method for obtaining gravity anomaly blank area data based on potential field characteristics. It solves the problems of poor stability, insufficient physical consistency and high dependence on measured data in the application of the existing gravity anomaly interpolation and estimation methods in areas with variable terrain.

[0009] To achieve the above objectives, the present invention provides a method for obtaining gravity anomaly blank area data based on potential field characteristics, the method comprising the following steps:

[0010] S1. Obtain gravity and topographic data of the exploration area and its surrounding area, preprocess and grid the data to generate standard grid data and mark blank areas;

[0011] S2. Perform interpolation calculations on the data and conduct quantitative evaluation to obtain reliable Bouguer gravity anomaly values ​​in the blank area;

[0012] S3. Perform mathematical analysis and filtering smoothing on the interpolation results to eliminate noise and obtain Bouguer gravity anomaly values ​​that conform to the physical continuity and smoothness of the gravity field.

[0013] S4. The processed Bouguer gravity anomaly is superimposed with the topographic gravity effect to obtain the complete spatial gravity anomaly field.

[0014] This method is based on existing Bouguer gravity anomaly data and digital elevation model data. It comprehensively utilizes the physical laws that Bouguer gravity anomalies are basically independent of terrain elevation and that spatial gravity anomalies are highly correlated with elevation. Based on Bouguer gravity anomaly data from neighboring areas, it performs interpolation calculations for Bouguer gravity anomalies in blank areas. Simultaneously, it introduces three interpolation methods: inverse distance weighting, ordinary kriging, and multiple regression models. Cross-validation is used to verify the accuracy of the three interpolation methods, quantitatively evaluate the calculation effect, and select the optimal interpolation model to ensure that the interpolation results of Bouguer gravity anomalies in blank areas have the minimum numerical deviation and the highest stability. The obtained interpolation results are mathematically analyzed and filtered to make the Bouguer gravity anomaly field physically conform to the continuity and analyzability characteristics of the potential field. The physically decoupled and numerically smoothed Bouguer anomaly field is superimposed with the gravity effect caused by terrain elevation to recover the physically complete spatial gravity anomaly field. It successfully achieves the filling and balancing of gravity anomaly values ​​in areas with complex terrain and sparse data.

[0015] Compared with the prior art, the advantages of the present invention are as follows:

[0016] (1) Based on the physical mechanism, this invention comprehensively utilizes the differences in the correlation between spatial gravity anomaly and Bouguer gravity anomaly in terms of terrain correlation, combines the strong correlation between spatial gravity anomaly and terrain elevation with the non-correlation of Bouguer gravity anomaly, breaks through the empirical assumptions of traditional mathematical interpolation, and establishes a physical constraint model. This invention can deduce the Bouguer gravity anomaly and spatial gravity anomaly in the blank area based solely on terrain elevation and neighborhood gravity information in the absence of measured gravity data, thereby realizing the restoration and balance filling of the gravity field.

[0017] (2) In the Bouguer gravity anomaly interpolation stage, this invention introduces a multi-algorithm parallel and cross-validation mechanism. At the same time, it uses the inverse distance weighting method, ordinary kriging method and multiple regression model for interpolation calculation. Through cross-validation, it calculates three indicators: root mean square error (RMSE), mean absolute error (MAE) and coefficient of determination (R²) to quantitatively evaluate the model performance. The model result with the smallest error and the highest goodness of fit is selected, which greatly improves the interpolation accuracy and the reliability of the results.

[0018] (3) In the analytical stage, the present invention innovatively proposes a dual processing scheme of "mathematical analysis + filtering and smoothing". Through analytical extension based on the Laplace equation, the Bouguer gravity anomaly field is made mathematically continuous and physically analytical. Then, the frequency domain Gaussian low-pass filter and the spatial domain Laplace smoothing operator are used to remove high-frequency noise and local discontinuities, and a smooth Bouguer anomaly field that conforms to the potential field characteristics is obtained. This "analysis-filtering-smoothing" joint method effectively suppresses high-frequency oscillations and boundary abrupt changes in traditional interpolation. Attached Figure Description

[0019] Figure 1 This is a flowchart of the method for obtaining gravity anomaly blank area data based on potential field characteristics according to the present invention;

[0020] Figure 2 This is a flowchart of the steps for calculating gravity anomalies in blank areas according to the present invention;

[0021] Figure 3 These are topographic maps of blank areas and Bouguer gravity anomaly distribution maps of adjacent areas, as described in this embodiment of the invention.

[0022] in, Figure 3 (a) shows the topographic map of the blank area and adjacent areas. Figure 3 (b) shows the distribution of Bouguer gravity anomalies in the adjacent areas of the blank area;

[0023] Figure 4 This is a comparison chart of Bouguer gravity anomaly interpolation results according to an embodiment of the present invention;

[0024] in, Figure 4 Figure (a) shows the interpolation result using the ordinary kriging method. Figure 4Image (b) shows the interpolation result using the inverse distance weighted method. Figure 4 (c) shows the interpolation results of the multiple regression model;

[0025] Figure 5 This is a Bouguer gravity anomaly diagram after mathematical analysis and filtering smoothing according to an embodiment of the present invention;

[0026] Figure 6 This is a topographic gravity effect diagram according to an embodiment of the present invention;

[0027] Figure 7 This is the final spatial gravity anomaly diagram of an embodiment of the present invention. Detailed Implementation

[0028] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. This embodiment takes the western edge of a basin as the research object. This area has drastic topography and sparse gravity survey lines, and is a typical representative of a region with variable topography and discontinuous gravity data. The method proposed in this invention for obtaining gravity anomaly blank area data based on potential field characteristics is used to estimate the Bouguer gravity anomaly and spatial gravity anomaly in the blank area of ​​this region and reconstruct the gravity field.

[0029] like Figure 1 , 2 As shown, this invention proposes a method for obtaining gravity anomaly blank area data based on potential field characteristics, specifically including:

[0030] Step S1: Obtain geological information such as slope and topographic curvature of the area to be explored, high-precision digital elevation model (DEM) data, and Bouguer gravity anomaly data of the area adjacent to the blank area.

[0031] Specifically, Bouguer gravity anomaly data is typically stored in .dat, .grd, or ASCII table format, containing latitude and longitude (L... i B i The data includes Bouguer gravity anomaly values ​​and elevation information of observation points; the digital elevation model can use publicly available terrain products, such as SRTM or ASTER DEM, and be uniformly projected to the CGCS2000 coordinate system.

[0032] Data preprocessing projects all data onto the same coordinate system to ensure consistent latitude and longitude. Then, it performs gridding to create a regular grid in the target blank area and performs quality control.

[0033] Specifically, based on spatial consistency constraints, any two observation points are subjected to dual screening based on horizontal distance and elevation difference. If the horizontal distance between any two observation points is less than a preset threshold and the elevation difference exceeds 3-4 times the empirical standard deviation, they are identified as outlier pairs. Observations with higher elevation accuracy and better terrain fit are retained, while redundant data is removed. If an observation point does not meet the outlier pair criteria, it is considered a normal observation point. If the horizontal distance between normal observation points is less than 50 m, they are identified as duplicate or neighboring observation points. Observations with higher elevation accuracy are retained first. If the relative difference between the elevation standard deviation or root mean square error of multiple observation points does not exceed 10%-20%, a representative value is generated using a weighted average or median method, and the merged event source point, weight, and standard deviation are recorded. For observation points that do not meet the outlier pair criteria and whose horizontal distance is greater than 50 m, they are directly retained as valid observation data without further processing. Then, based on the distribution of observation points, cells in the grid without observation coverage are identified as cells to be filled and marked as blank areas for subsequent interpolation calculations.

[0034] Specifically, a preset threshold is used to define the horizontal distance threshold for spatial proximity relationships between observation points. Its value must match the spatial density of observation points in the survey area and the resolution of the digital elevation model (DEM). In practice, this preset threshold can be set to 10-200 m, preferably 30-90 m, to ensure effective consistency in identifying spatially adjacent observation points under a unified gridding condition. This threshold is not a fixed constant but can be reasonably adjusted according to the data accuracy and terrain complexity of different survey areas.

[0035] This invention makes full use of existing Bouguer gravity anomaly data and digital elevation model data, eliminating the need for additional on-site measurements and significantly reducing data acquisition costs and manpower input.

[0036] Step S2: After obtaining the standardized data from Step S1, the interpolation of the Bouguer gravity anomaly in the blank area is calculated. To ensure reliability, three interpolation algorithms are used in parallel, including the inverse distance weighted method, the ordinary kriging method, and the multiple regression model.

[0037] Specifically, the inverse distance weighting method estimates Bouguer anomalies of blank points using the inverse distance of neighboring observations as weights. The calculation formula is as follows:

[0038] ,

[0039] in, The Bouguer gravity anomaly is at the point to be estimated. For the neighborhood observations, d(L0, L) i ) represents the distance between two points, and p is the weighted power exponent (usually taken as 1.5 to 2);

[0040] Ordinary kriging, while considering spatial autocorrelation, also provides an estimated variance to measure the reliability of the results. First, a spatial correlation model of Bouguer anomalies is established using an empirical semivariance function, calculated as follows:

[0041] ,

[0042] in, The spatial distance between two points. Let N(h) be the average semivariance of the distance interval h, and N(h) be the number of sample points at the distance interval h. Using the semivariance model, parameters such as range, sill, and platform are determined. A Kriging equation system is then established to solve for the weighting coefficients of each observation point. The calculation formula is:

[0043] ,

[0044] in, Sample x i With x j The semivariance between the values ​​is given by μ, a Lagrange multiplier used to ensure the unbiasedness of the estimate. Kriging interpolation has good accuracy in regions with high data stationarity. Bouguer outliers at blank points are calculated as a weighted average. The calculation formula is:

[0045] ;

[0046] The multiple regression model establishes a functional relationship between Bouguer anomalies and elevation, slope, and terrain curvature to obtain Bouguer anomaly values ​​for blank areas. The physical deduction and mathematical fitting are performed, and the model formula is as follows:

[0047] ,

[0048] Where H is elevation, S is slope, and a i These are the fitting coefficients. As a residual, this method has the dual advantages of topographic constraints and statistical fitting. Under the condition of sparse measured data, it makes full use of the spatial geomorphological information in DEM and maintains the physical rationality and geological continuity of the anomaly field.

[0049] To quantitatively evaluate the effectiveness of the three interpolation methods in calculating Bouguer gravity anomalies in blank areas, this invention employs cross-validation to systematically verify the interpolation accuracy.

[0050] Specifically, individual samples are removed one by one from the known set of observation points, and the predicted value of the point is recalculated by interpolation using the remaining sample points. The predicted value of the point is then compared with the measured value to obtain the interpolation error sample set. This process is repeated for all sample points to obtain complete error statistics.

[0051] To comprehensively evaluate interpolation accuracy and model stability, this invention uses three statistical indicators: root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²).

[0052] Specifically, the root mean square error (RMSE) measures the overall deviation between predicted and measured values, and is defined as:

[0053] ,

[0054] in, Let i be the predicted Bouguer outlier value for the i-th sample point. The value represents the corresponding measured Bouguer outlier. The smaller the RMSE, the lower the overall prediction error and the higher the stability of the model.

[0055] Mean absolute error is used to measure the average deviation of interpolation results, and is defined as:

[0056] ,

[0057] The mean absolute error (MAE) is less sensitive to outliers and reflects the average fit of the model on most sample points. The smaller the value, the closer the interpolation result is to the measured distribution.

[0058] The coefficient of determination measures the ability of predicted values ​​to explain the variance of measured data, and is defined as:

[0059] ,

[0060] in, R² represents the average Bouguer anomaly across all observation points. The value of R² ranges from 0 to 1. The closer the value is to 1, the better the model fits and the more consistent the interpolation result is with the actual distribution.

[0061] This invention executes three interpolation algorithms in parallel and introduces a cross-validation method to quantitatively evaluate the interpolation results. The interpolation result corresponding to the case with the smallest RMSE, the smallest MAE, and the largest R² is selected as the optimal estimation field of Bouguer gravity anomaly, ensuring that the interpolation results of Bouguer gravity anomaly in the blank area have the smallest numerical deviation and the highest stability.

[0062] Step S3: After obtaining the interpolation result from step S2, two processing stages, mathematical analysis and filtering smoothing, are performed to eliminate problems such as local discontinuities, abrupt changes, or high-frequency noise in the Bouguer gravity anomaly obtained by mathematical fitting interpolation, so that it meets the physical characteristics of the Bouguer gravity anomaly field conforming to the continuity and analyzability of the potential field.

[0063] Specifically, the mathematical analysis is based on the fact that the gravitational potential field satisfies Laplace's equation ▽ 2 The property of V=0 can be expressed mathematically in the spatial domain using an analytical continuation method based on Laplace constraints, as follows:

[0064] ,

[0065] Where Ω represents the study region, which is discretized into a finite difference form:

[0066] ,

[0067] An iterative relaxation algorithm is used to iterate through the initial interpolation results. The calculation formula is as follows:

[0068] ,

[0069] in, Let k be the relaxation factor and k be the number of iterations.

[0070] The filtering and smoothing are performed by using a Gaussian low-pass filter in the frequency domain to eliminate high-frequency noise and smooth the interpolated Bouguer gravity anomaly field. Perform a two-dimensional fast Fourier transform to obtain the frequency domain spectral function. A Gaussian low-pass filter is constructed to suppress high-frequency fluctuations, and its expression is:

[0071] ,

[0072] Where, k L and k B These are the transverse and longitudinal wavenumbers in the frequency domain, respectively, k c Given the cutoff wavenumber, the filtered spectrum function The formula is:

[0073] ,

[0074] The filtered Bouguer gravity anomaly field is obtained by restoring the field to the spatial domain through inverse Fourier transform. :

[0075] ,

[0076] The Laplace smoothing operator is added to the spatial domain, and the calculation formula is as follows:

[0077] ,

[0078] in, This is the smoothed Bouguer gravity anomaly field, where The smoothing coefficient is 0.05 to 0.2.

[0079] This invention introduces a two-stage method of mathematical analysis and filtering smoothing, which ensures that the interpolation results conform to the continuity of the gravitational potential field and the constraints of the Laplace equation, and eliminates local discontinuities caused by data inhomogeneity or algorithm fitting. At the same time, a Gaussian low-pass filter is applied in the frequency domain to remove high-frequency noise, and the Laplace smoothing operator in the spatial domain is used to achieve natural transition and boundary smoothing of the numerical field.

[0080] Step S4: Superimpose the Bouguer gravity anomaly field obtained in step S3 with the gravity effect caused by terrain elevation to obtain a physically complete spatial gravity anomaly field. The calculation formula is as follows:

[0081] ,

[0082] in, The spatial gravity anomaly is located at the target point (L,B). To analyze and filter the smoothed Bouguer gravity anomaly, For the Bouger plate correction item, For terrain corrections, the formula for calculating Bouguer plate corrections is as follows:

[0083] ,

[0084] Where G=6.674 10 -11 N m 2 / kg 2 , Assuming a uniform crustal density ( =2.67 10 3 kg / m 3 H(L,B) represents the topographic elevation above the reference surface. Bouguer plate correction uses an infinitely horizontal, homogeneous plate to approximately compensate for differences in gravity baselines caused by varying elevations; the reference surface is typically chosen as mean sea level; the topographic correction term... The theoretical expression for accurately calculating the gravitational differences exerted by local terrain features such as peaks and valleys on observation points is as follows:

[0085] ,

[0086] Where Δh(L′,B′) represents the superelevation of the terrain element relative to the reference plane, r represents the distance between the observation point (L0,B0) and the center of the terrain element, θ is the angle between the line connecting the centroid of the element to the observation point and the vertical direction, and the integration region S is the terrain area to be considered; this integral is difficult to solve directly analytically and is usually implemented using numerical discretization methods:

[0087] Specifically, the DEM is discretized into regular grids or polygonal prisms. Each element is considered as a cuboid prism with a thickness equal to the element height difference. The total topographic correction is obtained by summing the contribution of each element using the analytical expression or numerical approximation of the prism's gravitational attraction. In the case of oceans or bodies of water, seawater density should be used. 水 =1.025 10 3 kg / m 3 (and corresponding sea surface reference processing) to avoid systematic bias.

[0088] After completing the two corrections and superimposing them on the smooth Bouguer anomaly, the final spatial gravity anomaly in the blank area is obtained. .

[0089] This invention utilizes the systematic processing of existing Bouguer gravity anomaly and terrain data, fully leverages the significant differences in the correlation between spatial gravity anomalies and Bouguer gravity anomalies with terrain, establishes a physical constraint model suitable for gravity anomaly estimation in terrain-variable areas, and successfully reconstructs gravity anomalies in blank areas within complex terrain and sparse data regions based on physical laws.

[0090] To verify the effectiveness of the method for obtaining gravity anomaly blank areas based on potential field characteristics described in this invention, a specific implementation example and its effect are presented below, using an actual area on the western edge of a basin. The specific steps are as follows:

[0091] Step S1, as follows Figure 3 As shown in this embodiment of the invention, Bouguer gravity anomaly data and digital elevation model data are collected within the western edge of a basin, ranging from 101°–105°E to 29°–32°N. The Bouguer gravity anomaly data (in .grd and .dat formats) are from the WGM2012 global model, and the digital elevation model data are from the Shuttle RadarTopography Mission (SRTM). The digital elevation data is projected onto a unified latitude and longitude datum under the CGCS2000 coordinate system, and spatial matching and coordinate transformation are performed. Subsequently, gridding is performed to establish a unified Bouguer anomaly and elevation dataset. Based on the distribution of observation points, areas not covered by gravity anomaly values ​​are identified and marked as blank areas. The blank areas range from 102.5°–104°E to 30°–31°N.

[0092] Step S2, as follows Figure 4 As shown in the embodiments of the present invention, gravity anomalies in the blank area are calculated using the inverse distance weighting method, the ordinary kriging method, and the multiple regression model, respectively.

[0093] Specifically, for the inverse distance weighted method, a weighted power exponent is set. The results reflect the smooth field characteristics under the control of nearby measuring points;

[0094] For the ordinary kriging method, the model parameters are obtained by fitting the experimental semivariogram function: sill value C0 = 2.8mGal 2 The table surface value C = 5.6 mGal 2 The range a = 25 km is used to solve the weight matrix using a spherical model;

[0095] For the selection of elevation in the multiple regression model ,slope and terrain curvature Bouguer anomaly is the independent variable and the dependent variable, calculated using the following formula:

[0096] ,

[0097] Substituting these values ​​into the formula, we obtain the regression coefficients a1 = -0.042, a2 ​​= 0.018, and a3 = -0.006.

[0098] The root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²) were calculated using cross-validation. The results are as follows:

[0099] Inverse distance weighted method: RMSE=4.22 mGal, MAE=3.35 mGal, R²=0.80;

[0100] Ordinary Kriging: RMSE = 3.84 mGal, MAE = 3.01 mGal, R² = 0.89;

[0101] Multiple regression model: RMSE=8.13 mGal, MAE=6.27 mGal, R²=0.69.

[0102] In summary, Kriging interpolation outperforms the other two algorithms in both statistical accuracy and physical consistency. Therefore, the Kriging interpolation result is selected as the optimal estimation field for Bouguer gravity anomalies in the blank area.

[0103] Step S3, as follows Figure 5 As shown, in this embodiment of the invention, the Bouguer gravity anomaly field obtained by ordinary Kriging interpolation in step S2 is subjected to mathematical analysis and filtering smoothing.

[0104] Specifically, firstly, an analytical continuation is performed using an iterative relaxation method, with a relaxation factor of λ=0.2 in each step, for a total of 200 iterations, until the field function satisfies the Laplace condition. This process eliminates local discontinuities and artifacts. Subsequently, a Gaussian low-pass filter is applied in the frequency domain to effectively retain 90% of the power spectrum energy. After inverse transformation back to the spatial domain, a smooth Bouguer anomaly field with effectively suppressed high-frequency noise is obtained. To further improve edge transitions, a Laplace smoothing operator is applied in the spatial domain with a smoothing coefficient of α=0.1. After processing, the gradient of the Bouguer anomaly field is continuous, and the anomaly morphology is consistent with the terrain changes, possessing physical rationality.

[0105] Step S4, as follows Figure 6 , 7 As shown, in this embodiment of the invention, the smoothed Bouguer gravity anomaly is superimposed with the terrain gravity effect to obtain the spatial gravity anomaly in the blank area.

[0106] Specifically, the Bouguer board correction item is adopted. ,in , Using the DEM elevation as H(L,B), the topographic correction term was calculated using the discrete integral method of polygonal prisms, and finally the spatial gravity anomaly field of the western edge of the Sichuan Basin was obtained.

[0107] The embodiments of the present invention described above do not constitute a limitation on the scope of protection of the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the claims of the present invention.

Claims

1. A method for obtaining gravity anomaly blank area data based on potential field characteristics, characterized in that, Includes the following steps: S1. Obtain gravity and topographic data of the exploration area and its surrounding area, preprocess and grid the data to generate standard grid data and mark blank areas; S2. Perform interpolation calculations on the data and conduct quantitative evaluation to obtain reliable Bouguer gravity anomaly values ​​in the blank area; S3. Perform mathematical analysis and filtering smoothing on the interpolation results to eliminate noise and obtain Bouguer gravity anomaly values ​​that conform to the physical continuity and smoothness of the gravity field. The mathematical analysis is based on the fact that the gravitational potential field satisfies the Laplace equation ▽ 2 The property of V=0 can be expressed mathematically in the spatial domain using an analytical continuation method based on Laplace constraints, as follows: ， Where Ω represents the study region, which is discretized into a finite difference form: ， An iterative relaxation algorithm is used to iterate through the initial interpolation results. The calculation formula is as follows: ， in, Let k be the relaxation factor and k be the number of iterations. The filtering and smoothing are performed by using a Gaussian low-pass filter in the frequency domain to eliminate high-frequency noise and smooth the interpolated Bouguer gravity anomaly field. Perform a two-dimensional fast Fourier transform to obtain the frequency domain spectral function. A Gaussian low-pass filter is constructed to suppress high-frequency fluctuations, and its expression is: ， Where, k L and k B These are the transverse and longitudinal wavenumbers in the frequency domain, respectively, k c Given the cutoff wavenumber, the filtered spectrum function The formula is: ， By restoring the spatial domain through inverse Fourier transform, the smoothed Bouguer gravity anomaly field is obtained. : ， The Laplace smoothing operator is added to the spatial domain, and the calculation formula is as follows: ， in, The smoothed Bouguer gravity anomaly field. For smoothing coefficients; S4. The processed Bouguer gravity anomaly field is superimposed with the topographic gravity effect to obtain the complete spatial gravity anomaly field.

2. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 1, characterized in that, In step S1, the gravity and topographic data include the slope of the explored area, geological information on topographic curvature, high-precision digital elevation model data, and Bouguer gravity anomaly data of the area adjacent to the blank area. The Bouguer gravity anomaly data includes latitude and longitude, Bouguer gravity anomaly value, and elevation information of the observation point.

3. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 1, characterized in that, In step S1, data preprocessing projects all data onto the same coordinate system and performs gridding to establish a regular network in the target blank area, while simultaneously performing quality control, specifically including: Based on spatial consistency constraints, any two observation points are subjected to dual screening based on horizontal distance and elevation difference. If the horizontal distance between any two observation points is less than a preset threshold (10-200 m), and the elevation difference exceeds 3-4 times the empirical standard deviation, they are identified as outlier pairs. Observations with higher elevation accuracy and better terrain fit are retained, while redundant data is discarded. If an observation point does not meet the outlier pair criteria, it is considered a normal observation point. If the horizontal distance between normal observation points is less than 50 m, they are identified as duplicate or neighboring observation points. Observations with higher elevation accuracy are retained first. If the relative difference between the elevation standard deviations or root mean square errors of multiple observation points does not exceed 10%-20%, a representative value is generated using a weighted average or median method, and the merged event source point, weight, and standard deviation are recorded. For observations that do not meet the outlier pair criteria and whose horizontal distance is greater than 50 m... The observation points of m are directly retained as valid observation data without processing; then, according to the distribution of observation points, the cells in the grid without observation coverage are marked as cells to be filled and marked as blank areas for subsequent interpolation calculations.

4. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 1, characterized in that, In step S2, multiple interpolation algorithms are combined for parallel computation, and cross-validation is used for quantitative evaluation. The multiple interpolation algorithms include inverse distance weighting, ordinary kriging, and multiple regression models.

5. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 4, characterized in that, The inverse distance weighting method estimates Bouguer anomalies of blank points using the inverse distance of neighboring observations as weights. The calculation formula is as follows: ， in, The Bouguer gravity anomaly is at the point to be estimated. For the neighborhood observations, d(L0, L) i ) represents the distance between two points, and p represents the weighted exponent; The ordinary kriging method establishes a spatial correlation model of Bouguer anomalies using an empirical semivariance function, and the calculation formula is as follows: ， in, The spatial distance between two points. Let N(h) be the average semivariance of the distance interval h, and N(h) be the number of sample pairs at the distance interval h. The range, sump, and platform parameters are determined using a semivariance model, and the weighting coefficients for each observation point are solved by establishing the Kriging equations. The calculation formula is: ， in, Sample x i With x j The semivariance values ​​between the points are calculated as a weighted average, where μ is the Lagrange multiplier, and the Bouguer outliers of the blank points are calculated as a weighted average. The calculation formula is: ； The multivariate regression model establishes a functional relationship between Bouguer anomalies and elevation, slope, and terrain curvature, and is used to measure Bouguer anomaly values ​​in blank areas. The physical deduction and mathematical fitting are performed, and the model formula is as follows: ， Where H is elevation, S is slope, and a i These are the fitting coefficients. This is the residual.

6. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 5, characterized in that, In step S2, cross-validation is used to evaluate the accuracy of the three interpolation methods, specifically including: In the known set of observation points, individual samples are removed one by one, and the predicted value of the point is recalculated by interpolation using the remaining sample points. The predicted value of the point is then compared with the measured value to obtain the interpolation error sample set. After traversing all the sample points, the complete error statistics are obtained. The optimal interpolation model is selected by comprehensively quantifying three statistical indicators: root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²). 2 .

7. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 1, characterized in that, In step S4, the complete spatial gravity anomaly field is defined by the following formula: ， in, The spatial gravity anomaly is located at the target point (L,B). To analyze and filter the smoothed Bouguer gravity anomaly, For the Bouger plate correction item, This is a terrain correction item.

8. The method for obtaining gravity anomaly blank area data based on potential field characteristics according to claim 7, characterized in that, The formula for calculating the Bouger plate correction term is as follows: ， Where G=6.674 10 -11 N m 2 / kg 2 Uniform crustal density =2.67 10 3 kg / m 3 H(L,B) is the topographic elevation above the reference surface; The terrain correction item The theoretical expression is: ， Where Δh(L′,B′) represents the superelevation of the terrain unit relative to the reference plane, r represents the distance between the observation point (L0,B0) and the center of the terrain micro-element, θ is the angle between the line connecting the centroid of the micro-element to the observation point and the vertical direction, and the integration region S is the terrain range to be considered.

Images