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

Abstract

This invention discloses a three-dimensional inversion method and system for weak aeromagnetic anomalies in skarn-type iron ore. Through adaptive zoning management and search mechanisms, the system improves the characterization of alteration zones and weak magnetic anomalies, and enables more flexible retrospective adjustments when the results deviate significantly. Furthermore, the system incorporates a layered approach to different geological prior constraints, balancing local accuracy and global consistency across multiple periods and scales, thus providing a more reliable basis for predicting the potential location and morphology of skarn-type ore bodies.

Description

Technical Field

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

[0002] In the process of mineral resource exploration, the use of aeromagnetic data to locate and evaluate target ore bodies is a relatively mature technology. Many existing technologies often use multiple periods of aeromagnetic observations, combined with conventional spatial interpolation and coordinate transformation algorithms to generate magnetic distribution patterns. After eliminating or correcting anomalous measurement points through basic filtering and pulse denoising, a geological layer database is generally used to correct magnetic field coordinate offsets and reduce observation errors, thereby obtaining a basic grid model for subsequent analysis. Some solutions also combine surface geological sampling results to fuse multi-source measurement data to achieve early identification and preliminary assessment of potential mineralization zones.

[0003] When integrating multiple layers of magnetic information, some approaches incorporate borehole databases or regional geological profiles into the model and identify strata and structures based on established lithologic classifications. These approaches typically introduce geological constraints during the interpolation phase, attempting to assign magnetic anomalies to specific lithologic zones or alteration types. The interpolation results are then displayed using 3D visualization software, or geostatistical methods are used to estimate the magnetization intensity and likely distribution of ore bodies at unknown locations. Different interpolation techniques, such as covariance analysis or kriging interpolation, can, to some extent, reveal the spatial correlation between the magnetic field and geological features, thus providing additional support for inferring the likely extent of mineralization.

[0004] In subsequent data processing and anomaly identification, existing methods often use the magnitude and gradient distribution of magnetic anomalies to identify key exploration areas, then combine this with prior geological information to comprehensively determine the possible orientation of mineralization. If significant high-magnetic features are discovered within certain tectonic or alteration zones, these are typically considered priority exploration targets, supplemented by conventional optimization search processes to verify the consistency between magnetic anomalies and lithologic information. Ultimately, after meeting certain accuracy requirements, predictions based on the correspondence between magnetization intensity and geological structure are output, providing preliminary reference for further drilling deployment and orebody verification.

[0005] Current multi-period observation and three-dimensional interpolation methods based on aeromagnetic data often struggle to account for the deviations caused by different observation times, instrument accuracy, and data resolution during spatial fusion. In the identification of weak magnetic anomalies, existing technologies rely more on large-scale overall inferences. If local observations are incomplete or noise values ​​cannot be effectively corrected, potential anomalies are easily overlooked or over-smoothed. Faced with scattered borehole lithologic information and regional geological prior data, these technologies generally lack sophisticated adaptive constraint strategies, making it difficult to differentiate between different depths or different alteration-sensitive areas, resulting in inflexible or inaccurate predictions of possible ore bodies.

[0006] When integrated into subsequent geological interpretation and evaluation, most existing solutions are unable to dynamically correct for conflicts between their a priori constraints and measured magnetic anomalies. When inversion results for a local area significantly differ from known drillhole or surface measurement data, existing processes often lack an effective backtracking mechanism to reassess the model, making it difficult to accurately map deep ore bodies or multi-stage alteration zones using existing technologies. Summary of the Invention

[0007] The present invention provides a three-dimensional inversion method for weak aeromagnetic anomalies of skarn-type iron ore, which is used to solve the problem that the prediction of ore bodies in the prior art is not flexible or accurate enough.

[0008] A three-dimensional inversion method for weak aeromagnetic anomalies of skarn-type iron ore, comprising:

[0009] S1: Acquire aeromagnetic data from 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 the borehole lithology and alteration zone information;

[0010] 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-regions; three-dimensional gradient and differential operators are applied to each sub-region to extract field value mutation points and identify weak magnetic anomalies, which are then recorded in a local anomaly sequence. Constraints are introduced for each sub-region to generate and store results;

[0011] S3: For the possible mineralization space in the sub-area, a global search and group iterative update are performed based on the sparrow search algorithm to obtain a set of candidate solutions;

[0012] S4: Iterating the candidate solution set as the initial conditions of variational inference 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 applied to the sub-region during the iteration process; uncertainty is quantified by calculating the variance of the posterior distribution, and at the end of the iteration, a maximum posterior solution and uncertainty assessment results of the ore body position, morphology and magnetization intensity are output;

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

[0014] S6: Perform multi-scale fusion on the corrected sub-region results 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 Indicates 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 point after the unified coordinate system;

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

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

[0023] S15: mapping the drill hole lithology and geological log 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: reading the spatial position of the contact between the alteration zone and the surrounding rock on the three-dimensional grid and the prior matrix;

[0027] 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 multiple sub-regions with different attention levels;

[0028] S23: scheduling and storing local magnetic field data in the sub-region, 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;

[0029] S24: After completing the extraction of weak magnetic anomalies in the sub-area, a constraint condition set C is introduced based on the lithologic boundary and alteration zone information;

[0030] S25: completing search constraints for each sub-area based on the constraint condition set C, and performing parameter enhancement for a specific depth interval within the alteration zone;

[0031] S26: After completing the extraction of weak magnetic anomalies and search constraint configuration for each sub-area, unify the data packaging and generate storage results.

[0032] Optionally, the S3 further includes:

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

[0034] S32: Construct an allowable search domain D corresponding to each possible boundary, conduct multi-dimensional initial group settings for key physical quantities in skarn ore bodies and combine them into individual vectors;

[0035] S33: Applying geological heuristic factors to individuals during the iterative update process of the group;

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

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

[0038] Optionally, the S4 further includes:

[0039] S41: In the candidate solution set, a scoring function is performed on all solutions based on the subregion partition index, and the top several solutions with the highest performance in each subregion are selected as initial conditions for variational inference;

[0040] S42: Define a scoring function based on a comprehensive measure of magnetic field fitting and geological prior conformity, and use the comprehensive scoring function to screen candidate solutions in the sub-region to form an initial solution set Θ0 for variational reasoning.

[0041] Optionally, the S5 further includes:

[0042] The ore body position coordinates and magnetization intensity are extracted from the three-dimensional inversion results; compared with the measurement values ​​of known drill hole measurement points and close-range aeromagnetic data, the measured coordinates of all benchmark measurement points and the close-range measured magnetization intensity are recorded to form a comparison sample and the comprehensive error is calculated; iterative inversion is re-executed based on the results until the error meets the preset threshold, and the final three-dimensional distribution prediction of skarn-type iron ore is output.

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

[0044] Optionally, the S7 further includes:

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

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

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

[0048] S74: Determine the priority drilling locations and obtain several potential sections through interval classification;

[0049] S75. According to S71 to S74, a unified output data set is obtained, which includes magnetic susceptibility distribution data, ore body morphology and uncertainty matrix, priority drilling coordinate point set, and spatial index of classified potential areas.

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

[0051] A first acquisition unit is configured to acquire aeromagnetic data for 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;

[0052] The second acquisition unit is used to identify the contact position between the alteration zone and the surrounding rock based on the three-dimensional grid and the prior matrix and divide the sub-regions; apply the three-dimensional gradient and differential operators to each sub-region, extract the field value mutation points and identify the weak magnetic anomalies, and then record them in the local anomaly sequence; and generate and store the results after introducing constraints for each sub-region;

[0053] The third acquisition unit is used to perform global search and group iterative update based on the sparrow search algorithm for the possible mineralization space in the sub-area, and obtain a set of candidate solutions;

[0054] a fourth acquisition unit, configured to iterate the candidate solution set as an initial condition for variational inference 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 applied to the sub-region during the iteration process; uncertainty is quantified by calculating the variance of the posterior distribution, and at the end of the iteration, a maximum a posteriori solution for the ore body position, morphology, and magnetization intensity and an uncertainty assessment result are output;

[0055] a fifth acquisition unit, configured to compare the maximum a posteriori solution with known borehole measurement values ​​and short-range aeromagnetic data to calculate a comprehensive error in position and magnetization intensity; and re-execute iterative inversion based on the result until the error meets a preset threshold;

[0056] The sixth acquisition unit is used to perform multi-scale fusion on the corrected sub-region results and obtain a three-dimensional exploration model;

[0057] The seventh acquisition unit is used to extract the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid from the three-dimensional exploration model, and calculate the potential section based on multiple rounds of search probability and alteration zone distribution characteristics; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results and priority drilling coordinates.

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

[0059] The method of this invention improves the characterization of alteration zones and weak magnetic anomalies through adaptive zoning management and search mechanisms, and enables more flexible retrospective adjustments when results deviate significantly. Furthermore, the method also incorporates a layered approach to different geological prior constraints, balancing local accuracy with global consistency across multiple periods and scales, thus providing a more reliable basis for predicting the potential location and morphology of skarn-type ore bodies. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1This is a flow chart of a three-dimensional inversion method for weak aeromagnetic anomalies of skarn-type iron ore in an embodiment of the present invention;

[0061] Figure 2 This is a structural diagram of a three-dimensional inversion system for weak aeromagnetic anomalies of skarn-type iron ore in an embodiment of the present invention. DETAILED DESCRIPTION

[0062] The present invention will be further described in detail below with reference to the accompanying drawings and examples. It will be understood that the specific embodiments described herein are intended only to illustrate the present invention and are not intended to limit the present invention. It should also be noted that, for ease of description, the accompanying drawings only illustrate portions relevant to the present invention, not all structures.

[0063] It should be understood that in the various embodiments of this document, the size of the serial numbers of the above-mentioned processes does not mean the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this document.

[0064] The embodiment of the present invention provides a three-dimensional inversion method for weak aeromagnetic anomalies of skarn iron ore, such as Figure 1 As shown, the method includes:

[0065] S1: Acquire aeromagnetic data from multiple observation periods, perform data purification on the data, and then generate a 3D grid using a convolutional 3D interpolation method. A priori matrix is ​​constructed based on the borehole lithology and alteration zone information. Specifically, the observation year and instrument parameters are recorded for each acquired aeromagnetic data. Data purification includes normalizing the aeromagnetic data coordinates, projecting each observation coordinate onto a unified coordinate system using an affine transformation, filtering and removing measurement points with a high proportion of missing values ​​or significant impulse noise, and performing local neighborhood statistics correction on a small number of outliers. A measurement point is defined as each data point in the aeromagnetic data.

[0066] S2: Based on the three-dimensional grid and prior matrix, the contact location between the alteration zone and the surrounding rock is identified and divided into sub-regions. Three-dimensional gradient and differential operators are applied to each sub-region to extract field value mutation points and identify weak magnetic anomalies, which are then recorded in a local anomaly sequence. Constraints are introduced for each sub-region and stored results are generated. Specifically, the sub-regions are divided according to depth, alteration sensitivity, etc. The stored results include the weak magnetic anomalies and search constraint configuration for each sub-region.

[0067] S3: Based on the sparrow search algorithm, a global search and swarm iterative update are performed on the possible mineralization space within the sub-region to obtain a set of candidate solutions. Specifically, during the search process, geological heuristic factors are used to increase the search priority of alteration-sensitive areas, and jump operations are performed when individuals continue to converge. Dynamic alert factors are changed according to swarm diversity to improve the coverage of multi-peak regions, and finally a set of candidate solutions is obtained.

[0068] S4: The candidate solution set is used as the initial condition of variational inference to iterate and 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.

[0069] S5: Based on the maximum a posteriori solution, the combined errors of position and magnetization intensity are compared with known borehole point measurements and close-range aeromagnetic data to calculate them. Iterative inversion is then re-executed based on the results until the error meets a preset threshold. Specifically, if the calculated combined error is excessively large locally, it is determined to be due to a priori constraint conflict or insufficient search coverage. If the search is insufficient, the jump strength or warning factor of the Sparrow algorithm is increased. If the priori conflict is present, the penalty strength of the alteration degree-magnetic susceptibility coupling term is updated. Subsequently, the iterative inversion is re-executed until the error meets the preset threshold.

[0070] S6: Multi-scale fusion of the corrected sub-region results is performed to obtain a 3D exploration model. Specifically, in areas with high gradients or multi-phase overlap, a refined grid is used and the conjugate deviation metric is combined to correct the inversion parameters to ensure reasonable consistency with the original variational results. Only local updates are performed in areas with high confidence, while a global search is restarted in areas with low confidence, ultimately resulting in a unified and refined 3D exploration model.

[0071] 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 potential segments based on multi-round search probabilities and alteration zone distribution characteristics; outputting the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results, and preferred drilling coordinates. Specifically, after completing the fusion model, extracting the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid, and calculating potential segments based on multi-round search probabilities and alteration zone distribution characteristics; outputting the three-dimensional model and the corresponding ore body spatial position, uncertainty quantification results, and preferred drilling coordinates, thereby achieving the location and prediction of skarn-type iron ore concealed bodies.

[0072] The above-described embodiment of the present invention utilizes adaptive zoning management and search mechanisms to enhance the characterization of alteration zones and weak magnetic anomalies, enabling more flexible retrospective adjustments when results deviate significantly. Furthermore, by layering different geological prior constraints, it balances local accuracy with global consistency across multiple periods and scales, providing a more reliable basis for predicting the potential location and morphology of skarn-type ore bodies.

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

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

[0075] S12: Import the original data set D0 into the processing environment, determine the difference in data resolution and sampling density, and construct a two-dimensional affine transformation matrix A for normalization to project the measurement points in different coordinate systems into the same spatial reference system. i ,y i ,1) is obtained by the following formula to obtain the normalized coordinates (X i ,Y i ):

[0076]

[0077] Where A is a two-dimensional affine transformation matrix, x i with y i Indicates the horizontal and vertical coordinates of the i-th sampling point in the original data set, X i With Y i Represents the horizontal and vertical coordinates of the sampling point after the unified coordinate system; the affine transformation matrix A is solved based on the common reference point and fitted by the least squares method.

[0078] S13: Obtain the normalized data set D1, and perform validity checks on the magnetic field strength and coordinate information of each observation point; use a threshold-based filtering mechanism to remove data blocks where the proportion of missing values ​​exceeds a predetermined limit, and delete abnormal samples in the high-amplitude impulse noise area to obtain a preliminary purified data set D2; for field values ​​with a small number of jump points, correct them based on local neighborhood statistics to prevent individual outliers from causing excessive interference to the subsequent grid interpolation stage.

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

[0080]

[0081] Where (x, y, z) is the coordinate point in the target grid, (α, β, γ) is the original coordinate of the measuring point, K(x, y, z; α, β, γ) is the kernel function used to measure the degree of spatial correlation between the target point and the known measuring points, and the measuring point field value F(α, β, γ) represents the observation value of the measuring point.

[0082] The interpolation is completed to obtain a unified three-dimensional grid G0, and the results match the spatiotemporal distribution characteristics of multi-period aeromagnetic data.

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

[0084]

[0085] Where r represents the drillhole location, θ represents an azimuth parameter related to the tectonic zone direction, and z represents depth. This boundary function is used to map the boundary between the rock mass and the skarn alteration zone onto a three-dimensional grid.

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

[0087]

[0088] Among them, m k Represents the prior indicator information of the kth grid unit, and n is the total number of grid units.

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

[0090] S17: Merge the three-dimensional grid G0 obtained in the previous step with the prior matrix M to generate the output multi-period aeromagnetic three-dimensional grid G1 and the corresponding prior markers.

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

[0092] S21: Read the spatial position of the contact between the alteration zone and the surrounding rock on the three-dimensional grid G1 and the prior matrix M; define the boundary area in the three-dimensional coordinate space Used to identify the construction band range and define the indicator function I B (x, y, z) is 1, indicating that the grid point (x, y, z) is located in B, and 0 indicates other cases. For the grid point (x, y, z), the following part locking function Φ is constructed w (x,y,z) determines whether the grid point is in a structural area where ore bodies are likely to appear:

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

[0094] Among them, (u, v, w) represents the integrand coordinates, Ω represents the three-dimensional domain, and δ represents the time when I B The value is 1 when (u, v, w) = 1, otherwise it is 0. K represents the kernel correlation function, which is used to measure the spatial correlation between (x, y, z) and (u, v, w).

[0095] Get Φ w (x,y,z), if Φ w If (x, y, z) is greater than the threshold, it is judged that the location is affected by the structural position and is prone to skarn-type ore body distribution.

[0096] S22: Obtain potential structural parts, and perform partition management on the three-dimensional grid in combination with the prior matrix, splitting it into multiple sub-regions with different attention levels; specifically, the depth coordinate of the three-dimensional domain Ω is marked as z, and the alteration sensitivity identifier is marked as m i , determine whether grid cell i is marked as an alteration sensitive area. Define the depth segment set Z = {z0,z1,…,z k}, whether it belongs to the alteration sensitive area is used as the judgment condition within each segment to form a sub-area 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] Among them, k continues from 1 to m to represent the segment index, m i =1 means grid cell i has alteration sensitivity, m i =0 means no alteration sensitivity.

[0099] Through the above logical partitioning, the overall three-dimensional grid is divided into multiple sub-areas with different attention levels.

[0100] S23: Schedule and store the local magnetic field data in the sub-region, and record the magnetic field distribution in the region as F(x, y, z). Based on the differential operator and the three-dimensional gradient G(x, y, z), the rapid mutation position of the field value is characterized:

[0101]

[0102] The amplitude of G(x,y,z) is scanned. If it is greater than the threshold, it is determined to be a magnetic field gradient mutation point.

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

[0104] S24: After completing the weak magnetic anomaly extraction in the sub-region, based on the lithologic boundary and alteration zone information, a constraint condition set C is introduced; for each sub-region Γ k , define the search weight function ω k and the magnetization threshold τ k , written as the constraint vector Ψ k =(ω k ,τ k ). The vector varies with the sub-region location and lithologic type:

[0105]

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

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

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

[0109] Among them, η represents the amplification coefficient, h(z) represents the increasing function, and in z∈[z a ,z_b] range, and a larger amplification factor is applied to areas closer to the 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 completing the weak magnetic anomaly extraction and search constraint configuration for each sub-area, the data is packaged uniformly and the storage result O is generated, including L and C and all tuning parameters:

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

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

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

[0115] S31: Determine the boundary of the sparrow search algorithm in the three-dimensional space based on the coordinate information of the weak magnetic anomaly sequence L in the stored result and the depth constraint z; i ,y i ,z i ) and the depth constraint z∈[z min ,z max ], determine the boundaries of the sparrow search algorithm in three-dimensional space. It is agreed that each possible boundary corresponds to an allowed 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 It is obtained by intersecting the sub-region coordinate boundary output by S2 with the depth interval of the alteration zone.

[0118] S32: Construct an allowable search domain D corresponding to each possible boundary, carry out multi-dimensional initial group setting for key physical quantities in skarn ore bodies and combine them into individual vectors; specifically, the multi-dimensional initial group includes magnetization intensity M and volume estimation V, which are combined into individual vectors P j =(x j ,y j ,z j ,M j ,Vj ). In order to ensure the diversity of the initial population in the high-dimensional space, the initial distribution function f0 is introduced:

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

[0120] Among them, j represents the individual index, M max With V max This represents the upper limit of magnetization intensity and volume estimation, which is limited by embedding historical experience with weak magnetic anomalies output by S2. The coverage of f0 is appropriately expanded based on geological characteristics to meet the exploration needs of ore body size and magnetic parameter distribution.

[0121] S33: During the group iterative update process, the individual P j Applying geological heuristic factor Ψ a , increasing its search probability in the local alteration sensitive area. The current population is recorded as S(t), where t represents the number of iteration steps. Denotes the state of the jth individual at step t. Define the geological heuristic correction function H(P j ,z j ) is used to judge P j Priority within alteration sensitive areas:

[0122]

[0123] Among them, α is the adjustment parameter, R j Indicates P j The degree of overlap between the sub-region and the alteration information is obtained by C recorded at S2. j ,z j ) increases, the geological priority of the individual increases, resulting in a higher search probability in subsequent position updates.

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

[0125]

[0126] U~Uniform(0,1)

[0127] Where γ is the jump strength coefficient, sign represents the sign function, It represents the magnetic field gradient component extracted at S2, which is used to measure the change of the neighborhood field value. If it is judged that the fitness of an individual is not improved enough in several consecutive generations, the probability p is used to calculate the value of the magnetic field gradient component extracted at S2. j Enforcement of J(P j ). Through this jumping mechanism, individuals are prevented from falling into areas of weak abnormality and deficiency.

[0128] S35: Improve coverage of multimodal regions by varying a dynamic alert factor with population diversity; the candidate solution set is distributed in areas of concentrated magnetic anomalies and alteration-constrained regions. Specifically, to avoid missing potential solutions in a multimodal function environment, a dynamic alert factor λ(t) is designed that varies with the number of iterations and is related to the population diversity index θ(t).

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

[0130]

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

[0132] If the population gradually converges, causing θ(τ) to be significantly smaller than φ, the cumulative integral term becomes negative, causing λ(t) to gradually increase with each iteration. A larger alert factor increases the search radius or the tendency to jump, ensuring that there is still a certain probability of finding other feasible solutions in multi-peak regions.

[0133] S36: Execute multiple iterations, combine the jump and warning mechanisms of S34 and S35 to scan the feasible domain D. Finally, stop the search at the preset number of steps T and output the candidate solution set Each candidate solution Contains coordinates (x * ,y * ,z * ) and physical parameters (M * ,V * ), along with the fitness function value. The solutions in the ensemble are distributed in the magnetic anomaly concentration area and the alteration constraint area, reflecting the potential distribution relationship of skarn-type ore bodies as a whole.

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

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

[0136] S42: Define scoring function According to the comprehensive measurement of magnetic field fitting and geological prior conformity, Among them, F obs To observe the magnetic field, is the candidate solution P j After forward calculation, the magnetic field measure The matching degree in the geological prior matrix, ω1 and ω2 are weighted coefficients.

[0137] The comprehensive scoring function is used to screen candidate solutions that better meet the inversion requirements in the sub-region to form the initial solution set Θ0 of variational reasoning.

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

[0139] Let the magnetic susceptibility component be m(θ) and the alteration degree component be e(θ). When constructing the prior p(θ), the penalty term Epen(θ) is added:

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

[0141] Among them, g(e(θ)) is the empirical function of alteration degree and magnetic susceptibility, κ is the penalty intensity constant, E pen It amplifies exponentially as the difference between the two increases.

[0142] Integrating this penalty term into p(θ) and affecting p(θ,X) can suppress solutions that are inconsistent with the geological prior during the variational inference update process.

[0143] S44: Combine the sub-region index k divided in S2 with the variational distribution Associate, let the learning rate η k (t) Perform differential updates on the sub-region at the tth iteration.

[0144] The update rules are defined as follows:

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

[0146] Among them, αk is the magnification coefficient corresponding to sub-region k, GradVal(Γ k ,t) represents the variational inference step t for Γ k Statistics of the internal magnetic field gradient or loss gradient.

[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 converges, q(θ; φ * ) gives an approximate posterior distribution of θ, which is used to describe the location, morphology and magnetization of the ore body.

[0149] To quantify the certainty of each dimension, the uncertainty measurement function U(θ d ), for the d-th dimension parameter θ d The variance in is calculated:

[0150]

[0151] Among them, θ d represents the dth dimension of the parameter, and Eθ[·] represents the Expectation operation of distribution.

[0152] Will U(θ d ) results are visualized in three-dimensional space, with uncertainty shown through isosurface mapping.

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

[0154] Visualization of the maximum a posteriori solution for the morphology distribution and magnetization, along with uncertainty ranges.

[0155] The final output structure is represented as:

[0156]

[0157] in, is the maximum a posteriori estimate of the ore body configuration and magnetization characteristics, Indicates the uncertainty quantification of 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 drilling points and close-range aeromagnetic data, and record the measured coordinates of all benchmark measuring points The magnetization intensity measured at close range A comparison sample is formed and the comprehensive error is calculated.

[0160] Define the comprehensive error E total Express the total amount of difference between two types:

[0161]

[0162] Where n represents the number of measured points, and λ1 and λ2 are the weight coefficients of position deviation and magnetization deviation. total Used to quickly determine the overall error level.

[0163] S52: If E total If the local error value in the data is significantly larger than the average value, detailed error distribution information is collected for this sub-area, including indicator functions of ore body boundaries and alteration zones, magnetic anomaly peaks, and lithology revealed by drill holes.

[0164] S53: Construct a regional correlation metric R(Δ, Π) to distinguish whether the deviation conflicts with the variational prior or is caused by 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 the position.

[0166] If R is high, the inferred bias concentration does not match the prior constraints;

[0167] If R is low, it indicates that there is something wrong with 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(Δ, Π) is much larger than the average level in the sub-region, it can be judged that the deviation is mainly caused by the conflict of prior constraints. Otherwise, it is related to insufficient sparrow search.

[0170] S54: After determining that the search is insufficient, the deviation distribution function in the sub-region Γ is recorded as ζ(x, y, z). In order to improve the focus of the sparrow algorithm on the region, the warning factor Λ is cumulatively corrected.

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

[0172] Λ new =Λ old +∫∫∫ Γ ζ(x,y,z)dxdydz

[0173] Among them, Λ old Represents the warning value of the previous sparrow search algorithm, and the triple integral on the right reflects the total deviation in the sub-region. After the calculation is completed, return to the sparrow search in step S3 and use Λ new Increase the intensity of global jumps or key area searches to ensure that areas that were originally under-searched receive higher weights in re-iterations.

[0174] S55: When the deviation is found to be due to a priori conflict according to the correlation metric R(Δ, Π) of S53, the priori constraint correction process is initiated to update the penalty strength κ of the “alteration degree-magnetic susceptibility” correlation function and other hyperparameters.

[0175] Define a new penalty term

[0176]

[0177] in, represents the original prior penalty function, μ represents the adjustment coefficient, μδ prior Represents the prior mismatch measure in the sub-region. Go 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 observation data.

[0178] S56: After completing the modification of the sparrow search process or the prior constraints, the three-dimensional inversion is re-performed using the updated search mechanism to obtain a new result G3.

[0179] The method of S51 is then repeated to compare with the drilling and close-range measurement data, and the deviation values ​​of position and magnetization intensity are calculated to ensure that the new results are significantly improved in the sub-areas where the errors were previously serious.

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

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

[0182] The specific steps include:

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

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

[0185]

[0186] Among them, k traverses all sub-areas from 1 to K, B k (x, y, z) takes a larger value in sub-region k and tends to zero in the uncovered area.

[0187] The weighted splicing mentioned above ensures a smooth transition at the regional boundaries and avoids multiple solution conflicts within the same grid cell. The merged G4(x, y, z) inherits the optimized and corrected inversion results of different sub-regions.

[0188] S62: To highlight the weak abnormal signals in the large ore body cluster area, based on the statistical analysis of G4 (x, y, z), the magnetization intensity peak and the middle section of the alteration zone are mapped to the coordinate set R main , in which the spatial discretization accuracy is improved. Define the local grid scale function h k (x,y,z) to dynamically adjust the mesh size:

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

[0190] Where h0 is the initial grid size, α is the magnification factor, and GradVal(F(x,y,z)) represents the statistical value of the magnetic field gradient. The gradient is often larger in the main ore body area, resulting in a smaller h k (x,y,z) to achieve local fine segmentation.

[0191] S63: In the exploration of skarn-type ore bodies, multiple periods of mineralization will overlap in time and space to form a complex pattern, and it is necessary to perform deep conjugation processing on easily confused areas. Let the variational inference posterior distribution output in step S4 be q(θ; φ * ), and define the sub-region where multiple periods may overlap as Ω m .

[0192] In order to reduce the conflict between the observed data and the prior distribution, the new distribution q'(θ) is compared with the original distribution q(θ; * ) is denoted as D KL :

[0193]

[0194] If D KLIf the value is always too high, it means that the current description of the multi-period superposition area is seriously deviated from the original prior and subsequent results. In this case, the parameters of q'(θ) are gradually corrected to make it more integrated with the original variational results, thereby achieving 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 local updates in high-confidence areas to avoid wasting computing resources.

[0197] In the low confidence area, the global search is restarted to prevent missing potential ore bodies. Define the partition update function Γ update (x,y,z):

[0198]

[0199] Among them, δ represents the confidence threshold, Γ local and Γ global They represent the local update strategy and the global re-search strategy respectively.

[0200] With the help of this partition update mechanism, we can balance high precision and wide coverage and enhance the exploration of uncertain areas.

[0201] S65: When the magnetic susceptibility or ore body morphology found in the local grid is significantly different from the geological interpretation, a new prior weight W is introduced. prior (x, y, z) amplifies the geological constraints in the region, prompting the variational inference to converge to a reasonable solution faster.

[0202] Define the updated prior weights for:

[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, The significant increase helps constrain the local model during variational inference and enhances geological rationality.

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

[0207] In order to take into account both the large-scale overview and the small-scale details, the multi-scale fusion expression in the 3D exploration model is defined:

[0208]

[0209] in, Represents the output result of the lth scale, ω l is the fusion weight of the results at different scales, and L is the number of multi-scale layers. It strengthens the small-scale information of the concealed ore body and suppresses 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, from G final Extract the magnetic susceptibility distribution χ of each sub-region from (x,y,z) k (x,y,z), and perform unified index mapping according to the spatial position.

[0212] In order to preserve the contrast of low, medium and high magnetization features in the final output, the merging function M is defined as total (x,y,z):

[0213]

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

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

[0216]

[0217] Where N represents the number of grid cells, I mine It is used to identify the interface between the ore body and the surrounding rock, and U represents the confidence level of the location.

[0218] The matrix is ​​visualized to provide intuitive spatial structure and credibility distribution of the ore body.

[0219] S73: In the previous inversion and search iterations, several rounds of possible solutions are accumulated for each three-dimensional grid, and a probability evaluation value is obtained based on the variational distribution; in order to determine the priority drilling location in S7, the comprehensive probability P of multiple iterations is introduced joint(x, y, z) and obtain the probability evaluation value, which is synthesized using the following formula:

[0220]

[0221] Where T represents the total number of iterations, pt(x,y,z) represents the probability of the existence of the ore body at coordinates (x,y,z) after the tth iteration, and α is the amplification factor, which is used to emphasize the area consistently pointed to by multiple rounds of searches.

[0222] If P joint The higher (x, y, z) is, the more valuable it is to prioritize drilling there.

[0223] S74: Determine the priority drilling locations and obtain several potential sections through interval classification; define the potential metric D pot (x,y,z), the magnetic susceptibility M total , morphological indicator function I mine And the comprehensive probability P joint Combined 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, β3 are weighted coefficients. pot Size divides the area into high potential area, medium potential area and marginal potential area.

[0226] S75. Obtain a unified output data set O according to S71 to S74. final , which contains the magnetic susceptibility distribution data M total , ore body morphology and uncertainty matrix R, priority drilling coordinate point set, and spatial index of classified potential areas. The whole can be expressed as:

[0227]

[0228] in, Refers to the set of priority drilling locations, D pot [zones] indicates the distribution and weight of each potential zone.

[0229] The specific embodiment of the present invention also provides a 3D inversion system for weak aeromagnetic anomalies of skarn iron ore, such as Figure 2 As shown, the system includes:

[0230] The first acquisition unit 201 is configured to acquire aeromagnetic data for 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;

[0231] The second acquisition unit 202 is used to identify the contact position between the alteration zone and the surrounding rock based on the three-dimensional grid and the prior matrix and divide the sub-regions; apply the three-dimensional gradient and difference operators to each sub-region, extract the field value mutation points and identify the weak magnetic anomalies, and then record them in the local anomaly sequence; and generate and store the results after introducing constraints for each sub-region;

[0232] The third acquisition unit 203 is used to perform a global search and group iterative update based on the sparrow search algorithm for the possible metallogenic space in the sub-region, and finally obtain a candidate solution set;

[0233] A fourth acquisition unit 204 is configured to iterate the candidate solution set as the initial condition of variational inference 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 applied to the sub-region during the iteration process; uncertainty is quantified by calculating the variance of the posterior distribution, and at the end of the iteration, a maximum a posteriori solution for the ore body position, morphology, and magnetization intensity and an uncertainty assessment result are output;

[0234] A fifth acquisition unit 205 is configured to compare the maximum a posteriori solution with known borehole point measurements and short-range aeromagnetic data to calculate a comprehensive error in position and magnetization intensity; and to re-execute iterative inversion based on the result until the error meets a preset threshold.

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

[0236] The seventh acquisition unit 207 is used to extract the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid from the three-dimensional exploration model, and calculate the potential section based on multiple rounds of search probability and alteration zone distribution characteristics; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results and priority drilling coordinates.

[0237] An embodiment of the present invention further provides a computer-readable storage medium, 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 any of the methods described above.

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

[0239] 1. This embodiment of the present invention performs refined coordinate unification and three-dimensional fusion of multiple periods of aeromagnetic data. By establishing an affine transformation matrix and threshold filtering mechanism for the differential data generated during each observation period, coordinate projection and anomaly removal are simultaneously processed, preserving key measurement points and reliable field values ​​within the same three-dimensional grid. Subsequently, drillhole information and geological catalog content are mapped to the same coordinate system, laying the foundation for subsequent regional segmentation and prior constraints, ensuring that multi-source information remains closely linked during the fusion process.

[0240] 2. The entire 3D domain is divided into several sub-regions. We then focus on calibrating weak magnetic anomalies within potential tectonic zones, combining lithologic boundaries, alteration sensitivity, and spatial gradient characteristics. Furthermore, we construct adjustable search weights and amplification operators based on the depth range and known lithologic locations within the region to highlight field signals at key depths or sensitive areas. This allows us to detect concentrated strong magnetic anomalies while also fully capturing relatively weaker local anomalies within the same 3D grid, providing more targeted input data for inversion and optimization.

[0241] 3. During the global search and variational inversion phases, the iterative update and jump mechanism of the sparrow search algorithm is used to cover possible ore bodies, and geological heuristic factors are combined to dynamically increase attention to points near the alteration zone. To ensure that the final inversion is consistent with the measured data both locally and globally, this embodiment of the present invention proposes a backtracking correction step, performing differentiated iterations on sub-regions where deviations occur. When search coverage is insufficient, the search intensity is increased and the scan is re-scanned. If a priori conflicts are determined, the coupling penalty function is adjusted to align with geological laws.

[0242] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0243] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple 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 device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0245] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0246] The foregoing descriptions of specific exemplary embodiments of the present invention are for purposes of illustration and description. These descriptions are not intended to limit the invention to the precise forms disclosed, and it is apparent that many variations and modifications are possible in light of the foregoing teachings. The exemplary embodiments have been selected and described for the purpose of explaining the specific principles of the invention and their practical application, thereby enabling those skilled in the art to realize and utilize a variety of exemplary embodiments of the invention and various options 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 iron ore, characterized in that: The method comprises: S1: Acquire aeromagnetic data from 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 the 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-regions; three-dimensional gradient and differential operators are applied to each sub-region to extract field value mutation points and identify weak magnetic anomalies, which are then recorded in a local anomaly sequence. Constraints are introduced for each sub-region to generate and store results; S3: For the possible mineralization space in the sub-area, a global search and group iterative update are performed based on the sparrow search algorithm to obtain a set of candidate solutions; S4: Iterating the candidate solution set as the initial conditions of variational inference 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 applied to the sub-region during the iteration process; uncertainty is quantified by calculating the variance of the posterior distribution, and at the end of the iteration, a maximum posterior solution and uncertainty assessment results of the ore body position, morphology and magnetization intensity are output; S5: Based on the maximum a posteriori solution, the solution is compared with the measured values ​​of known borehole points and the close-range aeromagnetic data to calculate the comprehensive error of the position and magnetization intensity; and the iterative inversion is re-performed according to the result until the error meets the preset threshold. S6: Perform multi-scale fusion on the corrected sub-region results 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 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 ): Where A is a two-dimensional affine transformation matrix, x i with y i Indicates 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 point after the unified coordinate system; S13: Obtain the normalized data set D1, perform validity checks on the magnetic field strength and coordinate information of each observation point; remove data blocks with missing value ratios exceeding a predetermined limit, and delete abnormal samples in high-amplitude impulse noise areas to obtain a preliminary cleaned data set D2; S14: performing three-dimensional interpolation on the preliminary cleaned data set D2 to obtain a three-dimensional grid G0; S15: mapping the drill hole lithology and geological log 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 of the contact between the alteration zone and the surrounding rock on the three-dimensional grid and the prior 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 multiple sub-regions with different attention levels; S23: scheduling and storing local magnetic field data in the sub-region, 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 completing the extraction of weak magnetic anomalies in the sub-area, a constraint condition set C is introduced based on the lithologic boundary and alteration zone information; S25: completing search constraints for each sub-area based on the constraint condition set C, and performing 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-area, unify the data packaging and generate storage results.

4. The method according to claim 1, wherein Said S3 further comprises: S31: determining the boundary of the sparrow search algorithm in the three-dimensional space based on the coordinate information of the weak magnetic anomaly sequence L in the stored result and the depth constraint z; S32: Construct an allowable search domain D corresponding to each possible boundary, conduct multi-dimensional initial group settings for key physical quantities in skarn ore bodies and combine them into individual vectors; S33: Applying geological heuristic factors to individuals during the iterative update process of the group; S34: Perform a jump operation during the group iterative update process; S35: Improving the coverage of multi-peak areas by changing the dynamic alert factor with 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 Said S4 further comprises: S41: In the candidate solution set, a scoring function is performed on all solutions based on the subregion partition index, and the top several solutions with the highest performance in each subregion are selected as initial conditions for variational inference; S42: Define a scoring function based on a comprehensive measure of magnetic field fitting and geological prior conformity, and use the comprehensive scoring function to screen candidate solutions in the sub-region to form an initial solution set Θ0 for variational reasoning.

6. The method according to claim 1, characterized in that The S5 further includes: The ore body position coordinates and magnetization intensity are extracted from the three-dimensional inversion results; compared with the measurement values ​​of known drill hole measurement points and close-range aeromagnetic data, the measured coordinates of all benchmark measurement points and the close-range measured magnetization intensity are recorded to form a comparison sample and the comprehensive error is calculated; iterative inversion is re-executed based on the results until the error meets the preset threshold, and the final three-dimensional distribution prediction of skarn-type iron ore is output.

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

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

9. A three-dimensional inversion system for weak aeromagnetic anomalies of skarn iron ore, characterized by: The system comprises: A first acquisition unit is configured to acquire aeromagnetic data for 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; The second acquisition unit is used to identify the contact position between the alteration zone and the surrounding rock based on the three-dimensional grid and the prior matrix and divide the sub-regions; apply the three-dimensional gradient and differential operators to each sub-region, extract the field value mutation points and identify the weak magnetic anomalies, and then record them in the local anomaly sequence; and generate and store the results after introducing constraints for each sub-region; The third acquisition unit is used to perform global search and group iterative update based on the sparrow search algorithm for the possible mineralization space in the sub-area, and obtain a set of candidate solutions; a fourth acquisition unit, configured to iterate the candidate solution set as an initial condition for variational inference 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 applied to the sub-region during the iteration process; uncertainty is quantified by calculating the variance of the posterior distribution, and at the end of the iteration, a maximum a posteriori solution for the ore body position, morphology, and magnetization intensity and an uncertainty assessment result are output; a fifth acquisition unit, configured to compare the maximum a posteriori solution with known borehole measurement values ​​and short-range aeromagnetic data to calculate a comprehensive error in position and magnetization intensity; and re-execute iterative inversion based on the result until the error meets a preset threshold; The sixth acquisition unit is used to perform multi-scale fusion on the corrected sub-region results and obtain a three-dimensional exploration model; The seventh acquisition unit is used to extract the magnetic susceptibility distribution and ore body morphology indicator function of each three-dimensional grid from the three-dimensional exploration model, and calculate the potential section based on multiple rounds of search probability and alteration zone distribution characteristics; output the three-dimensional exploration model and the corresponding ore body spatial position, uncertainty quantification results 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