Applying Statistical Analysis and Machine Learning for Modeling the UCS from P-Wave Velocity, Density and Porosity on Dry Travertine

In the rock mechanics and rock engineering field, the strength parameter considered to characterize the rock is the uniaxial compressive strength (UCS). It is usually determined in the laboratory through a few statistically representative numbers of specimens, with a recommended minimum of five. The UCS can also be estimated from rock index properties, such as the effective porosity, density, and P-wave velocity. In the case of a porous rock such as travertine, the random distribution of voids inside the test specimen (not detectable in the density-porosity test, but in the compressive strength test) causes large variations on the UCS value, which were found in the range of 62 MPa for this rock. This fact complicates a sufficiently accurate determination of experimental results, also affecting the estimations based on regression analyses. Aiming to solve this problem, statistical analysis, and machine learning models (artificial neural network) was developed to generate a reliable predictive model, through which the best results for a multiple regression model between uniaxial compressive strength (UCS), P-wave velocity and porosity were obtained.


Introduction
The uniaxial compressive strength (UCS) is one of the most important parameters used in the rock mechanics and rock engineering field for design. It is widely used to characterize rock [1] and to classify rock masses in a quantitative way [2][3][4]. There are several methods to determine it in the field and in the laboratory, mainly provided by the ISRM Suggested Methods [5,6] and the ASTM [7][8][9] Standard Test Methods.
The determination of UCS in the laboratory can be time-consuming and, in some cases, unaffordable, due to the specimen preparation often requires relatively high-quality cores in order to fulfill geometric and surface specifications, as recommended by the ISRM [5] or ASTM [10]. This may be one of the reasons behind the research of correlations between this parameter and others come from non-destructive methods like P-wave velocity (v p ). The v p was found to be mainly dependent on multiple factors, with several correlations developed for density (ρ) and porosity (n) [11][12][13][14][15][16]. The UCS depends mainly on ρ and n [17], and therefore on v p , so many researchers have determined different where m dry is the dry mass (g), m sat is the saturated mass (g), m sub is the submerged mass (g), V v is the pore (void) volume (cm 3 ), V T is the bulk volume (cm 3 ), and ρ w is the water density (g·cm −3 ).

Ultrasonic Pulse Transmission Tests
A PROCEQ Pundit Lab+ apparatus was used to measure primary (v p ) and secondary wave (v s ). velocities. The test specifications are those of the method suggested by the ISRM [6]: 54 kHz bandwidth and direct configuration for the transducer pair (transmitter-receiver). The force applied to the transducers involved a stress less than 10 kPa. The programmed voltage was 500 V and the gain was 1x. The measure of v p was determined automatically by the machine, and v s was determined manually when both P and S waves were accoupled in the range of 1.4 to 2.3 times the P-wave time arrival. For good transmission between the transducers, the polished surfaces of the specimens were smeared with stiffer grease.

Uniaxial Compressive Strength Test
The method suggested by ASTM [9,10] was followed, by using a fully servo-controlled 300 t press (CONTROLS 50-C52Z00 + MCC8 50-C8422/M), which allows uniaxial compression strength tests to be carried out on specimens up to 150 mm in diameter. The loading speed was 0.2 MPa·s −1 so that the test lasted between 5 and 10 min.

Conventional Statistical Models
The effect of independent variables in the modeling of uniaxial compressive strength was studied by means of an experimental design [30,31] generating multiple linear regression models and a quadratic model, adjusted to the experimental tests to study the strength of travertine to variations of v p , n and ρ dry .
The methodology used consisted of compression tests, and 29 experimental tests were carried out, studying the effects of velocity of propagation, saturated density and porosity on the dependent variable. For the modeling and experimental design of regression models, the Python "Scipy" library and "statsmodels" module was used (using the ordinary least squares method), making it possible to investigate the linear effects, interactions and quadratic effects of the independent variables in the UCS. The experimental data were fitted first to a simple (v p -dependent) regression Appl. Sci. 2020, 10, 4565 4 of 14 model, and then through a multiple regression analysis [32] (incorporating n and ρ dry ) to a quadratic model, considering only those factors that help explain the model variability and that have statistical significance (p < significance level).
The general form of the experimental model is presented in Equation (4): where x i are the independent variables (v p , ρ dry , n), and the answer Y corresponds to the dependent variable (UCS). The determination coefficient (R 2 ), p-values and F-value indicate whether the model obtained is adequate to describe the uniaxial compressive strength for the set of values sampled [33,34].

Artificial Neural Networks
ANNs are supervised machine learning techniques, which determine associations between a known set of observations (i.e., training points) and different environmental variables to classify new and unknown data (test set). ANNs can approximate non-linear relationships and generalize complex systems from relatively imprecise information, and are robust in handling noise, overfittings, and outliers [35].
In an ANN system, the nodes are connected by means of synapses, this connection structure determines the behavior of the network. The most used structure is the multilayer perceptron [36], as presented in Figure 1, where x i represents the inputs, oi the outputs and σ the activation function. The neurons in the input layer depend on the information available to be classified, which is given by the variables: compression speed, ρ sat and n, while the output layer is given by the uniaxial compression resistance. The general form of the experimental model is presented in Equation (4): (4) where xi are the independent variables (vp, ρdry, n), and the answer Y corresponds to the dependent variable (UCS). The determination coefficient (R 2 ), p-values and F value indicate whether the model obtained is adequate to describe the uniaxial compressive strength for the set of values sampled [33,34].

Artificial Neural Networks
ANNs are supervised machine learning techniques, which determine associations between a known set of observations (i.e., training points) and different environmental variables to classify new and unknown data (test set). ANNs can approximate non-linear relationships and generalize complex systems from relatively imprecise information, and are robust in handling noise, overfittings, and outliers [35].
In an ANN system, the nodes are connected by means of synapses, this connection structure determines the behavior of the network. The most used structure is the multilayer perceptron [36], as presented in Figure 1, where xi represents the inputs, oi the outputs and σ the activation function. The neurons in the input layer depend on the information available to be classified, which is given by the variables: compression speed, ρsat and n, while the output layer is given by the uniaxial compression resistance. The mathematical expression that describes the output of each neuron is given by the activation function σ, where ωi, are the synaptic weights that weight the xi inputs, ω0 is the threshold and n is the total number of synaptic weights connected to the input of the neuron. Each neuron has its respective threshold level, and the hyperbolic tangent was used as a trigger function on all neurons in the network [37]. The learning method used was backpropagation, which follows the following stages: it starts when a sample of the process is captured by the input layer, then the activation value of each of the neurons is calculated progressively through each neuron or node through all the layers, from the input layer to the output one. Then, in the check stage, the activation value of the output neurons is compared with the expected output data and, if this error differs too much, the error is corrected in a distributed way by updating the weights of the neural network in a return stage, from the output layer to the input one. The whole process is repeated until the network adjusts the synaptic weights and the error value of the outputs is permissible [35].
The implementation of neural network models was developed by generating a model using the Keras neural network library. Keras models are defined as a sequence of layers, for which a "Sequential" model is created, adding layers one after another, until the requirements are met. Then, due to the limited number of samples, the models to be generated are 1 hidden layer, with variations in the number of neurons, which are completely connected by using the dense class. For the hidden layer, the ʺreluʺ activation function will be defined, and for the output layer, a ʺsigmoidʺ function. The mathematical expression that describes the output of each neuron is given by the activation function σ, where ω i , are the synaptic weights that weight the x i inputs, ω 0 is the threshold and n is the total number of synaptic weights connected to the input of the neuron. Each neuron has its respective threshold level, and the hyperbolic tangent was used as a trigger function on all neurons in the network [37]. The learning method used was backpropagation, which follows the following stages: it starts when a sample of the process is captured by the input layer, then the activation value of each of the neurons is calculated progressively through each neuron or node through all the layers, from the input layer to the output one. Then, in the check stage, the activation value of the output neurons is compared with the expected output data and, if this error differs too much, the error is corrected in a distributed way by updating the weights of the neural network in a return stage, from the output layer to the input one. The whole process is repeated until the network adjusts the synaptic weights and the error value of the outputs is permissible [35].
The implementation of neural network models was developed by generating a model using the Keras neural network library. Keras models are defined as a sequence of layers, for which a "Sequential" model is created, adding layers one after another, until the requirements are met. Then, due to the limited number of samples, the models to be generated are 1 hidden layer, with variations in the number of neurons, which are completely connected by using the dense class. For the hidden layer, the "relu" activation function will be defined, and for the output layer, a "sigmoid" function. Each ANN will be trained for a total of 1000 epochs, the number of neurons in the hidden layer will be defined in the next section, while that median absolute deviation (MAD), mean squared error (MSE), coefficient of determination (R 2 ) and accuracy (ACC) were used to evaluate the fit of the modeled networks. Accuracy is a performance measure defined how the ratio of correctly predicted observation over the total observations.

Analysis and Discussion of Results
The results from index properties (v p , v s , ρ dry , ρ sat and n) and strength (UCS) are shown in Table 1.
The mean values of ρ dry , ρ sat and n were 2.43 g·cm −3 , 2.47 g·cm −3 and 4.05%, respectively. The n values range from 0.7 to 8.1%, a fairly wide interval, due to the heterogeneity in the pore distribution and the travertine petrofabrics [38,39], which is in line with the UCS values, whose average value was 91.72 MPa but in the 25-29 test specimens the UCS values became 70 MPa on average to 120 MPa, coinciding with this decrease in n. The increase in dry and saturated primary velocity ∆v p was 369 m·s −1 on average, from which a clay mineral content ECC practically zero, has been deduced [40]. The average dry v p /v s ratio is 1.84, within the range 1.8-2 obtained by Soete et al. [38] for travertine and 1.6-2.3 by Schön [41], for consolidated sediments.
The dry v p values were correlated with the UCS of the rocks using the simple and multivariate regression methods. The experimental data set was used to generate empirical equations developed to estimate the UCS of travertine from the independent variables v p , ρ dry and n. In simple regression analyses, models were generated to individually relate the independent variables (v p , ρ dry and n) to the dependent one (UCS) by means of simple linear regression models, exponential models, potentials and quadratics through a logarithmic transformation, while in multivariate regression analyses, all variables considered in the sampling that are significant, and that contribute to explain the variability of the response were included. The analytical model, the coefficient of determination (R 2 ) and the p-value were calculated, to study and compare the goodness of fit of the models generated for each regression.

UCS Prediction from Selected Bibliography
As indicated by González et al. [18], the existing simple correlations for limestone from the literature were not suitable for the 'Antofagasta' limestone, which implied the development of a multiple correlation in which more variables (v p , ρ sat and n) were introduced to increase the fit of the experimental equation, and to be able to characterize the rock with easily measurable properties. This is due to the fact that the UCS depends both on the index properties (v p , ρ dry , n, I s(50) , R N , etc.), and on the petrofrabrics and internal pore distribution, as can be seen in the values of Table 2. Likewise, the UCS depends on the size of the specimen [1,42], as well as the v p [43,44], parameters that must be considered for the scale effect, when applying it to an engineering project.
* The authors do not explicitly indicate saturated or dry conditions, so it is assumed as dry (S r = 0).
In the case of travertine, the same occurs. Table 3 shows the different linear correlations that relate UCS with v p , (Equations (5)-(9)), with n (Equations (10) and (11)) and ρ dry (Equation (12)). Equation (13) corresponds to the multiple correlation that includes the three properties. Figure 2 shows how no correlation fits the experimental results (line 1:1), and Table 3 shows the mean values and relative errors. From these average values, it can be concluded that Equations (5) and (11) are the closest to 101.09 and 105.27 MPa, respectively (they have the lowest relative errors in absolute value). It should be noted that in the determination of these equations no travertine specimens were considered, only limestone, marble, dolomitic limestone, dolomite, marble and graveled limestone.
* The authors do not explicitly indicate saturated or dry conditions, so it is assumed as dry (Sr = 0)

Correlation and Univariate Regression Analyses
Developing Pearsonʹs correlation analysis (Figure 3), the positive correlation between UCS and vp,dry is highlighted, and the negative correlation between UCS and n. It is also observed that the UCS dependent variable does not have a strong linear correlation with the other sampled variables, especially with the ρsat, vs,dry and vs,sat variables, where there is no linear correlation (r = 0).

Correlation and Univariate Regression Analyses
Developing Pearson's correlation analysis (Figure 3), the positive correlation between UCS and v p,dry is highlighted, and the negative correlation between UCS and n. It is also observed that the UCS dependent variable does not have a strong linear correlation with the other sampled variables, especially with the ρ sat , v s,dry and v s,sat variables, where there is no linear correlation (r = 0). With the purpose of studying the relationship between the independent variables, vp,dry (hereinafter referred to as vp) and n, and the dependent variable (UCS), a simple regression model (Equations (14)-(19)) was set up, to model the compression strength as a function of the independent variable (considering the exponential and potential models, after log transformation). Additionally, the curvature is considered in the analytical model, to study its impact on the response variable With the purpose of studying the relationship between the independent variables, v p,dry (hereinafter referred to as v p ) and n, and the dependent variable (UCS), a simple regression model (Equations (14)- (19)) was set up, to model the compression strength as a function of the independent variable (considering the exponential and potential models, after log transformation). Additionally, the curvature is considered in the analytical model, to study its impact on the response variable (Equations (20) and (21)). The mathematical models generated as a function of v p and n and the respective goodness-of-fit statistics (R 2 ) are presented in Table 4. From the analysis of the observed data distributions in Figure 4, and considering the goodness-of-fit statistics presented in Table 4, the best model to explain the UCS as a function of v p is the quadratic model (R 2 = 0.6804), while the model as a function of n that presents a better fit also is the quadratic model (R 2 = 0.7573). On the other hand, although the adjusted models presented in Table 4 yielded considerably high coefficients of determination, these equations only include one independent variable. Therefore, an analysis that considers a greater number of independent variables that affect the response variable jointly, and that help explain the variability of the response variable has the potential to improve the percentage of variability explained by the analytical model.

Multivariate Regression Analyses
By developing the multiple regression fitting and based on the information obtained from the ANOVA analysis, the first-and second-degree multiple regression models (Equations (22) and (23), respectively), shown in Table 5, are able to represent the model in question for the set of sampled values. In addition, the linear effects and interactions of all factors must contribute greatly to explaining the experimental model. The ANOVA test of the model indicates that the model presented in Equation (22) is adequate to represent the UCS under the range of established parameters, and although there is no lack of fit of the model, and the value of R 2 (0.6788) validates it (v p and ρ dry do not an effect in the response for the single multiple model), the effect of the other variables does not exceed the fit of the univariate models shown in Equations (18) and (20). Additionally, the p-value of Equation (22) indicates that the model is statistically significant.
On the other hand, the ANOVA test of the second degree model (Equation (23)) indicates that the model is also adequate to represent the response variable as a function of the independent variables, and the goodness-of-fit indicator R 2 (0.8447) validates this. The predicted and experimental UCS values are compared in Figure 5, where it is shown that the graphs of observed vs. predicted data generally follow the 1:1-line and the p value of the normality test of residuals are greater than the level of significance (value p > 0.05), indicating a good fit of the models presented in the Table 5.
to represent the UCS under the range of established parameters, and although there is no lack of fit of the model, and the value of R 2 (0.6788) validates it (vp and ρdry do not an effect in the response for the single multiple model), the effect of the other variables does not exceed the fit of the univariate models shown in Equations (18) and (20). Additionally, the p-value of Equation (22) indicates that the model is statistically significant.
On the other hand, the ANOVA test of the second degree model (Equation (23)) indicates that the model is also adequate to represent the response variable as a function of the independent variables, and the goodness-of-fit indicator R 2 (0.8447) validates this. The predicted and experimental UCS values are compared in Figure 5, where it is shown that the graphs of observed vs. predicted data generally follow the 1:1-line and the p value of the normality test of residuals are greater than the level of significance (value p >> 0.05), indicating a good fit of the models presented in the Table 5.

Fitting of Artificial Neural Networks
With the aim of studying the fitting of the UCS by using models based on ANN's, different architectures were generated, and their fitting was studied using the indicators of goodness of fit indicated above. Three artificial neural networks architectures were considered for training and Figure 5. Predicted UCS vs. experimentally determined UCS for simple multiple regression (a) (Equation (22)) and second degree multiple regression (b) (Equation (23)).

Fitting of Artificial Neural Networks
With the aim of studying the fitting of the UCS by using models based on ANN's, different architectures were generated, and their fitting was studied using the indicators of goodness of fit indicated above. Three artificial neural networks architectures were considered for training and testing, with network architectures of 1-2, 1-3 and 1-4 hidden layer and neurons per hidden layer, respectively. The statistics of the percentage of relative errors between the predicted responses and the measured values have been summarized in Table 6. From the analysis of the goodness-of-fit statistics of the different ANN configurations presented in Table 6, it can be checked that the model presenting the best fit in is the simplest configuration (one hidden layer with two neurons), which is validated by the highest determination coefficient (0.6975), a high accuracy (85%) and low errors (MAD and MSE). The increase of the coefficient of determination through the training and testing periods is presented in Figure 6.  Table 6. From the analysis of the goodness-of-fit statistics of the different ANN configurations presented in Table 6, it can be checked that the model presenting the best fit in is the simplest configuration (one hidden layer with two neurons), which is validated by the highest determination coefficient (0.6975), a high accuracy (85%) and low errors (MAD and MSE). The increase of the coefficient of determination through the training and testing periods is presented in Figure 6.  Table 7 presents the results of both the regression analyses and the best ANN architecture. Regression models and their fit are presented, relating UCS (dependent variable) first with vp or n (linear, exponential, potential and quadratic models) and then including all rock parameters for the travertine stone studied (multiple regression and neural network models). Simple regression  Table 7 presents the results of both the regression analyses and the best ANN architecture. Regression models and their fit are presented, relating UCS (dependent variable) first with v p or n (linear, exponential, potential and quadratic models) and then including all rock parameters for the travertine stone studied (multiple regression and neural network models). Simple regression analyses found that UCS is related to v p and n, and the best fit was obtained with a quadratic model (R 2 = 0. 7573), dependent only on the porosity of the rock, while the neural network model is the one that presents the best fit for multivariate modeling, considering the totality of adjusted models. Neural network analyses are competitive with single and multiple linear regression analyses, but not with the quadratic multiple model (model that presents the best fit). However, due to the lack of sample data, the adjustment potential of the tool and its ability to abstract the modeled system is not evident.

Conclusions and Future Works
The combined effect of v p , ρ dry and n in UCS for travertine-saturated cylindrical specimens was studied. It is evident that v p and n have a synergistic effect in UCS (both individually and combined), while ρ dry does not have a statistically significant influence (linear, exponential, potential, and quadratic) on the dependent variable. This is explained by the fact that, for the same rock type, the density remains relatively constant, since the measured porosity only affects the external pores of the test specimen, which do not necessarily have to be interconnected with the internal ones. Therefore, v p may vary, but density does not, and, therefore, there is no correlation between them.
The experimental results were analyzed through simple regression analysis (or in its absence, after log transformation) to develop the analytical models that best adapt to UCS depending on v p (as developed by most authors). Subsequently, multiple regression analysis and an artificial neural network model (ANN) were performed, to obtain a predictive model of UCS from v p , n and ρ dry .
Simple regression analyses (using only v p or n as independent variable) showed a better R 2 value of 0.7573 for the univariate model, dependent only on porosity, while the univariate model dependent on v p turned out to be the quadratic adjustment, with an R 2 of 0.6804. On the other hand, a multiple regression analysis, considering interactions between factors and quadratic terms resulted in a better R 2 value (0.8447), which indicates a better adjustment when additional explanatory variables are considered. Multiple regression analyses were also able to identify ρ dry as a statistically non-significant independent variable for predicting UCS from ultrasonic pulse transmission tests. Finally, the ANN-based model presented a competitive fit with the single and multiple linear regression models in predicting UCS (R 2 = 0.6975), showing that UCS is directly dependent on the independent variables v p , ρ dry and n.
To continue this research line with travertine, the incorporation of a greater number of independent variables in the analysis and a larger number of samples is planned. This research aims to identify the variables that explain the experimental variability of UCS under sampling conditions, and to model the dynamics of the system using machine learning algorithms, such as neural networks or Bayesian networks [46]. The possible variables to be considered are the corrected point load index (I s(50) ), and the number of Schmidt's hammer rebounds (R N ). Further tests will be carried out to try to improve the analytical models.