For computational convenience, assume that remanent magnetization does not exist and that all magnetization is purely vertical. In this case, the magnetic forward modeling problem can be expressed as Equation 1 (Green, 1996): ∆ T = μ 0 κ H 4 π ζ − z arctan − x − ξ y − η r z − ζ | ξ 1 ξ 2 | η 1 η 2 | ζ 1 ζ 2 . (1) here, μ 0 is the magnetic permeability of free space, κ is the magnetic susceptibility, and H is the magnetization intensity of the magnetic medium. ξ 1 , ξ 2 are the X-axis coordinates of the underground source, and η 1 , η 2 are the Y-axis coordinates, ζ 1 , ζ 2 ​ are the Z-axis coordinates. x, y, z represent the coordinates of the observation point. In magnetic anomaly boundary extraction tasks, traditional methods (such as TILT and THG) typically rely on gradients, derivatives, or second-order edge enhancement functions based on the magnetic anomaly field to indirectly extract edge information. The results are often presented as continuous response images or pseudo-color intensity maps. These methods amplify abrupt changes in field intensity to indicate the probable locations of anomaly boundaries. Generally, traditional magnetic data boundary detection methods can be expressed as Equation 2: M B = arg ⁡ max Ω F ∆ T * K . (2) here, M B represents the boundary of the underground magnetic body, F (.) is the preprocessing mapping operator, K is the filter or differential operator, and Ω is the parameter space.

Unlike the aforementioned methods, DL methods adopt a supervised learning strategy, training neural networks with large amounts of ∆ T and M B data to learn the nonlinear mapping relationship between them. This gradually builds a mapping model from the original magnetic anomaly field to the target boundary, thereby enabling the direct prediction of M B . This mapping relationship can be expressed as Equation 3 (Yu and Ma, 2021): M B = Net ∆ T , θ , (3) where Net represents the neural network, M B ^ is the predicted magnetic body boundary by the neural network, and θ denotes the network parameters.

The network architecture proposed in this paper (TCNet) is shown in Figure 1a. Assuming the magnetic anomaly data ∆ T has a size of 1 × 32 × 32, where 1 is the number of channels. TCNet consists of 9 stages, each of which is a combination of a TF Block and a conv block. The specific structure of the TF Block is shown in Figure 1b. The TF Block first divides the raw ∆ T data into P × P patches and serializes them for subsequent processing. In this paper, P = 4. Then, each patch is transformed into a feature vector (token) with a dimension of d through linear projection, constructing a serialized feature representation with spatial correlations. Simultaneously, positional encoding is applied to embed the position information of each patch into the corresponding token. For the token at the t-th position in the sequence, the position vector is defined as Equation 4: p t i = sin t 10000 i − 1 d   i f   i % 2 = 0   e l s e   p t i = cos t 10000 i d , (4) here i represents the index of each magnetic anomaly point in the t-th token, and d refers to the dimension of the ∆ T after serialization. After position encoding, the magnetic anomaly data is sent into the Transformer encoder, where it undergoes layer normalization via the Norm layer, followed by multi-head self-attention computation. In the self-attention mechanism, each token is multiplied by the corresponding matrix to obtain q, k, and v. The self-attention mechanism can be expressed as Equation 5: α = softmax Q K T d k Attention Q , K , V = softmax Q K T d k V , (5) here, all the q, k, and v form the matrices Q, K, and V, respectively. All the α i , j form the matrix α, and N represents the dimension of the sequence. The output of the self-attention mechanism is the weighted sum of all Attention Q , K , V , and this global perception capability enables the model to overcome the local receptive field limitation of traditional convolutional operations when processing magnetic anomaly sequences. It simultaneously analyzes the spatial correlation characteristics between different underground source bodies, thereby constructing a global strong dependency relationship across regions, which provides ample spatial information support for the subsequent boundary reconstruction. To enhance the information extraction capability, the multi-head self-attention mechanism increases the number of Q, K, and V. Generally, Multi-head Self-attention can be expressed as Equation 6: Multihead Q , K , V = Contact h e a d 1 , h e a d 2 , ⋯ ⋯ , h e a d h W h e a d h = A t t e n t i o n Q W h Q , Q W h K , Q W h V , (6) here, W h Q , W h K , W h V are the weights learned by the linear projection layers. The output of the multi-head self-attention mechanism is combined with the input sequence through a residual connection, forming a residual compensation structure to alleviate the gradient vanishing problem. Layer normalization is then applied to ensure the stability of the feature distribution. Afterward, the output with the residual connection is passed to the MLP, which uses the ReLU activation function to connect two fully connected layers. High-order interaction modeling is achieved by first increasing and then decreasing the feature dimensions. Similarly, the output of the MLP is connected to the input sequence via a residual connection, followed by layer normalization. The mathematical process of the Transformer module can be expressed as Equation 7: z 0 = x p 1 E ; x p 2 E ; ⋯ ⋯ ; x p n E + E P O S z l ′ = M S A LN z l − 1 + z l − 1   l = 1 ⋯ L z l ′ = M S A LN z l − 1 + z l − 1   l = 1 ⋯ L   . (7) here, E P O S is the position vector, MSA stands for Multi-head Self-Attention, and LN represents layer normalization. After passing through the TF Block, the data reaches the Conv Block, where the local receptive field characteristics of the convolutional kernel focus on the detailed features between neighboring anomalies, enhancing the model’s ability to extract local detail information. The formula is as Equation 8: F = σ B N C o n v F ′ . (8) here, F ′ is the data processed by the TF Block, BN stands for BatchNorm, and σ is the ReLU activation function. Each stage repeats the above operation.

The cascaded design of self-attention and convolution in TCNet enables the model to leverage the Transformer for modeling the complex non-local interactions between source bodies, allowing the network to focus on the most relevant magnetic anomaly areas, thereby improving boundary identification accuracy when processing large-scale geological structures. At the same time, the convolutional layers facilitate local detail extraction, enabling the network to better reconstruct subtle boundaries when handling locally complex structures.

During the network training process, the parameters are continuously updated through backpropagation to solve the optimization problem shown in Equation 9. θ ^ = arg ⁡ min θ 1 N ∑ i = 1 N L M B i , N e t ∆ T i , θ . (9) where L represents the loss function, N is the number of samples in the training set, and M B i , ∆ T i denotes the i-th pair of samples in the training set. In magnetic data boundary detection, the loss function is used to measure the difference between the reconstructed source body boundary and the true underground source body boundary, guiding the network training. The most commonly used loss function is the MSE (Mean Squared Error) loss function, which is expressed as Equation 10: L M S E = 1 N ∑ i = 1 N ∥ M B i − M B ̂ i ∥ 2 , (10)

As can be seen from the above equation, the MSE tends to smooth the prediction results to reduce the overall loss, which may lead to inaccurate boundary reconstruction by the model. To address this issue, this paper introduces an edge-enhancing loss function that explicitly constrains the differences between the predicted values of boundary pixels and their neighboring background pixels, thereby improving the model’s sensitivity to boundary structures and enhancing the clarity of the network’s boundary predictions. formula is as Equation 11: L c o n t r a s t = ∑ max 0 , 0.5 − M B ̂ n − M B ̂ b , (11) here, M B ̂ n represents the network’s predicted background value, and M B ̂ b represents the network’s predicted boundary value. As shown in the formula, this loss function constructs a lower bound constraint (set to 0.5) on the difference between the values of boundary regions and their adjacent background regions. When the predicted difference is less than this threshold, the loss term is positive and generates gradients to drive the model to enlarge the response difference between the boundary and the background. When the difference exceeds 0.5, the loss automatically becomes zero, preventing excessive adjustments to regions that are already well distinguished. This mechanism not only enhances the discriminability of boundary regions but also effectively suppresses the tendency of boundary predictions to average out, thereby improving the clarity of the boundaries. Therefore, the final loss function is defined as Equation 12: L t o t a l = L M S E + α L c o n t r a s t (12)

α is the weighting coefficient, which needs to be determined through multiple experiments. For the optimal selection of α, this paper empirically tested nine values (0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.1, and 0.3) by training the network with each α. After extensive experimentation, it was found that the network’s boundary prediction results were best when α = 0.03.

In this paper, the underground space is divided into 32 × 32 × 16 cubes, with each cube sized 1.00 × 1.00 × 1.00 km. A heterogeneous dataset containing 30,000 samples is constructed using a random walk method and forward modeling, as shown in Figure 2. The model’s magnetization is set to 0.2 SI, and the background field is set to 0. The random walk method first divides the 3D space into a four-quadrant symmetrical structure. Starting points are randomly selected in each quadrant, and their movement is controlled to take s steps (where s∈[40,70]) in a random direction, dynamically generating irregular 3D geological models. This results in a highly random dataset that can approximate any underground source, forming the basis for training a well-optimized network. All models are set with a uniform magnetization of 0.2 SI and background field interference is eliminated. A regular dataset, which does not exist in the 2,000 training sets, is created as the test set and input into the network at a ratio of 14:1:1. During the network training phase, mini-batch gradient descent (batchsize = 32) combined with the Adam optimizer is used. The initial learning rate is set 2 × 10⁻⁴, with a dynamic decay mechanism implemented using a step decay strategy, where the learning rate is multiplied by 0.8 every 15 epochs. To prevent overfitting, a Dropout regularization layer with a probability of 0.2 is introduced, and the parameter space convergence is achieved after 100 training epochs. The loss function during the network training process is shown in Figure 3. The PC is configured with a 1 × 3.90 GHz Intel Xeon W-2245 processor, a 1 × 8 GB NVIDIA Quadro RTX 4000 GPU, and 192 GB of memory. The complete training of the network takes approximately 70 min.

Considering that this study focuses on deep learning-based underground anomaly boundary detection, the Intersection over Union (IOU) metric is used to more intuitively measure the overall accuracy of boundary region localization. As shown as Equation 13: IOU = M B ̂ ∩ M B M B ̂ ∪ M B . (13)

The intersection ratio helps us intuitively assess the quality of the boundary prediction results. Additionally, to further assess the pixel-wise differences between the model’s predicted maps and the ground truth label maps, this paper introduces Mean Squared Error (MSE) as an auxiliary evaluation metric. Its calculation formula is as Equation 14: MSE = 1 N ∑ i = 1 N M B ̂ − M B 2 (14)

This paper compares the boundary detection results of the proposed method, Swin-Transformer, U-net, THG, and TILT, which are presented in the following sections.

Model 1 consists of four horizontal prismatic bodies with a magnetization of 0.2 SI, all buried at a depth of 5 km. As shown in Figure 4a, the resulting magnetic anomaly is shown in Figure 4b, and its boundary in the horizontal projection plane is displayed in Figure 4c. Figures 4d–h show the boundary detection results for each method. It can be observed that the boundaries predicted by the deep learning methods are more focused and better reflect the true boundary positions compared to traditional methods. Moreover, the proposed method performs better than the Swin-Transformer and U-Net method, with boundaries that are clearer and more closely aligned with the true model boundaries. The IOU values and MSE errors are shown in Table 1.

Model 2 consists of two step-like models of the same size and depth, but with different horizontal positions, as shown in Figure 5a. The magnetic anomaly it produces is shown in Figure 5b. Its projection in the horizontal plane is shown in Figure 5c. Figures 6d–h show the boundary detection results for three methods. From the results, it can be seen that traditional methods still fail to intuitively display the boundary positions of the anomalous bodies. The Swin-Transformer method performs well in extracting the overall shape of stepped anomaly boundaries, particularly excelling in maintaining boundary continuity. However, since its feature extraction process focuses more on global structure modeling and lacks a mechanism for enhancing local fine-grained edge features, some boundary regions exhibit a certain degree of blurriness. The U-Net method accurately locates the boundary of the step-like anomalous bodies at shallow depths, but for deeper regions, the boundary detection results show discontinuity and false anomalies. The proposed method provides more continuous and clearer boundary detection results, with higher identification accuracy.

Model 3 consists of two prisms of different sizes and depths but with the same horizontal center position, as shown in Figure 6a. The magnetic anomaly it produces is shown in Figure 6b. Its projection in the horizontal plane is shown in Figure 6c. Figures 5d–h show the boundary detection results for three methods. From the figures, it is clear that traditional boundary detection methods still diverge and cannot accurately locate the specific boundary positions. In this model, the Swin-Transformer demonstrates strong capability in modeling the overall shape of large-scale anomaly boundaries, with good boundary continuity. However, it shows signs of boundary blurring in small-scale anomaly regions, indicating a limited ability to capture features in local, small areas. The U-Net boundary detection results are clearer, with better performance in detecting the edges of shallow anomalous bodies, but the detection results for deeper anomalous bodies are worse. The possible reason is that the network extracts less global information and fails to effectively model the complex relationships between large-area anomalies. The proposed method achieves the best detection accuracy, especially in accurately characterizing the horizontal boundaries of large, deep anomalies.

To evaluate the robustness of the proposed method, 5% Gaussian white noise was added to the test model, with a signal-to-noise ratio (SNR) of approximately 22 dB. To better illustrate the impact of noise on traditional methods, we also compared with the Theat method. Figure 7 shows the magnetic anomaly and boundary detection results of the test model with added random noise. From the figure, it can be seen that the boundary detection results of the proposed method still match the true boundaries well and maintain high accuracy, indicating strong noise resistance. In contrast, the boundary detection results of traditional methods are more blurred, demonstrating their sensitivity to noise.

To further demonstrate the role of each component of the proposed method, we performed boundary identification on the above models using TCNet constrained by L_MSE and the composite loss function L_total proposed in this paper, and evaluated the results using IOU. The results are shown in Table 2. It can be seen that L_total achieves higher boundary detection accuracy on all three models, indicating that it can guide the network to more effectively extract the boundary features of underground anomalies, thereby improving the boundary detection accuracy.

The Yushishan area is located in the northwestern part of China, in Akesai County, Gansu Province. It is situated at the confluence of the Altyn Mountains fault, the Qilian orogenic belt, and the Quanjie block, on the northern edge of the Qaidam Basin (Yu et al., 2015; Liu et al., 2007; Zhang et al., 2015; Wu et al., 2001; Xu et al., 2006), as shown in Figure 8. This region is an important structural window for understanding the tectonic evolution of northeastern Tibet and southern Tibet (Wang et al., 2018; Jiang et al., 2020). The area has undergone multi-stage tectonic, magmatic, and metamorphic processes, forming unique rock assemblages and metallogenic systems (Jiang et al., 2022). Tectonically, Yushishan is located in the western part of the Central Qilian block, with the Altyn northern block to the north, the Quanjie block to the south, the Hongliugou-Lapazhen mixed rock belt to the west, and the Central-Southern Qilian arc-basin system to the east (Yu et al., 2012; Jia et al., 2016; Yang et al., 2014). It preserves a complex geological record from the breakup of the Rodinia supercontinent in the late Neoproterozoic to the subduction and collision of the Paleo-Tethys Ocean in the Early Paleozoic (Xu et al., 1999). The left-lateral Altyn Fault, as a major transverse structure, dominates the regional tectonic pattern, facilitates the displacement of crustal blocks, and promotes the ascent of deep magmas/fluids, significantly affecting the mineralization channels (Wu et al., 2001).

The regional basement is composed of the Paleoproterozoic Daken Daban Formation (Pt₁D), primarily consisting of amphibolite facies metamorphic rocks, including gneiss, amphibolite, and marble. The overlying strata, previously attributed to the Mesoproterozoic Aoyougou Formation (Ch_a), have now been redefined as the Neoproterozoic Yushishan Formation (Pt₃ys), based on zircon U-Pb ages ranging from 790 to 843 Ma (Yang et al., 2012; Liu et al., 2022). This revision resolves a long-standing stratigraphic dispute and clarifies that the Yushishan Formation is composed of interlayered light rocks, marble, and amphibolite, with distinct layered features. It contrasts sharply with the shallow marine clastic-volcanic sedimentary sequence of the regional Duoruoer Formation (Zhou et al., 2022).

The total area of the study region is 1,590 square kilometers. The survey was conducted at a 1:50,000 scale following standard measurement procedures. The magnetic anomaly data is shown in Figure 9. To meet the input requirements of the network, the data dimensions were transformed to 32 × 32 using mean padding and interpolation. The processed data were then fed into the network trained on the dataset established in Section 2.4. After prediction, the output was cropped to match the size of the study area, as shown in Figure 10a. Figure 10b presents the boundary prediction results using TILT. It can be seen that the boundaries obtained using the proposed method are clearer than those produced by traditional methods and can reflect boundary details of some complex geological structures. The method clearly divides the study area into three regions (Regions I, II, and III, shown in Figure 10a). Region I corresponds to a narrow high-magnetic anomaly zone, displaying an elongated, bead-like or strip-like distribution with good continuity. The two boundary faults in this region correlate well with the known east-west major faults (F1 and F2 in Figure 11). Region II corresponds to the northern RTP anomaly. The western section shows a moderate high-magnetic anomaly, attributed to weakly magnetic Quaternary sediments and the limestone/marble of the Annanba Formation in the western Jixi region. The stable high-value anomaly here may originate from the strata in the Jixi region. The eastern section presents a narrow high-magnetic belt trending northeast, with discontinuous bead-like anomalies along the rock contact zone. The exposure of the Annanba Formation in the Jixi region and the intrusion of diorite/monocrystalline granite from the Permian suggest that these anomalies are related to the intrusion. Region III is located in the southern part of the study area, showing large-scale negative anomalies containing two anomaly groups. The high-value anomalies are speculated to originate from intrusions within the Paleoproterozoic Daken Daban Formation. Furthermore, the proposed method successfully identified two major faults (F1 and F2) and multiple minor faults (F3–F20 in Figure 11), showing strong spatial consistency. To further validate the effectiveness of the proposed method in boundary identification under real geological conditions, the IOU metric introduced in Section 3.1 was used to perform spatial overlap analysis between the model’s predictions and existing geological structure data (sourced from the 1:50,000 geological map and Yu et al., 2015). The results show an IOU of 0.85, indicating that the proposed method demonstrates good spatial alignment with known geological boundary locations. Combining the geological data, high-magnetic anomalies and their gradient areas, as well as faults and host rock alterations (potassic, sodic, and sericitic), provide a basis for delineating potential ore targets.