Comparison of seismic interferometry techniques for the retrieval of seismic body waves in CO₂ sequestration monitoring

Abstract

Ambient noise seismic interferometry performed by cross-correlation has been proven to be a potential cost-effective technique for geological studies. To improve the resolution of images created by interferometry, additional techniques using deconvolution and cross-coherence have been introduced. While all three methods have previously been evaluated using surface wave data for shear-wave imaging of the near surface, comparatively little study has been devoted to assess the three methods for the retrieval of body waves in reflection surveys for time-lapse application. Moreover, although the application of seismic interferometry to CO₂ sequestration by cross-correlation has been investigated by many researchers, to our knowledge, similar time-lapse studies have not been conducted using deconvolution and cross-coherence methods. We evaluate the three methods of cross-correlation, deconvolution and cross-coherence for the retrieval of phase information contained in virtual seismic records by applying seismic interferometry to synthetic data, using a model reservoir before and after CO₂ injection. By examining two approaches of regularization and smoothing factors to suppress spurious reflection events observed on the deconvolution and cross-coherence results, we note that both approaches provide similar results. We investigate noise effects by adding random noise independently at each geophone. Finally, we apply these techniques to field data recorded near the CO₂ storage site in Ketzin, Germany. For both our numerical and field data studies, we find that the cross-coherence technique retrieves the phase information of body-wave data more effectively than the cross-correlation and deconvolution techniques, and is less sensitive to uncorrelated noise from shallow sources.

interferometry, body wave, passive seismic, CO₂ sequestration

1. Introduction

Conventional ambient noise seismic interferometry (ANSI) refers to a passive seismic technique by which we can obtain the Green's function by cross-correlating noise data from two receivers (Schuster et al.2004; Wapenaar & Fokkema 2006). The ambient noise used in ANSI can be initiated by different sources such as vehicular traffic, construction, scattered body waves in the coda of teleseismic events and the dynamic change of deep petroleum reservoirs or CO₂ sequestration sites due to fluid movement and stress redistribution. ANSI has been proposed as an alternative (but cannot replace an active source in terms of imaging resolution) to active source seismic surveys for purposes such as pre-exploration evaluation, time-lapse monitoring and civil engineering, using a variety of noise sources and receiver geometries in both synthetic dataset and field studies. For example, Bakulin & Calvert (2004) and Bakulin et al. (2007) demonstrate the ANSI application to synthetic for time-lapse reservoir monitoring. Poletto et al. (2011) presented a seismic interferometry experiment using data acquired by an array of borehole receivers to obtain the reflectivity sequence. Using passive seismic noise data from a desert area, Draganov et al. (2007) obtained coherent reflection events that were comparable to the results of an active seismic survey. Cheraghi et al. (2015) imaged ore deposits in a crystalline rock environment using 300 hours of ambient noise generated by underground mining activities. Liu et al. (2016) and Zhou et al. (2015) demonstrated that using non-transient noise sources such as drilling noise with known positions one can retrieve reflection response of shallow geological structures. Additional studies have been reported that extract body waves from teleseismic signals and ambient noise for local and global velocity–structure determination (Schuster et al.2004; Roux et al.2005; Kumar and Bostock 2006; Tonegawa et al.2009; Nakata et al.2011; Olivier et al.2015; Nakata et al.2016; Cheraghi et al.2017; Mosher & Audet 2017). As a result of these studies, it is clear that ANSI has the potential to reduce the operational cost and environmental footprint that active seismic surveys have, and in particular may have application in time-lapse monitoring. Apart from the cost savings, it is worth mentioning that ANSI may provide images with lower frequencies than ever recorded in active sources.

Seismic interferometry has been used to retrieve surface waves to estimate the S-wave velocity for a wide range of applications (Larose et al., 2004, 2006; Ritzwoller et al.2005; Shapiro et al.2005; Nakata & Snieder 2012). For instance, Kimman & Trampert (2010) investigated various approximations used in ANSI and their effects on surface-wave extraction, and Halliday & Curtis (2009) studied seismic interferometry for scattered surface waves. By incrementing the aperture of the receiver, Quiros et al. (2016) were able to image deep formations from railroad noise data. Surface-wave retrieval has tended to yield more robust results than body-wave studies due to the larger amplitudes of surface waves for most geometries, largely due to the more rapid decay of body waves with distance (Xu et al.2012; Draganov et al.2013).

To enhance the resolution of both surface-wave and body-wave retrievals, several approaches have been proposed in both recording geometry and processing techniques. For instance, Mehta et al. (2007) used wavefield separation to overcome limitations of acquisition aperture that had resulted in artifacts. Mehta et al. (2008) used synthetic models to examine acquisition geometry, including allowable source spacing and source-receiver distributions. Harmon et al. (2010) demonstrated that the distribution of noise sources affects the accuracy of seismic interferometry.

Arogundade et al. (2016) analyzed several factors that can affect quality of the retrieved reflections including geophone interval, geophone depth and the presence of shallow noise source.

While cross-correlation is the conventional technique to estimate the inter-receiver Green's function required for ANSI, some other techniques including deconvolution and cross-coherence have been also proposed to retrieve body and surface waves from the ambient noise. Vasconcelos and Snieder (2008) suggested using deconvolution for interferometry, and Wapenaar et al. (2008) improve the deconvolution interferometry by introducing a multidimensional deconvolution method to overcome the limitation of source irregularity. In addition, Wapenaar et al. (2011) presented a systematic comparison between the conventional cross-correlation and proposed multi-dimensional deconvolution methods. Prieto et al. (2009) applied cross-coherence for the analysis of inelastic earth structure from seismic interferometry data. Nakata et al. (2014) applied multidimensional deconvolution to regional earthquake data to retrieve body waves. Ravasi et al. (2015) proposed multidimensional deconvolution using full-wave fields without wavefield separation. Nishitsuji et al. (2016) and Weemstra et al. (2016) used multidimensional deconvolution for noise data recorded in Malargue, Argentina. Although the three methods, cross-correlation, deconvolution and cross-coherence have been evaluated for both local and global imaging of subsurface structure from body waves (Nakata et al.2011; Nishitsuji et al.2016), no conclusive analysis has compared their results for reflection surveys for time-lapse application. In this study, we evaluate these three methods for the retrieval of body waves for time-lapse purposes from ANSI using finite-difference numerical modeling method. In addition, we show the application of the three methods for the retrieval of the body waves from field data. Our field data example can be considered as an extension of the work of Xu et al. (2012) that shows potential SNR improvements using alternatives to cross-correlation.

2. Methods

2.1. Cross-correlation method

Conventional seismic interferometry extracts Green's functions by cross-correlating wavefields excited by a noise source that has been observed by receiver pairs (Vasconcelos & Snieder 2008). Following Vasconcelos and Snieder (2008) and assuming a transient point source, with wavelet |$S( \omega )$| in the frequency domain, is emitted from a deep formation, the observed wavefields at receivers |${r_A}$| and |${r_B}$| are |$u( {{r_A},\, \omega } )$| and |$u( {{r_B},\, \omega } )$|⁠. In the Fourier domain, the cross-correlation operation is given by

$$\begin{equation} C\left( \omega \right) = {u^*}\,\left( {{r_A},\omega } \right)u\,\left( {{r_B},\omega } \right) \end{equation}$$

(1)

where the asterisk represents the conjugate operator. Using the Green's functions |$G( {{r_A},\ \omega } )$| and |$G( {{r_B},\ \omega } )$|⁠, we obtain

$$\begin{equation} u\left( {{r_A},\omega } \right) = S\left( \omega \right)G\left( {{r_A},\omega } \right)\!, \end{equation}$$

(2)

$$\begin{equation} u\left( {{r_B},\omega } \right) = \ S\left( \omega \right)G\left( {{r_B},\omega } \right)\!. \end{equation}$$

(3)

Cross-correlation can now be expressed as

$$\begin{equation} C\left( \omega \right) = {\left| {S\,\left( \omega \right)} \right|^2}\ G\,\left( {{r_A},\omega } \right){G^*}\,\left( {{r_B},\omega } \right)\!. \end{equation}$$

(4)

where |${| {S( \omega )} |^2}$| is the modulus of the autocorrelation of the source signature in frequency domain.

2.2. Deconvolution method

Snieder et al. (2006) first investigated deconvolution of waveforms for seismic interferometry. In the frequency domain, deconvolution is given by

$$\begin{equation} D\left( \omega \right) = \frac{{u\left( {{r_A},\omega } \right)}}{{\ u\left( {{r_B},\omega } \right)}}. \end{equation}$$

(5)

Applying the approach used in equations (2) and (3), deconvolution is equivalent to

$$\begin{equation} D\left( \omega \right) = \frac{{G\left( {{r_A},\omega } \right){G^*}\left( {{r_B},\omega } \right)}}{{\ {{\left| {G\left( {{r_B},\omega } \right)} \right|}^2}}}. \end{equation}$$

(6)

Note that in equation (6), the source function is eliminated in the deconvolution. This implies that the deconvolution result depends on the Green's functions between the source and two receivers, and not on the source signature.

2.3. Cross-coherence method

Aki (1957) proposed the spatial autocorrelation method, from which cross-coherence is derived. Cross-coherence, sometimes referred to as spatial coherency, applies a normalization algorithm to the spectral amplitudes of recorded signals. Prieto et al. (2009) retrieved the ambient noise Green's functions from spatial coherency. In the frequency domain, cross-coherence is given by

$$\begin{equation} H\left( \omega \right) = \frac{{u\left( {{r_A},s} \right){u^*}\left( {{r_B},s} \right)}}{{{\rm{\ }}\left| {u\left( {{r_A},s} \right)} \right|\left| {u\left( {{r_B},s} \right)} \right|}}. \end{equation}$$

(7)

Again using equations (2) and (3) applied to equation (7), we can rewrite the cross-coherence in terms of Green's functions as

$$\begin{equation} H\left( \omega \right) = \frac{{G\left( {{r_A},s} \right){G^*}\left( {{r_B},s} \right)}}{{{\rm{\ }}\left| {G\left( {{r_A},s} \right)} \right|\left| {G\left( {{r_B},s} \right)} \right|}}. \end{equation}$$

(8)

From equations (6) and (8), we note that by normalizing the amplitude, both deconvolution and cross-coherence minimize the effect of the amplitudes of the recorded noise. Indeed, both methods retrieve phase information. Particularly, cross-coherence normalizes the amplitudes across the frequency bandwidth. In other words, even if the undesired noises coming from the shallow formations have higher amplitudes than desired signals, because they usually have a broader bandwidth (due to their shorter paths to receivers and hence lower attenuation), they are balanced in deconvolution and cross-coherence. That is why cross-coherence is useful in the presence of strong undesired coherent noise.

3. ANSI numerical simulation

3.1. CO₂ sequestration time-lapse study

We create two velocity models in a two-dimensional acoustic layered medium with 1000 deep noise sources randomly placed within the deepest formation (figure 1). The two geological models shown in figure 1a and b are 10,000 m wide and 4000 m deep with 201 receivers regularly spaced on the top surface, each corresponding to a vintage of passive seismic survey, one before and one after CO₂ injection into a ridge-type structure near the center of the model at 1950 m depth. The velocity within the injected zone after CO₂ injection is modeled at 1000 m s⁻¹ in order to ensure a large effect, even though that value may be unreasonably low. Kim et al. (2015) show the P wave velocity in unconsolidated sediments of CO₂ storage site can be lower than 1000 m s⁻¹. The noise source distribution within the deep formation in this model is adapted from Boullenger et al. (2014), who demonstrate that this distribution yields optimal body-wave retrieval. In contrast, shallow body-wave noise sources will degrade the ANSI imaging by introducing ghosts, also called spurious events, compared to ‘deep noise’ that comes from below the reflector of interest and originates from, for instance, teleseismic signals or signals caused by reservoir activities. The noise sources for both baseline and repeat surveys are identical. The reservoir modeled for CO₂ storage will change impedance after CO₂ injection, as shown in figure 1b.

Figure 1.

An acoustic model of a buried ridge representing the location of a reservoir: 1000 noise sources, marked by black dots, are randomly distributed in the deep formation between 2700 and 4000 m depth, and 201 geophones are located on the free surface. (a) Before CO₂ injection. (b) After CO₂ injection, note the change in physical property at the top of the dome.

The ambient noise synthetic data is produced by a 2D finite-difference method using publicly available software developed by Thorbecke & Draganov (2011). The accuracy of the finite-difference program to calculate the Green's function in an acoustic homogeneous medium has been verified by comparison to an analytical Green's function, with a difference of less than 1% (Thorbecke & Draganov 2011). Examples of noise signals are presented in figure 2. Figure 2a shows 20 randomly chosen source noise signatures from among the 1000 noise sources used,as a time series with a sample rate of 0.008 s and varying durations. To suppress numerical dispersion, the maximum frequency of the source wavelet is set to be 30 Hz. Figure 2b and c shows details of one noise signature as a time series (part b) and amplitude spectrum (part c). We limit each source duration up to 10 s, frequencies extend up to 30 Hz and we use a total recording time of 120 s. These individual sources are randomly activated during the recording, and their signals typically overlap, creating a very long, complicated and random composite signal.

Figure 2.

(a) Examples of 20 of the random noise signals exhibiting varying source duration. (b) One source signature shown at a larger scale. (c) The amplitude spectrum of the source signature in (b) showing that the maximum frequency is 30 Hz.

To generate virtual shots gathers by seismic interferometry, we apply the three methods, cross-correlation, deconvolution and cross-coherence, to the noise data recorded by the receivers in our models before and after CO₂ injection. In the synthetic data, amplitude distribution of the noise sources is uniform meaning that all the sources have the same amplitude. Therefore, we did not apply any trace by trace normalization whereas normalization is necessary for field data because the noise distribution is not uniform. The outputs based on the three methods, which all involve the cross-correlation calculation in the numerator as shown in equations (4), (6) and (8) is a time series, that is symmetric about time zero, where we keep only the causal part. Because the source distribution is homogeneous in the model, the retrieved signals are redundant in the acausal part, and thus does not provide more information. In other words, the even illumination of the receiver array results in a more or less symmetric correlation function around time equal zero. However, this may not be the case in practice. In natural settings where illumination is often uneven, a careful illumination diagnosis would allow for a wiser choice of how to use or combine the causal and acausal parts of the correlation results. For the field data used in this study, we sum the casual and acausal retrieved image for all three methods. Figure 3 shows retrieved (virtual) common shot gathers using one single master trace, located at the middle of the recorded noise data (101th receiver, lateral position 0 in figure 1); there would be 200 other gathers in the full data set. We observe clearer reflection hyperbolae in the images retrieved from cross-correlation (figure 3a and d) and cross-coherence (figure 3c and f) compared to those obtained from deconvolution (figure 3b and e). The negative-polarity reflections are in black color. Using all 201 virtual common shot gathers extracted at each of the receivers, the stacked seismic data in figure 4 are obtained following common seismic data processing procedures. The stacking velocity is computed by picking the time-velocity pairs from the semblances of each of the CMP gathers. More details of the processing of the data can be found in Cao (2016).

Figure 3.

Examples of one retrieved common shot gather at the location of middle receiver (101st receiver among the total 201 receivers) after cross-correlation, deconvolution and cross-coherence, respectively, (a)–(c) before CO₂ injection and (d)–(f) their counterparts after CO₂ injection.

Figure 4.

(a)–(c) Interferometric images retrieved by cross-correlation before and after CO₂ injection and the difference between the two surveys of (a) and (b). (d)–(f) Interferometric images correspond to (a)–(c) produced by deconvolution. (g)–(i) by cross-coherence. An artificial reflection around 0.7 s is highlighted by the black ovals. The two-way travel time examples for major reflections are marked in (a).

Figure 4 shows the final, stacked images before and after CO₂ injection, and their difference (‘after’ minus ‘before’ injection) obtained from cross-correlation, deconvolution and cross-coherence. The prominent reflections retrieved at the two-way travel times of 0.533, 1.233, 1.783 and 2.117 s correspond to the reflections from interfaces in the geological model at 400 m (water bottom), 1100, 2200 and 2700 m. In addition, we can clearly identify at around 1.658 s, in the center of each section, corresponding to impedance changes at the top and bottom of the CO₂ reservoir that crests at 1950 m depth. In this case, with no random noise added, we note that the stacked sections and the difference results all recreate the original model, and that the cross-correlation method seems to be most noise-free. In particular, our cross-coherence results are in agreement with an experimental (laboratory-scale) study conducted by Draganov et al. (2016) who demonstrated the validity of cross-coherence for observing CO₂ injection. However, we point out that there are numerical artifacts in the images produced by the deconvolution and cross-coherence methods prior to differencing (such as those highlighted by the black ellipses). In addition to the numerical artifacts, the polarity reversed peg-leg multiples are also generated because the geophones are placed on the free surface.

3.2. Regularization factor and smoothing

In the usual case of the bandlimited recording of ambient noise data, deconvolution (equation 6) and cross-coherence (equation 8) can become numerically unstable when the recorded signal (in the denominators) goes to zero at one or more frequencies. This instability introduces some of the artifacts that we observe in figure 4d, e and f (for deconvolution) and 4g, h and i for cross-coherence. In order to overcome this instability, a parameter can be added providing non-zero values in the denominators. The value to be added can be controlled by a regularization factor (ε) referring to a percentage of the mean value of the spectral power of |$G( {{r_B},\ s} )$| at the reference receiver (Clayton & Wiggins 1976; Mehta et al.2007). The equation for regularized deconvolution interferometry is given by (Mehta et al.2007)

$$\begin{equation} D\left( \omega \right) = \frac{{G\left( {{r_A},s} \right){G^*}\left( {{r_B},s} \right)}}{{{{\left| {G\left( {{r_B},s} \right)} \right|}^2} + \varepsilon \left\langle {{{\left| {G\left( {{r_B},s} \right)} \right|}^2}} \right\rangle }}, \end{equation}$$

(9)

where |$\ {\langle| {G( {{r_B},\ s} )} |^2\rangle}$| is the average power spectrum of |$G( {{r_B},\ s} )$|⁠.

The expression for regularized cross-coherence interferometry is given by

$$\begin{equation} H\left( \omega \right) = \frac{{G\left( {{r_A},s} \right){G^*}\left( {{r_B},s} \right)}}{{\left| {G\left( {{r_A},s} \right)} \right|\left| {G\left( {{r_B},s} \right)} \right| + \varepsilon \left\langle {{{\left| {G\left( {{r_B},s} \right)} \right|}^2}} \right\rangle }}. \end{equation}$$

(10)

In this study, to obtain a reasonable value for ε, first we use cross-correlation as a benchmark because its output appears to contain no such artifacts, and we compare the results of deconvolution and cross-coherence with it. Then, we choose a small value for ε and subjectively assess whether the artifacts have been suppressed; if we observe that the artifacts persist, we iteratively increase ε until the artifacts become visually negligible. Figure 5 shows examples of a virtual shot record obtained from deconvolution and cross-coherence with varying values of|${\rm{\ }}\varepsilon $|⁠. We note that the signal quality (based on the similarity to the original cross-correlation virtual shot records of figure 3) improves as ε increases. However, as ε increases, the retrieved images look increasingly identical comparing figure 5c and f where ε applied is 10%.

Figure 5.

(a)–(c) Virtual shot gathers produced by deconvolution with the regularization factor ε equal to 0, 1 and 10%, respectively and (d)–(f) produced by cross-coherence.

The regularization parameter needs to be chosen carefully. If the value of ε is too small, the instability persists and the artifacts are not suppressed. On the other hand, if a too large value of ε is selected, the retrieved images from both methods become very similar to each other and to the cross-correlation results. Therefore, we need to find a value of ε that is optimal: reducing the artifacts while providing better images than cross-correlation.

We also investigate denominator spectral smoothing method to remove artifacts. This smoothing uses a moving window in the frequency domain to make the rough spectrum more smooth (Prieto et al.2009). A comparison of various smoothing techniques is beyond the scope of this paper. The smoothing approach used here is a Gaussian-weighted moving-average filter, resulting in better stopband attenuation than the simpler moving-average filter (Zhang & Fu 2016). The weighted values are calculated as shown in equation (11):

$$\begin{equation} G_i^{\prime} = \frac{{\mathop \sum \nolimits_{n = - k}^k {G_{i + n}}{W_{i + n}}}}{{\mathop \sum \nolimits_{n = - k}^k n}}, \end{equation}$$

(11)

where |${G_i}$| is the value at the center of the selected window of size |$2k + 1$|⁠; |${W_i}$| is the Gaussian weight for each of the neighborhood values in one window. And |$G_i^{\prime}$| represents the weighted values after Gaussian average.

The standard deviations and mean values used for calculating the Gaussian function values are updated as the smoothing window moves to each successive sample point. The values that significantly deviate from the mean are assigned a smaller weight than those values near the mean. Figure 6 shows the weighted-smoothing denominator power spectrum of one reference trace for deconvolution and cross-coherence methods with window lengths of 20 and 200 samples (0.16 and 1.6 s). The power spectrum is plotted in dB scale; recall from figure 2c that the noise signals are constrained to frequencies below 30 Hz, and have randomly assigned amplitudes at each frequency. As the window length for smoothing increases, the spectrum is increasingly smoothed.

Figure 6.

Raw (gray lines) versus smoothed (black lines) denominator power spectra for the window lengths of W_l = 20 (a and c) and W_l = 200 (b and d) for deconvolution (top panel) and cross-coherence (bottom panel), respectively.

Figure 7 shows a set of examples using denominator spectrum smoothing, with different moving-window lengths and meanwhile without any regularization factor, to obtain a single virtual shot gather using both deconvolution and cross-coherence methods. As we can see, the artifact around 0.6 s is suppressed as the smoothing window length increase. We note, once again, that very large values tend to result in images that are very similar between the two methods, and that continuing to increase the window length beyond some value (here about 100 to 200) does not noticeably improve the image.

Figure 7.

The smoothing results by different window lengths (W_l), (a) and (d) W_l = 0, (b) and (e) W_l = 20, and (c) and (f) W_l = 200 for deconvolution and cross-coherence, respectively. For all processing, ε|$\ = \ 0$|⁠.

For a comparison between the regularization and moving-window smoothing results, we selected traces 80, 100 and 120 to show the differences (trace 100 is adjacent to middle trace among total 201 traces) and normalize them using roughly optimal values (ε = 0.05) for each approach (figure 8). For spectral smoothing, we used values for window lengths of 20 and 200, and for regularization we used ε = 0.05. We chose the value ε = 0.05 by comparing figure 5b and c and figures 5e and f, respectively, where ε = 0.01 and 0.1, and by examining multiple simulations that showed ε = 0.05 is sufficient to suppress the artifact. We note that traces retrieved by the moving-window approach tend to be slightly noisier than those obtained by regularization. Therefore, in the remainder of this study, we apply regularization in our processing. However, one could consider trying a combination of the two approaches to suppress artifacts.

Figure 8.

(a) Retrieved seismic trace at the location of receiver 80 by cross-coherence method with smoothing window W_l = 20, W_l = 200 and regularization factor i = 0.05, respectively. (b) and (c) show the same for receivers 100 and 120.

3.3. Signal to noise ratio analysis

Ambient noise recorded by a geophone can be imagined as a superimposition of two (or more) noise types. One type is the noise coming from deeper layers (below the depth range being imaged) that is generated, e.g. from the reservoir's activities or from teleseismic events. Boullenger et al. (2014) demonstrated that this noise type (called ‘deep noise’ in the remainder of this paper) has a constructive contribution to ANSI imaging. The other type of noise considered here is that which does not have a constructive effect on the retrieval of seismic waves and, instead, is destructive in the process to retrieve useful seismic data; this is usually shallow in nature. The location of a noise source and the geometry of its ray path are the two primary factors determining whether noise has constructive or destructive effects on retrieving estimates of the inter-receiver reflection response. Wu (2016) showed that in contrast to deep noise, near-surface noise degrades ANSI imaging. Because near-surface noise will appear different at the different receiver locations, due to large differences in travel path, we will include the effect of this destructive noise type through the addition of Gaussian noise to the individual records made of deep noise. In this paper, we consider the deep noise data as coherent noise (due to its positive contribution to the ANSI) and the Gaussian noise as random noise. (Of course, all of it is ‘noise’ in the commonly used sense in that it is not purposely created in order to provide a known signal.)

To evaluate the efficacy of the cross-correlation, deconvolution and cross-coherence methods for the retrieval of seismic images, we consider different scenarios through the addition of white Gaussian noise to the records made of deep noise. Here we define ‘signal’ to ‘noise’ ratio (SNR) by

$$\begin{equation} SNR = \frac{{\sigma _{signal}^2}}{{\sigma _{noise}^2}}, \end{equation}$$

(12)

where |${\sigma _{signal}}$| and |${\sigma _{noise}}$| are, respectively, the standard deviation of the recorded coherent noise from the intended sources (coming from the deep formation) and random noise (Gaussian noise applied at each receiver, intended to simulate destructive shallow noise sources). Given the mean amplitude and standard deviation of the coherent noise for each of the 201 traces, random noise for each trace was generated, based on equation (12) to achieve selected SNR values. A new passive data set is created by adding the additional random noise to the coherent noise from the deep formation.

Because we add white random noise to the coherent noise, the spectrum is nearly white and the problem of artifacts, addressed in the previous section, does not appear. In our first case, we consider an acquisition site at which the receivers are uniformly affected by random noise. In this scenario, the SNR is assumed to be 1 for all the geophones. In other words, the amount of the added Gaussian random noise is the same as the deep-origin coherent noise for all receivers. Figure 9 shows the final ANSI stacked images that are obtained from the cross-correlation, deconvolution and cross-coherence methods. In this case, we observe that the cross-correlation and cross-coherence methods provide better images of the subsurface than the deconvolution method.

Figure 9.

Stacked interferometric images by (a) cross-correlation, (b) deconvolution and (c) cross-coherence methods. The signals were recorded by 201 receivers with the signal to noise ratio equal to 1.

In our second case, we assume that the SNR randomly varies from one geophone to another. This random noise simulation is reasonable as in real data acquisitions, geophones record different level of noise due mostly to their coupling and their locations relative to nearby sources of noise. To simulate this situation, we consider different SNR values for each recorded geophone, ranging from 0.1 to 1, with that value randomly varying from one geophone to another. From the results shown in figure 10, using the three processing methods, we observe that the cross-coherence method (figure 10c) presents a stronger image. For both cross-correlation and deconvolution methods, the reflections are contaminated by random noise at different receivers, drowning out the signal at the locations where the SNR is low.

Figure 10.

Stacked interferometric images by (a) cross-correlation, (b) deconvolution and (c) cross-coherence methods. Random single to noise ratio for the input data varies from 0.1 to 1.0.

In our final case, we consider an acquisition environment in which the SNR increases gradually across the line of acquisition. This is intended to simulate situations where, for instance, one end of the receiver array is close to a highway or to production facilities, degrading the data recorded by the receivers closest to the noise source. We assume the SNR linearly increases from left (SNR = 0.01) to right (SNR = 1). Again, by comparing the results (figure 11), we observe that cross-coherence has obtained a better image of the subsurface, distinguishing the reservoir structure. The image obtained by cross-correlation is poor throughout, and strong upward-curved features on left side of the image are probably caused by over-correction of NMO. Although the deconvolution method does bring out all the main reflections through most of the section, its resolution is not as good as it is in the images processed by cross-coherence.

Figure 11.

Stacked images by (a) cross-correlation, (b) deconvolution and (c) cross-coherence methods. 201 receivers with the signal to noise ratios increasing from 0.01 to 1.0 from left to right in each image. Oval identifies the location of the buried ridge in the model.

From our synthetic data simulations, we observe that cross-coherence is sufficiently robust to image the subsurface and less sensitive to the strong undesired surface noise than the other techniques. The artifacts produced by the deconvolution and cross–coherence methods in figure 4 are eliminated when bandlimited random noise has been added. We would expect that no regularization would be required if the recorded ambient noise had a wide and smooth frequency spectrum, a situation we did not model.

3.4. Influence of random noise from surface waves on time-lapse study

In this paper, we have already discussed the application of ANSI to seismic imaging using noise sources uniformly distributed in a deep subsurface without surface wave contamination, and have demonstrated that a limited amount of surface noise with random amplitude can serve to suppress artifacts through their spectral-whitening effect. Of course, this would not be the case for surface noise that is predominantly in a narrow frequency band, as many traffic or industrial noise sources would tend to be. We now move on to a more realistic situation for CO₂ sequestration time-lapse study in which both the base and repeat surveys are recorded by merging the deep noise and specific random noise. We use the results of the studies presented earlier in this paper, employing only the optimal image-enhancing parameters for each case, but now we include the time-lapse difference section as well. In doing so, we examine the possibility that poorly imaged individual ‘before’ and ‘after’ cases may still produce a well-imaged difference section.

First, we examine the case previously presented in figure 10 – that of a random distribution of noise across the surface receiver array (figure 12), for all three methods examined in this paper. In figure 13 we examine the case (as in figure 11) for linearly increasing random noise across the receiver array. From these two simulations, it becomes apparent that cross-coherence method maintains its advantage in the time-lapse study.

Figure 12.

Stacked sections generated from the passive noise data with random signal-to-noise ratio ranging from 0.01 to 1.0 between input traces. The first row of images (a–c) were processed by cross-correlation, while the second and third row of images (d–f), (g–i) Images processed by deconvolution and cross-coherence, respectively. The left panels (a, d and g) show the results of the base survey; middle panels (b, e and h) after CO₂ injection and the difference is shown in the right panels (c, f and i).

Figure 13.

Stacked sections generated by the passive noise data with increasing signal-to-noise ratio ranging from 0.1 to 1.0 from left to right in each section. The left panels (a, d and g) show the results of the base survey; middle panels (b, e and h) after CO₂ injection and the difference is shown in the right panels (c, f and i).

3.5. Low frequency noise source

One limiting factor in the use of ANSI for imaging the subsurface is the bandwidth of the deep noise sources. In the case of a low frequency deep noise source (i.e. 5 Hz), the retrieved image would be expected to have poor vertical and lateral resolution. To investigate the effect of bandwidth on ANSI, we model passive coherent noise (from the deep formation) with a maximum frequency of 5 Hz (details in figure 14), for comparison with the 30 Hz maximum used earlier. In this case, we do not add any additional surface random noise. Figure 15 shows the stacked images obtained from cross-correlation, deconvolution and cross-coherence. We observe that the temporal and vertical resolutions of the stacked images generated from deconvolution and cross-coherence are much higher than those from cross-correlation. While the response of cross-correlation is a function of the source spectrum |$S( \omega )$|⁠, the responses of the deconvolution and cross-coherence are solely controlled by the Green's functions and not the source spectrum. This demonstrates a strong advantage of the deconvolution and cross-coherence methods over the cross-correlation method.

Figure 14.

(a) Example of 20 random noise signatures with varying source duration. (b) One source signature shown on a large scale. (c) The amplitude spectrum of the source signature in (b) showing that the maximum frequency is 5 Hz.

Figure 15.

The stacked image of seismic interferometry produced by (a) cross-correlation, (b) deconvolution and (c) cross-coherence. The maximum frequency of source noise recorded from the deep formation is 5 Hz.

4. Field data

We now apply these three methods to field data acquired at the Ketzin CO₂ storage site (Germany) in February 2011 (Xu et al.2012). Unfortunately, this dataset does not include a time-lapse component, but we include the dataset in this paper because it provides a real-world test of our methods. The dataset consists of 63 traces from vertical geophones spaced at an interval of 24 m, sampling rate of 4 ms and recording time of 28.7 h. The vertical component of the ground motion recorded by the geophone was split evenly into 3348 ‘panels’ of 30 s each. The noise sources were apparently comprised of natural earthquakes, wind and vehicles, and observations suggested that surface waves were large in many time panels. Xu et al. (2012) indicate that the surface waves are present mostly below 14 Hz and the body waves are observed between 12 and 26 Hz.

4.1. Pre-processing

The raw data for one noise panel is shown in figure 16a in which we can identify particularly high amplitude traces. Although these are likely to be simply ‘bad’ traces, we treated them as if they were valid, but recorded with too high a gain in the event that they actually contained useful data. Before the cross-correlation, deconvolution and cross-coherence processing, root mean-square (RMS) normalization was performed to balance the energy between the traces. Other techniques can be applied: Bensen et al. (2007) summarized the workflow for seismic ambient noise data processing and compared a variety of normalization methods including one-bit, clipped waveform, automated event detection and running absolute mean normalization for ambient noise data pre-processing. Although the length of time window used for calculating the RMS for normalization can affect the results, such comparisons are beyond the scope of this paper, and a seemingly satisfactory window (30 s) was selected and used for the complete data set. After normalization, we obtained a frequency-wavenumber plot for the complete set of 3348 noise panels (figure 16c). Recognizing that most of the energy is contained in a limited frequency range, we applied a simple 2–30 Hz bandpass filter.

Figure 16.

(a) The raw data of 63 traces from one time panel. The traces of anomalous amplitude are marked by the red rectangle. (b) The noise signal data after normalization using the RMS method. (c) The frequency-wavenumber spectrum of the normalized noise data prior to filtering.

4.2. Seismic response retrieval for field data

Because the passive field recordings at Ketzin were made only following CO₂ injection, and no baseline seismic survey is available, we simply focus on the comparison of retrieved post-injection ANSI images obtained from cross-correlation, deconvolution and cross-correlation.

For storage of the data, the seismic data recorded in the Ketzin site had been divided into short-time panels. Because these panels are short, there is a very limited range of useful coherent noise in each panel, and reflections are not identifiable on the retrieved seismic response, regardless of technique used (figure 17). The regularization factor applied to deconvolution and cross-coherence is 1% of the average value of the denominator, a value found to be optimal in the synthetic cases. To improve the signal to noise ratio for the retrieved response, after first applying ANSI to each of the 3448 panels we then stack the 3448 outputs together for each virtual shot gather location. One such virtual shot gather, retrieved at a central location, is shown (figure 18a–c) as processed by the three methods. This can be compared with the common shot gather shown from synthetic data in figure 3; we recognize that the steeply dipping straight-line events must represent ground roll (with low frequency energy and velocity), and any flatter events that may be discerned could represent body waves. In order to further reject the surface waves and retain only the body waves, a 12–26 Hz bandpass filter, which is decided by trying different bandwidths, was applied (figure 18d–f) after a rough bandpass filter 2–30 Hz in the pre-processing. We see that the surface waves seem to have been largely removed, but many other straight-line events, probably trapped waves, cross the gathers. The left-side traces that appeared to contain only large-amplitude noise seem to remain contaminated, particularly for the cross-correlation and deconvolution methods.

Figure 17.

(a) (b) and (c) represent one single, unstacked, retrieved shot gather at the location of the middle receiver (32nd along the line of 63 receivers) using only single one ‘panel’ noise recording, through cross-correlation, deconvolution and cross-coherence, respectively.

Figure 18.

The virtual common shot gathers retrieved at location of receiver 30 after summing all 3348 panels by (a) cross-correlation, (b) deconvolution and (c) cross-coherence methods. A bandpass filter of 12–26 Hz is applied to extract body waves, and (d), (e) and (f) correspond to the filtered output from (a), (b) and (c), respectively. Traces that appear to contain only large-amplitude noise in the cross-correlation and deconvolution displays are indicated in the black ellipses.

Because most of the noisy panels are dominated by surface or trapped waves both before and after bandpass filtering, this data set requires additional processing beyond that which we applied to the synthetic data examples. To make the hyperbolic reflections more identifiable and taking advantage of the known reasonably flat structure based on a previous study of the area from both passive and active surveys (Xu et al.2012), we first constructed a single supergather by sorting all 63 common shot gathers (filtered to 2–30 Hz) to common offset gathers and merged the traces with common offsets (figure 19a, d and g). We then applied an f-k filter to reject the surface velocity range 400–1600 m s⁻¹ (figure 19b, e and h), followed by spectral whitening over 12–26 Hz to equalize the amplitudes of all frequencies (figure 19c, f and i), and finally a zero-phase bandpass filter (12–26 Hz) to remove inadvertently introduced noise outside the band of interest. As a result, the surface waves were largely removed from the supergather, whether derived from cross-correlation, deconvolution or cross-coherence. We observed that the body waves were apparently retrieved more effectively by cross-coherence, as observed cresting at 0.4 and 0.5 s in figure 19i. Note that the processing scheme of the real data is only slightly different from what was used for the synthetic data because the field data is contaminated with strong surface wave arrivals and non-uniform noise distribution. Therefore, we had to apply additional processing steps.

Figure 19.

The supergathers formed by merging by all 63 common shot gathers by (a) cross-correlation, (d) deconvolution and (g) cross-coherence method. The f-k domain (b), (e) and (h) is filtered to remove the highlighted areas, and applied to extract body waves. Following the f-k filtering step, spectral whitening was applied over 12–26 Hz to obtain the supergathers in (c), (f) and (i).

Our purpose in constructing the single survey-wide supergather in figure 19 was to identify the retrieved hyperbolic reflections by improving the SNR. Having done this, we recognize that there are, indeed, reflections present in the data. The following processing sequence was then performed on the individual virtual shot gathers in order to extract a stacked seismic section.

Figure 20 shows the final stack of the seismic sections obtained from the cross-correlation, deconvolution and the cross-coherence methods using the Ketzin field data. First, we identified 68 panels out of 3348 total recorded noise by visual inspection, representing about 0.56 out of 25 h recording time in which surface waves were not dominant and which contained apparently identifiable body waves. The method we used to visually identify the body waves from the noise panels is explained in Xu et al. (2012). An alternate method using illumination diagnosis for body-wave identification is explained in Vidal et al. (2014). We then applied cross-correlation, deconvolution and cross-coherence to each reduced-panel data set and merged all the virtual shot gathers retrieved at the same reference receiver. Finally, we applied seismic data processing as outlined in Table 1. The resulting stacked images are shown in figure 20. To avoid the very low signal to noise ratio of the stack data at the survey edges due to the low fold, we focus our interpretation only on the central part of the section (15–45 out of the total 61 common depth points (CDPs)).

Figure 20.

Comparison of the stacked section between (a) cross-correlation, (b) deconvolution and (c) cross-coherence. The arrows indicate possible reflecting events at around 0.2 and 0.5 s on the sections where they may be apparent.

Table 1.

Processing steps applied to the Ketzin virtual shot gathers for all the three methods.

Step	Processing
1	Sort to CMP
2	Velocity analysis
3	NMO
4	Stack
5	F-X deconvolution
6	Automatic gain control
7	Spectrum whitening
8	Bandpass filter (12–26 Hz)

A disappointing feature of these images is the ringiness (narrow-band reverberations). This was also observed on an active source data survey (Xu et al.2012) when filtered to match the ANSI frequency band. However, the active data we received, and display in figure 20d, has been filtered such that reflections, beyond the frequency range of the ANSI data, are visible. The deconvolution and cross-coherence results do not exhibit a significantly sharper wavelet than the cross-correlation results, as expected from equations and numerical results. Because we applied a bandpass filter to remove the surface waves, the data became bandlimited. Therefore, we applied spectral whitening to all the three techniques. We have not attempted to reduce the ringing in our ANSI results by additional wavelet processing techniques.

Based on geological data and the results of the active source seismic surveys (Xu et al.2012), we would expect to find major reflectors at the two-way travel times of 0.2, 0.5 and 0.9 s. In figure 20, we see that each processing method results in a different apparent identification of reflected events. All three methods (figure 20a–c) appear to retrieve the shallowest reflector (0.2 s) and the reflector at around 0.5 s, which matches well to the Tertiary and K2 horizon, respectively, according to Xu et al. (2012) and Boullenger et al. (2014). However, deconvolution and cross-coherence render slightly better images for both reflectors. The spectral whitening applied to the virtual shot gathers and final stacked data sharpened the reflections, but it also boosted the noise. Nevertheless, cross-coherence showed more stability at retrieving reflection events based on figure 20 if one compares the coherence of the retrieved signal at the shallow reflection near the red arrow and at the far offset. For the cross-correlation and deconvolution methods, the FX-deconvolution could also produce artifacts if the traces are not balanced in each virtual shot gather (e.g. figure 18a and b). Recall that we selected only those noise panels with identifiable body waves (about 30 min out of 25 h); if the recording time had been longer, the stacked result would have likely been improved.

5. Conclusion and discussion

We showed, using synthetic examples of ANSI, in addition to the cross-correlation technique previously applied in time-lapse monitoring studies, that the deconvolution and cross-coherence may in some cases provide more useful images than the cross-correlation approach. Modeling the noise sources used in the synthetic studies as bandlimited noise originating in the deep formation, as may occur during the process of CO₂ sequestration, the deconvolution and cross-coherence methods tend to generate numerical artifacts produced by the division in the frequency domain. These artifacts could be reduced using a regularization parameter or spectral smoothing, rendering reasonably good images from both deconvolution and cross-coherence.

In addition, we studied potentially destructive sources of noise. The worst of these are likely to be shallow sources because of the strong differences in raypaths to neighboring receivers. To simulate such noise, we added Gaussian noise to the recorded ambient noise, considering three different scenarios that we might face in the field; (i) all geophones record the same amount of destructive noise; (ii) some geophones are randomly more strongly affected by the noise and (iii) the level of destructive white noise gradually increases (decreases) along the line of acquisition. In all scenarios, cross-coherence resulted in better images. Both cross-correlation and deconvolution techniques were found to be more susceptible to noise from shallow sources.

Furthermore, we investigated the effect of bandwidth of the deep noise source. We noticed that deconvolution and cross-coherence provided greater temporal and lateral resolution because their outputs were controlled only by the Green's functions and not the source bandwidth.

The Ketzin field data contained extremely strong surface-wave contamination, presenting a new challenge. Without additional pre-processing, no image of value could be obtained. But using only those noise panels in which body waves could be visually observed over the noise provided a small suite of data from which some images could be obtained. Again, cross-coherence appeared to image the reflections better than the other techniques. Because cross-coherence was independent of the signal-noise source signature, this method appeared to be suitable for noise data with strong amplitude variations and complex signal to noise waveforms.

Acknowledgements

We thank to J. Thorbecke for providing his finite-difference wavefield modeling code; C. Juhlin and S. Lüth for sharing the Ketzin field data and W. Pennington for useful discussions and assistance in editing this paper. Most of conventional data processing was performed by Seismic Unix. We hereby express our gratitude to the Seismic Un*x project at the Center for Wave Phenomena, Department of Geophysics, Colorado School of Mines. We are grateful to P. Persaud for her valuable suggestions, which greatly improved our paper.

Conflict of interest statement. None declared.