For sulfur-containing oil products, the safety problems in production and processing are very serious. The FeS component contained in sulfur-containing oil products often reacts with air during storage and transportation, causing chemical reactions that lead to the release of a large amount of heat and cause many fire and explosion accidents [1-2].

To avoid the occurrence ofsuch malignant events, much research has been conducted in two aspects. On the one hand, many scholars have tried to analyze the reaction mechanism of spontaneous accidents and general reasons of producing T_max through chemical reaction theories. They proposed that the anti-corrosion coating on the inner wall of the storage tank would be partially peeled off and the iron in the coating would be exposed after the long service cycle. When the sulfur in the crude oil was exposed to a low temperature water-containing environment, the sulfurized rust was generated and the production facilities and tank walls was corroded through chemical and electrochemical modes [3-5]. Furthermore, the sulfurized rust underwent the oxidation reaction with air and released a large amount of heat to incur spontaneous combustion accidents [6-8]. Dou et al. [9-10] studied the chemical reaction mechanism of the self-heating process by analyzing the thermal decomposition kinetics of the sulfurized rust and obtained the thermal decomposition kinetic characteristic parameters and spontaneous combustion mechanism in different stages of the oxidation reaction of the sulfurized rust. With the same conditions in edible oil refining, Landucci [11-13] research ededible oil refining hazards using a thermodynamic model for the estimation of vapor phase composition in storage tanks as a function of operating condition, and this model provides a quick tool for preliminary assessment of hazards due to the formation of flammable mixtures in edible oil storage plants. In addition, the simulation experiments of the self-heating process of sulfurized rust under some special condition shave been conducted to explore the impacts of different factors. At the same time, many researchers adopted experiments to explore the implicit relationship between T_max and influence factors in the self-heating process. It was obtained from experiments that some factors, including the mass of sulfurized rust, operating temperature, water content of sulfurized rust, air flow rate through the safety valve, oxygen concentration in the safety valve, and pH, had an extremely important impact on the oxidation reaction of sulfurized rust [14-15]. Additionally, the researchers found that 70 ${}^{\text{o}}\text{C}$ is not only the critical temperature at which the chemical reaction takes place but also the temperature of the smoke point, at which SO₂is formed. Therefore, 70 ${}^{\text{o}}\text{C}$ is deemed as the threshold temperature for the self-heating process, and when the temperature exceeds 70${}^{\text{o}}\text{C}$, the fire or explosion accidents are triggered [16-18]. However, it is difficult to obtain ${{T}_{\max }}$ of the oxidation self-heating process through experiments because of the lack of accurate or theoretical relation ships between ${{T}_{\max }}$ and multiple uncertainties, namely influence factors. Thus, it is not realistic to quantify the detection and early warning for fire and explosion accidents [19-20].

On the other hand, the modeling method of machine learning was gradually applied to predict ${{T}_{\max }}$ of the oxidation self-heating process based on the relationship between ${{T}_{\max }}$ and the five main factors mentioned above. The Support Vector Machine (SVM) method was first used to build the ${{T}_{\max }}$ prediction model describing the relationship between the rising temperature and various uncertainties in the oxidation process of sulfurized rust based on the experimental dataset in literature [10]. ${{T}_{\max }}$ was obtained from the model based on the SVM method. Compared with the experimental results, the predictive ${{T}_{\max }}$ are within the error ranges. This shows that it is advantageous to study the relationship between ${{T}_{\max }}$ and its influence factors using the machine learning algorithm. In this paper, the Gaussian process regression(GPR) method is adopted to build the ${{T}_{\max }}$ prediction model. Compared with the SVM method, the GPR method has better adaptability to complex nonlinear problems with uncertainties and can make a probabilistic interpretation on the prediction result. Furthermore,the GPR method can meet the accuracy requirements needed in engineering applications and has achieved many successful cases. For example, Samulesson et al. used the GPR method to successfully address the monitoring and fault detection of waste water treatment processes [21]. The GPR method was also successfully applied in short-term wind speed forecasting, solar power forecasting, state-of-charge estimation for batteries, and fault detection [22-25], demonstrating its ability and application prospects to solve practical engineering problems.

In this paper, the GPR method is adopted to model the relationship between ${{T}_{\max }}$ and its influence factors during the oxidation self-heating process of sulfurized rust, and ${{T}_{\max }}$ is predicted using the built model. It is more helpful to solve the nonlinear ${{T}_{\max }}$ prediction problem in the oxidation self-heating process of sulfurized rust, and this will also provide stronger protection for oil and gas production and transport safety.

The remainder of this paper is organized as follows. Section 2 presents our methodology, including data preprocessing, the GPR method, and the model evaluation standard. Section 3 shares a case study on the oxidation self-heating process of the sulfurized rust based on the dataset [10]. Finally, Section 4 provides the conclusion.

To explore how ${{T}_{\max }}$ occurs in a self-heating process of sulfurized rust, regression methods are adopted. For this purpose, a set of experiments have been conducted at Nanjing Tech University. In the experiments, the five main influence factors mentioned above and ${{T}_{\max }}$ were monitored, and a dataset called set 1, which in cluded 85 samples, was obtained. Our goal is to find the implicit relation ship, as shown in Equation (1) [10], between ${{T}_{\max }}$ and the five factors: water content, mass of sulfurized rust, operating temperature in oil tank, air flow rate, and oxygen concentration in the safety valve.

Where ${{F}_{\text{wat}}}$ [unit:% as mass fraction]=water content of sulfurized rust, ${{M}_{\text{sr}}}$ [unit: g]=mass of sulfurized rust, ${{T}_{\text{ope}}}$ [unit: ${}^{\text{o}}\text{C}$]=operating temperature in oil tank, ${{R}_{\text{air}}}$ [unit: ml/min]=air flow rate through the safety value, and ${{C}_{\text{oxy}}}$[unit: % as volume fraction]=oxygen concentration in the safety value.

In Equation (1), the five factors can be considered as in dependent variables, and ${{T}_{\max }}$ is the dependent variable. The five factors are also considered as the five features to the GPR model, and the data in set 1 must be preprocessed before the GPR model is built.

During the training process of the GPR model to predict ${{T}_{\max }}$, in order to eliminate the in fluence of the different value ranges among five features on the prediction results, the Z-score normalization method is adopted to standardize the data involved in this paper. The standardized data will improve the training efficiency and prediction accuracy of the GPR model. For example, the feature ${{F}_{\text{wat}}}$ will be standardized through Equation (2) [26]:

Definition: ${{\mu }_{\text{wat}}}$ and ${{\delta }_{\text{wat}}}$ are the mean and standard deviation of feature ${{F}_{\text{wat}}}$ respectively. ${{F}_{\text{wat,}i}}$ is the i^th value of the feature ${{F}_{\text{wat}}}$ in the dataset involved in this paper. The other four features ${{M}_{\text{sr}}},\text{ }{{T}_{\text{ope}}},\text{ }{{R}_{\text{air}}},\text{ }{{C}_{\text{oxy}}}$, and ${{T}_{\max }}$ will conduct the same data pre-processing as ${{F}_{\text{wat}}}.$

The 85 samples in set1 are divided into the training set and the testing set, which are used to train and test the GPR model respectively. For obtaining a good GPR model, it is needed that all the minimum and maximum values of the five features and ${{T}_{\max }}$ in the dataset should be included in the training set. In addition, in order to prevent the GPR model from overfitting, the number of the samples included in the training set is supposed to be less than 80% of the total set [27-28]. Therefore, 60 samples in set1 (nearly 70%) are chosen to train the GPR model, and the remaining 25 samples are used to test the GPR model for the prediction of ${{T}_{\max }}$ in this paper. The additional 17 samples in literature [10], called set2,are used to further verify the generalization and effectiveness of the GPR model. We can see the change of ${{T}_{\max }}$ in different conditions from Figure 1.

Fora regression problem, it is usually assumed that there is a training set D with n observations expressed as

$D=\left\{ \left( {{\mathbf{x}}_{1}},{{y}_{1}} \right),\left( {{\mathbf{x}}_{2}},{{y}_{2}} \right),\cdots ,\left( {{\mathbf{x}}_{n}},{{y}_{n}} \right) \right\}$

Let ${{\mathbf{x}}_{i}}=\left[ {{x}_{i1}},{{x}_{i2}},\cdots ,{{x}_{id}} \right]$ and $i=1,2,\cdots ,n$, $\mathbf{X}={{\left[ {{\mathbf{x}}_{1}},{{\mathbf{x}}_{2}},\cdots ,{{\mathbf{x}}_{n}} \right]}^{T}}$, $\mathbf{y}={{\left[ {{y}_{1}},{{y}_{2}},\cdots ,{{y}_{n}} \right]}^{T}}.$ The task of the regression is to learn the mapping relationship $\left( f\left( \bullet \right)|\mathbf{X}\mapsto \mathbf{y} \right)$ between $\mathbf{X}$ and $\mathbf{y}$ based onthe training set D and to predict the most likely output value ${{f}_{*}}$ corresponding to the input point ${{x}_{*}}$ via the mapping relationship.

The above assume that the training set also holds true for the use of Gaussian processes(GPs) to solve regression problems. Generally, from the perspective of the function space, we define a GP to describe the distribution of function and make a Bayesian inference directly in the function space [29-30]. A GP is a set of limited random variables with joint Gaussian distribution and is determined by its mean function and covariance function completely. An $f\left( \mathbf{x} \right)$’s mean function and covariance function can be defined as

Therefore, a GP can be written as

Where $E \left( \bullet \right)$ denotes the operator of expectation.

In general, observation data ${{y}_{i}}$ can be described by the underlying function $f\left( {{\mathbf{x}}_{i}} \right)$, which is corrupted by noise $\delta $ for real-world regression problems. This can be defined as

Where $f\left( {{\mathbf{x}}_{i}} \right)$ and $\delta$ are two independent GPs, and $\delta$ is a GP with mean 0 and variance $\sigma _{n}^{2}.$

Therefore, the prior distribution of $\mathbf{y}$ can be expressed as

Meanwhile, the joint prior distribution of observations $\mathbf{y}$ and predictions ${{\mathbf{f}}_{\mathbf{*}}}$ is

Where ${{\mathbf{X}}_{\mathbf{*}}}$ denotes a ${{n}_{*}}\times d$ matrix of prediction set. For the input dataset $\mathbf{X}$ with $n$ samples and the prediction dataset ${{\mathbf{X}}_{\mathbf{*}}}$ with ${{n}_{*}}$ samples, ${{K}_{*}}=K\left( \mathbf{X,}{{\mathbf{X}}_{*}} \right)$ expresses a $n\times {{n}_{*}}$ matrix of covariance between the input and prediction samples. The same calculation rule is endowedwith $K=K\left( \mathbf{X,X} \right)$, ${{K}_{**}}=K\left( {{\mathbf{X}}_{\mathbf{*}}}\mathbf{,}{{\mathbf{X}}_{*}} \right)$,and $K_{*}^{T}=K{{\left( \mathbf{X,}{{\mathbf{X}}_{*}} \right)}^{T}}.$ ${{\mathbf{I}}_{n}}$ is the identity matrix. ${{\mathbf{f}}_{\mathbf{*}}}$ denotes the most likely output values corresponding the prediction dataset ${{\mathbf{X}}_{\mathbf{*}}}$.

Because the joint prior distribution of observations $\mathbf{y}$ and predictions ${{\mathbf{f}}_{\mathbf{*}}}$ is normal distribution, the posterior distribution of predictions ${{\mathbf{f}}_{\mathbf{*}}}$ can be directly calculated as

Where

${{\overline{\mathbf{f}}}_{\mathbf{*}}}$ and $\operatorname{cov}\left( {{\mathbf{f}}_{*}} \right)$ are respectively the mean and covariance of predictions ${{\mathbf{f}}_{\mathbf{*}}}$ corresponding to the prediction dataset ${{\mathbf{X}}_{\mathbf{*}}}$. Therefore, the posterior distribution of predictions ${{\mathbf{f}}_{\mathbf{*}}}$ will be obtained.

Various covariance functions can be selected when the GPR model is used to solve such problems. The most common covariance function is the squared exponential kernel:

Where $\theta =\left\{ l,\sigma _{_{f}}^{2},\sigma _{_{n}}^{2} \right\}$ is the parameter set of the GPR modeland is also called its hyper-parameters. $l$ is the length-scale, $\sigma _{_{f}}^{2}$ is the signal variance, and $\sigma _{n}^{2}$ is the noise variance. As we know, the simplest way to obtain the GPR model’s hyper-parameters $\theta $ is to maximize the log-likelihood function of the data set as shown in Equation (12):

The optimal hyper-parameters will be gotten through the above calculation. Then, Equations (9) and (10) will be used to calculate ${{\overline{\mathbf{f}}}_{\mathbf{*}}}$ and $\operatorname{cov}\left( {{\mathbf{f}}_{\mathbf{*}}} \right)$.

The following is the pseudo-code of the GPR algorithm according to the theory above:

Where a single prediction point ${{\mathbf{x}}_{\mathbf{*}}}$ is used in this pseudo-code, $\mathbf{\alpha }={{\left( K+\mathbf{\sigma }_{n}^{2}{{\mathbf{I}}_{n}} \right)}^{-1}}\mathbf{y}$, and $cholesky\left( K+\mathbf{\sigma }_{n}^{2}{{\mathbf{I}}_{n}} \right)$ denotes Cholesky decomposition.

To assess the prediction performance of the GPR model, we introduce a score function to quantitatively evaluate the performance of the GPR model that is relative to the training set, test set, and validation set. The score values indicate the robustness and accuracy of the GPR model. The closer the score value of the score function is to 1, the better the GPR model. The score values may be negative if the GPR model is bad enough. The specific expression of the score is

Where $u=\sum\limits_{l=1}^{N}{{{\left( {{T}_{\max ,\exp ,l}}-{{T}_{\max ,gpr,l}} \right)}^{2}}}$ and $v=\sum\limits_{l=1}^{N}{{{\left( {{T}_{\max ,\exp ,l}}-{{\overline{T}}_{\max ,\exp }} \right)}^{2}}}.$ ${{T}_{\max ,\exp ,l}}$ and ${{T}_{\max ,gpr}}_{,l}$ denote the l^th experimental value and the predictive value of ${{T}_{\max }}$ respectively. ${{\overline{T}}_{\max ,\exp }}$ denotes the mean of all ${{T}_{\max ,\exp ,l}}.$ $N$ denotes the number of samples in the training set, the test set, and the validation set respectively.

Meanwhile, the variable $\sigma $, defined as the mean relative error between ${{T}_{\max ,\exp }}$ and ${{T}_{\max ,gpr}}$ and shown in Equation (14), is also used to measure the performance of the GPR model, and the GPR model will demonstrate better performance if $\sigma $ has a lower value.

Where M denotes the number of all samples included in the training set, the test set, and the validation set.

When the GPR model is used to predict ${{T}_{\max }}$, the five features in the training set are inputted to the GPR model, and ${{T}_{\max }}$ is the model out put. According to the deviation between the model output and the maximum temperature measured by experiments, the GPR model is trained. In order to evaluate the performance of the trained GPR models, the data in the test set is inputted to the trained GPR models, and the obtained values of score and $\sigma$ will determine the best-suited GPR model. Further, the 17 samples in set2, called the validation set, are similarly used to validate the best GPR model through the values of score and $\sigma$. All the above work, including model training, data preprocessing, and model evaluation, are based on the programming language of python and relative machine learning package such as scikit-learn, pandas, and Numpy [26, 31].

For the training set and the test set in set1, the experimental ${{T}_{\max }}$ versus the predicted ${{T}_{\max }}$ by the GPR model are shown in Figure 2(a) and Figure 2(b)respectively. The experimental ${{T}_{\max }}$ versus the predicted ${{T}_{\max }}$ by the GPR model for set 2 is presented in Figure 2(c). The evaluation parameter values of $score $, $R$, and $\sigma$ of the GPR model are shown in Table 1.R denotes the correlation coefficient.

From the results shown in Figure 2(a), Figure 2(b), and Table 1, we can see that the GPR model shows excellent prediction performance on the training set and test set, especially on the training set, and the values of the score and R values are all over 0.9. Although the relative error of the test set is larger than the training set due to some poor fitting points, the predictive performance is still good for set1. As for the poor fitting points shown in Figure 2(b), the reason may be due to experimental mistakes or singularities [10]. No one can conquer the curse of original data for any regression model unless we eliminate the anomalous data in order to train a remarkable model.

For the validation set, the performance of the GPR model is not as good as set1. Not only is the score value in set 2significantly less than that in set1, but there are also several points that deviate sharply from the fitting curve, as shown in Figure 2(c). This indicates that the performance of the GPR model declines for set 2. The reasons may be due to factors such as experimental conditions, raw materials, and operation rules inset2.

Given the water content and the mass of sulfurized rust, operating temperature in the oil tank, air flow rate, and oxygen concentration in the safety valve, the prediction performance of ${{T}_{\max }}$ by the GPR model decreases sequentially on the training set, test set, and validation set. The prediction performance especially decreases on the validation set. However, the GPR model has high correlative values over 0.97 between ${{T}_{\max ,\exp }}$ and ${{T}_{\max ,gpr}}$, and the score values are very close to 1 on each dataset. For comparison, we also obtained the related results using the SVM model for the same problem, as shown in Table 2.

Evidentially, from the model evaluation results in Table 1 and Table 2, we can gain that the ${{T}_{\max }}$ prediction performance of the GPR model is superior to that of the SVM model in the score of score function, $R$ and the mean relative error between ${{T}_{\max ,\exp }}$ and ${{T}_{\max ,gpr}}.$ This proves the advantages of the GPR model in dealing with the problems of nonlinear uncertainties. It also verifies that the GPR method is more suitable for modeling observations with noise. Further, due to the more accurate predictions of ${{T}_{\max }}$ by the GPR model, it will be more favorable for sounding an alarm or initiating an automatic protection system for a fire or an explosion in an oil tank.

In this paper, the Gaussian process regression (GPR) method was applied to model the prediction of the maximum temperature (${{T}_{\max }}$) during the oxidation self-heating process of sulfurized rust. Five factors that have main effect son ${{T}_{\max }}$, such as water content, mass of sulfurized rust, operating temperature in the oil tank, air flow rate, and oxygen concentration in the safety value, are used as inputs of the GPR model. The validity and rationality of the GPR model are not only verified by the model parameter evaluation value based on the training set, test set, and validation set, but also compared with the SVM modelling results to reflect the excellent performance of the GPR. The GPR model achieves more accurate prediction of ${{T}_{\max }}$ than the SVM model, and it is more favorable for reducing the risk of fire and explosion of crude oil in production and transportation and preventing accidents.

The authors of this paper would like to thank Peng Hou and others in the State Key Laboratory of Explosion Science and Technology for providing an experimental space and technical support.