# Fractional integration of seismic wavelets in anelastic media to recover multiscale properties of impedance discontinuities

Stephan Ker, Yves Le Gonidec

► To cite this version:

Stephan Ker, Yves Le Gonidec. Fractional integration of seismic wavelets in anelastic media to recover multiscale properties of impedance discontinuities. Geophysics, 2018, 83 (2), pp.V61-V71. 10.1190/geo2017-0527.1. insu-02024817

HAL Id: insu-02024817

https://insu.hal.science/insu-02024817v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Fractional integration of seismic wavelets in anelastic media to recover multiscale properties of impedance discontinuities

Stephan Ker<fnref n="1" /> and Yves Le Gonidec<fnref n="2" />

## ABSTRACT

Multiscale seismic attributes based on wavelet transform properties have recently been introduced and successfully applied to identify the geometry of a complex seismic reflector in an elastic medium. We extend this quantitative approach to anelastic media where intrinsic attenuation modifies the seismic attributes and thus requires a specific processing to retrieve them properly. The method assumes an attenuation linearly dependent with the seismic wave frequency and a seismic source wavelet approximated with a Gaussian derivative function (GDF). We highlight a quasi-conservation of the Gaussian character of the wavelet during its propagation. We found that this shape can be accurately modeled by a GDF characterized by a fractional integration and a frequency shift of the seismic source, and we establish the relationship between these wavelet parameters and Q. Based on this seismic wavelet modeling, we design a time-varying shaping filter that enables making constant the shape of the wavelet allowing retrieval of the wavelet transform properties. Introduced with a homogeneous step-like reflector, the method is first applied on a thin-bed reflector and then on a more realistic synthetic data set based on an in situ acoustic impedance sequence and a high-resolution seismic source. The results clearly highlight the efficiency of the method in accurately restoring the multiscale seismic attributes of complex seismic reflectors in anelastic media by the use of broadband seismic sources.

## INTRODUCTION

In the framework of seismic analysis aiming at quantifying physical properties of seismic reflectors, recent methodologies have been

developed using multiscale decomposition of seismic traces (Cas-tagna et al., 2003), in particular, based on multifrequency wavelet analyses (Ker et al., 2013). The interest and scientific issue of such developments have been motivated by a significant increase of the seismic data frequency bandwidth related to recent improvements in seismic acquisition: broadband seismic acquisition (Soubaras and Dowle, 2010) and source design (Marsset et al., 2010) now allow assessing seismic information over several octaves, a prerequisite when imaging complex seismic structures. The potential of multi-frequency quantitative characterization of various kinds of seismic reflectors has been extensively demonstrated, including subsurface thin beds (Ker et al., 2011), gas-charged layers (Ker et al., 2014), and more recently thermocline-related oceanographic structures (Ker et al., 2015). The method is the so-called wavelet response (WR) introduced by Le Gonidec et al. (2002) that considers wavelet propagation as a natural extension of continuous wavelet transform (CWT) signal processing. In the WR method, the seismic sources consist of wavelet signals that are dilated versions of a common analyzing wavelet defined in the theoretical framework of the CWT: This wavelet family is used to probe the medium to perform a multiscale analysis of its acoustic impedance profile. Assuming the validity of the Born approximation, the seismic WR method is equivalent to computing the CWT of the in situ impedance profile (Le Gonidec et al., 2002). In particular, this allowed (1) extending the use of the CWT ridge functions as new multiscale seismic attributes (Ker et al., 2011), (2) providing an original processing sequence to correct the WR from the seismic source signature (Ker et al., 2012), and (3) merging seismic sources in the wavelet domain to improve seismic resolution capability (Ker et al., 2013).

But, the capacity of the WR in characterizing reflector geometry could be altered. For instance, Le Gonidec and Gibert (2006) highlight the impact of viscoelastic properties in the ridge function associated to a complex reflector. A similar effect is observed by Ker et al. (2014) with the WR associated to a half-space gas-charged medium. This suggests that the WR can be strongly affected by

Manuscript received by the Editor 9 August 2017; revised manuscript received 28 September 2017; published ahead of production 11 November 2017; published online 26 December 2017.

<footnote n="1">¹IFREMER, Géosciences Marines, Plouzané, France. E-mail: stephan.ker@ifremer.fr.</footnote>

<footnote n="2">²Université Rennes 1, Géosciences Rennes, Rennes Cedex, France. E-mail: yves.legonidec@univ-rennes1.fr.</footnote>

<footnote>© 2018 Society of Exploration Geophysicists. All rights reserved.</footnote>
seismic attenuation, a physical phenomenon that may affect seismic waves during their propagation and constitutes two main research topics. The first topic deals with attenuation as a quantitative parameter, i.e., the so-called quality factor *Q*, related to the physical properties of the medium, for instance, in seismology and in exploration seismology, in which attenuation can be related to fault density (Worthington and Hudson, 2000) or gas saturation (Morgan et al., 2012). The second topic deals with the impact of attenuation on the resolution of seismic images, i.e., so-called *Q* inverse filtering, which aims at correcting this effect to improve imaging of oil and gas reservoirs (Wang, 2008). To better address these issues, seismic attenuation mechanisms have to be identified, with a distinction between extrinsic and intrinsic phenomena. The former refers to attenuation generated by geometric spreading of the seismic wave and to scattering effects related to the presence of heterogeneities in the medium. The latter is related to the viscoelastic properties of the medium, in which the elastic energy of a seismic wave is gradually converted into heat and then dissipated (Toksöz and Johnston, 1981). This attenuation is inversely proportional to the quality factor *Q* defined as the wave energy divided by the energy loss per cycle. Seismic attenuation is a frequency-dependent phenomenon, and its effect during wave propagation is to modify the amplitude and the shape of the seismic waveform, related to dispersion (frequency-dependent wave velocity). As a consequence, anelasticity may make the WR different from the CWT of the acoustic impedance properties, preventing the use of ridge functions as accurate attributes to assess the geometry of seismic reflectors. The present paper deals with extended developments of the WR method to the case of anelastic media for which we establish an original framework to analyze seismic reflectors by the use of broadband seismic sources.

After a brief overview of the WR method established for elastic media, we introduce the modification of the WR when dealing with anelastic media characterized by an attenuation factor and a dispersive phase velocity. We consider seismic source wavelets defined as Gaussian derivative functions (GDFs), commonly used in seismic analysis. We show that the Gaussian character of such sources is preserved by integration during propagation in a Kolsky-Futterman (K-F) attenuating medium (Kolsky, 1956; Futterman, 1962). A description of a K-F medium, characterized by a linear dependency of the attenuation with the seismic frequency, is provided in this paper to briefly recall the key properties of such an anelastic medium. By introducing a dimensionless parameter, it is possible to establish the expression of the quality factor *Q* of the anelastic medium as a function of the GDF properties. Then, we consider a seismic propagation medium characterized by a *Q* factor constant with depth and an acoustic reflectivity sequence defined by either step-like or thin-bed discontinuities, and we illustrate how attenuation affects the WR analysis. Finally, we take advantage of the quasi-conservation of the Gaussian character of a GDF to design a time-varying shaping filter, an original approach used to compensate the WR for attenuation effects. We illustrate the approach on a synthetic data set involving an in situ acoustic reflectivity sequence embedded in an anelastic medium and a broadband seismic source signature.

## WAVELET RESPONSE OF AN ANELASTIC MEDIUM

The WR is a multiscale seismic probing of a medium based on the properties of the wavelet transform (Le Gonidec et al., 2002). Instead

of using a single seismic source, the WR method uses a family of source wavelets that are dilated or contracted versions of the same signal called the mother analyzing wavelet $\psi$. The WR is thus a collection of seismic traces obtained by a set of seismic sources forming a wavelet family that enables a time-scale decomposition of the acoustic impedance $p$ of the medium. The expression of the WR related to a plane-wave propagation in 1D medium is given by

$$R[\psi, p](t, a) = \mathcal{D}_a \psi(t) * r(t), \quad (1)$$

where $*$ stands for the convolution operator, $r(t)$ is the reflectivity sequence in the time domain, and $\mathcal{D}_a$ is the dilation operator controlled by the dilation $a$ that enables to contract or dilate the mother wavelet $\psi(t)$ according to $\mathcal{D}_a \psi(t) = (1/a) \psi(t/a)$. When dealing with one element of this wavelet family, we will consider the reference wavelet $\psi(t)$ for the sake of simplicity. The function $R$ is the wavelet transform of the reflectivity $r(t)$, and Le Gonidec et al. (2002) show that it is strictly equivalent to the wavelet transform of $p$ performed in the traveltime domain $t$ when the Born approximation is valid. The mother wavelet $\psi(t)$ has to satisfy some necessary conditions, i.e., an oscillating differentiable function with compact support and zero-mean value (Mallat, 1998). In the present work, $\psi(t)$ is a GDF, such as the Ricker wavelet commonly used in seismic imaging and discussed in the next section. The potential of GDF properties has been used by Ker et al. (2012) to design a multiscale processing of high-resolution seismic data, and it is now extended to study the propagation of multiscale seismic waves in anelastic media.

In an elastic homogeneous medium, the shape of a wavelet source signal $\psi(t)$ generated by a seismic emitter remains unchanged and the wavelet conserves the same features during its propagation at depth $z$, where the seismic wave is $\phi_0(z, t) \equiv \psi(t)$. When the medium is anelastic, the seismic wave experiences some changes, including amplitude decrease, dispersion, and phase rotation. As a consequence, the seismic wavelet becomes a time-depth-varying wave defined as

$$\phi(z, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \hat{\psi}(\omega) \hat{u}(z, \omega, t) d\omega, \quad (2)$$

where $\hat{\psi}(\omega)$ is the Fourier transform of $\psi(t)$ and $\omega$ is the angular frequency. The term $\hat{u}(z, \omega, t)$ is the Fourier transform of the displacement or the pressure (i.e., the waveform) which can be, in the framework of plane-wave propagation in 1D anelastic media, expressed as

$$\hat{u}(z, \omega, t) = \exp(-\beta z) \exp \left[ i\omega \left( t - \frac{z}{\gamma} \right) \right], \quad (3)$$

where $\beta$ is the attenuation of the wave amplitude and $\gamma$ is the phase velocity, which describes the velocity dispersion and the phase rotation of the wave. The WR is thus modified and, noted $\tilde{R}$, it can be described as a nonstationary process (Margrave, 1998), which can be expressed, based on the superposition principle, as

$$\tilde{R}[\psi, p](t, a) = \int_{0}^{\infty} \mathcal{D}_a \phi(z, t) r(z) dz. \quad (4)$$

In the following, $\phi(z, t)$ is noted $\phi_\tau(t)$, where $\tau$ stands for the time of propagation in the anelastic medium, and we are reminded that $\phi_0(t)$ corresponds to the elastic case.
As a result, attenuation effects preclude the use of the properties of the WR to study acoustic impedance discontinuities from multi-scale seismic attributes. The aim of the present paper is to develop a method to recover reliable properties of multiscale seismic attributes allowing quantitative description of the geometry of reflectors. The approach is based on few assumptions described below.

# QUASI-CONSERVATION OF THE GAUSSIAN CHARACTER OF A SEISMIC WAVELET IN ANELASTIC MEDIA

## Kolsky-Futterman anelastic medium

A first assumption deals with the modeling of the anelastic propagation medium. In numerous attenuation models, the attenuation parameter $\beta$ introduced in equation 3 may be expressed as a power law of the frequency $\beta \propto \omega^\delta$ (Wang, 2008). In particular, $\delta = 1$ corresponds to the Kolsky-Futterman model, which describes a linear frequency dependency (Kolsky, 1956; Futterman, 1962). This model is intensively used for attenuation analysis (Xintao et al., 2014) or designing an inverse $Q$ filter (Wang, 2008) because it has a high degree of similarity to several $Q$ models (Ursin and To-verud, 2002; Wang, 2008). In the K-F model, the attenuation $\beta$ and phase velocity $\gamma$ are defined by

$$ \beta = \frac{\omega}{2cQ}, \qquad (5a) $$

$$ \frac{1}{\gamma} = \frac{1}{c} \left( 1 - \frac{1}{\pi Q} \ln \frac{\omega}{\omega_h} \right), \qquad (5b) $$

where $c$ is the reference phase velocity and $\omega_h$ is a tuning parameter corresponding to the highest possible seismic frequency and ensuring compliance with the Kramers-Kröning dispersion relation (Wang and Guo, 2004).

When $Q$ tends to infinity, the attenuation $\beta$ tends to zero and the phase velocity $\gamma$ tends to the reference velocity; i.e., the medium is equivalent to an elastic medium: This is the reference case used in the present work. In the following, a value of $Q = 50$ is selected for numerical simulation purposes.

## Gaussian derivative seismic source signal

A second assumption is to assimilate a seismic source signal with a wavelet defined as a GDF, i.e., a derivative of order $\alpha$ of a Gaussian function (Ker et al., 2012). This assumption is commonly used in seismic analysis: The first ($\alpha = 1$) and second ($\alpha = 2$) derivatives of a Gaussian function correspond to the wavelet erroneously called minimum-phase Ricker (Hosken, 1988) used in seismic modeling, for instance, and to the so-called Ricker wavelet (Ricker, 1953), respectively. This second wavelet is a symmetric wavelet widely used to model seismic source signals in seismic imaging, attenuation estimation, and $Q$ inverse filtering (Wang, 2008, 2015a). Higher derivative orders $\alpha$ are used in acoustic logging to model sources with a high number of cycles (Heigl, 2007); we note that fractional derivative orders, which correspond to asymmetric GDF, have been recently used to better represent seismic signals, such as vertical seismic profile (VSP) data waveforms (Wang, 2015b) or to process large frequency bandwidth seismic data (Ker et al., 2013).

The general expression of a reference wavelet $\psi$ defined by a GDF to model a seismic source signal in the time domain is given by

$$ \psi(t) = \Omega \frac{d^\alpha}{dt^\alpha} \exp \left( -\frac{\omega_0^2 t^2}{4} \right). \qquad (6) $$

The parameter $\Omega$ is the amplitude, the derivative order $\alpha$ can either be an integer or be fractional, and the angular frequency $\omega_0$ is the reciprocal of the reference dilation defined by $a_0 = 2/\omega_0$. Note that a GDF can be written as the product of a Hermite polynomial and a Gaussian potential (Heigl, 2007; Ker et al., 2012); i.e., the amplitude spectrum $\hat{\psi}$ of the GDF $\psi$ is always derived from the Gaussian function and is expressed in the frequency domain by (Wang, 2015b)

$$ \hat{\psi}(\omega) = \Omega \frac{\omega^\alpha}{\omega_0^\alpha} \exp \left( -\frac{\omega^2}{\omega_0^2} \right). \qquad (7) $$

In seismic analysis, a seismic wavelet is commonly characterized by its peak frequency $f_p$, which is the frequency associated with the maximum of the amplitude spectrum; i.e., $f_p = \omega_p/2\pi$, where the peak angular frequency $\omega_p$ is defined by $(\partial\hat{\psi}/\partial\omega)|_{\omega_p} = 0$. For a GDF wavelet, $\omega_p = \omega_0\sqrt{\alpha/2}$ (Ker et al., 2012; Wang, 2015b), which means that $f_p$ depends on the derivative order $\alpha$ and the natural angular frequency $\omega_0$. In the following, $\alpha$ and $f_p$ will be the quantified wavelet parameters.

## Attenuated wavelet modeled by a GDF

In this section, we consider a Ricker seismic source, a GDF wavelet $\psi$ defined by a derivative order $\alpha = 2$ and a peak frequency $f_p = 180$ Hz typical of a small airgun source (Piété et al., 2013). We model the attenuated wavelet signal $\phi_\tau(t)$ according to equation 2, based on the K-F model with $Q = 50$ and measured at different times of propagation $\tau = [0; 25; 50; 75; 100]$ ms, where $\tau = 0$ s corresponds to $\phi_0(t) = \psi(t)$. The results, plotted in solid red lines in the time domain (Figure 1a1–1a5) and in the frequency domain (Figure 1b1–1b5), clearly highlight the dispersive character of the seismic wave during its propagation in an anelastic medium; the wavelet symmetry is broken and the frequency content is shifted to lower frequencies, which means that there is no conservation of $\alpha$ and $f_p$ when a seismic signal propagates in an anelastic medium. Note that recently, Wang (2015b) observes a similar behavior on VSP data and suggests attenuation as a plausible mechanism, without exploring the possible relationship with the wavelet-shape parameters.

The decrease of $\alpha$ and $f_p$ with $\tau$ suggests that the attenuated wavelet signal $\phi_\tau(t)$ may be an integrated version of the wavelet $\psi(t)$: This means that we may approximate $\phi_\tau(t)$ with an equivalent GDF wavelet $\phi'_\tau(t)$ defined by a derivative order $\alpha'$, a natural angular frequency $\omega'_0$ associated with a peak frequency $f'_p$, and an amplitude $\Omega'$. This approach is expressed by assuming

$$ \phi_{\tau}(t) \sim \phi'_{\tau}(t) = \Omega' \frac{d^{\alpha'}}{dt^{\alpha'}} \exp \left( -\frac{1}{4} \omega_0^{2/3} t^2 \right). \qquad (8) $$

To evaluate this assumption, we apply the method developed by Wang (2015b) to determine the $\phi'_\tau(t)$ parameters from the centroid frequency $\omega_m$ and the standard deviation $\omega_\sigma$ evaluated from the power spectrum of $\phi_\tau(t)$
$$ \left(\frac{1}{2\alpha'} + 1\right) \left(\frac{\Gamma(\alpha' + 1/2)}{\sqrt{\alpha'}\Gamma(\alpha')}\right)^2 = \frac{\omega_{\sigma}^2}{\omega_m^2}, \qquad (9a) $$

$$ \omega_0' = 2\sqrt{\frac{\omega_m^2 + \omega_{\sigma}^2}{1 + 2\alpha'}}, \qquad (9b) $$

where $\Gamma$ is the gamma function (Abramowitz and Stegun, 1972). For each $\phi_{\tau}(t)$ seismic wavelet plotted in solid red lines in Figure 1,

<figure><img src="image_3.png" /><figcaption>Figure 1. Normalized attenuated seismic wavelets $\phi_{\tau}(t)$, associated to a Ricker seismic source (GDF derivative order $\alpha = 2$ and peak frequency $f_p = 180$ Hz), based on the Kolsky-Futterman model with $Q = 50$ and for $\tau = [0; 25; 50; 75; 100]$ ms, respectively (solid red lines) and associated GDF models $\phi'_{\tau}(t)$ (dashed black line): representation in the time domain (a1-a5) and in the frequency domain (b1-b5).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. (a1-a5) Normalized GDF wavelet sources $\psi(t)$ associated to derivative orders $\alpha = [1; 2; 3; 4; 5]$ and peak frequency $f_p = 180$ Hz, respectively. (b1-b5) Normalized attenuated seismic wavelets $\phi_{\tau}(t)$ based on the Kolsky-Futterman model with $Q = 50$ and $\tau = 100$ ms (solid red curves) and associated GDF models $\phi'_{\tau}(t)$ (dashed black curves) in the time domain, and (c1-c5) in the frequency domain, in which the normalized solid black curve correspond to the nonattenuated wavelet $\phi_0(t)$ (equivalent to $\psi(t)$). Derivative order $\alpha'$ and peak frequency $f'_p$ stand for the GDF wavelet model.</figcaption></figure>

we determine the associated GDF approximation $\phi'_{\tau}(t)$ according to this method: The results, plotted in dashed black lines in Figure 1, show a very good agreement between $\phi_{\tau}(t)$ and $\phi'_{\tau}(t)$.

Because the derivative order $\alpha$ is a key parameter to model the shape of seismic source signatures (Heigl, 2007; Wang, 2015b) and to analyze complex discontinuities (Le Gonidec et al., 2002), it is of critical importance to evaluate the accuracy of the GDF approximation as a function of $\alpha$. We study how different GDF wavelet sources $\psi(t)$ with a similar peak frequency $f_p = 180$ Hz, but different integer derivative orders $\alpha = [1; 2; 3; 4; 5]$ (Figure 2a1-2a5, solid black curves) are modified after a propagation of $\tau = 100$ ms, corresponding to a depth of 150 m for a sound velocity of 1500 m/s. We model the nonattenuated wavelet $\phi_0(t)$, which is equivalent to $\psi(t)$ and used as a reference; the attenuated wavelet signal $\phi_{\tau}(t)$ according to equation 2; and the approximated attenuated wavelet signal $\phi'_{\tau}(t)$ defined as a GDF with the parameters given in equation 9a and 9b. The results are plotted in the time domain in Figure 2b1-2b5 and in the frequency domain in Figure 2c1-2c5: The normalized solid black, dashed black, and solid red curves correspond to $\phi_0(t)$, $\phi_{\tau}(t)$, and $\phi'_{\tau}(t)$, respectively.

As a key result, we show that a GDF seismic wavelet affected by a K-F attenuation model is an attenuated wavelet signal, which can be modeled by a specifically designed GDF $\phi'_{\tau}(t)$. The quasi-conservation of the Gaussian character means that attenuation effects can be approximated by a filter $F_{\tau}(t)$, which has Lévy alpha-stable properties (Ker et al., 2012) according to

$$ \phi'_{\tau}(t) = \psi(t) * F_{\tau}(t), \qquad (10) $$

where the wavelet source $\psi(t)$ depends on the parameters $(\Omega, \alpha, \omega_0)$ according to equation 6 and the filter $F_{\tau}(t)$ depends on the anelastic medium quality factor $Q$ and is function of the time of propagation $\tau$. The attenuated GDF wavelet $\phi'_{\tau}(t)$ depends on the parameters $(\Omega', \alpha', \omega_0')$, and the aim of the next section is to establish the explicit relationship between the physical parameter $Q$, the propagation time $\tau$, and the wavelet parameters $\omega_0, \alpha, \omega_0'$, and $\alpha'$.

## QUANTITATIVE RELATIONSHIP BETWEEN GDF PARAMETERS AND ATTENUATION

### Attenuation-dependent wavelet parameters

The WR method involves a wavelet family $D_a\psi(t)$ for which each member is defined by the initial set of parameters $(\Omega, \alpha, \omega_0)$, and it is modified during propagation within an anelastic medium characterized by a quality factor $Q$. To characterize how each member of the wavelet family is modified, we perform a quantitative
analysis of the variation of the set of parameters (Ω', α', ω'₀) associated with the GDF that models the attenuated wavelet signal. We assess the range of these attenuation-related parameters for a seismic propagation in an anelastic medium taking into account the wavelet family source signals. To do so, we introduce the dimensionless parameter Υ = τf_pQ⁻¹. For a range of Υ as wide as [10⁻³ – 30] describing nonattenuating to highly attenuated seismic signals, we determine (Ω', α', ω'₀) according to equation 9a and 9b for the five GDF analyzing wavelets introduced in Figure 2, i.e., for α = [1; 2; 3; 4; 5]. The variations of the amplitude ratio (Figure 3a, in dB) are weak for Υ < 10⁻¹, with a decreasing trend when Υ increases: For Υ > 1, the amplitude drops are more pronounced for high derivative orders α. The derivative order α' (Figure 3b) exhibits a hysteresis behavior: At low Υ values, which corresponds in particular to short time of propagation and/or low attenuation, the attenuation effects can be neglected and α' ~ α as expected. At high Υ values, α' tends toward much lower and fractional orders: For α = [1; 2; 3; 4; 5], we obtain α' ~ [0.2; 0.8; 1.3; 1.9; 2.4], i.e., α' ~ α/2 as a rough estimate, which shows the significant integrating effects. The ratio ω'₀/ω₀ (Figure 3c) also exhibits a hysteresis behavior: The curve starts with a value close to one at low Υ values, followed by a steep decrease and trend toward zero at high Υ values.

For a practical applicability of the method, the asymptotic behavior of the derivative order and the frequency shift of a seismic wave that propagates in a highly attenuating medium (Υ ≫ 0.1) is sensitive to noise conditions, as suggested by the strong decrease of the amplitude (Figure 3a). As a consequence, even if signal-to-noise ratio (S/N) issues affect not only this method but also many conventional geophysical processing methods, it may be unlikely to observe this asymptotic behavior if the S/N is not high enough.

## General Q expression from fractional integrated wavelet parameters

We now demonstrate how the quality factor can be evaluated from the parameters of the attenuated wavelet signal $\phi_{\tau}(t)$. The associated amplitude spectrum $\hat{\phi}_{\tau}(\omega)$ is first defined following equations 2 and 7 for a time of propagation $\tau$, which gives, based on a K-F attenuating medium:

$$ \hat{\phi}_{\tau}(\omega) = \Omega \frac{\omega^{\alpha}}{\omega_{0}^{\alpha}} \exp\left(-\frac{\omega^2}{\omega_{0}^2}\right) \exp\left(-\frac{\tau\omega}{2Q}\right). \quad (11) $$

By definition, the peak angular frequency $\omega_{\tau}$ of the attenuated seismic wavelet $\phi_{\tau}(t)$ can be obtained by solving $(\partial\hat{\phi}_{\tau}/\partial\omega) = 0$. As a result, the expression of Q is given by

$$ Q = \frac{\tau\omega_{\tau}\omega_{0}^{2}}{2(\alpha\omega_{0}^{2} - 2\omega_{\tau}^{2})}, \quad (12) $$

which is a generalization to any fractional GDF wavelet source of the result given by Zhang and Ulrych (2002) for a Ricker wavelet. This highlights the usefulness of approximating a seismic signal, the amplitude spectrum of which may be defined by several solutions of $(\partial\hat{\phi}_{\tau}/\partial\omega) = 0$, by a GDF wavelet characterized by a unique solution $\omega_{\tau}$. As shown in the previous section, the approach based on the centroid frequency and standard deviation of the seismic signal $\phi_{\tau}(t)$ makes this assessment more robust compared with peak-frequency detection (Tary et al., 2017). According to the

approximation of $\phi_{\tau}$ by a GDF wavelet $\phi'_{\tau}$, $\omega_{\tau}$ can be assimilated to $\omega'_{p} = \omega'_{0}\sqrt{\alpha'/2}$, resulting in the following approximation:

$$ Q \sim Q' = \frac{\tau\omega'_{0}\sqrt{\frac{\alpha'}{2}\omega_{0}^{2}}}{2(\alpha\omega_{0}^{2} - \alpha'\omega_{0}'^{2})}. \quad (13) $$

Interestingly, this can be rewritten $Q' = (\tau\pi f'_{p} f'_{p})/\alpha(f'_{p} - f'_{p})$, which means that the quality factor can actually be assessed by only determining the peak frequency $f'_{p}$ of the attenuated seismic wavelet.

The accuracy of the quality factor estimated from equation 13 is now assessed by modeling the attenuation of seismic wavelets previously presented for typical Q values in the range [5–250] and a seismic source wavelet peak frequency $f_p$ in the range [50–1000 Hz] to consider the propagation of low to very high resolution seismic sources in strongly to weakly attenuating media. In this $f_p$-Q domain, we determine the error $|Q - Q'|/Q$ (in %) for five derivative orders $\alpha = [1; 2; 3; 4; 5]$ of the GDF source wavelet (Figure 4). As a first observation, the distribution of errors shifts toward the minimum error value as $\alpha$ increases and on average, the error decreases with $\alpha$: Its highest value is larger than 11% for $\alpha = 1$, approximately 8.5% for $\alpha = 2$ (Ricker wavelet) and less than 4% for $\alpha = 5$. The decreasing tendency of the mean error with $\alpha$ is not linear, and a preliminary estimate suggests an exponential behavior. At low frequencies ($f_p < 100$ Hz), the error slightly increases with Q but remains very small; i.e., estimating Q from the GDF parameters is accurate at low frequencies for any Q values. Note that for $\alpha = 1$, the error is minimum for $Q < 50$, when it is maximum for $\alpha > 1$. However, a general trend shows that the error increases with $f_p$ and decreases with Q, remaining in a range lim-

<figure><img src="image_5.png" /><figcaption>Figure 3. Variations of the attenuated seismic wavelet parameters as a function of the adimensional parameter Υ for source wavelets defined by $\alpha = [1; 2; 3; 4; 5]$, respectively: (a) amplitude $\Omega'$ normalized by the reference source amplitude $\Omega_0$, (b) derivative order $\alpha'$, and (c) peak frequency $f'_p$ normalized by the reference source peak frequency $f_p$.</figcaption></figure>
ited to approximately 2% independently of $\alpha$. This means that the error is mainly controlled by $\alpha$ and $f_p$, and by $Q$ to a lesser extent. These results highlight the potential of the present approach in quantifying the quality factor from a GDF approximation: The efficiency is particularly improved by considering a seismic source wavelet defined by a high derivative order.

# SEISMIC ATTENUATION AFFECTING MULTISCALE REFLECTORS: QUANTITATIVE INTERPRETATION BASED ON THE WAVELET RESPONSE

## Preliminary considerations

We now consider seismic reflection propagation in an anelastic medium characterized by a homogeneous quality factor $Q$, using the K-F model, and an acoustic impedance profile. In a first approximation, we consider seismic reflections induced by density contrasts.

The analyzing wavelet $\psi(t)$ is an antisymmetric GDF of order $\alpha = 5$ composed of six extrema: As introduced in previous papers, this wavelet has been used to model high-resolution seismic sources (Ker et al., 2013). The wavelet family is composed of 48 elements, the peak frequency of which $f_p$ ranges between 100 and 1000 Hz. The WR is thus a collection of seismic wavelets of different peak frequencies, i.e., different dilations, acquired in normal incidence: The WR corresponds to a time-dilation representation of the wavelet amplitude. Following Ker et al. (2013), we consider the dilation

<figure><img src="image_6.png" /><figcaption>Figure 4. Error on the evaluation of the quality factor (equation 13) for anelastic media characterized by $Q$ in the range [5–250] and for seismic peak frequencies in the range [50–1000 Hz], for different shapes of the seismic source defined by a derivative order $\alpha$.</figcaption></figure>

defined by the wavelet-dominant period (Gesret et al., 2010), called the dominant dilation $a_d = \sqrt{\alpha/2}(\delta/\pi f_p)$, where $\delta \sim 1.75$ for $\alpha = 5$ (Ker et al., 2013). A ridge function is a multiscale seismic attribute, which represents the evolution of an extremum amplitude as a function of the dilation plotted in a log-log diagram: It is a sparse support of the wavelet decomposition of the seismic reflectivity, and that is why it has been proposed in seismic analysis as a multiscale seismic attribute (Ker et al., 2011).

To assess the effect of attenuation on multiscale seismic attributes, we compute $R$ (equation 1) and $\tilde{R}$ (equation 4 and associated equations 2 and 3), respectively, by setting $Q = \infty$ (elastic case) and 50 (anelastic case), respectively. According to the previous sections, we then approximate $\tilde{R}$ by $R'$ for which each element is modeled by a GDF wavelet. Two seismic reflectors are considered below as preliminary case studies, i.e., a step-like reflector and a thin bed.

## Step-like density contrast in an anelastic medium

The acoustic impedance profile of the medium is described by a step-like density contrast located at the time of 37.5 ms, i.e., a homogeneous Heaviside-like reflector geometry Figure 5a. For an elastic medium, $R$ is the WR of a homogeneous reflector already described by Le Gonidec et al. (2002). The WR is equivalent to the wavelet family, scaled by the reflection coefficient and normalized in Figure 5b, which points toward the position of the discontinuity. In that case, the number of ridge functions is equal to the number of extrema of the analyzing wavelet: As illustrated with the third extremum, the ridge function (Figure 5b, blue line) of the Heaviside discontinuity is a flat straight line (Figure 5d, blue line).

For an anelastic medium, the global structure of the WR $\tilde{R}$ looks at a first glance barely changed, still being a cone-like structure pointing toward the discontinuity (Figure 5c), but the third ridge function (Figure 5c, red curve) now differs from a simple straight line (Figure 5d, red curve). In addition, the last ridge function tends to disappear at small dilations, as highlighted when plotting the derivative order $\alpha'$ as a function of the dilation $a_d$ (Figure 5e, red line). At large dilations, $\alpha'$ tends to $\alpha = 5$, which characterizes the elastic case (Figure 5e, blue line), and decreases to four at small dilations, which means a symmetric wavelet with five extrema only. The dilation $a_0'$ is very similar to the reference dilation $a_0$ of the elastic case (Figure 5f, red and blue lines, respectively), with a weak divergence when $a_d$ decreases. Based on equation 13 expressed in the wavelet-transform formalism, we use the wavelet parameters extracted from $R'$ (Figure 5e and 5f) to assess $Q$, which gives a mean value of $51.7 \pm 0.5$ in very good agreement with the expected value of 50. It has to be noted that the error of 3.4% in $Q$ determination is in perfect agreement with results obtained in Figure 4 for $\alpha = 5$.

This simple case study illustrates the method to assess the quality factor of an anelastic medium by the use of the WR method, but also the misinterpretation of a seismic reflector geometry if only one source signal is considered, and reciprocally. This is even worse for complex seismic reflectors for which characteristic size quantification is not straightforward when it is embedded in an attenuating medium, as discussed in the following section.

## Thin bed in an anelastic medium

The geometry of a thin bed is a window function defined as two successive step-like discontinuities: The distance between them
corresponds to the window aperture, which is the characteristic size of the thin bed (Figure 6a). In the elastic case, the structure of the WR *R* is symmetric (Figure 6b) (for a full description, see Le Goinidec et al., 2002). At large dilations compared with the window aperture, the WR has a cone-like structure pointing toward the center of the discontinuity and is composed of five ridge functions: Note that the second and fourth ridge functions (Figure 6b, the solid and dashed blue lines, respectively) are similar in a log-log diagram (Figure 6d) and are characterized by a slope -1 at large dilations, where the window function is equivalent to a Dirac-like reflector. At small dilations, the WR is composed of two subconical structures: Each one points toward a step-like discontinuity composing the window function and is characterized by six ridge functions as discussed in the previous section. At intermediate dilations, the WR is characterized by a coalescence structure with an amplitude extremum related to seismic wave interferences between the edges of the internal structure of the reflector: The associated characteristic dilation is related to the characteristic size of the reflector.

In the anelastic case, the global symmetry of the WR is lost (Figure 6c), as shown by the two selected ridge functions (the solid

and dashed red lines, respectively) which strongly differ from each other when plotted in a log-log diagram (Figure 6e) and cannot be used to properly characterize the reflector geometry. Actually, the slope of the ridge function analysis is not -1 at large dilations, which would mean a reflector not equivalent to a Dirac-like discontinuity as expected. At small dilations, the slope is not zero, which would not allow us to conclude to internal step-like structures of the complex reflector. And more critically, the characteristic dimension of the reflector, which should be identified at intermediate dilations, cannot be assessed with confidence.

These analyses clearly highlight that a quantitative characterization of a seismic reflector embedded in an anelastic medium is not straightforward when the reflector has a complex structure. Reciprocally, assessing the quality factor of a complex seismic reflector is not straightforward because of the multifrequency dependency of the seismic wavelets reflected by the multiscale discontinuities of the propagating medium, including layering effects and interference patterns that may induce apparent attenuation and frequency tuning effects (Tary et al., 2017).

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) The WR associated to a step-like reflector in an (b) elastic and (c) anelastic medium: The blue and red lines, respectively, stand for the ridge functions plotted in (d) as a function of the dominant dilation $a_d$. The variations of the derivative orders $\alpha$ and $\alpha'$, and the variations of the dilations $a_0$ and $a'_0$ of the WR wavelets as a function of $a_d$ are plotted in (e and f), respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. (a) The WR associated to a thin-bed reflector in an (b) elastic and (c) anelastic medium: the blue and red lines, respectively, stand for ridge functions plotted in (d and e), respectively, as a function of the dominant dilation $a_d$. Note that for the elastic case (d), the two ridge functions (solid and dashed lines) are identical, but they differ for the anelastic case (e).</figcaption></figure>
# TIME-VARYING SHAPE FILTERING TO RECOVER MULTISCALE SEISMIC ATTRIBUTES

## Principles of the method

When embedded in an elastic medium, a seismic reflector can be identified in terms of morphology and characteristic dimensions by quantitative analyses of its multiscale seismic attributes based on the WR R properties. This can no longer be achieved when the medium is anelastic, as illustrated in the previous section in which the attenuated WR $\tilde{R}$ is obtained with $Q = 50$: Actually, each time-varying seismic wavelet related to $\tilde{R}$ can be accurately approximated by a GDF wavelet $\phi'_{\tau}(t)$ with a particular shape (defined by $\alpha'$), which deviates more from the shape of the source wavelet $\psi(t)$ the greater $\omega'_0$ and/or $\tau$ are (see Figure 3). We introduce a new method to compensate for this deviation to restore the properties of the WR analysis: The aim is to define a compensated WR $\tilde{R}$ for which each element $\tilde{\phi}_{\tau}(t)$ is a dilated version of a common GDF wavelet $\tilde{\psi}(t)$ defined by a derivation order $\tilde{\alpha}$ and a natural frequency $\tilde{\omega}_0$, associated to a peak frequency $\tilde{f}_p$ as previously discussed. The approach consists of a time-varying filtering method that enables the transformation of a time-varying seismic wavelet $\phi'_{\tau}(t)$ into a wavelet $\tilde{\phi}_{\tau}(t)$ of constant shape (given by $\tilde{\alpha}$), taking advantage of the mathematical properties of the GDF. The choice of the reference constant shape is based on the attenuated wavelet most affected by seismic attenuation, i.e., which corresponds to the smallest dilation of the wavelet source family and the largest propagation time $\tilde{\tau}$. Note that $\tilde{\tau}$ controls not only the derivative order but also the natural angular frequency $\tilde{\omega}_0$, i.e., the peak frequency $\tilde{f}_p$ of the wavelet $\tilde{\phi}_{\tau}(t)$. With this method, we do not aim at recovering the frequency content lost or reduced by seismic attenuation, such as Q inverse filtering or other resolution enhancement techniques would do (Wang, 2008; van der Baan, 2012), but rather at obtaining a constant shaping of wavelets, regardless of the global frequency content reduction of the seismic wavelets. In practice, considering $\tilde{\alpha}$ as an integer value facilitates the ridge functions analysis, but this is not mandatory.

We now define the so-called time-varying shape-filtering method. For each dilation and propagation time, the attenuated seismic wavelet has to be transformed into a dilated version of the analyzing wavelet $\tilde{\psi}$ involving (1) an amplitude correction term $\kappa$; (2) a fractional integration of order $\eta$ to remove the variations of the seismic wave derivative order; (3) a Gaussian term, where $\omega_s$ compensates for the frequency shift of the seismic wave; and (4) an exponential term, where $\tau_s$ compensates for the time shift. In the Fourier domain, the expression of the shaping filter is thus given by

$$ \hat{h}(\omega) = \kappa(i\omega)^{-\eta} \exp\left(-\frac{\omega^2}{\omega_s^2}\right) \exp(-i\omega\tau_s), \quad (14) $$

with

$$ \eta = \alpha' - \bar{\alpha}, \quad (15a) $$

$$ \omega_s = \sqrt{\bar{\omega}_0^2 - \omega_0'^2}. \quad (15b) $$

We perform the shape filtering by using a nonstationary convolution approach: To do so, we define the impulse response of the filter by using the definition of the fractional derivation of a Gaussian function

in the time domain (Caputo, 1967) and the Laplace convolution of causal functions Wang (2015b). According to this approach, we can express the impulse response in the time domain by

$$ h(t) = \begin{cases} 0 & \text{if } t < 0, \\ \kappa \frac{f^{\tau-1}}{\Gamma(\eta)} * \exp\left(-\frac{1}{4}\omega_s^2(t - \tau_s)^2\right) & \text{if } t \ge 0. \end{cases} \quad (16) $$

The numerical implementation of the shape filtering relies on the multiplication of (1) a matrix, in which each row contains the time-varying filter computed for each traveltime sample, by (2) a row vector containing the seismic trace.

In practice, the parameters $\tilde{\alpha}$ and $\tilde{f}_p$ of the constant-shape wavelet are defined for the attenuated seismic wavelet associated with the smallest dilation for a reference propagation time $\tilde{\tau}$ related to the depth position of the seismic reflector. This is done by performing on the attenuated seismic wavelets the fractional integration operation expressed in equation 15a. The parameters of the shaping filter are obtained from the relationships governing the wavelet parameters variation and the relationship describing the attenuation-induced time-shift variations.

## Application to the step-like and thin-bed seismic reflectors embedded in an anelastic medium

The WR analysis of the step-like reflector has been developed above (Figure 5b): The WR of the anelastic case is characterized by a frequency-dependent $\alpha'$, which tends to four at small dilations (Figure 5d, red line). According to the principles of the method described above, we determine the time-varying shape-filtered WR $\tilde{R}$ by considering an analyzing wavelet $\tilde{\psi}$ defined by $\tilde{\alpha} = 4$. Note that this arbitrary integer value corresponds to the less dilated wavelet ($\tilde{f}_p \sim 550$ Hz) for a reference time $\tau = 41.5$ms larger than the two-way traveltime of the seismic reflector. As a result, the method allows the recovery of the symmetric cone-like structure of the WR for the anelastic case, which is now defined by dilated wavelets with a constant-shape $\tilde{\alpha}$, in very good agreement with the WR performed with the same analyzing wavelet in the elastic case (Figure 7a, black and green curves, respectively): We then retrieve the correct morphology of the seismic reflector from the ridge function (Figure 7a, magenta line), which is now properly defined by a flat, straight line (Figure 7c, magenta line).

The efficiency of the time-varying shape filter is clearly demonstrated when assessing the characteristic dimension of a complex reflector in an anelastic medium, such as the thin bed. Similarly to the previous analysis performed on a step-like reflector, we determine the same time-varying filtering with $\tilde{\alpha} = 4$. The compensated WR $\tilde{R}$ has a global complex structure (Figure 7b, black curves) in very good agreement with the characteristic WR R of a window discontinuity (Figure 7b, green curves): In particular, the compensated ridge function (Figure 7b, black line) now allows an accurate estimation of the window aperture from the dilation associated with the maximum amplitude (Figure 7d). The compensation of the attenuation-induced shape effects enables the recovery of the shape of the ridge functions but involves a GDF derivative order $\tilde{\alpha} = 4$, which differs from $\alpha = 5$ of the seismic source wavelet: This explains the slight difference between the thin-layer ridge functions plotted in Figures 7d (magenta and black lines) and 6d (blue line), respectively. It has to be noted that until the GDF wavelet has several vanishing moments (which is related to the derivative order)
greater than the regularity of the reflector discontinuity, this has no consequence on the quantitative descriptions based on the multi-scale seismic attributes.

## Extended application to seismic data reflected by complex geologic reflectors in an anelastic medium

To extend the analysis to a more realistic seismic data set, we now consider the complex reflectivity sequence defined from the in situ impedance log introduced in Ker et al. (2011), embedded in a K-F anelastic medium defined by a quality factor *Q*. We model the synthetic data set by using (1) a real seismic source signature related to the deep-towed seismic device SYstème SIsmique Fond de mer/ Deep-towed Seismic System (SYSIF) (Marssset et al., 2010) and (2) the nonstationary convolution of the Green’s function. The SYSIF source signature consists in a broadband chirp signal (220–1050 Hz). Applying the source-correction developed by Ker et al. (2013) to compute a WR from a broadband seismic source signature, we associated GDF wavelets of derivative order $\alpha = 5$ for

<figure><img src="image_9.png" /><figcaption>Figure 7. The WR of (a) a step-like and (b) a thin-bed reflectors embedded in an anelastic medium, processed by the time-varying shaping filter (black curves), in comparison with the WR performed for an elastic medium (green curves). The black and magenta lines stand for the ridge functions extracted from the corrected WR and the elastic WR respectively plotted in (c and d) for the step-like and thin-bed reflectors.</figcaption></figure>

the SYSIF seismic source in the dilation range (logarithm values) [-6.6; -5.7], equivalent to a peak-frequency range [315–830 Hz].

As a first step, we perform *R* and $\tilde{R}$, which stand for the WR of the complex acoustic impedance profile embedded in an elastic medium (Figure 8a) and in an anelastic medium (Figure 8b), respectively. The two particular reflectors extensively studied by Ker et al. (2013) are identified by a solid (A) and a dashed (B) line, respectively. The associated ridge functions are plotted in a log-log diagram (Figure 8c and 8d, respectively, where the blue and red curves correspond to the elastic and anelastic cases, respectively). As introduced in the section presenting the WR background, the aim is not to describe these multiscale seismic attributes, but to highlight the efficiency of the time-varying shape filtering compensation in recovering the WR properties required to assess the geometry of the realistic complex seismic reflectors.

In a second step, we apply the time-varying shape filtering compensation to the realistic seismic data set associated with $\tilde{R}$. We still consider the derivative order $\tilde{\alpha} = 4$ for the sake of simplicity in the interpretation of the ridge function and the time $\tilde{\tau} = 41.5$ ms, which is larger than the two-way traveltimes of the analyzed reflectors A and B. The less dilated SYSIF wavelet corresponds to a reference

<figure><img src="image_10.png" /><figcaption>Figure 8. The WR of a synthetic seismic data, involving a seismic source and an in situ reflectivity sequence embedded in (a) an elastic and (b) an anelastic medium, in which the blue and red lines stand for the ridge functions of reflectors A (solid line) and B (dashed line), plotted in (c and d), respectively.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9. (a) The WR of a synthetic seismic data processed by the time-varying shaping filter (black curves), in comparison with the WR performed for an elastic medium (green curves) for which ridge functions (black and magenta lines, respectively) of reflectors A and B are plotted in (b and c), solid and dashed lines, respectively.</figcaption></figure>

peak frequency $\tilde{f}_p = 468$ Hz. Based on these parameters, the result is the compensated WR $\tilde{R}$, which is in very good agreement with the WR $R$ related to the reference case with $\alpha = 4$ (Figure 9a, the black and green curves, respectively). The agreement is particularly highlighted when plotting the ridge functions in a log-log diagram (Figure 9a, black and magenta lines, respectively) with an error in the reconstruction lower than 5% for reflector A (solid lines) and 10% for reflector B (dashed lines): The time-varying shape filtering approach is thus very accurate to allow further quantitative analyses of broadband SYSIF data of complex reflectors in anelastic media.

## CONCLUSION

We have studied the effects of the intrinsic attenuation on a seismic source signal during its propagation into an anelastic medium. The analysis considers a seismic source signal defined as a GDF wavelet of derivative order $\alpha$ and peak frequency $f_p$, and a K-F anelastic medium, whose attenuation is $\beta \propto \omega^\delta$, with $\delta = 1$ and the proportionality factor is related to the reciprocal of the quality factor $Q$, considered constant with depth in the present study. We have shown that a GDF source wavelet that propagates into this anelastic medium is modified according to a fractional integration operation and a frequency shift; i.e., the attenuated signal is equivalent to another GDF. We have also demonstrated the dependency between the wavelet parameters and $Q$, by establishing an expression of their relationship.

As a consequence of these attenuation-related effects, the multiscale seismic attributes analysis based on the WR of a seismic reflector embedded in an anelastic medium is not straightforward. In particular, the ridge functions cannot be used to characterize the geometry of the reflector. To go further into this study, in which a complex acoustic reflectivity induces tunning effects, for instance, and attenuation also modifies the shape of the seismic source wave-

let, we took advantage of the quasi-conservation of the Gaussian character of a GDF to compensate for these effects and to recover properly the quantitative properties of the multiscale seismic attributes. We designed a time-varying shaping filter based on fractional integration that enables to substitute the time-varying wavelet by a constant-shape wavelet. We demonstrated this approach first on a step-like and a thin-bed related reflectors, widely analyzed in the framework of the WR in elastic media, and then on realistic synthetic data modeled by using a deep-towed SYSIF source signature and in situ acoustic-impedance log. The multiscale seismic attributes are accurately corrected from their attenuation-related distortions, while spanning in a smaller dilation range. The proposed approach allows extending the use of multiscale seismic attributes to anelastic media enabling us to properly analyze the geometry of the seismic reflector.

These results have been performed for an anelastic medium characterized by $\delta = 1$ and $Q$ homogeneous in space. It is interesting to note that similar developments should be extended to other models, such as power-law attenuation model with $1 < \delta < 2$. The potential of the approach in developing attenuation analysis methods such as $Q$ estimation has been introduced in the present work from the variations of the GDF parameters and deserves to be tested with in situ seismic transmission measurements, such as VSP data. The interest also deals with seismic source signature deconvolution, by providing a simple model for the seismic wavelet based on only three parameters that describe the time-varying shape of the seismic wavelet and can be used to extend the ongoing work dedicated to seismic wavelet modeling.

## ACKNOWLEDGMENTS

The authors would like to thank A. Gesret and two anonymous reviewers for their valuable suggestions that contributed to improving the manuscript.

## REFERENCES

* Abramowitz, M., and I. A. Stegun, 1972, Handbook of mathematical functions: U.S. Department of Commerce, National Bureau of Standards, US Government Printing.

* Caputo, M., 1967, Linear models of dissipation whose $Q$ is almost frequency independent. Part II: Geophysical Journal of the Royal Astronomical Society, **13**, 529–539, doi: **10.1111/j.1365-246X.1967.tb02303.x**.

* Castagna, J., S. Sun, and R. W. Siegfried, 2003, Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons: The Leading Edge, **22**, 120–127, doi: **10.1190/1.1559038**.

* Futterman, W. I., 1962, Dispersive body waves: Journal of Geophysical Research, **67**, 5279–5291, doi: **10.1029/JZ067i013p05279**.

* Gesret, A., M. Laigle, J. Diaz, M. Sachpazi, and A. Hirn, 2010, The oceanic nature of the African slab subducted under Peloponnesus: Thin-layer resolution from multiscale analysis of teleseismic P-to-S converted waves: Geophysical Journal International, **183**, 833–849, doi: **10.1111/j.1365-246X.2010.04738.x**.

* Heigl, W. M., 2007, Computing Gaussian derivative waveforms of any order: Geophysics, **72**, no. 4, H39–H42, doi: **10.1190/1.2716624**.

* Hosken, J. W. J., 1988, Ricker wavelets in their various guises: First Break, **6**, 24–33, doi: **10.3997/1365-2397.1988002**.

* Ker, S., Y. Le Gonidec, and D. Gibert, 2012, Multiscale seismic attributes: Source-corrected wavelet response and application to high-resolution seismic data: Geophysical Journal International, **190**, 1746–1760, doi: **10.1111/j.1365-246X.2012.05601.x**.

* Ker, S., Y. Le Gonidec, and D. Gibert, 2013, Multiresolution seismic data fusion with a generalized wavelet-based method to derive subseabed acoustic properties: Geophysical Journal International, **195**, 1370–1383, doi: **10.1093/gji/ggt317**.

* Ker, S., Y. Le Gonidec, D. Gibert, and B. Marsset, 2011, Multiscale seismic attributes: A wavelet-based method and its application to high-resolution seismic and ground truth data: Geophysical Journal International, **187**, 1038–1054, doi: **10.1111/j.1365-246X.2011.05207.x**.
Ker, S., Y. Le Gonidec, L. Marié, Y. Thomas, and D. Gibert, 2015, Multi-scale seismic reflectivity of shallow thermoclines: Journal of Geophysical Research, **120**, 1872–1886, doi: **10.1002/2014JC010478**.

Ker, S., Y. Le Gonidec, B. Marsset, G. Westbrook, D. Gibert, and T. A. Minshull, 2014, Fine-scale gas distribution in marine sediments assessed from deep-towed seismic data: Geophysical Journal International, **196**, 1466–1470, doi: **10.1093/gji/ggt497**.

Kolsky, H., 1956, The propagation of stress pulses in viscoelastic solids: Philosophical Magazine, **1**, 693–710, doi: **10.1080/14786435608238144**.

Le Gonidec, Y., and D. Gibert, 2006, The wavelet response as a multiscale characterization of scattering processes at granular interfaces: Ultrasonics, **44**, 381–390, doi: **10.1016/j.ultras.2006.05.212**.

Le Gonidec, Y., D. Gibert, and J.-N. Proust, 2002, Multiscale analysis of waves reflected by complex interfaces: Basic principles and experiments: Journal of Geophysical Research, **107**, 2184, doi: **10.1029/2001JB000558**.

Mallat, S., 1998, A wavelet tour of signal processing, 2nd ed.: Academic Press Inc.

Margrave, G. F., 1998, Theory of nonstationary linear filtering in the Fourier domain with application to time-varying filtering: Geophysics, **63**, 244–259, doi: **10.1190/1.1444318**.

Marsset, T., B. Marsset, S. Ker, Y. Thomas, and Y. Le Gall, 2010, High and very high resolution deep-towed seismic system: Performance and examples from deepwater geohazard studies: Deep-Sea Research I, **57**, 628–637, doi: **10.1016/j.dsr.2010.01.001**.

Morgan, E. C., M. Vanneste, I. Lecomte, L. G. Baise, O. Longva, and B. McAdoo, 2012, Estimation of free gas saturation from seismic reflection surveys by the genetic algorithm inversion of P-wave attenuation model: Geophysics, **77**, no. 4, R175–R187, doi: **10.1190/geo2011-0291.1**.

Müller, G., 1983, Rheological properties and velocity dispersion of a medium with power-law dependence of Q on frequency: Journal of Geophysics, **54**, 20–29.

Piété, H., L. Marié, B. Marsset, Y. Thomas, and M.-A. Gutscher, 2013, Seismic reflection imaging of shallow oceanographic structures: Journal of Geophysical Research, **118**, 2329–2344, doi: **10.1002/jgrc.20156**.

Ricker, N., 1953, The form and laws of propagation of seismic wavelets: Geophysics, **18**, 10–40, doi: **10.1190/1.1437843**.

Skauvold, J., J. Eidsvik, and U. Theune, 2016, A parametric model for seismic wavelets — With estimation and uncertainty quantification: Geophysical Journal International, **205**, 796–809, doi: **10.1093/gji/ggw029**.

Soubaras, R., and R. Dowle, 2010, Variable-depth streamer: A broadband marine solution: First Break, **28**, 89–96.

Tary, J. B., M. van der Baan, and R. H. Herrera, 2017, Applications of high-resolution time-frequency transforms to attenuation estimation: Geophysics, **82**, no. 1, V7–V20, doi: **10.1190/geo2016-0022.1**.

Toksöz, M. N., and D. H. Johnston, 1981, Seismic wave attenuation: SEG, 45.

Ursin, B., and T. Toverud, 2002, Comparison of seismic dispersion and attenuation models: Studia Geophysica et Geodaetica, **46**, 293–320, doi: **10.1023/A:1019810305074**.

van der Baan, M., 2012, Bandwidth enhancement: Inverse Q filtering or time-varying Wiener deconvolution?: Geophysics, **77**, no. 4, V133–V142, doi: **10.1190/geo2011-0500.1**.

Wang, Y., 2008, Seismic inverse Q filtering: Blackwell, Oxford.

Wang, Y., 2015a, The Ricker wavelet and the Lambert W function: Geophysical Journal International, **200**, 111–115, doi: **10.1093/gji/ggu384**.

Wang, Y., 2015b, Generalized seismic wavelets: Geophysical Journal International, **203**, 1172–1178, doi: **10.1093/gji/ggv346**.

Wang, Y., and J. Guo, 2004, Modified Kolsky model for seismic attenuation and dispersion: Journal of Geophysics and Engineering, **1**, 187–196.

Worthington, M. H., and J. A. Hudson, 2000, Fault properties from seismic Q: Geophysical Journal International, **143**, 937–944, doi: **10.1046/j.1365-246X.2000.00315.x**.

Xintao, C., S. Wang, S. Yuan, J. Zhao, L. Sun, and X. Wei, 2014, Sparse reflectivity inversion for nonstationary seismic data: Geophysics, **79**, no. 3, V93–V105, doi: **10.1190/geo2013-0313.1**.

Zhang, C., and T. J. Ulrych, 2002, Estimation of quality factors from CMP records: Geophysics, **67**, 1542–1547, doi: **10.1190/1.1512799**.