Well Logging Based Lithology Identification Model Establishment Under Data Drift: A Transfer Learning Method

Abstract

Recent years have witnessed the development of the applications of machine learning technologies to well logging-based lithology identification. Most of the existing work assumes that the well loggings gathered from different wells share the same probability distribution; however, the variations in sedimentary environment and well-logging technique might cause the data drift problem; i.e., data of different wells have different probability distributions. Therefore, the model trained on old wells does not perform well in predicting the lithologies in newly-coming wells, which motivates us to propose a transfer learning method named the data drift joint adaptation extreme learning machine (DDJA-ELM) to increase the accuracy of the old model applying to new wells. In such a method, three key points, i.e., the project mean maximum mean discrepancy, joint distribution domain adaptation, and manifold regularization, are incorporated into extreme learning machine. As found experimentally in multiple wells in Jiyang Depression, Bohai Bay Basin, DDJA-ELM could significantly increase the accuracy of an old model when identifying the lithologies in new wells.

1. Introduction

Well logging data have grown dramatically over the past few decades due to the widespread deployment of oil wells and the rapid development of sensing technology. A large amount of data not only brings us more opportunities to understand the underground but also brings more great challenges to the interpretation of logging data [1,2]. Well logging provides an objective and continuous method with which to observe the properties of the rocks through which the drill bit passes and to describe the deposition process quantitatively. As a bridge between surface geophysical survey and subsurface geology, well logging is an effective and irreplaceable method to understand reservoir characteristics. As conventional reservoirs dry up, oil and gas companies are turning to unconventional exploration and development in shale and low-permeability reservoirs, posing more challenges for logging interpretation [3].

In traditional logging interpretation, lithology determination, porosity, and permeability calculations are performed by specialists in exploration and geology with specialized knowledge. Lithology identification is a fundamental problem in well logging interpretation and is of considerable significance in petroleum exploration engineering. It is the basis of reservoir parameter calculation (such as porosity, shale volume, permeability) and geological research (such as formation correlation, sedimentary modeling, favorable zone prediction, etc.). Due to the complexity of reservoir geological conditions, the uncertainty of exploration data and the inconsistency of expert experience, the results of lithology identification are mainly dependent on expertise. With the increasing diversity of logging data, the traditional logging interpretation methods which rely on human experience have some shortcomings and limitations. As a result, researchers are turning to more advanced data analysis methods for breakthroughs in lithology identification.

In recent years, with the rapid development of machine learning technology, its application in lithology identification has also attracted full attention [4]. In order to accurately and effectively determine the lithology type, a large amount of application work based on different machine learning began to emerge [5]. For example, G. Askari [6] used satellite remote sensing data to find the lithological information in Deh-Molla sedimentary succession by principal component analysis. Al-Anazi et al. [7] proposed a support vector machine-based classification feature method and selection method based on fuzzy theory to realize the recognition of potential features and the improvement of lithology recognition performance. Wang et al. [8] proposed a novel back propagation (BP) model by modifying the self-adapting algorithm and activation function, which is proven to be effective in predicting the lithologies of the Kela-2 gas field. The experimental results show that, compared with discriminant analysis and probabilistic neural network, support vector machines can identify different lithology of heterogeneous sandstone reservoirs more accurately. Xie et al. [9] evaluated five typical machine learning methods of naive Bayes, support vector machines, artificial neural networks, random forest and gradient tree enhancement based on the formation lithology identification data of the Danudui gas field and the Hangjinqi gas field. Experimental results show that the integrated method, including random forest and gradient tree enhancement, has a lower prediction error. At the same time, the gradient tree enhancement method has the highest accuracy compared with the other four methods. In order to reduce exploration uncertainty, Bhattacharya et al. [10] compared different methods of facies division and prediction of mudstone in the current conventional logging data, and applied them to the Devonian Bakken and Mamentago–Marsilus formations in North America. In this study, support vector machines (SVM), artificial neural networks (ANN), self-organizing mapping (SFM), and multi-resolution map based clustering (MRMC) are compared experimentally. Dev et al. [11] analyzed the data from China’s Danuci gas field and Hangjinqi gas field by the gradient lifting decision tree system (i.e., XGBoost, LightGBM, and CatBoost) to study formation lithology classification, and compared their performance with random forest, AdaBoost, and gradient boosting machines. Experiments show that LightGBM and CatBoost are the preferred algorithms for lithology classification by using well logging data. A large number of similar studies can be seen in [12,13,14,15,16].

In addition to these directly-applied studies, more and more scholars are focusing on how to improve existing machine learning tools to solve practical problems in lithology identification. Due to the different distribution of underground lithology, there is a severe class imbalance problem in the training data set. Furthermore, Deng et al. [17] introduced a borderline-stroke technique for dealing with unbalanced data, and the results showed that this method could effectively improve the classification accuracy of SVM, especially for the minority classes. There are many free hyper-parameters in machine learning algorithms, and the settings of these hyper-parameters will have significant influences on the performance of lithology identification. Therefore, Saporetti et al. [18] adopted an evolutionary parameter tuning strategy, and combined gradient boosting (GB) with differential evolution (DE) to achieve the optimization of super parameters, thereby making the lithology identification more stable. In the study of [19], the wavelet decomposition was used to construct multi-channel images of logging data, and then the lithology identification problem based on the logging curve was skillfully transformed into an image segmentation problem. Finally, the feasibility of this method was verified by the application in the Daqing oilfield. Aiming at the issue of data drift between different wells, Ao et al. [20] proposed a hybrid algorithm for lithology identification that combines the mean drift algorithm and the random forest algorithm in the prototype similarity space. It is pointed out that a more accurate lithology identification model can be obtained by transforming the classification problem into prototype similarity space. Li et al. [21] proposed a semi-supervised algorithm based on Generated Adversarial Network, which uses logging curves as the labeled data and seismic data as the unlabeled data. The model was trained by Adam algorithm and uses the discriminator to identify the lithologies. Compared with the experimental results of various supervised methods, the model can effectively use unlabeled data to achieve higher prediction accuracy. Similar work can be found in [22,23,24].

Although much work has been done, a practical problem, i.e., the well loggings from different wells are different in the probability distribution, is not taken into consideration. Hence, the model trained on old wells might not perform well on new wells. As illustrated in Figure 1, the phenomenon of data drift occurs between two wells, even when they are geospatially near. In particular, applying the model trained on well $W_b$ to well $W_a$ has many errors, as shown in Figure 1a. To suppress the data drift-induced accuracy decrease, we propose a transfer learning method named the data drift joint adaptation extreme learning machine (DDJA-ELM) to increase the accuracy of the old model applying to new wells. By incorporating the project mean maximum mean discrepancy, joint distribution domain adaptation, and manifold regularization into the extreme learning machine, we realize the knowledge transfer from old wells to new wells. As experimented in multiple wells in Jiyang Depression, Bohai Bay Basin, DDJA-ELM could increase the accuracy of an old model significantly when identifying the lithologies in new wells. In the remainder of the paper, Section 2 expatiates on the proposed DDJA-ELM, which is evaluated in Section 3. The last section concludes the paper.

2. Methodology

2.1. Notation

The dataset $\mathcal{D}={\left\{({\mathbf{x}}_i,{\mathbf{y}}_i)\right\}}_{i=1}^n$ is composed of the sample ${\mathbf{x}}_i={[x_1,\dots ,x_d]}^{\prime }\in {\mathbb{R}}^d$ with d dimension and the label ${\mathbf{y}}_i=[{\underbrace{0,\dots ,0,1}}_k,0,\dots ,0]\in \{0,1\}$ if ${\mathbf{x}}_i$ belongs to the k-th $(1\le k\le c)$ class, where n and c are the numbers of samples and classes, respectively. In well logging-based lithology identification, d generally means the number of loggings types, and c denotes the number of lithology types. A sample is composed of the logging values at a certain depth. Considering the problem of data drift, we use ${\mathcal{D}}_S$ and ${\mathcal{D}}_T$ to differentiate the datasets with drift; i.e., ${\mathcal{D}}_S$ and ${\mathcal{D}}_T$ represent the source dataset for training and the target dataset to predict, respectively. Specifically, the source dataset is labeled and the target dataset is unlabeled, thereby denoting ${\mathcal{D}}_S={\left\{({\mathbf{x}}_i^s,{\mathbf{y}}_i^s)\right\}}_{i=1}^{n_s}$ and ${\mathcal{D}}_T={\left\{{\mathbf{x}}_i^t\right\}}_{i=1}^{n_t}$ with $n_s$ and $n_t$ samples accordingly, and $n_s+n_t=n_{s+t}$, where ${\mathbf{x}}_i^s$ and ${\mathbf{x}}_i^t$ are source and target-dataset samples, respectively; and ${\mathbf{y}}_i^s$ is a source-dataset label of ${\mathbf{x}}_i^s$. By defining ${\mathbf{X}}_S=[{\mathbf{x}}_1^s;\dots ;{\mathbf{x}}_{n_s}^s]$, ${\mathbf{Y}}_S=[{\mathbf{y}}_1^s;\dots ;{\mathbf{y}}_{n_s}^s]$, and ${\mathbf{X}}_T=[{\mathbf{x}}_1^t;\dots ;{\mathbf{x}}_{n_t}^t]$, then the ${\mathcal{D}}_S=\left\{({\mathbf{X}}_S,{\mathbf{Y}}_S)\right\}$ and ${\mathcal{D}}_T=\left\{{\mathbf{X}}_T\right\}$.

2.2. Problem Definition and Formulation

In general, a classifier $f\left(x\right)$ trained on $\mathcal{D}$ classifies well on those samples that satisfied the independent identically distributed (i.i.d.) assumption. However, the classifier $f\left(x\right)$ trained on ${\mathcal{D}}_S$ might not perform well on ${\mathcal{D}}_T$ because the data drift means that the i.i.d. assumption does not hold; i.e., the knowledge learned from ${\mathcal{D}}_S$ cannot transfer to ${\mathcal{D}}_T$. In this case, it is necessary to add some new constraint conditions to achieve expected performance while learning the classifier f. Exactly, the expected performance can be concluded as the following constraint conditions: (i) minimizing the structural risk to avoid overfitting; (ii) minimizing the data drift between the source dataset and the target dataset; (iii) minimizing the prediction inconsistency within the target dataset. If considering the ELM as the basic classifier f, $f\left(\mathbf{x}\right)=f(\mathbf{x};\mathbf{\beta })$ where $\mathbf{\beta }$ is the output weight matrix of ELM.

According to the structural risk minimization (SRM) and regularization techniques [26], the classifier f can be denoted as

$$f=\underset{\mathbf{\beta }}{argmin}\ell (f(\mathbf{x};\mathbf{\beta }),\mathbf{y})+R(\mathbf{x};\mathbf{\beta }),$$

(1)

where the first term represents the empirical loss on samples (i.e., describing the fitness of applying the model to predict the training data); the second term represents the regularization term (i.e., representing the formulation of constraint conditions). Thus, combined with the constraint conditions mentioned above, the objective function is formulated mathematically as follows

$$f=\underset{\mathbf{\beta }}{argmin}\ell (f(\mathbf{x};\mathbf{\beta }),\mathbf{y})+\frac{1}{2}{∥\mathbf{\beta }∥}^2+\frac{\lambda }{2}\Omega ({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })+\frac{\gamma }{2}\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta }),$$

(2)

where the empirical loss term $\ell \left(f\right(\mathbf{x};\mathbf{\beta }),\mathbf{y})$ and the structural risk regularization term ${∥\mathbf{\beta }∥}^2$ compose the model accuracy and complexity of ELM. The term $\Omega ({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ indicates the extent of data drift between ${\mathcal{D}}_S$ and ${\mathcal{D}}_T$. Additionally, we introduce the manifold regularization term $\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })$ to improve the prediction consistency within the target dataset. $\lambda$ and $\gamma$ are the regularization parameters accordingly.

In the remainder of this section, we will expatiate on each term of the objective function.

2.2.1. ELM

Recent years have witnessed the development of a promising machine learning model; i.e., extreme learning machine (ELM). ELM is actually an artificial neural network with a single hidden layer, which was first proposed by Huang et al. [27] and found its application in many domains, such as robotic perception, hyperspectral image classification, lithology identification, and human activity recognition [28,29,30,31,32]. Compared with support vector machine and other artificial neural networks, ELM has significant superiority in generalization performance and training time. In addition, many variants for ELM have been investigated, including semi-supervised ELM, multi-kernel ELM, rough ELM, one-class ELM, etc. [33,34,35,36].

According to the objective function (2), we first describe the mathematical model of ELM as follows [37],

$$\begin{array}{cc}\hfill & \underset{\mathbf{\beta }\in {\mathbb{R}}^{z\times c}}{argmin}\ell (f(\mathbf{x};\mathbf{\beta }),\mathbf{y})+\frac{1}{2}{∥\mathbf{\beta }∥}^2=\frac{C}{2}\sum _{i=1}^{n_{s+t}}{∥{\mathbf{e}}_i∥}^2+\frac{1}{2}{∥\mathbf{\beta }∥}^2\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.j_i\mathbf{h}\left({\mathbf{x}}_i\right)\mathbf{\beta }={\mathbf{y}}_i-{\mathbf{e}}_i,i=1,\dots ,n_{s+t},\hfill \end{array}$$

(3b)

where ${\mathbf{e}}_i\in {\mathbb{R}}^c$ represents the error vector with respect to the i-th training sample, and ${\sum }_{i=1}^{n_{s+t}}{∥{\mathbf{e}}_i∥}^2$ is the sum of prediction errors. The tradeoff coefficient C is used to balance the contribution between the two terms. Since the target dataset ${\mathcal{D}}_T$ is unlabeled, $j_i=1$ if $i\le n_s$; otherwise, $j_i=0$; and ${\mathbf{y}}_i$ equals a zero vector if $i>n_s$. The output of hidden layer is equal to $\mathbf{h}\left({\mathbf{x}}_i\right)=[\mathbf{g}({\mathbf{x}}_i;{\mathbf{w}}_1,b_1),\dots ,\mathbf{g}({\mathbf{x}}_i;{\mathbf{w}}_z,b_z)],$ where ${\mathbf{w}}_j\in {\mathbb{R}}^z$; z is the number of hidden neurons; $b_j\in \mathbb{R}$ are the j-th random generated weight vector and bias constant; $\mathbf{g}(·)$ is a piecewise continuous nonlinear activation function, such as the sigmoid function $\mathbf{g}({\mathbf{x}}_i;{\mathbf{w}}_j,b_j)=1/\{1+exp(-{\mathbf{w}}_j^{\prime }·{\mathbf{x}}_i-b_j)\}$ or the Gaussian function $\mathbf{g}({\mathbf{x}}_i;{\mathbf{w}}_j,b_j)=exp(-b_j\parallel {\mathbf{x}}_i-{\mathbf{w}}_j{\parallel }^2)$. In this paper, the sigmoid function is used as the activation function.

2.2.2. Weighted ELM

Additionally, since each lithology in dataset $\mathcal{D}$ usually contains the samples of different amounts, appropriate weights should be assigned to each error vectors for the samples imbalance issue, so (3) is rewritten as

$$\begin{array}{cc}\hfill & \underset{\mathbf{\beta }\in {\mathbb{R}}^{z\times c}}{argmin}\ell (f(\mathbf{x};\mathbf{\beta }),\mathbf{y})+\frac{1}{2}{∥\mathbf{\beta }∥}^2=\frac{C}{2}\sum _{i=1}^{n_{s+t}}{\omega }_i{∥{\mathbf{e}}_i∥}^2+\frac{1}{2}{∥\mathbf{\beta }∥}^2\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.j_i\mathbf{h}\left({\mathbf{x}}_i\right)\mathbf{\beta }={\mathbf{y}}_i-{\mathbf{e}}_i,i=1,\dots ,n_{s+t},\hfill \end{array}$$

(4b)

where ${\omega }_i=\frac{{\nu }_i}{{\sum }_{i=1}^{n_{s+t}}{\nu }_i}$ denotes the weight, ${\nu }_i=\frac{1}{n_i^{\tau }}$ and $\tau$ is a constant, and $n_i$ is the number of samples belonging to $y_i$, thereby modifying the ELM to the weighted ELM (WELM).

Substituting the constraints (4b) into (4a) yields the equivalent unconstrained optimization problem and the matrix form:

$$\underset{\mathbf{\beta }\in R^{z\times c}}{argmin}\frac{C}{2}\mathbf{W}{∥\left(\mathbf{Y}-\mathbf{J}\mathbf{H}\mathbf{\beta }\right)∥}^2+\frac{1}{2}{∥\mathbf{\beta }∥}^2,$$

(5)

where $\mathbf{W}=\mathrm{blockdiag}({\mathbf{W}}_S,{\mathbf{0}}_{n_t\times n_t})\in {\mathbb{R}}^{n_{s+t}\times n_{s+t}}$; ${\mathbf{W}}_S\in {\mathbb{R}}^{n_s\times n_s}$ is a diagonal matrix with each item ${\omega }_i$; and $\mathbf{Y}=[{\mathbf{Y}}_S;{\mathbf{0}}_{n_t\times c}]\in {\mathbb{R}}^{n_{s+t}\times c}$, $\mathbf{H}=[\mathbf{h}\left({\mathbf{x}}_1\right);\dots ;\mathbf{h}\left({\mathbf{x}}_{n_{s+t}}\right)]\in {\mathbb{R}}^{n_{s+t}\times z}$, $\mathbf{J}=\mathrm{blockdiag}({\mathbf{I}}_{n_s\times n_s},{\mathbf{0}}_{n_t\times n_t})\in {\mathbb{R}}^{n_{s+t}\times n_{s+t}}$, where $\mathbf{0}$ and $\mathbf{I}$ are the zero matrix and the identity matrix with appropriate dimension, respectively. By setting the gradient of (5) over $\mathbf{\beta }$ to zero, we have

$$\mathbf{\beta }-C{\mathbf{H}}^{\prime }{\mathbf{J}}^{\prime }\mathbf{W}(\mathbf{Y}-\mathbf{J}\mathbf{H}\mathbf{\beta })=0,$$

(6)

where ′ represents the transpose of matrix. There are two forms of an optimal solution to the ${\mathbf{\beta }}^{*}$. When the number of the training data is larger or smaller than the number of hidden neurons, then $\mathbf{H}$ has more or fewer rows than columns, thereby resulting in an underdetermined or overdetermined least squares problem. The closed-form solution for (6) can be described as

$${\mathbf{\beta }}^{*}=\left\{\begin{array}{c}\hfill {\mathbf{H}}^{\prime }\mathbf{W}{\left(\mathbf{H}{\mathbf{H}}^{\prime }\mathbf{W}+\frac{1}{C}{\mathbf{I}}_{n_{s+t}}\right)}^{-1}\mathbf{Y},\mathrm{if}{n}_{s+t}<z\\ \hfill {\left({\mathbf{H}}^{\prime }\mathbf{W}\mathbf{H}+\frac{1}{C}{\mathbf{I}}_z\right)}^{-1}{\mathbf{H}}^{\prime }\mathbf{W}\mathbf{Y},\mathrm{if}{n}_{s+t}\ge z\end{array}\right.,$$

(7)

where ${\mathbf{I}}_{n_{s+t}}$ and ${\mathbf{I}}_z$ are the identity matrix of $n_{s+t}$ and z dimensions.

2.2.3. Data Drift Adaptation

In general, the data distribution properties can be used to statistically describe the correlation among the samples $({\mathbf{x}}_i,{\mathbf{y}}_i)\in \mathcal{D}$, thus the data drift-induced distribution discrepancy can be estimated by statistical criteria. In this paper, we utilize the projected maximum mean discrepancy (PMMD) criterion for distribution discrepancy measure [38]. Thus, the regularization term $\Omega ({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ with PMMD is formulated as

$${\Omega }_{\mathrm{P}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })={∥\frac{1}{n_s}\sum _{i=1}^{n_s}\mathbf{h}\left({\mathbf{x}}_i^s\right)\mathbf{\beta }-\frac{1}{n_t}\sum _{j=1}^{n_t}\mathbf{h}\left({\mathbf{x}}_j^t\right)\mathbf{\beta }∥}^2,$$

(8)

Further, (8) can be rewritten to a matrix form; that is,

$${\Omega }_{\mathrm{P}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })=\mathrm{tr}\left({\mathbf{\beta }}^{\prime }{\mathbf{H}}^{\prime }{\mathbf{M}}_{\mathrm{P}}\mathbf{H}\mathbf{\beta }\right),$$

(9)

where $\mathrm{tr}(·)$ denotes the trace of a matrix. The elements of matrix ${\mathbf{M}}_{\mathrm{P}}\in {\mathbb{R}}^{n_{s+t}\times n_{s+t}}$ are calculated by

$$\begin{array}{c}\hfill {\left({\mathbf{M}}_{\mathrm{P}}\right)}_{pq}=\left\{\begin{array}{cc}\frac{1}{n_s^2},\hfill & \mathrm{if}p,q\le n_s\hfill \\ \frac{1}{n_t^2},\hfill & \mathrm{if}p,q>n_s\hfill \\ -\frac{1}{n_s{n}_t},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(10)

According to (8), the samples ${\mathbf{x}}_i^s$ and ${\mathbf{x}}_j^t$ are transformed from the feature space to the mapping space by $\mathbf{h}(·)$, where the distribution discrepancy between the source and target datasets can be reduced by adjusting $\mathbf{\beta }$. Thus, the modified WELM named data drift adaptation WELM (DDA-WELM) classifier can be achieved by introducing the regularization term ${\Omega }_{\mathrm{P}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ to adapt the data drift.

However, the overall accuracy is increased with sacrificing the accuracies of some classes, because DDA-WELM is a method of reducing the distribution discrepancy between two datasets as a whole. Consequently, we modify the PMMD to propose the joint PMMD (JPMMD) that aims to reduce the distribution discrepancy between classes of two datasets, respectively. The regularization term $\Omega ({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ with JPMMD is formulated as

$${\Omega }_{\mathrm{J}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })=\sum _{k=1}^c{∥\frac{1}{n_s^{\left(k\right)}}\sum _{{\mathbf{y}}_i^s=k}\mathbf{h}\left({\mathbf{x}}_i^s\right)\mathbf{\beta }-\frac{1}{n_t^{\left(k\right)}}\sum _{{\tilde{\mathbf{y}}}_j^t=k}\mathbf{h}\left({\mathbf{x}}_i^s\right)\mathbf{\beta }∥}^2,$$

(11)

where ${\mathbf{y}}_i^s$ is the real label of ${\mathbf{x}}_i^s$ and ${\tilde{\mathbf{y}}}_j^t$ is the pseudo label of ${\mathbf{x}}_j^t$ which is generated by the DDA-WELM classifier; $n_s^{\left(k\right)}$ and $n_t^{\left(k\right)}$ are the number of source-dataset samples which belong to the class k and the number of target-dataset samples whose pseudo labels are k, respectively.

Similarly, (11) can be rewritten to the matrix form; that is,

$${\Omega }_{\mathrm{J}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })=\mathrm{tr}\left({\mathbf{\beta }}^{\prime }{\mathbf{H}}^{\prime }{\mathbf{M}}_{\mathrm{J}}\mathbf{H}\mathbf{\beta }\right),$$

(12)

where the elements of matrix ${\mathbf{M}}_{\mathrm{J}}\in {\mathbb{R}}^{n_{s+t}\times n_{s+t}}$ are computed by

$${\left({\mathbf{M}}_{\mathrm{J}}^{\left(k\right)}\right)}_{pq}=\left\{\begin{array}{cccc}\hfill & \frac{1}{{\left(n_s^{\left(k\right)}\right)}^2},\hfill & \hfill & \mathrm{if}p,q\le n_s\mathrm{and}{\mathbf{y}}_p^s,{\mathbf{y}}_q^s=k\hfill \\ \hfill & -\frac{1}{n_s^{\left(k\right)}{n}_t^{\left(k\right)}},\hfill & \hfill & \mathrm{if}p>n_s,q\le n_s\mathrm{and}{\tilde{\mathbf{y}}}_{p-n_s}^t,{\mathbf{y}}_q^s=k\hfill \\ \hfill & -\frac{1}{n_s^{\left(k\right)}{n}_t^{\left(k\right)}},\hfill & \hfill & \mathrm{if}p\le n_s,q>n_s\mathrm{and}{\mathbf{y}}_p^s,{\tilde{\mathbf{y}}}_{q-n_s}^t=k\hfill \\ \hfill & \frac{1}{{\left(n_t^{\left(k\right)}\right)}^2},\hfill & \hfill & \mathrm{if}p,q>n_s\mathrm{and}{\tilde{\mathbf{y}}}_{p-n_s}^t,{\tilde{\mathbf{y}}}_{q-n_s}^t=k\hfill \\ \hfill & 0,\hfill & \hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(13)

Thus, the data drift joint adaptation WELM (DDJA-WELM) classifier can be obtained by introducing the regularization term ${\Omega }_{\mathrm{J}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ and the DDA-WELM classifier to improve the WELM.

2.2.4. Manifold Regularization

Manifold regularization is widely used to improve the smoothness of predictions and suppress the classifier cutting through the high-density regions [39]. Specifically, we introduce the manifold regularization $\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })$ to assign smooth labels within the target dataset and make the classifier more adaptable to the target dataset.

The formulation of the manifold regularization $\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })$ is that

$$\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })=\frac{1}{2n_t^2}\sum _{i=1}^{n_t}\sum _{j=1}^{n_t}{a}_{i,j}{∥\mathbf{h}\left({\mathbf{x}}_i^t\right)\mathbf{\beta }-\mathbf{h}\left({\mathbf{x}}_j^t\right)\mathbf{\beta }∥}^2$$

(14)

where the similarity $a_{i,j}$ between the sample ${\mathbf{x}}_i^t$ and ${\mathbf{x}}_j^t$ is calculated by

$$\begin{array}{c}\hfill a_{i,j}=\left\{\begin{array}{cccc}\hfill & exp\left(\frac{-{∥{\mathbf{x}}_i^t-{\mathbf{x}}_j^t∥}^2}{4{\sigma }^2}\right),\hfill & \hfill & {\mathbf{x}}_i^t\in \mathcal{N}\left({\mathbf{x}}_j^t\right)\mathrm{or}{\mathbf{x}}_j^t\in \mathcal{N}\left({\mathbf{x}}_i^t\right),\hfill \\ \hfill & 0,\hfill & \hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(15)

where $\mathcal{N}\left({\mathbf{x}}_j^t\right)$ is the set of $\kappa$-nearest neighboring samples of ${\mathbf{x}}_i^t$ under the metric of Euclidean distance in feature space, and $\sigma >0$ is the width of Guassian kernel. Additionally, (14) can be rewritten to a matrix form

$$\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })=\mathrm{tr}\left({\mathbf{\beta }}^{\prime }{\mathbf{H}}^{\prime }\mathbf{L}\mathbf{H}\mathbf{\beta }\right).$$

(16)

where $\mathbf{L}=\mathrm{diag}({\mathbf{0}}_{n_s\times n_s},{\mathbf{L}}_T)\in {\mathbb{R}}^{n_{s+t}\times n_{s+t}}$, ${\mathbf{L}}_T=\mathbf{D}-\mathbf{A}$ is the target dataset Laplacian matrix, and $\mathbf{A}={\left[a_{i,j}\right]}_{n_t\times n_t}$, $\mathbf{D}$ is a diagonal matrix with diagonal elements $d_i=\sum _{j=1}^{n_t}{a}_{i,j}$.

Thus, the classifiers DDA-WELM and DDJA-WELM can be upgraded to the DDA-S2WELM and DDJA-S2WELM by introducing the manifold regularization term $\mathcal{M}({\mathcal{D}}_T;\mathbf{\beta })$. Moreover, the regularization term ${\Omega }_{\mathrm{J}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ with pseudo labels generated by the DDA-S2WELM classifier can further improve the data drift adaptation performance of the DDJA-WELM and DDJA-S2WELM classifiers.

2.3. Solution of Objective Function

The solution of the objective function (2) is introduced in this section. Especially, the regularization term $\Omega ({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ is divided into ${\Omega }_{\mathrm{P}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ with PMMD and ${\Omega }_{\mathrm{J}}({\mathcal{D}}_S,{\mathcal{D}}_T;\mathbf{\beta })$ with JPMMD so that we will discuss them in the followings, respectively.

2.3.1. Solution of Objective Function with ${\Omega }_{\mathrm{P}}$

By incorporating (5), (9), and (16) into (2), we have

$$f=\underset{\mathbf{\beta }\in {\mathbb{R}}^{z\times c}}{argmin}\frac{C}{2}\mathbf{W}{∥\mathbf{Y}-\mathbf{J}\mathbf{H}\mathbf{\beta }∥}^2+\frac{1}{2}{∥\mathbf{\beta }∥}^2+\frac{\lambda }{2}\mathrm{tr}\left({\mathbf{\beta }}^{\prime }{\mathbf{H}}^{\prime }{\mathbf{M}}_{\mathrm{P}}\mathbf{H}\mathbf{\beta }\right)+\frac{\gamma }{2}\mathrm{tr}\left({\mathbf{\beta }}^{\prime }{\mathbf{H}}^{\prime }\mathbf{L}\mathbf{H}\mathbf{\beta }\right).$$

(17)

By setting the gradient of f with respect to $\mathbf{\beta }$ to be zero

$$\nabla f=C{\mathbf{H}}^{\prime }{\mathbf{J}}^{\prime }\mathbf{W}(\mathbf{J}\mathbf{H}\mathbf{\beta }-\mathbf{Y})+\mathbf{\beta }+\lambda {\mathbf{H}}^{\prime }{\mathbf{M}}_{\mathrm{P}}\mathbf{H}\mathbf{\beta }+\gamma {\mathbf{H}}^{\prime }\mathbf{L}\mathbf{H}\mathbf{\beta }=0.$$

(18)

According to (18), the closed-form solution of the optimal $\mathbf{\beta }$ is

$${\mathbf{\beta }}^{*}=\left\{\begin{array}{c}\hfill {\mathbf{H}}^{\prime }\mathbf{W}{\left(\frac{1}{C}{\mathbf{I}}_n+\left(\mathbf{W}+\frac{\lambda }{C}{\mathbf{M}}_{\mathrm{P}}+\frac{\gamma }{C}\mathbf{L}\right)\mathbf{H}{\mathbf{H}}^{\prime }\right)}^{-1}\mathbf{Y},\mathrm{if}{n}_s+n_t<z\\ \hfill {\left(\frac{1}{C}{\mathbf{I}}_z+{\mathbf{H}}^{\prime }\left(\mathbf{W}+\frac{\lambda }{C}{\mathbf{M}}_{\mathrm{P}}+\frac{\gamma }{C}\mathbf{L}\right)\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\prime }\mathbf{W}\mathbf{Y},\mathrm{if}{n}_s+n_t\ge z\end{array}\right.$$

(19)

where z is the number of hidden neurons.

2.3.2. Solution of Objective Function with ${\Omega }_{\mathrm{J}}$

Similarly, incorporating (5), (12), and (16) into (2), and solving $\frac{\partial f}{\partial \mathbf{\beta }}=0$ yields the optimal $\mathbf{\beta }$ as follows

$${\mathbf{\beta }}^{*}=\left\{\begin{array}{c}\hfill {\mathbf{H}}^{\prime }\mathbf{W}{\left(\frac{1}{C}{\mathbf{I}}_n+\left(\mathbf{W}+\frac{\lambda }{C}{\mathbf{M}}_{\mathrm{J}}+\frac{\gamma }{C}\mathbf{L}\right)\mathbf{H}{\mathbf{H}}^{\prime }\right)}^{-1}\mathbf{Y},\mathrm{if}{n}_s+n_t<z\\ \hfill {\left(\frac{1}{C}{\mathbf{I}}_z+{\mathbf{H}}^{\prime }\left(\mathbf{W}+\frac{\lambda }{C}{\mathbf{M}}_{\mathrm{J}}+\frac{\gamma }{C}\mathbf{L}\right)\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\prime }\mathbf{W}\mathbf{Y},\mathrm{if}{n}_s+n_t\ge z\end{array}\right.$$

(20)

where z is the number of hidden neurons.

3. Experimental Verification

In this section, we conduct extensive experiments to verify the effectiveness of our method, using the well-logging data collected from multiple wells in Jiyang Depression, Bohai Bay Basin. The experimental datasets and settings will be described first. Then, the performance and impact of each regularization terms are shown in detail. The analysis of hyper-parameters sensibility is presented at last.

3.1. Experimental Settings

As shown in

Table 1. Dataset description.

Dataset	1				2				3
LOGS	AC, CAL, CNL, GR, RT, SP				AC, CAL, GR, R25, SP				AC, CAL, GR, R25, SP
Samples	LITH	Mu	Si	Co	LITH	Mu	Sa	Do	LITH	Mu	Sa	Sh
	Well A	3820	2048	452	Well D	4508	3076	1000	Well F	6072	2852	1644
	Well B	2164	1584	432	Well E	5000	1464	1112	Well G	6996	3404	1548
	Well C	2368	1452	400	-				-

[1] Logging curves. AC: acoustic log, CAL: caliper log, CNL: compensated neutron log, GR: gamma ray log, RT: true formation resistivity, SP: spontaneous potential log, R25: 2.5 m bottom gradient resistivity. [2] Lithology. Mu: mudstone, Si: siltstone, Co: conglomeratic sandstone, Sa: sandstone, Do: dolomite, Sh: shale.

, the experimental datasets are composed of three datasets collected from different regions in the Jiyang Depression, Bohai Bay Basin and contain 3, 2 and 2 wells, respectively (The wells in Figure 1 do not appear in the experimental datasets). Their relative positions are shown in Figure 2.

Table 2. Data drift description.

	A and B	A and C	B and C	D and E	F and G
MMD Class 1	0.107	0.182	0.281	0.109	0.090
MMD Class 2	0.230	0.067	0.277	0.238	0.251
MMD Class 3	0.203	0.026	0.384	0.189	0.601
MMD Overall	0.540	0.275	0.942	0.536	0.942

[1] MMD is short for maximum mean discrepancy, which can be found in [40]. [2] For (A,B), (A,C), and (B,C), classes 1–3 mean Mu, Si, and Co, respectively. For (D,E), classes 1–3 mean Mu, Sa, and Do, respectively. For (F,G), classes 1–3 mean Mu, Sa, and Sh, respectively.

describes the data drift statistically by maximum mean discrepancy. In one experiment, we set the well A as the training data and the well B as the testing data, thereby verifying the effectiveness of our method through A → B. In this case, the other experiments are denoted as A → C, B → A, B → C, C → A, C → B, D → E, E → D, F → G, and G → F, respectively. Considering the different ranges of measurements, the value of samples are transformed to [0, 1] by a min-max normalization method. In the experiments, we adopt the index Recall (i.e., the number of classify correctly divided by the number of samples) to calculate the classification accuracy, thus denoting the average recall (Macro-R) to represent the overall classification accuracy.

3.2. Experimental Results

Figure 3 exhibits the logging curves, core (i.e., label), classification performance before B → A transfer (i.e., using model trained on data B to directly predict data A), and classification performance after B → A transfer. It is observed that transfer learning could eliminate the data drift-induced accuracy loss. As shown in Figure 4, Figure 5, and Figure 6, the classification performances of applying the classifiers ELM, S2ELM, DDA-ELM, DDA-S2ELM, DDJA-ELM, and DDJA-S2ELM to datasets 1, 2, and 3 are presented, respectively. It can be obviously observed from Figure 4a–c,g–i, Figure 5a,c, and Figure 6a,c, that the accuracy of each class is gradually increased with the introducing of DDJA regularization terms, especially when it comes to the accuracy of Si and Sh that increase from 0 to more than 80 in Figure 4i, Figure 6a,c, respectively. Additionally, the ELM and S2ELM classifiers without DDJA regularization term are insufficient to the tasks with data drift.

Considering the ELM–based classifier with the random weights and biases parameters, we conduct multiple experiments by setting different random seeds to generate different weights and biases parameters which yield the experimental results in Figure 4d–f,j–l, Figure 5b,d and Figure 6b,d. According to the results, we have the following observations: (i) The overall Macro-R shows a trend of increasing step by step and the DDJA-S2ELM classifiers achieve the highest accuracy. Moreover, the Macro-R of DDJA-S2ELM are increased to 88% at least compared with the ELM on dataset 1 (B → C and C→ A), 80% on dataset 2 (D → E), and 69% on dataset 3 (F →G and G → F). Additionally, compared the ELM classifier with the DDJA-S2ELM classifier, the Macro-R of our method are increased 52% on dataset 1 (C→ B), 10% on dataset 2 (E→ D), and 42% on dataset 3 (F→ G). (ii) Although the Macro-R is not increased on dataset 1 (A → C) and dataset 2 (D → E), the standard deviation can be kept at a lower level. Thus, the stability can be improved by introducing the DDJA-S2ELM. (iii) Comparing the experimental results on dataset 3 with datasets 1 and 2, the ELM classifiers only achieve 27% (F→ G) and 35% (G→ F). The performance is significantly enhanced with our method regarding the data drift-induced accuracy decrease surges.

3.3. Parameters’ Sensitivity

In this section, we aim to analyze the sensitivity on the key hyper-parameters of the DDJA-S2ELM: the trade-off coefficient C, the contribution coefficient of DDJA regularization term $\lambda$, and the contribution coefficient of semi-supervised regularization term $\gamma$. By analyzing the Macro-R over different settings of hyper-parameters, the configuration of these parameters is given. To avoid repetition, we only show the results of B → C.

In the Figure 7, we set the hyper-parameter C range from 100 to 10,000,000. It can be seen from these figures in Figure 7 that the maximum and minimum Macro-R are increased first and then decreased with the increasing of C, indicating that a small C incurs under-fitting in classification and a high C causes the over-fitting. Additionally, when C is set too small or too big, the overall accuracies are preserved in the large variation range highest 94.9% in Figure 7a and lowest 41.1% in Figure 7a or the overall just more than 70% sightly in Figure 7f. According to Figure 7c,d with $C=1000$ and $C=10000$, respectively, the overall Macro-R are held steady highest 96.9% in Figure 7c and lowest 73.8% in Figure 7c and the accuracies are presented a trend of increasing gradually. Therefore, it is important to configure a moderate C first.

The influence of adjusting $\lambda$ and $\gamma$ are given by fixing the C. It can be seen from Figure 7 that the maximum Macro-R is achieved 96.9% in Figure 7c when set $\lambda /C=10,000$ and $\gamma /C=100$. Additionally, the Macro-R almost increase first and then decrease with $\gamma$ getting larger under fixed $\lambda$ showing in Figure 7a–c. Moreover, these results show that the maximum is achieved almost when setting the $\lambda$ larger than $\gamma$ by two to four orders of magnitude. Since the $\gamma$ control the contribution of semi-supervised regularization term that based on the manifold assumption, but the assumption is invalid when the dataset with data drift. We introduce the DDJA regularization term to suppress the data drift, so the coefficient $\lambda$ should be more than $\gamma$.

According to the analysis of setting C, $\lambda$, and $\gamma$ mentioned above, it can be concluded that: (i) C should be set a relatively larger range first to configure a moderate value; (ii) the setting of $\lambda$ should be bigger than $\gamma$ by two to four orders of magnitude. Here are some suggested settings: $C\in [{10}^3,{10}^6]$, $\frac{\lambda }{C}\in [{10}^3,{10}^4]$, $\frac{\gamma }{C}\in [{10}^2,{10}^3]$.

Furthermore, we use Sobol method (In our experiment, the python tool "SALib" is employed, which can be found at https://salib.readthedocs.io/en/latest/index.html). Ref. [41] to implement a global sensitivity analysis where all parameters are varied simultaneously over the entire parameter space. The sensitivities shown in

Table 3. Sensitivity analysis using Sobol method.

Param.	First-Order Sensitivity	Total Sensitivity		Param. 1	Param. 2	Second-Order Sensitivity
C	0.511	0.575		C	$\lambda$	0.074
$\lambda$	0.337	0.495		C	$\gamma$	0.008
$\gamma$	0.051	0.102		$\lambda$	$\gamma$	0.085

demonstrate that: (i) C and $\lambda$ contribute mainly compared with $\gamma$, (ii) C contributes more than $\lambda$ slightly, and (iii) the interactions between these parameters are weak, so they are relatively independent.

4. Conclusions

In this paper, we have investigated the well logging-based lithology identification under data drift, thus proposing a transfer Extreme Learning Machine method to handle it. According to the projected maximum mean discrepancy (PMMD) criterion and extreme learning machine, we introduce a new PMMD criterion and then propose the DDJA-ELM to minimize the data drift between the source dataset and the target dataset. Additionally, in order to improve the prediction consistency within the target dataset, the manifold regularization is introduced to promote the DDJA-ELM to the DDJA-S2ELM. Extensive experiments have validated the high and stable accuracy of our method.