Lithology Identification of Uranium ‐ Bearing Sand Bodies Using Logging Data Based on a BP Neural Network

: Lithology identification is an essential fact for delineating uranium ‐ bearing sandstone bodies. A new method is provided to delineate sandstone bodies by a lithological automatic classi ‐ fication model using machine learning techniques, which could also improve the efficiency of bore ‐ hole core logging. In this contribution, the BP neural network model for automatic lithology identi ‐ fication was established using an optimized gradient descent algorithm based on the neural network training of 4578 sets of well logging data (including lithology, density, resistivity, natural gamma, well ‐ diameter, natural potential, etc.) from 8 boreholes of the Tarangaole uranium deposit in Inner Mongolia. The softmax activation function and the cross ‐ entropy loss function are used for lithology classification and weight adjustment. The lithology identification prediction was carried out for 599 samples, with a prediction accuracy of 88.31%. The prediction results suggest that the model is effi ‐ cient and effective, and that it could be directly applied for automatic lithology identification in sandstone bodies for uranium exploration.


Introduction
Sandstone type uranium prospecting is a comprehensive exploration activity based on various exploration approaches and analytical methods [1][2][3], such as finding uranium sources, delineating uranium-bearing sandstone bodies, locating redox transitional zones, and checking material composition. The delineation of the uranium-bearing sandstone bodies is one of the most critical factors during the exploration of this type of deposit since the uranium orebodies are all hosted within the sandstone bodies. Furthermore, the correct lithology identification in the sandstone bodies is of decisive significance to the development of this type of uranium deposit.
Currently, manual logging by geologists and logging interpretation by geophysicists are the two main ways of lithology identification in the exploration of sandstone-type uranium deposits. In most cases, the logging interpretations are corrected by geological logging, or both are verified mutually [4,5]. Geological logging is commonly of high accuracy and strong reliability but low efficiency. In comparison, the geophysical interpretation is fast and of low cost but with lower accuracy and reliability. The approaches of identifying lithology from logging data include the cross-plot method [6], the curve overlap method [7], the statistical method [8], the structural characteristic parameter method, etc. For example, Wang et al. (2021) established a principal component cross plot of the Chang-8 reservoir of the Triassic Yanchang Formation in the Jiyuan area of the Ordos Basin, and the distribution in the characteristic space of different lithology is directly observed. The distinction boundary of lithology is obvious, and the distribution of lithology conforms to the law of sandstone oil-bearing level. Compared with the sensitive logging curve, the cross plot established with principal components F1 and F2 has advantages for lithology identification. Zhao (2020) has a good application effect on the identification of a low-resistivity gas reservoir in the Dongsheng Gas Field by using the logging curve overlap method and the four porosity difference/ratio method. Both methods are better than the cross-plot method. The interpretation coincidence rate of the conventional logging method for identifying low-resistivity gas reservoirs reaches 83%. Yukun et al. (2013) used the P-wave velocity and the density parameters for lithology identification through the Markov random field neighborhood system method. This method also has some limitations. The Markov random field prior model produces the identification result of a block, and the identified stratum needs to emphasize a continuous change.
With the successive development and improvement of machine learning technology, artificial neural network models represented by back propagation algorithm (BP) [9,10] and convolution neural network algorithm (CNN) [11][12][13] have been widely used in various fields. The CNN algorithm has many excellent performances in image classification and recognition, text classification, speech recognition, and unmanned driving recognition [14,15]. In contrast, the BP algorithm has been used in shallow learning and deep learning, and it is used in power equipment fault diagnosis, atmospheric environment monitoring, radar signal recognition, lithology identification, and other technical applications have performed well. Qu (2021) used the method based on the BP neural network model to identify the lithology of 873 rock samples from the Fauqi oilfield in Missan Province in southeast Iraq, which has a higher efficiency and accuracy than the conventional method. In particular, the recognition rate of dolomite and dolomitic sandstone has increased from 84.1% and 77.1% to 97.8% and 87.5%, respectively [16]. Ying (2015) estimated the relationship between gas content and natural gamma, acoustic moveout, density, and compensated neutron in five coalbed methane exploration wells in the Galilee Basin of Australia, by multiple linear regression and the BP neural network method. The interpretation accuracy of the BP neural network method is generally high, and the interpretation effect is affected by the number and the quality of input layer samples. The more samples, the stronger the regional representation, and the better the interpretation effect [17]. Taking the Yin'e Basin as an example, Xue et al. compared the accuracy of lithology identification by both cross plot and a BP neural network. The results showed that the coincidence rate of the prediction using a BP neural network was significantly improved [18]. The BP neural network model has achieved sound application results in the above practice-especially in the forward propagation process-through linear and nonlinear coupling, the input layer data can reach the output layer smoothly, but most of the references have less analysis of the backward propagation process. The first-order optimization methods such as stochastic gradient descent (SGD), AdaGrad, and Adam are easy to fall into a local minimum rather than a global minimum for weight parameter estimation [19][20][21][22][23][24]. The second-order optimization methods [25,26], such as the Newton method, quasi-Newton method, and conjugate gradient method have the disadvantages of high complexity, significant difficulty in constructing a Hessian matrix, a large amount of calculation, and a slow training speed when there are too many hidden layers.
The logging data of the sandstone-type uranium deposits are of dense logging points of which the sampling interval is 5 cm, with a large amount of data. Moreover, the analysis of the iterative adjustment of the weight parameters of backward propagation in lithology prediction is poor. The main contributions of this manuscript are as follows: ① The mathematical formulas of weight iteration are deduced in detail with examples in the geology domain, which provide theoretical support for geological staff to deeply understand the mathematical principle of a BP neural network model. ② An optimized gradient algorithm based on the AdaGrad method is proposed to improve the computational efficiency of the training model.
③ The codes to realize the lithology prediction are self-written, and it is convenient to analyze the factors affecting the lithology classification by tracking the trend of the estimated parameters.

Methodology
The realization of the BP neural network consists of two parts: forward propagation and backward propagation. Forward propagation refers to the sample vector X = [x1, x2, x3, …, xn] as the input layer, and each element of the sample vector is weighted (h = w T X+ b) and nonlinear activated (σ) by the hidden layer, and a prediction vector ( ) is finally output by the output layers ( Figure 1). Backward propagation refers to adjusting the errors calculated by comparing the predicted value of the output layer with the actual value through processing the training samples iteratively. The adjustment process is realized by calculating the weight through the loss function until the predicted value is infinitely close to the actual value.

Forward Propagation
Input by a BP neural network, Vector X has n samples, and each sample has m features. The weight value from the input layer to the hidden layer is set to w, and the activation function (σ) is used from the hidden layer to the output layer. Common activation functions include sigmoid, tanh, ReLu, softmax, etc. In uranium exploration, it is necessary to distinguish between the permeable layer and the non-permeable layer. The Sigmoid function is the best choice if there are only two classifications; if different types of lithology are identified according to the logging data and the lithology of the same section is mutually exclusive classification, the softmax function can realize the output of categories k with k > 2.
The linear regression formula from the input layer to the hidden layer: where h denotes the output value of the hidden layer, X denotes the input vector, w denotes the connection weight, b denotes the bias value, and its initial value is set to 1: where n denotes the number of samples, m denotes the feature of each sample, and k denotes the category of lithology. Namely: The softmax activation function formula of the hidden layer: where 1,2, ⋯ , , k is the category of lithology, and ℎ is the output of the neural node, is the n-th predicted value.

Backward Propagation
According to the difference between the predicted value and the actual value, from the output layer to the hidden layer, and then from the hidden layer to the input layer, the weight value of each layer is adjusted step by step to improve the prediction accuracy. The process is to adjust the weight data of each neuron by constructing a loss function and calculating the descending gradient. The common loss functions include the quadratic loss function [27], exponential loss function [28], hinge loss function [29], and cross-entropy loss function [30]. The derivation of a cross-entropy loss function is relatively simple in the process of gradient descent calculation, and there is no gradient dispersion. The loss function formula is based on the softmax function formula: where k denotes the category of lithology; denotes the i-th true value of training samples; and denotes the i-th predicted value. The optimization of the cross-entropy loss function is realized by calculating the descent gradient. The gradient descent formula is: where w is the weight matrix; n denotes the number of samples; j is the label of lithology; y is the true value of training samples; and p denotes the probability of the neural node of the output layer. 1{•} is the indicator function: To obtain the gradient of the loss function of the i-th output value at t-th iteration, the formula can be expressed in the following mathematical form, From Equations (3) and (6), and according to the chain rule, we have where 1 The in Equation (7) is calculated: Furthermore, applying Equations (7)-(10), we have

Optimized Gradient Descent Algorithm
In the process of updating the weights with backward propagation, the stochastic gradient descent method is the most commonly used method as a first-order optimization algorithm. The current optimization algorithms improved on the basis of the stochastic gradient descent algorithm include Momentum, AdaGrad, RMSProp, Adam, etc. In order to adapt to the characteristics of logging data of sandstone type uranium deposits with highly intensive sampling (sample interval of 5 cm), accelerate the convergence speed of the algorithm, and improve its stability, the "learning rate" is optimized based on AdaGrad algorithm. To increase the step size in the initial stage of optimization, the logarithmic operation is performed on the cumulative square gradient. Taking the function of as an example, while the initial values are set to (4,16) (Figure 2a,b), the logarithmic operation of the optimized algorithm tends to the minimum value of the function earlier (Figure 2b) after the same iteration of 5 steps, indicating that the optimized algorithm converges faster.

Regional Geology
The studied Tarangaole uranium deposit in the Ordos basin lies within the Yinshan-Liaohe uranium province ( Figure 3). Geographically, it is located in the northwest of the Ordos City in the northeast margin of the Ordos basin, extending in an east-west direction. Tectonically, the studied area is hosted within the Yimeng uplift in the north of the Ordos Basin [31].
The basement of the Ordos Basin has the characteristics of a "dual" basement, composed of the direct basement of the upper Paleozoic Permian (P2)-Mesoproterozoic (Pt2) age and the Paleoproterozoic (Pt1)-middle Archean (Ar2) crystalline basement. The cover rocks in the Ordos Basin are mainly Cenozoic and Mesozoic sediments, among which the Jurassic Zhiluo Formation (J2z) is the main uranium orebody-bearing prospecting target in the northeast of the Ordos Basin, followed by the Jurassic Yan'an Formation (J2a). The logging data selected in this manuscript are mainly from the Zhiluo Formation.

Deposit Geology
The strata revealed by drilling at the Tarangaole uranium deposit includes the Lower Cretaceous and the Middle Jurassic Zhiluo Formation. The ore-bearing sedimentary layers are dominated by the Lower Member of the Middle Jurassic Zhiluo Formation, which is further divided into the upper sub-member (J2z1-2) and the lower sub-member (J2z1-1). The water-retaining mudstone layers are only locally developed between the upper and the lower sub-member sandstones (Figure 4), and thus they are still classified as the same sand body as a whole. The sand body is mainly composed of green-gray medium-grained and coarse-grained sandstone, with some mudstone layers locally identified. The sandstone debris in the Lower Member of Zhiluo Formation at the Tarangaole deposit is mainly composed of the feldspar sandstone (77.76% in volume), the feldspar quartz sandstone (21.35% in volume), and a small proportion of quartz sandstone and graywacke sandstone, accounting for 0.78% and 0.11%, respectively. Feldspar sandstone accounts for a large proportion, while quartz sandstone and greywacke account for a small proportion, indicating that the Lower Member of Zhiluo Formation was deposited close to the source area, and the hydrodynamic conditions were relatively stable during deposition. Quartz is the main component of clasts, accounting for 68.22% of the total amount of clasts. There is no obvious alteration, and it is only partially replaced by carbonate in the strongly carbonated rocks.
Feldspar accounts for roughly 26.22% of the total debris. Potassium feldspar is the main component of feldspar, accounting for approximately two thirds of the total feldspar, followed by plagioclase. Most of the feldspar debris have not undergone obvious alteration, some feldspars are strongly clayed, the potassium feldspar is mainly kaolinized, and the plagioclase is hydromica, chlorite, and epidote (Figure 6a). Clayified feldspar accounts for about one third of the total feldspar.
Metamorphic rock clasts are mainly composed of quartzite and mica quartz schist (Figure 6d). The granite debris and the volcanic debris are observed as well. Spinel quartz, mica, iron minerals, quartz, feldspar, and heavy minerals with dissolution structures can be seen in volcanic debris.
The content of the mica varies greatly, from 2.0% to 5.0%, with an average of 3.66%, and up to 8-10% locally. Mica is biotite, with a small amount of muscovite. Biotite is brown biotite, followed by green biotite, which is mostly the product of brown biotite alteration. Biotite is often subjected to chloritization in varying degrees (Figure 6c) to become chlorite.
Organic debris are unevenly distributed in sandstone with content generally less than 0.50%; however, up to 5.0-10.0% in local sandstone, and the organic veinlets or strips are observed (Figure 6b).

Data Characteristics
The data used in this manuscript are the logging and the geological logging data from eight boreholes at the Tarangaole deposit. The logging data include 6 logging characteristic variables, such as density, resistivity, natural gamma, natural potential, well-diameter and gamma, with a total of 4578 sets of samples. The JHQ-2D digital comprehensive logging tool and the FD-3019 logging probe were used to carry out natural gamma and gamma logging. The sampling depth was between 470-500 m, with a sampling interval of 5 cm. The ZKC232-87 borehole is shown as the representative to check the consistency and the reliability of the characteristics of six logging variables (Table 1 and Figure 7). Different logging methods have different response characteristics to the same lithology. The more shale content in sandstone bodies, the higher the radioactivity intensity. Therefore, natural gamma logging is very sensitive to shale content. Gamma logging is the quantitative measurement of uranium and thorium content in sandstone bodies. Lithologic density logging classifies lithology by calculating the photoelectric absorption crosssection index of sandstone bodies. Spontaneous potential can divide permeable formation, determine its interface, and estimate the formation water resistivity. Well-diameter measurement can assist in judging mudstone layers and loose rock layers. Due to the well wall collapse caused by mud immersion and scouring during drilling, the actual welldiameter is larger than the bit diameter, resulting in well-diameter expansion.  Six lithology types were identified by geological logging (Figure 8), including conglomerate (CG), coarse-grained sandstone (CS), medium-grained sandstone (MS), finegrained sandstone (FS), calcareous sandstone (KS), mudstone (MD), etc. The content of medium sandstone was the highest (up to 38.14%), followed by coarse sandstone (36.98%), with minor mudstone and calcareous sandstone (2.94% in total) ( Table 2).  In the model training process, the logging data from seven boreholes were selected as training samples and the data from the other borehole data were selected as test samples.

Model Building
According to the data characteristics in Figure 7, the trends of natural gamma and gamma logging data are almost the same. Furthermore, the correlation coefficient between the natural gamma and the gamma ray is 0.9332, as shown in Figure 9. Both are strongly correlated, and they have the same influence weight on the training model. Only one of them is selected as one of the attribute data. Therefore, the natural gamma (GP), density (DEN), spontaneous potential (SP), well-diameter (Φ), and resistivity (ρ) of the logging data were selected as the input layer variables of the BP neural network; let X denotes variables, X = {GP, Φ, DEN, SP, ρ}. The nodes of the hidden layer network were set to be five, let H denote the value of node, H = w T X + b. The main lithology of the Tarangaole deposit revealed by drilling is conglomerate (CG), coarse-grained sandstone (CS), medium-grained sandstone (MS), fine-grained sandstone (FS), calcareous sandstone (KS), and mudstone (MD). The output layer variables were set to be six, let Y denote this variable, Y = {CG,CS,MS,FS,KS,MG}. From the hidden layer to the output layer, the classification was realized by the softmax function, and the w weight value was adjusted by the cross-entropy loss function. The optimization algorithm of O_AdaGrad proposed in this paper was used to accelerate the calculation of the w weight.

Model Training
The computing environment was based on Python 3.7, while the numpy module, the pandas module, the sklearn.neural_network module, the sklearn.model_selection model and the sklearn.metrics modelwere imported to perform the training and the testing data. The logging and the lithology data of seven boreholes, including zkc232-87, zkc236-87, zkc232-91, zkc264-87, zkc248-95, zkc252-87, and zkc264-79 in eight boreholes were selected as training data, and the zkc240-87 boreholes were selected as the verification data set. Due to density (DEN), resistivity (ρ), Natural gamma (GR), well-diameter (Φ), and spontaneous potential (SP) there was a large variation range and different measurement units. In order to eliminate the impact caused by the difference, logarithmic normalization was carried out for the training data and the verification data at the same time, as shown by this mathematical formula: (12) where X is the processed log data; is the raw log data; is the minimum value of the original logging data; and is the maximum value of the original logging data. Codes and labels were, respectively, set for six types of lithological data (Table 3), including conglomerate (CG), coarse-grained sandstone (CS), medium-grained sandstone (MS), fine-grained sandstone (FS), calcareous sandstone (KS), and mudstone (MD). The maximum number of iterations was set to be 10,000, the initial value of the weight index (w) was 1, the initial value of the weight decay index was 0.005, the bias value (b) was 1, and the initial value of the learning rate (η) was 0.1. The softmax function of formula (2) and the cross-entropy loss function of formula (3) were called to optimize and to classify the weight. The values of cross-entropy loss function were 1.195, 0.897, 0.797, 0.713, and 0.713 when iterating 100, 1000, 5000, 10,000, and 12,000 times, respectively. The value was still 0.713 (Figure 10a) when the number of iterations increased to 12,000, indicating that the optimal solution had been found and the weight update reached a reasonable position. The weights of five attribute values, such as DEN, ρ, GR, Φ, and SP were distributed on each type of lithology after 10,000 iterations, as shown in Figure 10b. After 100, 1000, 5000, 10,000, and 12,000 iterations, the training accuracy rates of 3979 training samples were 63.99%, 79.42%, 89.67%, 89.72%, and 89.73%, respectively. The comparison of the original lithology and the predicted lithology between 10,000 and 12,000 iterations, as shown in the Figure 10c, was almost the same. It indicated that the weight reached a reasonable value after 10,000 iterations.

Prediction Validation
The lithology of 599 tested samples from the borehole ZKC240-87 were predicted using the trained BP neural network model to verify the reliability of this model ( Figure 11). The confusion matrix was created based on the testing data ( Table 4). The overall prediction accuracy rate on the testing data was 88.31%. The average sensitivity of the testing data was 74.83%, and the sensitivity of the coarse-grained sandstone and the mediumgrained sandstone were 90.91% and 94.44%, respectively. The average specificity of the testing data was 92.23%, and the specificity of the conglomerate and the mudstone were both 100.00%. The average precision of the testing data was 79.34%. The results show that the coarse-grained sandstone and the medium-grained sandstone as "true positive" and the conglomerate and the mudstone as "true negative" are correctly predicted with high accuracy. The low prediction results of mudstone and fine-grained sandstone are attributed to the small training samples in the training model, indicating that the uniformity of data distribution of the training model has an important impact on the prediction results. The predicted lithology of 599 tested samples, which are identified by the trained BP neural network model, are almost the same as the core lithology ( Figure 12). The reliability of the BP neural network model is proved.

Conclusions
This paper proposes an optimized algorithm based on the AdaGrad method to improve the BP neural network model that is employed to predict the lithology of the sandstone-type uranium deposit. 1. Through comparative analysis of logging data, the parameters of natural gamma (GP), density (DEN), natural potential (SP), resistivity (ρ), and well-diameter (Φ) are selected to establish the BP neural network used for the lithology identification of the Tarangaole uranium deposit. A total of 4578 sample data were selected from 8 boreholes to train the BP neural network model, and 599 samples for data verification. 2. The prediction accuracy rate using a BP neural network model reaches 88.31% in this contribution, supporting the effectiveness of the BP neural network model in predicting lithology in the studied area. The average sensitivity and the average specificity of the testing data are 74.83% and 92.23%, respectively. 3. The loss function value, the weight value, and the prediction accuracy under different iteration times in the prediction process are reported in detail, which is favorable to the objective evaluation of the influencing factors of data changes in the process of training and prediction. Data Availability Statement: Not applicable.