Pattern Recognition of the Vertical Hydraulic Fracture Shapes in Coalbed Methane Reservoirs Based on Hierarchical Bi-LSTM Network

Pattern recognition of the hydraulic fracture shapes is very important and complex for the refracturing design of coalbed methane (CBM) wells. In this paper, we explore a new idea by regarding the pattern recognition process as understanding what a CBM reservoir “says” during hydraulic fracturing. *en we present a hierarchical Bidirectional LSTM (Bi-LSTM) network to recognize the pattern of hydraulic fracture geometry in CBM reservoirs. Inputting the wavelet denoised sequences of data to the presented network, we can extract the implicit features of the hydraulic sand fracturing operation curves and automatically combine them to make the classification of the fracture shapes. With this method, we can cope with the problems happened in early stage of the CBM field development such as the lack of monitoring wells and the information of rock mechanics. Moreover, the experiences of the engineers and the measured data are combinationally used, which can efficiently reduce the subjectivity and assist the engineers to make the refracturing design. *e validity of this method is verified by the testing data and comparing with the simulated results of Fracpro PT software.


Introduction
China has abundant coalbed methane (CBM) resources. Until 2018, there have been 26 CBM fields found in Qinshui, Hancheng, etc. e total proved reserve of CBM in China is 6345 × 10 8 m 3 , among which the economical recoverable reserve is 2537 × 10 8 m 3 [1]. e exploitation of CBM can not only help reduce the security risk during coal mining, but also can regulate the structure of the energy supply of China. Different from the geological condition of the coalbeds in US or Canada, the CBM reservoirs in China are of low porosity and low permeability, which means that the industrial exploitation of CBM has to be realized through reservoir stimulation. Hydraulic fracturing (HF) is a key technology to promote the development of tight reservoirs (e.g., CBM ) by enhancing production and recovery rates and improving their exploitation economics. It is implemented to increase reservoir permeability and productivity by creating and extending the fracture network with high-pressure fluid injection [2,3]. However, permeability and conductivity of the hydraulic fractures in coalbeds are easily reduced because of stress sensitivity or coal powder migration, which will lead to negative effect on the exploitation of CBM; see Figure 1. us it is necessary to implement the refracturing operation to create new fractures or extend the existed ones, so that the production of CBM can be improved.
Before refracturing, it is a key point to analyze the possible shape of the existed fractures, which has direct effect on the selection of fracturing techniques, scales, and proppant schedules. Now there are a lot of research studies focusing on the pattern recognition of fracture shapes. For example, the direction, length, height, and growth history of hydrofracture propagation can be obtained by microseismic monitoring [4]. However, the technique of microseismic has high requirement on the monitoring well and well site, and the measurement error of the fracture position sometimes is really high. e fracture azimuth, dip, and volume can be measured by surface tiltmeter [5], but the result is easily influenced by the position of the measuring instruments [6]. Moreover, the size of individual fracture cannot be obtained with this method in complex fracture network. e other evaluation techniques [7], such as the memory well temperature logging and the radioactive tracing, can only determine the vertical height of the fractures in near borehole zones, while the information of the fractures 1 m beyond the borehole in the horizontal direction is not able to be obtained. e Potentiometric method can measure the fracture azimuth and length exactly, but the result depends on the difference of the salinity between the fracturing fluid and the reservoir [8]. ese methods aforementioned can help monitor and evaluate the hydraulic fractures "directly", the application of which is limited by the expense of instruments and the disturbance from the reservoir properties.
"Indirect" method, the technique of pressure analysis, is of low cost and able to obtain more reservoir information. It is considered to be the most powerful and economical method to monitor and evaluate the hydraulic fractures. e traditional pressure analysis method is originated from the work of Nolte et al. [9][10][11] etc. is method can help determine the propagation of the fractures, based on the analysis of the slope information of the double-logarithmic curves of net pressure and time. Until now, there has been a lot of investigation on this method. Martins [12] developed and improved the fracturing pressure analysis technology, making it become a classic analytical in the fracturing pumping process. Crokett et al. [13] studied fracture pressure analysis method based on three-dimensional fracture extension mode. Benelkadi et al. [14] analyzed the fracturing pressure after stopping the pump. However, in practical situation, it is not effective to judge the shapes of the hydraulic fractures only by the analysis of the double-logarithmic curves. e essence of the above indirect analysis method and the related modification is to extract effective numerical features (slope, etc.) through the observation of the hydraulic fracture operating curves. With the extracted "meaningful" features, the operation supervisors can recognize the pattern of the fracture shapes empirically. e parameters during this process are often tuned by the supervisor, and thus the accuracy of the prediction depends on his own experience. For feature extraction, there are a lot of useful techniques in machine learning [15,16]. Along with the development of deep learning, there have been many powerful tools to extract the features (meaningful/meaningless) of different data types and apply them for pattern recognition [17]. For example, the convolutional networks are applied to extract the features of digital figures for image classification [18][19][20][21][22][23]. However, for the tasks involving sequential inputs, such as the speech recognition, it is often better use recurrent neural networks (RNNs), which can maintain the implicit historic information by the hidden units [24][25][26]. e RNNs have a defect that it is difficult to learn to store information very long because of the gradient exploding or vanishing [27]. To correct this, the long shortterm memory (LSTM) networks are presented to remember the inputs for a long time by using special hidden units [28]. LSTM and its modification (bidirectional LSTM, Bi-LSTM, etc.) have been applied effectively in the processing of speech or language [29][30][31][32].
A set of hydraulic sand fracturing operation curves for CBM reservoir consists of 3 time series: the tubing pressure, the proppant ratio, and the fluid rate; see Figure 2.
Regarding the above curves as the "Speech" from the CBM reservoir during fracturing process, the analysis of the hydraulic fracture operating curves can be taken as that we are trying to understand what the CBM reservoir "say". Inspired by this observation, we present a hierarchical Bi-LSTM network to recognize the pattern of hydraulic fracture geometry in CBM reservoirs; see Figure 3.

Complexity
In this framework, the input data is denoised by wavelet transform firstly, as we care more about the trends and the smoothed shapes of the operation curves. e hidden hierarchical layer consists of five Bi-LSTM layers [33][34][35], the output of which is the input of the classification layer. Here, the hierarchical Bi-LSTMs are regarded as encoders to extract the features. e advantage of doing this is that we can unfold the framework to analyze the hidden feature information explicitly: the type of the operation curves, sensitivity to proppant ratio, the level of pressure decline rate after stopping pump, and the manner of fracture propagation. With the presented method, the pattern of the fracture shapes in CBM reservoir can be recognized more rapidly and efficiently. is method can cope with the problems happened in early stage of the coalbed gas field such as the lack of monitoring wells and the information of rock mechanics. Moreover, the presented method does not require to fit the formation stress profile, so it also helps to improve the conventional analyzing methods for double-logarithmic curves.
e main body of the paper is organized as follows: in Section 2, the patterns of the fracture shapes in CBM reservoirs and the related influence factors (the features) are analyzed. In Section 3, the details of our method are separated into 2 parts, the data preparation and the hierarchical Bi-LSTM encoder for classification. In Section 4, the whole network is trained by the data of over 1000 wells of one coalbed gas field in China. e validation of the presented method is verified by comparing the predicted results with the simulated results of Fracpro PT.

Patterns of the Fracture Shapes in Coalbeds and the Related Influencing Factors
e orientation of the induced fractures in coal seam is dominated by the vertical principal stress and the horizontal principal stress of the reservoir mutually, which makes the main fracture plane vertical or horizontal. Moreover, disturbed by many subsurface uncertainties, it is almost impossible to predict the exact geometry of the fractures. Luckily, in reservoir stimulation engineering, the engineers care more about the pattern of the fracture shapes. e data of this work mainly comes from a coalbed gas field in Shanxi of China, where the hydraulic fractures are mainly vertical. us, in this paper, we mainly focus on the pattern recognition of vertical hydraulic fractures.

Patterns of the Fracture Shapes in Coalbeds.
Recognition of the classic fracture shape patterns is very important to fracture monitoring and evaluation. Before the refracturing operation for a CBM well, the fracturing engineers need to make the diagnosis empirically on the subsurface working condition according to the former fracture shape patterns; see Figure 4. Combining the information of drainage characteristics and the geological potential information of the CBM well, they can offer suitable suggestions for the fracturing design.
For example, there was a well, TS07-3D, in a coalbed field in China. Before the refracturing operation, we found that the production of this well is mainly controlled by reservoir permeability. Because of the jam of coalbed ash and that the pattern of the fracture shapes is High-Short type (the hydraulic fractures were not fully extended), the production rate of this well was very low after the first fracturing. With the expert system MDSCBM we complied, the suggestion for the design of refracturing operation for this well is obtained: (1) Implement the perforating and fracturing at the horizon near the top of the coal seam where the principal stress is minimal (2) e proppant schedule should be impulse sand adding, and the propping agent should be the sand of low density According to the field experience and simulation results, the shapes of the vertical hydraulic fractures in CBM reservoirs are mainly divided into 4 patterns in this paper. In Here, we should notice that the names of the four patterns in Figure 4 are just notations. e more important thing behind is the phenomenon that each name of the fracture patterns stands for.
(a) Short-Wide Fracture. e fluid loss is too much; thus, the fracture length cannot be extended normally, and the fluid is only used to expand the fracture width.
(b) Long-Narrow Fracture. e penetration of the fracture is good and the height of the fracture proper enough. is is the ideal fracture-making. (c) Multiple-Narrow Fracture. e natural fractures in this horizon are developed and the artificial fractures are multiple ones of small width. (d) High-Short Fracture. In this situation, the height of the fracture is seriously out of control and penetration of the fracture is relatively bad.

e Influencing Factors on the Pattern of Fracture
Shapes. e geometry of the hydraulic fractures in CBM reservoirs is influenced by many factors; it is almost impossible to predict exact shapes of the fractures after the hydraulic sand fracturing operation. Fortunately, the evaluation of the fractures does not require exact prediction and correct recognition of the pattern of fracture shapes sometimes is enough for the engineers to make the design of refactoring operation. According to the reservoir stimulation mechanics [2] and the field experience, the pattern of the fracture shapes is mainly influenced by the following features.

2.2.1.
e Types of Hydraulic Sand Fracturing Operation Curves. By the analysis of hydraulic fracture operation curves (variation of pumping pressure, etc.), we can determine qualitatively how difficult the fractures extend in the coal seam. e reason is that the variation of pumping pressure is the result of the compound effect of the ground stress field, the flow friction resistance in fractures, and the pressure in wellbore. In general, the hydraulic fracture operation curves can be divided into 4 types; see Figure 5.
Similarly, these 4 curve types represent four cases of the fracture propagation: Here, we should notice that it is easy to classify the typical operation curves. But in most cases, the type of the operation curves is discriminated subjectively by the fracturing engineers according to their experience. us, it is meaningful to make this classification by machine learning methods. Figure 6, once the proppant ratio (the black line) gets larger than 0, the fracturing operation turns to pump the sand-laden fluid (the 4 Complexity    Complexity yellow box). e ideal case is that the proppant (sand or ceramsite) is distributed uniformly on the fracture planes, and the fracture is effectively propped. Here, the sensitivity to proppant ratio reflects the fracture width and the sand distribution. Generally, if the fracture width is not enough, the proppant is easy to contact the fracture plane, and then fracture is sensitive to proppant ratio. e sensitivity to proppant ratio can be discriminated through the variation of tubing pressure. If the tubing pressure varies along with the variation of proppant ratio (positively correlated), the fracture is sensitive to proppant ratio. Otherwise, it is not sensitive. For example, in Figure 6 the tubing pressure increases along with the increasing of proppant ratio, so we can say that the fracture is sensitive to sand ratio. Similarly to the former subsection, the discrimination of the sensitivity to proppant ratio is often subjective and fuzzy; thus, it is also necessary to find a machine learning method to solve this problem.

Pressure Decline Rate after Stopping Pump.
e pressure drop curve is change curve of the bottom hole pressure or the well head pressure after stopping pump (see Figure 7, in the pink box). During this period, the fracturing fluid will leak off to the formation by the inner and outer pressure difference of the fractures. us, we can analyze the filtration of the formation according the pressure decline rate (see Figure 7).
In Figure 8, we can stiffly divide the pressure drop curves into 3 levels, the curves of high, moderate, and low decline rate, according to the field experience on the fractured CBM well. For example, by defining the pressure decline rate R d as the average pressure drop in 60 minutes after stopping pump, we take that if R d is less than 0.1 MPa·min − 1 , then the decline rate is low; if R d is between 0.1 and 0.15 MPa·min − 1 , then the decline rate is moderate; and if R d is larger than 0.15 MPa·min − 1 , then the decline rate is high.
As is known to all, the concepts "High", "Moderate", and "Low" are quite subjective and fuzzy, and it is not appropriate to classify them stiffly. e methods such as fuzzy comprehensive evaluation might be more suitable if we only focus on this point. However, the pressure decline level is just an intermediate (hidden) feature in our method; we can take the pressure drop curves as a sequential input of the related Bi-LSTM3 in Figure 2, the output of which is the numerical feature of the pressure drop curves and also the input of the classification layer.

Manner of Fracture Propagation.
e manner of fracture propagation is determined through the slope of the double-logarithmic curves of net pressure and time (lgP-lgt). Here, the slopes of lgP-lgt curves are obtained based on the classic PKN or KGD models [10,11].
During the fracturing operation, we can measure the wellhead pressure P surface . To get the net pressure P net in the fractures, we need to transform P surface to bottom hole pressure, and take the fluid column pressure in the wellbore P hyrostatic , the friction loss in the pipeline P fluidfic , the friction in perforation tunnels P pf , and the fracture closure stress (nearly equals to the minimal principal horizontal stress σ min ) into calculation. at is, One can see the Appendix for the details of the calculation of each term in (1). By (1), we can get the doublelogarithmic curve of net pressure to time (lgP-lgt) for each fracturing operation; see Figure 9.
According to Nolte's method [10], we know that the slope of lgP-lgt curve can reveal the efficiency of the fracture propagation during the fracturing operation. Based on the theory of reservoir stimulation mechanics [2] and the analysis on the historic materials of the discussed coalbed gas field in this paper, the manner of fracture propagation can be    divided into 6 types, which correspond to different intervals of the slope of the lgP-lgt curve; see Table 1.
Because of the uncertainties underground, the manner of fracture propagation may vary in different stages of the fracturing operation process. e pattern of the fracture shapes is influenced by the fracture propagation manners of the whole operation process. It is not easy to find suitable static features to characterize the lgP-lgt curve, and that is why we choose Bi-LSTM to analyze this kind of sequential data.

The Hierarchical Bi-LSTM Model to Recognize the Pattern of Fracture Shapes in CBM Reservoirs
e whole framework of the presented method in this paper can be divided into 3 steps: (1) data preprocessing; (2) feature extraction by the Bi-LSTMs, and (3) make the classification through the classification layer. e detail is listed in the following subsections.

Data Preparation.
First of all, the inputs of the Bi-LSTM layers in Figure 2 correspond to different data sequences, which should be extracted from the raw data as follows: (1) For the Bi-LSTM1, it is to extract the implicit features to demonstrate the types of hydraulic sand fracturing operation curves; thus, the input data composes of 3 sequences (the tubing pressure, the proppant ratio, and the injection rate), which are truncated at the stopping pump point (the injection rate is 0). (2) For the Bi-LSTM2, it is to discuss the correlation between the tubing pressure and the proppant ratio; thus, the input contains 2 sequences of tubing   pressure and proppant ratio in the sand-laden stage (the time interval is selected from the point that the prppant ratio is larger than 0, to the point that the injection rate goes to 0). (3) For Bi-LSTM3, the input is sequence of the tubing pressure after stopping pump (the length of the time interval is 60 minutes). (4) e lgP-lgt curve obtained by (1) is the input of Bi-LSTM4.
Generally, the raw data of the fracturing operation curves is sampled second by second, which is noisy and too dense to be the input of the deep network. us, as a start, the input sequences of the 4 Bi-LSTM layers should be filtered and smoothed to show clearer variation trends. In this paper, we employ the wavelet transform to denoise the data by choosing an appropriate wavelet basis [24]. en the smoothed sequences are sampled by minutes.

Classification through the Hierarchical Bi-LSTM
Networks. As a modification of RNN, the LSTM network is much more robust to gradient exploding or vanishing problems [28]. is network is capable of learning long-term information from large variety of sequential data. e basic structure of LSTM memory block contains 4 components: input gate, output gate, forget gate, and memory cell; see Figure 10. e gates and cell state of a LSTM unit are updated by the following formulae.
where x t , h t ∈ R d are the input and output of the network at time t. e matrices W and R with different subscripts are the weight matrices of the input part and the recurrent part respectively. b i , b o , b f , and, b c are the bias vectors with respect to the 3 gates and memory cell. σ, g, and φ are the activation functions (sigmoid or tanh). In this paper, we use sigmoid functions. ⊙ is the denotation of pointwise multiplication.
As a variant of LSTM, Bi-LSTM has been proved to be more effective as it can analyze information from the past and future [36]. e basic idea of Bi-LSTM is to use two parallel layers of hidden units to capture the past and future information. e two layers are concatenated to form the output. Comparing with LSTM, the output formula (7) should be adjusted to en (8) and (9) are concatenated to the final output: e first reason that we chose the hierarchical structure of the Bi-LSTM layers in this paper ( Figure 10) is that the pattern of the fracture shapes is the compound result of the 4 factors in Section 2. e data sequences related to the influence factors are interdependent with each other. It has been proved that the hierarchical Bi-LSTM network outperforms nonhierarchical ones to model the interdependencies of the sequences [29]. e second reason is that after the network is trained, we can analyze the implicit information (the influence factors of the pattern of the fractures) by unfolding the network [20,28]. e framework of the hierarchical Bi-LSTM is demonstrated in Figure 11.
In Figure 11, the fully connected layer is to embed the features into the neural network hidden layers, which is to combine all the extracted local information to get the global results in the next step. e final classification is implemented by using the softmax function [28].
where h i stands for the output of the network from the input x i , and j is the index of the label of the 4 patterns; thus, h j i means that h i is assigned with the label j. en the total loss function is defined as follows: where 1 · { } is the indicator function that equals one if h i � j and zero otherwise; h i is the true label of the input x i ; m is the number of the samples.

Dataset.
e proposed method is implemented on a dataset, which is composed of the hydraulic fracturing operation data of over 1101 CBM wells in Shanxi Qinshui Basin. e labels of the pattern of the fracture shapes are estimated according to the experience of the engineers on site and the simulation by Fracpro PT.

Training Details and Implementation.
As the max length of all the input sequences in Figure 9 is no more than 120 (2 hours), each Bi-LSTM layer in this paper consists of 240 units, the other training details are given by following the work of Li [33] and Sutskever [37]. e parameters of the LSTM units and the related neural network layers are initialized from a uniform distribution between [− 0.08, 0.08]. e stochastic gradient-based optimization based on 8 Complexity adaptive moment estimation (Adam) is employed in this paper to optimize the loss function [38]. Here, the learning rate is firstly fixed to be 0.001. e gradient is clipped when the norm exceeds a threshold of 1.

4.3.
Results. e pattern recognition of the fracture shapes in CBM reservoirs is a novel idea, and until now, no deep learning method is applied to this topic. us, the performance of our method is firstly tested by using the trained Input gate  Complexity network to make classification on the test data. In this paper, 80% of the data is applied to train the given hierarchical Bi-LSTM network and 20% is left for testing.
Secondly, we compared the performance of the LSTM network with different structures; see Table 2. Here, the abbreviation H-Bi-LSTM stands for the hierarchical Bi-LSTM network. From Table 2, we can see that the presented H-Bi-LSTM in this paper holds the highest accuracy of 81.8%, which is not that high as we hope (>95%). e reason is that the labeling of the training and testing data is done through the experience of the engineers and the Fracpro PT simulation, which has subjective error.
Finally, for some practical CBM well in China, inputting the data in Figure 12 and Table 3 into the presented network, we found that the predicted pattern of the fracture shapes after this fracturing operation is "High-Short".
Here, we can see that the simulated hydraulic fracture of this CBM well is indeed a High-Short fracture.

Conclusion
e hydraulic fracture propagation in CBM reservoirs is almost impossible to predict exactly, as it is a complicated nonlinear dynamic modeling problem, which is easily influenced by uncertain factors. In this paper, we explore a new idea to make the pattern recognition of the fracture shapes. In this idea, pattern recognition is regarded as understanding what the CBM reservoirs say during the hydraulic fracturing. e "words" of the CBM reservoirs are the related operation curves. en the methods of deep learning can be applied. e work of this paper is a part of the expert system, "Management decision software for low-productivity and low-efficiency coalbed methane Wells" (MDSCBM). Comparing with the conventional methods, the presented method in this paper takes full advantage of the nonlinear characterization capacity of deep networks and the processing capability of sequential information of LSTM, which provide a new perspective to predict the pattern of the fracture shapes. is method takes the experiences of the engineers and the measured data both into consideration, which can effectively assist the engineers on site to make the refracturing design. Appendix e calculation of the terms on the right side of formula (1) in the main body of this paper is listed as follows.
e fluid column pressure in wellbore is calculated through the following formula: P hyrostatic � ρ fluid gh, (A.1) in which ρ fluid is the density of the fracturing fluid, kg/m 3 ; g is the acceleration of gravity, m/s 2 ; h is the liquid height, m. e friction loss in pipeline P fluidfic is calculated by Lord DL's method [39]:    in which σ is the antifriction ratio, dimensionless; d is the inner diameter of the pipe, m; Q is the injection rate of fluid, m 3 /min; L is the length of the pipe, m; u is the average flow rage, m/s; G is the concentration of thicker, kg/m 3 ; P is the concentration of proppant, kg/m 3 . Next, the friction in perforation tunnels P pf is calculated by Willingham's method [40]: in which, Q is the injection rate of fluid, m 3 /min; ρ fluid is the density of the fracturing fluid, kg/m 3 ; d is the diameter of perforation tunnels, m; n is the number of perforation tunnels; C is the flow coefficient of perforation tunnels, C ≤ 8.9.

Data Availability
e data used to support the findings of this study have not been made available because they are confidential by the agreement with the CBM exploration company.

Conflicts of Interest
e authors declare that they have no conflicts of interest.