Efficient Identification of Crude Oil via Combined Terahertz Time-Domain Spectroscopy and Machine Learning

Abstract

The quality of crude oil varies significantly according to its geographical origin. The efficient identification of the source region of crude oil is pivotal for petroleum trade and processing. However, current methods, such as mass spectrometry and fluorescence spectroscopy, suffer problems such as complex sample preparation and a long characterization time, which restrict their efficiency. In this work, by combining terahertz time-domain spectroscopy (THz-TDS) and a machine learning analysis of the spectra, an efficient workflow for the accurate and fast identification of crude oil was established. Based on THz-TDS of 83 crude oil samples obtained from six countries, a machine learning protocol involving the dimension reduction of spectra and classification was developed to identify the geological origins of crude oil, with an overall accuracy of 96.33%. This work demonstrates that THz spectra combined with a modern numerical scheme analysis can be readily employed to categorize crude oil products efficiently.

1. Introduction

Crude oil is a complex mixture consisting of a variety of compounds [1]. Alongside liquid hydrocarbons like alkanes, cycloalkanes, and aromatic hydrocarbons, it also contains compounds containing nitrogen, sulfur, oxygen, and other elements, including hydrogen sulfide, quinoline, and phenols. Additionally, crude oil contains trace metals such as nickel, iron, and vanadium [2,3,4]. The composition and properties of crude oil are greatly influenced by the formation environment, resulting in variations in the quality and commercial value of crude oil from across the globe [5,6]. Rapid and accurate source identification of crude oil holds significant importance for crude oil trading and subsequent processing.

The existing methods for crude oil identification mainly include mass spectrometry methods such as gas chromatography–mass spectrometry (GC-MS), inductively coupled plasma–mass spectrometry (ICP-MS), Fourier transform–ion cyclotron resonance–mass spectrometry (FT-ICR MS), atmospheric pressure photoionization Fourier transform–ion cyclotron resonance mass spectrometry (APPI FT-ICR MS) [7,8,9,10], fluorescence spectroscopy [11,12,13], Raman spectroscopy [14,15], X-ray analysis [16,17], and infrared spectroscopy [18,19,20], etc. However, these methodologies exhibit distinct limitations. Mass spectrometry requires complex sample pretreatment, has a long experimental process, and yields complex fingerprint spectra [21], making it difficult to achieve classification analysis of multitudinous samples. Fluorescence spectroscopy fluctuates due to interference from scattered light, mutual interference between elements, and overlapping peaks. Raman spectroscopy suffers from weak Raman signal intensity and a low signal-to-noise ratio. X-ray spectroscopy is hazardous to human health and not extensively applicable. Currently, the most widely used methods in the field of crude oil detection and identification are attenuated total reflection infrared spectroscopy (ATR-IR) and Fourier transform attenuated total reflection infrared spectroscopy (ATR-FTIR), which are affected by uncertainties such as interfering molecules like resins in the obtained spectral range [22].

As a subfield of crude oil identification, crude oil classification most widely adopts the American Petroleum Institute (API) gravity method. Based on API gravity, crude oil can be classified into three categories: “light”, “medium”, and “heavy”. However, the use of standardized experiments to obtain data for crude oil analysis cannot be applied to geographical origin recognition. In order to achieve the faster and more efficient classification of crude oil, the combination of spectroscopic techniques with chemometrics or neural networks has attracted more and more scientific interest. As is shown in

Table 1. Applications of spectroscopic techniques: spectroscopic techniques combined with chemometrics or neural network methods in crude oil classification.

Type of Classification	Instrumental Methods	Chemometric Methods	Neural Networks	Ref.
Characteristics of trace elements	ICP-MS, GS-MS	CA 1	-	[9]
API gravity	FT-ICR MS	-	-	[10]
API gravity	ATR-FTIR	PCA, PLS-DA 2	-	[19]
Fuel and crude oil types	GS-MS	COW-PCA-LDA 3	-	[21]
Origin of crude oils	ATR-FTIR	PCA, SIMCA 4	-	[23]
Geographical origin	GS-MS	-	Kohonen self-organizing maps	[24]
Well and geographical origin	APPI FT-ICR MS	PCA, HCA 5	-	[25]
Certified reference materials	GC-MS	MPCA 6, PARAFAC 7	-	[26]

¹ CA: cluster analysis. ² PLS-DA: partial least squares-discriminant analysis. ³ COW-PCA-LDA: correlation optimized warping plus principal component analysis plus linear discriminant analysis. ⁴ SIMCA: soft independent modeling of class analogies. ⁵ HCA: hierarchical clustering. ⁶ MPCA: multi-way principal components analysis. ⁷ PARAFAC: principal factors analysis.

, the most commonly used methods for crude oil classification are a combination of mass spectrometry, ATR-FTIR, and chemometrics. However, inherent defects are inevitable in these methods, such as complex sample pretreatment and low sensitivity. In order to overcome these shortages and expand the application of spectroscopic techniques in the field of crude oil classification, this study employed terahertz time-domain spectroscopy (THz-TDS) to measure crude oil samples.

THz-TDS is a spectral measurement technique based on ultrafast femtosecond laser pulses [27]. Due to the high transmittance of THz waves in crude oil, it can measure the amplitude attenuation and time delay of the incident terahertz pulse without sample pretreatment, making it operationally simple. Additionally, the THz-TDS can be applied to analyze the compound structures and intermolecular interaction among organic molecules (such as hydrogen bonding, Van der Waals forces, dipole rotation, and vibrational transitions) [28,29,30,31]. It has been demonstrated that THz-TDS is applicable in qualitative and quantitative analysis of the crude oil properties [32,33], but its application in crude oil classification has not yet been explored.

THz-TDS with a frequency range from 0.2 to 2.5 THz was employed in this study to measure the terahertz time-domain spectra of 83 types of crude oil from six countries: Angola, Brazil, Saudi Arabia, Russia, Congo, and Iran. The refractive index spectra and absorption coefficient spectra of these crude oils were extracted and subjected to principal component analysis (PCA). Subsequently, the principal component of the refractive index spectra and absorption coefficient spectra were input into a convolutional neural network (CNN) to classify crude oils of different origins. The method of THz-TDS combined with PCA and CNN shows the new application of THz-TDS combined with machine learning in crude oil origin classification, which provides a new solution for the on-site classification of crude oil origin.

2. Experiment

2.1. Spectral Acquisition

The TDS-1008 system (BATOP GmbH, Jena, Germany, Figure S1), driven by a Mai Tai femtosecond laser (Spectra-Physics, Milpitas, CA, USA, with a wavelength of 780–980 nm, a pulse width of 80 fs, and an output power of 1.5–2.5 W) was used to conduct the THz-TDS measurements. The time domain signals were Fourier transformed to the frequency range of 0.2–2.5 THz, with a signal-to-noise ratio higher than 65 dB. The poly-4-methyl-1-pentene (TPX) cell adopted in the experiment had dimensions of 46 mm × 46 mm × 3 mm (length, width, and thickness, respectively). Every sample in the cell was placed into the THz-TDS chamber for measurement, and then when switching to a new sample, the previous sample was taken out, and the cell was cleaned and dried, and then injected with a new sample. The number of crude oil samples from different countries measured in the experiment is shown in Figure 1.

In addition, the refractive index and absorption coefficient of the sample are given by [34,35]

$$\mathrm{n}\left(\omega \right)=\left[\frac{\phi \left(\omega \right)c}{\omega d}\right]+1$$

(1)

$$\mathsf{\alpha }\left(\omega \right)=\frac{2}{d}\ln \left[\frac{4n\left(\omega \right)}{\rho \left(\omega \right){[n\left(\omega \right)+1]}^2}\right]$$

(2)

where $\omega$ is the angular frequency; $c$ is the velocity of light in vacuum; $d$ is the sample thickness; and $\phi \left(\omega \right)$ and $\rho \left(\omega \right)$ are the phase difference and amplitude ratio between the reference signal and the sample signal, respectively.

2.2. Principal Component Analysis (PCA)

PCA is a multivariate statistical analysis method that has been widely applied to deal with multivariable data, as it is capable of reducing the dimensionality of the data [33]. PCA transforms the original data into a set of linearly independent representations of each dimension via linear transformation, while retaining as much information present in the original dataset as possible [21]. By applying linear transformation, the original data variables are converted into a new set of variables known as principal components (PCs). Each principal component (PC) contains two types of information: scores and loadings. Scores represent the projection values of the original data on each PC, while loadings represent the weight information between the original variables and each PC [36,37]. The variance contribution rate is used to measure the importance of PCs, indicating the proportion of variance accounted for by a particular PC. The calculation formula for the variance contribution rate $R_i$ is as follows:

$$R_i=\frac{{\lambda }_i}{{\sum }_{k=1}^p{\lambda }_k}, i=1, 2, \dots ,p$$

(3)

where $R_i$ represents the eigenvalues (i.e., variances) obtained through eigenvalue decomposition of the covariance matrix of the original data and $i$ indicates the $i$-th PC. PCs with higher variance contribution rates have a greater degree of explanation for the original data. The cumulative variance contribution rate ${CR}_i$ refers to the proportion of the total variance of the first $i$ PCs, given by

$${CR}_i=\frac{{\sum }_{k=1}^i{\lambda }_k}{{\sum }_{k=1}^p{\lambda }_k}, i=1, 2, 3, \dots ,p.$$

(4)

2.3. Convolutional Neural Networks (CNNs)

CNNs, as one of the representative algorithms in deep learning, have been involved in remarkable achievements in various fields such as image recognition [38] and spectral recognition [39,40]. CNNs typically consist of convolutional layers, pooling layers, and fully connected layers. Generally, non-linear activation functions are added between the convolutional and pooling layers in CNNs to process the linear output data from the convolutional layers non-linearly, thereby enhancing the network’s representational capacity.

The two-dimensional CNN structure used in this study is illustrated in Figure 2. It comprises an input layer, two convolutional layers, two pooling layers, two fully connected layers, a SoftMax layer, and an output layer. The activation function applied in each convolutional layer is the exponential linear unit (ELU). This definition is shown in Equation (5):

$$ELU\left(x\right)=\left\{\begin{array}{c}\alpha \left(e^x-1\right), x<0\\ x, x\ge 0\end{array}\right.$$

(5)

where $\alpha$ is an adjustable parameter, usually with a value of 1.

The first convolutional layer consists of 6 convolutional kernels (or filters), while the second convolutional layer consists of 4 convolutional kernels. Both layers have a kernel size of 2 × 1 and a stride of 1. The input layer has a matrix size of 11 × 1 × 1, and the pooling layer has a matrix size of 2 × 1 and a stride of 1. The cross-entropy loss function (also known as the SoftMax function) is introduced to convert the model output into a classification probability, which measures the discrepancy between the output $q$ and the true labels $p$ for each class, and is given by

$$Loss=-\sum _{j=1}^m{p}_j\mathrm{l}\mathrm{o}\mathrm{g}\left(q_j\right)$$

(6)

where $m$ is the number of categories.

The Adam optimizer with an initial learning rate of 0.001 and a learning rate decay factor of 0.1 was used to train the model, and the trainable parameters were updated by minimizing the cross-entropy loss function in the training phase. The maximum number of iterations was set to 450, and at iteration 400, the learning rate was reduced to 0.0001. Furthermore, an L2 regularization parameter with a value of 0.0001 was introduced.

3. Results and Discussion

3.1. Terahertz Spectra of Different Origin Crude Oils

In Table S1, the origins and names of all crude oil samples are presented, along with their assigned numbers. Figure 3a,b displays the spectra of a subset of samples (two samples per origin) in the time and frequency domain, respectively. The complete THz spectroscopy for the crude oil samples is exhibited in Figures S2 and S3. In the figures, the solid black line represents the reference signal, while each color corresponds to a specific origin. The reference signal is the terahertz time-domain waveform of a terahertz wave passing through an empty sample cell. Solid and dashed lines of the same color distinguish different samples from the same origin.

Figure 3a shows the time delays and amplitude attenuations of the signals from different crude oil samples. The time delay and amplitude magnitude of the sample signals are influenced by the asphaltene content in the crude oil [22]. A lower asphaltene content results in a reduced time delay and amplitude attenuation transmitted terahertz pulse. The water content in the crude oil is another factor influencing the time delay and amplitude magnitude of terahertz waves [32]. Increased water content in the crude oil leads to greater time delay and amplitude attenuation in transmitted terahertz waves. Furthermore, the time delay and amplitude attenuation of the transmitted terahertz wave also fluctuate due to the complexity of crude oil composition and the presence of trace metals, sulfur, and nitrogen [40]. As shown in Figure 3b–d, the spectral features of different crude oil samples differ greatly in the range of 0.2–2.5 THz. Therefore, spectral data in the range of 0.2–2.5 THz were subjected to principal component analysis.

3.2. PCA-CNN Classification Model

In this study, PCA was adopted to investigate the refractive index spectrum and absorption coefficient spectrum of crude oil. The matrix was of the scale 83 × 1150 for both spectra, in which each row represented a sample and each column represented data points at various frequencies. PCs with a variance contribution rate exceeding 0.03% and the cumulative variance contribution rate of the total PCs more than 99.92% were chosen as the variables for the analysis. The variance contribution rates and cumulative variance contribution rates of the principal components of the refractive index spectra and absorption coefficient spectra are shown in

Table 2. Variance contribution rate and cumulative variance contribution rate of the first three PCs of the refractive index spectra.

Principal Component	Variance Contribution Rate/%	Cumulative Variance Contribution Rate/%
PC 1	99.688	99.688
PC 2	0.258	99.946
PC 3	0.033	99.979

and

Table 3. Variance contribution rate and cumulative variance contribution rate of the first eight PCs of the absorption coefficient spectra.

Principal Component	Variance Contribution Rate/%	Cumulative Variance Contribution Rate/%
PC 1	74.485	74.485
PC 2	22.301	96.786
PC 3	2.069	98.855
PC 4	0.647	99.503
PC 5	0.239	99.742
PC 6	0.099	99.841
PC 7	0.052	99.893
PC 8	0.035	99.928

, respectively.

By using PCA, the 2300 original spectral data including 1150 variables of the refractive index and absorption coefficient spectra separately were reduced to 11 variables including the first three PCs of the refractive index spectra and the first eight PCs of the absorption coefficient spectra. To eliminate the impact of differences in range of values between variables, the first three PCs of the extracted refractive index spectra and the first eight PCs of the absorption coefficient spectra were normalized. Subsequently, these normalized components were fed into a CNN for training and testing and the normalization formula is as follows:

$$X_{norm}=\frac{X-X_{min}}{X_{max}-X_{min}}$$

(7)

where X is the absolute value, $X_{min}$ is the minimum value, and $X_{max}$ is the maximum value of the dataset.

For countries for which the number of crude oil samples was less than or equal to 10, one sample was randomly chosen for testing, while the remaining samples were the training set. For countries with more than 10 crude oil samples, the training and test sets were divided at a ratio of 4:1. The k-fold cross-validation with k = 5 was implemented to validate the model. In this study, the training set data were divided into five mutually exclusive subsets utilizing fivefold cross-validation, among which one was selected for model verification and the remaining four were used for training in each iteration. This process was repeated five times, with a different subset acting as the validation set each time. The average classification accuracy of the five iterations was 97.14%. The confusion matrix of the test set, depicted in Figure 4, serves as a visual comparison between the actual category and the model prediction category. The figure reveals that only one sample was misclassified. The experiment was carried out 20 times and the average classification accuracy was 96.33%.

4. Conclusions

A methodology is proposed for categorizing the source of crude oil by employing THz spectroscopy in conjunction with PCA and CNN. To establish a THz spectroscopy dataset, THz-TDS technology was adopted to measure the spectra of 83 distinct crude oil samples from six countries. The PCA technique was then employed to reduce dimensionality and extract features from the refractive index spectra and absorption coefficient spectra of the crude oil samples, resulting in a combined PC dataset. Subsequently, a CNN was utilized for feature learning and classification recognition of the combined PC data. The experimental results demonstrate that the proposed approach, which combined THz-TDS, PCA, and a CNN, achieved a classification accuracy of 96.33% for distinguishing crude oil from different origins. This accomplishment confirms the feasibility and promising prospects of this method in the field of crude oil classification, extending the existing classification methods of crude oil.