Deep coverage area skarn type iron-rich ore target detection method based on gravity and magnetic data

Abstract

The invention relates to the technical field of iron ore detection, and discloses a deep coverage area skarn type iron-rich ore target detection method based on gravity and magnetic data. The method comprises the steps of collecting and preprocessing gravity and magnetic measurement data of a target area to obtain a standardized gravity and magnetic data set; based on the geological background of the analysis area, the structure and lithology distribution characteristics are extracted, and the metallogenic conditions are matched to screen the potential metallogenic target area. And performing multi-scale decomposition on the gravity and magnetic data of the target region, obtaining different depth level abnormal features, constructing a three-dimensional geologic structure model, and obtaining density and magnetic susceptibility distribution based on model inversion. Accordingly, the ore body is recognized, spatial distribution features are extracted, continuity and scale of the ore body are determined through spatial correlation analysis, and a final prediction result is output after geological constraint condition optimization is combined. The method is suitable for a complex environment with signal shielding and abnormal superposition in a deep coverage area, deep ore body information can be accurately captured, and a feasible path is provided for skarn type rich iron ore detection.

Description

Technical Field

[0001] This invention relates to the field of iron ore detection technology, specifically a method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data. Background Technology

[0002] my country's dependence on imported iron ore has long been high. Exploration and development of skarn-type rich iron ore deposits in deeply overburdened areas are of practical significance in alleviating the contradiction between resource supply and demand. Skarn-type rich iron ore deposits mostly form at the contact zone between the rock mass and carbonate rocks. However, these deeply overburdened areas are generally covered by thick Cenozoic strata. The geophysical signals of underground geological bodies weaken significantly after long-distance propagation, and mineral anomalies are often masked by large-scale magnetic rock anomalies, making identification and extraction extremely difficult.

[0003] Traditional detection methods often rely on single gravity and magnetic data or conventional processing techniques, which have many limitations. In the data processing stage, gravity and magnetic data often lack a unified standardized process, and the mixing of data from different sources and with varying precision leads to biases in subsequent analysis results. In regional geological background analysis, the extraction of structural features and lithological distribution is often fragmented, failing to form an organic correlation and making it difficult to accurately match the mineralization conditions of skarn-type iron deposits. The screening of potential mineralization target areas often relies on empirical judgment, lacking a systematic integration of gravity and magnetic data with mineralization patterns.

[0004] The core challenge in deep orebody detection in heavily overburdened areas lies in the separation of anomalies between geological bodies at different depths. Traditional analytical extrapolation and moving average methods are insufficient to effectively separate shallow interferences from deep mineral-induced anomalies. For instance, Fourier analysis can only obtain the overall signal spectrum and cannot capture local depth features, leading to biases in the judgment of orebody depth and morphology. Furthermore, the integration of 3D geological modeling and gravity / magnetic inversion is insufficient, and the lack of effective quality constraints during the inversion process easily leads to multiple solutions, resulting in discrepancies between the calculated density and magnetic susceptibility distributions and the actual geological conditions.

[0005] Existing methods often stop at morphological description after ore body identification, lacking quantitative means to analyze the spatial continuity and scale of ore bodies, making it difficult to determine whether dispersed anomalies belong to the same ore body or ore zone. In some ore clusters, the anomaly areas in the central part of the region are small in scale and weak in intensity due to ultra-deep overburden, which has prevented mineral exploration breakthroughs for a long time, confirming the limitations of traditional methods in the detection of skarn-type rich iron ore in deeply overburdened areas. Summary of the Invention

[0006] The purpose of this invention is to provide a method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, this invention provides a method for detecting skarn-type rich iron ore targets in deeply overburdened areas based on gravity and magnetic data, the method comprising: Gravity and magnetic measurement data of the target area are collected, and the gravity and magnetic measurement data are preprocessed to obtain a standardized gravity and magnetic dataset. Based on the standardized gravity and magnetic dataset, regional geological background analysis was performed to extract the structural features and lithological distribution features of the target area. Based on the aforementioned structural and lithological distribution characteristics, the mineralization conditions of skarn-type rich iron deposits are matched to screen potential mineralization target areas. The gravity and magnetic data of the potential mineralization target area are decomposed at multiple scales to obtain gravity and magnetic anomaly features at different depth levels. Three-dimensional geological modeling is performed by combining the gravity and magnetic anomaly characteristics at different depth levels to construct a geological structure model of the target area; Based on the geological structure model, gravity and magnetic data inversion is performed to calculate the density distribution and magnetic susceptibility distribution of the target area. Based on the density distribution and magnetic susceptibility distribution, the ore bodies of skarn-type rich iron ore are identified, and the spatial distribution characteristics of the ore bodies are extracted. Spatial correlation analysis was performed on the spatial distribution characteristics of the ore body to determine its continuity and scale. The spatial distribution characteristics of the ore body are optimized based on geological constraints, and the final ore body prediction results are output.

[0008] Preferably, the preprocessing of the gravity measurement data and magnetic measurement data to obtain a standardized gravity and magnetic dataset includes the following steps: The gravity measurement data are subjected to terrain correction and regional field correction to eliminate the influence of terrain undulation and regional background field. The magnetic measurement data are subjected to diurnal variation correction and magnetic declination correction to eliminate the influence of geomagnetic diurnal variation and magnetic declination; A gridded interpolation method was used to spatially interpolate the corrected gravity and magnetic measurement data to form regular grid data. The regular grid data is normalized so that the gravity measurement data and magnetic measurement data have the same numerical range, forming the standardized gravity and magnetic dataset.

[0009] Preferably, the step of performing regional geological background analysis based on the standardized gravity and magnetic dataset to extract the tectonic features and lithological distribution features of the target area includes the following steps: Gradient calculations are performed on the standardized gravity and magnetic dataset to extract gravity gradient and magnetic gradient features; Fracture structure identification is performed based on the gravity gradient and magnetic gradient characteristics to determine the fracture distribution in the target area; The distribution of the faults was verified by combining known geological data, and the main ore-controlling structures were screened. Based on the standardized gravity and magnetic dataset, lithological classification is performed to divide the target region into lithological units; Based on the superposition relationship between the lithological units and the main ore-controlling structures, favorable ore-forming lithological assemblages are extracted.

[0010] Preferably, the step of matching the mineralization conditions of skarn-type rich iron ore deposits based on the structural features and lithological distribution characteristics, and screening potential mineralization target areas, includes the following steps: Establish a metallogenic geological model for skarn-type rich iron deposits, and define the ore-controlling structural types and favorable lithological combinations; The main ore-controlling structures and favorable ore-forming lithological combinations are matched with the ore-forming geological model to calculate the ore-forming matching degree; A mineralization matching threshold is set, and areas with a matching degree higher than the threshold are selected as potential mineralization target areas.

[0011] Preferably, the step of performing multi-scale decomposition of the gravity and magnetic data of the potential mineralization target area to obtain gravity and magnetic anomaly features at different depth levels includes the following steps: The gravity and magnetic measurement data of the potential mineralization target area were decomposed into multiple scales using wavelet transform. Wavelet coefficients at different scales were extracted, corresponding to the gravity and magnetic anomaly characteristics of shallow, middle and deep regions, respectively; The wavelet coefficients at different scales are reconstructed to obtain gravity and magnetic anomaly distribution maps at different depth levels.

[0012] Preferably, the step of combining the gravity and magnetic anomaly features at different depth levels to perform three-dimensional geological modeling and construct a geological structure model of the target area includes the following steps: Based on the gravity and magnetic anomaly distribution maps at different depth levels, geological interfaces are identified to delineate the boundaries of strata and rock masses. By combining borehole data and geophysical logging data to constrain the geological interface, the modeling accuracy is improved. A three-dimensional interpolation algorithm is used to construct the density volume and magnetic susceptibility volume of the target area to form the geological structure model.

[0013] Preferably, the step of performing gravity and magnetic data inversion based on the geological structure model to calculate the density distribution and magnetic susceptibility distribution of the target area includes the following steps: The geological structure model was iteratively optimized using a constrained inversion algorithm to fit measured gravity and magnetic measurement data. Calculate the inverted density distribution and magnetic susceptibility distribution, and evaluate the reliability of the inversion results; Adjust the constraints based on the inversion error distribution to improve inversion accuracy.

[0014] Preferably, the step of identifying skarn-type rich iron ore bodies based on the density distribution and magnetic susceptibility distribution, and extracting the spatial distribution characteristics of the ore bodies, includes the following steps: Set density and magnetic susceptibility thresholds to screen out anomalous areas that conform to the characteristics of skarn-type iron-rich minerals; Spatial clustering analysis was performed on the aforementioned abnormal region to divide it into independent ore body units; The spatial location, morphology, and scale characteristics of the independent ore body units are extracted to form the spatial distribution characteristics of the ore body.

[0015] Preferably, the step of performing spatial correlation analysis on the spatial distribution characteristics of the ore body to determine its continuity and scale includes the following steps: Calculate the spatial distance and physical property similarity between adjacent ore body units; The connectivity of the ore body is assessed based on the spatial distance and physical property similarity. Merging highly connected ore body units forms a continuous ore body distribution.

[0016] Preferably, the optimization of the spatial distribution characteristics of the ore body based on geological constraints, and the output of the final ore body prediction result, includes the following steps: The spatial distribution characteristics of the ore body are corrected based on known mineralization information and geological patterns; Eliminate abnormal areas that do not conform to geological laws and optimize the ore body boundaries; Generate orebody prediction maps and mark the location, depth, and size of the orebody.

[0017] Compared with the prior art, the beneficial effects of the present invention are: By standardizing and preprocessing gravity and magnetic measurement data, raw data from different sources and with varying precision can be integrated, eliminating the influence of systematic errors and environmental interference. This gives the gravity and magnetic dataset a unified data benchmark and comparable properties, providing a consistent foundation for all subsequent geological analysis stages. The standardized dataset can more accurately reflect the differences in the physical properties of underground geological bodies, avoiding feature extraction biases caused by data heterogeneity.

[0018] Regional geological background analysis based on standardized data can extract structural features and lithological distribution features in a linked manner, clearly presenting their spatial relationship. The formation of skarn-type rich iron ore deposits is closely related to tectonic fracture zones and intrusive contact zones. This linked analysis model can accurately capture the combined characteristics of key geological elements for mineralization, making the screening of potential mineralization target areas more consistent with mineralization patterns and reducing the omission or misjudgment of target areas due to judgment based on a single feature.

[0019] The application of multi-scale decomposition technology enables in-depth hierarchical analysis of gravity and magnetic anomalies, effectively separating shallow surface soil interference, mid-level rock mass anomalies, and deep ore body anomalies to obtain anomaly features at different depth levels. This hierarchical analysis capability is suitable for complex scenarios with signal superposition in deeply covered areas. For example, for superimposed anomalies generated by shallow magnetic rock masses and deep iron ore bodies, the scale decomposition can clearly identify the anomaly morphology and distribution range of each, providing direct evidence for locating deep ore bodies.

[0020] A three-dimensional geological structure model constructed by combining multi-depth anomaly features can transform planar gravity and magnetic information into a three-dimensional geological framework, intuitively presenting the spatial relationship between strata, rock masses, and potential ore bodies. Gravity and magnetic inversion based on this model can fully utilize the spatial constraints of geological bodies, reduce inversion ambiguity, and make the calculated results of density and magnetic susceptibility distribution more consistent with actual geological conditions. Skarn-type rich iron ore has significantly higher magnetic susceptibility and specific density characteristics than the surrounding rocks. This precise physical property distribution data can directly provide clear identifiers for ore body identification and accurately delineate the spatial distribution morphology of the ore body.

[0021] Spatial correlation analysis can verify the correlation of dispersed orebody anomalies. By analyzing the consistency of the strike, occurrence, and physical properties of different anomalies, it can be determined whether they belong to the same orebody or continuous ore zone, thereby clarifying the extension range and scale characteristics of the orebody. Finally, combined with the optimization process of geological constraints, regional stratigraphic age, tectonic evolution history, and known borehole information can be incorporated to correct the initially identified orebody distribution characteristics, making the final output prediction results more consistent with geological evolution laws and mineralization reality. The entire method forms a complete logical chain from data processing to result output, adapting to the complex needs of skarn-type rich iron ore exploration in deeply overburdened areas, and can effectively capture information on deep concealed orebodies. Attached Figure Description

[0022] Figure 1 A comprehensive map of the detection results; Figure 2 A flowchart for gravity and magnetic data preprocessing and standardization; Figure 3 A flowchart for regional geological background analysis and feature extraction; Figure 4 This is a multi-scale anomaly analysis diagram. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] Please see Figure 1 This invention provides a method for detecting skarn-type rich iron ore targets in deeply covered areas based on gravity and magnetic data. The method includes: deploying gravity and magnetic measuring points in the target area; collecting raw data using high-precision gravimeters and magnetometers; ensuring uniform distribution of measuring points throughout the area during data collection to avoid data gaps; preprocessing the raw gravity and magnetic data after collection, including removing jump point errors and instrument drift effects to form a preliminary dataset; normalizing the data using a standardization algorithm to ensure gravity and magnetic values ​​are of the same order of magnitude, facilitating subsequent fusion analysis, and finally generating a standardized gravity and magnetic dataset; conducting regional geological background analysis based on this dataset, identifying fault zones and lithological boundaries by calculating gravity and magnetic gradient fields, and extracting structural features such as fault strike and lithological distribution features such as the extent of igneous bodies; matching the extracted features with skarn-type rich iron ore mineralization models, with matching factors including structural ore-controlling conditions and lithological assemblage types; and selecting potential mineralization target areas based on the matching degree score. Within the target area, wavelet transform or multi-resolution methods are applied to decompose gravity and magnetic data at multiple scales, separating shallow, mid-deep, and deep anomalous signals to obtain gravity and magnetic anomaly characteristics at different depth levels. Combining borehole logs and geophysical well logging data, a three-dimensional geological model is constructed, building a geological structure model that includes stratigraphic interfaces and rock mass morphology. Based on this model, constrained inversion algorithms such as least squares inversion are used to calculate the density and magnetic susceptibility distributions of the target area. Iterative optimization is employed during the inversion process to reduce fitting errors. Based on the inversion results, density and magnetic susceptibility thresholds are set to identify orebody anomaly zones, extracting the spatial location, morphology, and scale characteristics of the orebody. Spatial correlation analysis is performed on orebody units to assess the connectivity of adjacent units, merging continuous orebody units and determining the overall scale. Geological constraints such as mineralization regularities and known mineralization points are introduced to optimize orebody boundaries, generate orebody prediction maps, and output the final results.

[0025] Example 1: See Figure 2 In the topographic correction stage of gravity measurement data preprocessing, digital elevation model (DEM) data is required. This DEM data originates from a 1:50,000 scale DEM product provided by the surveying and mapping department. The coordinate system of the gravity measurement points and the DEM coordinate system both adopt the National Geodetic Coordinate System, and the elevation datum uses the National Elevation Datum. The elevation correction value for each gravity measurement point in the topographic correction calculation is calculated using the Bouguer plate model formula, which is:

[0026] in: This represents the terrain correction value, in milligal. It is the gravitational constant, with a value of ; It is the density of surface rocks, with values ​​ranging from... ; It is the actual elevation of the measuring point; This refers to the regional reference elevation. The mean sea level of the survey area is used as the reference surface during the calculation, and the elevation difference at each measuring point is obtained through interpolation using a digital elevation model. After topographic correction, regional field correction is performed. This regional field correction employs a frequency domain filtering method. A high-pass filter is designed with a cutoff wavelength of 15 kilometers. Fourier transform is used to convert the gravity data to the frequency domain, and after applying the filter function, an inverse transform is performed back to the spatial domain, achieving separation of the regional field from local anomalies.

[0027] In the preprocessing of magnetic measurement data, diurnal variation correction requires the establishment of a reference station. A high-precision proton magnetometer is selected as the reference station, with a sampling interval of 10 seconds to continuously record the diurnal variation of the total geomagnetic field intensity. The rover magnetometer is synchronized with the reference station in time, achieving a synchronization accuracy of 0.1 seconds. The diurnal variation correction is calculated by linear interpolation of the reference station value corresponding to each rover sampling time, subtracting the diurnal variation influence from the rover's raw data. Magnetic declination correction uses the International Geomagnetic Reference Field Model. The latitude and longitude coordinates of the survey area and the measurement date are input to calculate the magnetic declination value for the survey area. The magnetic declination correction formula is:

[0028] in: This is the corrected magnetic force value. It measures the magnetic force value. It is the magnetic declination value. This is the horizontal component. The corrected magnetic data is uniformly converted to the geographic North Pole direction.

[0029] The gridded interpolation process employed the Kriging interpolation algorithm. A spherical model was chosen as the variogram model in the Kriging interpolation algorithm, and the range parameter of the spherical model was determined based on the spatial correlation of the measurement points. Through experimental variogram fitting, the range value was found to be 800 meters. During interpolation, the search radius was set to 1.5 times the range, i.e., 1200 meters. The interpolation output grid size was set to 50 meters × 50 meters, with the grid direction parallel to the boundary of the survey area. Gravity and magnetic data were gridded separately, generating two regular grid data files in ASCII grid format, containing a grid header file and a data matrix.

[0030] The normalization process uses a minimum-maximum scaling algorithm, the formula of which is:

[0031] in: These are normalized values. These are the original grid data values. It is the minimum value of the grid data. This represents the maximum value of the grid data. Gravity and magnetic data are normalized separately so that the values ​​of both types of data are between 0 and 1. The normalized grid data is stored according to the coordinates of the measurement points, forming a standardized gravity and magnetic dataset. The dataset contains three fields: latitude and longitude coordinates, normalized gravity value, and normalized magnetic value.

[0032] During terrain correction, special attention must be paid to the impact of steep terrain. For areas with elevation differences greater than 300 meters, a precise terrain correction method is adopted, dividing the digital elevation model into smaller computational units. The terrain gravitational effect is calculated separately for each unit, and then the results are integrated and summed. For regional field correction, a Butterworth filter is used with an 8th-order filter to ensure a steep transition zone. In magnetic declination correction, the 13th generation International Geomagnetic Reference Field model is used, with model coefficients updated annually to ensure the accuracy of magnetic declination calculations. During Kriging interpolation, for sparse data areas, the number of search points is increased to a minimum of 8 to ensure interpolation stability. Before normalization, outlier removal is performed on the grid data using a 3x standard deviation rule to avoid extreme values ​​affecting the normalization results. Gravity measurement points are laid out in a regular grid pattern, with the point spacing determined according to the detection depth: 100 meters for shallow detection points and 200 meters for deep detection points. Each measurement point is measured three times repeatedly, and the average value is taken as the final observation value. Magnetic and gravity measurements were conducted simultaneously using an integrated gravity and magnetic measurement system to ensure spatiotemporal consistency of the data. Meteorological parameters, including temperature, air pressure, and humidity, were recorded during the measurements for instrument reading correction. Data preprocessing software employed a professional geophysical processing system, automating batch processing to reduce human error. For topographic correction calculations, surface rock density values ​​were based on regional geological survey results, with zoned density values ​​used for different lithological zones: 2.60 g / cm³ for sedimentary rock zones and 2.75 g / cm³ for igneous rock zones. A 30m × 30m resolution digital elevation model was used to ensure detailed topographic representation. In regional field correction, the cutoff wavelength was selected based on regional tectonic scale analysis, with the optimal cutoff wavelength determined through power spectrum analysis. During magnetic declination correction, normal field correction was performed simultaneously to remove the influence of the geomagnetic field gradient. After gridded interpolation, data quality was checked using methods including contour smoothness analysis and cross-validation to ensure interpolation accuracy. The normalized dataset underwent consistency testing using methods including correlation analysis and statistical distribution testing.

[0033] Gravimeter calibration uses an absolute gravity point as the benchmark, and instrument comparisons are performed periodically. Magnetometer calibration uses a standard magnetic generator, calibrated once per working day. Measurement data is transmitted to the processing center in real time for quality monitoring. Each step in the preprocessing process generates a quality report, including data processing parameters, correction statistics, and accuracy assessment. Standardized gravity and magnetic datasets are stored in a standard format for easy retrieval by subsequent analysis modules. Dataset management uses a database system for version control and tracking. Topographic correction accuracy is assessed through repeated measurements at checkpoints located in areas with typical topographic features. The regional field correction effect is verified through residual analysis; the residuals should exhibit a random distribution. Magnetic declination correction accuracy is verified through comparison with benchmark stations. Grid interpolation accuracy is evaluated through cross-validation, reserving 20% ​​of measurement points as verification points. Normalization processing effectiveness is checked using a data distribution histogram to ensure a reasonable data distribution. Tidal correction is considered during gravity measurement data acquisition, calculated using an Earth tidal model. Instrument drift correction is achieved through repeated observations of the base point, and a linear model is used for drift correction. Vehicle interference is considered in magnetic measurements, and a safe distance is maintained from vehicles during measurement. The diurnal variation correction base station is located far from electromagnetic interference sources to ensure data quality. The gridded interpolation algorithm has been optimized, and a block processing strategy is used when handling large amounts of data. Data standardization is performed before normalization to ensure that the data conforms to a normal distribution. Preprocessing results output include data files and quality report files for use in subsequent processing stages. The measurement coordinate system adopts the 2000 National Geodetic Coordinate System, and the elevation datum adopts the 1985 National Elevation Datum. Coordinated Universal Time (UTC) is used for time recording to avoid time zone confusion. Data storage uses both binary and text formats for dual backup to ensure data security. The processing log records the parameter settings and processing results for each step in detail, facilitating traceability and review.

[0034] Example 2: See Figure 3The regional geological background analysis begins with a standardized gravity and magnetic dataset, which includes normalized gravity and magnetic grid data in a 50m × 50m regular grid format. Gradient calculation employs the horizontal gradient modulus method. The Sobel operator, a 3×3 convolution kernel, is applied to the gravity grid data to calculate the horizontal gradient, performing convolution operations in both the east-west and north-south directions to obtain the gravity horizontal gradient components. Magnetic gradient calculation first involves polarization processing, using the dip and deflection parameters of the geomagnetic field in the survey area to convert the magnetic anomaly into an equivalent polar magnetic anomaly. Then, the vertical gradient is calculated using a frequency-domain vertical first-order derivative operator. The gradient calculation results generate a gravity horizontal gradient modulus map and a magnetic vertical gradient map, stored in a grid format. Fault structure identification is based on gradient features; high-value bands of the gravity horizontal gradient modulus correspond to density interfaces, and high-value bands of the magnetic vertical gradient correspond to magnetic susceptibility interfaces. The fracture identification algorithm employs the Canny edge detection algorithm, which includes four steps: Gaussian filtering, gradient magnitude calculation, non-maximum suppression, and dual-threshold detection. The Gaussian filtering standard deviation is set to two grid cells, non-maximum suppression uses a 3×3 neighborhood, the high threshold is set to the 70th percentile of the gradient value, and the low threshold is set to the 30th percentile. Identified fracture lines are vectorized, and their attributes include strike, dip, and length. Known geological data includes 1:50,000 geological maps and borehole lithological records. The digitized geological maps are overlaid with the fracture lines to verify their geological significance. The selection criteria for major ore-controlling structures include fracture scale, extension depth, and relationship with known mineralization. Faults extending greater than 2 kilometers and with a cutting depth exceeding 1000 meters are selected as major ore-controlling structures.

[0035] Lithological distribution features were extracted using a supervised classification method, with training samples drawn from known lithological outcrops on geological maps. Feature variables included normalized gravity, normalized magnetic, gravity gradient, and magnetic gradient. A support vector machine (SVM) classifier was used, with a radial basis function kernel, a penalty parameter C of 10, and a gamma parameter of 0.1. The classification results divided the survey area into lithological units such as granite, limestone, sandstone, and basic rocks, generating polygonal vector files for each unit. Spatial overlay analysis was performed between the lithological units and the main ore-controlling structures. This overlay analysis employed the overlay operation in a geographic information system (GIS), extracting the lithological contact zone within 500 meters on either side of the fault zone as a favorable ore-forming lithological assemblage. In the ore-forming condition matching stage, a skarn-type iron-rich ore-forming geological model was established. The ore-controlling structural type was defined as the intersection of a NE-trending fault and a NW-trending fault, and the favorable lithological assemblage was the contact zone between limestone and granite. The matching degree was calculated using a weighted superposition analysis method, with a structural matching degree weight of 0.6 and a lithological matching degree weight of 0.4. Structural matching degree was calculated based on the fault intersection angle and distance. For every 10-degree deviation of the intersection angle from the ideal 45-degree angle, the matching degree decreased by 0.1; for every 100-meter increase in distance, the matching degree decreased by 0.05. Lithological matching degree was calculated based on the similarity between the lithological combination and the ideal combination. A complete match scored 1 point, a partial match scored 0.5 points, and no match scored 0 points. The matching degree threshold was set at 0.7. Areas exceeding the threshold were used to generate potential mineralization target polygons, with a 500-meter buffer zone extended along the target area boundaries.

[0036] During gradient calculation, the gravity horizontal gradient is calculated using a 3×3 Sobel operator convolution kernel, with the east-west kernel set to [-1,0,1;-2,0,2;-1,0,1] and the north-south kernel set to [1,2,1;0,0,0;-1,-2,-1]. The magnetic vertical gradient is calculated in the frequency domain, multiplied by a function of wavenumber k after Fourier transform, and then inversely transformed back to the spatial domain. For fracture identification, the Canny algorithm's dual-threshold detection retains edge points above the high threshold, discards those below the low threshold, and retains edge points in between only if connected to the high-threshold edge. Verification using known geological data employs spatial consistency analysis, calculating the average distance between the identified fracture and the fracture on the geological map; fractures less than 100 meters away are considered verified. For lithology classification, each lithology category has at least 100 grid points, evenly distributed across the survey area. The support vector machine classifier is trained using 5-fold cross-validation, achieving a classification accuracy of 85% or higher before deployment. Morphological closure operations are performed on the boundaries of lithological units to eliminate small cavities and burrs. When extracting favorable mineralized lithological assemblages, in addition to the contact zones on both sides of the fault zone, favorable locations such as the depressions at the top interface of the rock mass and the turning points of the strata are also considered.

[0037] The metallogenic geological model parameters are set based on the characteristics of typical regional deposits, with NE-trending faults striking between 40-50 degrees and NW-trending faults between 310-320 degrees. The matching degree calculation is automated, inputting structural and lithological vector data and outputting a matching degree raster map. After target area screening, each target area undergoes manual review based on geochemical anomalies and remote sensing alteration information. Target area results are stored as a vector layer, with attribute tables including target area number, area, average matching degree, and main ore-controlling factors. Gradient maps are displayed using color gradient rendering, with gradient values ​​transitioning from blue to red from low to high. Fault line symbolization uses different colors to represent fault confidence: high-confidence faults are represented by solid red lines, medium-confidence by dashed orange lines, and low-confidence by dotted yellow lines. Lithological units are filled with geologically accepted color codes: red for granite, gray for limestone, and yellow for sandstone. A topographic base map is overlaid on the target area results map, with target area numbers and matching degree values ​​labeled. Quality control measures include verifying the accuracy of gradient calculations by checking the location of gradient extrema at known geological boundaries. Fault identification results are compared with manually interpreted results; a consistency rate of over 80% is considered acceptable. Lithological classification results are verified against borehole lithological records; a consistency rate exceeding 75% is considered acceptable. Sensitivity analysis is performed on matching degree calculations to analyze the impact of changes in weight parameters on the target area. The final results are organized into thematic atlases and databases for easy use in subsequent exploration work. The data flow is automated, from standardized gravity and magnetic data input to target area results output, generating intermediate results and quality reports at each stage. The processing software is developed using Python, employing the NumPy library for numerical calculations, the GDAL library for raster processing, and the Scikit-learn library for classification algorithms. Process logs record parameter settings and runtime status for each step, facilitating traceability and reproducibility. The entire implementation process adheres to geophysical data processing standards to ensure the reliability and repeatability of the results.

[0038] Example 3: Multi-scale decomposition processing is performed on gravity and magnetic grid data of potential mineralized target areas, with each target area defined as an independent processing unit. The wavelet transform method uses the Daubechies wavelet basis function, with an order of 4, providing tight support and appropriate smoothness. The number of decomposition layers is determined based on the depth of the target, consisting of four layers corresponding to shallow anomalies (0-300 meters), mid-shallow anomalies (300-800 meters), mid-deep anomalies (800-1500 meters), and deep anomalies (above 1500 meters). The wavelet decomposition algorithm uses the Mallat pyramid algorithm, performing two-dimensional discrete wavelet transforms on both gravity and magnetic grid data. The transform process includes row and column transformations, with each direction decomposing into low-frequency and high-frequency components. The first layer decomposition yields detail coefficients D1 and approximation coefficients A1; the second layer decomposes A1 into D2 and A2, and so on up to the fourth layer.

[0039] Wavelet coefficient extraction is performed for each decomposition layer, and the detail coefficients reflect local anomaly features at different scales. Coefficient thresholding uses a soft thresholding function, with the threshold λ calculated based on the noise level.

[0040] in: This represents the threshold of the j-th layer. is the standard deviation of the wavelet coefficients at the j-th level, and N is the total number of grid data points. Coefficient reconstruction is achieved through inverse wavelet transform, and the retained coefficients are reconstructed layer by layer to generate gravity and magnetic anomaly distribution maps corresponding to four scales. The distribution maps are stored in a grid format, with the grid size remaining constant at 50m × 50m.

[0041] 3D geological modeling is based on a multi-scale anomaly distribution map, and geological interface identification employs a moving window gradient extremum detection method. The window size is adjusted according to the anomaly scale, with a 3×3 window used for shallow anomalies and a 7×7 window used for deep anomalies. The interface identification algorithm calculates the gradient magnitude of each grid point and identifies the gradient maxima as geological boundary points. Stratigraphic boundaries are formed by connecting the gradient maxima of gravity anomalies, and rock mass boundaries are determined by the gradient maxima of magnetic anomalies. Spatial interpolation of boundary points generates continuous interfaces using the inverse distance weighted method, with a weighting index set to 2. Borehole data integration includes lithological records and depth data from all exploration boreholes within the survey area, and geophysical logging data includes density logging and magnetic susceptibility logging curves. Constraint processing uses a control point weighting method, assigning higher weights (10) to interface control points at borehole locations and a weight of 1 to areas without boreholes. Interface smoothing employs a minimum curvature algorithm to ensure the geometric rationality of geological interfaces. Density volume was constructed using multi-scale gravity anomaly inversion, employing a constrained least squares method with constraints derived from the interface model. Magnetization volume was constructed using a similar method, based on multi-scale magnetic anomaly inversion calculations. The 3D mesh model was configured with a horizontal grid spacing of 50 meters and a vertical grid spacing of 20 meters. Grid attribute assignments were determined based on the interface location, assigning each grid cell a density value, magnetic susceptibility value, and lithological code. The model output was in 3D volume data format, containing grid coordinate information and physical property values. Model validation was achieved through forward modeling, comparing the model response with the goodness of fit to the measured anomalies.

[0042] Before wavelet transform calculations, the grid data is extended at the boundary using a symmetric extension method to avoid boundary effects. The selection of decomposition layers is verified through spectral analysis to ensure that each scale corresponds to a clear depth range. During coefficient reconstruction, low-frequency components reflect the regional background field, while high-frequency components reflect local anomalies. In thresholding, the noise standard deviation is estimated using the highest frequency wavelet coefficients. The reconstructed anomaly distribution map is plotted using contour lines, with the contour line spacing dynamically adjusted according to the anomaly amplitude. During geological interface identification, gradient extrema are clustered to remove isolated points. The interface trend surface is fitted using a quadratic polynomial function, with the fitting residual controlled within the allowable range. Borehole data preprocessing includes coordinate unification and depth correction to ensure spatial consistency with gravity and magnetic data. Well logging curves are depth-matched and corrected for environmental impact to improve data quality. Density volume inversion sets density constraints: sedimentary rock density 2.4-2.7 g / cm³, igneous rock density 2.7-3.0 g / cm³. Magnetization volume inversion considers remanence and uses a total magnetization vector inversion method. When establishing the model grid system, the coordinate origin was set to the southwest corner of the survey area, and the elevation datum was set to sea level. Trilinear interpolation was used for grid attribute interpolation to ensure numerical continuity. Model visualization employed 3D rendering technology, with different colors representing different lithologies and transparency settings highlighting deep features. Model accuracy was evaluated through cross-validation, reserving 20% ​​of the borehole data as validation points. The wavelet basis function selection underwent comparative testing, with the Daubechies 4th order wavelet achieving a balance between preserving signal characteristics and computational efficiency. Due to the computationally intensive decomposition process, a block-based processing strategy was adopted, dividing the large grid into overlapping sub-regions for separate processing. A sparse matrix format was used for coefficient storage to reduce storage space requirements.

[0043] The reconstructed anomaly map undergoes a grid consistency check to ensure spatial alignment of data at all scales. The interface recognition algorithm is sensitive to noise; therefore, the anomaly map is preprocessed with median filtering in a 3×3 window. Boundary point connections utilize a minimum spanning tree algorithm to generate continuous, smooth boundary lines. In borehole-constrained weighted interpolation, the influence radius of control points is set to 200 meters; beyond this radius, conventional interpolation is used. The density volume inversion iteration count is set to 100, with a convergence tolerance of 0.001. Magnetization volume inversion considers the geomagnetic field direction, undergoing polarization processing before inversion. The 3D model is constructed using a layered structure, building layer by layer from shallow to deep. The thickness of each layer is adjusted according to geological conditions: sedimentary layers are 100-300 meters thick, and bedrock layers are 500-1000 meters thick. The model grid exceeds one million elements and is managed using an octree structure to improve access efficiency. The model output format supports import from common geoscience software, such as GOCAD and Surfer. Quality control includes wavelet decomposition reconstruction error checking; an error norm of less than 5% is considered acceptable. The interface recognition results were compared with the geological map; a boundary position deviation of less than 100 meters was considered acceptable. The 3D model underwent volume calculation checks to ensure the conservation of physical quantities. The final model was reviewed by experts, who corrected any unreasonable geological structures.

[0044] See Figure 4 Multi-scale decomposition techniques were used to demonstrate the gravity and magnetic anomaly characteristics at different depth levels. The top image shows the distribution of shallow anomalies (0-800 meters), while the bottom image shows the characteristics of deep anomalies (above 1500 meters). Shallow anomalies mainly reflect the influence of near-surface geological bodies, including weathered layers, shallow rock masses, and small structures. These anomalies typically exhibit high spatial frequency and strong local variation characteristics. Deep anomalies reveal the distribution pattern of basement structures and large rock masses, showing a broad and gentle regional anomaly pattern. Multi-scale analysis can effectively separate the contributions of anomalies at different source depths. Several local high anomaly centers can be clearly identified in the shallow anomaly map, which may be related to shallow mineralization. Deep anomalies reveal the regional tectonic framework, providing important evidence for understanding the three-dimensional structure of the mineralization system. By comparing the anomaly characteristics at different scales, the spatial distribution and genetic relationships of ore bodies can be better constrained.

[0045] Example 4: Gravity and magnetic data inversion is conducted based on a geological structure model, which is an initial model of density and magnetic susceptibility volumes stored in a three-dimensional grid. A damped least squares inversion algorithm is selected for the constrained inversion algorithm. The objective function consists of a data fitting term and a model constraint term. The data fitting term calculates the sum of squared residuals between the measured gravity data and the model's forward response; the same sum of squared residuals is calculated for the magnetic data. The model constraint term includes smoothing constraints and boundary constraints. Smoothing constraints are achieved by minimizing the gradients of adjacent grid values ​​in the model, while boundary constraints utilize interface information provided by the geological structure model to control the range of property variations. The inversion iteration uses the conjugate gradient method, updating the three-dimensional density and magnetic susceptibility values ​​in each iteration. The iteration stopping condition is set when the residual decrease rate is less than one ten-thousandth or when the maximum number of iterations (200) is reached. During the inversion process, the forward calculation uses a cuboid element stacking method, treating each grid cell as a cuboid with uniform properties, calculating its gravity and magnetic effects. The main diagonal elements of the data fitting weighting matrix are the reciprocals of the standard deviation of the observed data, while the off-diagonal elements are zero. Model constraint weights are determined using the L-curve method to balance data fitting with model smoothness. The inversion results are output as updated 3D density and 3D magnetic susceptibility distributions, stored as 3D arrays. Reliability assessment calculates the posterior standard deviation for each grid point, obtained through resolution matrix analysis. The inversion error distribution map displays the spatial distribution characteristics of the residuals, identifying high-error regions as areas of poor data quality or insufficient model parameters. Constraint adjustments are made specifically for high-error regions, either by increasing smoothing constraint weights or introducing prior information constraints.

[0046] Referring to Table 1, orebody identification was performed based on the inverted density and magnetic susceptibility distributions. The density threshold was set to 2.8 g / cm³, and the magnetic susceptibility threshold was set to 0.05 SI units. Threshold filtering traversed each grid cell, marking cells with a density value greater than 2.8 and a magnetic susceptibility value greater than 0.05 as candidate anomalous cells. Spatial clustering analysis employed the DBSCAN algorithm, with the neighborhood radius eps set to 150 meters and the minimum number of points min_samples set to 8 grid cells. The clustering results divided connected candidate anomalous cells into independent orebody cells, each assigned a unique number. Orebody spatial distribution feature extraction included calculating the centroid 3D coordinates, volume, average density, and average magnetic susceptibility of each orebody cell. Morphological features were obtained by calculating the minimum enclosing ellipsoid of the orebody cell; the major axis, minor axis, and vertical axis of the ellipsoid characterized the orebody's extension direction.

[0047] Table 1: Threshold Table for Ore Body Identification Parameters Parameter type lower threshold upper limit of threshold unit Applicable conditions density value 2.80 3.50 g / cm³ Skarn-type rich iron ore magnetic susceptibility 0.05 1.00 SI Skarn-type rich iron ore volume 1000 1000000 m³ Ore body size screening Aspect Ratio 0.1 10.0 Dimensionless Ore body morphology screening Density distribution data preprocessing includes background trend removal, which is obtained through 3D polynomial fitting with a polynomial order of 3. Residual field calculation is performed on the magnetic susceptibility distribution data; the residual field is the measured field minus the regional field. Threshold settings are based on known ore deposit statistics: a density threshold of 2.8 g / cm³ corresponds to the typical density value of skarn-type iron ore, and a magnetic susceptibility threshold of 0.05 SI corresponds to the enhanced magnetic susceptibility caused by magnetite mineralization in skarn. In spatial clustering analysis, the neighborhood radius of the DBSCAN algorithm is determined based on three times the grid spacing, with a minimum number of points ensuring the ore body has a certain scale. Morphological opening operations are performed on the ore body unit boundaries to eliminate small protrusions and voids. The ore body feature extraction algorithm calculates the geometric center coordinates of each ore body unit, obtained by averaging the unit position coordinates. Volume calculation is performed by multiplying the total number of grid units belonging to the ore body by the unit volume. The minimum bounding ellipsoid is calculated using principal component analysis, with the principal axis of the ellipsoid corresponding to the eigenvector direction. The calculation of attitude elements includes strike, dip angle, and plunging angle. The strike is the azimuth angle of the projection of the major axis of the ellipsoid onto the horizontal plane, the dip angle is the angle between the major axis and the horizontal plane, and the plunging angle is the projection angle of the major axis onto the vertical plane. The spatial distribution characteristics of the ore body are stored in vector format, including the ore body number, vertex coordinates, and attribute table.

[0048] Before inversion calculations, the 3D mesh model is merged, combining small-scale meshes into larger-scale meshes for inversion. The horizontal scale is merged from 50 meters to 100 meters, and the vertical scale from 20 meters to 40 meters. Forward modeling is accelerated using the Fast Fourier Transform algorithm, converting spatial domain convolution into frequency domain multiplication. During iteration, preconditioning techniques are used to improve the condition number of the coefficient matrix, with the Jacobian preconditioner selected. Residual monitoring displays the fitting error curve in real time, automatically adjusting the step size factor when the fitting error decreases slowly. Resolution analysis is performed on the inversion results, calculating the diagonal elements of the model's resolution matrix; areas with element values ​​greater than 0.7 are considered to have reliable resolution. During orebody identification, multi-parameter fusion analysis is performed; only areas where the overlap between density anomalies and magnetic susceptibility anomalies is greater than 70% are identified as mineralized anomalies. False anomaly removal utilizes geological structure model information to eliminate isolated anomalies located within known rock masses or sedimentary layers. Orebody unit classification is based on physical and morphological characteristics, categorized into types such as dense massive orebody, disseminated orebody, and banded orebody. The generated maps include a planar distribution map, a cross-section map, and a three-dimensional map of the ore body. The planar map indicates the ore body number and burial depth, the cross-section map shows the ore body morphology and stratigraphic relationship, and the three-dimensional map shows the spatial distribution of the ore body.

[0049] The quality control process includes comparing the inversion results with the borehole core density measurements; a density error of less than 0.1 g / cm³ is considered acceptable. The magnetic susceptibility inversion results are compared with ground magnetic survey data; a morphological consistency greater than 80% is acceptable. The orebody identification results are overlaid with known orebody boundaries; a boundary position deviation of less than 50 meters is acceptable. The entire process is automated, from 3D model input to orebody feature output, with intermediate results saved in a standard format. The processing log records detailed information such as inversion parameters, iteration count, and fitting error for subsequent analysis and optimization. Gravity and magnetic data inversion calculations are performed on a high-performance computing cluster, using the MPI parallel algorithm to accelerate large-scale 3D inversion. Each computing node processes a sub-region, and data exchange between nodes is achieved through message passing. The orebody identification algorithm optimizes memory usage and employs a block processing strategy for large-scale 3D data. A database system is established for results data management, supporting spatial query and statistical analysis functions. The final output includes an orebody spatial distribution feature file and a quality control report. The feature file contains complete 3D geometric information and physical property parameters, while the report includes a processing flow description and accuracy assessment results.

[0050] Example 5: Spatial Correlation Analysis focuses on the spatial distribution characteristics of ore bodies, which include the three-dimensional coordinates, physical property parameters, and geometric morphology information of each independent ore body unit. The spatial distance between adjacent ore body units is calculated using the Euclidean distance formula, determining the straight-line distance between the centroids of the ore body units. The distance calculation range is set to a search radius of 500 meters; ore body units exceeding this distance are considered to have no spatial correlation. Physical property similarity calculation includes two indicators: density similarity and magnetic susceptibility similarity. Density similarity is calculated using the relative error method, and magnetic susceptibility similarity is calculated using the correlation coefficient method. In the similarity weighting, density similarity has a weight of 0.6, and magnetic susceptibility similarity has a weight of 0.4. A comprehensive scoring model is established for ore body connectivity assessment. The comprehensive score is the weighted sum of the spatial distance score and the physical property similarity score. The distance scoring function uses a linear decay model, decreasing the score by 0.2 for every 100-meter increase in distance. The physical property similarity scoring function uses a piecewise linear function. A similarity greater than 0.8 scores 1, 0.6-0.8 scores 0.8, 0.4-0.6 scores 0.6, and less than 0.4 scores 0. A connection threshold is set to 0.75; ore body unit pairs with a comprehensive score higher than this threshold are marked as connectable unit pairs. Ore body unit merging uses a hierarchical clustering algorithm, starting with each ore body unit as an independent cluster and progressively merging connectable unit pairs. Merging is determined based on the comprehensive score between unit pairs, processed in descending order of score. The merging operation performs 3D bounding box fusion and recalculates the geometric parameters of the merged ore body. After continuous ore bodies are generated, morphological regularization is performed, using a convex hull algorithm to calculate the minimum convex polyhedron boundary. The distribution results of continuous ore bodies are stored as a new set of ore body units, with each continuous ore body assigned a unique identifier.

[0051] The geological constraints integrate known mineralization information and regional metallogenic regularity data. Known mineralization information includes surface mineralization outcrops, borehole mineralization sections, and geochemical anomaly data. Mineralization information is digitized into spatial point layers, and spatial correlation analysis is performed with the orebody distribution. Regional metallogenic regularities are summarized into three categories: ore-controlling structural styles, lithological assemblage types, and mineralization zoning characteristics. These rules are formalized as conditional statements, such as "If the orebody is located at the intersection of a NE-trending fault and a NW-trending fault, and the surrounding rock is a contact zone between limestone and granite, then the mineralization favorability is increased." The spatial distribution characteristics of the orebody are corrected using a rule-based reasoning mechanism, establishing a mapping relationship between geological rules and orebody parameters. Correction includes orebody boundary adjustment, orebody morphology optimization, and false anomaly removal. Boundary adjustment is based on the known spatial distribution of mineralization points; the orebody boundary is expanded in densely mineralized areas and contracted in sparsely mineralized areas. Morphology optimization maintains the geological rationality of the orebody, avoiding sharp protrusions or depressions that do not conform to geological regularities. False anomaly removal targets areas where physical property anomalies contradict geological background, such as high magnetic anomaly areas appearing within known sedimentary strata. Ore body prediction maps are generated using computer-aided mapping technology, and the maps include planar views, cross-sections, and three-dimensional views. Figure 3 The map includes several views. The plan view is projected onto a geographic coordinate system, marking the ore body's boundary lines, isopyres, and elevation information. The cross-sectional view vertically cuts through the ore body, showing its contact relationship with the strata. The 3D view is rendered using virtual reality technology, enabling interactive visualization of the ore body's spatial morphology. Map annotations adhere to mineral exploration map legend standards, and color symbols conform to industry conventions.

[0052] Spatial distance calculation is optimized using a KD-tree data structure to accelerate nearest neighbor search and reduce computational complexity. Data standardization is performed before physical property similarity calculation to eliminate the influence of dimensions. Spatial autocorrelation analysis is introduced into connectivity assessment, and the Moran index is calculated to verify the degree of orebody clustering. The orebody merging algorithm handles special cases, such as when multiple orebody units form a chain connection, using a minimum spanning tree algorithm to determine the merging order. A knowledge base system is established based on geological constraints, comprising a rule base and a fact base. The rule base stores formalized expressions of mineralization laws, while the fact base stores known geological and mineral data. The inference engine employs a forward chain reasoning strategy, applying rules to derive new conclusions from known facts. Uncertainty evaluation is performed on the corrected results, calculating the confidence index for each orebody unit. Orebody prediction maps are output in both vector and raster formats; vector format is used for further editing and analysis, while raster format is used for printing and publication. Map metadata records information such as production time, coordinate system, and data source. The map layout design follows cartographic principles, with complete main map, legend, and scale elements. 3D visualization enables dynamic cutting, supporting the generation of profiles in any direction. The quality control process includes multiple checkpoints. During spatial correlation analysis, the accuracy of distance calculations is checked, and distance values ​​are verified through random sampling. In the orebody connectivity assessment phase, the scoring calculation logic is checked, and the rationality of weight settings is verified. During the geological constraint application phase, the execution of rules is checked to ensure all rules are correctly triggered. In the map generation phase, the completeness of elements is checked to avoid missing important annotations.

[0053] A version control mechanism was established throughout the implementation process to record the changes and basis for each revision. Technical documentation detailed the analysis methods, parameter settings, and processing results. Results data were stored in a standardized format for easy integration with other geoscientific data. The final orebody prediction results comprised spatial data, attribute data, and metadata, forming a complete prediction results system. The spatial analysis algorithm was optimized, employing a block-based computation strategy for large datasets. A traceability mechanism was implemented for geological constraint inference, recording the basis for each inference step. Map generation was automated, requiring no manual intervention from data input to map output. An expert review meeting was organized for results acceptance, and final modifications and improvements were made based on the reviewers' comments.

[0054] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0055] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data, characterized in that, Includes the following steps: Gravity and magnetic measurement data of the target area are collected, and the gravity and magnetic measurement data are preprocessed to obtain a standardized gravity and magnetic dataset. Based on the standardized gravity and magnetic dataset, regional geological background analysis was performed to extract the structural features and lithological distribution features of the target area. Based on the aforementioned structural and lithological distribution characteristics, the mineralization conditions of skarn-type rich iron deposits are matched to screen potential mineralization target areas. The gravity and magnetic data of the potential mineralization target area are decomposed at multiple scales to obtain gravity and magnetic anomaly features at different depth levels. Three-dimensional geological modeling is performed by combining the gravity and magnetic anomaly characteristics at different depth levels to construct a geological structure model of the target area; Based on the geological structure model, gravity and magnetic data inversion is performed to calculate the density distribution and magnetic susceptibility distribution of the target area. Based on the density distribution and magnetic susceptibility distribution, the ore bodies of skarn-type rich iron ore are identified, and the spatial distribution characteristics of the ore bodies are extracted. Spatial correlation analysis was performed on the spatial distribution characteristics of the ore body to determine its continuity and scale. The spatial distribution characteristics of the ore body are optimized based on geological constraints, and the final ore body prediction results are output.

2. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 1, characterized in that, The preprocessing of the gravity and magnetic measurement data to obtain a standardized gravity and magnetic dataset includes the following steps: The gravity measurement data are subjected to terrain correction and regional field correction to eliminate the influence of terrain undulation and regional background field. The magnetic measurement data are subjected to diurnal variation correction and magnetic declination correction to eliminate the influence of geomagnetic diurnal variation and magnetic declination; A gridded interpolation method was used to spatially interpolate the corrected gravity and magnetic measurement data to form regular grid data. The regular grid data is normalized so that the gravity measurement data and magnetic measurement data have the same numerical range, forming the standardized gravity and magnetic dataset.

3. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 2, characterized in that, The regional geological background analysis based on the standardized gravity and magnetic dataset, extracting the tectonic features and lithological distribution features of the target area, includes the following steps: Gradient calculations are performed on the standardized gravity and magnetic dataset to extract gravity gradient and magnetic gradient features; Fracture structure identification is performed based on the gravity gradient and magnetic gradient characteristics to determine the fracture distribution in the target area; The distribution of the faults was verified by combining known geological data, and the main ore-controlling structures were screened. Based on the standardized gravity and magnetic dataset, lithological classification is performed to divide the target region into lithological units; Based on the superposition relationship between the lithological units and the main ore-controlling structures, favorable ore-forming lithological assemblages are extracted.

4. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 3, characterized in that, The process of matching the mineralization conditions of skarn-type rich iron ore deposits based on the structural and lithological distribution characteristics, and screening potential mineralization target areas, includes the following steps: Establish a metallogenic geological model for skarn-type rich iron deposits, and define the ore-controlling structural types and favorable lithological combinations; The main ore-controlling structures and favorable ore-forming lithological combinations are matched with the ore-forming geological model to calculate the ore-forming matching degree; A mineralization matching threshold is set, and areas with a matching degree higher than the threshold are selected as potential mineralization target areas.

5. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 4, characterized in that, The process of performing multi-scale decomposition of the gravity and magnetic data of the potential mineralization target area to obtain gravity and magnetic anomaly features at different depth levels includes the following steps: The gravity and magnetic measurement data of the potential mineralization target area were decomposed into multiple scales using wavelet transform. Wavelet coefficients at different scales were extracted, corresponding to the gravity and magnetic anomaly characteristics of shallow, middle and deep regions, respectively; The wavelet coefficients at different scales are reconstructed to obtain gravity and magnetic anomaly distribution maps at different depth levels.

6. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 5, characterized in that, The process of combining gravity and magnetic anomaly features at different depth levels to perform three-dimensional geological modeling and construct a geological structure model of the target area includes the following steps: Based on the gravity and magnetic anomaly distribution maps at different depth levels, geological interfaces are identified to delineate the boundaries of strata and rock masses. By combining borehole data and geophysical logging data to constrain the geological interface, the modeling accuracy is improved. A three-dimensional interpolation algorithm is used to construct the density volume and magnetic susceptibility volume of the target area, forming the geological structure model.

7. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 6, characterized in that, The process of inverting gravity and magnetic data based on the geological structure model to calculate the density distribution and magnetic susceptibility distribution of the target area includes the following steps: The geological structure model was iteratively optimized using a constrained inversion algorithm to fit measured gravity and magnetic measurement data. Calculate the inverted density distribution and magnetic susceptibility distribution, and evaluate the reliability of the inversion results; Adjust the constraints based on the inversion error distribution to improve inversion accuracy.

8. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 7, characterized in that, The identification of skarn-type rich iron ore bodies based on the density and magnetic susceptibility distributions, and the extraction of spatial distribution characteristics of the ore bodies, includes the following steps: Set density and magnetic susceptibility thresholds to screen out anomalous areas that conform to the characteristics of skarn-type iron-rich minerals; Spatial clustering analysis was performed on the aforementioned abnormal region to divide it into independent ore body units; The spatial location, morphology, and scale characteristics of the independent ore body units are extracted to form the spatial distribution characteristics of the ore body.

9. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 8, characterized in that, The spatial correlation analysis of the spatial distribution characteristics of the ore body to determine its continuity and scale includes the following steps: Calculate the spatial distance and physical property similarity between adjacent ore body units; The connectivity of the ore body is assessed based on the spatial distance and physical property similarity. Merging highly connected ore body units forms a continuous ore body distribution.

10. The method for detecting skarn-type rich iron ore targets in deep-covered areas based on gravity and magnetic data according to claim 9, characterized in that, The optimization of the spatial distribution characteristics of the ore body based on geological constraints, and the output of the final ore body prediction result, includes the following steps: The spatial distribution characteristics of the ore body are corrected based on known mineralization information and geological patterns; Eliminate abnormal areas that do not conform to geological laws and optimize the ore body boundaries; Generate orebody prediction maps and mark the location, depth, and size of the orebody.

Images