Introduction
The Boom Clay in Belgium is composed of clays and generally has embedded boulders and/or thin layers of septaria (according to the Merriam Webster dictionary, a septarium is ‘a concretionary nodule usually of limestone or clay ironstone intersected within by cracks filled with calcite, barite, or other minerals’). The formation of these concretions (lithification) started early in the diagenesis process near the sediment–water interface and before compaction occurred (De Craen et al., Reference De Craen, Wang, Van Geet and Moors2004). The septaria generally are made of ellipsoidal-shaped carbonate concretions, whose sizes range from a few decimetres up to a couple of metres (see Fig. 1). The seismic response of these thin layers is characterised by a set of parallel reflection events, which, depending on the frequency bandwidth, look like continuous horizontal events or alignments of diffraction events (Fig. 2) (see Hemerijckx et al., Reference Hemerijckx, Carpentier, De Schrijver, Henriet and Heldens1983; his Fig. 2).
A calcareous septarium embedded in Boom Clay. The width is approximately 60 cm (from Vis & Verweij, Reference Vis and Verweij2014).
A. In-line section showing a septaria level; B. Close-up with single diffraction events; C Horizon slice at 0.2 ms below that level; the green arrows refer to events identified in B; D. Enlarged section showing typical blotchy pattern possibly suggesting concretion distribution (from Missiaen et al., Reference Missiaen, Versteeg and Henriet2002).
The septaria may constitute a problem regarding the integrity of the Boom Clay and therefore it is essential to locate these geobodies, preferably with seismic waves. However, the fact is that the presence of high attenuation often degrades the seismic images, resulting in a low resolution of targets. In this case the conventional migration may not image the geobodies properly.
Recently, this problem has been reconciled with an advanced migration with attenuation compensation. For example, Wang & Guo (Reference Wang and Guo2004) presented a migration algorithm incorporating inverse Q filtering to improve the imaging resolution, although the imaging technique is limited to 1D velocity and attenuation models. Zhang et al. (Reference Zhang, Zhang and Zhang2010) proposed a reverse-time migration (RTM) scheme for compensating attenuation and phase dispersion effects. Their viscoacoustic wave equation is based on a constant-Q model, that is, attenuation is considered to be approximately linear with frequency. Dutta & Schuster (Reference Dutta and Schuster2014) used a linearized 2D viscoacoustic wave equation based on the Zener model, written in the particle velocity–stress formulation. The method is adapted from conventional least-squares migration and reconstructs the Earth's reflectivity image from the recorded wavefield under the Born approximation. Zhu (Reference Zhu2014) implemented a time-reverse modelling approach (or reverse-time focusing) based on a viscoacoustic wave equation which explicitly separates attenuation and dispersion following a constant-Q model based on fractional derivatives (Carcione, Reference Carcione2010). Later this idea was formulated in an RTM algorithm (Zhu et al., Reference Zhu, Harris and Biondi2014). Because of possible instability problems, attenuation compensation is necessary to stabilise the migration. Sena et al. (Reference Sena, Stoffa and Sen2006) introduced a plane-wave approximation in the split-step Fourier technique in the frequency domain to stabilise the algorithm in the presence of attenuation. Zhu et al. (Reference Zhu, Harris and Biondi2014) and Sun et al. (Reference Sun, Zhu and Fomel2015) adapted a low-pass filter to suppress the high-frequency noise during the compensation.
In this work, we propose to use a post-stack RTM (e.g. Baysal et al., Reference Baysal, Kosloff and Sherwood1983) taking into account the attenuation effects to recover the migration amplitude of septaria geobodies. We use the simple exploding-reflector concept to back-propagate the seismic waves. To extrapolate source and receiver wavefields with attenuation compensation, the equation used for back-propagation of the signals has the same velocity but a negative Q factor. We solve a lossy wave equation that is dispersionless. In the migration process, the recorded seismogram is used as a time-dependent boundary condition. The seismic trace is applied at each receiver in reversed time and the propagation goes back in time until the origin time, where the best focusing occurs. The reflectors are thus imaged. We also apply a low-pass filter to suppress the high-frequency noise.
The dispersionless wave equation
Wave motion (P waves) in the (x, z) plane of an acoustic medium is governed by the following equation:
$$\begin{equation}
{\partial _x}({\rho ^{ - 1}}{\partial _x}p) + {\partial _z}({\rho ^{ - 1}}{\partial _z}p) = {(\rho {c^2})^{ - 1}}{\partial _{tt}}p + s(x,z,t)
\end{equation}$$
(e.g. Carcione Reference Carcione2015, Chapter 9), where p is the pressure field, ρ is the mass density, c is the wave velocity, s is a source term, ∂ i is the partial derivative with respect to the spatial variable xi (x or z) and ∂ tt is the second-order partial derivative with respect to the time variable t. This equation can be rewritten as
$$\begin{eqnarray}
&&\!\!\!{\partial _t}\left( {\begin{array}{*{20}{l}} p\\ {\dot p} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - \gamma }&\ \ 1\\ {\rho {c^2}[{\partial _x}({\rho ^{ - 1}}{\partial _x}) + {\partial _z}({\rho ^{ - 1}}{\partial _z})]}&\ \ { - \gamma } \end{array}}\! \right)\!\left( {\begin{array}{*{20}{l}} p\\ {\dot p} \end{array}} \right) + \left( {\begin{array}{*{20}{l}} 0\\ {{s^\prime }} \end{array}} \!\right)\!,\nonumber\\
&&
\end{eqnarray}$$
where
$\dot p = {\partial _t}p$
(if γ = 0), s′ = ρc
2
s and γ = 0 yields equation (1). The introduction of γ implies wave attenuation, as we shall see below.
Eliminating the variable
$\dot p$
in equation (2), we obtain
$$\begin{equation}
{\partial _{tt}}p = \rho {c^2}[{\partial _x}({\rho ^{ - 1}}{\partial _x}p) + {\partial _z}({\rho ^{ - 1}}{\partial _z}p)] - 2\gamma {\partial _t}p - {\gamma ^2}p + {s^\prime }.
\end{equation}$$
We analyse the anelastic properties related to equation (3) in Appendix A. As shown here, a negative quality factor (implicit in equation (3)) is required to back propagate the wavefield for recovering the amplitudes in reverse-time imaging (Zhu, Reference Zhu2014). Thus, the reverse-time imaging equation is written as
$$\begin{eqnarray}
{\partial _{tt}}p = \rho {c^2}[{\partial _x}({\rho ^{ - 1}}{\partial _x}p) + {\partial _z}({\rho ^{ - 1}}{\partial _z}p)] + 2\gamma {\partial _t}p - {\gamma ^2}p + {s^\prime }.\nonumber\\
&&
\end{eqnarray}$$
In Appendix B we derive finite-difference versions of these equations to perform forward and backward propagation.
An exploding-reflector version of equation (3) can be obtained by assuming the impedance ρc constant in the whole space. In this way, multiples are avoided. An approximation to a stacked section is obtained in a single experiment by halving the velocities and by initiating the sources at time zero on all the reflecting boundaries. With this model, the recorded surface time section approximates the stacked or zero-offset section (Carcione et al., Reference Carcione, Böhm and Marchetti1994, Reference Carcione, Feliciangeli and Zamparo2002). In a forward simulation, the strength of each source at the interfaces is proportional to the normal-incidence reflection coefficient. Equation (3), for the exploding-reflector approach, becomes
$$\begin{equation}
{\partial _{tt}}p = c[{\partial _x}(c{\partial _x}p) + {\partial _z}(c{\partial _z}p)] - 2\gamma {\partial _t}p - {\gamma ^2}p + {s^\prime }.
\end{equation}$$
Examples
First, we consider the model displayed in Fig. 3, where septaria layers and boulders are embedded in homogeneous Boom Clay. The upper thin layer is composed of septaria bodies of 1 m size separated with a period of 4 m. The period is 10 m for the lower layer. Moreover, there are isolated septaria at different depths. The seismic velocity of the Boom Clay (homogeneous) is c = 1.8 km/s and the attenuation parameter is Q 0 = 20 on the left-hand side of the model (x < 100 m) and Q 0 = ∞ on the right-hand side. We perform an exploding-reflector simulation based on equation (4) using a Ricker wavelet with a peak frequency fp = 200 Hz. The waveform is
$$\begin{equation}
w(t) = \left( {u - \frac{1}{2}} \right)\exp ( - u),\;\;u = {\left[ {\frac{{\pi (t - {t_s})}}{T}} \right]^2},
\end{equation}$$
where T = 1/fp is the period of the wave and we take ts = 1.4T.
Septaria bodies and layers embedded in homogeneous Boom Clay. The upper layer is composed of septaria bodies of 1 m size separated with a period of 4 m. The period is 10 m for the lower layer. Moreover, there are isolated septaria at different depths. The Q factor of the Boom Clay at the left-hand side is 20 (surrounded by a box with dashed grey lines), while the right-hand zone is lossless.
We consider a numerical mesh with nx = nz
= 231 grid points and a grid spacing dx = dz = 1 m. Absorbing boundaries of size 24 grid points are implemented at the sides of the mesh. The algorithm for forward modelling, based on the Fourier pseudospectral method, uses a time step of 0.1 ms to propagate the wavefield 2000 steps. Fig. 4 shows the synthetic seismogram, where it is clear that the signal has been highly attenuated on the left-hand side. On the right-hand side the response of the layered sets of septaria is continuous and single geobodies cannot be distinguished. This is a Fresnel-zone effect (see Appendix C). According to the equations reported in Appendix C and the properties used here, the Fresnel radii before and after migration at the peak frequency are
$R = 2.12\sqrt z $
and R = 2.25 m, respectively. The first and second layers of septaria have z = 70 m and 115 m, implying that before migration the Fresnel radii are 18 m and 23 m, respectively. The geobodies therefore cannot be discriminated in the unmigrated section, since their separations are 4 and 10 m, respectively. The migration without (a) and with (b) Q compensation is displayed in Fig. 5, where it is clear that the amplitude has been recovered (the narrow left vertical zone is within the absorbing boundary, implemented to avoid wraparound). Moreover, the imaging collapses the energy such that the single septaria in the layers can be discriminated, since the horizontal Fresnel radius after migration is 2.25 m, as we have seen above. The two septaria at the left-hand side at a depth of 150 m have been recovered, even if the seismic responses of these geobodies cannot be seen in Fig. 4. A gain compensation could in principle recover amplitudes, but for events whose signal-to-noise ratio is high, since this method enhances the noise and the signal. On the other hand, Q compensation enhances the signal compared to the noise. Examples with noise are given at the end of this section.
Synthetic seismogram computed with the exploding-reflector method, corresponding to the model shown in Fig. 3.
Imaging by RTM without (A) and with (B) Q compensation, corresponding to the data shown in Fig. 4.
Fig. 6 shows a model where the medium is inhomogeneous, with a high-loss layer embedded in the Boom Clay. The septaria geobodies are denoted by black-diamond shapes. The simulation parameters are the same as in the previous example. The Ricker wavelet is used as a source, with a peak frequency of 200 Hz, and the time step is 0.08 ms. The synthetic seismograms are displayed in Fig. 7, where it can be observed that the septaria below the layer are attenuated with respect to those above the layer (Fig. 7B). The migration without (a) and with (b) Q compensation is displayed in Fig. 8. It can be seen that the septaria geobodies located below the high-loss layer are not well imaged in Fig. 8A, while they have been recovered using compensated RTM imaging in Figure 8B. Compared to the reference image in Fig. 8C, all the diffractors are clearly delineated.
Septaria bodies embedded in inhomogeneous Boom Clay.
Synthetic seismograms computed with the exploding-reflector method, corresponding to the model shown in Fig. 6. A. Lossless; B. Lossy; C. Lossy with noise (S/N = 10 dB); D. Lossy with noise (S/N = 5 dB).
In order to validate the robustness of the algorithm with random noise, we corrupted the seismograms with strong random noise (see Figs 7C and 7D). The diffractions from septaria bodies are not easy to identify. During compensation, a low-pass Tukey filter with a cut-off frequency of 350 Hz and taper ratio 0.1 is applied. The cut-off frequency has been determined by identifying the noise level in the seismic data spectrum. The imaging results are presented in Figs 9A and 9B. The background appears noisy but both migrated sections clearly show the septaria geobodies.
Conclusions
Conventional seismic processing is designed to restore amplitudes by a time-dependent scaling, but this procedure introduces amplitude and phase errors to the reconstruction of reflectors since it is not based on the geological model. Moreover, this simple scaling cannot recover the correct location of the interfaces, therefore seismic migration with loss compensation is the correct method to recover a true amplitude wavefield. Here, we have presented a seismic migration methodology to compensate for attenuation loss effects in the post-stacked domain. The approach is based on the exploding-reflector concept, where back-propagation is performed by taking a negative quality factor. The viscoacoustic equation used for the propagation is dispersionless, that is, the phase velocity is frequency-independent.
The imaging method was applied to image septaria bodies and thin layers present in the Boom Clay in Belgium. The first example considers a homogeneous medium in terms of velocity but inhomogeneous regarding the attenuation properties. The amplitudes of the septaria embedded in a high-loss zone are fully recovered and an analysis based on the Fresnel-zone concept indicates that the resolution is greatly improved after migration. The second example considers a fully inhomogeneous medium and isolated septaria at random locations. It can be seen that the amplitudes of the septaria geobodies located below a high-loss thick layer have been recovered and the geobodies correctly imaged, even in the presence of different levels of noise.
The results indicate that this approach can effectively improve the resolution and quality of images, particularly at and beneath high-attenuation zones.
Acknowledgments
T. Zhu is grateful to the Jackson Postdoctoral Fellowship at the University of Texas at Austin for his financial support.
Appendix A: Phase velocity and attenuation and quality factors
A plane-wave analysis of equation (3) in the homogeneous case, without the source term, implies [∂
x
(ρ−1∂
x
p) + ∂
z
(ρ−1∂
z
p)] = ρ−1Δp → −ρ−1
k
2
p, ∂
tt
p → −ω2
p and ∂
t
p → iωp, where Δ is the Laplacian, k is the wavenumber, ω is the angular frequency and
${\rm{i}} = \sqrt { - 1} $
. With this substitution, we obtain
$$\begin{equation}
{c^2}{k^2} = {\omega ^2} - {\gamma ^2} - 2i\gamma \omega .
\end{equation}$$
Then, the complex velocity is
$$\begin{equation}
v = \frac{\omega }{k} = \frac{{\omega c}}{{\omega - {\rm{i}}\gamma }}
\end{equation}$$
(Carcione, Reference Carcione2015). From this quantity, the phase velocity and attenuation factor are
$$\begin{equation}
{v_p} = {\left[ {{\rm{Re}}\left( {\frac{1}{v}} \right)} \right]^{ - 1}} = c
\end{equation}$$
and
$$\begin{equation}
\alpha = - \omega {\rm{Im}}\left( {\frac{1}{v}} \right) = \frac{\gamma }{c},
\end{equation}$$
respectively (e.g. Carcione, Reference Carcione2015) Both properties are frequency independent, meaning that there is no velocity dispersion.
On the other hand, defining the quality factor as the total energy density divided by the dissipated energy density (see equations (2.124), (3.131) and (3.132) of Carcione (Reference Carcione2015)), we have
$$\begin{equation}
Q = \frac{{{\rm{R}}{{\rm{e}}^2}(v)}}{{{\rm{Im}}({v^2})}} = - \frac{{{\rm{R}}{{\rm{e}}^2}(k)}}{{{\rm{Im}}({k^2})}} = \frac{\omega }{{2\gamma }}.
\end{equation}$$
The value of γ can be obtained from the quality factor Q 0 at a given reference frequency ω0 = 2πf 0, as
$$\begin{equation}
\gamma = \frac{{{\omega _0}}}{{2{Q_0}}}.
\end{equation}$$
Appendix B: Finite-difference discretisation
Assuming constant density, we discretise the Laplacian. Space and time are discretised as x = ih, z = jh and t = ndt, where h and dt are the respective cell sizes. A suitable fourth-order accurate representation is
$$\begin{eqnarray}
\Delta {p_{i,j}} = \frac{1}{{12{h^2}}}[ - 60{p_{i,j}} + 16({p_{i + 1,j}} + {p_{i,j + 1}} + {p_{i - 1,j}} + {p_{i,j - 1}})\nonumber\\
&& -\, ({p_{i + 2,j}} + {p_{i,j + 2}} + {p_{i - 2,j}} + {p_{i,j - 2}})]
\end{eqnarray}$$
(e.g., Abramowitz & Stegun, Reference Abramowitz and Stegun1964, p. 885), otherwise, the Fourier pseudospectral method can be used (e.g., Baysal et al., Reference Baysal, Kosloff and Sherwood1983; Carcione, Reference Carcione2015).
The time discretisation is
$$\begin{equation}
{\partial _{tt}}p \to \frac{{p_{i,j}^{n + 1} - 2p_{i,j}^n + p_{i,j}^{n - 1}}}{{d{t^2}}},\quad {\partial _t}p \to \frac{{p_{i,j}^{n + 1} - p_{i,j}^{n - 1}}}{{2dt}}.
\end{equation}$$
Forward propagation is performed from equation (3) as
$$\begin{eqnarray}
p_{i,j}^{n + 1} = \frac{1}{{1 + \epsilon }}\big[ (2 - {\epsilon ^2})p_{i,j}^n - (1 - \epsilon )p_{i,j}^{n - 1} \nonumber\\
&&+\, {c^2}d{t^2}\Delta {p_{i,j}} + d{t^2}{s^\prime } \big],
\end{eqnarray}$$
where ε = γdt.
From equation (15), back-propagation is performed as
$$\begin{eqnarray}
p_{i,j}^{n - 1} = \frac{1}{{1 - \epsilon }}\left[ {(2 - {\epsilon ^2})p_{i,j}^n - (1 + \epsilon )p_{i,j}^{n + 1} + {c^2}d{t^2}\Delta {p_{i,j}}} \right].\nonumber\\
&&\end{eqnarray}$$
Note that this equation is obtained from equation (15) if we consider a negative Q 0 (or negative ε), therefore back-propagation with equation (16) implies antidamping and recovery of the signal amplitude.
In the case of the exploding-reflector equation (4), the discretised equations are similar to equations (15) and (16), with the following substitution:
$$\begin{eqnarray}
&&{c^2}d{t^2}\Delta {p_{i,j}} \to cd{t^2}{\Delta _c}{p_{i,j}},\;\;\;{\Delta _c} = {\partial _x}c(x,z){\partial _x} + {\partial _z}c(x,z){\partial _z}p.\nonumber\\
&&\end{eqnarray}$$
In the migration process, the seismogram is a time-dependent boundary condition in equation (16). The time step dt is equal to the sample rate of the data. The seismic trace is applied at each receiver in reverse time and the propagation goes back in time until the origin time, where the best focusing occurs. The reverse modelling sums the energy of all receivers, enhancing the signal-to-noise ratio. The imaging condition is that of Gajewski & Tessmer (Reference Gajewski and Tessmer2005), that is, the origin times of the events are given by the time where maximum focusing (maximum amplitude) occurs.
To prevent the growing of high-frequency noise, we use a low-pass Tukey filter in the wavenumber domain (Tukey, Reference Tukey1967; Zhu, Reference Zhu2014). The cut-off wavenumber is calculated from the cut-off frequency based on the maximum phase velocity of the model. A suitable cut-off frequency is estimated by identifying the noise in the spectrum of the observed data.
Appendix C: The Fresnel zone before and after migration
The Fresnel zone is a measure of the horizontal resolution. Geobodies smaller than the Fresnel zone usually cannot be detected using seismic waves, but migration highly improves the resolution. Since the vertical resolution is λ/4, where λ is the wavelength, the horizontal resolution before migration, represented by the Fresnel radius, is
$$\begin{equation}
R = \sqrt {{{\left( {{z^2} + \frac{\lambda }{4}} \right)}^2} - {z^2}} \approx \sqrt {\frac{{\lambda z}}{2}},
\end{equation}$$
where z is the reflector depth and λ = c/fp (Elmore & Heald, Reference Elmore and Heald1969).
Migration is a downward continuation of the seismic energy from the receivers to the reflectors such that the theoretical limit is obtained for z = 0, that is, after migration the Fresnel radius is
$$\begin{equation}
R = \frac{\lambda }{4}.
\end{equation}$$
Open access
