Same-amplitude strong magnetic anomaly source discrimination method

Abstract

The invention discloses a same-amplitude strong magnetic anomaly source discrimination method, which comprises the steps of generating a mineral combination consistency graph, an aviation magnetic method characteristic grid set and an aviation gravity characteristic grid set in a contact zone working domain and a unified grid space reference, and obtaining a skarn compatibility weight graph and a residual magnetism direction feasible domain graph; by taking the skarn compatibility weight map and the residual magnetism direction feasible region map as constraints, embedding a gravity-magnetic vector magnetization joint inversion target function to obtain a density body, a magnetic susceptibility body and a residual magnetism direction body; the method comprises the following steps: dividing same-amplitude strong magnetic units into skarn ore bodies or basic / ultra-basic rock veins through a classification discrimination rule, calculating data residual consistency and prior consistency, and generating a source discrimination graph and a credibility graph; carrying out geometric consistency and homodromous change inspection along the direction of the contact zone, and screening out a priority zone; and generating a hole arrangement reference graph and a suggested drilling position and depth interval in the contact zone working area based on the priority zone. According to the scheme, the availability and deployment efficiency of the source judgment result can be improved.

Description

Technical Field

[0001] This invention relates to the field of geophysical applied geology, and in particular to a method for identifying the source of strong magnetic anomalies of the same amplitude. Background Technology

[0002] The contact zone between carbonate host rocks and intermediate-acidic intrusive bodies is a key control area for skarn mineralization. Hydrothermal-metastrophic processes readily form assemblages of garnet, pyroxene, epidote, and iron oxides, among which magnetite and other minerals tend to cause significant magnetic anomalies. Airborne magnetometry and gravity can reflect the magnetic susceptibility and distribution density of the physical body, while satellite / airborne hyperspectral imaging can identify mineral spectral absorption characteristics, providing remote sensing evidence for mineralization indicators. However, within the contact zone, the complex remanent magnetization direction and the superposition of multiple magmatic-tectonic phases result in diverse magnetic anomaly polarities and amplitudes. Strong magnetic fields of the same amplitude are common, potentially originating from economically valuable skarn ore bodies or from unmineralized basic to ultrabasic dikes. Furthermore, significant differences in spatial resolution and coverage between multi-source data mean that raster alignment and coordinate unification are prerequisites for pixel-level fusion. Connectivity and normal bandwidth along the contact zone strike also directly affect the expression and selection of anomalies. In addition, the coordinated changes in density and magnetic susceptibility, as well as the anisotropic structure along the strike, are also important characteristics that characterize the continuity of mineralization; if these are not reflected, the interpretation often deviates from the actual engineering requirements.

[0003] Existing technologies typically employ two approaches: ① Physical field interpretation primarily based on airborne magnetics / gravity. This involves enhancing boundaries and shape through processes such as real-time polarization (RTP), oblique derivatives, analytical signals, and downward extension. Density or magnetic susceptibility is then estimated using conventional L2 / L1 canonical gravity and magnetic inversion, often employing scalar magnetization assumptions and weak geometric constraints. Soft masks are sometimes applied to geological boundaries. ② Mineral identification primarily based on remote sensing spectroscopy. This utilizes characteristic absorption bands of minerals such as calcite, dolomite, and iron oxides to generate alteration mineral distribution maps, which are then spatially overlaid with magnetic anomalies for analysis. However, differences in spatial resolution and coverage exist between multi-source data. Raster alignment and coordinate unification are prerequisites for pixel-level fusion. Existing methods often remain at the level of empirical interpretation or simple machine learning classification, failing to transform mineral assemblage information into spatial constraints during the inversion process. The above-mentioned schemes either use geological-remote sensing evidence as a posteriori interpretation label, or only introduce general smoothing / sparse constraints in the inversion, making it difficult to transform mineral assemblage evidence and remanent magnetic direction information along the contact zone into computable constraints in the inversion; at the same time, the segmented consistency test along the contact zone and the engineering output (preferred zone, borehole depth) often do not form a closed loop.

[0004] In summary, existing technologies suffer from common shortcomings in contact zone scenarios: they lack scenario-based constraints for distinguishing the sources of strong magnetic anomalies of the same amplitude; they cannot incorporate the compatibility of mineral assemblages and contact zone geometric control into the spatial weights of the inversion; and they cannot convert the remanent magnetic direction information implied by the anomaly polarity / offset into a pixel-level directional feasible region and force its projection during iteration. The inversion solutions are easily affected by non-uniqueness and extra-domain extensions; classification often relies on amplitude thresholds or empirical superposition; and there is a lack of verifiable mechanisms for continuity and unidirectional changes along the contact zone, making it difficult to stably output discrimination results that can be directly used for priority zone labeling and aperture design. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying the source of strong magnetic anomalies of the same amplitude.

[0006] To achieve the above objectives, the present invention provides a method for identifying the source of strong magnetic anomalies of the same amplitude, comprising the following steps:

[0007] S1: Based on the geological boundary and tectonic line of the contact zone between carbonate rocks and intrusive bodies, and integrating the combined coverage of airborne magnetics, airborne gravity, airborne hyperspectral and satellite hyperspectral methods, a working domain for the contact zone and a unified grid spatial benchmark are constructed.

[0008] S2: Based on the contact zone working domain and unified grid spatial reference, spectral interpretation of airborne hyperspectral and satellite hyperspectral is performed to identify garnet, pyroxene, epidote and iron oxide and distinguish them from carbonate rock indicators, generate a mineral assemblage consistency map, and extract airborne magnetic and airborne gravity features to generate airborne magnetic feature grid sets and airborne gravity feature grid sets.

[0009] S3: Based on the mineral assemblage consistency map, the airborne magnetic characteristic grid set, and the airborne gravity characteristic grid set, a multi-flow heterogeneous map causal reasoning model is constructed, and the causal direction is transmitted from source type to physical properties to observation anomalies to determine the skarn compatibility weight map and the feasible domain map of remanent magnetization direction.

[0010] S4: Using the skarn compatibility weight map as spatial weight and the remanent magnetic direction feasible domain map as directional constraint, embed the gravity and magnetic vector magnetization joint inversion objective function to determine the density volume, magnetic susceptibility volume and remanent magnetic direction volume;

[0011] S5: Based on the density volume, magnetic susceptibility volume and remanent direction volume, skarn compatibility weight map and remanent direction feasible domain map, the same strong magnetic units are divided into skarn ore bodies or basic or ultrabasic dikes by classification discrimination rules. The data residual consistency and prior consistency are calculated, and the source discrimination map and confidence map are generated.

[0012] S6: Based on the source discrimination map, confidence map, and mineral assemblage consistency map, conduct geometric consistency and unidirectional variation tests along the contact zone direction to generate the priority zone;

[0013] S7: Priority band superimposed source discrimination map and confidence map. Based on the spatial relationship between density volume, magnetic susceptibility volume and remanent direction volume, a hole layout reference map is given in the working domain of the contact zone, and suggested drilling positions and depth ranges are generated.

[0014] In one embodiment of the present invention, S1 uses the geological boundary and tectonic line of the contact zone between carbonate rocks and intrusive bodies as a reference, and integrates the combined coverage of airborne magnetics, airborne gravity, airborne hyperspectral, and satellite hyperspectral to construct a working domain for the contact zone and a unified grid spatial reference, including:

[0015] Based on geological maps and geophysical interpretation results, the contact zone boundary lines and tectonic lines are extracted and unified into the same coordinate system;

[0016] A buffer zone is set along the normal direction with the contact zone boundary line as the center to form the initial contact zone working area;

[0017] The common coverage area and the initial contact zone working domain are spatially intersected to obtain the contact zone working domain;

[0018] Establish a unified grid spatial reference within the contact zone working domain.

[0019] In one embodiment of the present invention, S2 involves spectral interpretation of airborne hyperspectral and satellite hyperspectral data to identify garnet, pyroxene, epidote, and iron oxides, distinguishing them from carbonate rock indicators, and generating a mineral assemblage consistency map, including:

[0020] Diagnostic absorption feature identification and mineral indicator index calculation are performed on the reflectance sequences of airborne hyperspectral and satellite hyperspectral data, and the mineral indicator weights are output in pixels.

[0021] For garnet, pyroxene, epidote and iron oxide, respectively, diagnostic band combinations are selected, mineral indicator indices are calculated and corresponding mineral indicator maps are output;

[0022] Convert each mineral indicator map into a pixel-level indicator weight according to a unified dimension;

[0023] A carbonate rock indicator mask was constructed, and the indicator weights of each pixel were separately calibrated inside and outside the carbonate rock indicator mask to obtain the mineral indication results distinguished by carbonate rock indicators.

[0024] Based on the mineral indication results distinguished by carbonate rock indicators, a mineral assemblage consistency map was determined.

[0025] In one embodiment of the present invention, S3 constructs a multi-flow heterogeneity map causal reasoning model based on the mineral assemblage consistency map, the airborne magnetic characteristic grid set, and the airborne gravity characteristic grid set, and determines the skarn compatibility weight map and the remanent magnetization direction feasible domain map along the causal direction from source type to physical properties to observation anomalies, including:

[0026] Within a unified raster space reference, establish source type nodes, physical property nodes, observation anomaly nodes, and contact zone geometry nodes for each pixel;

[0027] At the source type node, the prior probability of the source type being skarn is set based on the mineral assemblage consistency map and contact zone geometry. At the property node, prior intervals for density increment, magnetic susceptibility increment, and remanence compatibility with skarn are set. The observation likelihood is calculated at the observation anomaly node using airborne magnetic and airborne gravity characteristics, and the posterior of the property node and the source type node is updated from bottom to top along the directed edges. Consistency propagation is performed between pixels through the spatial adjacent edges controlled by the contact zone geometry node to obtain the posterior distribution of the source type and the posterior of the property related to the remanence direction.

[0028] Using the posterior probability of skarn as the source type as the pixel-level weight, the pixel-level weights are used to output a skarn compatibility weight map within a unified raster space reference.

[0029] For each pixel, the feasible domain map of the remanent magnetization direction is determined based on the characteristics and source type of airborne magnetics and the a posteriori of physical properties.

[0030] In one embodiment of the present invention, the objective function for the joint inversion of gravity and magnetic vector magnetization includes gravity observation residual terms, magnetic field observation residual terms, structural constraint terms, and orientation constraint terms.

[0031] S4 uses the skarn compatibility weight map as spatial weight and the remanent magnetic direction feasible region map as directional constraint, embedding a joint inversion objective function of gravity and magnetic vector magnetization to determine the density volume, magnetic susceptibility volume, and remanent magnetic direction volume, including:

[0032] Spatial weight coefficients are generated at the pixel level using the skarn compatibility weight map, and these spatial weight coefficients are used for weighted structural constraint terms and parameter boundary constraints.

[0033] The direction allowable set is generated at the cell level using the feasible domain map of the remanent magnetization direction, and the direction allowable set is written into the direction constraint term;

[0034] Within a unified grid spatial reference, the density volume, magnetization volume, and remanent orientation volume are initialized. The weighted sum of gravity and magnetic field observation residuals, adjusted by spatial weights, is iteratively minimized, and structural constraints are applied. After each update, the remanent orientation volume is projected at the pixel level onto the feasible remanent orientation region map within the allowable azimuth and dip angles given by the orientation constraints. Iteration continues under the mask control of geometric hard constraints until the joint inversion objective function of gravity and magnetic vector magnetization converges. The resulting output is a skarn-compatible weighted map and a feasible remanent orientation region. Figure 1 Density volume, magnetic susceptibility volume, and remanent direction volume.

[0035] In one embodiment of the present invention, S5 classifies strong magnetic units of the same amplitude into skarn ore bodies or basic or ultrabasic dikes based on density volume, magnetic susceptibility volume and remanent magnetic direction volume, skarn compatibility weight map and remanent magnetic direction feasible domain map, using classification discrimination rules, including:

[0036] Based on the density volume, magnetic susceptibility volume, and remanent direction volume, the unit mean and quantile statistics of strong magnetic units of the same amplitude are calculated within a unified grid space reference, and the coincidence rate of density index, magnetic susceptibility index, remanent direction volume and remanent direction feasible domain map is obtained.

[0037] Based on the skarn compatibility weight diagram, the average value of each element is calculated within the element.

[0038] Based on the density index, magnetic susceptibility index, and the coincidence rate and unit average of the remanent direction volume and the feasible domain diagram of the remanent direction, the classification and discrimination rules are determined.

[0039] According to the classification and discrimination rules, the same strong magnetic units are divided into skarn ore bodies or basic or ultrabasic dikes.

[0040] In one embodiment of the present invention, according to classification rules, strong magnetic units of the same amplitude are divided into skarn ore bodies or basic or ultrabasic dikes, including:

[0041] When the average value of the unit in the skarn compatibility weight map is greater than the first threshold and the coincidence rate between the remanent magnetic direction body and the remanent magnetic direction feasible domain map is greater than the second threshold, and the magnetic susceptibility index reaches the third threshold, then the strong magnetic unit with the same amplitude is judged as a skarn ore body.

[0042] When the average value of the skarn compatibility weight map is less than the first threshold and the density index and magnetic susceptibility index are both greater than the fourth threshold, and the remanent magnetic direction volume matches the remanent magnetic direction feasible domain map less than the second threshold, then the same strong magnetic unit is judged as a basic or ultrabasic dike.

[0043] In one embodiment of the present invention, S5 calculates the consistency of data residuals and prior knowledge, and generates a source discrimination graph and a confidence graph, including:

[0044] For each strong magnetic unit of the same amplitude, gravity prediction and magnetic field prediction are obtained by forward modeling density volume, magnetic susceptibility volume and remanent direction volume, and compared with the corresponding observations to obtain the unit-level consistency score of gravity residual and magnetic field residual, which is defined as data residual consistency.

[0045] Based on the source type label, calculate the average value of the skarn compatibility weight map in the corresponding unit and the consistency rate between the remanent magnetic direction volume and the remanent magnetic direction feasible domain map, and combine them into prior consistency.

[0046] Write the source type label into the corresponding cell of the source discrimination graph;

[0047] The consistency of data residuals and the consistency of prior knowledge are combined into a confidence score according to a preset weight and written into the corresponding cell of the confidence map.

[0048] In one embodiment of the present invention, step S6 involves performing geometric consistency and unidirectional variation checks along the contact zone direction based on the source discrimination map, confidence map, and mineral assemblage consistency map, including:

[0049] Using the contact zone geometry as a linear reference, a linear index is established along the contact zone geometry within a unified grid space reference, and the contact zone geometry is discretized into several segments of equal length.

[0050] Within each segment, the pixel proportion of skarn orebody tags is extracted from the source discrimination map, the segment-level average value is calculated from the mineral assemblage consistency map, the segment-level average value is calculated from the confidence map, and the angle between the main direction of the skarn source pixel and the geometric orientation of the contact zone is calculated within the corresponding segment.

[0051] The geometric consistency condition is defined as the angle between the main direction and the geometric orientation of the contact zone not exceeding the first threshold and the proportion of skarn ore body labels not less than the fifth threshold.

[0052] The condition for consistent change is defined as the skarn orebody label ratio, mineral assemblage consistency segment-level average value, and credibility segment-level average value of the changes in two adjacent segments being consistent, and the three not exhibiting contradictory abrupt changes between adjacent segments.

[0053] In one embodiment of the present invention, S6 generates a priority band, including:

[0054] Segments that simultaneously satisfy both the geometric consistency condition and the same direction of change condition are marked as consistent segments;

[0055] Connectivity analysis is performed on the set of consistent segments according to the contact zone geometry, and the maximum continuous sequence formed by adjacent consistent segments is extracted. Each maximum continuous sequence is defined as a continuous skarn source segment.

[0056] For each continuous skarn source segment, the sequence length is calculated and compared with a first length threshold. Only sequences with a length not less than the first length threshold are retained.

[0057] Within the retained sequence, verify whether the segment-level average value of the confidence map and the segment-level average value of the mineral assemblage consistency map are not lower than the sixth threshold and the seventh threshold, respectively. Then, mark the continuous skarn source segment that passes the verification as a segment-level priority zone.

[0058] Centered on the linear labeling of the segment-level priority zone, a set of strip-shaped pixels that are consistent with the normal range of the contact zone working domain are selected within the unified raster spatial reference as the spatial coverage of the priority zone. The set of strip-shaped pixels corresponding to each continuous skarn source segment is assigned a priority zone identifier code, and the priority zone is output within the contact zone working domain.

[0059] This invention discloses a method for identifying sources of strong magnetic anomalies of the same amplitude. The method first constructs a unified grid and common coverage based on the contact zone geometry as the main axis. Then, it defines strong magnetic units of the same amplitude on the joint high amplitude of observation and prediction, and classifies them using unit statistics from density volume, magnetic susceptibility volume, and remanent direction volume as a priori. Finally, it verifies the geometric consistency and unidirectional variation along the contact zone and outputs a preferred zone and borehole reference map. This scheme connects the spatial framework, prior construction, inversion solution, and segment-level quality control in a closed loop. It no longer relies on a single amplitude threshold for decision-making, but instead uses joint evidence of material properties, direction, priors, and residual consistency to constrain misjudgments. Under hard geometric constraints, it avoids model volume out-of-bounds errors, ensures traceability of results along the direction, and outputs engineering results including suggested borehole locations and suggested depth ranges, improving the usability and deployment efficiency of the source identification results.

[0060] This invention uses the skarn compatibility weight map obtained through causal reasoning as spatial weights, jointly adjusting observation residuals, anisotropic structure regularization, and parameter limits; it uses the remanent magnetization direction feasible domain map as pixel-level directional constraints and performs projection after each parameter update; it implements geometric hard constraints and strike / normal differential constraints within the contact zone working domain, and introduces density-magnetic susceptibility boundary coordination. Unlike joint inversion using equal weights or uniform regularization, this method employs source weight-driven spatial variable constraints and directional feasible domain projection. The improved technical effect is that, without relying on a preset magnetization direction, it improves the identifiability of remanent magnetization direction and magnetic susceptibility, reduces the non-uniqueness and extra-domain expansion risk of vector magnetization inversion, and ensures that the results maintain reasonable continuity along the contact zone direction and are geometrically compatible with mineralization, providing a stable physical property basis for distinguishing the genesis of strong magnets of the same amplitude.

[0061] This invention organizes mineral assemblage consistency, aeromagnetic / gravity characteristics, and contact zone geometry based on the causal orientation of source type, physical properties, and observed anomalies. It combines this with strike-adjacent propagation to obtain the source posterior, and derives a skarn compatibility weight map and a feasible remanent magnetic direction map. The feasible direction region is obtained through joint screening using polarity morphology and anomaly peak-valley shift patterns. The technical effect of this method is that it quantifies remote sensing mineralogical and geometric evidence into spatial weights and directional constraints that can directly drive inversion. This differs from traditional interpretations and classifications based solely on amplitude or boundary enhancement. In low-latitude and strong remanent magnetic conditions, it provides more reliable directional and source priors, reduces confusion between bodies with similar amplitudes, and provides a quantifiable constraint basis for subsequent inversion and classification. Attached Figure Description

[0062] Figure 1 This is a flowchart illustrating the method for identifying the source of strong magnetic anomalies of the same amplitude provided by the present invention. Detailed Implementation

[0063] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings, but it should be understood that the scope of protection of the present invention is not limited to the specific embodiments.

[0064] Unless otherwise expressly stated, throughout the specification and claims, the term "comprising" or its variations such as "including" or "comprises" shall be understood to include the stated elements or components without excluding other elements or other components.

[0065] like Figure 1 As shown, a method for identifying the source of a strong magnetic anomaly of the same amplitude according to a preferred embodiment of the present invention includes the following steps:

[0066] S1: Based on the geological boundary and tectonic line of the contact zone between carbonate rocks and intrusive bodies, and integrating the combined coverage of airborne magnetics, airborne gravity, airborne hyperspectral and satellite hyperspectral methods, a working domain and a unified grid spatial benchmark for the contact zone are constructed as a coordinate framework for subsequent data and prior knowledge.

[0067] Specifically, geological boundaries and tectonic lines include intrusive contact boundaries, carbonate rock boundaries, major faults, shear zones, and fold traces, and are represented by vector lines from geological maps and geophysical interpretations to define the geometric characteristics of the contact zone.

[0068] The common coverage is the spatial intersection area of ​​airborne magnetics, airborne gravity, airborne hyperspectral and satellite hyperspectral in the same coordinate system, and is limited to the area where each data valid pixel exists at the same time.

[0069] The unified raster spatial reference includes coordinate system, projection, cell size, starting point and row and column numbers. The raster coverage is consistent with the contact zone working domain, which is used to ensure that multi-source data and priors are aligned at the cell level.

[0070] The coordinate framework used for subsequent data and prior knowledge refers to mapping the mineral assemblage consistency map, skarn compatibility weight map, remanent magnetic direction feasible domain map, density volume, magnetic susceptibility volume, and remanent magnetic direction volume to a unified raster spatial reference to ensure spatial consistency.

[0071] Based on the geological boundary and tectonic line of the contact zone between carbonate rocks and intrusive bodies, a working domain and a unified grid spatial reference are constructed for the contact zone, including:

[0072] S1.1. Using geological maps and geophysical interpretation results as sources, extract the contact zone boundary lines and structural lines and unify them into the same coordinate system to obtain linear elements for geometric constraints;

[0073] S1.2. A buffer zone is set up along the normal direction with the boundary line of the contact zone as the center. The width of the buffer zone is determined according to the mineralization control range of the skarn in the region and is continuously connected along the direction of the contact zone to form the initial contact zone working area.

[0074] S1.3 Perform spatial intersection calculations on the common coverage of airborne magnetics, airborne gravity, airborne hyperspectral and satellite hyperspectral and the initial contact zone working domain, eliminate non-common coverage areas to obtain the final contact zone working domain. The final contact zone working domain is used to limit the range of causes of strong magnetic anomalies of the same amplitude and to provide domain boundaries for subsequent geometric hard constraints.

[0075] S1.4. Establish a unified raster spatial reference within the final contact zone working domain, with the raster long axis laid out along the contact zone direction. The pixel size is selected according to the dominant scale in the common coverage, and the starting point is fixed at the corner of the contact zone working domain, forming a raster framework with unique row and column numbers. This aligns the rasterization results of airborne magnetism, airborne gravity, airborne hyperspectral, and satellite hyperspectral to the same pixel, thus providing unified coordinate support for the subsequent generation of mineral assemblage consistency maps and the construction of skarn compatibility weight maps and remanent magnetic direction feasible domain maps. It also provides a unique solution domain and pixel-level constraints for the joint inversion of gravity and magnetic vector magnetization, enabling the discrimination of strong magnetic anomalies of the same amplitude to be carried out within the contact zone working domain and avoiding extra-domain expansion.

[0076] S2. Based on the contact zone working domain and a unified grid spatial reference, spectral interpretation of airborne hyperspectral and satellite hyperspectral is performed to identify garnet, pyroxene, epidote and iron oxides and distinguish them from carbonate rock alteration. Masking, resampling and normalization are used to generate a mineral assemblage consistency map and extract airborne magnetic and airborne gravity features.

[0077] Specifically, spectral interpretation refers to the diagnostic absorption feature identification and mineral indicator index calculation of reflectance sequences of airborne hyperspectral and satellite hyperspectral data within the contact zone working domain and a unified grid spatial reference, outputting mineral indicator values ​​in pixels.

[0078] Distinguishing from carbonate rock indicators refers to constructing a carbonate indicator mask based on the diagnostic absorption characteristics of calcite and dolomite, and distinguishing pixels inside the carbonate indicator mask from pixels outside the carbonate indicator mask, in order to define the alteration zone of skarn.

[0079] Extracting airborne magnetic and airborne gravity features involves calculating derived quantities such as amplitude, gradient, and boundary enhancement from airborne magnetic and airborne gravity data under a unified grid spatial reference, forming a feature grid for subsequent processing.

[0080] The identification of garnet, pyroxene, epidote, and iron oxides, and their differentiation from carbonate rock indicators, is specifically as follows: Within the contact zone working domain and a unified grid spatial reference, diagnostic band combinations are selected for garnet, pyroxene, epidote, and iron oxides, respectively. Mineral indicator indices are calculated, and corresponding mineral indicator maps are output. Each mineral indicator map is converted into a pixel-level indicator weight according to a unified dimension. A carbonate rock indicator mask is constructed, and the indicator weights of each mineral are calibrated both inside and outside the carbonate rock indicator mask. The weights inside the carbonate rock indicator mask are used to express the skarn alteration probability within the carbonate rock, while the weights outside the carbonate rock indicator mask are used to identify the non-carbonate rock background. The resulting mineral indicator, differentiated by carbonate rock indicators, serves as input for generating a mineral assemblage consistency map.

[0081] The generation of the mineral assemblage consistency map aims at the compatibility of the skarn mineralization process. Specifically, it uses the indicator weights of garnet, pyroxene, epidote, and iron oxide, distinguished by carbonate rock indicators, as inputs. Pixel-level weight coefficients are set according to the co-occurrence relationships of mineral assemblages in skarn mineralization, ensuring that the weights of iron oxides, garnet, and pyroxene are higher than those of epidote. A weighted consistency score is calculated for each pixel and normalized to the range of zero to one, resulting in the mineral assemblage consistency map. This map only expresses the strength of compatibility with skarn mineralization and does not indicate mineral content or thickness. Furthermore, the mineral assemblage consistency map is screened for connectivity based on the geometric orientation of the contact zone working domain, retaining only the high-consistency bands continuous along the contact zone as input for S3 to support the construction of spatial weights in subsequent multi-flow heterogeneous graph causal inference and vector magnetization joint inversion.

[0082] S3, using the mineral assemblage consistency map, airborne magnetic features, and airborne gravity features output by S2, along with the contact zone geometry output by S1, as input, constructs a multi-flow heterogeneous graph causal reasoning model, and propagates it along the causal direction from source type to physical properties to observational anomalies, outputting a skarn compatibility weight map and a feasible domain map of remanent magnetization direction.

[0083] Specifically, the multi-flow heterogeneous graph causal reasoning model includes node types, edge types, causal directions, and propagation rules. Among them, node types include mineral assemblage consistency nodes, aeromagnetic feature nodes, aerogravity feature nodes, and contact zone geometric nodes. Edge types include spatial adjacency edges of similar nodes and cross-layer dependency edges of dissimilar nodes.

[0084] The causal propagation along the source type to the property to the observation anomaly refers to setting directed edges in the multi-flow heterogeneous graph causal inference model, with the source type node pointing to the property node and the property node pointing to the observation anomaly node, and using this direction for prior setting and posterior update.

[0085] The skarn compatibility weight map is a raster map output with the contact zone working domain as the scope and a unified raster spatial reference as the pixel unit. The pixel value represents the weight of the pixel's compatibility with the skarn source type.

[0086] The feasible domain map of remanent magnetization direction is a directional constraint map output with the working domain of the contact zone as the range and a unified grid space reference as the pixel unit. The pixel value represents the allowable range of the remanent magnetization azimuth and tilt angle of that pixel.

[0087] The multi-source fusion reasoning model is a multi-source fusion reasoning route aimed at identifying the source of strong magnetic anomalies of the same amplitude. Its construction steps are as follows: taking the mineral assemblage consistency map, airborne magnetic features, and airborne gravity features output by S2, along with the contact zone geometry output by S1, as input, source type nodes, physical property nodes, and observation anomaly nodes are established for each cell within a unified grid spatial reference, and contact zone geometry nodes are set to control spatial correlation; within the cell, mineral assemblage consistency nodes, airborne magnetic feature nodes, and airborne gravity feature nodes are connected to physical property nodes by cross-layer dependent edges, and then physical property nodes are connected to observation anomaly nodes; spatial adjacency edges are set between cells along the direction of contact zone geometry, and the weight along the contact zone direction is increased to reflect the continuity of mineralization, thereby forming a graph structure that can perform source type reasoning and direction constraint generation within the contact zone working domain.

[0088] The causal propagation along the path from source type to physical properties to observed anomalies involves setting the prior from top to bottom and correcting the posterior from bottom to top. Specifically, at the source type node, the prior probability of the source type being skarn is first set based on the mineral assemblage consistency map and contact zone geometry. At the physical property node, prior intervals for density increment, magnetic susceptibility increment, and remanence compatibility with skarn are set. Then, the observation likelihood is calculated at the observed anomaly node using airborne magnetic and airborne gravity characteristics, and the posterior of the physical property node and the source type node is updated from bottom to top along the directed edges. Simultaneously, consistency propagation is performed between pixels through the spatial adjacent edges controlled by the contact zone geometry node, ensuring that the posterior remains reasonably continuous along the contact zone direction. This yields the posterior distribution of the source type and the physical property posterior related to the remanence direction.

[0089] The skarn compatibility weight map is a generation route that maps the posterior of the source type to spatial weights. Its specific technical route can be understood as follows: the posterior probability of the source type being skarn obtained by causal direction transmission is used as the pixel-level weight, and this weight is output as a skarn compatibility weight map within a unified raster spatial reference; high weights are maintained for pixels that are geometrically continuous along the contact zone, and the weights are reduced for pixels that do not meet the mineral assemblage consistency requirement, so as to form a skarn compatibility weight map that can be directly used as S4 as spatial weight.

[0090] The feasible domain map of remanent magnetization direction is the generation route constrained by observed anomaly features for the range of remanent magnetization azimuth and dip angles. Specifically, for each pixel, a set of discrete azimuth and dip angle candidates are constructed by combining the polarity, anomaly offset, and dip angle derivative sign pattern in the aeromagnetic features with the source type and posterior properties. The consistency score between the candidate direction and the observed anomaly is evaluated, and the azimuth and dip angle ranges with scores exceeding the threshold are selected as the allowable range for that pixel. The consistency of adjacent pixels is then screened along the contact zone geometry, and the output is the feasible domain map of remanent magnetization direction, which is used by S4 for directional constraints.

[0091] S4. Using the skarn compatibility weight map and remanent magnetic direction feasible domain map output from S3 as input, the skarn compatibility weight map is used as spatial weight, and the remanent magnetic direction feasible domain map is used as directional constraint to embed into the gravity and magnetic vector magnetization joint inversion objective function. Geometric hard constraints are implemented with the contact zone working domain and the unified grid spatial reference to solve the density volume, magnetic susceptibility volume and remanent magnetic direction volume.

[0092] Specifically, spatial weight refers to the weight coefficient field set for pixels within a unified raster spatial reference. It is used to adjust the relative influence of the observation residual, structural constraint, and directional constraint terms in the inversion objective function on that pixel. Its source is the pixel value of the skarn compatibility weight map.

[0093] Directional constraints refer to the allowable range of azimuth and tilt angles set at the pixel level based on the feasible domain map of remanent magnetization direction, which are used to limit the value of the remanent magnetization direction in vector magnetization.

[0094] The objective function for the joint inversion of gravity and magnetic vector magnetization is constructed using gravity and magnetic field data as observations and density volume, magnetic susceptibility volume, and remanent direction volume as unknowns. It includes gravity observation residuals, magnetic field observation residuals, structural constraints, and directional constraints.

[0095] Geometric hard constraints refer to using the contact zone working domain as the solution domain boundary, fixing the model parameters of pixels outside the contact zone working domain as background values ​​and not participating in the solution, thereby prohibiting the model volume from expanding beyond the boundary.

[0096] A density volume refers to a voxelized density distribution field under a uniform grid space reference. A magnetic susceptibility volume refers to a voxelized magnetic susceptibility distribution field under a uniform grid space reference. A remanent direction volume refers to a voxelized remanent direction parameter field under a uniform grid space reference.

[0097] The objective function for the joint inversion of gravity and magnetic normal vector magnetization is an embedded approach that directly connects the skarn compatibility weight map and the feasible domain map of remanent magnetization direction. Specifically, it involves generating spatial weight coefficients at the pixel level using the skarn compatibility weight map, and using these coefficients to weight structural constraints and parameter bounds. This allows pixels with high spatial weights to have stronger structural synergy between density and magnetic susceptibility volumes and to broaden the range of values ​​for magnetic susceptibility and remanent magnetization amplitudes. Conversely, it tightens the feasible range of magnetic susceptibility and remanent magnetization for pixels with low spatial weights and reduces the synergistic strength between density and magnetic susceptibility volumes. At the pixel level, it generates a direction allowance set using the feasible domain map of remanent magnetization direction, and writes this set into the direction constraint term. Penalties are imposed on remanent magnetization directions that deviate from the allowance set, and the remanent magnetization direction volume is projected back into the allowance set after each parameter update. This guides the vector magnetization to converge to a solution consistent with the mineralization process during the solution process, directly addressing the business pain point of difficulty in distinguishing the sources of strong magnetic anomalies of the same amplitude.

[0098] The implementation of geometric hard constraints is to limit the model support range by defining the working domain of the contact zone and the unified grid space reference. Specifically, a solution mask is set within the working domain of the contact zone, and only the pixels within the mask participate in the calculation of the objective function and parameter updates. Fixed boundary conditions are set on the model parameters at the mask boundary to keep the parameters at the boundary as background values ​​and prevent the model volume from crossing the boundary. Anisotropic structural constraints are set inside the mask along the direction of the contact zone to strengthen the continuity along the direction and suppress unreasonable extension of the normal, thereby ensuring that the inverted volume is geometrically controlled by the contact zone and does not produce out-of-domain artifacts.

[0099] Solving for the density volume, magnetic susceptibility volume, and remanent orientation volume involves a combination of joint updating and directional projection. Specifically: the density volume, magnetic susceptibility volume, and remanent orientation volume are initialized within a unified grid spatial reference. The weighted sum of gravity and magnetic field observation residuals is iteratively minimized under spatial weight adjustment, while structural constraints are applied to enhance the co-location and boundary consistency of the density and magnetic susceptibility volumes in high-weight regions. After each update, the remanent orientation volume is projected at the cell level onto the feasible remanent orientation region map within the given azimuth and dip allowable ranges based on the directional constraints. Iteration continues under the mask control of geometric hard constraints until the objective function converges, outputting a skarn-compatible weighted map and a feasible remanent orientation region. Figure 1 The density volume, magnetic susceptibility volume, and remanence direction volume are used as inputs to S5 for source discrimination.

[0100] S5. Using the density volume, magnetic susceptibility volume, and remanent direction volume output from S4 as input, and combining the skarn compatibility weight map and remanent direction feasible domain map output from S3, the same strong magnetic units are divided into skarn ore bodies or basic or ultrabasic dikes according to the classification discrimination rules. The data residual consistency and prior consistency are calculated, and the source discrimination map and confidence map are generated.

[0101] Specifically, a common-amplitude strong magnetic unit refers to a connected region within the contact zone working domain and the unified grid spatial reference, based on the common high-amplitude region of the amplitude map observed by airborne magnetic method and the magnetic field prediction amplitude map obtained by forward modeling of the magnetic susceptibility volume and the remanent direction volume, and divided according to the amplitude threshold and pixel connectivity, which serves as the spatial object for classification.

[0102] The classification and discrimination rules refer to a set of rules that take strong magnetic units of the same amplitude as objects, and output source type labels based on the unit statistics of density volume, magnetic susceptibility volume and remanent direction volume, combined with the unit statistics of skarn compatibility weight map and remanent direction feasible domain map, according to preset thresholds and logical relationships.

[0103] Data residual consistency refers to the degree of agreement between the residuals of gravity prediction and magnetic field prediction obtained by forward modeling density volume, magnetic susceptibility volume and remanent direction volume relative to the observed data within the same strong magnetic cell. It is used to verify the consistency between the classification results and the observations.

[0104] Prior consistency refers to the degree of conformity between the source type judgment of a strong magnetic unit and the skarn compatibility weight map and the feasible domain map of the remanent magnetic direction within that unit. It is used to verify the consistency between the classification results and prior constraints.

[0105] The source discrimination map is a spatial representation that uses strong magnetic units of the same amplitude as basic objects within a unified raster spatial reference, and assigns source type labels such as skarn ore bodies or basic or ultrabasic dikes.

[0106] A confidence map is a spatial representation that maps data residual consistency and prior consistency to confidence values ​​of strong magnetic cells of the same amplitude within a unified raster spatial reference.

[0107] The classification of strong magnetic units of the same amplitude into skarn ore bodies or basic or ultrabasic dikes is driven by a combination of S4 and S3 outputs. Specifically, using the density volume, magnetic susceptibility volume, and remanent direction volume from the S4 output as inputs, the unit mean and quantiles of strong magnetic units of the same amplitude are calculated within a unified grid spatial reference to obtain the conformity rate of density index, magnetic susceptibility index, and remanent direction volume and remanent direction feasible region map. Simultaneously, the average value and connectivity index of the skarn compatibility weight map from the S3 output are calculated within the unit. The classification discrimination rule is set as a logical combination, specifically, when the unit average value of the skarn compatibility weight map is greater than a first threshold and the remanent direction... If the agreement rate between the volume and the feasible domain map of the remanent magnetic direction is greater than the second threshold, and the magnetic susceptibility index reaches the third threshold, then the same-amplitude strong magnetic unit is classified as a skarn ore body. If the average value of the unit in the skarn compatibility weight map is less than the first threshold, and the density index and magnetic susceptibility index are both greater than the fourth threshold, and the agreement rate between the volume in the remanent magnetic direction and the feasible domain map of the remanent magnetic direction is less than the second threshold, then the same-amplitude strong magnetic unit is classified as a basic or ultrabasic dike. The first, second, third, and fourth thresholds are preset thresholds used to directly respond to scenarios where the source of the same-amplitude strong magnetic field is difficult to distinguish, so that the classification relies on the joint evidence of physical properties and prior knowledge rather than just the magnitude of the amplitude.

[0108] The calculation of data residual consistency and prior consistency, generating source discrimination maps and confidence maps, is a spatial product output by jointly scoring observational and prior consistency. Specifically, for each strong magnetic unit of the same amplitude, gravity and magnetic field predictions are obtained through forward modeling using density volume, magnetic susceptibility volume, and remanent direction volume, and compared with corresponding observations to obtain unit-level consistency scores for gravity and magnetic field residuals, defined as data residual consistency. Based on the given source type labels, the average value of the skarn compatibility weight map within the unit and the consistency rate between the remanent direction volume and the remanent direction feasible domain map are calculated and combined to form prior consistency. The source type labels are written into the corresponding units of the source discrimination map, and the data residual consistency and prior consistency are combined according to preset weights to form a confidence score and written into the corresponding units of the confidence map. The source discrimination map and confidence map are output with the contact zone working domain as the range and a unified raster spatial reference as the pixel unit, and are used as inputs for S6 for subsequent geometric consistency and co-directional change verification.

[0109] S6. Based on the source discrimination map, confidence map, and mineral assemblage consistency map, conduct geometric consistency and co-directional variation tests along the contact zone direction, mark continuous skarn source segments as priority zones and mark inconsistent segments as low priority segments to generate priority zones.

[0110] Specifically, the priority zone is a strip-shaped result marked along the geometric direction of the contact zone within a unified grid spatial reference. It consists of continuous skarn source sections that meet the tests of geometric consistency and unidirectional variation, and is used to indicate priority working sections.

[0111] A continuous skarn source segment refers to a continuous band-shaped portion in space formed by a set of adjacent pixels along the geometric direction of the contact zone in the source discrimination map, all of which are labeled as skarn ore bodies, and which simultaneously meets the preset threshold conditions in the confidence map and the mineral assemblage consistency map.

[0112] Inconsistent segments refer to segments in a set of pixels distributed along the geometric direction of the contact zone that do not meet the geometric consistency condition or the condition of unidirectional change.

[0113] Low-priority segments refer to the annotation results of inconsistent segments, and are used to annotate low-priority segments outside the priority band.

[0114] The geometric consistency and co-directional variation test along the contact zone's geometric orientation is a segment-level test route based on contact control and multi-source covariance evidence. Specifically: First, using the contact zone geometry as a linear benchmark, a linear index is established along the contact zone's geometric orientation within a unified raster space, and the contact zone geometry is discretized into several segments of equal length. Second, within each segment, the pixel proportion of skarn orebody tags is extracted from the source discrimination map, the segment-level average is calculated for the mineral assemblage consistency map, and the segment-level average is calculated for the confidence map. Within that segment, the angle between the principal direction of the skarn source pixels and the contact zone's geometric orientation is calculated. Third, the geometric consistency condition is defined as the principal direction. The angle between the contact zone and the geometric trend does not exceed the first threshold, and the proportion of skarn orebody labels is not less than the fifth threshold. The condition of unidirectional change is defined as the trend of the proportion of skarn orebody labels, the average value of mineral assemblage consistency at the segment level, and the average value of credibility at the segment level are consistent between two adjacent segments, and there are no contradictory abrupt changes between the three between adjacent segments. Finally, segments that simultaneously meet the geometric consistency condition and the unidirectional change condition are marked as consistent segments, and those that do not meet the condition are marked as inconsistent segments, so as to avoid misclassifying isolated or abnormal segments that deviate from the geometric control of the contact zone as continuous mineralized segments, and directly respond to the business pain point that it is difficult to stably indicate the source of strong magnetic anomalies of the same amplitude along the contact zone.

[0115] The annotation approach, which prioritizes continuous skarn source segments and denotes inconsistent segments as low-priority segments, is based on segment-level connectivity and confidence screening. The specific technical approach can be understood as follows: First, connectivity analysis is performed on the set of consistent segments according to the geometric orientation of the contact zone. The largest continuous sequence composed of adjacent consistent segments is extracted, and each largest continuous sequence is defined as a continuous skarn source segment. Second, the sequence length of each continuous skarn source segment is calculated and compared with a first length threshold, retaining only sequences whose length is not less than the first length threshold. Third, within the retained sequences, the segment-level average value of the confidence map and the segment-level average value of the mineral assemblage consistency map are verified to ensure they are not lower than the sixth and seventh thresholds, respectively; those that do not meet these thresholds are removed. Finally, the verified continuous skarn source segments are marked as priority segments, and inconsistent segments are marked as low-priority segments. The annotation results are maintained to ensure continuous representation along the geometric orientation of the contact zone, avoiding disruption of the continuity of the priority working area due to local fluctuations.

[0116] The priority zone is generated by constructing a strip-shaped spatial product based on the line annotation results. Specifically, it involves selecting a set of strip-shaped pixels within a unified raster spatial reference that is consistent with the normal range of the contact zone working domain, centered on the linear annotation of the priority zone. Each set of strip-shaped pixels corresponding to a continuous skarn source segment is assigned a priority zone identifier code, while the set of strip-shaped pixels corresponding to low-priority segments is assigned a low-priority identifier code. The priority zone is then output within the contact zone working domain.

[0117] S7. Using the priority band output by S6 as input, superimpose the source discrimination map and confidence map output by S5, and refer to the spatial relationship of the density volume, magnetic susceptibility volume and remanent direction volume output by S4, give a hole layout reference map in the working domain of the contact zone, and generate suggested drilling positions and depth ranges as the final result.

[0118] Example

[0119] S1. Construct the contact zone working domain and a unified grid space reference.

[0120] This step addresses the difficulty in distinguishing the sources of strong magnetic anomalies of the same magnitude by prioritizing the establishment of contact zone geometry, contact zone working domain, and a unified raster spatial reference. This ensures that subsequent mineral identification, causal reasoning, and joint inversion under the control of the contact zone are all conducted within the same spatial reference and pixel unit. The input to step S1 is regional geological data and surface / shallow contact zone alignment information. The output is contact zone geometry, contact zone working domain, and unified raster spatial reference, forming a mask and linear index along the line.

[0121] Define a spatial reference frame and a set of symbols. The projected coordinate reference frame is defined as the symbolic spatial reference frame, denoted as: This is used to unify the coordinate representation of all subsequent data. The contact zone geometry is represented by line objects, defined as a set of parameterized curves:

[0122]

[0123] in, For contact zone geometry, In a space reference frame The following is based on the arc length parameter The coordinates of the contact point in the positioning system. Here is the arc length parameter along the contact zone. Let be the total arc length of the contact zone. To characterize the geometric control along the direction and normal, tangential and normal descriptions are introduced: For contact strip in The tangential unit vector at that point To and Orthogonal normal unit vector.

[0124] The contact zone working domain is defined as a normal strip-shaped region surrounding the contact zone geometry, used to limit the solution range and cell set for subsequent calculations. The contact zone working domain is denoted as:

[0125]

[0126] in, For contact zone working area, For spatial reference frame The location point in the middle, Location point To contact zone geometry The shortest distance function, This refers to the normal range (half-width) of the contact strip's working domain. The typical normal width of the mineralized alteration zone and the engineering deployment width are used as the basis to ensure that the subsequent continuity test and joint inversion along the contact zone direction are strictly geometrically controlled and do not extend into invalid areas outside the domain.

[0127] Establish a unified raster space benchmark so that all outputs can be uniformly expressed in pixels.

[0128] The unified grid space reference is denoted as:

[0129]

[0130] in, To unify the grid space benchmark, The origin of the grid in the spatial reference frame coordinates in and The grid is respectively in and Pixel size in direction, and The grid is respectively in and Number of pixels in the direction.

[0131] Within a unified raster space reference, the set of pixels is denoted as:

[0132]

[0133] in, For a set of pixels, For the first Line 1 Columnar image index.

[0134] To connect the contact band working domain to the cell set, a mask function is defined:

[0135]

[0136] in, For masking functions, Represents a pixel Located in the contact zone working area Inside, This indicates that the pixel is located outside the working domain. Mask function S2 to S4 are used to limit the calculation range, and S3 to S6 are used for connectivity screening and consistency verification along the contact zone direction.

[0137] To support subsequent segment-level analysis along the contact zone direction, a linear index is established along the line, denoted as:

[0138]

[0139] in, For the linear index set of contact zone geometry, For linear index numbering, This represents the total number of index segments.

[0140] Index and contact band parameters The relationship is denoted as:

[0141]

[0142] in, For the first The position of the arc length parameter corresponding to each index. This is the index step size along the contact zone geometry. (Through...) In S6, geometric consistency and unidirectional change tests can be performed directly according to the contact zone geometry, and the strict implementation of segment-level connectivity and length thresholds can be guaranteed.

[0143] In summary, the output of step S1 is: contact zone geometry. Contact zone working area Unified raster space benchmark Pixel set Mask function With contact band linear index .

[0144] In S2, the contact strip is the working area. With a unified grid space reference Spectral interpretation and mineral assemblage consistency map generation are performed based on the input.

[0145] In S3, the contact zone geometry Set the causal direction and spatial adjacency for the input;

[0146] In S4, a mask function Implement geometric hard constraints;

[0147] In S5 and S6, the contact strip is linearly indexed. Perform unidirectional change tests and priority zone markings along the geometric direction of the contact zone to ensure that subsequent determination of the source of strong magnetic anomalies of the same amplitude is always carried out under a unified spatial reference and contact control.

[0148] S2. Construction of mineral assemblage consistency maps and airborne magnetic feature grid sets, and airborne gravity feature grid sets.

[0149] This step, constrained by the contact zone geometry, contact zone working domain, and unified raster spatial reference established in S1, involves spectral interpretation and generation of mineral assemblage consistency maps. This ensures that the differentiation of sources of similarly strong magnetic anomalies is based on mineral evidence and carbonate indicators compatible with the skarn mineralization process. Input data is in a spatial reference frame. Unified expression within the contact zone working domain mask function Pixel-level calculations are performed within a limited range, and the output is a mineral assemblage consistency map and an airborne magnetic feature raster set and an airborne gravity feature raster set for subsequent processing.

[0150] Under a unified grid space benchmark With cell set Within this range, spectral interpretation is performed to obtain mineral indicator weights. Mineral indicator indices are calculated for garnet, pyroxene, epidote, and iron oxides, and their dimensions are standardized to the interval [0,1], denoted as: For garnet, the pixel-level indicator weight, For pyroxene, the pixel-level indicator weight, As a pixel-level indicator weight for epidote, Pixel-level indicator weights for iron oxides; where each This indicates the degree to which the diagnostic features of the corresponding mineral match in that pixel. A carbonate indicator mask is constructed based on the diagnostic absorption characteristics of calcite and dolomite, denoted as: ,in This indicates that the pixel is in the carbonate indicator area. This indicates that the pixel is in a non-carbonate background. The above weights and masking are only applied to the contact zone working domain. Intrinsic participation in calculations, through a mask function Implement domain-specific restrictions.

[0151] Combining contact zone geometry With the set of linear indices along the line In the contact zone working area The consistency diagram of endogenous mineral assemblage, denoted as The calculation formula is as follows, used to simultaneously express evidence of carbonate indicator differentiation, mineral co-occurrence relationships, and connectivity along the contact zone:

[0152]

[0153] in: The pixel value of the mineral assemblage consistency map; This is the mask function (1 inside the domain, 0 outside the domain); For carbonate indicator mask; For mineral weighting coefficients, satisfying This is used to reflect the dominant co-occurrence relationship between iron oxides and garnet and pyroxene in skarn mineralization; These are the weighting coefficients for co-occurrence terms; For mineral co-occurrence terms, the mineralization compatibility is expressed by the intensity of paired co-occurrence of iron oxides with garnet / pyroxene; This is a normalization constant that normalizes the values ​​within the parentheses to [0,1]. The inhibition coefficient for non-carbonate background is used in... This reduces the consistency value; The connectivity factor is along the geometry of the contact zone; For pixels The set of neighborhoods along the contact zone geometry, with the neighborhood defined by the arc length radius. limited; This represents the number of neighboring pixels. This is an indicator function that takes the value 1 if the condition is true and 0 otherwise. This represents the threshold for the co-occurrence of mineral pairs. The nearest arc length parameter position of the pixel in the contact zone geometry; For the contact zone geometry in arc length parameters The coordinates of the location; The coordinates of the pixel center are based on a unified raster space reference. Calculated; This is the total arc length of the contact zone; The coordinates of the grid origin; Pixel size; This is a cell row and column index, and all the above calculations are performed within the contact strip working domain. Implemented internally.

[0154] The above formula uses carbonate indicator masks for differentiation, mineral-dominated co-occurrence terms to enhance compatibility expression, and connectivity factors along the contact zone geometry to screen out isolated high values. This serves the scenario of differentiating the sources of strong magnetic anomalies of the same amplitude, and provides the spatial weights and directional constraints of subsequent causal inference and joint inversion with the physical evidence basis of the mineralization process.

[0155] Airborne magnetic and airborne gravity features are extracted and used to generate feature grids, which serve as input for subsequent steps. The set of airborne magnetic feature grids is denoted as... ,in To establish a unified grid space benchmark The magnetic field characteristic grids (amplitude grids, gradient grids, boundary enhancement grids) calculated below; the set of airborne gravity characteristic grids is denoted as... ,in To establish a unified grid space benchmark The gravity feature grid (such as amplitude grid, gradient grid, boundary enhancement grid) is calculated below.

[0156] The above and All are in the contact zone working area Inner through mask function It is limited to being used as an observation anomaly input and spatial feature input when constructing the multi-stream heterogeneous graph causal inference model in step S3.

[0157] The output of step S2 includes a mineral assemblage consistency diagram. Airborne magnetic feature grid set With airborne gravity feature grid set Among them, the mineral assemblage consistency diagram In the contact zone working area With a unified grid space reference The strength of skarn mineralization compatibility is expressed in pixels, and the geometry along the contact zone is also considered. The connectivity has been incorporated into the consistency calculation; airborne magnetic features and airborne gravity features are stored in the form of feature grids for use in step S3 with contact zone geometry. We will jointly construct causal directions and spatial adjacencies, and continue to advance the causal transmission from source type to physical properties to observational anomalies.

[0158] S3. Generation of multi-flow heterogeneous graph causal reasoning model and skarn compatibility weight map and remanent magnetization direction feasible region map.

[0159] This step involves the contact zone geometry already constructed in step S1. Contact zone working area Unified raster space benchmark With cell set Under the constraints, combined with the mineral assemblage consistency map output in step S2 Airborne magnetic feature grid set With airborne gravity feature grid set A multi-flow heterogeneous graph causal reasoning model is constructed and propagated along the causal direction of source type → physical property → observation anomaly, outputting a skarn compatibility weight map and a feasible domain map of remanent magnetization direction.

[0160] Among them, the source type node uses skarn ore bodies and basic or ultrabasic dikes as candidate labels; the physical property node is described at the pixel level by parameters such as density increment, magnetic susceptibility increment, and remanence direction (azimuth and dip); the observation anomaly node uses airborne magnetic and airborne gravity characteristics as observations; spatial correlation is determined by contact zone geometry. Compared with the geometric connectivity factor along the contact zone calculated in step S2 Control, so that the inference results are within the contact zone working domain Maintain reasonable continuity along the direction.

[0161] In pixels Above, using a mineral assemblage consistency diagram The prior strength is set as the source type of skarn ore body, using an airborne gravity feature raster set. Density increment compatibility is scored using aeromagnetic feature grid sets. The compatibility between magnetic susceptibility and direction is scored under candidate remanent magnetization directions, and the geometric connectivity factor along the contact band is used to determine the compatibility. Spatial propagation was performed to obtain the pixel values ​​of the skarn compatibility weight map and the directional allowable range of the remanent magnetization direction feasible domain map.

[0162] The calculation uses the following formula (based on a unified raster space reference). Contact zone working area Internal computation, and only for mask functions (pixel execution)

[0163]

[0164] in: Pixel values ​​for the skarn compatibility weight map; The feasible region map for the remanent magnetization direction in pixels The direction allowed set on; For the mask function (derived from step S1, defining the working domain of the contact zone) ); The connectivity factor along the contact zone geometry (derived from step S2, used to suppress isolated high values ​​and enhance continuity along the direction). These are the weighting coefficients for prior knowledge, gravitational compatibility, and magnetic compatibility; To create a mineral assemblage consistency diagram A function mapped to the a priori strength of the skarn ore body source type (e.g., a monotonic mapping) ,in (This refers to the prior sensitivity coefficient). For the first Pixel values ​​of an airborne gravity feature; This is a compatibility scoring function for gravity characteristics and the density increment range of skarn ore bodies. This represents the a priori interval for the density increment of the skarn ore body. The weights of each gravitational characteristic; For the first Pixel values ​​of aeromagnetic features; Magnetic characteristics in the range and direction of magnetic susceptibility increment The compatibility scoring function is as follows: This represents the a priori interval for the magnetic susceptibility increment of the skarn ore body. The weights of each magnetic characteristic; For the azimuth candidate set (e.g., discrete in...) (angle set) For the candidate set of tilt angles (e.g., discrete in...) (angle set) These are the pixel-level azimuth and tilt angles in the direction of remanence, respectively. For directional consistency gate functions, when The value is 1 if it matches the polarity morphology, anomalous offset, and sign pattern of the dip derivative in the characteristics of airborne magnetics; otherwise, it is 0. This is the feasible threshold for the direction, used to filter the feasible range of remanent magnetization directions; Let be the normalization constant, so that Normalize to the interval [0,1].

[0165] The above calculations are based on the working area of ​​the contact strip. For the scope, within a unified raster space reference With cell set Executed pixel-by-pixel, and by the mask function. The solution domain is limited. The node types in the multi-flow heterogeneous graph causal inference model are source type nodes, property nodes, and observation anomaly nodes. The edge types are directed edges from source type nodes to property nodes and from property nodes to observation anomaly nodes, and these edges are geometrically connected along the contact zone between pixels. The direction of propagation is determined by spatial adjacency to ensure consistency; causal direction propagation is first determined by the source type node. Define the a priori strength of the skarn ore body, and then use the density increment range of the physical property nodes. Magnetic susceptibility increment range Candidate set of remanent magnetization directions The airborne gravity and airborne magnetic characteristics of observed anomaly nodes are scored for compatibility. Finally, the observation compatibility scores are summarized from bottom to top along the directed edges and then processed by the connectivity factor. Propagation yielded a skarn compatibility weight map. Feasible region diagram with remanent magnetization direction .

[0166] The output of this step includes: a skarn compatibility weight map in pixels. The pixel value serves as the spatial weight for step S4; and the feasible region map of remanent magnetization direction is plotted in pixels. Its direction allows for the set Directional constraints used in step S4.

[0167] S4. Joint inversion objective function of skarn compatibility weight and remanent magnetic direction constraint embedded gravity and magnetic vector magnetization.

[0168] Under the constraints of the contact zone geometry, contact zone working domain, and unified grid space reference established in step S1, step S4 takes the skarn compatibility weight map and remanent magnetic direction feasible domain map output in step S3 as input, directly writes the pixel-level spatial weight and directional constraints into the gravity and magnetic vector magnetization joint inversion objective function, and implements geometric hard constraints with a mask function. The density volume, magnetic susceptibility volume, and remanent magnetic direction volume are solved only in the contact zone working domain, which aims to solve the business pain point of difficulty in distinguishing the source of strong magnetic anomalies of the same amplitude.

[0169] Specifically, a unified grid space reference is adopted. Pixel set With mask function In the contact zone working area The solution is performed at the cell level; the contact zone geometry is represented by a set of parametric curves. ,in Arc length parameter The coordinates of the location Let be the total arc length, and be at position. The tangential and normal unit vectors at point A are denoted as follows: and .

[0170] The skarn compatibility weight map output by S3 is denoted as... , used as spatial weight; the feasible region graph of the remanent magnetization direction is denoted as This represents the allowable set of remanent magnetization azimuth and tilt angles at the pixel level. The unknowns to be solved include density volume. Magnetic susceptibility Azimuth angle of the remanent magnetization direction body and tilt angle .

[0171] Observational data are recorded as gravity observations under a unified grid spatial reference. With magnetic field observation The corresponding forward modeling operator is gravity forward modeling. Forward modeling with magnetic field The forward modeling of the magnetic field takes into account the geomagnetic field direction vector. As known environmental parameters.

[0172] In the structural constraints, an anisotropic term along the geometric direction of the contact zone is introduced, and mapped through the nearest arc length parameter of the contact zone geometry. The pixel center coordinates ; using spatial gradient operator Construct a tangential component by approximating the gradient between the density volume and the magnetic susceptibility volume at the pixel center. , and normal components , and with weighting coefficients and Adjust the constraint strength along the direction and normal, while using the coordination coefficient. The boundary between the density volume and the magnetic susceptibility volume is consistent.

[0173] Directional constraints are measured by set distance. Characterizes the deviation of the pixel-level orientation from the allowable set. The observation residual weighting coefficient is denoted as... and The total weight of the structural constraints is denoted as The directional constraint weight is denoted as Under geometrically hard constraints, for pixels outside the working domain of the contact zone. The density and susceptibility volumes are directly fixed as background values. and It does not participate in the solution process to prevent the model from expanding out of bounds.

[0174] Under the above unified expression, the objective function for the joint inversion of gravity and magnetic vector magnetization is calculated only within the working domain of the contact zone, and is written as:

[0175]

[0176] in: The joint inversion objective function; For the first Line 1 Column image index; This is a mask function (takes 1 only within the working domain of the contact strip, and 0 outside the domain); The pixel values ​​of the skarn compatibility weight map (importance weights due to prior source type); and These are the pixel values ​​from gravity and magnetic field observation data, respectively. For density volume Gravity predictions obtained through forward modeling; For the magnetic susceptibility body Remanent magnetization direction body In the direction of the geomagnetic field Magnetic field prediction obtained by forward modeling under constraints; Weights for observed residuals; For structural constraint weights; These are the structural constraint coefficients along the geometric direction of the contact zone and the normal direction, respectively; It is the coordination coefficient between the density volume and the magnetic susceptibility volume along the direction, used to strengthen the consistency of their boundaries; To establish a unified grid space benchmark Spatial gradient operator; Geometry from the cell center to the contact zone The location of the nearest arc length parameter; and These are the tangential and normal unit vector fields of the contact zone, respectively. For directional constraint weights; Residual magnetization direction Allowable set of feasible region graphs to the direction of remanent magnetization The set distance metric (e.g., minimum angular distance) is used to penalize deviations from the allowed set. Geometric hard constraints are implemented through... It is explicitly reflected in the objective function and fixed at the out-of-domain cells. , To prevent extraterritorial expansion.

[0177] During the solution process, minimize Obtain compatibility weights with skarn Figure 1 The spatial weighted response and feasible region of remanent magnetization direction Figure 1 The directional constraint response makes the density volume Magnetic susceptibility Within the working domain of the contact zone, maintain reasonable continuity and consistent boundaries along the geometric direction of the contact zone, and ensure the remanent magnetization direction is consistent. Converging at the pixel level to the allowable set The final output density volume Magnetic susceptibility Remanent magnetization direction body It is expressed in pixels under a unified raster space reference and is used as input to step S5 for the generation of source discrimination map and confidence map.

[0178] S5. Generation of classification and source discrimination diagrams and credibility diagrams.

[0179] The density volume output in step S4 Magnetic susceptibility Remanent magnetization direction body As input, combined with the skarn compatibility weight map output from step S3. Feasible region diagram with remanent magnetization direction The contact strip working area established in step S1 Unified raster space benchmark With cell set Under constraints, the system completes the division of strong magnetic anomalies of the same amplitude, determines the source type according to the classification rules, and calculates the consistency of data residuals and prior knowledge, ultimately generating a source discrimination map and a confidence map. To ensure that the classification serves scenarios where the source of strong magnetic anomalies of the same amplitude is difficult to distinguish, both classification and consistency quantification are directly driven by physical properties and prior maps.

[0180] In the contact zone working area Internally based on a unified grid space standard Constructing strong magnetic units of the same amplitude. Selecting a set of airborne magnetic feature grids. The characteristic grid representing the magnetic field amplitude is denoted as . ; and the magnetic field forward modeling in step S4 Obtain the amplitude of the predicted field, denoted as... ,in For the reason and The magnetic field amplitude obtained through forward modeling. Based on a preset amplitude threshold. With cell connectivity, within the contact zone working domain, the following conditions are met:

[0181] or The pixels are clustered to obtain a set of strong magnetic units with the same amplitude, denoted as . ,in For the first A set of pixels in units, The number of units is specified. The above construction ensures that the classification objects are strongly magnetically connected regions within the working domain along the contact zone.

[0182] For each strong magnetic unit of the same amplitude Calculate prior and physical property statistics and residual indices at the unit level to support classification rules and consistency scores.

[0183] set up This is the unit average of the skarn compatibility weight map; let... The unit average of the mineral assemblage consistency diagram (derived from step S2);

[0184] set up , These are the unit averages for magnetic susceptibility volume and density volume, respectively;

[0185] Let the direction coincidence rate be... ,in The indicator function is set to 1 if the remanent magnetization direction of a pixel falls within the allowable set of the feasible region map of the remanent magnetization direction, and 0 otherwise. To quantify the consistency of data residuals, cell-level gravity residuals and magnetic field residuals are defined as follows:

[0186]

[0187] in, and To unify gravity and magnetic field observations under a grid space benchmark, and For the forward modeling of gravity and magnetic field in step S4, and is the residual normalization constant.

[0188] To quantify prior consistency, a prior interval is introduced. and (Derived from the prior physical properties in step S3), and adopted , and The joint expression is used as a measure of prior conformity.

[0189] Based on the above quantization, classification rules are applied to each strong magnetic field unit of the same amplitude within the contact strip working domain, and unit-level confidence scores are generated. The formula is as follows, used to simultaneously output source type discrimination indicators and confidence scores, thereby distinguishing strong magnetic sources:

[0190]

[0191]

[0192] in: This is an indicator for identifying the source type of skarn ore bodies. A value of 1 indicates that the ore body is a skarn ore body, and a value of 0 indicates that the ore body is a basic or ultrabasic dike. Assess credibility at the unit level; The threshold for skarn compatibility weight and directional consistency rate; The threshold values ​​for gravity and magnetic field residuals; For consistency fusion weights; The sensitivity coefficient for residual consistency; is the normalization constant for prior consistency.

[0193] The first term in the formula is a logical combination of classification rules, ensuring that the rules only apply if they are a priori compatible. , , , ) and observation compatibility ( , The first term is only output when both conditions are met; the second term is a weighted composite of data residual consistency (exponential decay term) and prior consistency (three-factor joint term) to obtain the unit-level confidence level. .

[0194] The cell-level results are written into the source discrimination map and confidence map under the unified raster spatial benchmark.

[0195] Specifically, for each pixel Assignment source type label diagram skarn ore body Otherwise, assign it as a basic or ultrabasic dike; simultaneously assign a confidence level map. The above drawings are only applicable to the contact strip working area. Inner through mask function The output ensures that source identification and credibility expression are strictly controlled by contact zone geometry. Control. The generated source discrimination map. With credibility graph As input to step S6, geometric consistency and unidirectional change checks are performed on the contact strip geometry to extract the preferred strip.

[0196] S6. Perform geometric consistency and unidirectional change checks along the contact zone geometry and generate the preferred zone.

[0197] This step is based on the source discrimination map and confidence map output in step S5, and on the contact zone geometry established in step S1. Contact zone working area Unified raster space benchmark With cell set Under the constraints, geometric consistency and unidirectional variation tests were carried out along the geometric direction of the contact zone, and continuous skarn source sections were marked as priority zones, while inconsistent segments were marked as low priority segments.

[0198] Input includes source discrimination graph With credibility graph (Derived from step S5), Mineral Assemblage Consistency Diagram Connectivity factor along the contact zone geometry (Derived from step S2), and mask function (Derived from step S1, used to define the working area of ​​the contact strip) ).

[0199] To construct segment-level analysis along the contact zone geometry, a linear index set of the contact zone is used. The lower section center arc length parameter ,in For the index step size along the contact zone geometry, Let the number of segments be defined; define the first segment as... The set of strip pixels of a segment is ,in For pixels The location of the nearest arc length parameter on the contact zone geometry ensures that the segment-level calculations are strictly performed along the contact zone geometry.

[0200] To simultaneously perform geometric consistency and unidirectional change checks in a single calculation and output a segment-level consistency determination, this implementation method adopts the following formula:

[0201]

[0202] in: For the first The consistency of a segment is determined by setting a value of 1 to indicate that the segment satisfies the geometric consistency and unidirectional change tests, and setting a value of 0 to indicate that the segment is an inconsistent segment. In the first The percentage of skarn orebody labels in the source identification diagram within the segment; The segment-level average of the credibility graph; The segment-level average of the mineral assemblage consistency diagram; The segment-level average of the connectivity factor along the contact zone geometry; These represent the changes in proportion, confidence, and consistency between adjacent segments, used for testing for changes in the same direction. This is a set of geometric consistency thresholds, which are applied to source proportion, confidence, mineral assemblage consistency, and connectivity factor, respectively. It is a tolerance coefficient that changes in the same direction, used to relax minor contradictions between adjacent segments; This is an indicator function; it takes the value 1 if the condition inside the parentheses is true, and 0 otherwise. For the first The number of pixels in a segment; Pixel labels for the source discrimination map (skarn ore body or basic or ultrabasic dike); These are the pixel values ​​of the credibility map; The pixel value of the mineral assemblage consistency map; is the connectivity factor along the geometry of the contact zone.

[0203] The above formula uses four segment-level thresholds to uniformly perform geometric consistency judgment and uses the product sign of the three quantities of adjacent segments to uniformly perform the same-direction change test. It directly serves the stable screening along the geometric direction of the contact zone and avoids misclassifying isolated high values ​​or abnormal segments that are inconsistent with the direction as continuous mineralized segments.

[0204] Based on the above segment-level determination, a linear index set along the contact zone is established. To satisfy Connectivity analysis is performed on the segments to extract the longest contiguous sequence consisting of adjacent identical segments. To avoid short segments affecting the engineering usability of the preferred band, a length threshold is introduced. (Based on the number of segments), only those with a length not less than [a certain value] are retained. The continuous sequence was labeled as a continuous skarn source segment, and a priority zone identifier was assigned along the geometric strike of the contact zone; accordingly, The set of segments is marked as inconsistent segments and is marked as low priority segments in the output.

[0205] To convert segment-level annotations into strip-shaped spatial products, each segment is represented by a set of strip-shaped pixels. For coverage, in the contact strip working area Output priority band within the range.

[0206] S7. Generate borehole location and depth ranges based on the borehole layout reference map.

[0207] This step takes the priority band output from step S6 as input, and the contact band geometry established in step S1. Contact zone working area Unified raster space benchmark With cell set Under the constraints, the source discrimination map output by step S5 is superimposed. With credibility graph And refer to the density volume output in step S4 Magnetic susceptibility Remanent magnetization direction body Based on the spatial relationship, a hole layout reference map is generated within the working domain of the contact zone, and suggested drilling positions and depth ranges are output as the final result.

[0208] Based on the linear labeling of the priority zone, in the linear index set of the contact zone Below, the set of strip pixels for each segment. Based on this, pixels marked as skarn ore bodies in the source discrimination map were screened, and the confidence map was used as the quantitative basis for borehole priority. Within each priority zone, points were selected using the minimum spacing principle along the geometric direction of the contact zone: First select those that satisfy skarn ore bodies and A set of candidate pixels with a value not lower than a preset threshold is selected, and local peak points along the direction are used as candidate drilling point locations. Minimum spacing constraints are applied along the geometric direction of the contact zone to avoid overly dense hole distribution. For continuous segments with longer priority zone lengths, linear indexing of the contact zone is used. Distribute the material evenly to ensure reasonable coverage along the direction and good accessibility for construction.

[0209] Reference density volume Magnetic susceptibility Remanent magnetization direction body Based on the spatial relationships, a suggested depth range is given for each candidate drilling point.

[0210] Specifically, this refers to the set of pixels in the strip where the candidate drilling point is located. Within the local neighborhood, the high-value connected regions of the combined magnetic susceptibility volume and density volume are consistent with the boundary (derived from the along-direction cooperative constraint in step S4), and candidate drilling points are preferentially located in... and The common high-altitude connection point; subsequently, combined with the forward modeling of the magnetic field. With gravity forward model Based on the fitting results at this location, and according to the preset engineering drilling depth range and the observation response main contribution depth window, the recommended top and bottom depth intervals are determined. This interval ensures that it covers the main contributing layer of the physical property volume and aligns with the allowable set of the remanent magnetization direction volume. Maintain consistency to avoid insufficient crossing due to being too shallow or separation from the target connected body due to being too deep.

[0211] The above-mentioned recommended depth range is determined by the high-value area of ​​the prior and confidence map of the skarn ore body as the source type. This directly solves the scenario where it is difficult to distinguish the source of the same strong magnetic anomaly, and ensures that the drilling depth design is consistent with the source identification and the spatial relationship of the physical properties.

[0212] Finally, in a unified raster space benchmark Contact zone working area Internal output orifice reference map. The orifice reference map displays candidate drill point locations based on the priority zone, overlaid with a source discrimination map. With credibility graph The base map is used to mark the suggested depth range and the corresponding linear index segment number of the contact zone for each candidate drilling point. This is to facilitate on-site deployment and project verification.

[0213] The terminal results include:

[0214] Candidate drill point locations (coordinates calculated from the cell center coordinates using a unified raster space reference) (Given)

[0215] The recommended depth range (top and bottom depths) and the priority zone marking of the drilling point are used for subsequent engineering implementation.

[0216] The foregoing description of specific exemplary embodiments of the invention is for illustrative and explanatory purposes. These descriptions are not intended to limit the invention to the precise forms disclosed, and it will be apparent that many changes and variations can be made in accordance with the foregoing teachings. The exemplary embodiments were chosen and described in order to explain the specific principles of the invention and its practical application, thereby enabling those skilled in the art to implement and utilize various different exemplary embodiments of the invention, as well as various different choices and variations. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A method for identifying the source of a same-amplitude high magnetic anomaly, characterized in that, The method comprises the following steps: S1: taking the geological boundary and structural line of the contact zone between the carbonate rock and the intrusive body as the reference, fusing the common coverage of the airborne magnetic method, the airborne gravity, the airborne hyperspectrum and the satellite hyperspectrum, and constructing the contact zone working domain and the unified grid space reference; S2: based on the contact zone working domain and the unified grid space reference, performing spectral interpretation on the airborne hyperspectrum and the satellite hyperspectrum, identifying the garnet, pyroxene, epidote and iron oxide and distinguishing from the carbonate rock indication, generating a mineral combination consistency map, and extracting the airborne magnetic method features and the airborne gravity features to generate an airborne magnetic method feature grid set and an airborne gravity feature grid set; S3: according to the mineral combination consistency map, the airborne magnetic method feature grid set and the airborne gravity feature grid set, constructing a multi-flow heterogeneous graph causal reasoning model, and transmitting along the causal direction from the source type to the physical property to the observation anomaly to determine a skarn phase compatibility weight map and a remanent magnetic direction feasible region map; S4: taking the skarn phase compatibility weight map as a space weight and the remanent magnetic direction feasible region map as a direction constraint, embedding a gravity and magnetic vector magnetization joint inversion target function to determine a density body, a magnetic susceptibility body and a remanent magnetic direction body; S5: according to the density body, the magnetic susceptibility body and the remanent magnetic direction body, the skarn phase compatibility weight map and the remanent magnetic direction feasible region map, dividing the same amplitude strong magnetic unit into a skarn ore body or a basic or ultrabasic rock vein through a classification discrimination rule, calculating data residual consistency and prior consistency, and generating a source discrimination map and a reliability map; S6: according to the source discrimination map and the reliability map, the mineral combination consistency map, performing geometric consistency and same direction change test along the contact zone direction to generate a priority zone; S7: superimposing the priority zone on the source discrimination map and the reliability map, giving a hole distribution reference map in the contact zone working domain based on the spatial relationship of the density body, the magnetic susceptibility body and the remanent magnetic direction body, and generating a recommended drilling position and depth interval.

2. The method of claim 1, wherein the same-amplitude strong magnetic anomaly source discrimination method is characterized by, The S1 takes the geological boundary and structural line of the contact zone between the carbonate rock and the intrusive body as the reference, fuses the common coverage of the airborne magnetic method, the airborne gravity, the airborne hyperspectrum and the satellite hyperspectrum, and constructs the contact zone working domain and the unified grid space reference, comprising: extracting the contact zone boundary line and the structural line from the geological map and the earth explanation result and unifying them to the same coordinate system; setting a buffer zone along the normal direction centering on the contact zone boundary line to form an initial contact zone working domain; performing a spatial intersection operation on the common coverage and the initial contact zone working domain to obtain the contact zone working domain; establishing the unified grid space reference in the contact zone working domain.

3. The method of claim 1, wherein the same-amplitude strong magnetic anomaly source discrimination method is characterized by, In the S2, the spectral interpretation is performed on the airborne hyperspectrum and the satellite hyperspectrum, the garnet, the pyroxene, the epidote and the iron oxide are identified and distinguished from the carbonate rock indication, and the mineral combination consistency map is generated, comprising: performing diagnostic absorption feature identification and mineral indication index calculation on the reflectivity sequence of the airborne hyperspectrum and the satellite hyperspectrum to output the mineral indication weight in units of pixels; A diagnostic waveband combination is selected for each of garnet, pyroxene, epidote and iron oxide, a mineral indicator index is calculated and a corresponding mineral indicator map is output; Each mineral indicator map is converted into a pixel-level indicator weight according to a unified dimension; A carbonate rock indicator mask is constructed and each pixel-level indicator weight is labeled inside and outside the carbonate rock indicator mask to obtain a carbonate rock indicator distinguished mineral indicator result; According to the carbonate rock indicator distinguished mineral indicator result, a mineral combination consistency map is determined.

4. The method of claim 1, wherein the same-amplitude strong magnetic anomaly source discrimination method is characterized by, S3 constructs a multi-flow heterogeneous graph causal reasoning model according to the mineral combination consistency map, the set of airborne magnetic feature grids and the set of airborne gravity feature grids, and transmits in the causal direction from source type to physical property to observation anomaly to determine a skarn compatible weight map and a remanent magnetization direction feasible region map, including: A source type node, a physical property node and an observation anomaly node and a contact zone geometry node are established for each pixel in the unified grid space reference; A prior probability of the source type being skarn is set on the source type node according to the mineral combination consistency map and the contact zone geometry, and prior intervals of the density increment, the magnetic susceptibility increment and the remanent magnetization strength compatible with skarn are set on the physical property node; the observation likelihood is calculated on the observation anomaly node according to the airborne magnetic feature and the airborne gravity feature, and the posterior of the physical property node and the source type node is updated from bottom to top along the directed edge; the consistency propagation is performed between pixels through the spatial adjacency edge controlled by the contact zone geometry node to obtain the posterior distribution of the source type and the physical property posterior related to the remanent magnetization direction; The posterior probability of the source type being skarn is taken as a pixel-level weight, and the pixel-level weight is output in the unified grid space reference to obtain the skarn compatible weight map; On each pixel, the remanent magnetization direction feasible region map is determined according to the airborne magnetic feature and the source type and physical property posteriors.

5. The method of claim 1, wherein the same-amplitude strong magnetic anomaly source discrimination method is characterized by, The gravity and magnetic vector magnetization joint inversion objective function includes a gravity observation residual term, a magnetic field observation residual term, a structure constraint term and a direction constraint term; S4 embeds the skarn compatible weight map as a spatial weight and the remanent magnetization direction feasible region map as a direction constraint into the gravity and magnetic vector magnetization joint inversion objective function to determine a density body, a magnetic susceptibility body and a remanent magnetization direction body, including: A spatial weight coefficient is generated at the pixel level according to the skarn compatible weight map, and the spatial weight coefficient is used for the weighted structure constraint term and the parameter boundary constraint; A direction allowed set is generated at the pixel level according to the remanent magnetization direction feasible region map, and the direction allowed set is written into the direction constraint term; The density body, the magnetic susceptibility body and the remanent magnetization direction body are initialized in the unified grid space reference, the weighted sum of the gravity observation residual term and the magnetic field observation residual term under the adjustment of the spatial weight is iteratively minimized, and a structural constraint term is applied; after each update, the remanent magnetization direction body is projected to the given azimuth and inclination allowable range of the remanent magnetization direction feasible region at the pixel level according to the direction constraint term; the iteration is continued under the control of the mask of the geometric hard constraint until the gravity-magnetic vector magnetization joint inversion objective function converges, and the density body, the magnetic susceptibility body and the remanent magnetization direction body consistent with the skarn compatible weight map and the remanent magnetization direction feasible region map are output.

6. The method of claim 1, wherein the same-amplitude strong magnetic anomaly source discrimination method is characterized by, The S5 divides the same amplitude strong magnetic unit into a skarn ore body or a basic or ultrabasic rock vein according to the density body, the magnetic susceptibility body and the remanent magnetization direction body, the skarn compatible weight map and the remanent magnetization direction feasible region map, and comprises the following steps. According to the density body, the magnetic susceptibility body and the remanent magnetization direction body, the statistical quantities such as the cell mean value of the same amplitude strong magnetic unit are calculated in the unified grid space reference, and the density index and the magnetic susceptibility index and the coincidence rate of the remanent magnetization direction body and the remanent magnetization direction feasible region map are obtained. According to the skarn compatible weight map, the cell average value in the unit is obtained. According to the density index and the magnetic susceptibility index and the coincidence rate of the remanent magnetization direction body and the remanent magnetization direction feasible region map, the cell average value, the classification discriminant rule is determined. According to the classification discriminant rule, the same amplitude strong magnetic unit is divided into a skarn ore body or a basic or ultrabasic rock vein.

7. The method of claim 6, wherein the step of determining the source of the same-amplitude high magnetic anomaly is characterized by, The S5 divides the same amplitude strong magnetic unit into a skarn ore body or a basic or ultrabasic rock vein according to the density body, the magnetic susceptibility body and the remanent magnetization direction body, the skarn compatible weight map and the remanent magnetization direction feasible region map, and comprises the following steps. When the cell average value of the skarn compatible weight map is greater than a first threshold value, the coincidence rate of the remanent magnetization direction body and the remanent magnetization direction feasible region map is greater than a second threshold value, and the magnetic susceptibility index reaches a third threshold value, the same amplitude strong magnetic unit is judged to be the skarn ore body; When the cell average value of the skarn compatible weight map is less than the first threshold value, the density index and the magnetic susceptibility index are greater than a fourth threshold value at the same time, and the coincidence rate of the remanent magnetization direction body to the remanent magnetization direction feasible region map is lower than the second threshold value, the same amplitude strong magnetic unit is judged to be the basic or ultrabasic rock vein.

8. The method of claim 1, wherein the same amplitude magnetic anomaly source discrimination method is characterized by, The S5 calculates the data residual consistency and the prior consistency, generates a source discrimination map and a reliability map, and comprises the following steps. For each same amplitude strong magnetic unit, the gravity prediction and the magnetic field prediction are obtained by using the density body, the magnetic susceptibility body and the remanent magnetization direction body, and are compared with the corresponding observation to obtain the cell level coincidence score of the gravity residual and the magnetic field residual, which is defined as the data residual consistency; According to the source type label, the average value of the skarn compatible weight map in the corresponding unit and the coincidence rate of the remanent magnetization direction body and the remanent magnetization direction feasible region map are calculated, which are combined as the prior consistency; The source type label is written into the corresponding unit of the source discrimination map; The data residual consistency and the prior consistency are combined into a reliability score according to a preset weight and are written into the corresponding unit of the reliability map.

9. The method of claim 1, wherein the same amplitude magnetic anomaly source discrimination method is characterized by, The S6 is consistent with the source discrimination map, the reliability map, and the mineral assemblage consistency map in the contact zone direction. The geometric consistency and the same direction change test include: Taking the contact zone geometry as a linear reference, a linear index is established along the contact zone geometry in the unified grid space reference, and the contact zone geometry is discretized into several equal-length segments; In each segment, the proportion of skarn ore body tags in the source discrimination map is extracted, the segment-level average value of the mineral assemblage consistency map is calculated, the segment-level average value of the reliability map is calculated, and the angle between the main direction of the skarn source pixels and the contact zone geometry is calculated; If the angle between the main direction and the contact zone geometry is not more than the first angle threshold value, and the proportion of skarn ore body tags is not less than the fifth threshold value, it is determined as the geometric consistency condition; If the proportion of skarn ore body tags, the change trend of the mineral assemblage consistency segment-level average value and the reliability segment-level average value of the adjacent two segments are consistent, and the three do not have mutual contradictory mutations between the adjacent two segments, it is determined as the same direction change condition.

10. The method of claim 9, wherein the step of determining the source of the same-amplitude high magnetic anomaly is performed by the steps of: determining a first magnetic anomaly source; and determining a second magnetic anomaly source. The S6 generates a priority zone, including: Labeling the segments that meet the geometric consistency condition and the same direction change condition as consistent segments; Performing connectivity analysis on the consistent segment set according to the contact zone geometry, extracting the maximum continuous sequence composed of adjacent consistent segments, and defining each maximum continuous sequence as a continuous skarn source segment; Calculating the sequence length of each continuous skarn source segment and comparing it with the first length threshold value, and only retaining the sequences with a length not less than the first length threshold value; In the retained sequence, whether the segment-level average value of the reliability map and the segment-level average value of the mineral assemblage consistency map are not less than the sixth threshold value and the seventh threshold value, respectively, is reviewed. The continuous skarn source segment that passes the review is labeled as a segment-level priority zone; Taking the linear label of the segment-level priority zone as the center, selecting a strip pixel set consistent with the normal range of the contact zone working domain as the spatial coverage of the priority zone in the unified grid space reference, assigning a priority zone identification code to each continuous skarn source segment corresponding to the strip pixel set, and outputting the priority zone in the contact zone working domain range.

Images