Gravity anomalies are closely related to the density and spatial distribution of the earth’s internal matter. Heck et al. (2007) proposed that in the spherical coordinate system the gravity can be expressed as: g = G ρ ∭ v r ′ 2 r − r ′ ⁡ cos ⁡ ψ cos ⁡ φ ′ d r ′ d φ ′ d λ ′ l 3 (1) where l is the Euclidean distance between the observation point P (φ, λ, r) and the source point Q (φ′, λ′, r′); and ψ is the angle between the position vectors of P and Q; G is the gravitational constant, ρ is the density, φ is the latitude, λ is the longitude, and r is the radial distance. l = r ′ 2 + r 2 − 2 r ′ r ⁡ cos ⁡ ψ (2) cos ⁡ ψ = ⁡ sin ⁡ φ ⁡ sin ⁡ φ ′ + ⁡ cos ⁡ φ ⁡ cos ⁡ φ ′ ⁡ cos λ ′ − λ (3) For u = r ′ 2 r − r ′ ⁡ cos ⁡ ψ cos ⁡ φ ′ l 3 (4) we have g = G ρ ∭ v u d r ′ d φ ′ d λ ′ (5)

Density ρ can be expressed as a nonlinear function of variables such as rock mineral composition m, content p, burial depth h, rock formation time t, etc. ρ = v m , p , h , t (6)

Inserting Eq. 6 into Eq. 5 yields g = G ⋅ v m , p , h , t ⋅ ∭ v u d r ′ d φ ′ d λ ′ (7)

The above integral equation can be abstracted as the following nonlinear functional: v = f u , g (8) where u is related to the spatial location, g is the gravity anomaly, and v is related to the rock properties. This nonlinear functional defines a nonlinear mapping relationship from the spatial location and gravity anomaly to lithology.

According to the universal approximation theorem (Cybenko, 1989; Hornik et al., 1989), a feedforward neural network with a sufficient number of hidden units and a nonlinear activation function can approximate any Borel function from one finite-dimensional discrete space to another with arbitrary accuracy. In other words, a deep convolutional neural network can be used to approximate the mapping defined in Eq. 8.

In general, a convolutional neural network of depth n can be represented as: y = f n f n − 1 ⋯ f 2 f 1 x ; θ = f x ; θ (9)

The data sample x (gravity anomaly) is fed into a cascading n-layer nonlinear transform network to obtain the desired output y (lithologic distribution). The parameter θ is the learning parameter of this nonlinear transformation. We use an optimizing the algorithm to find θ so that the neural network can maximize the approximation of the mapping defined in Eq. 8.

The parametric model y = f (x; θ) defines a Conditional probability distribution p (y | x; θ). We use the principle of maximum likelihood to estimate it. The maximum likelihood estimator for θ is then defined as θ M L = argmax θ p mod ⁡ e l Y | X ; θ = argmax θ ∏ i = 1 m p mod ⁡ e l y i | x i ; θ (10) where X = x 1 , ⋯ , x m is a set of m examples. The p_model (y | x; θ) is a parametric family of probability distributions indexed by θ and it maps any configuration x to a real number estimating the true probability p_data (y | x).

The product of many probabilities is prone to numerical underflow, which is not easy to calculate. We observe that the logarithm of the likelihood does not change its arg max, but conveniently transform the product into a summation: θ M L = argmax θ ∑ i = 1 m log ⁡ p mod ⁡ e l y i | x i ; θ (11)

Dividing by m, we obtain the expectation with respect to the empirical distribution define by the training data as the estimation criterion. θ M L = argmax θ E x ∼ p ^ d a t a ⁡ log ⁡ p mod ⁡ e l y | x ; θ (12)

Deep neural network learning is estimating the parameter θ using the principle of maximum likelihood. The essence of this optimization problem is to maximize the log-likelihood, that is, to minimize the negative log-likelihood, which is equivalent to minimizing the cross entropy between the empirical distribution defined by the training set and probability distribution defined by model (Goodfellow et al., 2016).

The cost function is given by J θ = − E x , y ∼ p ^ d a t a ⁡ log ⁡ p mod ⁡ e l y | x ; θ (13)

In order to enhance the generalization ability of the neural network and avoid overfitting during the optimization, we add a parameter regularization term to the cost function to obtain a new objective function: J ∼ θ ; x , y = − E x , y ∼ p ^ d a t a ⁡ log ⁡ p mod ⁡ e l y | x ; θ + α Ω θ (14)

The problem of minimizing the objective function is a nonconvex optimization problem. This means that we cannot accurately obtain the global optimal solution of the problem. Therefore, deep neural network training uses an iterative, gradient-based optimization method to obtain a local optimal solution that makes the objective function sufficiently small. The stochastic gradient descent (SGD) algorithm is employed to solve the above nonconvex optimization problem.

We decompose the cross-entropy cost function as a sum over training examples of some per-example loss function. J θ = − E x , y ∼ p ^ d a t a ⁡ log ⁡ p mod ⁡ e l x i , y i , θ = − 1 m ∑ i = 1 m log ⁡ p mod ⁡ e l x i , y i , θ (15)

For these additive cost function, the gradient of the cross-entropy cost function is: ∇ θ J θ = − 1 m ∑ i = 1 m ∇ θ ⁡ log ⁡ p mod ⁡ e l x i , y i , θ (16)

The above gradient is an expectation that we can approximately estimate using a small set of samples. On each step of the SGD algorithm, we can sample a minibatch of examples B = x 1 , ⋯ , x m ′ drawn randomly from the training set. The estimate of the gradient is formed as g r a d i e n t = − 1 m ′ ∇ θ ∑ i = 1 m ′ log ⁡ p mod ⁡ e l x i , y i , θ (17)

The stochastic gradient descent algorithm then follows the estimated gradient downhill: θ ← θ − ε ⋅ g r a d i e n t (18) where ε is the learning rate.

The stochastic gradient descent algorithm is sometimes very slow or unreliable in the learning process. The method of stochastic gradient descent with momentum (SGD with momentum) is designed to accelerate learning (Robbins and Monro, 1951). SGD with momentum can avoid training into saddle points (Lee et al., 2016) and improve network generalization performance (Hardt et al., 2015; Wilson et al., 2017). SGD scales the gradient uniformly in all directions to determine the descending step size, which can be particularly harmful to ill-conditioned problems. Therefore, SGD needs to frequently modify the learning rate according to the actual situation. To address this issue, adaptive methods such as Adam (Kingma and Ba, 2015), Adagrad (Duchi et al., 2011), and RMSprop (Tieleman and Hinton, 2012) have been proposed that adaptively correct the learning rate during training.

Although the convergence speed and generalization ability of Adam and other adaptive methods are better than SGD in the initial stage of training, their performance in the convergence part has stagnated. A more natural strategy is to use the Adam algorithm to initialize the training, which allows the model to converge quickly and then convert to the SGD with momentum when appropriate (Keskar and Socher, 2017).

The satellite Bouguer gravity data are downloaded from the website (https://bgi.obs-mip.fr/data-products/grids-and-models/wgm2012-global-model/) with the resolution of 2′×2′ and the range of E65^°-95^°, N40^°-55^°. We obtained the regional gravity anomaly by upward continuing the satellite Bouguer gravity anomaly (Figure 1A) to 10 km. We subtract the regional field from the total field to obtain the residual gravity anomaly (Figure 1B). Considering that the resolution of the satellite Bouguer gravity anomaly data used in this paper is 2′×2′, the matching 1:5000000 Asia-European geological map (Figure 2) is used as a priori information for data annotation. In order to accurately depict the boundary contours of the geologic body and accurately classify it by stratum, we annotate the training data at the pixel level. The labeling map adopts the index map mode. The specific categories and index values are shown in Table 1.

In order to take into account the semantic segmentation accuracy of both small and large targets and provide more global semantic information to the network, we use a multi-scale sliding window clipping method (Table 2) to clip the data of the entire region. In order to preserve the details of the remaining gravity anomalies, all image sizes are 2048×2048 pixels. A total of 61 samples with bedrock outcrops less than 10% of the total area of the whole map were used as a test set. There are 475 effective samples participating in network training, 15% are randomly selected as the verification set, and the remaining are the training set.

The feature encoder uses a convolution operation to encode the input 3-channel RGB image into a high-dimensional feature map with 1/8 pixel size and 96 channels (Oktay et al., 2018). It consists of three convolution modules, each of which contains a convolution with the kernel size 3×3 and stride = 2, followed by a batch normalization layer (BN) (Ioffe et al., 2015) and a rectification linear unit (ReLU).

In semantic segmentation task, spatial location information, contextual semantic information, and receptive fields are crucial for segmentation accuracy. The increase of network depth can obtain better contextual semantic information. The skip connections (He et al., 2015; 2016) can enrich spatial location information. The network depth, the size of the convolution kernel, and the position of the skip connections will affect the receptive field of the feature map used for segmentation. Increasing the width of the network can improve the receptive field of the feature map, at the same time, the number of network parameters also increase dramatically. With limited data sets, it is difficult to train a large network with strong generalization capabilities. It is a challenging task to make the network have rich contextual semantic information, precise spatial location information, and sufficient receptive field, which is the ultimate goal of network design.

With the increase of GPU computing power, neural architecture search algorithms are more and more widely used. The Efficientnet (Tan et al., 2019), MobilenetV2 (Sandler et al., 2018), Auto-DeepLab (Liu et al., 2020) are all designed using neural architecture search algorithms. They have achieved good performance in image classification, detection, segmentation and other tasks. We adopt the neural architecture search algorithm (Zoph and Quoc, 2016; Brock et al., 2017; Bender et al., 2018; Pham et al., 2018; Gong et al., 2019) to design multi-resolution feature extraction module. Inspired by Chen, (2019), a multi-resolution feature extraction module (MRFEM) is obtained by searching directly on the target dataset. In order to take into account the experience of manual design and the flexibility of neural architecture search, we designed the search space consists of 5 operators (Table 3).

The zoomed convolution, proposed by Chen, (2019), reduces the size of the input feature map by bilinear downsampling, followed by a convolution operation, and finally restores the output to the original input size by bilinear upsampling. This special design enjoys a lower calculation amount and 2 times larger receptive field compared to standard convolution.

During the neural network search, we fix the network depth of a multi-resolution branch and simultaneously search the paths of three resolution branches with the sampling rates of 1/8, 1/16, and 1/32. For each layer of a single branch, the expansion rate of the network width can be any value in {2, 4, 6, 8}. The operation type can be any one in the search space, and the position of the skip connections can be selected arbitrarily. Since we are more concerned with the accuracy of network segmentation, we use weighted bootstrapped cross-entropy loss function in the search process.

In order to provide the maximum receptive field with global context information, we use Attention Refinement Module (Yu et al., 2018) to refine the features of each stage. Attention Refinement Module (ARM) employs global average pooling to remove the redundant information of the feature map (Figure 3B). It extracts the global context semantic information through the convolution operation with a kernel size of 1×1. The attention vector is calculated through the sigmoid function which is merged into the output to guide the feature learning during the training process.

The features of the three paths are different in level of feature representation. The 1/8 resolution branch represents the relatively macroscopic and superficial semantic information, while the 1/16 and 1/32 resolution branches focus on high-level semantic information such as microscopic details and inter-pixel relationships. Therefore, we cannot simply sum up these features. Moreover, we also need to introduce spatial position information for each resolution branch to achieve precise pixel positioning. We employ the skip connection to integrate the original image information into the output feature map, forming a new feature map that contains rich semantic information and accurate spatial location information. Therefore, we employ a specific Feature Fusion Module (Figure 3C) (Yu et al., 2018) to fuse these features.

We first concatenate the output features of each branch and then utilize the batch normalization to balance the scales of the features. Next, we use global pooling to transform the concatenated feature to a feature vector. We compute a weight vector through the convolution operation with a kernel size of 1×1. The weight vector can re-weight the features, which amounts to feature selection and combination.

The output of the semantic segmentation module is a matrix with the same size as the original image, and the number of channels is equal to the number of categories. Each element of the matrix stores the category of the current pixel. A bilinear up-sampling of the feature map after feature fusion is required to restore its size to the original map. We bilinear upsample the 1/8 resolution feature map to 1/4 size, and cascade it with the original 1/4 resolution feature map through a long skip connection. The feature map is sequentially processed with the deformable convolution (Dai et al., 2017) with a kernel size of 3×3, batch normalization layer (BN) and rectification linear unit (ReLU). By doing the same operation on the feature map with 1/4 and 1/2 resolution, the final feature map with the same size as the original pixel is obtained. The network output is converted to pixel classification results by a 1×1 convolution layer.

We trained all the models directly on the target data set without pre-training. We initialized the network weight parameters by random initialization. We use 4 GPUs for parallel training, with each card having a batch size of 4 and the training period is 600 epochs. The specific training details are as follows.

For the problem of bedrock prediction in covered areas, which is the focus of this article, it can be formulated as a semantic segmentation task. The metrics employed to evaluate performance include Mean Pixel Accuracy (MPA) and Intersection over Union (IoU). Given the category imbalance present in our dataset, to more objectively assess the network’s predictive efficacy, in addition to the aforementioned evaluation criteria, we use Mean Pixel Accuracy (MPA) and Frequency Weighted Intersection over Union (FWIoU) as the evaluation metrics in order to evaluate the network performance more objectively.

Mean Pixel Accuracy is the average ratio for all categories of samples between the number of correctly classified pixels and the total number of pixels. MPA = 1 k ∑ i = 0 k − 1 p i i ∑ j = 0 k p i j (20)

Where k is the total number of categories; p_ij is the number of pixels of class i but is predicted to be class j.

FWIoU is a weighted summation on the IOU_i of each category and its weight w_i. Where w_i is calculated according to the frequency of each class in the data set. IoU = area g r o u n d t r u t h ∩ a r e a p r e d i c t i o n area g r o u n d t r u t h ∪ a r e a p r e d i c t i o n (21) FWIoU = ∑ i = 0 k − 1 w i · IoU i (22)

Under the aforementioned training strategy and initial conditions, the proposed GGMNet behaves stably during the training process (Figure 4). At the beginning of training, the Adam algorithm adaptively adjusts the learning rate, the network converges rapidly, the loss function curve decreases linearly (Figure 4A), and the frequency-weighted cross-parallel curve rises rapidly (Figure 4B). When the loss function curve tends to flatten, the optimization algorithm switches to the SGD algorithm with momentum. After many iterations of the network, the FWIoU curve becomes smooth which means that the network performance is almost saturated. It is time to terminate network training.

The residual Bouguer gravity anomalies in the Junggar Basin and its surrounding areas are feed into the trained GGMNet for prediction. Figure 6 is a visualization of the prediction results. In the northern margin of the Junggar Basin, there is a local high gravity anomaly in the south of Fuhai-Fuyun and north of Kelamayi-Mulei. The GGMNet predicts the coexistence of Carboniferous and Devonian strata. The North-East trending Carboniferous and Devonian strata are also widespread in the bedrock outcrops at the periphery of the basin. The local low gravity anomaly is predicted to be granite. In the south of Karamay-Mulei, there is a local negative anomaly area, which is predicted to be Proterozoic, Cambrian strata and Granitic acidic intrusive rocks.

Through deep neural network prediction, the distribution of concealed formations or rock bodies throughout the entire Junggar Basin has been ascertained. Within the basin, the major north-northwest trending Kalamayi-Mulei-Hami fault serves as a boundary; to its north, Carboniferous and Devonian strata predominate, while to its south, a mixture of Proterozoic, Cambrian, and Ordovician strata dominate, with acid granite intrusions occurring along the stratigraphic interfaces. The prediction indicates that in the Fuhai area along the northern margin of the Junggar Basin, primarily Carboniferous and Devonian strata are present, intermixed with granitic intrusions. This outcome is consistent with the deep geological structures revealed by traditional geophysical methods.

Studies on the crystalline basement of the Junggar Basin suggest the possible existence of a continental crustal crystalline basement, with its lower portion consisting of strata predating the Neoproterozoic and an upper section characterized by widely distributed Devonian and Carboniferous folded basements. This aligns well with the distribution of concealed geological bodies beneath the Junggar Basin’s cover as predicted by our deep neural network model, indicating a high degree of reasonableness in our predictive results. This concurrence signifies that the predictions made in this study offer a credible depiction of the basin’s sub-surface geology.