Abstract
In order to avoid the occurrence of caprock integrity damage and gas escape due to injection pressure overrun in CO2 sequestration, an optimized back propagation (BP) neural network model based on Monte Carlo simulation (MCS) and an improved lion swarm optimization (ILSO) algorithm is proposed for the maximum sustainable injection pressure prediction. Firstly, a hydromechanical model is constructed to simulate the damage changes of the reservoir caprock during injection by ABAQUS. Secondly, in view of the uncertainties of formation parameters that could lead to deviations between the model calculation results and actual geological condition, the MCS method is used, and then, the probability distribution interval of the maximum injection pressure with high probability of caprock failure under different formation parameters is obtained by MATLAB. This effectively reduces the uncertainty and improves the calculation accuracy. Finally, based on the numerical simulation results, the maximum injection pressure prediction model is constructed. Aiming at the problems of the sensitivity of the BP neural network to initial weights and its poor convergence, tent chaotic mapping and difference mechanism are introduced to improve the LSO algorithm. Following this, the neural network is optimized by ILSO algorithm whose superiority is verified through 8 benchmark functions, and the maximum injection pressure is effectively predicted according to porosity, permeability, and other parameters. Experimental results show that, compared to the other three optimization methods, the ILSO_BP model has a faster convergence speed, higher prediction accuracy, and stability, which can provide a powerful guide for the safe injection of CO2 and efficient sequestration management.
1. Introduction
High CO2 concentration in the atmosphere can lead to a series of problems such as global warming and extreme climate, thus posing a serious threat to the sustainable development of the ecological environment. Geological CO2 sequestration is one of the most scientific and effective means for controlling the CO2 concentration in the atmosphere, which is of great significance to mitigate the greenhouse effect, improve oil and gas recovery, and promote the development of low-carbon economies [1]. The injection of large amounts of CO2 into the geological reservoir will break the equilibrium in the reservoir and have an important impact on the mechanical stability of the caprock [2] (Figure 1). Once a leak occurs, it will likely cause a serious environmental disaster. Thus, it is necessary to quantitatively characterize the extent of the caprock damage during the injection. Many complex factors affect the caprock integrity in CO2 sequestration projects, such as the injection rate and material properties [3–5]. Among them, the injection pressure is the key external factor triggering the geomechanical effect, whereas the mechanical properties of rocks are the objective internal factor affecting the sealing performance. Moreover, the injection pressure directly affects the efficiency, storage capacity, and potential hydraulic fracture risk. Therefore, exploring an efficient and accurate prediction and evaluation method for calculating the maximum sustainable injection pressure (hereafter referred to as the maximum injection pressure) is a particularly important basis for ensuring the safety of the CO2 sequestration project, improving the storage efficiency, and reducing the cost.
Several scholars have carried out studies on the maximum injection pressure in CO2 sequestration, including theoretical analysis methods [6, 7], laboratory tests or gas injection tests in shallow sites [8], and numerical simulations [9]. Among them, the numerical simulation method can efficiently make up for the shortcomings of other methods (Figure 2) and has been widely used in many aspects such as potential assessment of CO2 sequestration sites [1], analysis of the formation pressure change during injection [10], leakage risk assessment [11], and analysis of the influence of injection parameters in CO2 enhanced oil recovery [12]. Geological CO2 sequestration is a complex thermal-hydrological-mechanical-chemical coupling process, and the constitutive models which are used for simulating this process are mainly of three types: elastic constitutive models, plastic constitutive models, and discontinuous media models [13]. The impact of different constitutive models on the simulation results of the caprock integrity analysis is different. Related researches show that the key to ensuring the integrity of the caprock is to make sure that the reservoir pressure does not exceed the breakthrough pressure of the caprock [14], whereas in engineering application, they often increase the injection pressure to improve the CO2 storage efficiency. The maximum safe injection pressure in the traditional method is obtained by multiplying the formation breakthrough pressure by a safety factor (0.8-0.9). However, this method can only roughly determine a range and cannot accurately predict whether the injection pressure will exceed the critical value and lead to the risk of rupture [15, 16]. In addition, the mechanical and permeability properties of reservoir caprock are obtained from the core test data, and porosity and permeability, as important parameters affecting the sequestration performance, present a complex multiparameter, nonlinear mapping relationship with the maximum injection pressure [17]. However, the measurement and sampling errors and the complexity of the geological environment lead to the incompleteness and some uncertainty in the logging data [18, 19]. The mean value of the formation parameters in the field is generally used in the model calculation, which leads to the deviations between the calculated results and the actual geological conditions. This, to some extent, reduces the reliability of the model and increases the risk of the sequestration project. Therefore, it is necessary to perform an uncertainty analysis on the relevant formation parameters to accurately predict the maximum injection pressure. However, at present, there are relatively few studies on this aspect.
Artificial neural networks (ANNs) have unique advantages in the analysis and processing of massive, high-dimensional big data and can efficiently characterize the complex nonlinear mapping relationships between the input and output data. As a result, ANNs have been widely used in mining the potential laws and predicting the unknown data in the field of oil and gas development, e.g., the prediction of the fracturing capacity of tight sandstone reservoirs based on the genetic algorithm-Elman neural network [20], the remaining oil plane distribution prediction based on machine learning [21], and the seismic reservoir prediction based on a convolutional neural network [22]. In particular, ANNs have also been widely used in CO2 sequestration engineering, e.g., deep learning-based reservoir mechanical response at the storage site [23], predictive characterization of the subsurface porosity by fusing ANNs [24], and ANN-based trap effectiveness prediction in CO2 sequestration in saline aquifers [25]. Related studies have achieved better results in the fine characterization of reservoirs, evaluation of CO2 storage safety, and prediction of the sequestration potential. However, the traditional ANN still has the problems of slow convergence speed and low accuracy, etc. Swarm intelligent algorithms have gained the attention of many researchers because of their fast convergence speed and better global search ability, and the lion swarm optimization (LSO) algorithm is a new type of intelligent algorithm that exhibits better global convergence and robustness compared to the particle swarm optimization (PSO) algorithm, grey wolf optimization (GWO) algorithm, etc. [26]. Therefore, in order to overcome the shortcomings of the back propagation (BP) neural network with a slow convergence speed and tendency to fall into a local optimum, this study combines the method of chaotic tent mapping and the difference evolution (DE) mechanism to improve the LSO algorithm and then proposes a CO2 maximum injection pressure prediction model based on BP neural networks optimized by the improved LSO (ILSO) algorithm.
This study was aimed at analyzing the variation characteristics of the formation stress field during CO2 injection through numerical simulation and constructing a model based on optimized neural networks for predicting the maximum injection pressure so as to improve the calculation efficiency and prediction accuracy. The contributions and innovations of the study are as follows. (1)For the study area, a hydromechanical coupled model has been constructed using the ABAQUS software based on the ideal elastoplastic constitutive model, the stress field distribution, and Mohr-Coulomb criterion. Then, the damage evolution law of the reservoir caprock under different injection pressures has been simulated and analyzed. Further, the maximum injection pressure corresponding to the initial damage of caprock has been obtained(2)Monte Carlo simulation (MCS) is used for analyzing the uncertainty of the formation parameters affecting the maximum injection pressure. Combining the numerical simulation results, the fitting of the relationship between permeability, elastic modulus, and Poisson’s ratio with the maximum injection pressure is performed by MATLAB so as to determine the probability distribution interval of the maximum injection pressure under different parameters, which effectively reduces the uncertainty and improves the accuracy of the model calculation(3)To address the problem of the LSO algorithm being prone to falling into local minimum at the later stage of the optimization search, the chaotic tent mapping strategy and the DE mechanism have been combined in this study, and the superiority of the ILSO algorithm has been verified by eight benchmark functions. Experimental results show that the ILSO algorithm enhances the uniformity and traversal of population distribution and improves the global search capability(4)The ILSO_BP model is proposed. Based on numerical simulation results, the optimized BP neural network has been constructed, and the ILSO algorithm has been used for optimizing the initial weights and thresholds of the neural network for improving the convergence speed and training accuracy. Subsequently, the proposed model is applied to predict the maximum injection pressure by porosity, permeability, and other related parameters. Taking mean square error (MSE), root mean square error (RMSE), and the mean absolute error (MAE) as the evaluation indices, the ILSO_BP model has been compared with the standard BP neural network model, the BP neural network model optimized by the LSO (LSO_BP) and GWO (GWO_BP), respectively. Results show that the model proposed achieves better performance in predicting the maximum injection pressure in CO2 sequestration, with a faster convergence speed and higher prediction accuracy compared to the other three methods
2. Numerical Simulation of Maximum Injection Pressure
ABAQUS is widely used for simulations in engineering fields for its powerful data analysis ability [27, 28]. In this study, based on the geological formation and hydrogeological conditions of the study area, the damage changes of the caprock and the maximum injection pressure has been analyzed by ABAQUS. In this section, we present the geomechanical model and numerical simulation.
2.1. Caprock Failure Criterion
Studies have shown that pore pressure has a direct effect on the deformation and damage of porous media [29]. The caprock is an elastoplastic body, which will break and weaken when the load reaches its yield strength, and its failure is controlled by effective stress [30]. When the reservoir is injected with a fluid or a gas, the general reduction in the effective stress might lead to the damage of the caprock and reduce its sealing capacity, thus triggering the risk of gas leakage [31]. When the effective stress is lower than the tensile strength, tensile failure will occur. Prediction of the tensile failure is relatively easy due to the relatively low tensile strength of most sedimentary rocks, which is usually assumed to be zero. In this study, based on the Mohr-Coulomb criterion and the maximum tensile stress theory, we have defined the caprock shear failure threshold, , as given by Equation (1), where and are the rock cohesion and the internal friction angle, respectively; and are the maximum and minimum effective principal stresses, respectively. When , the caprock belongs to stable state, whereas when , shear failure occurs. The smaller the value of , the more serious is the damage of the caprock.
In this paper, it is assumed that the reservoir skeleton deformation follows the Terzaghi effective stress principle, as shown in Equation (2), where is the effective stress, is the total stress, is the pore fluid pressure, is the Kronecker tensor, and is the Biot coefficient.
The stress balance equation is shown in Equation (3), where is the volume component.
The fluid flow in the reservoir is governed by Darcy’s law. Taking the influence of gravity and capillary force in the flow process into consideration, the gas-water flow equation can be deduced from the law of mass conservation.
The input parameters of the equations that ABAQUS solves include density, elastic modulus, Poisson’s ratio, permeability, and porosity. The output parameters include pore pressure, stress, and the shear failure threshold of the cap, namely, , etc.
2.2. Geomechanical Model
2.2.1. Geological Survey and Conceptual Model
Taking a subsurface aquifer in a block as the study subject, fully utilizing the measured data of cores, it is known that the cover layer is thickly bedded green and brown mudstone with a maximum single-layer thickness of 50 m, a stable lateral distribution, and good seal. The reservoir is buried at a depth of 1050 m, with a cumulative thickness of 100 m, dominated by fine sandstone, with an average porosity of 19% and an average permeability of m2, which has good physical properties. Considering the sealing safety, it is necessary to carry out the mechanical integrity analysis of the caprock for the injection process. The numerical model adopts the flow-solid coupling seepage model for single-phase fluids in porous media, and the three-dimensional situation of straight-well injection has been simulated by an axisymmetric method centered on the injection well according to the requirements for numerical calculations. The conceptual model is shown in Figure 3(a).
(a)
(b)
The constructed model extends vertically from the surface to a depth of 2500 m and includes four typical layers, namely, the overburden, the caprock, the reservoir, and the bedrock. No fault through the caprock is found at the top of the overburden, and the injection point is located in the middle and lower part of the reservoir at the symmetry axis. Based on the delineation criteria for finely portraying the vertical thickness and planar spreading scale of the reservoir, as well as the time and efficiency of the simulation calculation, the model has been meshed using an eight-node hexahedral displacement cell and a pore pressure cell. The grid has been divided equally according to the reservoir thickness in the vertical direction. A dense to sparse grid has been used in the horizontal direction, with a small-sized grid near the injection well and a gradually coarsening large-sized grid in the distant area, as shown in Figure 3(b).
2.2.2. Material Properties and Boundary Conditions
Without considering the heterogeneity, the reservoir caprock parameters were set to the general values of the study area, as shown in Table 1.
The internal friction angle, the maximum horizontal principal stress, and the minimum horizontal principal stress given in Table 1 are mainly used for the calculation of the in Equation (1) to judge whether the shear failure occurs in the cap. The maximum horizontal principal stress and minimum horizontal principal stress were derived from the test results of the differential strain in situ stress. The value of the internal friction angle came from the reference [32]. Ignoring the effect of boreholes, the reservoir was injected with gas at constant pressure, the inner boundary of the caprock was nonpermeable, the outer boundary was a constant pressure boundary, and the normal displacement around the model was fixed.
The detailed initial boundary conditions are as follows. The mechanical boundary is the upper boundary without constraint; the lower boundary is normal displacement constraint; the left boundary is axially symmetrical, and the right boundary is a normal constraint. The permeable boundaries are as follows. The upper boundary is the seepage boundary; the lower boundary is undrained; the left boundary is undrained, and the right boundary is a closed boundary. The initial conditions also include the initial pore ratio, the initial geostress field, and the initial pore pressure given before the model calculation.
2.3. Numerical Simulation of the Maximum Injection Pressure
Based on the ideal elastic model, a numerical simulation study of gas injection in the reservoir was carried out to analyze the damage change of the caprock under different injection pressures. For the sake of clear presentation, only the visualization results of the analysis area ( for the reservoir and for the caprock) have been presented.
2.3.1. Simulation Results Based on ABAQUS
The distribution of the initial ground stress was simulated before gas injection to bring the model to equilibrium. Subsequently, the reservoir was injected with CO2 at constant pressure, the simulation time was set to 1 year, and 10000 steps were set to calculate and obtain the pore pressure, horizontal stress, and vertical stress distribution before and after injection, as shown in Figures 4–6, respectively. For a clear representation of the temporal and spatial variations before and after the injection, the variation curve of the effective stress at the measuring point at the reservoir caprock interface with injection time is shown in Figure 7.
(a) Before injection
(b) After injection
(a) Before injection
(b) After injection
(a) Before injection
(b) After injection
Comparing the results of the temporal and spatial changes of the stress field and formation pore pressure in the reservoir caprock from the beginning of the simulation to the end of one year of injection, the following conclusions can be drawn. (1)The injection of CO2 has a significant effect on the reservoir pore pressure, which directly affects the effective stress of the formation and the mechanical integrity of the caprock. With CO2 injection, the formation pore pressure gradually increases, the effective stress significantly decreases, and the maximum stress reduction appears near the reservoir injection zone(2)After injection, the magnitude of vertical stress reduction was greater than that of horizontal stress reduction. One year after the simulated injection, the distribution trend of the pore pressure is consistent with the effective stress of the formation(3)With the increase of injection time, the vertical effective stress and horizontal effective stress at the measuring point gradually decrease, and the change rate of effective stress within one month of injection is large and then gradually decreases and tends to be stable(4)Under the influence of CO2 injection, the amplitude of the effective stress change at the bottom of the caprock is significantly smaller than that of the reservoir, and the maximum stress change occurs at the intersection of the reservoir and the caprock
2.3.2. Analysis of Reservoir Caprock Damage
Varying the injection pressure at the bottom of the well, the shear failure threshold, , defined in Section 2.1, was used to assess the damage to the reservoir caprock during CO2 injection. The value of was represented by SDV5 in the constructed model. Figure 8 shows the SDV5 nephogram of the damage occurred in the reservoir caprock in the equilibrium state before injection and 1 year after the injection at different bottom hole injection pressures.
(a) Before injection
(b) Injection MPa
(c) Injection MPa
As can be seen from Figure 8, when the injection pressure at the bottom of the well is 22.529 MPa, the SDV5 at the reservoir caprock junction is 1.001, and the caprock is still in a stable state. As the injection pressure increases to 22.530 MPa, the SDV5 reaches 1, and the caprock starts to undergo shear damage. In summary, with the increase of the injection pressure, the caprock gradually began to be damaged from the initial stable state, which could produce the risk of gas leakage. From the simulation results, it is known that the bottom-of-the-well maximum injection pressure for maintaining the integrity of the caprock is 22.529 MPa.
3. Uncertainty Analysis of Formation Parameters
When using ABAQUS to simulate the maximum injection pressure, the parameters used in the calculation model are empirical or average values of the study area, but measurement and sampling errors and the complexity of the geological environment can lead to some uncertainty in physical and rock mechanics parameters. For example, permeability is related to porosity, and neither has a fixed value but varies within a certain range depending on the actual geological conditions. However, a fixed value is used in the model calculations, which leads to uncertainty in the simulation. In order to obtain results that are closer to the actual geological conditions, the MCS method was used to perform the uncertainty analysis of the relevant formation parameters affecting the maximum injection pressure.
3.1. Uncertainty Analysis Methods
3.1.1. MCS
MCS is a numerical calculation method based on probabilistic statistical theory, which can accurately analyze the effects of various uncertainty factors to solve the probability of occurrence of a certain event or the expected value of a certain random variable. It is a relatively accurate method for the probability of failure calculation of multidimensional problems, which has been widely used in risk simulation analysis in many fields such as hydropower machinery [33] and geological engineering [34]. In this study, we use Python and MCS to perform an uncertainty analysis of the factors influencing the maximum injection pressure. The analysis contains four basic steps. (i)Describe the distribution type of parameters and construct a probabilistic model to make the solution correspond to certain characteristics of the random variables of the model(ii)Perform a MCS sampling based on the characteristics of the probabilistic model and the distribution of the random variables(iii)Conduct simulation experiments based on the constructed model to find the stochastic solution of the problem(iv)Statistically analyze the results of the simulation experiments to obtain a probabilistic solution as well as the cumulative probability
3.1.2. Determination of the Probability Density Function
Statistical information in MCS sampling shows that geological data are triangularly or normally distributed in most cases [35, 36]. The triangular distribution requires the determination of the minimum , maximum , and most probable values , as shown in Equation (4). In this study, the elastic modulus and Poisson’s ratio obeyed the triangular distribution as known from the field reports and core data of the study area.
If the random variable obeys a Gaussian distribution with an expectation value of and variance , denoted , it is modeled using a normal distribution with the probability density function given in Equation (5). In this study, it is known that the porosity and permeability conform to the normal distribution.
3.1.3. Prediction Model for the Failure Probability of Caprock
To simulate the probability distribution of the maximum injection pressure when the initial damage of the caprock occurs, the range of variation of each formation parameter was obtained by interpreting the logging data, and its probability density function was determined. Further, the influence of different values in the variation range of each parameter on the calculation results of the maximum injection pressure was analyzed by numerical simulation. The relationship between the maximum injection pressure and the formation parameters was established by polynomial fitting. On this basis, the prediction model for the failure probability of caprock was constructed, as shown in
In Equation (6), is the failure probability of caprocks, is a function related to the maximum injection pressure, and is the probability density function. Since it is difficult to calculate the failure probability using analytical methods, the MCS and the Latin hypercube sampling method were used for dividing each dimension of the 3D model uniformly into intervals that do not overlap. A point was randomly selected from each interval in each dimension to obtain subsamples, with each sample point denoted as , where is the number of sensitive parameters, and each sensitive parameter in the sample point has its own distribution law. Different sensitive parameters were used as model uncertainties to conduct a large number of statistical simulation tests to obtain the probability type of the maximum injection pressure distribution when the caprock failure occurs, namely, .
In the model calculation, we first calculate the sum of caprock failure probability of each sensitive parameter corresponding to the sample point, namely, . Then, take the mean value as the failure probability of the maximum injection pressure during the initial rupture under different parameters; thus, the distribution model could be constructed. In this paper, the simulation times is 10000 and there are three sensitive parameters; namely, the value of is set to 100 and is 3 in our experiment.
3.1.4. Uncertainty Analysis of Model Parameters
Based on the above model, the uncertainty analysis of the influencing factors of maximum injection pressure was carried out by performing the following steps. (i)The uncertainty of the physical parameters and rock mechanics parameters of the reservoir-caprock were analyzed, and the sensitive parameters affecting the maximum injection pressure were determined. In this study, we focused on three parameters, namely, permeability, elastic modulus, and Poisson’s ratio(ii)The influence law of the uncertainty of different parameters on the maximum injection pressure was analyzed, and polynomial fitting was done using MATLAB. In this study, the fit was characterized by the coefficient of determination () and the RMSE. The closer the is to 1, the better is the fit; the closer the RMSE is to 0, the better is the fit(iii)The probability density functions of different parameters and the number of simulations were determined, and the MCS-based parameter uncertainty analysis was implemented using Python. In this study, we know that the permeability conforms to a Gaussian distribution, whereas the elastic modulus and the Poisson’s ratio obey the triangular distribution(iv)The probability distribution, the cumulative probability graph, and the statistical results were plotted, and the probability distribution interval of the maximum injection pressure based on the influence of different parameters was obtained
3.2. Uncertainty Analysis of the Maximum Injection Pressure
The maximum injection pressure during CO2 injection is influenced by several factors. This section presents the uncertainty analysis of the formation parameters using MCS and MATLAB based on the actual geological condition.
3.2.1. Distribution of Formation Parameters
From the logging data and core analysis data, it is known that the third section of the Quantou group of the study area is a “medium-porosity, medium-permeability” reservoir, and its porosity and permeability distribution is shown in Figure 9, with the elastic modulus of 0.93–1.01 GPa and Poisson’s ratio of 0.23–0.36.
3.2.2. Analysis of the Effect of Different Parameters on the Maximum Injection Pressure
To investigate the effect of the parameter uncertainty on the calculation results, the influence of different values on the maximum injection pressure was analyzed using the MCS method within the range of three parameters, namely, permeability (), elastic modulus (), and Poisson’s ratio (). Further, the relationship curves between the maximum injection pressure and the corresponding parameters were plotted (as shown in Figure 10).
(a)
(b)
(c)
The maximum injection pressure versus , , and was obtained by fitting them using MATLAB as shown in Equations (7)–(9), respectively. The values of the three relations are 0.995, 0.9862, and 0.9956, respectively, and the RMSEs are 0.1147, 0.0073, and 0.0791, respectively, all of which are well-fitted. These relational equations can be used for predicting the maximum injection pressure at different values of , , and in order to fully assess the risk of caprock damage during injection.
3.2.3. Probability Distribution of the Maximum Injection Pressure
Based on the corresponding fitted relational equations and the probability density functions of , , and taken randomly, the calculation results with certain probability distribution were obtained by applying MCS sampling. The interpretation of the logging data showed that the reservoir permeability obeys a normal distribution (136.3193 and 137.8410), the elastic modulus obeys a triangular distribution (0.93, 0.97, and 1.01), and Poisson’s ratio also obeys a triangular distribution (0.23, 0.29, and 0.36). After determining the probability density functions of the three parameters, the maximum injection pressure was defined as the prediction parameter, and the probability distribution and the cumulative probability distribution of the statistical maximum injection pressure were calculated by simulation using the direct sampling method and setting the number of sampling times to 10000, as shown in Figure 11.
Since several formation parameters have uncertainties, the maximum injection pressure uncertainty analysis should be based on the distribution of each parameter for a comprehensive assessment of the caprock failure probability. As can be seen from Figure 11, the maximum injection pressure intervals corresponding to the 5% to 95% confidence level for , , and are (22.5344 MPa, 26.0737 MPa), (22.9343 MPa, 23.0171 MPa), and (22.9138 MPa, 24.6648 MPa), respectively, whereas the maximum injection pressure simulation results are 22.529 MPa, which is lower than the value corresponding to the confidence level of 5% for the three parameters, indicating that the possibility of the caprock failure is relatively small. The uncertainty analysis of the formation parameters can provide a more comprehensive and accurate assessment of the maximum injection pressure and provide a powerful guide for the safe injection management of the CO2 geological sequestration.
4. BPNN Based on ILSO Algorithm
The method of fixing other parameters and analyzing only the effect of a single factor on the maximum injection pressure has some limitations because of the correlation between certain formation parameters. In contrast, ANNs do not require a prior determination of mathematical equations for mapping relationships between inputs and outputs and have powerful nonlinear approximation capabilities. Thus, the BPNN was used to predict the maximum injection pressure in this study. To address the problems of the sensitivity of the BPNN to the initial weight and threshold, its slow convergence speed, and tendency to fall into a local optimum, the LSO algorithm was improved by combining the chaotic tent mapping with the DE mechanism, and the ILSO algorithm was applied to optimize the network weights and threshold, in order to improve the global search ability, convergence speed, and prediction accuracy of the model.
4.1. ILSO Algorithm
4.1.1. LSO Algorithm
The LSO algorithm is a swarm intelligent optimization algorithm proposed on the basis of the collaborative hunting strategies performed by lions [37]. In the LSO algorithm, the search range becomes smaller in the late stage of the optimization for the lion population individuals, thus making the algorithm prone to early convergence and local optimum, which reduces the diversity and the global search ability. To this end, an improved chaotic tent mapping sequence was introduced to initialize the population for its randomness, ergodicity, and regularity which can be used for improving the diversity and uniform ergodicity of the population distribution, inhibit the algorithm from falling into a local optimum, and enhance its global search capability.
Additionally, the mechanism of lioness position update in the traditional LSO algorithm exists only within the lioness population, and the effective information of the lion king and cubs might be wasted to a large extent and affect the convergence speed of the algorithm. Hence, it was improved so that when the lioness updates her position, she learns from the lion king and other cubs, in addition to communicating with the lioness group to enhance the global search capability. Further, a DE mechanism was introduced in the lioness position update to improve the accuracy while enhancing the local search capability of the algorithm.
4.1.2. Population Initialization Based on Chaotic Tent Mapping
A tent map involves segmented linear mapping, which is also a two-dimensional chaotic mapping with a uniform distribution function and good correlation in its parameter range. It is expressed by Equation (10). A tent map has uniformity and randomness, which can make the algorithm easily escape from the local optimal solution so that the diversity of the population can be maintained and the global search ability can be improved [38]. However, there are small cycles and unstable cycle points in its sequence. Thus, a random variable was introduced into the original chaotic tent mapping, and the improved chaotic tent mapping expression is given by Equation (11).
In Equation (11), is the number of particles within the sequence and is a random number ranging between [0,1]. In this study, we use the improved chaos tent mapping to initialize the population. By introducing the random variables, not only does it maintain the randomness, ergodicity, and regularity of the chaotic tent mapping but can also efficiently avoid the iterations from falling into small and unstable periodic points, thus solving the problem of falling into local minima during the training process and consequently strengthening the global search capability of the algorithm.
4.1.3. Lioness Position Update Based on DE
DE algorithm is an efficient global optimization algorithm with a simple structure, easy implementation, fast convergence, and robustness [39]. The original DE algorithm includes three operations of mutation, hybridization, and selection [40]. This study mainly uses the mutation operation of the DE algorithm to optimize the position update strategy of the lioness. Assuming that denotes the th particle in the th iteration, the mechanism of the mutation operation is expressed by Equations (12) and (13), where represent the different individuals in the th generation, is the scaling factor, and is the variance operator. On constant updating, gradually converges to , which can avoid the loss of the optimal solution of the parent generation and preserve it for the offspring.
In a D-dimensional target search space, a group of lions is formed, among which only one is a male lion, and the rest are lionesses with a certain proportion of adults. The position of the th lion is denoted as , . Different types of lions move their positions differently. The lion king moves in a small area at the best food to secure his privilege, and his position is updated according to Equation (14). The lioness requires the collaboration of other lionesses and her position is updated according to Equation (15). The process of the lioness position update communicates only with the lioness group and ignores the valid information of the lion king and cubs. To enhance the global search capability of the algorithm, the lioness position update strategy was improved by using the DE mechanism, as shown in Equation (16).
In Equations (14)–(16), denotes the position of the current lion in the th iteration; is the optimal position in the th generation; is the historical optimum in the th iteration of the current lioness ; is the historical optimum in the th iteration of a lioness randomly selected in the population of lionesses other than the lioness ; is a random number generated based on the normal distribution ; is a scaling factor with the same meaning as in Equation (13); the probability factor, , is a uniform random number generated based on the uniform distribution ; and is the lioness position update perturbation factor, which is defined by
In Equation (17), is the maximum number of iterations; is the current number of iterations; , where and are the maximum and minimum values, respectively, on each dimension. The perturbation factor, , makes the move step of the lioness change continuously with the number of iterations, first surveying with a larger step, followed by a change in the step size from large to small, and finally searching finely with a small step that converges to zero.
4.1.4. Performance Test of ILSO Algorithm
To verify the performance of the ILSO algorithm, eight benchmark functions were selected for the experiments, as shown in Table 2, where is the independent variable of the test function, is the independent variable dimension, represents the range of variables, and is the global optimal position. The search history of the ILSO algorithm on these 8 functions is shown in Figure 12, which demonstrates the global best position of each function; the historical positions of the lion king, lioness, and lion cub in each iteration; and the convergence curve of the optimal fitness.
Additionally, the performance of the ILSO algorithm was compared with PSO, GWO, and LSO algorithms. In these experiments, the same parameter settings were used for all algorithms; the population size was 50, the maximum number of iterations was 6000, and 50 independent runs on each test function, and the mean value (Mean), standard deviation (Std), and best value (Best) were used as the evaluation metrics. Experimental results are presented in Table 3, where represents the characteristic of functions, including unimodal separable (US), unimodal nonseparable (UN), multimodal separable (MS), and multimodal nonseparable (MN), and is the global optimal value. These eight functions are all representative, among which F5 to F8 are multimodal functions, F7 has several local extremes, and F8 is complex due to its uneven search surface, which makes the search process more complicated.
Experimental results in Table 3 and Figure 12 show that the ILSO algorithm can make the search process more refined and obtain accurate results by introducing the chaotic tent mapping strategy and the DE mechanism in the face of the more complicated search process and the difficulty of finding the optimal value. In addition, when the change of the search value falls into a long-time stagnation, the search value can still continue to fall, which indicates that the improved algorithm has the ability to jump out of the local extremes. Compared to the other three algorithms, the ILSO algorithm proposed in this study has better accuracy and stability.
4.2. ILSO_BP Model
To address the problem of the random initialization of weights and thresholds in BPNN which leads to poor generalization ability and its tendency of falling into local optimal value, the ILSO algorithm was used for optimizing the initial weights and thresholds in order to improve the convergence speed and accuracy.
4.2.1. BPNN
The BPNN is a multilayer feed forward neural network trained by the back propagation of errors. It uses gradient descent and search techniques to continuously adjust the weights and threshold, which minimizes the MSE between the actual output and the desired output. The method has a strong nonlinear mapping ability, high self-learning, self-adaptivity, and generalization ability. It has been widely used in the field of geological engineering for lithology identification [41], logging interpretation [42], and reservoir prediction [43]. The traditional BPNN is a local search optimization method, and the network weights are gradually adjusted by following the direction of local improvement, which easily leads the algorithm to fall into local extremes, thus leading to the failure of network training. In addition, because the search process of the BPNN uses the gradient descent method, when the objective function to be optimized is very complex, the algorithm will be inefficient and the convergence speed is often slow. Therefore, many researchers optimize BPNN by swarm intelligent algorithms in order to improve their convergence speed and global search capability and avoid falling into local optima.
4.2.2. Implementation of ILSO_BP Model
The process of the ILSO_BP model is shown as Figure 13.
The main idea of the ILSO_BP model is to introduce chaotic tent mapping and the DE mechanism in the LSO algorithm to improve the population initialization distribution and the lioness position update method, respectively, and apply the ILSO algorithm to the optimization of the initial weights and thresholds of the BPNN to enhance the global search capability, realize a more extensive and uniform search, and achieve fast convergence speed and high accuracy.
5. Prediction of the Maximum Injection Pressure
5.1. Model Parameter Setting
Experience shows that a three-layer BPNN can approximate a function of arbitrary complexity with arbitrary accuracy. Thus, in this study, a three-layer BPNN was used for predicting the maximum injection pressure for CO2 sequestration, and the initial weights and thresholds of the network were generated and trained by the ILSO_BP model. The input parameters are the four causeless variables affecting the maximum injection pressure, namely, porosity, permeability, elastic modulus, and Poisson’s ratio, and the output is the maximum injection pressure. According to Kolmogorov’s theorem, the number of nodes in the hidden layer was determined to be four, the transfer function of the intermediate layer neurons was Tanh, the initial learning rate was set to 0.01, the Adam algorithm was used for gradient update, and the number of iteration steps was set to 3000 and 6000, respectively. The population was set to 50, the proportion factor, , of adult lions was 0.2, and the value interval of weights and thresholds was [–5, 5].
5.2. Experimental Results and Analysis
The 140 sets of sample data obtained from the above ABAQUS simulation were selected, and the training and test sets were divided in the ratio of 8 : 2. MSE, RMSE, and MAE were used as evaluation metrics for model performance, as shown in Equations (18)–(20), where is the true value, is the predicted value, and is the number of samples. The proposed ILSO_BP model was compared with BP, LSO_BP, and GWO_BP methods. The iterative curves of the training errors of the four methods are shown in Figure 14. The final training error and testing error of each method are given in Table 4.
(a) 3000 times
(b) 6000 times
From Figure 14 and Table 4, it can be seen that the ILSO_BP model exhibits the lowest training errors as well as testing errors, a higher prediction accuracy compared to the other three algorithms, and the fastest convergence rate.
In addition, the errors between the actual values and the predicted values of the 28 data in the test set were counted, as shown in Figure 15. As can be seen, compared to the other three methods, the experimental results of the ILSO_BP model are all satisfactory, with fewer predicted outliers, and the overall error with the actual value is within a small range. The method proposed uses the ILSO algorithm based on chaotic tent mapping and DE mechanism in the optimization of the BPNN, which solves the defects of slow convergence and low training accuracy and achieves a higher calculation speed and accuracy in prediction of the maximum injection pressure.
(a) True value and predicted value
(b) Error
6. Conclusion
A coupled hydromechanical model was constructed using ABAQUS to simulate the mechanical changes of the caprock during CO2 injection, and the maximum injection pressure was predicted based on the ideal elastoplastic model, Mohr-Coulomb, and tensile failure criteria. Considering the influence of the uncertainty of formation parameters such as porosity and permeability on the reliability of the computational model, the uncertainty analysis of the relevant parameters was carried out using the MCS method to achieve a comprehensive assessment of the probability of caprock damage. Based on simulation results, a prediction model was constructed. The LSO algorithm was improved to optimize the structure of the BPNN by combining the chaotic tent mapping with the DE mechanism to address the problems of slow convergence of the BPNN and its tendency to fall into local optimal solutions. Experimental results show that the ILSO_BP model achieves better performance in the comprehensive evaluation of the maximum injection pressure for CO2 geological storage, and the convergence speed and prediction accuracy are higher than those of other three methods. Thus, the method proposed in this work can provide useful guidance for the injection management of geological CO2 storage projects, improve safety and storage efficiency, and reduce the risk of gas leakage. However, the actual geological situation is much more complicated, and further research is required to evaluate the integrity and sealing of the caprock during CO2 injection more comprehensively and accurately. In particular, the presence of faults in the sealing zone during injection will be an important research topic in the future. In addition, improving the model calculation speed further will also be the next work.
Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Nos. 42002138 and 42072166), the Postdoctoral Scientific Research Development Fund of Heilongjiang Province (No. LBH-Q20073), the Natural Science Foundation of Heilongjiang Province (LH2019F042), and the Excellent Young and Middle-aged Innovative Team Cultivation Foundation of Northeast Petroleum University (No. KYCXTDQ202101).