Monitoring the Geometry Morphology of Complex Hydraulic Fracture Network by Using a Multiobjective Inversion Algorithm Based on Decomposition

: The fracturing technique is widely used in many ﬁelds. Fracture has a greater impact on the movement of ﬂuids in formations. Knowing information about a fracture is key to judging its effect, but detailed information about complex fracture networks is difﬁcult to obtain. In this paper, we propose a new method to describe the shape of a complex fracture network. This method is based on microseismic results and uses the L-system to establish a method for characterizing a complex fracture network. The method also combines with decomposition to construct a new method called the multiobjective fracture network inversion algorithm based on decomposition (MOFNIAD). The coverage of microseismic monitoring results and the degree of ﬁtting of production data are the two objective functions of the inversion fracture network. The multiobjective fracture network inversion algorithm can be optimized to obtain multiple optimal solutions that meet different target weights. Therefore, this paper established a multischeme decision method that approached the ideal solution, sorting technology and AHP to provide theoretical guidance for ﬁnding a more ideal fracture network. According to the error of microseismic monitoring results, we established two cases of fracture to verify the proposed method. Judging from the results of the examples, the fracture network ﬁnally obtained was similar to actual fractures.


Introduction
With the depletion of global conventional energy, unconventional energy and renewable energy are gradually attracting attention [1,2]. Tight oil and gas reservoirs are an important part of unconventional energy. There are more than 900 tight oil and gas reservoirs in the United States, and more than 40,000 production wells are exploited. The porosity of tight reservoirs is less than 0.1, and the permeability of tight reservoirs is less than 0.1 mD. It is very difficult to exploit such extremely low-permeability reservoirs [3][4][5]. For this reason, horizontal well technology and fracturing technology are widely used in the exploitation of low-permeability reservoirs [6][7][8].
Fracturing technology was first applied in the field of oil and gas development in the 1930s. Successful fracturing measures can increase the production of reservoirs by up to 10 times, which is very attractive for oilfield development. Thousands of fracturing projects have been implemented every year [9]. Fracturing technology is not only applied to increase production in oil and gas fields, but also widely used in other fields such as the extraction of geothermal energy [10], the treatment of solid hazardous wastes [11], the measurement of ground stress [12], the repair of soil and groundwater layers [13], and so on. Information on fracturing fractures in these fields is key to evaluating the effect of fracturing measures.
Fracturing monitoring technology is an indispensable part of the fracturing process. The purpose of using fracturing monitoring technology is to evaluate the effect of current fracturing and to guide future reservoir development plans [14,15]. At present, fracture information is mainly obtained through fracturing monitoring technology, which is mainly divided into two categories: direct fracturing monitoring technology such as microseismic monitoring, tracer monitoring, etc., [16,17] and indirect fracturing monitoring technology such as well test analysis, production data analysis, etc., [18,19]. Core CT scanning technology has further confirmed that fractures have complex structures [20][21][22]. Both the distribution of natural fractures and the ground stress affect the complexity of fractures. Secondary fractures during the fracturing process are connected to natural fractures to form a complex fracture network [23]. Describing the shape of a complex fracture network has become a challenge.
Microseismic monitoring technology is widely used in reservoir fracturing [24,25]. However, its results can only monitor the approximate location of fractures and cannot accurately predict the shape of the fracture network [26]. Therefore, based on the general information on fractures provided by the microseismic results (because the previous research on microseismic technology is rich enough, this paper does not introduce it too much, as it used only microseismic monitoring data) [27][28][29][30], microseismic point sets were taken as the constraint of the fracture network, and the Delaunay triangulation method was used to divide those points into the possible network structures of the fracture network. Then, the L-system was used to generate the fracture network and establish the fracture network simulated model. A multiobjective optimization algorithm based on decomposition was used to obtain the inverse solution of the fracture network morphology [31,32]. The two objectives were the coverage of microseismic monitoring results and the degree of fitting to production data [33][34][35][36][37][38]. There are some similarities between fracture parameter inversion and history matching in oilfields, and these methods could be used for reference [39]. In essence, these challenges all belong to the category of optimization problems; therefore, we found some ideas from optimization methods [40][41][42][43]. The degree of fitting of the microseismic monitoring results was represented by the number of microseismic points covered by the generated fracture network, and the degree of fitting of the production data was represented by the absolute value of the difference between the production data of the simulated fracture network and the real fracture network.

Embedded Discrete Fracture Model (EDFM)
The commonly used methods to simulate fluid flow in reservoirs include the doublemedium model [44,45], discrete fracture model [46,47], and embedded discrete fracture model [48]. In the double-medium model, the ideal state is suitable only for reservoirs that have great connectivity between fractures. In the discrete fracture model, the reservoir is divided by an unstructured grid. When simulating the flow of some complex fracture models, the closer to the fracture and the larger the model calculation, the less sufficient the convergence becomes.
The embedded discrete fracture model overcomes the shortcomings of the doublemedium model and discrete fracture model. It not only highlights the flow characteristics of a single fracture but also greatly reduces the amount of calculation. Therefore, this paper selected the embedded discrete fracture module (HFM) in MRST (MATLAB Reservoir Simulation Toolbox) to simulate the production of fractured reservoirs [49,50].
EDFM has two main characteristics; one is the method of meshing, and the other is the method of coupling the fracture to the matrix [51][52][53][54]. The embedded discrete fracture model uses a structured grid to divide the reservoir matrix, while for complex fractures, it uses an unstructured grid. Figure 1 illustrates fractures using EDFM, and Figure 2 illustrates the results of the fracture in the above model.

The Delaunay Triangulation Method
The microseismic technique is the most commonly used monitoring method for fractured reservoirs. In this section, the microseismic point set obtained from the monitoring was taken as the constraint of the fracture network. The Delaunay triangulation method [55] was used to divide the microseismic points into possible fracture network structures.
The microseismic points were the set of subdivision points in the monitoring area. For reservoir areas not detected by microseismic monitoring, some regular point sets were defined, which are shown as the red points in Figure 3. The regular points and microseismic points outside the monitoring area constituted a set of points to be subdivided as a constraint of the fracture network inversion.  The irregularly distributed points in the middle are the microseismic points detected, and the microseismic areas are the areas where fractures were observed. The area of the peripheral regular point is the area where the fractures may not be detected by microseismic monitoring. The network structure after subdivision was regarded as a possible form of the fracture network, based on which a complex fracture network generation method was constructed in this paper.

Complex Hydraulic Fracture Network Generation Method Based on the L-System
Current studies have shown that fractures have the characteristics of tree-like fractals, and the L-system [56] is often used to simulate the shape of fractal trees. Therefore, a fracture network generation rule based on L-system was established.
To increase the complexity of the fracture network, we established a random mechanism of L-system character substitution in which a random number c was defined. Different character replacement rules were defined within the range of different random numbers. In each iteration, the random numbers were generated again. The character replacement rules were as follows: Figure 4 shows the fracture network morphology obtained through different iterations based on the previously established subdivided reservoir structure. It can be seen from the figure that the range and complexity of the fracture network increased with the increase of the number of iterations. It should be noted here that the fracture network generated by a high number of iterations is not necessarily larger than that generated by a low number of iterations, because repeated fracture segments appear in the final result generated by the fracture network. Therefore, repeated fracture segments in the fracture network were removed in the processing of the final fracture network generated, as shown in Figure 4f. However, a high number of iterations had a higher probability of generating a fracture network with a higher degree of complexity, and such results were easier to extend beyond the microseismic monitoring results in the reservoir. In choosing the next point in the process of seismicity, we considered the same fracture produced by the microseismic point on most of the relatively close, and we considered that microseismic points generated by different fracture segments may be far. We therefore established our model based on the choice of the distance criterion; i.e., points close to the current candidate microseismic point were selected during generation as being at greater risk of being fracture points; meanwhile, faraway points were less likely to be fracture points based on the candidate microseismic point. In the actual operation, the Euclidean distance between each candidate point and the current microseismic point was calculated, the reciprocal was taken, and this reciprocal of the distance was normalized into the probability for the selected candidate point.

The Bayesian Objective Function
Bayesian theory is a statistical prediction method that introduces the prior information of the problem into the method [57]. In the petroleum field, it is often used to construct the objective function in the historical fitting method. In the history-matching method, the sedimentary facies of the reservoir are used as the prior information of reservoir geology inversion [58]. Bayesian theory can be summarized as follows: Assuming that the problem is subject to a Gaussian distribution with a priori information as the mean value, then the probability p(m) of the problem is: When the parameter m is given, the probability distribution function of the observed data is the probability distribution function of the observation error: where m represents the desired parameter, which in this paper was the candidate option of subspace, namely the microseismic monitoring results and the production data; m pr represents prior information; C m represents the covariance matrix of prior information, and its diagonal is the variance of the desired parameter; C D represents the error of observed data; d obs represents real observation data; g(m) represents the forward modeling results corresponding to different parameters, which in this paper were the seismic data of the reservoir under the production of different fracture networks detected by microseismic monitoring. According to Bayesian theory, under the premise that the true model is that which actually occurred, the probability distribution function p(m|d obs ) of the request parameter was as follows: Networks of greater probability demand more closely to the true model parameters. Therefore, the inversion of the target maximizes the probability p(m|d obs ), so the objective function was simplified as: In the optimization process, the parameter that minimizes the objective function f (m) should be obtained. When the objective function reaches the minimum, it can not only make the observation data fit best, but also conform to the prior model.

Multiobjective Fracture Network Inversion Algorithm Based on Decomposition
A general multiobjective minimization problem that aims to simultaneously optimize more than one conflicting objective can be formulated as where x = (x 1 , x 2 , . . . , x d ) T is the optimization vector, X is the feasible region in ddimensional decision space, m is the number of objective functions, and F: Rd → Rm is an objective vector. Multiobjective algorithms are mainly used to solve optimization problems with more than one objective function, where objective functions are conflicting objectives [59][60][61][62]. Using a multiobjective optimization algorithm, solutions with different target ratios can be optimized. For a multiobjective problem with two objective functions, the feasible solution that satisfies the multiobjective problem constitutes the feasible region of the multiobjective problem. The feasible region is usually an irregular closed plane, and the boundary curve inside the feasible region near the ideal side of the multiobjective problem is called the Pareto front. The solutions on this front are optimal and nondominated. The dominance relationship is defined as where ≺ refers to the dominance relationship and x 1 is said to dominate x 2 . The optimization objective of the multiobjective algorithm is to find the Pareto front in the feasible region, and then select the required solution from the Pareto front according to different needs. Figure 5 illustrates the Pareto front for a biobjective problem. The two objectives used in this paper were the degree of fit between the fracture network and the microseismic detection results and the degree of fit between the fracture network and the real fracture production data during the production process. Because of the error of microseismic monitoring, the change of reservoir pressure during the production process, the ambiguity of fracture inversion, and other factors, the two objective functions were directly proportional in a nonstrict sense. The multiobjective optimization algorithm is roughly divided into two parts, one to maintain the diversity of the multiobjective problem solution and the other to maintain the uniformity of the Pareto front. In this paper, the L-system-based fracture network production method was used to maintain the diversity of solutions, and the decomposition method was used to maintain the uniformity of Pareto front solutions.
The decomposition method transforms the multiobjective problem into multiple scalarization optimization problems according to different objective weights so that the multiobjective algorithm can solve a finite number of solutions on the Pareto front using the decomposition method. The Chebyshev method was used as the decomposition method: where z i * represents the optimal objective function value of the objective i in all solutions and z * is a collection of the optimal values of each objective function.
The principle of the Chebyshev method is to minimize the maximum value of each target and optimal value under different weights. In the optimization process, the optimized solution continuously approaches the Pareto front around the weight vector. The Pareto front solution obtained by this method depends on the choice of weights.
The multiobjective fracture network inversion algorithm established in this section combines the fracture network generation method with the Chebyshev decomposition method. The algorithm allocates the number of subproblems and evenly distributes the weight of different objectives in each subproblem, sets the number of neighborhoods for each subproblem, and assigns neighborhoods according to the Euclidean distance of each subproblem. The neighborhood is set to reduce the amount of calculation for updating the subproblem in each iteration. The optimal solution in each subproblem is updated based on the information in the neighborhood, and each subproblem continuously approaches the optimal solution on the Pareto front through the Chebyshev method. For the calculation of the objective function, this paper took the degree of fit of the microseismic monitoring results and the degree of fit of the production data as its two objective functions. The number of microseismic sources covered by different fracture networks was used to indicate the degree of fit of the microseismic monitoring results. The embedded discrete fracture model was used to simulate the flow and production results of different fracture networks, and the absolute value of the difference between the production data of the calculated fracture networks and the real production data was calculated to indicate the degree of fit of the production data. The maximum number of iterations was used as the condition for the algorithm to stop. Figure 6 illustrates the algorithm flowchart.

Multiple Criteria Decision Making
Through the established algorithm, we could optimize to obtain a variety of solutions that met the proportions of different attributes. The multiattribute decision-making method was needed for the decision making of these schemes in different situations.
The TOPSIS (technique for order preference by similarity to an ideal solution) helps decision makers to choose the best solution in different situations by comparing the gap between different schemes and the most ideal scheme and the worst scheme. This method makes full use of the information of different attributes to help decision makers quickly find the best solution, and there are fewer restrictions on attributes. Figure 7 illustrates the principle of this method. The maximization of the properties of the problem is in the ideal state. Compared with scheme 2, scheme 1 is closer to the positive ideal solution, so scheme 1 is more applicable than scheme 2.  The analytic hierarchy process (AHP) is a multiple-criteria decision-making process that is used to set priorities among different attributes. The AHP has been widely used to reflect the importance, or relative weight, of the factors associated with priorities. The AHP method consists of three basic steps: first, the problem is broken down and structured into a hierarchy of subproblems; second, the data are collected and measured through pairwise comparisons of the attributes; and finally, the priority weights of factors or items in each level are calculated. A judgment or comparison consists of a numerical representation of a relationship between two elements that share a common parent. The set of all relative comparisons is then reported in a square matrix in which the elements are compared with each other. Figure 8 illustrates a schematic diagram of the hierarchy structure in AHP. In this paper, considering the advantages and disadvantages of the ranking method that approaches the ideal solution, a multiattribute decision-making method that combined the TOPSIS and AHP was established [63].
To reduce the judgment error of the TOPSIS, AHP analysis was used for the sorting method of approximation to the ideal solution scheme. The weights of attributes in the AHP were used to construct the hierarchical structure of the model and to set the target for calculating the weights of every scheme as close as possible to the ideal solution layer in the sorting method solution properties. In this paper, the objective function was used in the multiobjective fracture inversion algorithm, and the criterion layer was defined as the index. The latter was selected to make it easy to judge the weight of the scheme layer, so that the decision maker can more accurately select a scheme without rich experience and with minimal use of subjective judgment. Figure 9 illustrates the flowchart of the established decision-making method. Multiple fracture network schemes optimized by the multiobjective fracture network inversion algorithm were taken as the input of candidate schemes. The judgment criteria and the objective function used in the multiobjective fracture network inversion algorithm were defined as the inputs of criteria and scheme attributes.

Reservoir Model
A 600 × 600 m low-permeability, heterogeneous oil reservoir model was established. This oil reservoir model was developed by using the five-point method. Four injection wells are located in the four corners of the reservoir, and one fractured production well is located in the center of the reservoir. Figure 10 shows the well locations in the reservoir and the permeability of the reservoir. The modeled reservoir has adopted fracturing and stimulation measures. The simulated reservoir has been produced for 5 years after fracturing. The embedded discrete fracture model was divided into 68-step simulations. The reservoir simulation was carried out with constant production bottom flow pressure. It was assumed that the formation contained only oil and water two phases. The specific reservoir parameters are shown in Table 1. The difference between case 1 and case 2 lay in the injection rate of the injection well. The EDFM was used to simulate the production of the fractured reservoir.

Case 1
For the shape of the real fracture, the method of using the L-system to characterize fracture morphology in the references was used to randomly generate a real fracture. Microseismic points with smaller errors reflected the approximate distribution of fracture. Smaller interference was added that conformed to the normal distribution, and regular point division as added.
Applying the above generation method to independently generate fracture networks of different shapes on both sides with the production well in the center, the actual fracture, shown in Figure 11, was generated. A point was randomly generated in the generated real fracture segment to indicate the location of the microseismic point, while the midpoint of each fracture segment represented the distribution of the microseismic points generated by the simulation. To consider the errors in the microseismic monitoring process, disturbances that conformed to the normal distribution were added to the generated microseismic point distribution to simulate the set of microseismic points in the fracturing area. The distribution of simulated microseismic point sets and the distribution of microseismic point sets after adding errors are shown in Figure 12. The cross marks represent the initial generated microseismic point sets, and the dots represent the distribution of microseismic point sets after adding random interference. Based on the generated microseismic points, regular points were added to the external area, and the points existing in the reservoir were Delaunay triangulated. The interval of adding regular points was 40 m. The shape after division was as shown. The red dots indicate the added regular points, and the black dots indicate the distribution of the set of microseismic points after adding the error.  After dividing the microseismic points and the added regular points in the reservoir, using the coverage of the microseismic monitoring results and the fitting degree of the production data as the objective function, the MOFNIAD method proposed in this paper was used to invert the fracture morphology. The coverage of the microseismic monitoring results was expressed by the number of microseismic points covered by the generated fracture network, and the fitting degree of the production data used the absolute value of the difference between the production data of the simulated and real fracture networks. For the production data of the real fracture model, the embedded discrete fracture model introduced previously was used to simulate the production of the real fracture network. The saturation field and production data of the reservoir at a certain time after fracturing are shown in Figure 13. The number of multiobjective subproblems was 50, and the maximum iteration number of the multiobjective inversion algorithm was 100. The Pareto front obtained by the final multiobjective inversion is shown in Figure 14. In this example, due to the small error of microseismic monitoring results, the number of microseismic points covered by the actual fracture morphology was highest, and the fitting degree of the microseismic monitoring results was basically proportional to the fitting degree of the production data. Therefore, the Pareto front was narrow, and only three schemes that were similar to the real fracture were optimized. Figures 15-17 show comparisons of the fracture morphology and production data of these three schemes.    To further evaluate the similarity between the fracture network and the real model in the optimization scheme theoretically, this paper defined the allowable error range. The fracture inversion method is a fracture network constructed on the basis of the Delaunay triangle. If the fracture segment obtained by the inversion belonged to the triangle where the real fracture segment was located, the fracture segment of the inversion was considered as within the acceptable error range. Figure 18, in which the solid line is the real fracture segment and the dotted line is the fracture segment within the acceptable error range, shows the acceptable error range. Figure 19 shows the distribution range of fracture segments within the allowable error of the theoretical optimal model fracture network.   Figure 20 shows the error area evaluation results based on the fracture networks of schemes 1, 2, and 3, respectively. The red line represents the coincident fracture segment with the theoretical optimal fracture network, the orange line represents the allowable error, and the green line represents that the error was larger than allowed. The error evaluation results of the three schemes based on the theoretical optimal fracture network are shown in Figure 21, where the red and orange lines' meanings are the same as those in Figure 20. The blue line represents the fracture segments in the theoretically optimal fracture network that were not obtained by inversion. The fracture network in scheme 1 was basically within the allowable error range. However, judging from the evaluation results based on the theoretical optimal fracture network, many fracture segments were not within the allowable error range. In scheme 3, the fracture network had more fracture segments within the allowable error range of the theoretical optimal fracture network, but the number of redundant fracture segments was also large. The fracture network in scheme 2 is between scheme 1 and scheme 3. Most of the three optimized fracture networks were within the allowable error range and were basically close to the theoretical optimal fracture network form.

Case 2
Research has shown that in the technology of microseismic monitoring, the process of picking up microseismic signals, obtaining velocity models, constructing forward models, constructing objective functions, and selecting optimization algorithms brings errors to the source location; these errors can be as much as 100 m. For this reason, it is necessary to consider the large error of microseismic monitoring results. Case 1 used microseismic monitoring results with a small error to verify the fracture network inversion method proposed in this paper. Case 2 considered another case where the error between microseismic monitoring results and real fracture morphology is large. The feasibility of applying MOFNIAD was considered. In Case 2, the morphology of the actual fracture was basically the same as that in Case 1, but the strike of the fracture was adjusted. Meanwhile, the microseismic points close to the fractured well were added with less interference, while the microseismic points set far away from the fractured well were added with more interference. Figure 22a shows the real fracture network. Figure 22b shows the comparison of the distribution of microseismic points without interference and after interference. The blue dots represent points without interference, and the black dots represent points after interference.  Figure 23 shows the obtained Pareto front. In this case, the microseismic monitoring results have a large error, so the coverage degree of the microseismic monitoring results and the fit degree of the production data are two conflicting objectives. The greater the coverage of microseismic results, the greater the gap between the fracture network and the real fracture network, and the more obvious the difference in the degree of fitting of production data. When the microseismic monitoring results have a poor fit, the production data may fit better with the actual production data. The Pareto front obtained in this case was wider than that in the previous section, and the optimal solution had a larger range and more numbers on the two optimization objectives. The schemes were named in the order of the abscissa on the Pareto frontier. The model with the best production data fit was the first scheme, and the tenth scheme had the best coverage of the microseismic monitoring results. In this case, there was no theoretical optimal model that dominated all solutions. When the inversion fracture is very close to the real fracture, the obtained production data has little difference from the real production data. In other words, the objective function of production data fit is smaller. This scheme has a dominant relationship with the solution within a small range of production data. Therefore, when the number of iterations is large enough, the scheme that is theoretically closest to the real fracture shape is the scheme with the smallest abscissa on the Pareto front (the first scheme). The higher the degree of microseismic coverage, the greater the possibility of deviation from the true fracture shape. Figure 24 shows the fracture morphology in the first and second schemes. In schemes 1 and 2, due to the better fit of the microseismic monitoring results, the fracture network shape extended to the microseismic monitoring area, slightly deviating from the real fracture shape. Figure 25 shows the fracture shapes in schemes 3 and 10. The fracture network form in scheme 10 deviates significantly from the real fracture network form.  From the comparison between the production data in Figures 26 and 27, since the area covered by the fracture in scheme 10 was larger, the simulated oil production was higher than the real production data.  Thus, it can be seen that the scheme theoretically closest to the real fracture network was scheme 1; it fits the production best. The allowable error range defined in the previous case was applied to evaluate the similarity between fracture morphology and the real fracture network in scheme 1. It can be seen from Figures 28 and 29 that most of the fracture segments of the fracture network in scheme 1 were within the allowable error range.  The optimization results of this case had a scheme that was far from the real model, such as scheme 10 with a larger coverage of microseismic results. Therefore, a multischeme decision-making method was used to choose a better scheme.
When the two objectives of optimization are proportional and the number of algorithm iterations is large enough, the optimized Pareto frontier is narrow and there are fewer solutions. If the two objectives are conflicting, as in this case, when the number of iterations of the algorithm is sufficient, the optimized Pareto front range is wider and there are more solutions. The following figure shows a comparison of the Pareto front (triangular mark) obtained from the case of microseismic monitoring with small errors and the Pareto front (dots) from the case of microseismic monitoring with large errors discussed in this section. Therefore, according to the multischeme decision-making method established in this article, one can first judge the weighting tendency of the optimization objective based on the optimized Pareto front. If the Pareto front range is wide, the relationship of the optimized targets is obviously in conflict. When using the analytic hierarchy process to determine the target weight, the accuracy or attribute confidence should be more inclined to the production data. In practice, the error in the production data monitoring results is small.
The above section shows a case of a multischeme decision-making method. It inclined toward production data in terms of accuracy and attribute confidence. The weights of the fitting degree of production data and the degree of microseismic coverage calculated by the AHP method were (0.78, 0.22), respectively. The best solution presented was scheme 1.

Conclusions
This paper proposes a new fracture network inversion method by using a decompositionbased multiobjective fracture network inversion algorithm. According to the error conditions of microseismic monitoring results, two different test cases were established. The MOFNIAD was used in these cases and obtained similar schemes to the real fracture network shape. If there were many optimization result schemes, the method was inclined to produce data using multiple decision methods. The results of the cases verified the characteristics and inversion effects of the fracture inversion method proposed in this paper.