https://orcid.org/0000-0002-9556-1243
SUMMARY
In this study, seismic inversion is carried out using a genetic algorithm (GA) as well as a hybrid genetic algorithm (HGA) approach to optimize the objective function designed for the inversion. An HGA is a two steps coupled process, where a local optimization algorithm is applied to the best model obtained from each generation of the GA. The study aims to compare the qualitative as well as the quantitative performance of both methods to delineate the reservoir zone from the non-reservoir zone. Initially, the developed algorithm is tested on synthetic data followed by its application to real data. It is found that the HGA for synthetic data is providing more accurate and high-resolution subsurface information as compared with the conventional GA although the time taken later is less as compared with the former methods. The application to real data also shows very high-resolution subsurface acoustic impedance information. The interpretation of the impedance section shows a low impedance anomaly zone at (1055–1070) ms time interval with impedance ranging from (7500 to 9500) m s−1*g cc−1. The correlation between seismic and well data shows that the low impedance zone is characterized as a clastic glauconitic sand channel (reservoir zone). In seismic inversion using an HGA, one can delineate the areal extent of the reservoir zone from the non-reservoir zone more specifically as compared to the GA-derived impedance. The convergence time of HGA is 4.4 per cent more than GA and can be even more for larger seismic reflection data sets. Further, for a more detailed analysis of the reservoir zone and to cross-validate inverted results, an artificial neural network (ANN) is applied to data, and porosity volume is predicted. The analysis shows that the low impedance zone interpreted in inversion results are correlating with the high porosity zone found in ANN methods and confirm the presence of the glauconitic sand channel. This study is important in the aspect of qualitative as well as quantitative comparison of the performance of the GA and HGA to delineate sand channels.
1 INTRODUCTION
The objective of any geophysical experiment is to understand how a physical property of the subsurface leads to a particular observation and it explains why a certain type of observation is obtained. The process of solving the observed data for the possible physical model is an inverse problem and the process is known as inversion. All the geophysical problems are inverse and nonlinear. To solve these inverse problems, the most significant part is the construction of objective functions judiciously. Once the objective function is constructed, one needs an optimization technique to minimize/maximize the error to get the desired solution (model). An optimization technique requires a fitness function relatable to the objective function with several properties to facilitate the performance of the algorithm.
There are two types of optimization techniques used in geophysical problems, the first is local optimization [e.g. hill-climbing algorithm, Nelder-mead algorithm, quasi-Newton methods (QNMs), etc.] and the second is the global optimization method [e.g. genetic algorithm (GA), simulated annealing, etc.]. Most of the seismic inversion problems are solved by local optimization methods because of their fast convergence to the solution. But the biggest drawback of local optimization is that, it depends on starting guess model and hence it converges to the local optimum solution many times (Mallick 1995; Hejazi et al. 2013). If the solution will not converge to the global optimum, a false solution will be found. It has been two decades; since researchers have started using global optimization techniques to overcome the problem encountered in the local optimization case. The global optimization methods are more time-consuming as compared to the local optimization methods but provide the best solution to the problem (Sen & Stoffa 2013; Datta & Sen 2016). Recently, a new algorithm is developed to combine local optimization along with global optimization methods to converge toward solutions faster. This technique is called the hybrid optimization technique. In this study, two types of global optimization methods, namely GA and hybrid genetic algorithm (HGA) are utilized to demonstrate their capability to estimate subsurface acoustic impedance from seismic reflection data. In the hybrid optimization technique, we have used QNMs as local optimization methods under GA as a global optimization technique.
The global optimization algorithm finds the global minimum or maximum of an objective function. Many global optimization methods are generally stochastic, which uses global information about cost function to update their current state. They are not sensitive to the initial model this is how they escape from local minima even with a poor initial guess. A GA is a search-based optimization technique developed by John Holland (1973). It has an analogy with Darwin's evolutionary theory and is based on the principle of ‘survival of the fittest’. In the GA, first, the model parameters are encoded in form of a bit string. Further, from the sample, space population is generated and from it, individual fitness values are calculated and based on this fitness value, one individual solution is selected. After the selection process, two genetic processes, crossover and mutation are performed in a specific order to make a better solution (Chunduru et al. 1997; Maurya et al. 2020).
In an HGA, both local and global algorithms are used to optimize the objective function for faster convergence. Porsani et al. (1993) described it as a coupled two-step procedure in which a local algorithm is applied to the best model that is obtained from each generation of the global optimization algorithm, that is, GA in this case. Like a local optimization algorithm, it does not require a good initial model and a good solution can be obtained even with a poor starting model. It is less computationally efficient than global or local optimization schemes (Romero et al. 2000).
This study aims to compare the GA and HGA algorithm and their performance qualitatively as well as quantitatively to extract subsurface acoustic impedance models from seismic reflection data. To cross-validate inverted results, the artificial neural network (ANN) technique is also used for seismic data using inverted impedance as one input, and porosity volume is predicted.
The ANN approach, which has its origins in the hard-rock mining industry, is extremely beneficial in the exploration of geophysics for reservoir characterization (Bosch et al. 2010). There are many types of ANN techniques is available but in this study probabilistic neural network (PNN) technique is adopted as the method is quite simple and easy to apply to geophysical problems. In classification and pattern recognition applications, PNNs offer a scalable alternative to traditional backpropagation neural networks (Pramanik et al. 2004). They do not necessitate the huge forward and backward calculations that ordinary neural networks necessitate. When applied to a classification task, these networks use the concept of probability theory to reduce misclassifications (Chambers et al. 1994; Hampson et al. 2001).
This study focuses on two main issues: the first is to develop an algorithm for GA and HGA by modifying fitness functions developed by Ma (2002) and demonstrated their performance for synthetically generated data. Second, the performance of the algorithm to invert seismic reflection data is compared to the find best solution. To cross-verify inverted results, the PNN tool is adopted to predict porosity in the interwell region and compared inverted impedance with predicted porosity sections.
2 THE STUDY AREA: BLACKFOOT FIELD, CANADA
This study uses post-stack seismic data along with well-log data from the Blackfoot oil field, Alberta, Canada. The seismic and well-log data are recorded in 1995 by Consortium for Research in Elastic Wave Exploration Seismology (CREWES). The oil field is located about 70 km southeast of Strathmore city, Alberta, Canada (Maurya et al. 2020). In the Blackfoot oil field, the 3C-3D seismic survey has been done in two overlying patches. The first patch is the clastic Glauconitic channel and the second patch is targeted at reef-prone Beaver hill lake carbonates. Both patches are hydrocarbon producing formations since 1996. This study uses seismic reflection data along with well-log data from the first patch that targeted the clastic Glauconitic channel. Compare to the second patch of reef-prone Beaver hill lake carbonates, the first patch is of primary importance because from the geological information of the study area it consists of a clastic Glauconitic sand channel of lower Cretaceous age which represents a sediment-filled incised valley (Andrichuk & Wonfor 1954).
The Geology of the Blackfoot area is shown by a simplified stratigraphic column by Maurya & Singh (2015). Within the study area, the Manville Group of the Lower Cretaceous period rests unconformably over the Mississippian period carbonates of the Shunda formation. Shunda formation is swallowing up-section, fluvial and estuarine glauconitic sediments were deposited during subsequent regression. Glauconitic members in the Manville formation consist of medium to fine grain size quartz sandstone in Strathmore, Alberta. The Ostra codsbeds underlying the Glauconitic are made up of brackish water shales, argillaceous, fossil ferrous limestone and thin quartz sandstones and siltstones (Bell 1990). The Sunburst Member consists of sheet sandstones made up of sublithic and quartz arenites. The Detrital Beds make up the basal part of the Manville Group which has diversified lithology containing chert pebbles, lithic sandstone, siltstone and abundant shale.
In the Blackfoot field area, the Glauconitic sandstone is encountered at a depth of approximately 1550 m and the valley-filled sediments vary from 0 m to over 35 m in thickness. The Glauconitic Member is subdivided into three units. The lower and upper members are consisting of quartz sandstones with an average porosity, while the middle member is tightlithic sandstone. The individual members range in thickness from 5 to 20 m. Hydrocarbon reservoirs are found in structural and stratigraphic traps where porous channel sands pinch out against non-reservoir regional strata or low-porosity channel sediments. The primary hydrocarbon at the Blackfoot Field is oil, although gas may also be present in the upper member. Fig. 1 shows the location of the study area along with the clastic Glauconitic channel trend. The geological map and stratigraphic division of the Blackfoot oil field are presented in Figs 2(a) and (b), respectively.
Location of Blackfoot field, Alberta, Canada. The recorded seismic data is highlighted by blue colour marked with patches 1 and 2. The sand channel trend is shown with a dotted line, whereas well-log data are shown with a star (*).
(a) Depicts geological map and (b) stratigraphic division of the study area.
One inline ranging from 0 to 50 traces with two-way traveltime from 900 to 1100 ms time interval is chosen for the analysis. As the adopted methodology is more time-consuming and hence analysis of small data sets is beneficial to demonstrate the algorithm's performance. In the Blackfoot oil field, 13 well logs are present along with 3-D post-stack seismic reflection data, and the location of the well log in the seismic section is given in Table 1.
Location of wells in the seismic section.
| Well name | Inline | Crossline |
|---|---|---|
| 01–08 logs | 12 | 41 |
| 01–17 logs | 65 | 40 |
| 04–16 logs | 63 | 28 |
| 05–16 logs | 80 | 29 |
| 08–08 logs | 26 | 41 |
| 09–08 logs | 38 | 41 |
| 09–17 logs | 96 | 40 |
| 11–08 logs | 38 | 59 |
| 12–16 logs | 96 | 27 |
| 13–16 logs | 110 | 17 |
| 14–09 logs | 55 | 17 |
| 16–08 logs | 51 | 39 |
| Well name | Inline | Crossline |
|---|---|---|
| 01–08 logs | 12 | 41 |
| 01–17 logs | 65 | 40 |
| 04–16 logs | 63 | 28 |
| 05–16 logs | 80 | 29 |
| 08–08 logs | 26 | 41 |
| 09–08 logs | 38 | 41 |
| 09–17 logs | 96 | 40 |
| 11–08 logs | 38 | 59 |
| 12–16 logs | 96 | 27 |
| 13–16 logs | 110 | 17 |
| 14–09 logs | 55 | 17 |
| 16–08 logs | 51 | 39 |
Location of wells in the seismic section.
| Well name | Inline | Crossline |
|---|---|---|
| 01–08 logs | 12 | 41 |
| 01–17 logs | 65 | 40 |
| 04–16 logs | 63 | 28 |
| 05–16 logs | 80 | 29 |
| 08–08 logs | 26 | 41 |
| 09–08 logs | 38 | 41 |
| 09–17 logs | 96 | 40 |
| 11–08 logs | 38 | 59 |
| 12–16 logs | 96 | 27 |
| 13–16 logs | 110 | 17 |
| 14–09 logs | 55 | 17 |
| 16–08 logs | 51 | 39 |
| Well name | Inline | Crossline |
|---|---|---|
| 01–08 logs | 12 | 41 |
| 01–17 logs | 65 | 40 |
| 04–16 logs | 63 | 28 |
| 05–16 logs | 80 | 29 |
| 08–08 logs | 26 | 41 |
| 09–08 logs | 38 | 41 |
| 09–17 logs | 96 | 40 |
| 11–08 logs | 38 | 59 |
| 12–16 logs | 96 | 27 |
| 13–16 logs | 110 | 17 |
| 14–09 logs | 55 | 17 |
| 16–08 logs | 51 | 39 |
3 METHODOLOGY
In geophysical data inversion, local and global optimization techniques are frequently utilized to estimate subsurface rock properties. Each algorithm has its own set of benefits and drawbacks and this study presents an approach for merging the two algorithms so that one can overcome their flaws while maximizing the benefits of both. To address issues of geophysical importance, the study merged a local optimization technique based on QNMs with a global GA approach as a part of a hybrid optimization technique. Global optimization algorithms find the global minimum or maximum of an objective function and do not depend on starting model, whereas local optimization methods need a good starting model to get optimum results. In the following sections, the GA and HGA are described briefly.
3.1 Genetic algorithm
John Holland invented the GA in 1970 as a search-based optimization approach (Holland 1970). It is founded on the notion of ‘survival of the fittest,’ which is analogous to Darwin's evolutionary theory. The model parameters are first represented as a bit string in the GA. Individual fitness scores are derived using the fitness function from the space population generated from the sample. The population's two best individuals are chosen based on their fitness scores. Following the selection procedure, the selected individuals are subjected to crossover and mutations to improve their solutions (Chunduru et al. 1997; Maurya et al. 2020). This process is described briefly in the following subsections.
(i) Fitness function
Initially, a GA generates a random solution/population of 200 individuals to the problem and this is called the initial population. A fitness score based on fitness function (eq. 1) is used to assign to each solution. Based on these fitness scores, a selection procedure was applied.
(ii) Selection
After determining the fitness of each solution in the initial population, the first, genetic process, selection is adopted. Individual models are paired using the selection procedure based on their fitness scores. The model's fitness value is utilized to determine whether they are selected or not. When compared to models with low fitness values, models with a very high value have a far higher probability of being chosen. This is because the more fit models will be chosen for reproduction more frequently. The model's fitness value is directly connected to the selection probability (Goldberg & Holland 1988). The selection can be done in many ways, that is, roulette wheel selection, rank selection, tournament selection, etc., and this study uses the roulette wheel selection procedure. The roulette wheel selection is based on the probability estimated from their fitness score and the individuals with a larger probability have a greater chance to be selected. This selection procedure is more frequent in geophysical problem as provide fast convergence and less computational time. In selection; individuals with good fitness values based on the probability of selection are selected and passed to the next process, that is, crossover.
(iii) Crossover
After the models have been chosen and coupled, the crossover is the next genetic operator used in this process. The two selected individuals/models are called parent are paired and operations are performed in which a crossover scheme is used to exchange the bits of a binary string of parents to generate a new offspring. It transmits some information between paired individuals/models, resulting in the generation of new models (Boschetti et al. 1996; Velez-Langs 2005; Maurya et al. 2020). Single-point crossover and multipoint crossover are the two forms of crossover commonly employed in geophysical problems. The single-point crossover's main premise is that the model's one-bit position is chosen at random using a uniform probability distribution. Following that, all the bits on the right side of the selected bit are exchanged between two models, resulting in a new model. In the multipoint crossover procedure, on the other hand, a two-bit location (or more) is chosen at random and all bits between the chosen bit (every alternate bit segment in case of more than two crossover points are selected) are exchanged with each other. In this study, the multipoint crossover is adopted as the method is the best suitable technique for this study. After the crossover process, the next genetic operator, mutation is adopted.
(iv) Mutation
The mutation is the last step of the GA and it involves random perturbing in the solution. This is used to make randomness/diversity to the solution. In this step, a bit is chosen randomly and flipped with another bit which is again chosen as randomly in the same offspring. In most cases, the mutation process takes place during the crossover step. After the completion of these three processes, the new models are compared to the previous generation and accepted based on an updated probability (Chunduru et al. 1997; Aleardi & Ciabarri 2017; Maurya & Singh 2019; Maurya et al. 2019). The procedure is repeated until convergence is reached. A flowchart of the GA technique adopted in this study is presented in Fig. 3.
Flowchart of the adopted methodology of GA. Three processes of GA namely selection, crossover and mutation are also presented in the figure.
3.2 Hybrid optimization
A hybrid optimization algorithm uses both local as well as global optimization algorithms at the same time to optimize an objective function. Porsani et al. (1993) describe it as a coupled two-step procedure in which a local algorithm is applied to the best model that is obtained from global optimization. These global optimization algorithms do not require a good initial model and a good solution can be obtained even with a poor starting model. In this study, local optimization based on the QNMs is chosen over other local optimization methods along with the GA as a global optimization technique. These approaches are primarily used to obtain the global minimum of a twice-differentiable function. For vast and complex nonlinear problems, QNMs have distinct advantages over other local optimization. However, depending on the sort of QNM utilized and the situation in which it is applied, these methods are not ideal and can have significant limitations. Despite this, except for very basic problems, QNMs are often worthwhile to use in the field of geophysics.
To optimize twice-differentiable functions, a variety of QNMs are applied. The particular methods they utilize to approximate the inverse Hessian matrix (H-1) that is used in the full Newton's method are different between them. The main advantage of QNMs is that they do not require the iterative calculation of the inverse Hessian matrix. As a result, the different types of QNMs are greatly reliant on the approximation utilized. The most basic approximation utilizes the same inverse Hessian value for each iteration. This is not ideal, but it serves as an example (Chunduru et al. 1997).
Initially, a global optimization technique is applied to data and optimum results are estimated. Mathematically, the fitness value between observed and modelled data should be close to zero but practically its not possible due to the limitation of computing power, cost and time. Generally, 1000–5000 iterations are chosen for the optimization depending on computing power although the global optimization needs mathematically infinite iterations to converge completely. To reduce iterations, the hybrid optimization technique developed in this study is very useful. First, the GA technique is adopted and a close solution to the problem is expected but considering less number of iterations, the exact solution to the problem is not possible to achieve. To overcome this problem, local optimization based on QNM is adopted after GA. The benefit of using QNM is that the local optimization technique needs a good starting model to perform better. The solution found in GA is used as the initial point to the QNM-based optimization and hence converges exactly (Krahenbuhl & Li 2009; Zhao & Sen 2021). In this study, 1500 iterations are selected for global optimization (GA) and 28 iterations are chosen for local optimization techniques as stopping criteria. The detailed working principle of the hybrid optimization technique is shown through a flowchart in Fig. 4.
A flowchart shows the working mechanism of the HGA. To reduce the error between modelled and actual data, local optimization is applied to data after GA.
The biggest drawback of local optimization is to converge at local minima if the starting points are not close to the global optimum solution and the biggest advantage of it is to converge towards to solution very quickly. Global optimization needs mathematically infinite time to completely converge to the optimum solution and that is not physically possible hence hybrid optimization is chosen to get better results in a short time. Using local optimization after global optimization is to choose starting points close to the optimum solution so that local optimization does not fall to the local minima.
In this study, GA and HGA algorithms were first tested on the synthetically generated data to optimize parameters and thereafter, applied to the real seismic reflection data from the Blackfoot oil field, Canada. After, qualitative and quantitative comparisons of both algorithms, reservoir delineation from inverted results are performed. To confirm this reservoir, porosity volume is predicted using an ANN technique. The predicted results validate the reservoir as well as cross-verify the inverted impedance found from GA and HGA methods. The detailed discussion about the ANN is as follows.
3.3 Artificial neural network
ANN, systems often known as neural networks, are mathematical models of theorized mind and brain activity. They have shown to be effective at recognizing patterns and problems that can be solved via practice. ANN is neutrally inspired brain and behaviour models; that is, neurological models are based on how the brain processes information. A neuron or a processing element (PE) or a node, which is a multi-input, single-output element, is the basic building unit of the ANN (Hajian & Styles 2018).
Under an ANN, there are many algorithms available, but in this study, the PNN is adopted as a prediction tool because the methods work more accurately and are less time taking processes (McCormack 1991; Peters et al. 2019). The PNN is a mathematical interpolation scheme that uses neural network architecture for its implementation. The methods analyse a variety of attributes estimated directly or indirectly from the seismic data the and best attributes are selected for further analysis (Bosch et al. 2010). Thereafter, these selected attributes are cross-plotted with the desired well-log property (i.e. porosity in this case). A best nonlinear fit gives the relation between porosity and seismic attributes that are used for predicting porosity in the seismic section (Bhatt & Helle 2002). The entire process is carried out separately for all the sample points and finally, porosity volume is generated (Masters 1995; Mahmood et al. 2017).
4 RESULTS AND DISCUSSIONS
A post-stack seismic inversion using a GA as well as an HGA is carried out in this study. An HGA is a coupled two-step process in which the first global optimization algorithm, that is, GA is run to generate the best possible solution for an objective function followed by a local optimization algorithm, that is QNM is applied to the best model which is obtained from each generation of the GA. The algorithm is tested on synthetic and real data and inverted results along with well-log information are interpreted. Synthetic data sets are generated by convolving earth reflectivity obtained from well 01 to 08 to build the initial subsurface model. The inversion results are analysed in conjunction with well logs and anomalous impedance zone depicting as a glauconitic sand channel and its areal extent of the reservoir within the Blackfoot field is obtained as the final output.
In this study, seismic inversion has been performed in two steps. In the first step, the inversion of synthetic data is carried out and inverted results are compared with a true solution to optimize inversion parameters. Thereafter, in the second step, real seismic data are inverted and subsurface acoustic impedances are estimated which is interpreted for the anomalous zone. These steps are performed for global optimization as well as for hybrid optimization techniques. The detailed discussions are as follows.
4.1 Synthetic data: algorithm testing
Using eqs (5)–(7), a synthetic seismogram is generated from well 08–08 and plotted in Fig. 5. The first track in Fig. 5 shows well-log impedance, track 2 depicts calculated reflectivity (using eq. 7) and track 3 shows generated synthetic trace. Track 3 has three traces, the first trace corresponds to 0 per cent noisy data, the second trace corresponds to 10 per cent noisy data and the third trace contains 20 per cent white Gaussian noise to data. The percentage of noise means that the signal-to-noise ratio (SNR) is an input argument that uses dB units to represent the ratio of signal-to-noise powers. In this study, we read the 10 per cent white Gaussian noise power as signal power 1 and noise power 0.1, resulting in a 10 dB setting for the SNR input. Similarly, for the 20 per cent noise case, we have used 0.2 as noise power and 1 unit as signal power. The idea to add noise to the data is to test the algorithm for noisy data along with noise-free data. Further, tracks 4, 5 and 6 of Fig. 5 show inverted results which are compared with original curves from well-log data for 0 per cent noise, 10 per cent noise, and 20 per cent noisy data, respectively. The red and green curves (Fig. 5) show inverted impedance from GA and HGA, whereas the blue curve is the original impedance from well-log data. The lower bound, upper bound and initial model is shown by black lines on each track. From the figure, one can note that the inverted curves from GA as well as HGA are following very well with original impedance from well-log data even for noisy data also. This is also evidenced by the correlation coefficient as it is 0.99 for almost all cases. Other statistical parameters are shown in Table 2. These statistical parameters are depicting that the algorithm is performing well for synthetic data cases.
Comparison of inverted impedance and well-log impedance. Track 1 depicts well-log impedance, track 2 is reflectivity, track 3 is synthetic trace with 0 per cent, 10 per cent, and 20 per cent noise in the data, and tracks 4, 5 and 6 depict a comparison of inverted impedance with actual impedance for 0 per cent, 10 per cent and 20 per cent noise in the data.
Statistical parameters for synthetic data case.
| Parameters | Well logs (A) | Hybrid (B) | Difference (A −- īB) | GA (C) | Difference (A − C) |
|---|---|---|---|---|---|
| Min | 3206 | 3306 | 3.1 | 3031 | 5.5 |
| Max | 17 000 | 16 930 | 0.4 | 17 070 | 0.4 |
| Mean | 7835 | 7878 | 0.5 | 7904 | 0.9 |
| Median | 7770 | 7778 | 0.1 | 7863 | 1.2 |
| Mode | 3206 | 3364 | 4.9 | 3031 | 5.5 |
| Std. | 1794 | 1803 | 0.5 | 1811 | 0.9 |
| Range | 13 800 | 13 570 | 1.7 | 14 040 | 1.7 |
| Parameters | Well logs (A) | Hybrid (B) | Difference (A −- īB) | GA (C) | Difference (A − C) |
|---|---|---|---|---|---|
| Min | 3206 | 3306 | 3.1 | 3031 | 5.5 |
| Max | 17 000 | 16 930 | 0.4 | 17 070 | 0.4 |
| Mean | 7835 | 7878 | 0.5 | 7904 | 0.9 |
| Median | 7770 | 7778 | 0.1 | 7863 | 1.2 |
| Mode | 3206 | 3364 | 4.9 | 3031 | 5.5 |
| Std. | 1794 | 1803 | 0.5 | 1811 | 0.9 |
| Range | 13 800 | 13 570 | 1.7 | 14 040 | 1.7 |
Statistical parameters for synthetic data case.
| Parameters | Well logs (A) | Hybrid (B) | Difference (A −- īB) | GA (C) | Difference (A − C) |
|---|---|---|---|---|---|
| Min | 3206 | 3306 | 3.1 | 3031 | 5.5 |
| Max | 17 000 | 16 930 | 0.4 | 17 070 | 0.4 |
| Mean | 7835 | 7878 | 0.5 | 7904 | 0.9 |
| Median | 7770 | 7778 | 0.1 | 7863 | 1.2 |
| Mode | 3206 | 3364 | 4.9 | 3031 | 5.5 |
| Std. | 1794 | 1803 | 0.5 | 1811 | 0.9 |
| Range | 13 800 | 13 570 | 1.7 | 14 040 | 1.7 |
| Parameters | Well logs (A) | Hybrid (B) | Difference (A −- īB) | GA (C) | Difference (A − C) |
|---|---|---|---|---|---|
| Min | 3206 | 3306 | 3.1 | 3031 | 5.5 |
| Max | 17 000 | 16 930 | 0.4 | 17 070 | 0.4 |
| Mean | 7835 | 7878 | 0.5 | 7904 | 0.9 |
| Median | 7770 | 7778 | 0.1 | 7863 | 1.2 |
| Mode | 3206 | 3364 | 4.9 | 3031 | 5.5 |
| Std. | 1794 | 1803 | 0.5 | 1811 | 0.9 |
| Range | 13 800 | 13 570 | 1.7 | 14 040 | 1.7 |
Further, cross-plots are generated between inverted impedance and original impedance and presented in Fig. 6. Fig. 6(a) compares results for 0 per cent noisy data, Fig. 6(b) compares for 10 per cent noisy data, and Fig. 6(c) compares for 20 per cent noisy data. From the figure, one can note that, the scatter points fall very close to the best-fitted line for both cases, GA (black solid) and HGA (black circle) which again depicts that the inverted results are very close to the original impedance.
Cross-plot between inverted impedance and original impedance for (a) 0 per cent noise, (b) 10 per cent noise and (c) 20 per cent noisy data.
In this study, the aim was to reduce errors between observed and modelled traces, and this is performed in two ways here one using a GA and the second, using an HGA. The variation of error with iterations is presented in Fig. 7. Fig. 7(a) depicts error variation for the GA, whereas Fig. 7(b) depicts error for the HGA. From the figure, one can note that major error minimization has been done within 50 iterations and thereafter error minimization is very small. From this error plot, one can also note that the local optimization as part of an HGA does not reduce much error. This is because the synthetic data are perfectly modelled and hence GA has already reached the optimum solution and no space is remaining for a local optimization to reduce error.
Variation of error with iterations for (a) GA and (b) HGA.
The biggest challenge in inversion based on global optimization is convergence time as these methods take more time as compared to the local optimization techniques only. In this study, by using Intel (R) Core (TM) i7-11 700@2.50 GHz with 8GB RAM, the x64-based processor time taken to invert a single synthetic trace (shown in Fig. 5) is 1926, 1939, 1950, 1976, 1991 and 2004 s for GA with 0 per cent, 10 per cent, 20 per cent noisy data and HGA with 0, 10 and 20 per cent noisy data, respectively. This shows that hybrid optimization (HGA) is even more time-consuming as compared with global optimization (GA) methods. After getting satisfactory results from synthetic data, analysis is carried forward to the real seismic reflection data from Blackfoot field, Canada, and results are found as follows.
4.2 Real seismic data
After the successful application of GA and HGA to synthetic data, real seismic data from the Blackfoot field is used to invert the acoustic impedance, and the results are presented in Figs 8–13. Fig. 8 shows the seismic section at inline 41 and the study aims to invert this seismic section into subsurface acoustic impedance using GA and HGA algorithms discussed in this paper. The benefit of inversion is to transform seismic amplitude into subsurface acoustic impedance which shows layer property, whereas the former shows only interface property. This layer property can provide fluid and lithology information, whereas the seismic section can provide only interface information. Further, the acoustic impedance information is always present in the well log but the problem is that these well logs are very scattered and hence 2-D information cannot be extracted from the well-log only. Fig. 8 highlights different horizons along with high amplitude contrast which may be due to the presence of shale and sand interface.
Seismic reflection data from the Blackfoot field, Canada. Different horizons along with high amplitude anomalies are highlighted. This section shows 50 traces only along with 200 ms two-way traveltime (900–1100 ms) that will be inverted for acoustic impedance.
Low-frequency acoustic impedance model generated from interpolation of well-log impedance in seismic section. A high-cut filter of 10–15 to 50–80 Hz is applied to remove a higher frequency in the data.
Process of seismic to well tie. Track 1 is density log, track 2 depicts P- and S-wave velocities, track 3 is statistical wavelet, track 4 is synthetic trace plotted 5 times, track 5 shows composite trace plotted 5 times and track 6 depicts seismic reflection data near to well location.
(a) Cross-section of inverted impedance generated using GA and (b) HGA. Well-log impedance from well 08–08 is also plotted on this cross-section for comparison purposes only. A good match along with a low impedance zone is highlighted by an ellipse.
(a) Comparison of Blackfoot seismic reflection data (black) with reconstructed synthetic (red) from inverted impedance generated using GA. A correlation of 0.64 shows quite good match. (b) Depicts the difference between Blackfoot seismic and reconstructed synthetic. (c) Comparison of Blackfoot seismic reflection data with reconstructed synthetic from inverted impedance derived using HGA. In this case, a correlation of 0.79 is estimated shows clearly higher matching as compared with GA case. (d) Depicts the difference between Blackfoot seismic and reconstructed synthetic traces.
Variation of error with iterations. The black lines show error for the GA, whereas the black and red lines together show variation of error for the HGA.
For the inversion of seismic reflection data, the low-frequency acoustic impedance model plays a very important role which is generated by the interpolation of well-log acoustic impedance in the seismic section using the seismic horizon as a guide. To generate a low-frequency AI model, a high-cut filter of 10,15, 60 and 80 Hz is used to remove higher frequencies from the data to get a low-frequency model. This step aims to supply low-frequency information to the inverted results as the inversion of seismic reflection data contains a higher frequency of more than 10 Hz and hence lacks a lower frequency which is very important for deeper information. The low-frequency acoustic impedance section is presented in Fig. 9 which is generated from interpolation of well-log acoustic impedance and high-cut filter.
This study uses seismic as well as well-log data together as a part of inversion input and these two data sets are always in a different domain (time and depth), and hence they need to convert into a single domain (time) before their application. Hence, depth-to-time conversion is performed in this study, and related results are presented in Fig. 10. In this figure (Fig. 10), 6 tracks are created and the first track shows the density log with depth from the well 08–08, track 2 shows P- and S-wave velocities from the well 08–08, track 3 shows the wavelet which is extracted from seismic reflection data, track 4 shows synthetic trace generated from velocity and density plotted in track 1 and 2. This single synthetic trace is plotted 5 times in track 4 to look like seismic data. Further, track 5 of Fig. 10 shows composite traces extracted from seismic data near the well location (well 08–08) and again plotted 5 times to look-like data. In last, track 6 of Fig. 10 shows seismic reflection data near the well location and highlights a single seismic trace (red colour) which is extracted and plotted against synthetic data (in track 5) to compare. The idea is that the extracted seismic trace should be very close to the synthetic trace generated from the well log as both are coming from the approximately same location and the same lithology in both data is expected. In the initial stage, both synthetic and seismic trace is not matching with each other because synthetic data comes from a well log and this well log is in-depth whereas seismic data are in time. Stretching and squeezing have been applied to the synthetic data to match higher peaks from the synthetic trace to the seismic trace. These higher peaks are expected to come from the same lithology and the moment when these peaks match with each other can provide time and depth relation which can be used at any point of time during inversion to transform one domain into another.
Inversion is performed trace by trace and once completed inversion of the entire seismic trace using GA and HGA, results are plotted against two-way traveltime and shown in Figs 11(a) and (b), respectively. From Fig. 11, one can note very high-resolution subsurface acoustic impedance in both cases (GA and HGA) with slightly better results in the HGA. In both figures (Fig. 11), well-log acoustic impedances are also plotted at a well location to compare inversion results. From the figure (Fig. 11), one can note that the inverted impedance matches very nicely with well-log impedance and shows good inversion results in both cases (GA and HGA). The interpretation of the inverted section shows a low impedance zone (7500–9500 m s−1*g cc−1) near to 1055–1065 ms time interval in both sections. This low impedance zone is cross-validating with amplitude contrast (Fig. 8) and interpreted as a sand channel.
Further, to examine inverted impedance, a synthetic seismic gather is reconstructed from the inverted impedances and plotted beside the real seismic reflection data. Mathematically, these two data sets (seismic and reconstructed synthetic) should match each other as one can say that we tried to minimize the error between these two data sets as a part of inversion. The related results are presented in Fig. 12. Fig. 12(a) depicts reconstructed synthetic gather from the GA-based impedance section (using eqs 4–7) plotted in red colour along with the real seismic section from the Blackfoot field shown in black colour. From the figure (Fig. 12a), one can easily note that the reconstructed synthetic and Blackfoot seismic are matching quite nicely. The average correlation between these two data sets is calculated to be 0.64 which is relatively high and depicts good matching between reconstructed synthetic and Blackfoot seismic data sets. Further, the misfit between the Blackfoot seismic and reconstructed synthetic section is estimated and given in Fig. 12(b). From Fig. 12(b), one can note that the misfit energy is insignificant at most of the places and again demonstrated good matching between the reconstructed synthetic and Blackfoot seismic and hence good inversion results. The reconstructed synthetic is also generated for HGA-derived impedance and plotted in Figs 12(c) and (d). Fig. 12(c) shows a comparison of the reconstructed synthetic section (red colour) from HGA-derived impedance and Blackfoot seismic reflection data (black colour). From Fig. 12(c), it is noted that the reconstructed and Blackfoot seismic are matching very nicely with each other. The average correlation coefficient between the reconstructed synthetic section and the Blackfoot seismic section is estimated and found to be 0.79 which is a relatively high value and depicts a good match between the two data sets. Further, the misfit between the reconstructed synthetic and Blackfoot seismic is estimated and presented in Fig. 12(d). The visual examination of the misfit section depicts that there are no larger peaks present and hence indicates a good match between the reconstructed synthetic and Blackfoot seismic section. Further, if one compares the reconstructed synthetic from GA-derived impedance and the reconstructed synthetic from HGA-derived impedance, a better match of HGA with Blackfoot seismic and least misfit is found and demonstrates that the performance of HGA is better as compared with the GA algorithm.
The process used here for the inversion is that a synthetic seismogram is generated based on an initial guess model, and the error between this synthetic and Blackfoot seismic data is estimated. Thereafter, GA and HGA are used one by one to reduce this error by changing the initial model and at the moment when this error becomes small enough, the modified model is our final solution. The GA and HGA minimize error via iterations and the variation of error with iterations is presented in Fig. 13. In this study, it is decided to use 1500 iterations for the global optimization technique and 28 iterations for the local optimization technique. The less number of iterations is chosen for local optimization because of its very fast convergence of it. The HGA uses global optimization (GA) along with local optimization and this local optimization depends on starting points. And if the starting point is very close to the global optimum solution (which is provided by GA), then this local optimization reduces error very rapidly and converges to the global optimum solution. Further, each trace has its error variation with iteration and a small alteration can be expected among the error plot of different traces but overall errors are decreasing with iteration exponentially. From the figure (Fig. 13), one can note that the error decreases rapidly within 0–25 iterations, and thereafter gradual decrease occurs. If more iterations is chosen, more error minimization can be achieved but considering large data sets and convergence time, an optimum number of iterations, that is, 1500 in our case is chosen. In Fig. 13, the black lines are plotted for error minimization estimated for the GA however, local optimization error variations are presented by red lines. One can note that when applying local optimization after global optimization, the convergence is very fast however if one uses only the local optimization technique for inversion, this convergence cannot be so rapid. This is due to the selection of starting model for the local optimization as it depends largely on it. The application of local optimization just after global optimization provides starting points in the vicinity of the optimum solution and hence local optimization performs minimization very rapidly and gives the true solution.
Further, the time consumption by GA and HGA is very important to consider as mathematically the global optimization needs infinite time to completely converge to the optimum solution and physically it is not possible. Hence the stopping criterion depends on time consumption by these two methods (GA and HGA). The comparison of time taken by GA and HGA is estimated and found that 132 and 128 s convergence time is taken by HGA and GA, respectively to invert one composite trace. Further, 6540 and 6829 s convergence time is taken by the GA and HGA to invert 50 traces presented in Fig. 8. One can note that the convergence time of HGA for 50 traces is 5400 s more as compared with GA techniques. All these inversions are performed on a machine with the specification of Intel (R) Core (TM) i7-11 700@2.50 GHz with 8GB RAM, x64-based processor. The Blackfoot seismic data is recorded from 300 to 1300 ms two-way traveltime and the inversion is carried out only from 900 to 1100 ms (only 200 ms) time interval because the sand channel is present in this interval only. This small time interval is chosen to reduce the further time of convergence. In some cases, seismic reflation data are recorded more than 1 min two-way traveltime and hence one can imagine how much time going to takes these methods to get good results. Hence, one needs to understand convergence time before choosing data segments to be inverted, and methods to apply for inversion. Table 3 shows a quantitative comparison of convergence time for synthetic as well as real data cases. From Table 3, it has been noted that 4.4 per cent more convergence time is required for HGA as compared with GA to invert just 50 traces. This time difference can be even more for larger data sets.
Quantitative comparison of convergence time for synthetic as well as seismic data.
| Data | GA (in s) | HGA (in s) | Difference (in per cent) |
|---|---|---|---|
| Synthetic with 0 per cent noise | 1926 | 1976 | 2.60 |
| Synthetic with 10 per cent noise | 1939 | 1991 | 2.68 |
| Synthetic with 20 per cent noise | 1950 | 2004 | 2.77 |
| Composite trace | 128 | 132 | 3.13 |
| All seismic trace | 6540 | 6829 | 4.42 |
| Data | GA (in s) | HGA (in s) | Difference (in per cent) |
|---|---|---|---|
| Synthetic with 0 per cent noise | 1926 | 1976 | 2.60 |
| Synthetic with 10 per cent noise | 1939 | 1991 | 2.68 |
| Synthetic with 20 per cent noise | 1950 | 2004 | 2.77 |
| Composite trace | 128 | 132 | 3.13 |
| All seismic trace | 6540 | 6829 | 4.42 |
Quantitative comparison of convergence time for synthetic as well as seismic data.
| Data | GA (in s) | HGA (in s) | Difference (in per cent) |
|---|---|---|---|
| Synthetic with 0 per cent noise | 1926 | 1976 | 2.60 |
| Synthetic with 10 per cent noise | 1939 | 1991 | 2.68 |
| Synthetic with 20 per cent noise | 1950 | 2004 | 2.77 |
| Composite trace | 128 | 132 | 3.13 |
| All seismic trace | 6540 | 6829 | 4.42 |
| Data | GA (in s) | HGA (in s) | Difference (in per cent) |
|---|---|---|---|
| Synthetic with 0 per cent noise | 1926 | 1976 | 2.60 |
| Synthetic with 10 per cent noise | 1939 | 1991 | 2.68 |
| Synthetic with 20 per cent noise | 1950 | 2004 | 2.77 |
| Composite trace | 128 | 132 | 3.13 |
| All seismic trace | 6540 | 6829 | 4.42 |
In this study, an ANN technique is also adopted to cross-verify inverted results in the inter-well region. This ANN is a well-established method and can be implemented using Hampson Russell software suits to cross-validate inversion results. The results found using ANN are discussed in the following section.
4.3 Artificial neural network
An ANN is a tool used in the field of geoscience to predict various petrophysical parameters in inter-well regions. The results of the ANNs can be used to cross-verify the inverted impedance section as the porosity always compliment acoustic impedance and hence can be used to cross-verify results found from seismic inversion. Further, this predicted porosity section can also strengthen the interpretation of seismic reflection data. There are several types of neural networks available but in this study PNN is chosen as the method is more accurate and takes less computational time and expertise to implement. In PNN, the number of seismic attributes is estimated from seismic reflection data and the best number of attributes is selected based on their correlation with the porosity log. In this study, more than 20 seismic attributes like a seismic envelope, cosine instantaneous amplitude, cosine instantaneous phase, etc. are estimated. These attributes are nothing but the combination of the real and imaginary parts of complex traces
List of attributes used for the PNN.
| Target | Attribute | Training error | Validation error |
|---|---|---|---|
| Porosity | 1/HGA_impedance | 5.33 | 5.56 |
| Porosity | Filter 25/30–35/40 | 5.20 | 5.38 |
| Porosity | Amplitude envelope | 5.13 | 5.33 |
| Porosity | Instantaneous phase | 5.05 | 5.23 |
| Porosity | Amplitude envelope | 4.88 | 5.20 |
| Porosity | X-coordinate | 4.78 | 5.18 |
| Porosity | Integrated absolute amplitude | 4.72 | 5.12 |
| Target | Attribute | Training error | Validation error |
|---|---|---|---|
| Porosity | 1/HGA_impedance | 5.33 | 5.56 |
| Porosity | Filter 25/30–35/40 | 5.20 | 5.38 |
| Porosity | Amplitude envelope | 5.13 | 5.33 |
| Porosity | Instantaneous phase | 5.05 | 5.23 |
| Porosity | Amplitude envelope | 4.88 | 5.20 |
| Porosity | X-coordinate | 4.78 | 5.18 |
| Porosity | Integrated absolute amplitude | 4.72 | 5.12 |
List of attributes used for the PNN.
| Target | Attribute | Training error | Validation error |
|---|---|---|---|
| Porosity | 1/HGA_impedance | 5.33 | 5.56 |
| Porosity | Filter 25/30–35/40 | 5.20 | 5.38 |
| Porosity | Amplitude envelope | 5.13 | 5.33 |
| Porosity | Instantaneous phase | 5.05 | 5.23 |
| Porosity | Amplitude envelope | 4.88 | 5.20 |
| Porosity | X-coordinate | 4.78 | 5.18 |
| Porosity | Integrated absolute amplitude | 4.72 | 5.12 |
| Target | Attribute | Training error | Validation error |
|---|---|---|---|
| Porosity | 1/HGA_impedance | 5.33 | 5.56 |
| Porosity | Filter 25/30–35/40 | 5.20 | 5.38 |
| Porosity | Amplitude envelope | 5.13 | 5.33 |
| Porosity | Instantaneous phase | 5.05 | 5.23 |
| Porosity | Amplitude envelope | 4.88 | 5.20 |
| Porosity | X-coordinate | 4.78 | 5.18 |
| Porosity | Integrated absolute amplitude | 4.72 | 5.12 |
In this study, total of six attributes are selected for the final analysis that is directly calculated from seismic data, and one attribute is taken, that is, inverted impedance derived from HGA methods. Further, a cross-plot between these attributes with well-log porosity has been generated and the best-fitting polynomial provides a nonlinear relationship between selected attributes and desired well petrophysical parameters (porosity in this case). Finally, this relationship is used to predict porosity at each sample point of seismic reflection data. This analysis is performed in two steps, first, a single composite trace is extracted from seismic reflection data and porosity prediction is performed using equations derived by PNN. Second, porosity volume is generated by using the entire seismic section as input to the PNN method. The reason to predict porosity for the composite trace (in the first step) is that the well log contains a porosity log and the predicted log can be compared with actual porosity from the well log. Fig. 14 shows a comparison of actual porosity with the predicted porosity estimated from composite trace near well 08–08. In Fig. 14(a), the solid line shows predicted (modelled) porosity, whereas the dotted line shows well-log (actual) porosity. From the figure (Fig. 14a), one can note that the modelled curve is matching nicely with the actual curve although the actual porosity has a larger frequency range present. Further, a cross-plot between actual and modelled porosity is generated and plotted in Fig. 14(b). This figure (Fig. 14b) also shows good prediction results as most of the scatter points are close to the best-fitting line.
Comparison of modelled porosity with well-log porosity (a) against depth and (b) depicts a cross-plot between actual and predicted porosity. The best-fitting line (solid line) is also plotted on this scatter plot to compare predicted results.
Thereafter, porosity volume is generated by predicting porosity in the entire seismic section and one cross-section at inline 41 is presented in Fig. 15. Fig. 15 shows the variation of porosity with depth in the y-direction and porosity with offset in the x-direction. The green colour (in Fig. 15) shows a lower porosity value, whereas the pink colour is used for a higher porosity value. From Fig. 15, it is noted that the porosity varies from 10 to 20 per cent in the study region with lower porosity in most of the parts and higher porosity concentrated near 1055–1065 ms two-way traveltime. This high porosity zone is nicely correlated to the low acoustic impedance zone interpreted in inverted results found from GA (Fig. 11a) and HGA (Fig. 11b) techniques. The porosity cross-correlation with acoustic impedance validates inversion results and hence good algorithm performance. This high porosity zone and low acoustic impedance zone found between 1055 and 1065 ms time intervals are interpreted as a sand channel.
A cross-section of predicted porosity with high porosity zone (pink colour) is highlighted by an ellipse. This high porosity varies from 17 to 19 per cent.
5 CONCLUSIONS
To delineate the Glauconitic sand channel located in Manville formation of a lower cretaceous geological period of Blackfoot field Alberta, Canada, seismic inversion using a GA as well as an HGA is performed on post-stack seismic reflection data. The presence of well logs 08–08 at trace location 42 earth models are used as prior information in seismic inversion. In this study, the methodology for the inversion of seismic reflection data from GA and HGA is developed and applied to synthetic as well as real data. In an HGA, a global optimization technique, that is, GA along with a local optimization, that is, based on QNMs is used. The idea is that local optimization largely depends on the starting model, whereas global optimization does not depend on starting model. So the benefit of using both (local as well as global optimization) methods is that the final results found from GA can be used as input or starting models for the local optimization methods and hence one can reduce error very rapidly and can reach the optimum results. After the development of the methodology, the algorithm is tested with synthetically generated data, and thereafter, the inversion of the composite trace along with the entire seismic reflection data is inverted for acoustic impedance. The analysis shows that the HGA performs superior as compared to the GA alone although the HGA take 4.4 per cent more time for their convergence. The interpretation of inverted sections depicts the presence of a sand channel at 1055–1065 ms two-way traveltime. In this interval, low acoustic impedance ranging from 7500 to 9500 m s−1*g cc−1 is estimated. To cross-validate inversion results, porosity prediction is performed using an ANN technique. The results from ANN is correlating very well with the inverted acoustic impedance found from GA and HGA. The low acoustic impedance (7500–9500 m s−1*g cc−1) is well correlated with the high porosity (17–20 per cent) zone at 1055–1065 ms two-way traveltime and is interpreted as a sand channel. The study shows that both GA and HGA a much time taking process but HGA is providing more accurate and high-resolution subsurface information as compared with GA only. The study suggests that, before applying these methods to any data, their convergence time plays a very important role as most of the effort and cost is taken by the algorithm to converge towards the best solution.
ACKNOWLEDGEMENTS
The authors acknowledge the CGG Geo software for providing Hampson Russell software and data. One of the authors gratefully acknowledges the funding agency, UGC-BSR (M-14–0585) and IoE BHU (Dev. scheme no. 6031B) for their financial support.
DATA AVAILABILITY
The seismic and well-log data used in this study are provided by CGG Geo Software Pvt. Ltd This data can be accessed by contacting https://www.geosoftware.tech/.
CONFLICT OF INTEREST
The authors declare that they have no conflict of interests.
