An adaptive prediction method for mechanical properties deterioration of sandstone under freeze–thaw cycles: a case study of Yungang Grottoes

Due to the location of the Yungang Grottoes, freeze–thaw cycles contribute significantly to the degradation of the mechanical properties of the sandstone. The factors influencing the freeze–thaw cycle are classified into two categories: external environmental conditions and the inherent properties of the rock itself. Since the parameters of rock properties are inherent to each rock, the effect of rock properties on freeze–thaw degradation cannot be investigated by the control variates method. An adaptive multi-output gradient boosting decision trees (AMGBDT) algorithm is proposed to fit nonlinear relationships between mechanical properties and physical factors. The hyperparameters in the GBDT algorithm are set as variables, and the Sequential quadratic programming (SQP) algorithm is applied to solve the hyperparameter optimization, which means finding the maximum Score. The case study illustrates that the AMGBDT algorithm can precisely determine the effect of each independent factor on the output. The patterns of mechanical properties are similar when the number of freeze–thaw cycles and porosity are used as variables separately and when both are used simultaneously. The uniaxial compressive strength decay rate is positively correlated with the number of freeze–thaw cycles and porosity. The modulus of elasticity is negatively correlated with the number of freeze–thaw cycles and porosity. The results show that the number of freeze–thaw cycles is the main factor influencing the freeze–thaw cycling action, and the porosity is minor. In addition, the fitting accuracy of the AMGBDT algorithm is generally higher than neural networks (NN) and random forests (RF). Studying the influence of porosity and other rock properties on the freeze–thaw cycle will help to understand the failure mechanism of rock freeze–thaw cycles.


Introduction
The Yungang Grottoes, one of the three major grotto groups in China and world-famous stone sculpture art, was inscribed on the UNESCO World Heritage List in 2001. Therefore, Yungang Grottoes has a significant artistic and cultural value. Yungang Grottoes (113°20′ E, 40°04′ N) is located in a continental monsoon semi-arid climate. The annual average temperature in the study area is 7-10 °C, the monthly average temperature in January is − 11.4 °C, and the monthly average temperature in July is 23.1 °C. The freezing period is up to six months, with a standard frozen depth of 1.5 m.
The major stratigraphy of Yungang Grottoes is the Mesozoic Jurassic stratigraphy, mainly composed of fluvial and lacustrine deposits [1]. Yungang Grottoes Open Access *Correspondence: liuchenchen19@mails.ucas.ac.cn 1 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China Full list of author information is available at the end of the article is distributed in a vast lenticular miscellaneous sandstone. The lithology of the sandstone is mainly calcareous cemented sandstone, and gravel is often visible at the bottom of the sandstone's lenticular body. The main components of the sand are quartz, feldspar, and rock debris. The main components of the conglomerate are quartz breccia, rock clasts, and a small portion of mudstone clumps. Due to the geographical location of the Yungang Grottoes and the environment, the effects of weathering caused by freeze-thaw cycles cannot be ignored.
A number of experiments to investigate the effect of freeze-thaw cycling tests on the mechanical properties of sandstone were conducted [2,3], such as porosity [4], the number of freeze-thaw cycles [5][6][7]. The water content of the rock and the number of freeze-thaw cycles affect the degree of deterioration by freeze-thaw cycles. The main reasons are as follows. (1) The main cause of freeze-thaw damage in rocks is the repeated generation and dissipation [8] of frost heaving forces. The number of freeze-thaw cycles defines the number of times the frost heaving forces are repeatedly generated and dissipated, which means the more times, the more serious the freeze-thaw damage. (2) Pore water influences rock damage through three mechanisms: phase swelling, hydrostatic pressure, and formation of ice lenses or wedges [9], phase swelling, hydrostatic pressure, and formation of ice lenses or wedges. The mechanisms indicate a positive correlation between the pore water content of rocks and damage to saturated rocks. The larger the initial porosity of the rock, the greater the water content at saturation, resulting in greater tensile stresses in the pore walls caused by volume expansion due to water-ice phase changes. The initial porosity affects the freeze-thaw deterioration of rocks in part.
However, most experiments [10,11] focused on the effect of external environmental conditions on the freeze-thaw deterioration of sandstone. The mechanism of freeze-thaw damage in rocks shows that the initial porosity of rocks also affects freeze-thaw damage. It is impossible to experimentally achieve the control of variables such as intrinsic rock parameters, which means that it is hard to explore the degree of influence of each variable when multiple variables act. In this paper, with the help of the developed AMGBDT algorithm, we investigate the changes in mechanical properties of sandstone when initial porosity and the number of freeze-thaw cycles are adopted as variables.
The selection of hyperparameters has a significant impact on the fitting performance of GBDT. In this paper, the AMGBDT algorithm is proposed to complete the GBDT hyperparameter tuning task, and a prediction model is given based on this algorithm. The prediction model is based on a series of freeze-thaw cycling experimental data to quantify the extent of porosity effects on rock deterioration under freeze-thaw cycles. The prominent advantage of this method is that the degree of contribution of the rock's inherent parameters to the deterioration can be calculated analytically. This paper provides a considerably more efficient and costeffective means of studying the failure mechanism of freeze-thaw cycles.
The remainder of this paper is organized as follows. The methodology section introduces the principles of the GBDT algorithm and the SQP algorithm and the construction framework of the AMGBDT-SQP algorithm. The comparison between the fitted curves of experimental data and the AMGBDT algorithm is introduced in the case study when the number of freeze-thaw cycles and the porosity are applied as independent variables separately. The discussion section discusses three methods to demonstrate the GBDT algorithm prediction accuracy by the algorithm Score's evaluation index. The fifth part outlines conclusions.

Methodology
High-accuracy predictions can be obtained by these nonparametric machine learning (ML) methods, including classification and regression trees (CART) [12], support vector machines (SVM) [13][14][15], and k-nearest neighbor algorithm (KNN) [16,17]. The above methods are supervised learning methods, and target variables need to be prepared for the dataset [18]. The potential relationships in the experimental data set can be captured effectively by these ML methods, with good predictive performance but hard to interpret [19].
Ensemble learning has been a popular machine learning method in recent years. It refers to a classification system that uses the idea of bagging [20], stacking [21], or boosting [22] to build and combine multiple weak base learners to complete the learning task. It can obtain superior generalization performance than a single learner [23]. Compared with other ML methods, GBDT can interpret the interactions between input variables and predictive models and also identify the relative importance of key factors through ensemble learning [24]. The algorithms combining weak learners with decision trees come in two relatively primary forms: RF and GBDT.
The GBDT algorithm differs significantly from the traditional boosting algorithm in that each computation is designed to reduce the residuals from the previous one. The model is built in the direction of the decreasing gradient of the residuals to eliminate the residuals. The GBDT algorithm is a non-parametric model. The number of parameters is uncertain before training, so the fully grown decision tree has a large degree of freedom, resulting in fitting the training data to the maximum extent. The gradient boosting method for integrating multiple decision trees can well solve the problem of overfitting.
Tree-based integration methods have been widely used in the field of geotechnical engineering in recent years. Due to the successful applications of the GBDT algorithm in many research fields [25][26][27], we argue that the GBDT algorithm will be reasonably competent for the work of discovering the relationship between the variables and responses.
To construct a hybrid model that accurately predicts the mechanical property deterioration of sandstone after freeze-thaw cycles, the AMGBDT method is proposed by combing the GBDT algorithm and the SQP algorithm into one framework. The relationship between independent variables and the response is learned in the proposed model by applying the GBDT algorithm as a regressor. The wide range of hyperparameter variations in the GBDT algorithm results in difficult manual adjustment of the hyperparameters. Therefore, the SQP algorithm is applied as an optimizer to find the optimal parameters of the GBDT algorithm, increasing the model's adaptability.

GBDT algorithm
The GBDT algorithm is an algorithm that classifies or regresses data by using an additive model (i.e., a linear combination of basis functions) and continuously reducing the residuals generated by the training process. The bottom layer of the algorithm is mainly based on regression trees and gradient descent algorithms in function space. The algorithm not only possesses the advantages of solid interpretability of the tree model, effective handling of mixed types of features, scale-invariance (no need to normalize the data), and robustness to missing values but also has the advantages of reliable predictive power and good stability.
Compared to its successor XGB (Extreme gradient boosting) and LGB (Light gradient boosting machine) algorithm, the GBDT algorithm only requires the loss function to be first-order derivable, and both convex and non-convex functions are applicable. However, the XGB/ LGB algorithm requires a stricter loss function, first-and second-order derivable, and a strictly convex function. Therefore, the GBDT algorithm was adopted to explore the effects of porosity and the number of freeze-thaw cycles on the deterioration of rock mechanical properties.
Data description. Before building the GBDT model, the freeze-thaw cycle experimental data are given as follows: where the ftt i refers to freeze-thaw cycle times of the i th sample; p i refers to the initial porosity of the i th sample; v i represents the longitudinal wave velocity of the i th sample;sl i is the loss rate of uniaxial compressive strength,e i refers to the elastic modulus after the deterioration of the i th sample,ml i refers to the loss rate of mass, N is the number of samples.
After inputting the training dataset, the GBDT algorithm first constructs an initial regressor which can be defined as follows: When initializing the base regressor, c adopts the value of the mean of all training sample labels.
Following the construction of an initial decision tree, the GBDT algorithm performs multiple iterations, where each iteration produces a new base regressor. Each base regressor is fitted with the previous regressor's residuals to minimize the current round's residuals.
The final total regressor is obtained by weighting and summing base regressors obtained from each round of training. The framework of the AMGBDT-SQP algorithm is shown in Fig. 1.
Aiming at the problem of the loss function fitting method, Freidman [28] proposed to fit an approximation of the loss of the current round with the negative gradient of the loss function, which in turn fits a CART regression tree. The specific implementation flow of the algorithm is as follows: For the number of iterations t = 1, 2, …, T.
(1) For i = 1, 2, …, N , calculate the negative gradient of the loss function in this iteration.
where L(•) is the loss function and is constructed by the least-squares method in this paper.
(2) The t th regression tree can be obtained by fitting a CART regression tree to the residuals generated at Where J is the number of leaf nodes of t th regression tree. (3) For the leaf node region j = 1, 2, …, J , calculate the best-fit value for each region.
Finally, a strong regressor is obtained.

SQP algorithm
The adaptive algorithm allows the model to automatically select the optimal parameters for any reasonable data set, resulting in the best prediction of the model, thus enabling the model to be more adaptive. The objective functions are all nonlinear functions of the independent variables when the hyperparameters in the model are the independent variables and need to be solved using nonlinear programming algorithms. SQP algorithm is one of the most effective methods in current algorithms for solving small and medium-scale nonlinear optimization problems. Its main idea is to utilize a number of quadratic programming to (6) sequentially approximate the original nonlinear programming problem [29,30]. The number of trees (M) and the learning rate of a model (µ) are the hyperparameters that significantly impact the performance of the GBDT model. Therefore, the M and the µ are considered as the independent variables of the SQP algorithm. The number of trees M refers to the number of iterations or base models in the additive model. Underfitting or overfitting is prone to occur when the value of M is too large or too small. The learning rate µ refers to the weight reduction coefficient of each CART regressor tree, and the value is greater than 0 and less than 1. Generally speaking, the smaller the µ is, the more trees are needed to fit a model. M is inversely related to µ. Therefore, M and μ should be adjusted simultaneously to optimize the model prediction.
The M and μ parameters are interconnected, causing the optimization problem a combined, constrained optimization problem. The mathematical model for parameter optimization is defined as follows.
where f (Z) is the coefficient of determination Score.

Fig. 1 Flowchart of the proposed work
The objective function of the nonlinearly constrained problem at the current iteration point Z k is reduced to a quadratic function on the variable S via Taylor expansion.
Adopt the optimum solution S * as the next search direction S k of the current problem, and use a onedimensional search along the above direction to obtain Z k+1 . An approximate solution to the original problem.
SQP algorithm has an excellent performance in convergence, computational efficiency, and boundary search. Besides, it is based on sound mathematical theory support. However, the relationship between the objective function and the independent variables is relatively complex for model parameter optimization problems. The objective function may have more than one extreme value point, rendering it easy to fall into a local optimum solution applying this algorithm. A large number of initial (11) points are generated and computed to find local optimal solutions to obtain the highest global Score. The local relative optimal solution can be approximated as the global optimal solution with a sufficiently large number of initial points.

Data description
The data used in this paper is from the existing test data recorded in the subproject III of Yungang Grottoes conservation project [1]. The sampling locations for the freeze-thaw cycle experiments and the photograph of part of the specimens are shown in Fig. 2. This section investigates the effect of freeze-thaw cycles and initial porosity on mechanical property deterioration. The statistical properties of the experimental data are listed in Table 1. The standard deviation is adopted to quantify the extent of the variation of data. The formula for standard deviation is as follows: where σ is the symbol for standard deviation and the X k is the symbol for the mean value.

Preprocessing of the experiment data
The number of freeze-thaw cycles and initial porosity vary in values and magnitudes. The number of freezethaw cycles varies from 10 to 35, while the range of initial porosity is from 3.32% to 8.2%. The range of values for the number of freeze-thaw cycles is more extensive than the initial porosity, making the contour of the loss function steep. When seeking the optimal solution using gradient descent, it is quite possible to take a Zig-Zag route (vertical contour line), as Fig. 3a shows, resulting in many iterations to converge. Figure 3b shows that after normalizing the data, the optimal solution finding process becomes significantly smoother and easier to converge to the optimal solution in the right way. Therefore, to eliminate the differences in magnitude and range of values between different features, the maximum minimization method is used for data (13) normalization in this paper. The specific formula for the maximum-minimum normalization process is defined as The normalization process allows the preprocessed data to be limited to [0,1], eliminating the undesirable effects caused by differences in magnitude and range and speeding up the training speed.

Prediction results of the AMGBDT algorithm
The normalized training set data are input into the GBDT model, and the SQP method is applied to optimize the GBDT parameters. The optimal parameters are obtained as M = 353.60919; µ = 0.87306678 . The optimal parameters are input into the GBDT model to obtain the prediction model of mechanical property deterioration of sandstone.

The influence of the number of freeze-thaw cycles
In this experiment, the results of the uniaxial compressive strength test of rock samples after freeze-thaw cycles were analyzed to determine a correlation between the number of freeze-thaw cycles and the decay of mechanical indexes of rock masses. The experimental data of the same group of samples were averaged to reduce the error caused by the dispersion of the experimental data. The experimental data were fitted to obtain fitted curves regarding the number of freeze-thaw cycles and mechanical indices.
(a) The decay rate of uniaxial compressive strength.
It is seen from Fig. 4a that the predicted curves of the model generally match with the fitted curve However, the threshold range of uniaxial compressive strength in the prediction curve of the GBDT model is slightly smaller than the threshold range of the test curve. It can be seen that, among the independent variables, the number of freeze-thaw cycles significantly influences the uniaxial compression strength. (b) The elastic modulus. Similarly, the two predicted curves have the same pattern of variation, with the elasticity modulus showing a decreasing trend with the increase of the number of freeze-thaw cycles. However, the magnitude of the change in the AMGBDT model fitted curve is slightly smaller than the experimental. In addition, the experimental curve fit's effectiveness is reduced due to the discrete nature of the experimental data points. The robustness of the AMGBDT model to missing values makes the dispersion of data points have less impact on the fitting effect.
The variation law of the fitted curve of the GBDT algorithm is basically the same as that of the fitted curve of the test. However, since the number of freeze-thaw cycles is set as the independent variable in the AMGBDT model only, the law of the number of freeze-thaw cycles and the decay of mechanical indexes obtained by the AMGBDT model is more accurate.

The influence of the initial porosity
The initial porosity is an inherent parameter of the rock, whose extent of influence in freeze-thaw cycles cannot be determined by experimental means. In other words, it is hard to obtain the experimental fitted curve of initial porosity and deterioration of mechanical properties. In the AMGBDT model, the number of freeze-thaw cycles as a fixed value, the strength decay rate is positively correlated with initial porosity, and the modulus of elasticity is negatively correlated with initial porosity. According to the magnitude of the curve, the deterioration caused by initial porosity as a variable is less than that caused by the number of freeze-thaw cycles. It was demonstrated that the number of freeze-thaw cycles had a more significant effect on the uniaxial compressive strength decay rate and elastic modulus than initial porosity. Figure 5 provides the exact changes in elastic modulus and uniaxial compressive strength of the sandstone after freeze-thaw cycles when the initial porosity is used as the independent variable.

Discussion
Different modeling methods can result in different impacts in mining the nonlinear mapping relationship between the independent vector X and the dependent vector Y. There are many studies on the application of GBDT in geological disciplines. In addition, RF, NN performs well in fitting highly nonlinear relationships. Therefore NN, RF, and the AMGBDT model are selected for fitting prediction. The mutual constraints and interactions between neurons in the neural network algorithm enable the entire network to exhibit a nonlinear mapping from the input state space to the output state space. Therefore, the NN algorithm is suitable for solving problems with complex environmental information and unclear inference rules. The structure of the deep neural network used is shown in Fig. 6. The number of units in the input and output layers is 3 and 3, respectively, and the number of layers in the hidden layer is 3. The connection between layers is  A random forest is composed of many decision trees in a randomly selected subset of training data [31], where each tree in the forest is voted on and decided jointly at decision time [32]. It performs well on a number of datasets and can handle high-dimensional data [33]. In addition, RF does not select features and can be computed in parallel to speed up training. RF exhibits some fault tolerance for training data and is an effective method for estimating missing data. Moreover, RF suits multi-classification problems and detects the interactions between variable features and the degree of importance. The number of trees in the random forest is 100.
Furthermore, the AMGBDT model's Score calculated by Eq. (10) is chosen as the evaluation index. The Score of the three models is shown in Fig. 7. A higher value of Score indicates a better prediction. The comparison in Fig. 7 illustrates that with the number of freeze-thaw cycles and initial porosity as independent variables, the Score value of the AMGBDT model is closer to 1 than the remaining two algorithms, which means that the AMGBDT fits best. Therefore, the AMGBDT algorithm with multiple outputs has been proven to be the most accurate.

Conclusions
To study the effects of freeze-thaw cycles and initial porosity changes on the deterioration of mechanical properties of sandstone under freeze-thaw cycles, an AMGBDT model was developed for the study in this paper. Some limitations should be noted. The initial porosity is adopted as the independent variable when we explore the extent of influence of pore space inside the rock in the freeze-thaw damage. Due to the low cost and accessibility of the measurement of connected porosity and the considerable contribution in the freeze-thaw cycle deterioration, the initial porosity in this paper refers to the connected porosity. The frost damage mechanism of sandstone [34] demonstrates that both closed pores and connected pores contribute to rock freeze-thaw damage. The role of connected pores during freeze-thaw cycles may be amplified to some extent because no data on closed pores were measured in our study.
The conclusions are as follows.
(1) Based on the theory of SQP and GBDT, the AMG-BDT algorithm is proposed to allow its application to explore the effects of freeze-thaw cycling on the macroscopic deterioration patterns of sandstones in various locations. (2) The comparison of the prediction curves of the AMGBDT model with the curves fitted to the experimental data yields the following conclusions. When the initial porosity is fixed and the number of freeze-thaw cycles is the variable, the predicted curve of the AMGBDT model follows the same trend as the fitted curve of the test. However, the change in the decay rate of strength and elasticity modulus is less than the change range corresponding to the experimental curve. When the number of freeze-thaw cycles is fixed, and initial porosity is the variable, the changes in strength decay rate and modulus of elasticity are small and positively correlated. In summary, for dense rocks like sandstone, both initial porosity and the number of freeze-thaw cycles exacerbate the deterioration of mechanical properties, and the number of freeze-thaw cycles is the dominant factor. (3) The NN and RF are selected to perform the same training prediction. Comparing the Score of three models reveals that the model built by the AMG-BDT algorithm performs optimally, which illustrates the feasibility of the proposed method.
The model constructed in this paper can be further applied to the multi-field coupling analysis of rock weathering to obtain the degree of contribution of the influencing factors of each action field to the deterioration of the final mechanical properties. The influence ranking weights of each factor at the time of coupling are obtained, and the mechanism of multi-field coupling can be further studied.