Probabilistic Fisher discriminant analysis based on Gaussian mixture model for estimating shale oil sweet spots

Abstract

The delineation of shale oil sweet spots is a crucial step in the exploration of shale oil reservoirs. A single attribute such as total organic carbon (TOC) is conventionally used to evaluate the sweet spots of shale oil. This study proposes a probabilistic Fisher discriminant approach for estimating shale oil sweet spots, in which the probabilistic method and Gaussian mixture model are incorporated. Statistical features of shale oil facies are obtained based on the well log interpretation of the samples. Several key parameters of shale oil are projected to data sets with low dimensions in each shale oil facies. Furthermore, the posterior distribution of different shale oil facies is built based on the classification of each shale oil facies. Various key physical parameters of shale oil facies are inversed by the Bayesian method, and important elastic properties are extracted from the elastic impedance inversion (EVA-DSVD method). The method proposed in this paper has been successfully used to delineate the sweet spots of shale oil reservoirs with multiple attributes from the real pre-stack seismic data sets and is validated by the well log data.

Graphical abstract

Keywords

probabilistic Fisher discriminant analysis / sweet spots / shale-oil facies / Bayesian inversion

1 Introduction

Shale oil is typically composed of both free oil and adsorbed oil, and is mainly retained or stored in the organic-rich source rocks (Raji et al., 2015; Wang et al., 2015). The probabilistic Fisher discriminant approach aims to delineate the shale oil sweet spots from the measured seismic data sets and well-log data. Total organic carbon (TOC), maturity, high Young’s modulus, and low Poisson’s ratio are the most important indicators of sweet spots in the shale reservoir (Zong et al., 2012b; Ouadfeul and Aliouane, 2016; Verma et al., 2016; Chen et al., 2018). The selection of proper reservoir inversion parameters is crucial in the accurate characterization of the shale reservoir. The four most important reservoir physical properties used to characterize unconventional reservoirs are the total organic carbon (TOC), porosity, pyrolysis S1, and clay content (Jarvie, 2012; Gorynski et al., 2019; Chen et al., 2020; Luo and Zong, 2020). This information together with fracturing-related features can be used to find and characterize the sweet spots of shale oil reservoirs.

The content of oil and gas is a key indicator of shale oil sweet spots. Particularly, the S1 value in pyrolysis of organic matter can be used to estimate the shale oil content in the source rock approximately and delineate the sweet spots of the shale oil reservoir (Song et al., 2013). In the pyrolysis of organic matter, gas is continuously released from the samples. Espitalie et al. (1977a) defined the first peak (S1) in the pyrolysis as the hydrocarbon content in the rock samples. Moreover, the S1 to TOC ratio, referred to as oil saturation index (OSI), can be used to determine the potential producibility of oil (Jarvie, 2012). The mobility of shale oil is affected by the properties of the retained oil and the light hydrocarbon component (Bao et al., 2016).

The pyrolysis S1 is related to the mobility of shale oil. Gorynski et al. (2019) used pyrolysis S1 data to assess the mobility of tight oil. The free or adsorbed hydrocarbon in the source rock can be characterized quantitatively by pyrolysis S1 from 3D pre-stack seismic data sets (Luo and Zong, 2019). Chen et al. (2020) used a large amount of Rock-Eval pyrolysis data in the Junggar Basin of China to evaluate tight oil resources. Pyrolysis S1 and fracture pressure gradient (FG) are key properties of any unconventional shale reservoir and the abundance of total organic carbon is the critical parameter to delineate geological sweet spots (Robert et al., 2015).

The pre-stack seismic inversion has been widely used to estimate TOC, clay content, and porosity. Ogiesoba and Hammes (2014) identified TOC-rich zones from 3D post-stack seismic attributes in the Eagle Ford Shale quantitatively. Yin et al. (2014, 2015) proposed a pragmatic approach, which combines the rock physics model and Bayesian inversion to directly inverse physical parameters of reservoirs from the elastic impedance of three different angles. Wang (2017) introduced a method for the pre-stack seismic inversion of shale gas sweet spots based on Yin’s works, and quantitatively inversed TOC, clay content and porosity. Ouadfeul and Aliouane (2016) inverted TOC using the multilayer perceptron neural network from three-dimensional seismic data in the Barnett Shale Formation. Yu et al. (2020) employed a rock-physics-based Bayesian method to inverse the total organic carbon in the Dongying Basin, China.

The Young’s modulus and Poisson’s ratio are directly related to the brittleness of rocks (Rickman et al., 2008). Yin et al. (2004) inversed Poisson’s ratio, P-wave impedance, and S-wave impedance from pre-stack seismic data sets based on elastic impedance. The shale oil sweet spots have relatively high Young’s modulus, but low Poisson’s ratio (Sena et al., 2011). Zong et al. (2012a) derived the YPD-Zoeppritz function, which uses the correlation among Young’s modulus, Poisson’s ratio, P-wave and S-wave velocity, and density to estimate Young’s modulus and Poisson’s ratio. The elastic impedance (EI) has been widely used in the inversion of reservoir attributes (Connolly, 1999). In the YPD-EI function (Zong et al., 2013), Young’s modulus, Poisson’s ratio, and density can be directly inverted from elastic impedance variation with incident angle. Huang (1984) proposed a model for the prediction of fracture pressure and considered several attributes such as overburden stress, rock strength, and Poisson’s ratio. Ma et al. (2020) utilized an AVAZ inversion method based on azimuth seismic data and derived the approximate correlation of the reflection coefficient, which is characterized by Young’s modulus, Poisson’s ratio, density, and fracture weaknesses parameters. Ma et al. (2020) calculated fracture pressure by an improved shale anisotropy fracture pressure model, which is based on Huang’s model (1984). Hu et al. (2015) assumed isotropic pressure coupling and estimated fracture pressure gradient (FG) quantitatively using the Poisson’s ratio, effective stress, and pore-pressure gradient.

Fisher discriminant analysis (Fisher, 1936), also known as Fisher’s LDA, was initially used to distinguish the two-species problem. It has been used in many fields such as text classification (AbuZeina and Al-Anzi, 2018), face recognition (Wang et al., 2016; Wu et al., 2017), and other data classification applications (Jiang and Chen, 2020; Guo et al., 2021). The FDA needs to find an optimal projection by maximizing the between-class distance and minimizing the within-class distance simultaneously. It is a supervised technique that requires class labels of the training data. Several effective Fisher discriminant functions have been proposed. However, the most commonly used discriminant functions are the ratio trace criterion (Rao, 1948), trace ratio criterion (Wang et al., 2007), and vectors ratio criterion (Foley and Sammon, 1975). Multiple inversed features of the shale oil reservoir can be used with Fisher’s discriminant analysis for the classifications and estimation of the shale oil sweet spots.

This study proposes a PFDA technique in which the probabilistic method is incorporated into the Fisher discriminant analysis to delineate the sweet spots of shale oil. This method could use multiple attributes from shale oil reservoirs to distinguish the sweet spots. The discriminate criterion was converted to a Gaussian mixture probability model, which was in turn expanded in the Fisher discriminant analysis method. Additionally, Young’s modulus, Poisson’s ratio, and density from real pre-stack seismic data sets were inverted based on the method proposed by Zong et al. (2012a, 2013). Inversion was performed for TOC, S1, clay content, and porosity using rock physics principles in the Bayesian method. With these inversion attributes of shale oil reservoirs, the sweet spots could be delineated by our PFDA method and the introduced method was applied to real seismic data and validated with the well logs.

2 Methods

2.1 Fisher discriminant analysis with probabilistic approach for sweet spots of shale oil

The main idea in the Fisher discrimination analysis is to project all types of data in a certain direction to obtain the highest separation distance among all projections. Fisher discrimination analysis attempts to find an optimal low-dimensional space to store the most discriminative information. FDA and its developments have been successfully applied in face recognition, speech recognition, and many other fields. However, Fisher discriminant analysis (FDA) focuses on a linear problem and due to the nonlinear separability of the shale oil sweet spots of different types, the FDA alone may not be able to achieve the delineation of the sweet spot of shale oil. Bayesian methods are often adopted to solve these linear restrictions. In this study, the conditional probability methods and the Fisher discriminant analysis (PFDA) are integrated and combined.

The information obtained from the lithofacies, petrophysical curves, and sedimentologic information was employed to categorize shale oil facies. Some reservoir attributes such as total content of organic matter, the thickness of source rock, maturity of kerogen, mobility of shale oil, and pyrolysis S1 may be used to distinguish the sweet spots in the delineation of shale oil. Given that there are m classes and n samples in different facies G_k, \mathrm{1}\le k\le m:

(1)\{\begin{array}{l}{G}_{\mathrm{1}}:{\mathit{x}}_{\mathrm{1}}^{\left(\mathrm{1}\right)},{\mathit{x}}_{\mathrm{2}}^{\left(\mathrm{1}\right)},\mathrm{...},{\mathit{x}}_n^{\left(\mathrm{1}\right)}\hfill \\ G_{\mathrm{2}}:{\mathit{x}}_{\mathrm{1}}^{\left(\mathrm{2}\right)},{\mathit{x}}_{\mathrm{2}}^{\left(\mathrm{2}\right)},\mathrm{...},{\mathit{x}}_n^{\left(\mathrm{2}\right)}\hfill \\ \hfill \mathrm{\cdots \cdots \cdots \cdots }\hfill \\ G_m:{\mathit{x}}_{\mathrm{1}}^{\left(m\right)},{\mathit{x}}_{\mathrm{2}}^{\left(m\right)},\mathrm{...},{\mathit{x}}_n^{\left(m\right)}\hfill \end{array},

where {\mathit{x}}_j{}^{\left(k\right)},\mathrm{1}\le j\le n, denote output features of a sample from the supervised facies. N=\mathrm{3}n, denotes the total samples of all classes.

Shale oil resources are organic-rich mudstones that have generated oil that is stored in the organic-rich mudstone intervals or continuous organic-lean intervals. The formation of shale oil resources can be categorized into three different types by their dominant organic matter and lithological composition. The three types are (G₁) organic-rich mudstones with high maturity, (G₂) organic-lean intervals of sand, and (G₃) organic-rich mudstones with low maturity. Several reservoir parameters influence the quality of the shale oil resources. Shale oil resource systems vary considerably from one system to another (Jarvie, 2012), but the best resources share some commonalities such as marine shales with excellent TOC values, high thermal maturity, high clay content, and relatively high mobility of shale oil. Generally, the high mobility of shale oil is associated with high free hydrocarbon content in the shale, which corresponds with a high pyrolysis S1 in the shale. In particularly, each type of shale oil facies in our study area has six main characteristics which enable the delineation and indication of the sweet spots of shale oil resources.

Based on the classification of G_{\mathrm{1}}, G_{\mathrm{2}}, G_{\mathrm{3}}, let {\mathit{x}}_j{}^{\left(k\right)}=\left[J_{\mathrm{1}},J_{\mathrm{2}},\mathrm{...},J_{\mathrm{6}}\right]denote the characteristic matrix of a sample. J_{\mathrm{1}},J_{\mathrm{2}},\mathrm{...},J_{\mathrm{6}} represent TOC, pyrolysis S1, clay content, Young’s modulus, porosity, and Poisson’s ratio, respectively.

The target of fisher discriminant analysis is to find the projection matrix \mathit{W}, by maximizing the between-class scatter distance, and meanwhile minimizing the within-class scatter distance (Wu et al., 2017), which can be formulated as

(2)\arg \underset{\mathit{W}}{\max }{\displaystyle \frac{{\mathit{W}}^{\mathrm{T}}{\mathbf{S}}_b\mathit{W}}{{\mathit{W}}^{\mathrm{T}}{\mathbf{S}}_{\omega }\mathit{W}}},

where {\mathit{S}}_w and {\mathit{S}}_b denote as the within scatter matrix and the between scatter matrix, respectively, which can be defined as

(3){\mathit{S}}_w={\displaystyle \sum _{k=\mathrm{1}}^m{\displaystyle \sum _{j=\mathrm{1}}^n({\mathit{x}}_j^{\left(k\right)}-{\mathit{\mu }}_k){({\mathit{x}}_j^{\left(k\right)}-{\mathit{\mu }}_k)}^{\mathrm{T}}}},

(4){\mathit{S}}_b={\displaystyle \sum _{k=\mathrm{1}}^m({\mathit{\mu }}_k-\mathit{\mu }){({\mathit{\mu }}_k-\mathit{\mu })}^{\mathrm{T}}},

where {\mathit{\mu }}_k is the mean vectors of atoms in the kth class G_k, and \mathit{\mu }={\displaystyle \frac{\mathrm{1}}{N}}{\displaystyle \sum _{k=\mathrm{1}}^m{\displaystyle \sum _{j=\mathrm{1}}^n{\mathit{x}}_j^{\left(k\right)}}} denotes the mean vector of all samples from all classes.

Specifically, the optimal solution (Wang et al., 2016) of Eq. (2) can be written as follows:

{\mathit{W}}^{*}=\arg \min -{\mathit{W}}^{\mathrm{T}}{\mathbf{S}}_b\mathit{W}

(5)s.t.{\mathit{W}}^{\mathrm{T}}{\mathbf{S}}_{\omega }\mathit{W}=\mathrm{1}.

We can use the Lagrange multiplier method to get Eq. (6):

(6){\mathit{S}}_b\mathit{W}=\lambda {\mathit{S}}_{\omega }\mathit{W}.

If the {\mathit{S}}_{\omega } is nonsingular, we can obtain the matrix form:

(7){\mathit{S}}_{\omega }^{-\mathrm{1}}{\mathit{S}}_b\mathit{W}=\mathit{W}\Lambda ,

where \mathit{W}=\left[w_{\mathrm{1}},\mathrm{...},w_q\right] the eigenvectors; \Lambda =\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[{\lambda }_{\mathrm{1}},\mathrm{...},{\lambda }_q\right],{\lambda }_{\mathrm{1}}\ge {\lambda }_{\mathrm{2}},\mathrm{...},\ge {\lambda }_q\ge {\lambda }_{q+\mathrm{1}}=\mathrm{0}, is a diagonal matrix of eigenvalues; q is the rank of the matrix of {\mathit{S}}_{\omega }^{-\mathrm{1}}{\mathit{S}}_b.

when we get the optimal projection direction, the raw data can then be projected onto a new hyperplane:

(8)y_j^{\left(k\right)}={\mathit{W}}^{\mathrm{T}}{x}_j^{\left(k\right)},

where y_j^{\left(k\right)} are the samples after the projection of the raw data; \mathit{W} is the optimal projection direction.

Assume that variables yin each type conform to the Gaussian mixture model with finite components N_d. Especially the variables y in each type G_jcan be generally sampled by Monte Carlo simulation and their distributions could be expressed as

(9)p\left(y|G_j\right)={\displaystyle \sum _{d=\mathrm{1}}^{N_d}{\pi }_d}N(y\mathrm{;}{\mathit{\mu }}_j^d,{\mathbf{\Sigma }}_j^d),

where N is the Gaussian distribution, {\mathit{\mu }}_j^d represents the vector of means in each type G_j; {\mathbf{\Sigma }}_j^d and {\pi }_d are the covariance matrix and the weights of each Gaussian distribution in the type G_j, respectively. Then the prior distribution (Grana and Rossa, 2010) can be written as

(10)p\left(y\right)={\displaystyle \sum _{n=\mathrm{1}}^{N_d}{\pi }_n}N(y\mathrm{;}{\mathit{\mu }}^n,{\mathbf{\Sigma }}^n),

where {\pi }_n represents the weights of each Gaussian distribution of all samples; {\mathbf{\Sigma }}^n and {\mathit{\mu }}^n are the covariance matrix and the vector of means of all samples. Then, the probability of each type G_jcan be obtained

(11)p\left(G_j\right)={\displaystyle \frac{p\left(y\right)}{p\left(y|G_j\right)}}={\displaystyle \frac{{\displaystyle \sum _{n=\mathrm{1}}^{N_d}{\pi }_n}N(y\mathrm{;}{\mathit{\mu }}^n,{\mathbf{\Sigma }}^n)}{{\displaystyle \sum _{d=\mathrm{1}}^{N_d}{\pi }_d}N(y\mathrm{;}{\mathit{\mu }}_j^d,{\mathbf{\Sigma }}_j^d)}}.

The type G_j of given sample y by comparing the probability of each category:

(12)G_j\to \arg \underset{j}{\max }P\left(G_j\right).

The algorithm used for solving PFDA is as follows:

2.2 Pyrolysis S1, TOC, clay content and porosity estimation from elastic impedance using Bayesian inversion

This study adopts Bayesian inversion to predict the pyrolysis S1, total organic carbon, clay content, and porosity from the elastic impedance data sets. In the Bayesian method, the posterior distribution of model attributes p\left(\mathit{m}|\mathit{d}\right)can be acquired through calculation of the prior distribution p\left(\mathit{m}\right)and the likelihood function p\left(\mathit{d}|\mathit{m}\right) (Grana et al., 2017):

(13)p\left(\mathit{m}|\mathit{d}\right)={\displaystyle \frac{p\left(\mathit{d}|\mathit{m}\right)p\left(\mathit{m}\right)}{p\left(\mathit{d}\right)}},

where \mathit{d}=\left[E{I}_{\mathrm{1}},E{I}_{\mathrm{2}},E{I}_{\mathrm{3}}\right] and \mathit{m}=\left[S\mathrm{1},\mathrm{T}\mathrm{O}\mathrm{C},\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{y},\mathrm{p}\mathrm{o}\mathrm{r}\right] are the elastic impedance and the model parameters, respectively. The statistical rock-physics model can be written as

(14)\left[E{I}_{\mathrm{1}},E{I}_{\mathrm{2}},E{I}_{\mathrm{3}}\right]={\mathbf{f}}_{\mathrm{R}\mathrm{P}\mathrm{M}}(S\mathrm{1},\mathrm{T}\mathrm{O}\mathrm{C},\mathrm{c}\mathrm{l}\mathrm{a}\mathrm{y},\mathrm{p}\mathrm{o}\mathrm{r})+\epsilon ,

where \epsilon is the error that represents the difference between model predictions and real data; and {\mathbf{f}}_{\mathrm{R}\mathrm{P}\mathrm{M}} is the rock-physics model.

The posterior distribution is affected by observed data \mathit{d} only through p\left(\mathit{d}|\mathit{m}\right), which is regarded as a function of \mathit{m}, referred to as likelihood function (Li et al., 2020):

(15)p\left(\mathit{d}|\mathit{m}\right)={\displaystyle \prod _{j=\mathrm{1}}^{j}p\left({\mathit{d}}_j|\mathit{m}\right)}.

We assume that the prior distribution is a Gaussian mixture model (GMM) with finite components N_c, which can be expressed as:

(16)p\left(\mathit{m}\right)={\displaystyle \sum _{n=\mathrm{1}}^{N_c}{\pi }_n}N(\mathit{m}\mathrm{;}{\mathit{\mu }}_{\mathit{m}}^n,{\mathbf{\Sigma }}_{\mathit{m}}^n),

where N is the Gaussian distribution; {\mathit{\mu }}_{\mathit{m}}^n, {\mathbf{\Sigma }}_{\mathit{m}}^n and {\pi }_n represent the vector of prior means, the covariance matrix and the weights of each Gaussian distribution, respectively.

Finally, the posterior distribution can be represented by a mixture Gaussian distribution:

(17)p\left(\mathit{m}|\mathit{d}\right)={\displaystyle \sum _{k=\mathrm{1}}^{N_c}{\lambda }_{\mathit{m}|\mathit{d}}^k}N(\mathit{m}\mathrm{;}{\mathit{\mu }}_{\mathit{m}|\mathit{d}}^k,{\mathbf{\Sigma }}_{\mathit{m}|\mathit{d}}^k),

where {\lambda }_{\mathit{m}|\mathit{d}}^kare mixture coefficients of the posterior, {\mathbf{\Sigma }}_{\mathit{m}|\mathit{d}}^k are the covariance matrices of each Gaussian component, and {\mathit{\mu }}_{\mathit{m}|\mathit{d}}^kare the conditional means. The petrophysical parameters can be attained by maximizing the posterior distribution:

(18)m^{*}\to \underset{m}{\arg \max }p\left(\mathit{m}|\mathit{d}\right).

2.3 Elastic impedance inversion for Young’s modulus and Poisson’s ratio

YPD function was used to estimate Young’s modulus and Poisson’s ratio and the properties were extracted from the elastic impedance inversion using the EVA-DSVD method proposed by Zong et al. (2013):

(19)EI\left(\theta \right)=A_{\mathrm{0}}{\left({\displaystyle \frac{E}{E_{\mathrm{0}}}}\right)}^{a\left(\theta \right)}{\left({\displaystyle \frac{\sigma }{{\sigma }_{\mathrm{0}}}}\right)}^{b\left(\theta \right)}{\left({\displaystyle \frac{\rho }{{\rho }_{\mathrm{0}}}}\right)}^{c\left(\theta \right)}.

If we take the logarithm on both sides of Eq. (19), then the Eq. (19) can be rewritten as

(20)\ln {\displaystyle \frac{EI\left(\theta \right)}{A_{\mathrm{0}}}}=a\left(\theta \right)\ln \left({\displaystyle \frac{E}{E_{\mathrm{0}}}}\right)+b\left(\theta \right)\ln \left({\displaystyle \frac{\sigma }{{\sigma }_{\mathrm{0}}}}\right)+c\left(\theta \right)\ln \left({\displaystyle \frac{\rho }{{\rho }_{\mathrm{0}}}}\right),

where A_{\mathrm{0}}, {\sigma }_{\mathrm{0}}, {\rho }_{\mathrm{0}}, and E_{\mathrm{0}} are the statistic average of elastic impedance, Poisson’s ratio, density and Young’s modulus, respectively; a\left(\theta \right), b\left(\theta \right), and c\left(\theta \right) are the coefficients associated with incident angles:

(21)a\left(\theta \right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\sec }^{\mathrm{2}}\theta -\mathrm{2}k{\sin }^{\mathrm{2}}\theta ,

(22)b\left(\theta \right)={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\sec }^{\mathrm{2}}\theta {\displaystyle \frac{(\mathrm{2}k-\mathrm{3}){(\mathrm{2}k-\mathrm{1})}^{\mathrm{2}}}{k(\mathrm{4}k-\mathrm{3})}}+\mathrm{2}k{\sin }^{\mathrm{2}}\theta {\displaystyle \frac{(\mathrm{2}k-\mathrm{1})}{(\mathrm{4}k-\mathrm{3})}},

(23)a\left(\theta \right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\sec }^{\mathrm{2}}\theta ,

Poisson’s ratio, density, and Young’s modulus can be determined from elastic impedances using EVA-DSVD method and Eq. (20).

Figure 1 illustrates the summarized workflow described thus far.

3 Results

The method proposed is applied to a shale oil area in southern China. The well logs in Fig. 2 show organic content, pyrolysis S1, clay content, porosity, Poisson’s ratio, and Young’s modulus. Based on the lithofacies, petrophysical curves and available sedimentologic information, the distribution of shale oil facies can be divided into three different shale oil facies (G_{\mathrm{1}},G_{\mathrm{2}},G_{\mathrm{3}}), as shown in Fig. 3. The classification is carried out by three steps: 1) Sand and clay facies were distinguished based on the lithofacies; 2) The distribution of G_{\mathrm{3}} was obtained by differentiating the facies containing a TOC smaller than 1.8 and the clay content below 0.36; 3) Then, the combination of porosity and pyrolysis was used as the means to classify the remaining part, where S1>7 mg/g and porosity>0.065 are indicators of the facies, and the rest is of G_{\mathrm{2}} type. The mean and standard deviation of various parameters are listed in Fig. 3, where variation in the statistical parameters is evident for different shale oil facies. For instance, G_{\mathrm{1}} has a relative high value of TOC, pyrolysis S1 and Poisson’s ratio compared with G_{\mathrm{2}} and G_{\mathrm{3}}. Besides, it has a lower value of Young’s modulus. The porosity of G_{\mathrm{1}} and G_{\mathrm{2}} is relatively lower than G_{\mathrm{3}}, but G_{\mathrm{1}} is slightly higher than G_{\mathrm{2}}.

Figure 4 displays the cross plot of the key evaluation parameters of the shale oil facies. Intersection analysis was performed on the pairwise plots of the features to primarily examine the influence of each parameter on the classification of the oil shale facies. The results show that each pair of the parameters can simply distinguish the facies of G_{\mathrm{1}}, G_{\mathrm{2}}, and G_{\mathrm{3}}. However, some overlap areas are detected in the cross-plots, which may cause errors in the characterization of shale oil sweet spots. The rock physics analysis results demonstrated in Figs. 4(b)–4(f), and Fig. 4(m) suggest that high-quality shale oil face G_{\mathrm{1}} has a relatively high value of TOC, S1, clay content, and Poisson’s ratio. Figure 4(a) and Fig. 4(i) indicate that high Poisson’s ratio and low Young’s modulus are the significant features of G_{\mathrm{1}}. In theory, high porosity is beneficial for a shale oil face, but there is no obvious difference in porosity in each facies according to Fig. 4(e), Fig. 4(j), Fig. 4(n), and Fig. 4(h).

After using Fisher dimension reduction, the original shale oil key parameters, which have six dimensions (TOC, Clay, S1, Por, E, σ), are turned to a new group of two-dimensional data sets as shown in Fig. 5, where G_{\mathrm{1}}, G_{\mathrm{2}}, and G_{\mathrm{3}} are depicted by circles of different colors. The coordinate “dimension 1” and “dimension 2” in Fig. 5 can be obtained by function 9. The optimal projection calculated by function 3 is

(24)\mathit{W}={\left[\begin{array}{l}\begin{array}{llll}\mathrm{0.1716}\hfill & \mathrm{0.0442}\hfill & -\mathrm{0.2896}\hfill & \begin{array}{ll}\mathrm{0.0216}\hfill & \begin{array}{ll}-\mathrm{0.00096}\hfill & \mathrm{0.9403}\hfill \end{array}\hfill \end{array}\hfill \end{array}\hfill \\ \begin{array}{llll}\mathrm{0.1073}\hfill & \mathrm{0.0668}\hfill & \mathrm{0.3668}\hfill & \begin{array}{ll}\mathrm{0.0246}\hfill & \begin{array}{ll}-\mathrm{0.0163}\hfill & \mathrm{0.9212}\hfill \end{array}\hfill \end{array}\hfill \end{array}\hfill \end{array}\right]}^{\mathrm{T}}.

Figure 5 demonstrates that the Fisher method can effectively realize the reduction of high-dimensional data, and furtherly separate different types of facies. Figure 5 shows a more integrated result, where all the previous factors are accounted for in generating the cross plot. After dimension reduction, the six factors of shale oil are converted to two-dimensional data sets in Fig. 5. There is some overlap between G_{\mathrm{1}} and G_{\mathrm{2}} due to the initial classification errors, but their core distributions are different. Gaussian mixture model in Eq. (9) can be used to rebuild the distribution of each type based on the new two-dimensional data.

The optimal projection \mathit{W} obtained in Eq. (24) was applied to the inversion results of shale oil key attributes depicted in Fig. 6 and two-dimensional data sets were generated. Figures 6(a)–6(f) are inversion results of clay fraction, total organic carbon content, Young’s modulus, pyrolysis S1, Poisson’s ratio, and porosity, respectively. Figure 7(a) is the delineation of sweet spots from TOC results shown in Fig. 6(b), where TOC>2.1 is the pre-determined criterion for the good shale oil sweet spots. Fig. 7(b) displays the sweet spots of shale oil reservoirs based on the classification of shale oil facies, determined by the PFDA method. The results depicted in Fig. 7(b) are validated by the well logs, where the final result obtained by the FPDA method is in good agreement with the well log. Although, G_{\mathrm{1}} exhibits a smaller distribution range in Fig. 7(b), the best sweet spots shown in Fig. 7(a) are consistent with the best shale oil face G_{\mathrm{1}} in Fig. 7(b).

4 Conclusions

A PFDA method was proposed to differentiate the shale oil sweet spots in this study, in which the probabilistic method and Gaussian mixture model are combined with Fisher’s method. This method was successfully used to delineate the sweet spots of shale oil reservoirs from the real pre-stack seismic data sets and is validated by the well log data. Additionally, a new workflow for the evaluation of shale oil sweet spots was introduced. Six key attributes of the shale oil sweet spots were selected based on the rock physics analysis, and these features were inversed by the Bayesian and EVA-DSVD methods. The PFDA was finally used to utilize more comprehensive and low-dimensional data to effectively delineate the sweet spots.

The accuracy of the PFDA method is relatively sensitive to the prior classifications. While significant errors in the inversion process can adversely affect the accuracy of the results, minimal error would not have a significant impact on the accuracy of the results generated by the PFDA method.

References

[1] AbuZeina D, Al-Anzi F S (2018). Employing fisher discriminant analysis for Arabic text classification. Comput Electr Eng, 66: 474–486

[2] Bao Y S, Zhang L Y, Zhang J G, Li J Y, Li Z (2016). Factors influencing mobility of Paleogene shale oil in Dongying Sag, Bohai Bay Basin. Oil & Gas Geo, 37(3): 408–414

[3] Chen S, Zhao W Z, Zeng Q C, Yang Q, He P, Gai S H, Deng Y (2018). Quantitative prediction of total organic carbon content in shale-gas reservoirs using seismic data: a case study from the Lower Silurian Longmaxi Formation in the Chang Ning gas field of the Sichuan Basin, China. Interpretation (Tulsa), 6(4): SN153–SN168

[4] Chen J Q, Pang X Q, Wang X L, Wang Y X (2020). A new method for assessing tight oil, with application to the Lucaogou Formation in the Jimusaer depression, Junggar Basin, China. AAPG Bull, 104(6): 1199–1229

[5] Connolly P (1999). Elastic impedance. Leading Edge (Tulsa Okla), 18(4): 438–452

[6] Espitalie J, Laporte J L, Madec M, Marquis F, Leplat P, Paulet J, Boutefeu A (1977a). Rapid method for source rock characterization and for determination of their petroleum potential and degree of evolution. Oil & Gas Sci Techn Revue, 32(1): 23–42

[7] Fisher R A (1936). The use of multiple measurements in taxonomic problems. Annals of Eugenics, 7: 179–188

[8] Foley D H, Sammon J W (1975). An optimal set of discriminant vectors. IEEE Trans Comput, C-24(3): 281–289

[9] Gorynski K E, Tobey M, Enriquez D, Smagala T, Dreger J, Newhart R (2019). Quantification and characterization of hydrocarbon-filled porosity in oil-rich shales using integrated thermal extraction, pyrolysis, and solvent extraction. AAPG Bull, 103(3): 723–744

[10] Grana D, Della Rossa E (2010). Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion. Geophys, 75(3): O21–O37

[11] Grana D, Fjeldstad T, Omre H (2017). Bayesian Gaussian mixture linear inversion for geophysical inverse problems. Math Geosci, 49(4): 493–515

[12] Guo R Y, Bai Y Q, Li C N, Shao Y H, Ye Y F, Jiang C Z (2021). Reverse nearest neighbors Bhattacharyya bound linear discriminant analysis for multimodal classification. Eng Appl Artif Intell, 97: 104033

[13] Hu R, Vernik L, Nayvelt L, Dicman A (2015). Seismic inversion for organic richness and fracture gradient in unconventional reservoirs: Eagle Ford Shale, Texas. Leading Edge (Tulsa Okla), 34(1): 80–84

[14] Huang R Z (1984). A model for predicting formation fracture pressure. J China U Petrol (Nat Sci), 8: 335–347

[15] Jiang C R, Chen L H (2020). Filtering-based approaches for functional data classification. Wiley Interdiscip Rev Comput Stat, 12(4): 1–15

[16] Jarvie D M (2012). Shale resource systems for oil and gas: part 2—shale-oil resource systems. AAPG Mem, 97: 89–119

[17] Luo K, Zong Z Y (2020). Pyrolysis S1 discrimination of shale gas with Pre-stack seismic inversion. In: Annual Meeting of Chinese geoscience union 2020, Beijing

[18] Li K, Yin X Y, Zong Z Y (2020). Facies-constrained pre-stack seismic probabilistic inversion driven by rock physics. Scientia Sinica Terrae, 50(6): 832–850

[19] Luo K, Zong Z Y (2019). Sweet spots discrimination of shale gas with Pre-stack seismic inversion. In: Annual Meeting of Chinese Geoscience Union 2019, Beijing, 53: 29–31

[20] Ma L, Lin Z L, Hu H F, Zhou D (2020). Seismic prediction method of fracture pressure in a shale formation. In: SEG International Exposition and 90th Annual Meeting, 1068–1072

[21] Ouadfeul S, Aliouane L (2016). Total organic carbon estimation in shale-gas reservoirs using seismic genetic inversion with an example from the Barnett Shale. Leading Edge (Tulsa Okla), 35(9): 790–794

[22] Ogiesoba O, Hammes U (2014). Seismic-attribute identification of brittle and TOC-rich zones within the Eagle Ford Shale, Dimmit County, South Texas. J Pet Explor Prod Technol, 4(2): 133–151

[23] Raji M, Gröcke D R, Greenwell C, Cornford C (2015). Pyrolysis, porosity and productivity in unconventional mudstone reservoirs: free and adsorbed oil. In: Unconventional Resources Technology Conference, 11: 270–279

[24] Rickman R, Mullen M J, Petre J E, Grieser W V, Kundert D (2008). A practical use of shale petrophysics for stimulation design optimization: all shale plays are not clones of the Barnett Shale. In: SPE Technical Conference and Exhibition

[25] Rao C R (1948). The utilization of multiple measurements in problems of biological classification. J R Stat Soc B, 10(2): 159–193

[26] Sena A, Castillo G, Chesser K, Voisey S, Estrada J, Carcuz J, Carmona E, Hodgkins P (2011). Seismic reservoir characterization in resource shale plays: “sweet spot” discrimination and optimization of horizontal well placement. In: 81st Annual International Meeting, SEG, Expanded Abstracts, 1744–1748.

[27] Song G Z, Lin Y, Lu S F (2013). Resource evaluation method for shale oil and its application. Earth Sci Front, 20(4): 221–228

[28] Verma S, Zhao T, Marfurt K J, Devegowda D (2016). Estimation of total organic carbon and brittleness volume. Interpretation (Tulsa), 4(3): T373–T385

[29] Wang H J (2017). Prediction of geological dessert in shale gas with pre-stack seismic inversion. Dissertation for the Master’s Degree. Qingdao: China University of Petroleum (East China)

[30] Wu L, Shen C H, Hengel A (2017). Deep linear discriminant analysis on fisher networks: a hybrid architecture for person re-identification. Pattern Recognit, 65: 238–250

[31] Wang S, Feng Q, Javadpour F, Xia T, Li Z (2015). Oil adsorption in shale nanopores and its effect on recoverable oil-in-place. Int J Coal Geol, 147–148: 9–24

[32] Wang H, Yan S C, Xu D, Tang X O, Huang T (2007). Trace ratio vs. ratio trace for dimensionality reduction. In: Proceedings of the Conference on Computer Vision and Pattern Recognition. Minneapolis, USA: IEEE, 1–8

[33] Wang Z, Ruan Q, An G (2016). Facial expression recognition using sparse local Fisher discriminant analysis. Neurocomputing, 174: 756–766

[34] Yin X Y, Cui W, Zong Z Y, Liu X J(2014). Petrophysical property inversion of reservoirs based on elastic impedance. Chinese J Geophys, 57(12): 4132–4140

[35] Yu J Q, Yu Z Q, Mao Z Q,Gao G, Luo K, Lei T, Zong Z Y(2020). Prediction of total organic carbon content in source rock of continental shale oil using pre-stack inversion. Geophys Prospect Petrol, 59(5): 823–830

[36] Yin X Y, Zong Z Y, Wu G C (2015). Research on seismic fluid identification driven by rock physics. Sci China Earth Sci, 58(2): 159–171

[37] Yin X Y, Yuan S H, Zhang F C (2004). Rock elastic parameters calculated from elastic impedance. CPS/SEG Technical Program Expanded Abstracts

[38] Zong Z Y, Yin X Y, Wu G C (2013). Elastic impedance parameterization and inversion with Young’s modulus and Poisson’s ratio. Geophysics, 78(6): N35–N42

[39] Zong Z Y, Yin X Y, Zhang F,Wu G C (2012a). Reflection coefficient and pre-stack seismic inversion with Young’s modulus and Poisson’s ratio. Chinese J Geophys, 55(11): 3786–3794

[40] Zong Z Y, Yin X Y, Wu G C (2012b). Pre-stack inversion for rock brittleness indicator in gas shale. In: Annual meeting of Chinese Geoscience Union 2012. Beijing, 17: 495

RIGHTS & PERMISSIONS

Higher Education Press