SVR has incomparable advantages in data mining of small sample data sets. Considering that there may be a small amount of data in some data sets after the construction of sub-data sets, this paper chooses SVR as the basic data mining algorithm for TOC logging interpretation. The basic concept of SVR method is to project the input data into a higher dimension by kernel function, so as to find a hyperplane to establish a regression function. For a given data set {(x ₁ ,y ₁ ),……, (x _l ,y _l ) }, where x _i ∈R ⁿ is the input data, y _i ∈R ¹ is the target output value, and the SVR estimation function is: f x = w T · Φ x + b , (2) where w and b are hyperplane parameters, Φ(x) denotes the eigenvectors after x projection. The standard form of SVR for solving hyperplane parameters is (Vapnik, 1998): min w , b ,   ξ ,   ξ *   1 2 w T w + C ∑ i = 1 l ξ i + ξ i * . (3a)

受制于.

Subject to w T · Φ x i + b − y i ≤ ε + ξ i , (3b) y i − w T · Φ x i − b ≤ ϵ + ξ i * , ξ i ,   ξ i * ≥ 0 ,   i = 1 , ⋯ , l , where C is the penalty coefficient or regularization parameter, and ε, ξ,ξ*∈R are slack variables introduced to penalize the fitting function. Eq. 1 can be transformed into a dual problem to solve, and the original problem is transformed into its corresponding Lagrangian function form, and by minimizing: min α , α *   1 2 ∑ i = 1 l ∑ j = 1 l α i − α i * α j − α j * K x i , x j + ε ∑ i = 1 l α i + α i * + ∑ i = 1 l y i α i − α i * . (4a)

Subject to ∑ i = 1 l α i − α i * = 0 , (4b) 0 ≤ α i ,   α i * ≤ C ,   i = 1 , ⋯ , l , where α _i = (α ₁ , α ₂ ,…α _l ) is the Lagrange multiplier, and K (x _i ,x _j ) is the kernel function. The final regression equation is: f x = ∑ i = 1 n α i − α i * K x i , x j + b . (5)

In this paper, polynomial kernel function, radial basis kernel function (RBF) and sigmod kernel function are selected to explain TOC respectively, so as to optimize the best kernel function type.

XGboost was first proposed by Chen and Guestrin (2016). It is a machine learning algorithm that relies on the boosting principle and explores weak learners to comprehensively predict. This is mainly due to its well-known high prediction accuracy. Its basic principle is to generate a sub-classifier to fit the prediction residuals of the previous sub-classifiers, thereby continuously reducing the residuals between the true value and the predicted value, and finally integrating all sub-classifiers to give the final prediction result. The expression is: y ^ i = ∑ k = 1 K f k x i   ,   f k ∈ F . (6)

Among them, y ^ i is the calculated predictive value of the i th sample; K denotes the number of decision trees; f _k denotes the kth submodel; x _i represents the input feature of the ith sample; F represents the set of sub-classifiers. In the Xgboost sub-classifier, the Classification and Regression Tree is usually selected. In the Xgboost algorithm, the objective function is composed of a loss function and regularization parameters. The expression is: L φ = ∑ i l y ^ i , y i + ∑ k Ω f k ,   w h e r e   Ω f = γ T + 1 2 λ w 2   . (7)

Among them, l y ^ i , y i is the residual between the predicted value y ^ i and the target value y _i ; f _k is the function expression of the k sub-classifier; Ω f k is the penalty term of the model complexity, which can be used to smooth the final learned weights to avoid overfitting. XGboost is trained iteratively to obtain an approximation of L(φ). Assuming that the sub-classifier trained in the t iteration is f _t , after the t iteration, the objective function can be expressed as: L t = ∑ i = 1 n l y ^ i t − 1 + f t x i , y i + Ω f k . (8)

Equation. 8 can be further optimized using second-order approximation: L t ≅ ∑ i = 1 n l y ^ i t − 1 , y i + g i f t x i + 1 2 h i f t 2 x i + Ω f k , (9) where, g _i and h _i are the first-order and second-order partial derivatives (gradients) of l, respectively, where g i = ∂ y ^ t − 1 l y ^ i t − 1 , y i 和 h i = ∂ y ^ t − 1 2 l y ^ i t − 1 , y i . After taking them into Eq. 9 and removing the constant term, we can obtain: L ∼ t = ∑ i = 1 n g i f t x i + 1 2 h i f t 2 x i + Ω f t . (10)

Define I _j ={i | q (x _i ) = j} as an instance set of leaf node j. Eq. 10 is rewritten by extending Ω to: L ∼ t = ∑ i = 1 n g i f t x i + 1 2 h i f t 2 x i + γ T + 1 2 λ ∑ j = 1 T w j 2 = ∑ j = 1 n ∑ i ∈ I j g i w j + 1 2 ∑ i ∈ I j h i + λ w j 2 + γ T . (11)

Therefore, the objective function is transformed into a function of the first and second partial derivatives of the loss function l, the leaf node weight, and the number of leaf nodes. In the case of fixed tree structure q(x) the optimal weight w _j *of leaf node j can be calculated by the following formula: w j * = − ∑ i ∈ I j g i ∑ i ∈ I j h i + λ . (12)

The optimal solution formula of the objective function is as follows: L ∼ t q = − 1 2 ∑ j = 1 T ∑ i ∈ I j g i 2 ∑ i ∈ I j h i + λ + γ T . (13)

XGboost iteratively adds branches to construct sub-classifiers on the initial leaf nodes through a greedy algorithm to determine the optimal tree structure of the CART tree. Suppose there is a leaf node, I _L and I _R are instances of the left and right nodes after the node is branched. Let I = I _L ∪I _R , then the loss after branching is reduced to: L s p l i t = 1 2 ∑ i ∈ I L g i 2 ∑ i ∈ I L h i + λ + ∑ i ∈ I R g i 2 ∑ i ∈ I R h i + λ − ∑ i ∈ I g i 2 ∑ i ∈ I h i + λ − γ . (14)

If L _split is greater than 0, the objective function decreases after the leaf node is split into two leaf nodes, so as to determine the node segmentation. On this basis, XGboost is optimized by feature pre-ranking, quantile approximation, and parallel lookup to quickly find the nearest split point.

In the SVM and XGboost algorithm, there are a large number of hyper-parameters, which will affect the final prediction results. Therefore, it is necessary to use hyper-parameter optimization algorithm to determine which hyper-parameter system the SVM and XGboost algorithm can achieve the best prediction results. This study used Genetic Algorithm to optimize hyperparameters.

Genetic algorithm was first proposed by Holland (1973). It is a parallel stochastic optimization algorithm developed from the simulation of natural genetic mechanism and biological evolution theory. The genetic algorithm starts with a set of randomly generated parameters to be optimized, which is called the initial population, where each parameter pair is called an individual. Genetic algorithm encodes each individual in series to form chromosome, and determines the fitness function according to the optimization objective to calculate the fitness of each individual. Several individuals with high fitness values are selected from the initial population, and the chromosomes encoded by these individuals are crossed and mutated to form a new generation of individual populations. Then the fitness of each individual in the new population is calculated, and the above operations are performed repeatedly until the target value or the maximum number of iterations satisfying the fitness is met. In the iterative process, the genetic algorithm can preserve the individuals with good fitness values and eliminate the individuals with poor fitness. The new population not only inherits the information of the previous generation, but also is superior to the previous generation. Through continuous iteration, the parameters can be optimized.

Based on core observation, according to the difference of sedimentary structure and structure, five lithofacies are developed in the shale of Yanchang Formation, which are argillaceous shale facies (AS), siltstone/very fine sandstone facies (SS), tuff facies (TUF), argillaceous shale facies with silty lamina (ASLS) and argillaceous shale facies with tuff lamina (ATLS). The above five rock facies are easy to identify at the core scale. AS are black, grayish black, fine particles (Figure 3A), and do not develop or develop a small amount of silty or tuffaceous layers; SS is mainly gray and grayish white, and a very small amount of grayish black argillaceous bands are developed (Figure 3B); TUF is grayish yellow, easily broken (Figure 3C), relatively homogeneous, and basically does not develop other lithologic layers; ASLS are mainly gray-black argillaceous shale, with a large number of gray-white and gray silty layers distributed inside. The thickness of these layers is generally millimeter and centimeter (Figure 3D), and the cumulative thickness of silty layer accounts for 20%–50%. The main body of ATLS is black argillaceous shale, with a large number of yellow or grayish yellow tuffaceous laminae distributed inside. The laminae thickness is generally in the millimeter and centimeter levels (Figure 3E). The cumulative thickness of the tuffaceous layer accounts for 20%–50%. Based on the above principles, a columnar distribution map of rock facies in 12 wells was drawn.

In the TOC frequency distribution diagram of different lithofacies shown in Figure 4, there are differences in the distribution range of TOC in different lithofacies. The TOC of TUF and SS is low, mainly distributed below 2.0 wt%, (Figures 4A,B )and the TOC distribution of AS and ASLS is medium (Figures 4C,D). The average values are 5.03 wt% and 6.39 wt% The above four lithofacies have no obvious exhibit skewed distribution with a long tail, and the data have good balance. In the Yanchang Formation shale, the TOC of the ATLS is generally high, and the numerical distribution range is from 3 wt% to 22 wt% (Figure 4E). The frequency distribution guidance diagram of the lithofacies shows a weak skew distribution. Compared with Figure 1C, the proportion of data greater than 8 wt% is all increased, and the imbalance of data is weakened.

The relative proportion of different rock facies in the shale section of the Yanchang Formation is the main reason for the unbalanced distribution of TOC data in Figure 1C. In the shale interval, AS has the highest proportion of thickness, which can account for the total thickness of the shale interval. Secondly, the TOC distribution characteristics of ASLS samples are similar to those of AS, which causes the overall TOC data to be concentrated in the TOC intervals of the above two lithofacies, while in other distribution intervals, especially in the high-value TOC interval of ATLS, there are fewer samples, resulting in unbalanced distribution of data in Figure 1C. The data imbalance of TOC sub-data set of rock facies obtained by classification is reduced, which is helpful for learning algorithm to obtain more accurate prediction model.

As mentioned above, the relationship between formation logging response and TOC is affected by other formation parameters, including mineral composition, elemental composition, and organic matter type. Table 2 shows the differences in mineral composition and organic matter types between different rock facies.

From the perspective of mineral composition, AS has the characteristics of high clay mineral content, low felsic content and medium pyrite content. TUF and SS are characterized by low clay mineral content, low pyrite and high felsic content. The difference is that TUF has high carbonate content and SS has high carbonate content. ASLS has obvious transitional characteristics between AS and SS, that is, clay mineral content, carbonate mineral content, clay mineral content and pyrite are all at a medium level; the main characteristic of ATSL is the highest content of pyrite, which can reach 8.9% on average, and other minerals are at a medium level.

Through pyrolysis data, it can be seen that S1 and S2 are higher in the three rock phases of ATSL, AS and ASLE, with an average value of more than 3.5 mg/g and 9.5 mg/g. The S1 and S2 values of SS are lower, with an average value of 2.15 mg/g and 4.2 mg/g. The S1 and S2 values of TUF are the lowest, and S1 and S2 are below 1 mg/g. From the S1/TOC × 100 index, the SS value is the highest, reaching an average of 173 mg/g TOC, followed by ASLS (an average of 116 mg/g TOC), The average value of AS and ATSL is about 65 mg/g TOC, and TU is the lowest, only 26 mg/g TOC. S1/TOC is often used to evaluate the oil content in rocks (Jarvie, 2008). Considering that shale oil may adsorb/dissolve in kerogen, the higher S1/TOC value is generally considered to be a higher content of movable oil, that is, the higher oil content in pores (Li et al., 2015).

In addition, Qiu et al. (2014) and Akhtar et al. (2018) studied the geochemical characteristics of tuff layers in the Yanchang Formation of the Ordos Basin and found that tuff layers generally have high U and Th contents. In a comparative study, Lu (2020) and Yin et al. (2017) found that the layers of siltstone or silty lamina in Zhangjiatan shale often have low U and Th content, while argillaceous shale has relatively high U and Th content.

Figure 5 shows the logging response distribution of each rock facies in different TOC intervals, and the trend line is drawn by the connection of 50th perecentile point in each interval. Because tuff generally has the characteristics of hole enlargement (as shown in Figure 3F, YY22 well 1310 m), the logging response value has great uncertainty, so the relevant data of tuff are not drawn. It can be seen from Figure 5 that there are great differences in the logging response trend lines between different TOC intervals in different rock facies, which also shows that the differences in mineral composition, element composition and oil content of different rock facies will affect the relationship between logging response and TOC.

From Figure 5, the relationship between TOC and logging response in different lithofacies can be classified into two categories:

One is that the trend lines are similar in direction but not coincident, as shown in Figures 5A–F. In Figures 5A, F, because ATLS has the highest U and Th content, it often has higher GR and U logging values under the same TOC conditions as other lithofacies. Similar SS and ASTL have lower U and Th than AS, which makes it have lower GR and U values. The difference in mineral composition may be the main reason for the inconsistency of the trend lines in Figures 5B, C. For example, the high density of pyrite and carbonate makes ASTL and SS lithofacies have higher density values under the same TOC, and similar minerals also make ASTL and SS have lower acoustic time difference. The difference in oil content caused the non-coincidence of the trend line in Figure 5D. Oil has a higher H+ content than kerogen, resulting in SS and ASLS with higher S1/TOC under the same TOC. Higher neutron porosity values, on the contrary, ASTL neutron porosity is low.

The second is the difference in the direction of the trend line, as shown in Figure 5E, the resistivity logging response distribution in different TOC intervals. The obvious feature is that the trend line of ATLS lithofacies is not obvious, and even in some TOC intervals, the resistivity decreases with the increase of TOC. The high content of pyrite in ATLS may be a key factor in this phenomenon, which also causes the resistivity of ATLS to be generally lower than that of other rock phases under the same TOC. Secondly, the low content of clay minerals with good conductivity and high oil content also lead to higher resistivity of SS and ASLS than AS.

It can be seen that the introduction of rock facies in the analysis of TOC and logging response relationship can better understand the relationship between TOC and different logging responses, so that the relationship is less affected by shale formation factors such as mineral composition and element composition.

In this study, four SVR models of rock facies were constructed, which were AS, SS, ATLS and ASLS. For three purposes, the prediction model of tuff facies (TUF) was not constructed: 1) the phenomenon of borehole enlargement is obvious in this lithofacies, and the quality of logging data is poor; 2) The proportion of tuff facies in the Yanchang Formation reservoir is low, and the number of TOC test samples is small. 3) The TOC content of the lithofacies is generally low and the values are concentrated (Figure 4A). Datas are derived from FY1, YY18, WY1, YY12, YY2, YY27, YY28, and B36 wells. The total number of data is 412, of which AS, SS, ATLS and ASLS are 214,55,71 and 72 respectively. The above data are randomly assigned to supervised training data sets and validation sets at a ratio of 1:4.

For the need of comparison, this study also constructed a prediction model under the uniform regression fitting mode. The same as the above data, the supervised data did not contain the relevant samples of tuff facies, and the supervised training data set and verification set were obtained from the corresponding data sets of the above four lithofacies. The supervised data include seven kinds of data such as AC, DEN, GR, Rt, PE, Th/K, U/Th, and TOC. The data normalization is carried out by the following formula: y i ′ = y max − y i y max − y min . (23)

Among them, the logging data uses the same maximum and minimum values in the above five supervised data sets. Since the TOC of the samples in SS is much lower than that of the other rock facies, in order to ensure the final prediction accuracy, the TOC of the SS supervised data set is normalized to [0,3], and the TOC of the remaining four supervised data sets is normalized to [0,30], which is normalized to [0,30] in the uniform regression fitting model.

In this study, SVR is used as the basic algorithm, and genetic algorithm is used to optimize the hyper-parameters in SVR. The optimized parameters include kernel function type and its key parameters (Table 3). The fitness function of the genetic algorithm is the cross-validation MSE of the training data (using the K-fold cross-validation method, K = 3). The genetic algorithm uses the bidding model to select the optimal individual (the selection ratio is 0.2). After crossover and mutation operations, a new generation of population is formed. The number of populations in each generation is 150, and the number of iterations is 200. Table 3 shows the range of hyperparameter optimization.

SVR and genetic algorithm are implemented based on libsvm (https://www.csie.ntu.edu.tw/∼cjlin/libsvm/) and geatpy package (http://geatpy.com/index.php/quickstart/).

Figure 6 shows the fitness curves of TOC interpretation model (Figures 6A–D) and homogeneous regression interpretation model (Figure 6E) based on rock facies classification and regression under genetic algorithm optimization. Table 4 lists the optimal parameters of the above optimization process and their corresponding cross validation MSE. It can be found from Figure 6 and Table 4 that the RBF kernel function obtains higher cross-validation accuracy in both interpretation models, that is, the cross-validation MSE is the smallest, which is better than the other two kernel functions. The optimal cross MSE obtained by the RBF kernel function in the homogeneous regression interpretation model is 1.3. In the TOC interpretation model based on rock facies classification regression, the optimal cross MSE in the tuffaceous/clay interbedded shale data set is 2.18, and the optimal root mean square error of other data sets is less than 0.7.

In order to better evaluate the generalization ability and prediction accuracy of the above model, the above model is applied in the corresponding validation set, and the TOC in each validation set data is predicted respectively. Figure 7 is the projection plot of the predicted TOC value and the measured TOC value of different prediction models. Figure 7A is the projection plot of the prediction results of the uniform regression interpretation model in its response data set, and Figure 7B is the projection plot of the prediction results of the classification regression interpretation model based on lithofacies in their respective validation data sets. Figure 7A shows that the overall prediction effect of the uniform regression interpretation model is poor. As shown in Table 4, MSE is 3.23, RMSE is 1.80, and MAPE is 31.65. It is worth noting that in the interval of TOC>9% in Figure 7A, the predicted value of the prediction model seriously deviates from the true value, and the relative error of the predicted value of individual data points is far more than 25%. In comparison, the classification regression interpretation model based on rock facies shown in Figure 7B has better prediction performance. The distribution of data points is closer to the baseline of real TOC equal to predicted TOC. The MSE, RMSE and MAPE in the evaluation indexes are 0.91, 0.95 and 19.56, respectively. The prediction ability is greatly improved compared with the uniform regression model.

Through the difference of prediction performance of the two interpretation models in different rock facies, the reason of poor prediction performance of homogeneous prediction model can be further analyzed. Figure 7C shows the relationship between real and predicted TOC in different lithofacies under the homogeneous regression model. It can be seen from the figure that although the overall prediction performance is not good, the homogeneous regression model has higher prediction performance on AS and ASLS, and each evaluation index is at a lower level. The reason for the great decrease of prediction performance is ATLS and SS. The evaluation indexes MSE, RMSE and MAPE of ATLS prediction results are 11.73, 3.42 and 25.14 respectively. Considering the low TOC characteristics of SS, it may be more appropriate to use MAPE for evaluation, but MAPE = 93.23 is obviously beyond the acceptable range. The reason for the above characteristics can be attributed to the fact that the uniform regression model makes the learning algorithm pay more attention to the learning of data features in the high distribution density interval. The SVR algorithm learns more feature information from the 3% ∼ 8% interval shown in Figure 1C for training, so that the evaluation index MSE is minimized. In this case, the prediction model will extract features from the two lithofacies of AS and ASLS. Figure 5 shows that the characteristics of TOC and logging in different lithofacies have certain differences, which results in the decline of prediction performance in ATLS and SS.

The classification regression prediction model based on lithofacies solves the problem of poor prediction performance of single prediction model in ATLS and SS, making the prediction performance close to AS and ASLS. As shown in Table 4, the prediction indexes of ATLS and SS have been greatly improved (Figure 7D). The MSE, RMSE and MAPE of ASLS are 2.24.1.49 and 16.37 respectively, and SS is 0.20,0.37 and 37.84 respectively. At the same time, the prediction performance of AS and ASLE has also been improved. For example, the MSE, RMSE and MAPE of AS are reduced to 0.70,0.84 and 15.76, and ASLS is reduced to 0.49,0.70 and 15.20, respectively.

The above prediction results show that the classification regression interpretation based on lithofacies can obtain better TOC prediction accuracy in the case of data imbalance and multi-stratigraphic factors.

The basis for the application of TOC interpretation model with better prediction performance is lithofacies interpretability. In this study, the lithofacies delineated in the coring section of 8 wells and their corresponding logging response values are used as supervised data sets, and the lithofacies logging recognition model is established by XGBoost algorithm.

The supervised data set includes seven kinds of logging data and lithofacies labels such as AC, DEN, GR, U, CAL, Th and Rt. The logging data uses the original data, and the lithofacies are coded according to TUFof 1, ATLS of 2, AS of 3, ASLS of 4, SS of 5. The XGboost algorithm is based on the XGboost package (https://github.com/d mlc/xgboost). The base model type is gbtree, the learning task is multi: softmax, and the learning objective is mlogloss. In XGboost training, K-fold cross-validation is used to obtain the cross-validation mlogloss value to determine the optimal number of iterations of the tree in XGboost, where K = 7, the maximum number is 200, and the optimal number of iterations is output after 30 iterations without performance improvement. Hyperparameters such as learning rate, max _ depth, min _ child _ weight, gamma and sub _ sample are optimized by genetic algorithm. The fitness function is set to cross mlogloss. The optimization range and coding method are shown in table 5.

Figure 8A shows the optimal and average fitness curves in the genetic algorithm optimization process, where the optimal mlogloss is 0.019, and the optimal parameters shown in Table 5 are determined. The maximum number of times obtained by using this parameter is 179. Using the prediction model obtained by this parameter training, 556 data points in the coring sections of WY1 and FY3 were identified, and the confusion matrix shown in Figure 8B was drawn (Figure 8C).

The prediction accuracy of the model in the prediction set, the prediction accuracy and recall rate of each class are calculated by the confusion matrix (Table 6). It can be seen that the accuracy and recall rate of the four rock phases of TUF, ATLS, AS and SS are greater than 0.80, and the prediction accuracy and recall rate of ASLS are low, which can still reach about 0.75. The overall prediction accuracy is 0.86, which shows that the prediction model has high prediction accuracy and lays the foundation for the application of LBCRM model in actual drilling.

Based on lithofacies prediction model and LBCRM method, the TOC content of shale can be predicted. Firstly, the lithofacies prediction model is used to identify the lithofacies types of shale section, and the response TOC interpretation model is used to predict TOC for each lithofacies type. Figure 9 shows the lithofacies identification and TOC interpretation results of the shale depth section of 1,250 ∼ 1,340 m in Well FY3. For comparison, the figure also shows the TOC calculated by a single regression model. It can be seen from Figure 9 that the TOC interpretation results obtained by the LBCRM, TOC interpretation model are closer to the measured values. Especially at 1,277 ∼ 1,283 m, LBCRM successfully explained that in addition to the high TOC content in this depth section, the single regression model explained it as lower TOC.