Estimating reservoir zone from seismic reflection data using maximum-likelihood sparse spike inversion technique: a case study from the Blackfoot field (Alberta, Canada)

Seismic inversion involves extracting qualitative as well as quantitative information from seismic reflection data that can be analyzed to enhance geological and geophysical interpretation which is more subtle in a traditional seismic data interpretation. Among many approaches that have been made to improve interpretation of post-stack seismic data, a great effort has been made to use maximum likelihood (ML), sparse spike inversion (SSI) along with multi-attribute analysis (MAA) aimed to increase the resolution power of interpreting seismic reflection data and mapping into the subsurface lithology. These methods are applied to the Blackfoot seismic reflection data to estimate reservoir. The methods were first applied to the composite trace close to well locations and were inverted for acoustic impedance (AI). The results depict that the inverted AI matches very well with the well log AI. The statistical analysis demonstrates good performance of the algorithm. Thereafter, the entire seismic section was inverted to acoustic impedance section. The analysis of the inverted impedance section shows an anomaly zone in between 1060 and 1075 ms time and characterize it as reservoir. Further, the multi-attribute analysis is performed to estimate porosity and density in the inter-well region. The inverted porosity section shows a high porosity anomaly and a low density anomaly in between 1060 and 1075 ms time intervals which corroborated well with the low impedance zone and confirm the presence of a reservoir.


Introduction
Seismic inversion is a procedure that helps extract the high-resolution subsurface model of the physical characteristics of rocks and fluids from low-resolution seismic reflection data with the integration of well log data (Krebs et al. 2009). In the petroleum industry, the seismic inversion techniques have been widely used as a tool to locate hydrocarbon-bearing strata in the subsurface from the seismic data (Bosch et al. 2010;Maurya and Singh 2015). There are many approaches available for seismic inversion such as pre-stack simultaneous inversion, elastic impedance inversion, recursive inversion, model-based inversion, sparse spike inversion which includes linear programming and maximum-likelihood inversion, colored inversion, and band-limited inversion. (Russell and Hampson 1991). In the present study, maximum-likelihood sparse spike inversion (MLSSI) technique is used to estimate rock property in the subsurface using seismic reflection data from the Blackfoot field Alberta, Canada. The motivation to use MLSSI technique among others is that the technique is robust and fast to compute the subsurface model (Russell and Hampson 1991;Veeken and Silva 2004;Jackiewicz 2009;Maurya and Sarkar 2016). Other methods like colored inversion are also fast techniques, but this gives an average variation of the rock parameters, whereas MLSSI techniques give quantitative variation which corroborate well with the well log parameters (Russell and Hampson 1991;Maurya and Sarkar 2016).
Although the maximum-likelihood sparse spike inversion is very popular nowadays, this technique is not used by the Blackfoot datasets to characterize the reservoir. Maurya and Singh (2015) used maximum-likelihood sparse spike inversion along with linear programming sparse spike inversion technique to compare them. The present work is an extension of Maurya and Singh's (2015) paper.
The maximum-likelihood sparse spike inversion technique used the theory of maximum-likelihood deconvolution (MLD) which aims to generate reflectivity series from the seismic data (Russell 1988;Russell and Hampson 1991;Maurya and Sarkar 2016;Maurya et al. 2018a, b). The maximum-likelihood deconvolution algorithm came into the existence in 1982 and published by Kormylo and Mendel (1983) and Chi et al. (1984). The application of the MLD algorithm developed by Kormylo and Mendel to the real data is very complex and hence a modification was given by Hampson and Russell in 1985 which was easy to use. Hampson and Russell (1985) described that the MLD could be used to estimate reflectivity series from the broadband seismic reflection data and then the reflectivity can be transformed to the acoustic impedance which is more meaningful.
In the present study, the theory developed by Hampson and Russell (1985) is used to estimate rock property from the seismic reflection data. The algorithm is applied to the seismic section to estimate the variation of impedance and characterize the reservoir. Thereafter, for the confirmation of the reservoir, porosity and density volumes are predicted in the inter-well region using multi-attribute analysis technique. The theory of multi-attribute analysis was developed by Hampson et al. (2001). In this methodology, the aim is to find an operator, that is linear or nonlinear, that can predict well logs from neighboring seismic data. These methods choose to analyze not only low-resolution seismic reflection data but also the seismic attributes. As these attributes are generally nonlinear, it increases the prediction power of the multi-attribute analysis tool. Using this technique, porosity and density volumes are predicted from the seismic reflection data and used to characterize the prospective (reservoir) zone. The analysis is performed in Hampson Russell software and the results are presented using Matlab programming language.

The Blackfoot study area
The Blackfoot seismic survey was conducted to record data in two overlying patches. The first patch is targeted to the clastic Glauconitic channel of lower Cretaceous age which represents sediment-filled incised valley, while the second patch is targeted to the reef prone Beaver hill lake carbonates (Lawton et al. 1996). In this study, the seismic data related to the first patch, i.e., Glauconitic channel is used for the analysis. The reservoir in the study area is produced from the Glauconitic sand formation at a depth of 1550 m (~ 1060 ms), where Glauconitic sandstones and shales are incised into the regional Mannville stratigraphy valleys (Maurya and Singh 2018a, b).
The Blackfoot field is located about 70 km southeast from the Strathmore city, Alberta, Canada (Dufour et al. 1998). The 3C-3D (three dimensional-three component recording) data are recorded by Consortium for Research in Elastic Wave Exploration Seismology (CREWES) in 1995. The initial aim for this survey was to evaluate the effectiveness of integrated PP and PS survey for improved hydrocarbon survey (Lawton et al. 1996). The first patch is about 2.7 × 3.8 km in dimension. The data is recorded using 708 sources and 690 receivers with 140-fold at the center of spread. The other recording parameters are: source and receiver interval is 60 m, source line interval is 210 m, receiver line interval is 255 m, the number of source lines is 12 and the total number of receiver lines is 15. The bandwidth of the seismic reflection data is 12-90 Hz (Lawton et al. 1996;Margrave et al. 1998;Wood and John 1992).
The primary target of the survey was to see the amplitude changes caused by the Glauconitic Member of the Mannville Group. The Glauconitic member consists of shales and sands of lacustrine (Lacustrine deposits) and channel origin (river stream has carried sediment into the basin). Within the channel, the sediments are subdivided into three units corresponding to three phases of valley incision with different qualities of sand deposited. These three units may not be encountered everywhere in the area (Miller et al. 1995). The detailed stratigraphic column of the Blackfoot area is discussed in Margrave et al. (1998) (Maurya and Singh 2018b). Figure 1a shows the location of first and second patch and Fig. 1b depicts the stratigraphic division of the Blackfoot region, Alberta, Canada. Post-stack seismic reflection data and thirteen wells from the Glauconitic region are used for the characterization of the reservoir. Figure 2 shows the location of wells in the area.

Maximum-likelihood inversion
Maximum-likelihood sparse spike inversion is used to estimate reservoir property from seismic reflection data of the Blackfoot field, Alberta, Canada. The workflows used in this study are as follows: first, a couple of seismic horizons are picked near to 930 and 1040 ms two-way travel time. These horizons serve as guides to interpolate well log property in inter-well regions (Kormylo and Mendel 2014;Zhang and Castagna 2011). Thereafter, wavelet is extracted from the seismic as well as well log data. The steps are as follows: (i) extract the analysis window in the seismic data/well log data, (ii) taper the start and end of the window with a taper length, (iii) calculate the autocorrelation of the data window, (iv) calculate the amplitude spectrum of the autocorrelation, (v) take the square root of the autocorrelation spectrum. This approximates the amplitude spectrum of the wavelet, (vi) add manually the desired phase, (vii) take the inverse FFT (fast Fourier transform) to produce the wavelet, (viii) sum this with other wavelets calculated from other traces in the analysis window (Maurya and Singh 2018a, b). This desired wavelet is used to generate synthetic seismograms from the well log data which is used for well correlation. Further, time to depth conversion is performed. These can be done in one of two ways, first, using check shot data (if available) and second, manually. In the present study, check shot data is not available, hence the time to depth conversion is performed manually. This manual conversion is performed for all the available wells. For this, one particular well is selected and synthetic seismograms are generated by the convolution of extracted wavelet and Earth's reflectivity series from the well log data (Krebs et al. 2009). Thereafter, these synthetic traces are plotted along with the seismic traces near to that well and the matching of higher amplitudes is noticed. The seismic amplitudes are stretched and squeezed to match them properly. The condition when all the major amplitudes match, gives the relation between time and depth which is used to convert well log data into time domain (Helgesen et al. 2000). This is repeated for all the wells. The process of time to depth conversion is important, since it improves the correlation of events on the well log data and the events on the seismic data. Thereafter, initial impedance model of low frequency is generated. This low-frequency model is important, because the seismic data have limited bandwidth and does not contain lower frequency band, hence to get broadband spectrum of the inverted results, this low-frequency model is added to the inverted acoustic impedance results (Li and Speed 2004;Banihasan et al. 2006;Maurya and Singh 2015). After pre-processing, seismic inversion can be performed to estimate subsurface rock properties such as impedance, density, and porosity velocity.
Inversion is performed in two steps. In the first step, maximum-likelihood deconvolution is employed to estimated reflectivity series from the seismic reflection data. The maximum-likelihood deconvolution assumes that the Earth's reflectivity series are a combination of large spikes superimposed with smaller Gaussian events in the background (Russell and Hampson 1991;Mendel 2012). These large spikes are geologically indicating the depositional gap and hence they are important.  The objective function can be written for reflectivity and noise which may be minimized to produce the reflectivity and wavelet combination consistent with the statistical assumption (Russell 1988;Li and Speed 2004). The objective function E is given as where r k is reflection coefficient at the kth interface, m is the total number of reflections, L is the total number of samples presents, N is the square root of noise variance, n k is the noise at kth sample/interface, is likelihood that a given sample has reflection and R is the square root of reflectivity variance.
In the second step, the reflectivity series are transformed into the acoustic impedance which is more meaningful compared to the reflectivity series. This can be achieved using two approaches, the first approach is that the acoustic impedances can be estimated using the direct relationship between reflectivity and impedance (Velis 2006(Velis , 2007Maurya and Sarkar 2016). However, this relationship has a limitation, i.e., if the data contains significant noise, then the transformation cannot be performed properly as the relationship does not include the noise component. A second approach is used in this study to generate acoustic impedance model from the reflectivity section which included the noise component and hence gives complete spectrum (Russell 1988;Wang et al. 2006;Zhang et al. 2016). The equation is as follows: where (i) is the known impedance trend, r(i) is Earth's reflectivity series, n(i) is the noise in the input trend, i is the location of sample points and H(i) is a factor that would be The inversion is first tested to the one composite trace close to the well location and cross verified by the well log parameters. After getting satisfactory results, the entire seismic volume is inverted for acoustic impedance. Thereafter, porosity and density are estimated using multi-attribute analysis technique which analyzes a number of attributes generated directly/indirectly from the seismic data.

Multi-attribute analysis
There are numerous approaches available to estimate porosity in the inter-well region such as Cokriging (Doyen 1988;Doyen et al. 1996), single attribute analysis (Maurya et al. 2018a, b), multi-attribute analysis (Maurya et al. 2018a, b), probabilistic neural network (Maurya et al. 2018a, b), relationship derived from the well logs (Maurya and Singh 2018a, b), multilayer feed forward neural network (Leite and Vidal 2011), and fuzzy logic (Kadkhodaie-Ilkhchi et al. 2009), but the multi-attribute analysis is very simple, easy to use and has a fast computation process (Doyen 1988). Hence, we selected it in this study to estimate the porosity and the density. The multi-attribute analysis has more than one attribute at a time to estimate the relationship. The multi-attribute analysis uses internal attributes which are estimated directly from the seismic data and external attribute which is derived from the seismic inversion technique (Haas and Dubrule 1994;Eskandari et al. 2004). The best attributes are selected on the basis of their correlation with the target logs. Thereafter, the selected attributes are cross-plotted corresponding to the well log property at a particular time sample to estimate the relationship between attributes and target well log property which is used to estimate porosity and density at that time sample (Russell et al. 1997;Hampson et al. 2001). The same process has been repeated at each time sample and porosity and density is estimated (Pramanik et al. 2004). Figure 3 shows the flow chart of the maximum-likelihood inversion and multiattribute analysis technique. This can be shown mathematically as follows.
Let three attributes A 1 , A 2 and A 3 be the best for the prediction of porosity which shows low error and high correlation coefficient with well log porosity. A cross-plot has been generated between these attributes and well log porosity. Further, a straight line may fit the regression which can be represented as where w 0 , w 1 , w 2 and w 3 are weights and is desired porosity. These weights can be estimated by minimizing the mean (4) (t) = w 0 + w 1 A 1 (t) + w 2 A 2 (t) + w 3 A 3 (t), Fig. 3 Flowchart of seismic inversion and multi-attribute analysis methods squared prediction error between attributes and target log, i.e., porosity (Hampson et al. 2001;Leiphart and Hart 2001) The solution of four weights can be written (Hampson et al. 2001) as: where i is the location of sample point. The mean squared error (Eq. 5) calculated using the derived weights (Eq. 6) constitutes goodness-of-fit measure for the transform (Hampson et al. 2001). A similar expression can be derived for prediction of density in inter-well region. (5)

Seismic inversion
In this study, seismic inversion is performed to estimate rock property from seismic reflection data with integration of well log data to characterize reservoir from the Blackfoot field, Canada and the results are shown in Figs. 4,5,6,7 and 8. Figure 4 shows the inversion results for the composite trace near to the well locations. In Fig. 4, the black solid curve depicts well log acoustic impedance, whereas the inverted acoustic impedances are shown in blue solid color. The analysis is performed for all the composite traces near to well locations, but the results are shown for four wells only  for simplicity. From Fig. 4, it is noticed that the inverted impedances are agreeing well (average correlation > 0.93) with the well log impedances. The well log data are not recorded from surface to bottom in all wells and hence in some wells, the data are only available between 1000 and 1100 ms time intervals, e.g., Fig. 4c. Thereafter, a cross-plot of the inverted impedance with actual impedance from all the wells are generated and shown in Fig. 5. The distribution of scatter points shows the quality check of inversion results and hence the performance of the algorithm. It is noticed that the inverted data points are very close ( r 2 = 0.59 ) to the best fit curve and hence depicts a good performance of the algorithm. However, for quantitative comparison, the error between inverted and actual impedance is estimated and shown graphically in Fig. 6. Figure 6a shows variation of correlation coefficient, whereas Fig. 6b depicts synthetic relative error (SRE) and Fig. 6c shows variation of root-mean-square (RMS) error with wells. It is noticed that the correlation coefficients vary from 0.91 to 0.94 with a mean value of 0.93, SRE varies from 0.19 to 0.23 with a mean value of 0.21 and RMS error varies from 1000 to 2100 m/s*g/cc with a mean value of 1300 m/s*g/cc. These values indicate that the algorithm successfully searched the optimum solution.
Thereafter, the algorithm is applied to the post-stack seismic data from the Blackfoot field, Alberta, Canada with the aim of characterizing reservoir. The seismic section is inverted for acoustic impedance and shown in Fig. 7 (inline 118). The special features are highlighted in the inverted impedance section. The analysis of inverted impedance section shows a low impedance anomaly ranging 6000-9000 m/s*g/cc within 1060-1075 ms time interval which may be due to the presence of reservoir.
Further, the amplitude spectrum of inverted synthetic and Blackfoot seismic is compared and shown in Fig. 8 to see the frequency match between them. The figure depicts that for a lesser frequency, the amplitude is matching accurately with the seismic amplitude, whereas for a larger frequency (> 80 Hz), the amplitude spectrum shows a small deviation from the seismic amplitude, however, the overall correlation is good (~ 96%) between them. Thus, it can be concluded that the MLSSI algorithm distorts the amplitude spectrum particularly for larger frequency during its implementation, however, the amount of distortion is small.

Multi-attribute analysis
For the confirmation of reservoir, the other important petrophysical parameters, i.e., porosity and density is predicted in the inter-well region using multi-attribute analysis (MAA) technique and using MLSSI-derived impedance as an external attribute and seismic-derived attributes as internal and important results are shown in Figs. 9, 10, 11 and 12. Figure 9 shows the result of multi-attribute analysis to predict porosity and density from composite trace near to well locations. Figure 9a, b depicts comparison of predicted and actual porosity values for traces near to well 01-17 and well 08-08, respectively, whereas Fig. 9c, d depicts comparison of predicted density with actual density for traces near to well 01-17 and 08-08, respectively. From Fig. 9, it is noticed that the predicted porosity and density agree well with the actual porosity and density, respectively. The Fig. 9 Comparison of predicted porosity (blue) and actual porosity (black) for the a well 01-17, b well 08-08, and comparison of predicted density (green) with actual density (red) for c well 01-17, and d well 08-08 correlation coefficients are estimated to be 0.82 and 0.86 for porosity and density case, respectively, which indicate good values and hence a good performance of the algorithm.
A cross-plot of the predicted values with the actual values from the well log is generated to see the quality check of the predicted results and is shown in Fig. 10. Figure 10a shows the cross-plot of predicted porosity versus actual porosity, whereas Fig. 10b shows cross-plot of predicted density versus actual density for all the 13 wells available in the study area. It is noticed that the distribution of scatter points is very close to the best fit line for both cases and depicts that the inverted results are very close to the actual values and hence it point towards a good performance of the algorithm. This can be also verified by the r 2 values which is estimated for both the cases (porosity and density) and found to be 0.56 and 0.58, respectively.
The validation error is estimated for porosity and density case and shown in Fig. 11. It is noticed that the validation error varies from 1.4 to 6% for both the cases which indicate good prediction results.
After getting satisfactory results from the composite trace, the entire seismic section is inverted first for porosity and then for density sections. Figure 12a depicts predicted porosity section, whereas Fig. 12b depicts predicted density section at inline 118 only for simplicity. Although the entire seismic volumes are inverted for impedance, porosity and density, the inline 118 is shown here for simplicity. It is noticed that there is a high porosity anomaly ranging 15-20% and a low density anomaly ranging 2.0-2.3 g/ cc between 1060 and 1075 ms time intervals are present which corroborated well with the low impedance anomaly and confirm the presence of reservoir. The reservoir zone is highlighted by the ellipse in both (porosity and density) the sections.
Further, a cross-plot among inverted impedance, porosity and density volumes is generated and shown in Fig. 13. This cross-plot helps to derive a relationship among the inverted parameters which can be used for quick estimation of these parameters in the Blackfoot region. Figure 13a depicts crossplot of predicted porosity versus predicted density with color bar as inverted impedances, whereas Fig. 13b depicts crossplot of inverted impedance versus predicted density with color bar as predicted porosity and Fig. 13c depicts crossplot of inverted impedance versus predicted porosity with color bar as predicted density. The reservoir zone is characterized by low acoustic impedance (6000-9000 m/s*g/cc), low density value (2-2.35 g/cc), and high porosity (15-20%) and highlighted by ellipse in Fig. 13. The relations among where ρ is the density, AI is the acoustic impedance and ϕ is the porosity. These relations can be used to estimate of density, porosity and/or acoustic impedances if any of these quantities are known. The highlighted zone in the cross-plots are mapped into the seismic sections and shown in Fig. 14. Figure 14a highlights the low impedance zone, Fig. 14b highlights the high porosity zone and Fig. 14c highlights the low density zone in between 1060 and 1075 ms time interval in the seismic section at inline 118. These zones correspond to the reservoir. The predicted reservoir is also noticeable from the well log data at same location and strengthens the prediction results. The location of same reservoir is also interpreted by Maurya and Sarkar (2016), . The reservoir varies from NE to SW direction in the subsurface.

Conclusions
The maximum-likelihood sparse spike inversion (MLSSI) technique is utilized in the present study to estimate the reservoir zone from the seismic reflection data from the Blackfoot field, Alberta, Canada. The aim of this study is to transfer invisible information from the band-limited seismic Fig. 13 Cross-plot of a predicted porosity versus predicted density, b predicted density versus inverted impedance, and c predicted porosity versus inverted impedance seismic section centered 1065 ms with 100 ms width. The reservoir zone is highlighted by the ellipse. The relationships between inverted properties are shown on top of the respective figures reflection data to the blocky structure of the subsurface and hence characterize the reservoir. Thereafter, the density and porosity volumes are estimated in inter-well region using multi-attribute analysis (MAA) technique for the confirmation of reservoir zone. In the first step, MLSSI technique is applied to the composite trace near to well locations and inverted for acoustic impedances. The results depict that the inverted impedance is agreeing very well with the actual impedance for almost all wells. The correlation is estimated to be 0.93, RMS error is 1220 m/s*g/cc and SRE is 0.35 which reveals the good performance of the algorithm. Thereafter, in the second step, MLSSI technique is applied to the entire seismic section to estimate distribution of acoustic impedance in the subsurface. A low impedance zone ranging 6000-9000 m/s*g/cc is noticed in between 1060 and 1075 ms two-way travel time and interpreted as the reservoir zone. In the third step, multi-attribute analysis is employed and porosity along with density distributions in the subsurface are estimated for composite trace near to well locations. The results show good agreement with the actual values from well logs. The correlation coefficients are estimated to be 0.82 and 0.86 for porosity and density case, respectively, which indicate good values and hence good performance of the algorithm. Thereafter, in the fourth step, the MAA technique is applied to the entire seismic volume to estimate porosity and density volumes of the region. The predicted porosity section shows high porosity anomaly ranging 15-20% in between 1060 and 1075 ms and inverted density section shows low density anomaly ranging 2.0-2.235 g/cc in between 1060 and 1075 ms two-way travel time and confirm the presence of reservoir zone. These anomaly zones are well corroborated with low impedance zone. We can not determine if the reservoir zone is productive or not from this analysis, as only 10-20% gas in the reservoir shows tremendous change in the seismic pattern and hence in the inverted sections. The reservoir varies from North-East to South-West direction in the subsurface.