An Unsupervised Clustering Method for Selection of the Fracturing Stage Design Based on the Gaussian Mixture Model

Abstract

In order to further improve the efficiency and economic benefits of multi-stage fracturing of unconventional oil and gas horizontal wells, it is urgently needed to conduct comprehensive reservoir quality evaluation research on the whole horizontal well section. Firstly, based on logging data, focusing on reservoir quality and completion quality, and comprehensively considering key factors such as reservoir physical property indexes and fracability indexes, a subjective and objective coupled evaluation model of the entropy weight method (EWM) and the analytic hierarchy process (AHP) without bias is established to obtain the composite reservoir quality index. Then, unsupervised gaussian mixture model (GMM) clustering algorithms are used to classify the reservoir comprehensive quality index and finally four grades of fracturing stages are established. Taking shale oil well A and B of the Permian Lucaogou Formation in Jimsar Sag, Junggar Basin, as examples, the comprehensive reservoir quality evaluation and clustering model training, testing, and prediction were carried out. By comparing the clustering results with the actual fracturing stages and oil production, it is found that the evaluation results obtained by the GMM clustering algorithms based on the coupled evaluation model of EWM and AHP can identify the good fracturing grades. The algorithm can also predict the fracturing grades of other wells in the same block. It proves the accuracy of the method proposed in this paper and provides a favorable technical basis for determining the placement of multi-cluster fracturing perforation.

1. Introduction

Currently, most unconventional reserves are characterized by a low grade, low and ultra-low permeability resources, and high commercial development costs [1]. Reservoir optimization evaluation and multi-stage fracturing are key technologies to reduce costs and increase efficiency in unconventional oil and gas development [2,3,4]. However, there are many indicators involved in reservoir optimization evaluation and inaccurate quantification [5,6], and when in multi-stage perforation fracturing, unreasonable selection of the fracturing stages leads to low production after fracturing [7,8]. Therefore, to define the basis of multi-stage and multi-cluster division of the horizontal wells, and to solve the problems of low accuracy and poor operability in the evaluation of reservoir quality and completion quality, we need to establish a set of comprehensive methods of evaluation based on well logging and mud logging data; this will realize being comprehensive and unify the evaluation parameters, as well as self-adaptive multi-stage optimization, to meet the demands of practical and efficient engineering applications.

Previously, domestic and foreign scholars have carried out some studies on reservoir evaluation. For example, Roberto et al. [9,10] considered the reservoir quality composed of hydrocarbon-filled porosity, permeability, pore pressure, organic carbon content, and organic maturity, and the completion quality composed of fracture height growth and fracability. Zou et al. [11] focused on the matching evaluation of the “six characteristics” of hydrocarbon source, lithology, physical property, brittleness, hydrocarbon potential, and stress anisotropy, and proposed eight evaluation criteria for unconventional oil and gas enrichment areas. Scholars at home and abroad have also carried out relevant studies on the optimization design of the fracturing stage. Lei et al. [12] proposed a geological engineering integrated reservoir reconstruction design platform, and the pre-fracture analysis module selected the fracturing well stage through comprehensive evaluation of reservoir quality and completion quality. Liu et al. [13] proposed stratigraphic modeling and evaluation technology for horizontal wells, integrating anisotropic reservoir quality evaluation and engineering quality evaluation methods such as in situ stress and brittleness to identify the sweet spot section of the horizontal wells and fracability, and finally optimized the horizontal well segmentation by integrating the stress and lithologic interval. Zhang et al. [14] used the unsupervised clustering algorithms based on Euclidean distance to cluster the reservoir according to reservoir flow and geomechanical parameters, and identified the fracability area of the fractured section of the reservoir so as to ensure the effectiveness of the perforation fracturing and improve the efficiency of the perforation fracturing. Slocombe et al. [15] compared fracture stages divided by engineering and geometric methods and found that effective stimulation was achieved by grouping formations with similar stress in engineering fracture stages. Cipolla et al. and Atanayev et al. [16,17,18] divided the well into several stages based on reservoir quality and completion quality by subjectively assigning quality evaluation and setting thresholds. Ejofodomi and Salah [19,20], based on Cipolla’s work, considered an anisotropic one-dimensional geomechanical model of a reservoir for completion quality, which improved the fracturing efficiency. Starting from the sweet spot evaluation of horizontal well engineering, Li et al. [21] comprehensively considered the factors of reservoir fracability indexes, physical property indexes, the fracture-induced stress field, and fracture pressure profile to divide the fracturing stages. Xia et al. [22] used the grey correlation analysis and multi-correlation coefficient method to evaluate the reservoir based on the characteristic quality of source rock, reservoir physical quality, and completion quality.

To sum up, for reservoir evaluation, foreign scholars usually use reservoir quality and completion quality to evaluate the reservoirs quantitatively through the subjective weighting method, with reservoir quality corresponding well with the reservoir physical properties such as high gas richness and high organic content, and completion quality corresponding to the reservoir’s good fracability; in turn, domestic scholars evaluates the six characteristics of source rock quality, completion quality, and reservoir physical quality to evaluate a reservoir by the objective weighting method. There are many evaluation indicators considered comprehensively by scholars at home and abroad, but it is difficult to evaluate reservoir quantitatively by unifying them, and the subjective and objective weighting methods are not combined to comprehensively evaluate a reservoir. As for the optimization design of fracturing stages, there are mainly geometric and engineering staging methods, where geometric staging means equal spacing of the perforating clusters. Because of the low production of the fracturing stages caused by geometric staging, this paper mainly discusses the engineering staging of the fracturing stages based on quality design. Most scholars divide the fracturing stages by subjectively setting thresholds based on the integration of geology and engineering, which has disadvantages such as strong subjectivity and poor operation. Some scholars considered the clustering algorithm based on Euclidean distance to improve the efficiency of the fracturing staging. This method has the problem of poor clustering, a result of non-spherical data sets. Therefore, this paper puts forward three quality evaluation methods—the reservoir quality, completion quality, and composite quality of the reservoir—and evaluates the three qualities of the reservoir based on the unbiased subjective and objective, coupled evaluation models of the entropy weight method and analytic hierarchy process. Then, unsupervised machine learning GMM clustering algorithms were used to classify the well segments and identify the fracturing, which provides a key technical support for multi-cluster fracturing and perforating.

2. Methodology

2.1. Calculation of Three Quality Indexes Based on Logging Data

Reservoir quality and completion quality are the keys to characterize the physical properties and fracability of the reservoir, and provide a basis for the later optimization of multi-stage and multi-cluster fracturing perforation; so, this article is based on well logging, such as the natural gamma ray, natural potential, density, and interval transit time, to extract the reservoir quality, completion quality, and composite quality index. The specific method is as follows.

2.1.1. Calculation of Reservoir Quality Index

The reservoir quality index is the main basis for logging the interpretation classification and fracturing fault selection, and the physical property of the reservoir is the main influencing factor, including the shale content, porosity, permeability, and oil saturation.

(1) For the calculation of effective porosity, the shale content is obtained according to the natural gamma ray and natural potential logging curves, and then the effective porosity is calculated by the Wylie formula.

$$\varphi =\frac{DT-\Delta t_{mp}}{T_f-\Delta t_{mp}}\cdot \frac{1}{C_p}-V_{sh}\left(\frac{T_{sh}-\Delta t_{mp}}{T_f-\Delta t_{mp}}\right)$$

(1)

where $\varphi$ is the effective porosity, %; $DT$ is the interval transit time logging value of the target layer; $\Delta t_{mp}$ is the interval transit time of rock skeleton, $\mathsf{\mu }\mathrm{s}/\mathrm{m}$; $T_f$ is the interval transit time of formation fluid, $\mathsf{\mu }\mathrm{s}/\mathrm{m}$; $V_{sh}$ is the shale content of formation, decimal; and $T_{sh}$ is the interval transit time of shale, $\mathsf{\mu }\mathrm{s}/\mathrm{m}$.

(2) The permeability calculation, based on the known porosity of Equation (1), Zeng [23] proposed a statistical empirical formula to obtain the permeability.

(3) The calculation of oil saturation is obtained by the Archie formula [24].

2.1.2. Calculation of the Completion Quality Index

Completion quality mainly considers the fracability of the reservoir and is used to characterize the difficulty of effective fracturing, including six evaluation indicators: the brittleness index, fracture toughness index, horizontal stress difference index, vertical stress difference index, natural fracture index, and bedding development index.

(1) Many scholars at home and abroad have done a lot of research on the calculation of the brittleness index. This paper mainly refers to the mineral content brittleness index method proposed by Jin et al. [25], the mechanical brittleness index method proposed by Rickman et al. [26], and the energy-based brittleness index method proposed by Ning Li et al. [27].

$$F_{I1}={\beta }_{1}B{I}_M+{\beta }_{2}B{I}_{YB}+{\beta }_{3}B{I}_E$$

(2)

where $B{I}_M$ is mineral brittleness index, dimensionless; $B{I}_{YB}$ is the mechanical brittleness index of rock, dimensionless; $B{I}_E$ is brittleness index of energy evolution, dimensionless; and ${\beta }_1, {\beta }_2, {\beta }_3$ are the weight values of each evaluation index to evaluate the brittleness index, dimensionless.

(2) The fracture toughness index is an important factor affecting the difficulty of the fracturing and reflects the ability of maintaining the fracture extension after fracture formation in the fracturing process. The calculation formula in this paper refers to the fracture toughness formula determined by Jin Yan et al. [28,29] using the multiple nonlinear regression method.

$$F_{I2}=\frac{KI{C}_{\max }-KIC}{KI{C}_{\max }-KI{C}_{\min }}$$

(3)

where $KIC$ is the fracture toughness, $\mathrm{Mpa}\cdot \sqrt{\mathrm{m}}$; and $KI{C}_{\max }, KI{C}_{\min }$ are the maximum and minimum fracture toughness, respectively, $\mathrm{Mpa}\cdot \sqrt{\mathrm{m}}$.

(3) The horizontal stress difference index is the difference between the maximum horizontal principal stress and the minimum horizontal principal stress, which is mainly affected by tectonic movement. Based on the transverse isotropic model [30], the calculation formula is as follows:

$$F_{\mathrm{I}3}=\{\begin{array}{l}0\Delta {\sigma }_h>{\sigma }_{h\max }\\ \frac{\Delta {\sigma }_{h\max }-\Delta {\sigma }_h}{\Delta {\sigma }_{h\max }-\Delta {\sigma }_{h\min }}\Delta {\sigma }_{h\max }\ge \Delta {\sigma }_h\ge \Delta {\sigma }_{h\min }\\ 1\Delta {\sigma }_{h\min }>\Delta {\sigma }_h\end{array}$$

(4)

where $\Delta {\sigma }_h$ is the horizontal stress difference value, $\Delta {\sigma }_h={\sigma }_H-{\sigma }_h$, $\mathrm{Mpa}$; $\Delta {\sigma }_{h\max }$ is the corresponding threshold of horizontal stress difference when forming hydraulic fractures; $\mathrm{Mpa}$; and $\Delta {\sigma }_{h\min }$ is the corresponding threshold of horizontal stress difference when the ideal hydraulic fractures is formed, $\mathrm{Mpa}$.

(4) The vertical stress difference index is the difference between vertical stress and the minimum horizontal principal stress. Its calculation formula is as follows:

$$F_{\mathrm{I}4}=\{\begin{array}{l}0\Delta {\sigma }_v\le \Delta {\sigma }_{v\max }or\Delta {\sigma }_v\le \Delta {\sigma }_{v\min }\\ 1-\frac{\Delta {\sigma }_v-\Delta {\sigma }_v^{\prime }}{\Delta {\sigma }_{v\max }-\Delta {\sigma }_v^{\prime }}\Delta {\sigma }_{v\max }>\Delta {\sigma }_v\ge \Delta {\sigma }_v^{\prime }\\ 1\Delta {\sigma }_v^{\prime }>\Delta {\sigma }_v>\Delta {\sigma }_{v\min }\end{array}$$

(5)

where $\Delta {\sigma }_v$ is the vertical stress difference, $\Delta {\sigma }_v={\sigma }_v-{\sigma }_h$, $\mathrm{Mpa}$; $\Delta {\sigma }_{v\max }$ is the corresponding vertical stress difference threshold when simple vertical fractures are formed by fracturing, $\mathrm{Mpa}$; $\Delta {\sigma }_{v\min }$ is the vertical stress difference threshold corresponding to the formation of simple bedding fractures during fracturing, $\mathrm{Mpa}$; and $\Delta {\sigma }_v^{\prime }$ is the optimal vertical stress difference value for forming the complex fractur network, $\mathrm{Mpa}$.

(5) For the natural fracture index calculation, the harmonic mean method is adopted to consider the influence of the natural fracture length, natural fracture density, and angle between the natural fracture orientation and horizontal maximum principal stress direction on fracability. The specific formula is as follows:

$$F_{I5}=\frac{3}{\left(\raisebox{1ex}{$(L_{\max }-L_{\min })$}\!\left/ \!\raisebox{-1ex}{$(L_i-L_{\min })$}\right.+\raisebox{1ex}{$({\rho }_{\max }-{\rho }_{\min })$}\!\left/ \!\raisebox{-1ex}{$({\rho }_i-{\rho }_{\min })$}\right.+\raisebox{1ex}{$45$}\!\left/ \!\raisebox{-1ex}{$\left(45-\left|45-{\theta }_i\right|\right)$}\right.\right)}$$

(6)

where $L$ is the normalized natural fracture length, $\mathrm{m}$; $L_{\max }$ is the maximum natural fracture length, $\mathrm{m}$; $L_{\min }$ is the minimum natural fracture length, $\mathrm{m}$; $\rho$ is the normalized natural fracture density line$/\mathrm{m}$; ${\rho }_{\max }$ is maximum fracture density, line$/\mathrm{m}$; ${\rho }_{\min }$ is minimum fracture density, line$/\mathrm{m}$; ${\theta }_i$ is the angle between the natural fracture and horizontal maximum principal stress direction, degree; and $\theta$ is the angle between normalized natural fracture orientation and horizontal maximum principal stress direction.

(6) Calculation of the bedding development index.

$$F_{\mathrm{I}6}=\{\begin{array}{ll}0& {\rho }_i\ge {\rho }_{\max }\hspace{0.17em}or{\rho }_i\le {\rho }_{\min }\\ 1-\frac{{\rho }_i-{\rho }^{\prime }}{{\rho }_{\max }-{\rho }^{\prime }}& {\rho }_{\max }>{\rho }_i\ge {\rho }^{\prime }\\ 1-\frac{{\rho }^{\prime }-{\rho }_i}{{\rho }^{\prime }-{\rho }_{\min }}& {\rho }^{\prime }>{\rho }_i>{\rho }_{\min }\end{array}$$

(7)

where ${\rho }_i$ is the bedding density; ${\rho }^{\prime }$ is the optimal bedding density; ${\rho }_{\max }$ is the bedding density threshold corresponding to the fracture morphology dominated by forming bedding; and ${\rho }_{\min }$ is the corresponding bedding density threshold when forming a fracture morphology dominated by hydraulic fractures.

2.1.3. Calculation of the Composite Quality Index

(1) Normalization of the original data

Due to the different physical meanings of various indicators for reservoir evaluation, normalization should be carried out before calculating the composite quality index, which can eliminate the influence of differences in the magnitude of absolute values of different parameters on the evaluation effect. The specific methods are as follows.

Suppose there is $n$ evaluation parameter indicators, and each evaluation parameter indicator has $m$ evaluation objects; then, the original data of the corresponding indicator of the evaluated object can be expressed as $R={\left(r_{ij}\right)}_{m\times n}$ by the following matrix, and the calculation is as follows:

Positive indictors:

$$s_{ij}=\frac{r_{ij}-\min \left\{r_{ij},\cdots ,r_{mj}\right\}}{\max \left\{r_{1j},\cdots ,r_{mj}\right\}-\min \left\{r_{1j},\cdots ,r_{mj}\right\}}$$

(8)

Negative indictors:

$$s_{ij}=\frac{\max \left\{r_{ij},\cdots ,r_{mj}\right\}-r_{ij}}{\max \left\{r_{1j},\cdots ,r_{mj}\right\}-\min \left\{r_{1j},\cdots ,r_{mj}\right\}}$$

(9)

where $s_{ij}$ is the parameters of the $i$th evaluation object of the $j$th evaluation indicator after normalization, dimensionless; and $r_{ij}$ is the parameter of the $i$th evaluation object of the $j$th evaluation indicator of original data ($i=1,2,\dots ,m$; $j=1,2,\dots ,n$).

(2) Calculation of the reservoir quality index:

$$Q_{RQ}={\displaystyle \sum _{i=1}^L{a}_i{G}_{Ii}}$$

(10)

where $Q_{RQ}$ is reservoir quality index, dimensionless; $G_{Ii}$ is the parameter value after normalization of the $i$th reservoir quality evaluation indicator, dimensionless; $a_i$ is the weight of the $i$th indicators, $L=4$.

(3) Calculation of the reservoir quality index:

$$Q_{CQ}={\displaystyle \sum _{i=1}^H{\alpha }_i{F}_{Ii}}$$

(11)

where $Q_{CQ}$ is completion quality index, dimensionless; and ${\alpha }_i$ is the weight of the $i$th indicator of completion quality, dimensionless, $H=6$.

(4) The composite quality index takes into account both reservoir quality and completion quality, and selects the best completion quality interval for fracturing. Then, the composite quality index of reservoir can be calculated as follows:

$$Q_{CQI}=W_1{Q}_{RQ}+W_2{Q}_{CQ}$$

(12)

where $Q_{CQI}$ is the composite quality index of reservoir, dimensionless; and $W_1$ and $W_2$ are the weight values of $Q_{RQ}$ and $Q_{CQ}$, dimensionless.

2.2. Coupled Evaluation Model of the Entropy Weight Method and Analytic Hierarchy Process

Based on the three quality indexes for the evaluation of the reservoir physical property and fracability, the core problem is the determination of the indicator weight. At present, the main weight calculation methods include subjective, objective, and combined weight methods [31]. The subjective weight method converts the subjective preferences of decision-makers into actual evaluation values, reflecting the leading role of the decision-makers, while the objective weight method completely relies on objective data information, completely abandoning the subjective initiative of the decision-makers and lacks flexibility in the decision-making process. Due to the complexity of unconventional oil and gas reservoirs, the subjective and objective weighting methods are often needed to comprehensively evaluate the reservoir, and finally optimize the subjective and objective weighting value. Therefore, this paper proposes a coupled evaluation model using the entropy weighting method and analytic hierarchy process.

2.2.1. The Entropy Weight Method (EWM)

EWM is based on the difference between the same indicators in the original data as the basis for weight determination [32]. The greater the difference between the evaluation indicators, the smaller the entropy value, and the more information this indicator contains and transmits, the greater the weight of the information obtained. Such methods are less affected by subjective factors, and different samples will have different weight values. For unconventional reservoirs with strong heterogeneity and large differences in physical properties, EWM can be combined with a variety of reservoir evaluation indicators to comprehensively weight the samples, so that the weights of different indicators are independent of each other, and only reflect the characteristics of strong or weak heterogeneity of the reservoir evaluation indicators, which can better realize the comparison of differences among the sample indicators. This method determines the weight value according to the relationship between the amount of information of the sample data and the probability of random events, which has a high calculation efficiency, and can be used to solve the weight value of the evaluation indicators more objectively.

Calculate the entropy weight of the $j$th evaluation indicator:

$$p_{ij}=\frac{s_{ij}}{{\displaystyle \sum _{i=1}^m{s}_{ij}}}$$

(13)

$$u_j=\frac{1+\frac{1}{\ln m}{\displaystyle \sum _{i=1}^m{p}_{ij}}\ln p_{ij}}{-{\displaystyle \sum _{j=1}^n\frac{1}{\ln m}{\displaystyle \sum _{i=1}^m{p}_{ij}}\ln p_{ij}}}$$

(14)

where $p_{ij}$ is the proportion of the $i$th evaluation object of the $j$th evaluation indicator, dimensionless; and $u_j$ is the entropy weight of the $j$th evaluation indicator, dimensionless, and the weight vector is $U={\left(u_1,u_2,\cdots ,u_j\right)}^T$.

2.2.2. The Analytic Hierarchy Process (AHP)

AHP is a subjective weight method to specify the weight value according to the judgment of the decision-makers. For the establishment of the layer structure of different evaluation indicators, it is decomposed into elements, and the elements are divided into several groups according to their properties, forming different levels. AHP avoids the deficiency of the objective weighting method without subjective preference. The weight setting of each layer will eventually directly or indirectly affect the result, and the influence degree of each factor in each layer on the result is quantified, very clear, and explicit.

(1) The establishment of the judgment matrix

The judgment matrix represents the importance of elements at a certain layer relative to elements at the next layer, and nine scales are introduced [33]. If there are $n$ evaluation parameter indicators, and each evaluation parameter indicator has $m$ evaluation objects, the judgment matrix is established $A={\left(a_{ij}\right)}_{n\times n}$.

(2) The weight calculation (summation method)

$$v_i=\frac{{\displaystyle \sum _{j=1}^n{v}_{ij}^{\prime }}}{{\displaystyle \sum _{j=1}^n{v}_j^{\prime }}}$$

(15)

where $v_{ij}^{\prime }$ is the normalization of the $j$th column vector of the judgment matrix, namely, $a_{ij}/{\displaystyle \sum _{i=1}^n{a}_{ij}}$; $v_i$ is the weight value of the $j$th evaluation parameter indicator, and the weight vector is $V={\left\{v_1,v_2,\cdots ,v_i\right\}}^T$.

(3) For the consistency test calculation [33], when $CR<0.1$, the judgment matrix is considered acceptable; otherwise, the judgment matrix should be modified appropriately.

2.2.3. Establishment of the Coupled Evaluation Model Entropy Weight Method and Analytic Hierarchy Process

In order to consider the shortcomings of the subjective and objective weights and minimize the total deviation between the evaluation results and the evaluation results under the subjective and objective weights, this paper adopts the quadratic programming method [34] to combine the subjective and objective weights without deviation. Establish the objective function:

$$\{\begin{array}{l}\min \hspace{0.17em}F(w)={\left[S(U-W)\right]}^T\left[S(U-W)\right]+{\left[S(V-W)\right]}^T\left[S(V-W)\right]\\ =4\left[\frac{1}{2}{W}^T(S^{T}S)W-{(\frac{1}{2}{S}^{T}SU+\frac{1}{2}{S}^{T}SV)}^{T}W\right]+b\\ s.t\hspace{0.17em}CW=1,W\ge 0\end{array}$$

(16)

where $W$ is the comprehensive weight value vector, $W={\left(w_1,w_2,\cdots ,w_n\right)}^T$; $S$ is the matrix after the normalization of the original data, $S={\left(s_{ij}\right)}_{m\times n}$; C is the $n\text{-}order$ unit row vector; and $b$ is an arbitrary constant.

Construct the Lagrange function for the objective function, getting

$$\left[\begin{array}{cc}{S}^{T}S& -C^T\\ -C& 0\end{array}\right]\left[\begin{array}{c}W\\ \lambda \end{array}\right]=\left[\begin{array}{c}\frac{1}{2}{S}^{T}SU+\frac{1}{2}{S}^{T}SV\\ -1\end{array}\right]$$

(17)

where $\lambda$ is the Lagrange multiplier. Let $H=S^{T}S$, $B=-S^{T}SU/2-S^{T}SV/2$, the comprehensive weight value vector of the solution $W=-QO+P^T$, where, $Q=H^{-1}-H^{-1}{C}^T{(C{H}^{-1}{C}^T)}^{-1}C{H}^{-1}$, $O={\left(C{H}^{-1}{C}^T\right)}^{-1}C{H}^{-1}$.

2.3. Gaussian Mixture Model (GMM) Clustering Algorithms

The classical k-means clustering algorithm often ends up with the local optimization of the SSE function of the clustering error, and is sensitive to noise and outliers [35,36]. It is not suitable for discovering non-spherical clusters or clusters with large size differences, and only applies to data sets with continuous attributes. However, the Gaussian mixture model (GMM) clustering algorithm is more flexible in processing data shapes than the k-means clustering algorithms [37,38], and the data set can be any non-spherical shape. GMM accurately quantifies a data set through a Gaussian probability density function by dividing the data set into several models based on the Gaussian probability density function. GMM is also the probability model with the fastest learning speed, which builds an appropriate linear GMM distribution combination model by fitting the input data set, so as to achieve the purpose of unsupervised clustering. Therefore, the characteristics of the spheroidal data set in this paper are suitable for clustering using the probability density model. GMM clustering algorithms are used to classify the reservoir composite quality, which makes the reservoir composite quality index of the same fracturing stage similar, and ensures that the reservoir quality and completion quality characteristics of the same clustering category are less different, while the reservoir quality and completion quality of different clustering categories are significantly different.

2.3.1. Standardization of Original Data

As the reservoir data changes with the depth of the formation, such original data sets fluctuate greatly, which will have a great impact on the clustering results. Therefore, in order to ensure that the data set is relatively stable and obeys a normal distribution, the Z-score method is adopted to standardize the composite quality index of the reservoir. Given $m$ samples and $n$ indicators, $X={\left(x_{ij}\right)}_{m\times n}$ is the original data matrix, and the standardized transformation is

$$x_{ij}^{\prime }=\frac{x_{ij}-{\overline{x}}_j}{S_j}$$

(18)

where ${\overline{x}}_j$ is the average value of the $j$th indicator; and $S_j$ is the variance of the $j$th indicator.

2.3.2. GMM Clustering Algorithms

The sample data set $X^{\prime }=\left\{x_{11}^{\prime },x_{11}^{\prime },\cdots ,x_{mn}^{\prime }\right\}$ of a known reservoir composite quality index conforms to a $K$ gaussian distribution; then, the GMM is as follows [39]:

$$p\left(x|\theta \right)={\displaystyle \sum _{k=1}^K{\alpha }_k{\mathcal{N}}_k}\left(x|{\theta }_k\right)$$

(19)

where ${\theta }_k=\left({\mu }_k,{\mathsf{\Sigma }}_k\right)$, and the unit gaussian distribution ${\mathcal{N}}_k\left(x|{\theta }_k\right)$, with mean ${\mu }_k$ and covariance matrix ${\mathsf{\Sigma }}_k$, is called a component of the GMM, where ${\alpha }_k$ is a mixed parameter, and is the weight of the $k$ gaussian distributions and represents prior probability:

$${\displaystyle \sum _{k=1}^K{\alpha }_k}=1,0\le {\alpha }_k\le 1$$

(20)

The probability density function of ${\mathcal{N}}_k\left(x|{\theta }_k\right)$ is

$${\mathcal{N}}_k\left(x|{\theta }_k\right)=\frac{1}{{\left(2\pi \right)}^{d/2}{\left|{\mathsf{\Sigma }}_k\right|}^{1/2}}{e}^{-1/2\left[{\left(x-{\mu }_k\right)}^T{\mathsf{\Sigma }}_k^{-1}(x-{\mu }_k)\right]}$$

(21)

2.4. Flowchart of the Identification of the Fracturing Grades

This paper is based on the logging or mud logging data to obtain the reservoir physical property parameters and fracability parameters, using the coupled evaluation model of EWM and AHP, for reservoir quality, completion quality, and composite quality evaluation; for the composite reservoir quality index for the division of the GMM clustering categories and identifying the grades of fracturing, the specific process is shown in Figure 1.

3. Case Study and Discussion

The horizontal section of well A selected in this paper is 4020.0~5725.0 m, and the reservoir porosity is 5.02~38.67%; the average porosity is 11.0%, the oil saturation is 45.0~100.0%, and the average oil saturation is 81.9%. The casing of the oil reservoir, with an outer diameter of 139.7 mm, is used for well completion. In the well section of 4020.0~5725.0 m, there are 15 stages, 145 clusters, a 981 m effective reconstruction section length, average cluster spacing of 6.8 m, and average segment spacing of 65.4 m. The interpreted horizontal section of well B is 4052.0~5718.0 m, the interpreted reservoir porosity is 0.12~31.3%, the average porosity is 7.7%, the oil saturation is 0.29~99.00%, and the average oil saturation is 63.9%. There are 103 clusters in 35 stages within the well section 4052.0~5718.0 m. The effective reconstructed section length is 915.5 m, the average cluster spacing is 15.8 m, and the average segment spacing is 48.3 m. After processing the cleaning data and outliers, this paper firstly selects the logging data of 4072.0–4472.0 m in well A to calculate the three quality indexes through the coupled evaluation model of EWM and AHP. Then the GMM clustering model was trained and classified by the composite quality index, and the clustering results were compared with the actual fracturing and oil production of the well section. Finally, the trained clustering model was used to predict and classify the reservoir quality of 4032.0–4506.0 m in well B and identify the grades of fracturing stages.

3.1. Analysis of the Coupled Evaluation Model of EMW and AHP

In this paper, the evaluation indicators of reservoir quality and completion quality were obtained based on logging data of well A, and the weight index of each indicator of reservoir quality and completion quality were obtained by EMW and AHP, respectively. Finally, the optimal comprehensive weight value was obtained by quadratic programming method.

As can be seen from

Table 1. Each indicator weight of reservoir quality.

Indicator	The Weight of Each Indicator
	EWM	AHP	Quadratic Programming
$V_{sh}$	0.39	0.055	0.2225
$\varphi$	0.209	0.1178	0.1634
$K$	0.213	0.2634	0.2382
$S_o$	0.187	0.5638	0.3754
$Q_{RQ}$	0.397	0.4500	0.4235

, when calculating the weight of each evaluation indicator of reservoir quality, the AHP is based on the knowledge of the decision-maker on reservoir oil content, and is the most important, so the weight of oil saturation is the largest. According to the heterogeneity of the reservoir, the entropy weight of the shale content is small, so the weight of information entropy weight is high. Finally, the quadratic programming method relies on the weight information and decision vector, limits the proportion balance of each index weight by minimizing the deviation, synthesizes the advantages of both subjective and objective weight assignment methods without preference, and finds the optimal solution of the comprehensive weight.

As can be seen from

Table 2. Each indicator weight of completion quality.

Indicator	The Weight of Each Indicator
	EWM	AHP	Quadratic Programming
$F_{I1}$	0.203	0.1818	0.1927
$F_{I2}$	0.122	0.4545	0.2885
$F_{I3}$	0.395	0.2727	0.333
$F_{I4}$	0.281	0.0909	0.1857
$Q_{CQ}$	0.603	0.550	0.5765

, the brittleness index, fracture toughness index, horizontal stress difference index, and vertical stress difference index were selected to evaluate the completion quality. When calculating the weight of each evaluation parameter of completion quality, the AHP firstly considers the fracture toughness index according to the fracability of the reservoir, so its weight is the largest. According to the entropy weight of the completion quality indicator, the weight value of the horizontal stress difference index is the largest, and the entropy weight of the fracture toughness index increases with the change of well depth, so its weight is the least.

Figure 2a and Figure 3a are the scatter plots of the subjective and objective weighting method and the quadratic combined weighting method for evaluating the three quality indexes, respectively. It can be seen that the quadratic programming method of optimization theory, as an optimal algorithm for solving nonlinear programming problems, can obtain the optimal solution of the combined weighting deviation function. For the weight evaluation results obtained based on the subjective and objective weighting method, the linear weighted grouping method (the objective coefficients were set as 0.6, 0.5, and 0.4 respectively) and the quadratic programming combined weighting method were used to combine the weights, and the results are shown in Figure 2b and Figure 3b. As can be seen from the figure, the robustness of the linear weighted grouping method is poor, and the weight changes of some evaluation indicator fluctuate greatly with different values of the objective coefficients. This is because it is easy to produce weight preference when a fixed value is used to roughly represent the importance of a certain weighting method, so that the importance of some evaluation indicators to the comprehensive weight value is artificially enhanced or reduced.

3.2. Determination of Clustering Number Value

The “Silhouette coefficient” is used to evaluate the clustering effect of GMM, so as to determine the clustering number $k$ value, which is selected according to the change in the silhouette coefficient with $k$ by combining the two factors of cohesion and separation degree. The silhouette coefficient is between −1 and 1, and the closer the distance between the samples in the cluster is, the farther the distance between the samples in the cluster is, the larger the overall silhouette coefficient is, and the better the clustering effect is. Figure 4 is a schematic diagram of $k$ determined by the silhouette coefficient method. As the value of $k$ increases, the silhouette coefficient increases first and then decreases. When the value is 5, the maximum silhouette coefficient is 0.993, which is the value in the red circle. Therefore, it can be seen that the GMM clustering algorithm is the best to divide the composite quality index into five clustering categories.

3.3. Cluster Simulation Experiment

After the standardized treatment of the composite quality index of well A, the training of the GMM clustering model and the classification of the cluster were carried out. Finally, the distribution boxplots of each evaluation indicator that constituted the three quality indexes under different clustering categories were drawn, respectively.

According to the distribution boxplot of the reservoir quality index in Figure 5, category 1 has a high shale content and oil saturation, and relatively low permeability and porosity; that is, the reservoir of Category 1 has the best physical properties. The comparison between Category 2 and Category 3 shows that the differences in shale content, permeability, porosity, and oil saturation are relatively small, indicating that the reservoir quality of Category 1 and Category 2 are similar. In Categories 4 and 5, the shale content is low, but the permeability, porosity, and oil saturation are high, indicating that the reservoir quality is better in Category 4.

Figure 6 shows the distribution boxplot of the completion quality indicators. It can be seen that Category 1 has a high brittleness index and a low fracture toughness index, horizontal stress difference index, and vertical stress difference index, and the reservoir has a good fracability; that is, the completion quality of Category 1 is the best. A comparison between Categories 2 and 3 shows that there are small differences in the brittleness index, fracture toughness index, horizontal stress difference index, and vertical stress difference index between Category 2 and Category 3, indicating that the completion quality of Category 2 and Category 3 is similar. The comparison between Categories 4 and 5 shows that Category 4 has the lowest brittleness index, but the fracture toughness index and horizontal stress difference index are lower than those of Category 5; so, the completion quality of the Category 4 reservoir is better than that of the Category 5 reservoir.

It can be seen from boxplot Figure 7 that the reservoir quality index and completion quality index of Category 1 are high, indicating that the composite quality index of Category 1 is the highest. The quality index differences between Category 2 and Category 3 were minimal by comparative analysis; so, Category 2 and Category 3 were classified into the same category. Compared with Category 4, Category 5 has a lower reservoir quality index, but it has a higher completion quality index and composite quality index. Therefore, the Category 5 reservoir is the least fracability and less prone to fracturing than Category 4. Figure 8a shows the gaussian density distribution of the clustering category of well A, indicating the distribution relationship among the five clustering categories and verifying the accuracy of the GMM clustering model.

Comprehensive analysis shows that Category 1 is the Grade 1 of fracturing grades; Categories 2 and 3 are Grade 2; Category 4 and Category 5 are Grade 3 and Grade 4, respectively. Typically, the last two grades are poor quality and can be selected according to the needs of the fracturing operation.

3.4. Cluster Test Experiment

Through the cluster simulation test, well A is clustered into four fracturing grades. Grade 1 and Grade 2 is the preferred fracturing, and Grade 3 and Grade 4 are selected according to the fracturing requirements. The fracturing grade of the cluster are compared with the fracturing stages and oil production of the well.

As shown in Figure 9, the section from 4024.0 to 4472.0 m in well A was divided into two grades, Grade 1 and Grade 4, based on reservoir quality index and completion quality. For reservoir quality (RQ), Grade 1 (blue) indicates the best reservoir properties, Grade 4 (red) indicates the worst reservoir properties, and for completion quality (CQ), Grade 1 (blue) indicates the best reservoir fracability and Grade 4 (red) indicates the worst reservoir fracability. RQ and CQ are weighted into composite reservoir quality (CQI), with Grade 1 (blue) representing the highest fracturing grade, Grade 2 (green) representing the second fracturing grade, Grade 3 (yellow) and Grade 4 representing the third and fourth fracturing grade, respectively, and the fracturing grade is low and not recommended. It can be seen that a fracturing grade with a high reservoir quality index does not necessarily have a high completion quality index. This is because although the evaluation of the reservoir quality is good, the corresponding fracturing effect is poor. Therefore, the comprehensive quality index was selected to classify the well segments. This method not only integrated the physical properties of reservoir quality and the fracability of the completion quality, but also avoided the subjective uncertainty caused by setting thresholds to classify the fracturing stages.

According to the composite quality index, as shown in

Table 3. Comparison of the test and actual fracturing stages in well B.

Fracturing Stages of GMM		Actual Fracturing Stages and Production
Fracturing stages/m	Grade	Fracturing stages/m	Oil production m3/d
4024.0~4089.9	2	4024.0~4088.0	2.536
4089.9~4109.0	4	Null	0
4109.0~4166.1	1	4109.0~4166.0	4.331
4166.1~4223.1	4	Null	0
4223.1~4307.3	2	4223.0~4307.0	2.499
4307.1~4351.1	3	Null	0
4351.2~4404.1	2	4351.0~4404.0	1.915
4404.1~4420.1	4	Null	0
4420.1~4472.0	2	4420.0~4472.0	1.896

, 4109.0–4166.1 m is the first-grade fracturing; 4024.0–4089.9 m, 4223.1–4307.3 m, 4351.2–4404.1 m, and 4420.1–4472.0 m are the second-grade fracturing. The 4307.1–4351.1 m was the third-grade fracturing, and the 4089.9–4109.0 m, 4166.1–4223.1 m and 4404.1–4420.1 m were the fourth-grade fracturing. By comparing the fracturing grades divided by GMM clustering model with the actual fracturing stages and oil production, it can be seen that Grade 1 and Grade 2 are exactly the fracturing grades required by the actual fracturing, while Grade 3 and Grade 4 have low evaluation grades and poor fracturing benefits, and these two grades are not used as the fracturing in the actual construction. The reservoir composite quality index of Grade 1 fracturing is the highest, and the oil production of the fracturing stage is also the highest at 4.331 m³/d after fracturing. Figure 10a shows the daily oil production of each fracturing stage under actual operation. Comprehensive analysis shows that the clustering result of well A is consistent with the actual fracturing stage and oil production, which confirms the rationality of the GMM clustering method selected in this paper.

3.5. Cluster Prediction Test

The trained GMM clustering model was used to predict the clustering category of fracturing grades and identify the fracturing grades of the 4032–4506 m reservoir composite quality index in well B, as shown in Figure 11. It can be seen from Figure 11 that the fracturing Grade 1 identified by clustering exactly corresponds to the fracturing stage with high production in actual construction, while Grade 3 and Grade 4 are the stages with no production or no construction in actual construction.

Table 4. Comparison of the prediction and actual fracturing stages in well B.

Fracturing Stages of GMM				Actual Fracturing Stages and Production
Fracturing stages/m	Grade	Fracturing stages/m	Grade	Fracturing stages/m	Oil production m3/d	Fracturing stages/m	Oil production m3/d
4032.0~4057.0	4	4281.8~4302.5	4	4032.0~4058.0	0	Null	0
4057.0~4075.4	3	4302.5~4327.3	2	Null	0	4303.0~4327.0	0.493
4075.4~4101.9	2	4327.3~4342.5	3	4077.0~4103.0	0.907	Null	0
4101.9~4120.8	3	4342.5~4370.0	1	Null	0	4345.0~4371.0	2.876
4120.8~4146.8	1	4370.0~4389.5	3	4122.0~4147.0	3.698	Null	0
4146.8~4165.9	4	4389.5~4410.9	3	Null	0	4391.0~4416.0	0
4165.9~4192.9	1	4410.9~4424.1	4	4166.0~4193.0	6.116	Null	0
4192.9~4211.1	3	4424.1~4462.5	4	Null	0	4434.0~4460.0	0
4211.1~4237.5	2	4462.5~4480.8	4	4211.0~4237.0	3.561	Null	0
4237.5~4255.9	3	4480.8~4506.0	2	Null	0	4483.0~4506.0	0.539
4255.9~4281.8	2			4256.0~4282.0	5.509

shows the comparison between the fracturing grades predicted by the GMM clustering model and the actual fracturing stages and oil production. It can be seen that the fracturing grades predicted by the clustering model are highly consistent with the actual fracturing stages. The analysis of the predicted results shows that the fracturing grades of 4032.0–4057.0 m, 4389.5–4410.9 m, and 4424.1–4462.5 m predicted by the GMM clustering model are 4, 3, and 4, respectively, and the fracture grades are low and difficult to fracture. However, in the actual fracturing, 4032.0–4058.0 m, 4391.0–4416.0 m, and 4434.0–4460.0 m were used as the fracturing stages. The oil production was 0 after fracturing (the gaussian density distribution of clustering categories and oil production of well B are shown in Figure 8b and Figure 10b, respectively), which indicates that the GMM clustering model can filter out the fracturing grades with low or no production and has good robustness.

4. Conclusions

In this paper, the logging data of the natural gamma ray, natural potential, density, and interval transit time are used to excavate various reservoir information, and the reservoir quality, completion quality, and composite quality are considered comprehensively, and the three quality indexes are accurately extracted. According to the logging data of shale oil wells A and B of the Permian Lucaogou Formation in Jimsar Sag, Junggar Basin, the reservoir was evaluated, and the evaluation results were used to classify the clustering categories and identify the fracturing grades. Then, the results of wells A and B were compared with the actual fracturing stage and oil production, and the following conclusions were obtained:

(1) A set of the EWM and AHP coupled evaluation model is established, which comprehensively considers the initiative of the subjective decision-maker and the characteristics of objective data when evaluating reservoir quality, completion quality, and composite quality, and integrates the subjective and objective weights impartially through a constrained nonlinear programming method of optimization theory. The evaluation results show that the coupling evaluation model can reflect the comprehensive quality of the reservoir more effectively and has a certain practical value.

(2) The GMM clustering algorithm is used to classify the reservoir composite quality, which can effectively identify the fracturing grade of the reservoir, and divide the fracturing grade into four grades, namely, I > II > III > IV. The convergence rate of the GMM clustering model is fast, and the efficient EM algorithm is adopted, which has a good clustering effect on the non-spherical data set with the change of well depth.

(3) According to the comparison analysis of the fracturing grades identified by the GMM cluster model with the actual fracturing stages and oil production of wells A and B, it is found that the first fracturing grade has the highest oil production, the second fracturing grade is lower than the first fracturing grade, while the third and fourth fracturing stages have not been constructed or have no oil production The accuracy of the GMM model proposed in this paper is verified, which provides a favorable basis for determining the perforation location of multi-cluster fracturing.