Three-dimensional inversion method and system for weak aeromagnetic anomalies of skarn type iron ore

Abstract

The invention discloses a skarn type iron ore weak aeromagnetic anomaly three-dimensional inversion method and a skarn type iron ore weak aeromagnetic anomaly three-dimensional inversion system, which improve the depiction capability of alteration zones and weak magnetic anomalies through a self-adaptive partition management and search mechanism, and realize more flexible backtracking adjustment when the result deviation is large. Meanwhile, hierarchical processing is carried out on different geological prior constraints, so that local precision and global consistency can be considered in a multi-stage or multi-scale environment, and a more reliable prediction basis is provided for the potential position and form of the skarn type ore body.

Description

Technical Field

[0001] The present invention relates to the field of geological exploration, and particularly to a three-dimensional inversion method and system for weak airborne magnetic anomalies of skarn-type iron ore deposits. Background Art

[0002] In the process of mineral resource exploration, using airborne magnetic survey data to locate and evaluate target ore bodies is a relatively mature technical system. Many existing technologies often adopt multi-phase airborne magnetic observations and cooperate with conventional spatial interpolation and coordinate transformation algorithms to generate the magnetic force distribution pattern. After removing or correcting abnormal measurement points through basic filtering and impulse noise reduction means, generally, a geological layer database is used to correct the magnetic field coordinate offset and reduce the observation error, so as to obtain a basic grid model for subsequent analysis. Some solutions will also combine the surface geological sampling results to perform multi-source fusion on the measurement data to achieve early identification and preliminary evaluation of potential ore-forming sections.

[0003] When integrating multi-layer magnetic information, some technical solutions will incorporate the borehole database or regional geological section into the model and identify the strata and structures according to the existing lithological classification. These implementation methods usually introduce geological constraints in the interpolation stage and attempt to correspond the magnetic anomalies to specific lithological zones or alteration types. Subsequently, the interpolation results are displayed through three-dimensional visualization software, or geostatistical methods are used to estimate the magnetization intensity at unknown positions and the possible distribution of ore bodies. Different interpolation techniques, such as covariance analysis or Kriging interpolation, can reflect the spatial correlation between the magnetic field and geological features to a certain extent, thus providing an auxiliary basis for inferring the possible mineralization range.

[0004] In subsequent data processing and anomaly identification, existing methods often use the magnetic anomaly amplitude and gradient distribution to determine the key exploration areas, and then combine geological prior information to comprehensively judge the possible mineralization directions. If significant high-magnetic features are found in some tectonic zones or alteration zones, they are usually regarded as priority exploration targets, and a conventional optimization search process is supplemented to verify the consistency between the magnetic field anomalies and lithological information. Finally, after meeting certain accuracy requirements, a prediction result based on the correspondence between the magnetization intensity and the geological structure is output, providing a preliminary reference for the next drilling deployment and ore body verification.

[0005] Current multi - period observation and three - dimensional interpolation methods based on airborne magnetic survey data often struggle to balance the deviations caused by different observation times, instrument accuracies, and data resolutions during spatial fusion. In the process of identifying weak magnetic anomalies, existing techniques rely more on large - scale overall inferences. In cases where local observations are incomplete or noise values cannot be effectively corrected, potential anomalies are likely to be overlooked or over - smoothed. Facing scattered borehole lithology information and regional geological prior data, these techniques generally lack refined adaptive constraint strategies and are difficult to differentially process different depths or different alteration - sensitive areas, resulting in inflexible or inaccurate predictions of potential ore bodies.

[0006] When integrating into subsequent geological interpretation and evaluation processes, most existing solutions are unable to dynamically correct the contradictions between their prior constraints and measured magnetic anomalies. When there are significant differences between the inversion results of local areas and the data of known boreholes or surface measurement points, the original process usually lacks an effective backtracking mechanism to re - evaluate the model, making it difficult for ore bodies at greater depths or multi - period superimposed alteration zones to be accurately corresponded in existing techniques. Summary of the Invention

[0007] The present invention provides a three - dimensional inversion method for weak airborne magnetic anomalies of skarn - type iron ore, aiming to solve the problem of inflexible or inaccurate prediction of ore bodies in existing technologies.

[0008] A three - dimensional inversion method for weak airborne magnetic anomalies of skarn - type iron ore includes:

[0009] S1: Obtain airborne magnetic data for multiple observation periods, perform data purification processing on the airborne magnetic data, generate a three - dimensional grid using three - dimensional interpolation in convolution form, and construct a prior matrix based on borehole lithology and alteration zone information;

[0010] S2: Identify the contact positions between alteration zones and surrounding rocks based on the three - dimensional grid and prior matrix and divide sub - regions; Apply three - dimensional gradient and difference operators to each sub - region, extract field value mutation points, identify weak magnetic anomalies, record them in the local anomaly sequence, and generate storage results after introducing constraint conditions for each sub - region;

[0011] S3: For the possible ore - forming space within the sub - region, perform global search and population iterative update based on the sparrow search algorithm, and obtain a candidate solution set;

[0012] S4: Use the candidate solution set as the initial condition of variational inference for iteration to generate a three - dimensional inversion result; Define a coupling penalty function between alteration degree and magnetic susceptibility, and adopt a differential learning rate or amplification factor for the sub - region during the iteration process; Quantify the uncertainty by calculating the variance of the posterior distribution, and output the maximum a posteriori solution of the ore body position, shape, and magnetization intensity and the uncertainty evaluation result at the end of the iteration;

[0013] S5: Based on the maximum a posteriori solution, the measured values of known borehole measurement points and close-range aeromagnetic data are compared to calculate the comprehensive errors of position and magnetization intensity; and iterative inversion is re-executed according to the result until the error meets a preset threshold;

[0014] S6: Perform multi-scale fusion on the corrected results of each sub-area and obtain a three-dimensional exploration model;

[0015] S7: extracting the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid from the three-dimensional exploration model, and calculating the potential section based on multiple rounds of search probability and alteration zone distribution characteristics; outputting the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results and priority drilling coordinates.

[0016] Optionally, step S1 includes:

[0017] S11: numbering the aeromagnetic data according to the observation time sequence and forming an original data set D0;

[0018] S12: Import the original data set D0 into the processing environment, determine the difference between data resolution and sampling density; agree that any observation point is at the original coordinate (x i ,y i ,1) is obtained by the following formula to obtain the normalized coordinates (X i ,Y i ):

[0019]

[0020] Where A is a two-dimensional affine transformation matrix, x i With y i represents the horizontal and vertical coordinates of the i-th sampling point in the original data set, X i With Y i Indicates the horizontal and vertical coordinates of the sampling points in the unified coordinate system;

[0021] S13: Obtain the normalized data set D1, perform validity check on the magnetic field strength and coordinate information of each observation point; remove data blocks with missing values exceeding a predetermined limit, and delete abnormal samples in the high-amplitude impulse noise area to obtain a preliminary purified data set D2;

[0022] S14: performing three-dimensional interpolation on the preliminary cleansed data set D2 to obtain a three-dimensional grid G0;

[0023] S15: mapping the borehole lithology and geological logging information into the coordinates of the three-dimensional grid, and mapping the boundary between the rock mass and the skarn alteration zone into the three-dimensional grid;

[0024] S16: Combining the distribution of ore bodies in the confirmed area, forming the prior matrix.

[0025] Optionally, step S2 includes:

[0026] S21: Read the spatial positions where the altered zone contacts the wall rock on the three-dimensional grid and the prior matrix;

[0027] S22: Obtain potential structural parts, and perform zoning management on the three-dimensional grid in combination with the prior matrix, and split it into multiple sub-regions with different degrees of attention;

[0028] S23: Schedule and store the local magnetic field data on the sub-regions, and based on the differential operator and the three-dimensional gradient, characterize the rapid mutation positions of the field values and local weak magnetic anomalies;

[0029] S24: After extracting the weak magnetic anomalies within the sub-regions, introduce the constraint condition set C according to the lithological boundary and altered zone information;

[0030] S25: Complete the search constraints for each sub-region based on the constraint condition set C, and perform parameter strengthening for specific depth intervals within the altered zone;

[0031] S26: After completing the extraction and search constraint configuration of the weak magnetic anomalies in each sub-region, uniformly encapsulate the data and generate the storage result.

[0032] Optionally, S3 further includes:

[0033] S31: Determine the boundary of the sparrow search algorithm in three-dimensional space based on the coordinate information of the weak magnetic anomaly sequence L and the depth constraint z existing in the storage result;

[0034] S32: Construct a permitted search domain D corresponding to each possible boundary, and perform multi-dimensional initial population setting for the key physical quantities in the skarn-type ore body and combine them into individual vectors;

[0035] S33: Apply a geological heuristic factor to the individuals during the population iterative update process;

[0036] S34: Perform a jump operation during the population iterative update process;

[0037] S35: Improve the coverage of multi-peak regions by changing the dynamic vigilance factor with the population diversity; the candidate solution set is distributed in the magnetic anomaly concentration area and the altered constraint area.

[0038] Optionally, S4 further includes:

[0039] S41: In the candidate solution set, based on the sub-region division index, perform a scoring function on all solutions, and select the top several solutions with the highest performance in each sub-region as the initial conditions for variational inference;

[0040] S42: Define a scoring function to comprehensively measure according to the magnetic field fitting degree and geological prior conformity. Screen candidate solutions within the sub-region through the comprehensive scoring function to form the initial solution set Θ0 of variational inference.

[0041] Optionally, the S5 further includes:

[0042] Extract the ore body position coordinates and magnetization intensity from the three-dimensional inversion result; compare with the measured values of known borehole measuring points and near-distance aeromagnetic data, record the measured coordinates of all benchmark measuring points and the magnetization intensity measured at close range, form a comparison sample and calculate the comprehensive error; re-perform iterative inversion according to the result until the error meets the preset threshold, and then output the final three-dimensional distribution prediction of skarn-type iron ore.

[0043] Optionally, the S6 further includes: Merge the updated inversion results of each sub-region into the global three-dimensional model through weighted stitching.

[0044] Optionally, the S7 further includes:

[0045] S71: Extract the magnetic susceptibility distribution of each sub-region of the three-dimensional grid from the three-dimensional exploration model, and perform unified index mapping on it according to the spatial position;

[0046] S72: After obtaining the integrated magnetic susceptibility distribution, synchronously output the ore body shape and uncertainty range; calculate the shape indicator function and uncertainty on the spatial coordinates of each three-dimensional grid based on the posterior distribution accumulated in the previous variational inference;

[0047] S73: In the previous inversion and search iteration, accumulate several rounds of possible solutions for each three-dimensional grid, and obtain the probability evaluation value based on the variational distribution;

[0048] S74: Determine the priority drilling positions, and obtain several potential sections through interval grading;

[0049] S75: Obtain a unified output data set according to S71 to S74, which includes magnetic susceptibility distribution data, ore body shape and uncertainty matrix, priority drilling coordinate point set, and spatial index of the graded potential areas.

[0050] A three-dimensional inversion system for weak aeromagnetic anomalies of skarn-type iron ore in the present invention, the system includes:

[0051] The first acquisition unit is used to acquire aeromagnetic data of multiple observation periods, perform data purification processing on the aeromagnetic data, generate a three-dimensional grid by using three-dimensional interpolation in convolution form, and construct a prior matrix based on borehole lithology and alteration zone information;

[0052] A second acquisition unit, configured to identify the contact positions between the altered zones and the surrounding rocks based on the three-dimensional grid and the prior matrix and divide sub-regions; apply three-dimensional gradients and difference operators to each sub-region, extract the field value mutation points, identify weak magnetic anomalies, record them in the local anomaly sequence, and generate storage results after introducing constraint conditions for each sub-region;

[0053] A third acquisition unit, configured to perform global search and population iterative update on the possible ore-forming spaces within the sub-regions based on the sparrow search algorithm, and obtain a candidate solution set;

[0054] A fourth acquisition unit, configured to use the candidate solution set as the initial condition of variational inference for iteration to generate a three-dimensional inversion result; wherein, a coupling penalty function between the alteration degree and the magnetic susceptibility is defined, and different learning rates or amplification factors are adopted for the sub-regions during the iteration process; quantify the uncertainty by calculating the variance of the posterior distribution, and output the maximum posterior solution of the ore body position, shape and magnetization intensity and the uncertainty evaluation result at the end of the iteration;

[0055] A fifth acquisition unit, configured to compare with the measured values of known borehole measurement points and near-distance airborne magnetic data based on the maximum posterior solution, and calculate the comprehensive error of the position and magnetization intensity; re-perform iterative inversion according to the result until the error meets the preset threshold;

[0056] A sixth acquisition unit, configured to perform multi-scale fusion on the results of each sub-region after correction and obtain a three-dimensional exploration model;

[0057] A seventh acquisition unit, configured to extract the magnetic susceptibility distribution and the ore body shape indication function of each three-dimensional grid from the three-dimensional exploration model, and calculate the potential sections based on the multi-round search probability and the distribution characteristics of the altered zones; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification result and preferred drilling coordinates.

[0058] A computer-readable storage medium in the present invention, wherein the computer-readable storage medium stores one or more programs, and the one or more programs can be executed by one or more processors to implement the method described in any one of the above.

[0059] The method of the present invention improves the characterization ability of altered zones and weak magnetic anomalies through an adaptive partition management and search mechanism, and realizes more flexible backtracking adjustment when the result deviation is large. At the same time, it also performs hierarchical processing of different geological prior constraints, so as to take into account both local accuracy and global consistency in multi-phase or multi-scale environments, thereby providing a more reliable prediction basis for the potential position and shape of skarn-type ore bodies. Description of the Drawings

[0060] Figure 1It is a flowchart of a 3D inversion method for weak airborne magnetic anomalies of skarn-type iron ore in an embodiment of the present invention;

[0061] Figure 2 It is a structural diagram of a 3D inversion system for weak airborne magnetic anomalies of skarn-type iron ore in an embodiment of the present invention. Detailed implementation manners

[0062] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the present invention, rather than limiting the present invention. Additionally, it should be noted that for the sake of description, only parts related to the present invention are shown in the accompanying drawings, rather than all the structures.

[0063] It should be understood that in various embodiments herein, the magnitudes of the serial numbers of the above processes do not mean the order of execution. The order of execution of each process should be determined by its function and internal logic, and should not constitute any limitation to the implementation process of the embodiments herein.

[0064] An embodiment of the present invention provides a 3D inversion method for weak airborne magnetic anomalies of skarn-type iron ore. As Figure 1 shown, the method includes:

[0065] S1: Obtain airborne magnetic data for multiple observation periods, generate a 3D grid by using 3D interpolation in convolution form after data purification processing of the airborne magnetic data, and construct a prior matrix based on borehole lithology and alteration zone information; specifically, record the observation year and instrument parameters of the obtained airborne magnetic data. The data purification processing includes: performing coordinate normalization processing on the airborne magnetic data, projecting each observation coordinate to a unified coordinate system by using affine transformation, filtering and removing the measuring points with too high missing value ratio or significant pulse noise, and performing local neighborhood statistical correction on a small number of outlier points. A measuring point is each data point in the airborne magnetic data.

[0066] S2: Identify the contact positions between alteration zones and surrounding rocks and divide sub-regions based on the 3D grid and the prior matrix; apply 3D gradient and difference operators to each sub-region, extract the field value mutation points and identify weak magnetic anomalies, and record them in the local anomaly sequence, and generate a storage result after introducing constraint conditions for each sub-region; specifically, divide sub-regions according to depth, alteration sensitivity, etc. The storage result includes the weak magnetic anomalies and search constraint configurations of each sub-region.

[0067] S3: For the possible ore-forming space within the sub-region, perform global search and population iterative update based on the sparrow search algorithm, and obtain a candidate solution set. Specifically, during the search process, use the geological heuristic factor to improve the search priority of the alteration-sensitive area, and perform a jump operation when the individual continuously converges; improve the coverage of multi-peak regions by changing the dynamic vigilance factor with the population diversity, and finally obtain the candidate solution set.

[0068] S4: Iterate using the candidate solution set as the initial condition for variational inference to generate a three-dimensional inversion result; wherein, define a coupling penalty function between the alteration degree and magnetic susceptibility, and adopt a differential learning rate or amplification factor for the sub-region during the iteration process; quantify the uncertainty by calculating the variance of the posterior distribution, and output the maximum a posteriori solution and the uncertainty evaluation result of the ore body position, shape, and magnetization intensity at the end of the iteration.

[0069] S5: Based on the maximum a posteriori solution, compare with the measured values of known borehole measurement points and near-range aeromagnetic data, and calculate the comprehensive error of position and magnetization intensity; re-perform iterative inversion according to the result until the error meets the preset threshold. Specifically, in the calculated comprehensive error, if the local error is too large, then judge that the reason is prior constraint conflict or insufficient search coverage; if it is insufficient search, then enhance the jump intensity or warning factor of the sparrow algorithm; if it is a prior conflict, then update the penalty intensity of the alteration degree - magnetic susceptibility coupling term. Subsequently, re-perform iterative inversion until the error meets the preset threshold.

[0070] S6: Perform multi-scale fusion on the results of each corrected sub-region to obtain a three-dimensional exploration model. Specifically, adopt encrypted grid meshing and combine conjugate deviation metrics to correct the inversion parameters in high-gradient or multi-phase superposition areas to ensure reasonable consistency with the original variational results; only perform local updates on areas with high confidence, while restart global search for areas with low confidence to finally obtain a unified and refined three-dimensional exploration model.

[0071] S7: Extract the magnetic susceptibility distribution and ore body shape indication function of each three-dimensional grid from the three-dimensional exploration model, and calculate the potential section based on the multi-round search probability and the distribution characteristics of the alteration zone; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification result, and preferred drilling coordinates. Specifically, after completing the fusion model, extract the magnetic susceptibility distribution and ore body shape indication function of each three-dimensional grid, and calculate the potential section based on the multi-round search probability and the distribution characteristics of the alteration zone; output the three-dimensional model and the corresponding ore body spatial position, uncertainty quantification result, and preferred drilling coordinates, so as to realize the positioning and prediction of concealed bodies of skarn-type iron ore.

[0072] The above embodiments of the present invention improve the ability to depict alteration zones and weak magnetic anomalies through an adaptive zoning management and search mechanism, and achieve more flexible backtracking adjustment when the result deviation is large. At the same time, through the hierarchical processing of different geological prior constraints, both local accuracy and global consistency can be taken into account in multi-phase or multi-scale environments, thus providing a more reliable prediction basis for the potential position and shape of skarn-type ore bodies.

[0073] In a preferred embodiment of the present invention, step S1 includes:

[0074] S11: Number the aeromagnetic data according to the observation time series T to form the original data set D0, so as to match the data source and the initial measurement conditions in the subsequent steps.

[0075] S12: Import the original data set D0 into the processing environment, determine the differences in data resolution and sampling density, and for the problem of inconsistent coordinate systems, construct a two-dimensional affine transformation matrix A for normalization to project the measurement points in different coordinate systems onto the same spatial reference system. It is agreed that any observation point in the original coordinates (x i , y i , 1) obtains the normalized coordinates (X i , Y i ) through the following formula:

[0076]

[0077] where A is the two-dimensional affine transformation matrix, x i and y i represent the horizontal and vertical coordinates of the i-th sampling point in the original data set, and X i and Y i represent the horizontal and vertical coordinates of the sampling point after being in the unified coordinate system; the affine transformation matrix A is solved according to the common reference points and is fitted by the least squares method.

[0078] S13: Obtain the normalized data set D1, and perform validity checks on the magnetic field intensity and coordinate information of each observation point; based on the threshold-based filtering mechanism, eliminate the data blocks with the missing value ratio exceeding the predetermined limit, and delete the abnormal samples in the high-amplitude pulse noise area to obtain the preliminary purified data set D2; for the field values with a small number of jump points, correct them based on the local neighborhood statistical values to prevent individual outliers from causing too much interference in the subsequent grid interpolation stage.

[0079] S14: Perform three-dimensional interpolation on the preliminary purified data set D2 to obtain the three-dimensional grid G0; specifically, denote the grid function obtained after interpolation as Let the spatial region be Ω, and perform a triple integral in the form of convolution of the kernel function K(x, y, z; α, β, γ) and the measured point field value F(α, β, γ):

[0080]

[0081] where (x, y, z) are the coordinate points in the target grid, (α, β, γ) are the original coordinates where the measured points are located, K(x, y, z; α, β, γ) is the kernel function, which is used to measure the spatial correlation degree between the target point and the known measured points, and the measured point field value F(α, β, γ) represents the observed value of this measured point.

[0082] The interpolation is completed to obtain a unified three-dimensional grid G0, and the result matches the spatio-temporal distribution characteristics of multi-phase aeromagnetic data.

[0083] S15: Map the borehole lithology and geological logging information to the coordinates of the three-dimensional grid. Let the set of borehole spatial coordinates be R, and the borehole lithology indicator function b(r) be used to represent the lithology category. Denote b(r)=1 as the lithology within the skarn alteration zone, and b(r)=0 as other rock masses. Use the following boundary function B(r,θ,z) to characterize the boundary area distribution:

[0084]

[0085] where r represents the borehole location, θ represents the azimuth parameter related to the direction of the tectonic zone, and z represents the depth. Through this boundary function, the boundary between the rock mass and the skarn alteration zone is mapped into the three-dimensional grid.

[0086] S16: Combine the ore body distribution in the confirmed area to form the prior matrix M. M is composed of the prior labels m i of each grid cell of the three-dimensional grid. If the i-th grid cell falls into the skarn or can be inferred as the ore body boundary area, then m i is 1, otherwise m i is 0:

[0087]

[0088] where m k represents the prior index information of the k-th grid cell, and n is the total number of grid cells.

[0089] M is in a sparse diagonal form, which is convenient for the constraint algorithm to quickly query whether each grid cell is consistent with the ore body prior.

[0090] S17: Combine the three-dimensional grid G0 obtained in the previous step with the prior matrix M to generate the output multi-phase aeromagnetic three-dimensional grid G1 and the supporting prior markings.

[0091] In a preferred embodiment of the present invention, step S2 includes:

[0092] S21: Read the spatial position where the alteration zone contacts the surrounding rock on the three-dimensional grid G1 and the prior matrix M; assume that a boundary area is defined in the three-dimensional coordinate space for identifying the scope of the tectonic zone, and define the indicator function I B (x,y,z) as 1, indicating that the grid point (x,y,z) is within B, and 0 indicating otherwise. For the grid point (x,y,z), use the following structural part locking function Φ w (x,y,z) to determine whether the grid point is in the structural area prone to ore bodies:

[0093] Φ w (x, y, z) = ∫∫∫ Ω δ(I B (u, v, w))·K(x - u, y - v, z - w)dudvdw

[0094] where (u, v, w) represents the integration coordinates, Ω represents a three - dimensional domain, δ takes the value 1 when I B (u, v, w) = 1 and 0 otherwise, and K represents the kernel correlation function used to measure the spatial correlation between (x, y, z) and (u, v, w).

[0095] After obtaining Φ w (x, y, z), if Φ w (x, y, z) is greater than the threshold, it is determined that this position is affected by the tectonic part and is prone to the distribution of skarn - type ore bodies.

[0096] S22: Obtain potential tectonic parts, manage the partition of the three - dimensional grid in combination with the prior matrix, and split it into multiple sub - regions with different degrees of attention; specifically, record the depth coordinate of the three - dimensional domain Ω as z and the alteration sensitivity identifier as m i , and determine whether the grid cell i is marked as an alteration - sensitive area. Define the depth segmentation set Z = {z0, z1, …, z k}, and within each segmentation, use whether it belongs to the alteration - sensitive area as a discriminant condition to form a sub - region set {Γ1, Γ2, …, Γ m}:

[0097] Γ k = {(x, y, z) | z ∈ [z k-1 , z k ), m i = 1} ∪ {(x, y, z) | z ∈ [z k-1 , z k ), m i = 0}

[0098] where k ranges from 1 to m representing the segmentation index, m i = 1 indicates that the grid cell i has alteration sensitivity, and m i = 0 indicates that it does not have alteration sensitivity.

[0099] Through the above logical partitioning, the overall three - dimensional grid is split into multiple sub - regions with different degrees of attention.

[0100] S23: Schedule and store the local magnetic field data on the sub - regions, and record the magnetic field distribution in this region as F(x, y, z). Characterize the rapid mutation positions of the field values based on the difference operator and the three - dimensional gradient G(x, y, z):

[0101]

[0102] Perform a scan on the amplitude of G(x, y, z). If it is greater than the threshold value, it is determined as a magnetic field gradient mutation point.

[0103] Further, statistically analyze the peak value of F(x, y, z) near the mutation point. If the peak value exceeds the preset threshold, it is determined as a magnetization intensity spike and recorded in the local weak magnetic anomaly sequence L.

[0104] S24: After the extraction of weak magnetic anomalies within the sub-region is completed, according to the lithology boundary and alteration zone information, introduce the constraint condition set C; for each sub-region Γ k , define the search weight function ω k and the magnetization intensity threshold τ k , written as the constraint vector Ψ k =(ω k , τ k ). The vector varies with the sub-region position and lithology type:

[0105]

[0106] Among them, σ k represents the coincidence degree between the lithology boundary and Γ k , ρ k represents the alteration sensitivity coefficient, and α and β represent hyperparameters used to adjust the search intensity and magnetization threshold.

[0107] S25: Based on the constraint condition set C, complete the search constraint for each sub-region. Perform parameter enhancement for a specific depth interval within the alteration zone to improve the detection ability of local anomalies. Denote the range where the alteration zone is located at the depth coordinate z ∈ [z a , z b as R ab . In this interval, linearly amplify the three-dimensional gradient G(x, y, z) and the magnetization intensity F(x, y, z). The amplification operator Λ(γ) performs the following transformation on the input variable γ:

[0108] Λ(γ)=γ×(1 + η·h(z))

[0109] Among them, η represents the amplification coefficient, and h(z) represents an increasing function that rises in segments within the range z ∈ [z a , z_b], and a larger amplification coefficient is used for the area closer to this interval.

[0110] By introducing the amplification operator Λ within this depth interval, weak signals can be highlighted and signal attenuation caused by smoothing filtering can be avoided.

[0111] S26: After the weak magnetic anomaly extraction and search constraint configuration for each sub-region are completed, the data is uniformly encapsulated and the storage result O is generated, which includes L and C as well as all tuning parameters:

[0112] O = {L, C, Ψ1, Ψ2, …, Ψ m , Λ, MetaInfo}

[0113] Among them, L represents the weak magnetic anomaly sequence obtained in S23, C represents the set of constraints related to the lithological boundary, Ψ k is the weight and threshold vector of each sub-region, Λ is the amplification operator, and MetaInfo includes the coordinate index, sub-region range, and grid identifier related to this step.

[0114] In a preferred embodiment of the present invention, S3 further includes:

[0115] S31: Based on the coordinate information and depth constraint z of the weak magnetic anomaly sequence L existing in the storage result, the boundary of the sparrow search algorithm in three-dimensional space is determined; based on the coordinate information (x i , y i , z i ) in L and the depth constraint z ∈ [z min , z max , the boundary of the sparrow search algorithm in three-dimensional space is determined. It is agreed that each possible boundary corresponds to an allowable search domain, denoted as D, which represents the set of all feasible coordinates:

[0116] D = {(x, y, z) | x min ≤ x ≤ x max , y min ≤ y ≤ y max , z min ≤ z ≤ z max}

[0117] Among them, x min , x max , y min , y max , z min , z max are obtained by intersecting the sub-region coordinate boundaries output by S2 with the alteration zone depth interval.

[0118] S32: Construct an allowable search domain D corresponding to each possible boundary, and conduct multi-dimensional initial population setting for the key physical quantities in the skarn-type ore body and combine them into an individual vector; specifically, the multi-dimensional initial population includes the magnetization intensity M and the volume estimate V, and they are combined into an individual vector P j = (x j , y j , z j , M j , Vj )。To ensure the diversity of the initial population in the high-dimensional space, an initial distribution function f0 is introduced:

[0119] P j ~f0(D × [0, M max × [0, V max )

[0120] where j represents the individual index, M max and V max represent the upper limit values of the magnetization intensity and volume estimation, which are limited by embedding the weak magnetic anomalies output by S2 into historical experience. Appropriately expand the coverage range of f0 in combination with geological characteristics to meet the exploration requirements for the ore body scale and magnetic parameter distribution.

[0121] S33: During the population iterative update process, apply the geological heuristic factor Ψ j to increase the search probability of the individual P a in the locally altered sensitive area. Denote the current population as S(t), where t represents the iteration step number, represents the state of the j-th individual at the t-th step. Define the geological heuristic correction function H(P j , z j ) to judge the priority of P j in the altered sensitive area:

[0122]

[0123] where α is the adjustment parameter, and R j represents the overlap degree between the sub-region where P j is located and the alteration information, which is obtained through C recorded at S2. When H(P j , z j ) increases, the geological priority of this individual is improved, so that it can obtain a higher search probability in the subsequent position update.

[0124] S34: Perform a jump operation during the population iterative update process. Specifically, during the iteration process, some individuals lack sufficient field value excitation in the area with sparse weak magnetic anomalies and quickly tend to local convergence, which is not conducive to the global optimum. To overcome this limitation, perform a jump operation J(P j ):

[0125]

[0126] U ~ Uniform(0,1)

[0127] where γ is the jump intensity coefficient, and sign represents the sign function, Denote the magnetic field gradient component extracted at S2, which is used to measure the change of neighborhood field values. If it is determined that the fitness improvement of an individual is insufficient within several consecutive generations, then with probability p j Enforce J(P j ). Through this jump mechanism, avoid the individual being trapped in a region with weak anomalies and deficiencies.

[0128] S35: Improve the coverage of multimodal regions by changing the dynamic vigilance factor with population diversity; the candidate solution set is distributed in the magnetic anomaly concentration area and the alteration constraint area. Specifically, to avoid missing potential solutions in a multimodal function environment, design a dynamic vigilance factor λ(t), which changes with the number of iteration steps and is related to the population diversity index θ(t).

[0129] Let θ(t) be the weighted sum of variances of each dimension in the population at step t, and it is calculated as:

[0130]

[0131] where λ(0) represents the initial vigilance value, μ is the weight coefficient, φ is the threshold parameter, and θ(τ) represents the population diversity at the τ-th iteration.

[0132] If the population gradually converges, resulting in θ(τ) being significantly less than φ, the cumulative integral term is negative, causing λ(t) to gradually increase with the increase of iterations. The larger the vigilance factor, the larger the search radius or the greater the tendency to jump, ensuring that there is still a certain probability of finding other feasible solutions in the multimodal region.

[0133] S36: Perform multiple rounds of iteration, and scan the feasible region D by combining the jump and vigilance mechanisms of S34 and S35. Finally, stop the search at the preset number of steps T and output the candidate solution set where each candidate solution contains coordinates (x * , y * , z * ) and physical property parameters (M * , V * ), and is accompanied by the fitness function value. The solutions in the set are distributed in the magnetic anomaly concentration area and the alteration constraint area, and overall reflect the potential distribution relationship of skarn-type ore bodies.

[0134] In a preferred embodiment of the present invention, the S4 further includes:

[0135] S41: In the candidate solution set , based on the sub-region division index k, execute the scoring function Score(P j ) for all solutions P j , and select the top several solutions with the highest performance in each sub-region as the initial conditions for variational inference; denote the selected solution set as Θ0 = {θ 10 , θ20 ,…, θ n0}。

[0136] S42: Define the scoring function According to the comprehensive measurement of the magnetic field fitting degree and the geological prior compliance degree, where F obs is the observed magnetic field, is the magnetic field after the forward calculation of the candidate solution P j , measures the matching degree in the geological prior matrix, and ω1 and ω2 are weighting coefficients.

[0137] Filter the candidate solutions that are more suitable for the inversion requirements in the sub-region through the comprehensive scoring function to form the initial solution set Θ0 of variational inference.

[0138] S43: When constructing p(θ), to strengthen the description of the coupling relationship between the alteration degree and the magnetic susceptibility, define the correlation function R(θ), and increase the penalty when the magnetic susceptibility and the alteration zone intensity deviate from the geological law.

[0139] Denote the magnetic susceptibility component as m(θ) and the alteration degree component as e(θ). Add a penalty term Epen(θ) when constructing the prior p(θ):

[0140] E pen (θ) = exp(κ{[m(θ) - g(e(θ))] 2})

[0141] where g(e(θ)) is the empirical correspondence function between the alteration degree and the magnetic susceptibility, κ is the penalty intensity constant, and E pen is exponentially amplified as the difference between the two increases.

[0142] Integrate this penalty term into p(θ) and affect p(θ, X), which can suppress the solutions that do not conform to the geological prior during the variational inference update process.

[0143] S44: Associate the sub-region index k divided in S2 with the variational distribution and let the learning rate η k (t) perform differential update on this sub-region at the t-th iteration.

[0144] The update rule is defined as follows:

[0145] η k (t + 1) = η k (t) × [1 + α k GradVal(Γ k , t)]

[0146] where αk is the magnification factor corresponding to the sub-region k, and GradVal(Γ k , t) represents the statistic of the internal magnetic field gradient or loss gradient for Γ k at the t-th step of variational inference.

[0147] If the weak magnetic anomaly is significant and a more sensitive inversion is required, then GradVal(Γ k , t) takes a high value, thereby increasing η k (t + 1).

[0148] S45: When asymptotically converging, q(θ; φ * ) gives the approximate posterior distribution of θ, which is used to describe the position, shape, and magnetization intensity of the ore body.

[0149] To quantify the certainty of each dimension, an uncertainty measurement function U(θ d ) is defined, and the variance in the d-th dimension parameter θ d in is calculated as follows:

[0150]

[0151] where θ d represents the d-th dimension of the parameter, and Eθ[·] represents the expectation operation for the distribution.

[0152] Visualize the result of U(θ d ) in three-dimensional space, and display the uncertainty through isosurface mapping.

[0153] S46: When the variational distribution is updated until the end of the iteration, the final parameter posterior is obtained, and it is combined with the three-dimensional grid coordinates to generate the three-dimensional inversion result G2.

[0154] Visualize the shape distribution and magnetization intensity with the maximum a posteriori solution, and attach the uncertainty range.

[0155] The final output structure is expressed as:

[0156]

[0157] where is the ore body configuration and magnetization characteristics of the maximum a posteriori estimation, represents the quantification of uncertainty in each dimension, and CoordRef represents the coordinate mapping in the three-dimensional grid.

[0158] In a preferred embodiment of the present invention, the S5 further includes:

[0159] S51: Extract the ore body position coordinates (x i, y i , z i ), and the magnetization M i ; Compare with the measured values of known borehole measuring points and near-range aeromagnetic data, and record the measured coordinates of all reference measuring points and the magnetization measured at close range to form a comparison sample and calculate the comprehensive error.

[0160] Define the comprehensive error E total to express the total amount of two types of differences:

[0161]

[0162] where n represents the number of measured points, and λ1 and λ2 are the weight coefficients of the position deviation and the magnetization deviation. E total is used to quickly judge the overall error level.

[0163] S52: If the local error value in E total in several sub-regions is significantly larger than the average value, then collect detailed error distribution information for this sub-region, including the ore body boundary and the alteration zone indication function, the magnetization anomaly peak value, and the lithology revealed by the borehole.

[0164] S53: Construct the regional association metric R(Δ, Π) to distinguish whether this deviation conflicts with the variational prior or results from the sparrow search not fully covering the sub-region.

[0165] Define Δ(x, y) to represent the position error at the coordinate (x, y), and Π(x, y) to represent the prior matching difference corresponding to this position.

[0166] If R is high, it is inferred that the deviation concentration does not match the prior constraint;

[0167] If R is low, it indicates that there is a deficiency in the search range or the convergence process:

[0168] R(Δ, Π) = ∫∫ Ω Δ(x, y)Π(x, y)dxdy

[0169] where Ω represents the projection range of the sub-region on the horizontal plane. If the value of R(Δ, Π) in this sub-region is much larger than the average level, it can be judged that the deviation mainly results from the prior constraint conflict, otherwise it is related to the sparrow search deficiency.

[0170] S54: After determining that it is a search deficiency, record the deviation distribution function in the sub-region Γ as ζ(x, y, z). To improve the focusing degree of the sparrow algorithm on this region, accumulate and correct the vigilance factor Λ.

[0171] Define the following formula to describe the new vigilance factor Λ new :

[0172] Λ new = Λ old + ∭ Γ ζ(x, y, z) dx dy dz

[0173] Among them, Λ old represents the alert value of the previous sparrow search algorithm, and the triple integral on the right reflects the total deviation within the sub-region. After the calculation is completed, return to the sparrow search in step S3, and use Λ new to increase the global jump or the search intensity in the key area, ensuring that the areas with insufficient original search get higher weights in the re-iteration.

[0174] S55: When it is found according to the correlation measure R(Δ, Π) in S53 that the deviation is due to a prior conflict, start the prior constraint correction process, and update the penalty intensity κ of the "alteration degree - magnetic susceptibility" correlation function, as well as other hyperparameters.

[0175] Define a new penalty term

[0176]

[0177] Among them, represents the original prior penalty function, μ represents the adjustment coefficient, and μδ prior represents the prior mismatch measure within the sub-region. Trace back to step S4 and perform the variational inference update process again to ensure that the new prior function can match the geological laws and the observed data.

[0178] S56: After completing the correction of the sparrow search process or the prior constraint, use the updated search mechanism to perform three-dimensional inversion again to obtain a new result G3.

[0179] Subsequently, repeat the method in S51 to compare with the borehole and near-distance measurement data, and calculate the deviation values of the position and magnetization intensity, ensuring that the new result is significantly improved in the sub-regions with serious previous errors.

[0180] Until the deviation meets the preset threshold, then output the final three-dimensional distribution prediction of skarn-type iron ore.

[0181] In a preferred embodiment of the present invention, the S6 further includes: merging the updated inversion results of each sub-region into the global three-dimensional model through weighted stitching.

[0182] In the specific steps, it includes:

[0183] S61: In step S5, secondary backtracking correction is performed on several sub-regions to obtain the updated inversion result G k (x, y, z) of each sub-region.

[0184] To merge into the global three-dimensional model, define a weight function Bk (x, y, z) represents the effective coverage area of sub-region k, and the merged three-dimensional model G4(x, y, z) is obtained using the following formula:

[0185]

[0186] where k traverses all sub-regions from 1 to K, and B k (x, y, z) takes larger values within sub-region k and approaches zero in non-covered areas.

[0187] Through the above weighted stitching, smooth transition at the region junction is ensured, and multi-solution conflicts within the same grid cell are avoided. The merged G4(x, y, z) inherits the inversion results optimized and corrected in different sub-regions.

[0188] S62: To highlight the weak anomaly signals within the large ore body aggregation area, based on the statistical analysis of G4(x, y, z), map the magnetization intensity peak and the middle section of the alteration zone to the coordinate set R main , and improve the spatial discretization accuracy therein. Define the local grid scale function h k (x, y, z) to dynamically adjust the grid meshing granularity:

[0189] h k (x, y, z) = h0 × [1 + αGradVal(F(x, y, z))]

[0190] where h0 is the initial grid size, α is the amplification factor, and GradVal(F(x, y, z)) represents the magnetic field gradient statistical value. The gradient is often large within the main ore body area, thus obtaining a smaller h k (x, y, z) to achieve local fine meshing.

[0191] S63: In the exploration of skarn-type ore bodies, multi-stage mineralization will form a complex pattern by superimposing in space and time, and depth conjugate processing needs to be carried out on the easily confused areas. Let the variational inference posterior distribution output by step S4 be q(θ; φ * ), and define the sub-region that may have multi-stage overlap as Ω m .

[0192] To reduce the conflict between the observed data and the prior distribution, the conjugate deviation metric between the new distribution q'(θ) and the original distribution q(θ; φ * ) is denoted as D KL :

[0193]

[0194] If D KLIf the value is always on the high side, it indicates that the current description of the multi-period superposition area seriously deviates from the original posterior. In this case, gradually correct the parameters of q'(θ) to make it more convergent with the original variational result, so as to achieve the purpose of conjugate refinement.

[0195] S64: Through multiple rounds of optimization in steps S3 and S4, the uncertainty evaluation C(x, y, z) at each grid point (x, y, z) can be obtained.

[0196] Perform small-range local updates in high-confidence regions to avoid wasting computing resources.

[0197] In low-confidence regions, restart the global search to prevent the omission of potential ore bodies. Define the partition update function Γ update (x, y, z):

[0198]

[0199] where δ represents the confidence threshold, Γ local and Γ global represent the local update strategy and the global rescan strategy respectively.

[0200] With the help of this partition update mechanism, it is possible to balance high precision and wide coverage and enhance the exploration of uncertain areas.

[0201] S65: When significant differences are found between the magnetic susceptibility or ore body morphology and the geological interpretation in the local grid, introduce a new prior weight W prior (x, y, z) to amplify the geological constraints in this area and promote the variational inference to converge faster to a reasonable solution.

[0202] Define the updated prior weight as:

[0203]

[0204] where χ(x, y, z) is the current magnetic susceptibility, Δ(·) is the matching deviation with the geological interpretation in this area, and ρ is the adjustment coefficient.

[0205] When Δ(χ) becomes larger, significantly increases, helping to constrain the local model during the variational inference process and strengthening geological rationality.

[0206] S66: After completing grid meshing, conjugate processing, and uncertainty partition search, three-dimensional models at different scales are obtained respectively

[0207] To balance the large-scale overview and small-scale details, define the multi-scale fusion expression in the three-dimensional exploration model:

[0208]

[0209] Among them, represents the output result of the l-th scale, and ω l is the fusion weight of different scale results, and L is the number of multi-scale layers. Strengthen the small-scale information of the concealed ore body and suppress noise in a large range.

[0210] In a preferred embodiment of the present invention, the S7 further includes:

[0211] S71: Extract the magnetic susceptibility distribution of each sub-region of the three-dimensional grid from the three-dimensional exploration model, and perform unified index mapping on it according to the spatial position; specifically, extract the magnetic susceptibility distribution χ final (x, y, z) of each sub-region from G k (x, y, z), and perform unified index mapping on it according to the spatial position.

[0212] To retain the contrast of low, medium, and high magnetization characteristics in the final output, define the merging function M total (x, y, z):

[0213]

[0214] Among them, k from 1 to K represents the sub-region index, and Ω k (x, y, z) is the binary weight function indicating sub-region k, represents the point-by-point multiplication operation.

[0215] S72: After obtaining the integrated magnetic susceptibility distribution, synchronously output the ore body shape and the uncertainty range; based on the posterior distribution q(θ; φ * ) accumulated in the previous variational inference, calculate the shape indicator function I mine (x, y, z) and the uncertainty U(x, y, z) for each three-dimensional grid in the spatial coordinates; define the following three-dimensional matrix formula:

[0216]

[0217] Among them, N represents the number of grid cells, and I mine is used to identify the interface between the ore body and the surrounding rock, and U represents the confidence level at that position.

[0218] Visualize this matrix to provide an intuitive spatial structure and credibility distribution of the ore body.

[0219] S73: In the previous inversion and search iteration, accumulate several possible solutions for each three-dimensional grid, and obtain the probability evaluation value based on the variational distribution; to determine the preferred drilling position in S7, introduce the comprehensive probability P of multiple iterative combinations joint(x, y, z) and obtain a probability evaluation value, which is specifically synthesized by the following formula:

[0220]

[0221] Among them, T represents the total number of iterations, pt(x, y, z) represents the probability of the existence of the ore body at the coordinate (x, y, z) after the t-th iteration, and α is an amplification factor used to emphasize the area consistently pointed to by multiple rounds of searches.

[0222] If P joint (x, y, z) is relatively high, it is more valuable to conduct priority drilling here.

[0223] S74: Determine the priority drilling location, and obtain several potential sections through interval grading; define the potential measure D pot (x, y, z), and fuse the magnetic susceptibility M total , the morphological indicator function I mine and the comprehensive probability P joint into an overall score:

[0224] D pot (x, y, z) = β1[M total (x, y, z)] + β2[I mine (x, y, z)] + β3[ln(P joint (x, y, z))]

[0225] Among them, β1, β2, and β3 are weighting coefficients. According to the size of D pot , the area is divided into high-potential areas, medium-potential areas, and marginal potential areas.

[0226] S75. Obtain the unified output dataset O final according to S71 to S74, which includes magnetic susceptibility distribution data M total , the ore body morphology and uncertainty matrix R, the set of priority drilling coordinate points, and the spatial index of the graded potential areas. It can be generally expressed as:

[0227]

[0228] Among them, refers to the set of priority drilling locations, and D pot [zones] represents the distribution and weights of each potential section.

[0229] The specific embodiment of the present invention also provides a three-dimensional inversion system for weak airborne magnetic anomalies of skarn-type iron ore, as Figure 2 shown, the system includes:

[0230] The first acquisition unit 201 is configured to acquire aeromagnetic data of multiple observation periods, perform data purification processing on the aeromagnetic data, generate a three-dimensional grid by using three-dimensional interpolation in a convolutional form, and construct a prior matrix based on borehole lithology and alteration zone information;

[0231] The second acquisition unit 202 is configured to identify the contact positions between alteration zones and surrounding rocks and divide sub-regions based on the three-dimensional grid and the prior matrix; apply three-dimensional gradients and difference operators to each sub-region, extract field value mutation points, identify weak magnetic anomalies, record them in a local anomaly sequence, and generate a storage result after introducing constraint conditions for each sub-region;

[0232] The third acquisition unit 203 is configured to perform global search and population iterative update on the possible ore-forming space within the sub-region based on the sparrow search algorithm, and finally obtain a candidate solution set;

[0233] The fourth acquisition unit 204 is configured to use the candidate solution set as the initial condition of variational inference to perform iteration to generate a three-dimensional inversion result; define a coupling penalty function between the alteration degree and magnetic susceptibility, and adopt a differential learning rate or amplification factor for the sub-region during the iteration process; quantify the uncertainty by calculating the variance of the posterior distribution, and output the maximum a posteriori solution of the ore body position, shape, and magnetization intensity and the uncertainty evaluation result at the end of the iteration;

[0234] The fifth acquisition unit 205 is configured to compare with the measured values of known borehole measuring points and nearby aeromagnetic data based on the maximum a posteriori solution, and calculate the comprehensive error of the position and magnetization intensity; re-perform iterative inversion according to the result until the error meets a preset threshold;

[0235] The sixth acquisition unit 206 is configured to perform multi-scale fusion on the results of each sub-region after correction and obtain a three-dimensional exploration model;

[0236] The seventh acquisition unit 207 is configured to extract the magnetic susceptibility distribution and ore body shape indication function of each three-dimensional grid from the three-dimensional exploration model, and calculate potential sections based on multi-round search probabilities and alteration zone distribution characteristics; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification result, and preferred drilling coordinates.

[0237] An embodiment of the present invention further provides a computer-readable storage medium, which stores one or more programs, and the one or more programs can be executed by one or more processors to implement the method described in any one of the above.

[0238] The method and system according to the embodiments of the present invention have the following technical effects:

[0239] 1. The embodiments of the present invention perform refined coordinate unification and three-dimensional fusion on multi-phase aeromagnetic data. For the differential data generated in each observation period, by establishing an affine transformation matrix and a threshold filtering mechanism, synchronous processing of coordinate projection and anomaly rejection is achieved, and key measurement points and reliable field values are retained within the same three-dimensional grid. Subsequently, the borehole information and geological logging content are mapped to the same coordinate system, laying a foundation for subsequent regional zoning and prior constraints, and keeping multi-source information closely related during the fusion process.

[0240] 2. The entire three-dimensional domain is divided into several sub-regions. Combining lithological boundaries, alteration sensitivity, and spatial gradient characteristics, weak magnetic anomalies within potential tectonic zones are highlighted. At the same time, according to the depth range and known lithological location within the region, an adjustable search weight and amplification operator are constructed to highlight the field value signals in key depths or sensitive areas. Thus, on the same three-dimensional grid, both strongly concentrated magnetic anomalies can be discovered, and local relatively weak anomaly effects can be fully captured, providing more targeted input data for inversion and optimization.

[0241] 3. In the global search and variational inference inversion stage, the sparrow search algorithm's iterative update and jumping mechanism are used to cover possible ore body areas, and the attention to points near the alteration zone is dynamically enhanced by combining geological heuristic factors. To ensure that the final inversion is consistent with the measured data both locally and globally, the embodiments of the present invention propose a backtracking correction step to perform differential iteration on sub-regions with deviations; when the search coverage is insufficient, the search intensity is enhanced and re-scanned; if it is determined that there is a prior conflict, the coupling penalty function is adjusted to conform to geological laws.

[0242] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, system, or computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memory, CD-ROM, optical memory, etc.) containing computer-usable program code.

[0243] The present application is described with reference to the flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or block in the flowchart and / or block diagram, and the combination of processes and / or blocks in the flowchart and / or block diagram, can be realized by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate for realizing the process Figure 1one or more processes and / or blocks Figure 1 means for the functions specified in one or more blocks

[0244] These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means that implement the functions in the process Figure 1 one or more processes and / or blocks Figure 1 specified in one or more blocks

[0245] These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, whereby the instructions executed on the computer or other programmable apparatus provide steps for implementing the functions in the process Figure 1 one or more processes and / or blocks Figure 1 specified in one or more blocks

[0246] The foregoing description of specific exemplary embodiments of the invention has been presented for purposes of illustration and example. These descriptions are not intended to limit the invention to the precise forms disclosed, and it is evident that many modifications and variations are possible in light of the above teaching. The purpose of selecting and describing the exemplary embodiments was to explain the particular principles of the invention and its practical application so as to enable those skilled in the art to implement and utilize the invention in its various different exemplary embodiments as well as various different selections and modifications. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A three-dimensional inversion method for weak aeromagnetic anomalies of skarn-type iron ore, characterized in that, The method comprises: S1: Acquire aeromagnetic data of multiple observation periods, perform data purification on the aeromagnetic data, generate a three-dimensional grid using three-dimensional interpolation in the form of convolution, and construct a priori matrix based on borehole lithology and alteration zone information; S2: Based on the three-dimensional grid and the prior matrix, the contact position between the alteration zone and the surrounding rock is identified and divided into sub-areas; the three-dimensional gradient and differential operators are applied to each sub-area, the field value mutation points are extracted and the weak magnetic anomalies are identified and recorded in the local anomaly sequence, and the storage results are generated after the constraint conditions are introduced for each sub-area; S3: For the possible mineralization space in the sub-area, global search and group iterative update are performed based on the sparrow search algorithm to obtain a set of candidate solutions; S4: using the candidate solution set as the initial condition of variational reasoning to iterate to generate a three-dimensional inversion result; wherein, a coupling penalty function between the degree of alteration and the magnetic susceptibility is defined, and a differentiated learning rate or amplification factor is used for the sub-region during the iteration process; the uncertainty is quantified by calculating the variance of the posterior distribution, and the maximum posterior solution and uncertainty assessment results of the ore body position, morphology and magnetization intensity are output at the end of the iteration; S5: Based on the maximum a posteriori solution, the measured values of known borehole measurement points and close-range aeromagnetic data are compared to calculate the comprehensive errors of position and magnetization intensity; and iterative inversion is re-executed according to the result until the error meets a preset threshold; S6: Perform multi-scale fusion on the corrected results of each sub-area and obtain a three-dimensional exploration model; S7: extracting the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid from the three-dimensional exploration model, and calculating the potential section based on multiple rounds of search probability and alteration zone distribution characteristics; outputting the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results and priority drilling coordinates.

2. The method according to claim 1, characterized in that, Step S1 includes: S11: numbering the aeromagnetic data according to the observation time sequence and forming an original data set D0; S12: Import the original dataset D0 into the processing environment to determine the difference between the data resolution and the sampling density; it is agreed that for any observation point, the normalized coordinates (X i , Y i ) are obtained through the following formula under the original coordinates (x i , y i , 1): where A is a two-dimensional affine transformation matrix, and x i and y i represent the horizontal and vertical coordinates of the i-th sampling point in the original dataset, and X i and Y i represent the horizontal and vertical coordinates of the sampling point after being in the unified coordinate system; S13: Obtain the normalized data set D1, perform validity check on the magnetic field strength and coordinate information of each observation point; remove data blocks with missing values exceeding a predetermined limit, and delete abnormal samples in the high-amplitude impulse noise area to obtain a preliminary purified data set D2; S14: performing three-dimensional interpolation on the preliminary cleansed data set D2 to obtain a three-dimensional grid G0; S15: mapping the borehole lithology and geological logging information into the coordinates of the three-dimensional grid, and mapping the boundary between the rock mass and the skarn alteration zone into the three-dimensional grid; S16: Combining the distribution of ore bodies in the confirmed area, forming the prior matrix.

3. The method according to claim 1, characterized in that, Step S2 includes: S21: reading the spatial position where the alteration zone contacts the surrounding rock on the three-dimensional grid and the priori matrix; S22: Obtaining potential structural parts, and performing partition management on the three-dimensional grid in combination with the prior matrix, and splitting the three-dimensional grid into a plurality of sub-regions with different attention levels; S23: scheduling and storing local magnetic field data in the sub-area, and characterizing the rapid mutation position of the field value and the local weak magnetic anomaly based on the differential operator and the three-dimensional gradient; S24: After extracting the weak magnetic anomalies within the sub-region, introduce the constraint set C based on the lithological boundary and alteration zone information; S25: Perform search constraints on each sub-region based on the constraint set C, and execute parameter enhancement for a specific depth interval within the alteration zone; S26: After completing the extraction of weak magnetic anomalies and search constraint configuration for each sub-region, uniformly encapsulate the data and generate the storage result.

4. The method according to claim 1, wherein The S3 further includes: S31: Determine the boundaries of the sparrow search algorithm in three-dimensional space based on the coordinate information and depth constraint z of the weak magnetic anomaly sequence L existing in the storage result; S32: Construct a permitted search domain D corresponding to each possible boundary, and conduct multi-dimensional initial population setting for the key physical quantities in the skarn-type ore body and combine them into individual vectors; S33: Apply a geological heuristic factor to individuals during the population iterative update process; S34: Perform a jump operation during the population iterative update process; S35: Improve the coverage of the multi-peak region by changing the dynamic warning factor with the population diversity; the candidate solution set is distributed in the magnetic anomaly concentration area and the alteration constraint area.

5. The method according to claim 1, wherein The S4 further includes: S41: In the candidate solution set, execute a scoring function for all solutions based on the sub-region division index, and select the top several solutions with the highest performance in each sub-region as the initial conditions for variational inference; S42: Define that the scoring function is comprehensively measured according to the magnetic field fitting degree and geological prior conformity, and screen candidate solutions within the sub-region through the comprehensive scoring function to form the initial solution set Θ0 of variational inference.

6. The method according to claim 1, wherein The S5 further includes: Extract the ore body position coordinates and magnetization intensity from the three-dimensional inversion result; compare with the measured values of known borehole measuring points and near-range aeromagnetic data, record the measured coordinates of all reference measuring points and the near-range measured magnetization intensity, form a comparison sample and calculate the comprehensive error; re-execute the iterative inversion according to the result until the error meets the preset threshold, and then output the final three-dimensional distribution prediction of the skarn-type iron ore.

7. The method according to claim 1, characterized in that The S6 further includes: Merge the updated inversion results of each sub-region into the global three-dimensional model through weighted stitching.

8. The method according to claim 1, characterized in that The S7 further includes: S71: Extract the magnetic susceptibility distribution of each sub-region of the three-dimensional grid from the three-dimensional exploration model, and perform unified index mapping on it according to the spatial position; S72: After obtaining the integrated magnetic susceptibility distribution, synchronously output the ore body morphology and uncertainty range; calculate the morphological indicator function and uncertainty for each three-dimensional grid in the spatial coordinates based on the posterior distribution accumulated in the previous variational inference; S73: In the previous inversion and search iteration, accumulate several rounds of possible solutions for each three-dimensional grid, and obtain the probability evaluation value based on the variational distribution; S74: Determine the priority drilling positions, and obtain several potential sections through interval grading; S75: Obtain the unified output data set according to S71 to S74, which includes magnetic susceptibility distribution data, ore body morphology and uncertainty matrix, priority drilling coordinate point set, and spatial index of the graded potential areas.

9. A 3D inversion system for weak airborne magnetic anomalies of skarn-type iron deposits, characterized in that, The system includes: The first acquisition unit is configured to acquire aeromagnetic data for multiple observation periods, perform data purification processing on the aeromagnetic data, generate a three-dimensional grid by using three-dimensional interpolation in a convolutional form, and construct a prior matrix based on borehole lithology and alteration zone information; The second acquisition unit is configured to identify the contact positions between alteration zones and surrounding rocks and divide sub-regions based on the three-dimensional grid and the prior matrix; apply three-dimensional gradient and difference operators to each sub-region, extract field value mutation points, identify weak magnetic anomalies, record them in a local anomaly sequence, and generate a storage result after introducing constraint conditions for each sub-region; The third acquisition unit is configured to perform global search and population iterative update based on the sparrow search algorithm for the possible ore-forming space within the sub-region, and obtain a candidate solution set; The fourth acquisition unit is configured to use the candidate solution set as the initial condition of variational inference to perform iteration to generate a three-dimensional inversion result; define a coupling penalty function between alteration degree and magnetic susceptibility, and adopt a differential learning rate or amplification factor for the sub-region during the iteration process; quantify the uncertainty by calculating the variance of the posterior distribution, and output the maximum a posteriori solution of the ore body position, shape and magnetization intensity and the uncertainty evaluation result at the end of the iteration; The fifth acquisition unit is configured to compare the maximum a posteriori solution with the measured values of known borehole measurement points and nearby aeromagnetic data, and calculate the comprehensive error of position and magnetization intensity; re-perform iterative inversion according to the result until the error meets a preset threshold; The sixth acquisition unit is configured to perform multi-scale fusion on the results of each sub-region after correction to obtain a three-dimensional exploration model; The seventh acquisition unit is configured to extract the magnetic susceptibility distribution and ore body shape indication function of each three-dimensional grid from the three-dimensional exploration model, and calculate potential sections based on multi-round search probabilities and alteration zone distribution characteristics; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification result and priority drilling coordinates.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores one or more programs, and the one or more programs can be executed by one or more processors to implement the method according to any one of claims 1 to 8.

Images