Stabilized inverse Q filtering algorithm

Cand.Real. Knut Sørsdal, University of Oslo

Seismic inverse Q-filtering is a data-processing technology for enhancing the resolution of seismic images. I have written some papers on the subject and on other subjects in seismic theory. URL’s to my different websites are:

Research website http://bki.net/mx/forskmidt.html

University of Oslo: http://folk.uio.no/knutsor

A html-version of this paper: http://bki.net/ricc/xtra/inverseQfiltering.html

In this paper I discuss a stabilization method in inverse Q-filtering. The algorithm is taken from the book “Seismic inverse Q-filtering” by Yanghua Wang (2008).

An inverse Q filter consists of two components, for amplitude compensation and phase correction. Whereas the phase component is unconditionally stable, the amplitude compensation operator is an exponential function of the frequency and traveltime, and including it in inverse Q filtering may cause instability and generate undesirable artefacts in seismic data. We will come back to that later.

First we will discuss inverse Q filtering from the point of Fourier transform. Inverse Q-filtering is a seismic data-processing technique for enhancing the image resolution. When a seismic wave propagates through the earth media, because of the anelastic property of the subsurface material, the wave energy is partially absorbed and the wavelet is distorted. As a consequence, it is generally observed in a seismic profile that the amplitude of the wavelet decreases and the width of the wavelet gradually lengthens along a path of increasing traveltime. An inverse Q filter attempts to compensate for the energy loss, correct the wavelet distortion in terms of the shape and timing, and produce a seismic image with high resolution.

An inverse Q filter includes two components, for amplitude compensation and phase correction. If an inverse Q filter considers phase correction only, it is unconditionally stable. When Q is constant in the medium, phase-only inverse Q filtering can be implemented efficiently as Stolt's (1978) wavenumber-frequency domain migration (Hargreaves and Calvert, 1991; Bano, 1996). This algorithm corrects the phase distortion from dispersion but neglects the amplitude effect. The amplitude-compensation operator is an exponential function of the frequency and traveltime, and including it in the inverse Q filter may cause instability and generate undesirable artefacts in the solution. Therefore, stability is the major concern in any scheme for inverse Q filtering.

It is desirable to have a stable algorithm for inverse Q filtering which can compensate for the amplitude effect and correct for the phase effect simultaneously, and does not boost the ambient noise. This chapter develops such a stable algorithm for inverse Q filtering. The implementation procedure is based on the theory of wavefield downward continuation (Claerbout, 1976). Within each step of downward continuation, the extrapolated wavefield, which is inverse Q filtered, is estimated by solving an inverse problem, and the solution is stabilized by incorporating a stabilization factor. In the implementation, the earth Q model is assumed to be a one-dimensional (1-D) function varying with depth or two-way traveltime.

1.1 Basics of inverse Q filtering

This section summarizes the basics of inverse Q filtering, to explain the new development in the following sections.

Inverse Q filtering can be introduced based on the 1-D one-way propagation wave equation

(1.1)

where U(r,w) is the plane wave of radial frequency w at travel distance r, k(w) is the wavenumber and i is the imaginary unit. Reflection seismograms record the reflection wave along the propagation path r from the source to reflector and back to the surface. Thus, in the following derivation we assume that the plane wave U(r,w) has already been attenuated by a Q filter through travel distance r.

Equation (1.1) has an analytical solution given by

U(r+∆r,w) = U(r,w) exp[ik(w)∆r] (1.2)

Considering only the positive frequencies, the complex wavenumber k defined by equation (3.27) in Wang(2008) becomes

(1.3)

where γ=(πQ_r)^-1 , and c_r and Q_r are the phase velocity c(w) and the Q(w) value, respectively, at an arbitrary reference frequency. The tuning parameter w_h is related to the highest possible frequency of the seismic band. Substituting this complex-valued wavenumber k(w) into solution (1.2) produces the following expression:

(1.4)

Replacing the distance increment ∆r by traveltime increment ∆τ = ∆r/c_r, equation (1.4) is expressed as

(1.5)

This is a basic inverse Q filter, in which the two exponential operators compensate and correct for, respectively, the amplitude effect (i.e. the energy absorption) and the phase effect (i.e. the velocity dispersion) of the earth Q filter.

Note that in derivation of equation (1.5), we assume that the reference velocity c_r approximately equals to the group velocity:

(1.6)

Equation (1.5) is the basis for inverse Q filtering in which downward continuation is performed in the frequency domain on all plane waves. The sum of these plane waves gives the time-domain seismic signal,

This summation is referred to as the imaging condition, as in seismic migration. Equations (1.5) and (1.7) must be applied successively to each time sample with sampling interval ∆r, producing u(τ) at each level.

1.2 Numerical instability of inverse Q filtering

Wang tests the basic downward continuation scheme described above on a group of noise-free synthetic seismic traces, to show the numerical instability of inverse Q filtering. He build a synthetic trace by

(1.8)

where S(w) is the Fourier transform of a source waveform s(t), defined as the real-valued Ricker wavelet,

(1.9)

with the dominant radial frequency w₀=100π (i.e. 50 Hz). The Ricker wavelet is a symmetric, noncausal wavelet. He use this symmetric wavelet to conveniently examine the phase change visually after the inverse Q filter is applied.

Given the travel distance r = τc(w₀) for traveltime τ of the plane wave with dominant frequency w₀, the kr term used in equation (4.8) can be expressed as

(1.10)

so that it is independent of the velocity c(w₀). In the example of Wang, he considered the signal to consist of a sequence of Ricker wavelets with τ = 100, 400, . . . , 1900 ms (increment 300 ms). (We refer to Figure 4.l.a in his book. The figure is also put in the notes of this article) It shows five synthetic traces with different Q values (Q = 400, 200, 100, 50 and 25) constant with depth in each case.

Figure 4.1.b (in the same book) displays the result of applying the inverse Q filter to the synthetic signals. For traces with Q = 400 and 200, the process restores the Ricker wavelet with correct phase and amplitude. However, there are strong artefacts as the Q value decreases and the imaging time increases, even though the input signal is noise-free. The appearance of noise in the output signal is a natural consequence of the downward continuation procedure: a plane wave is attenuated gradually, and beyond a certain distance the signal is below the ambient noise level, but the amplification required to recover the signal amplifies the ambient noise. In the data-noise-free case here, the background noise comes from the machine errors relative to working precision. This creation of strong artefacts is referred to as the numerical instability of the inverse Q filter.

In equation (1.5), the amplitude compensation operator is

(1.11)

If we set A(w) = 1, equation (1.5) becomes an inverse Q filter for phase correction only:

(1.12)

This phase-only inverse Q filter is unconditionally stable, as shown in Figure 4.1.c.in Wang’s book. Thus a restricted stability condition may be stated as

(1.13)

required for every sample interval ∆τ.

In practice, the artefacts caused by numerical instability can be suppressed by a low-pass filter. Based on experiments, Wang has found the following empirical stability condition:

˂=e (1.14)

That is, the (accumulated) exponent of the amplitude factor is not greater than 1. Equation (1.14) thus suggests an upper limit on the frequencies that can be included in the amplitude compensation,

(1.15)

where τ = ∑∆τ and w_q is a time-varying frequency limit.

Wang now test this low-pass filter with the following three methods:

1) both phase and amplitude compensation are truncated at frequency w_q, with cosine-square tapering;

2) phase compensation is performed on the full frequency band but with amplitude compensation only within the band limit (0, w_q );

3) both phase and amplitude are compensated on the full frequency band with the amplitude operator defined by

(1.16)

The results of these three schemes applied to the synthetic data set of Figure 4.l.a. in Wang’s book are shown in Figures 4.2.a-c in Wang’s book, respectively. We can see from Figure 4.2 that (b) shows an improvement over (a), but (c) is superior to both as the amplitudes for the traces Q < 200 are better compensated. In each case the numerical instability evident in Figure 4.1b (in Wang) is successfully suppressed. We may consider the third test scheme (Figure 4.2c) as a type of gain-limited inverse Q filter, but just for application to a seismic data set free of acquisition noise. For practical application, a gain limit should be set also according to the noise level in the seismic data set.

Figure. 1. The earth Q-filter and the inverse Q filter. a) Synthetic traces show the effect of the earth Q-filter with Q=400, 200, 100, 50 and 25. b) The inverse Q filtering (compensating for both phase and amplitude) result, which clearly indicates the numerical instability. C) The phase-only inverse Q filtered result. (Figure from Wang(2008))

2.1.Notes to Wang by Cand Real Knut Sørsdal, University of Oslo,Norway

I take out from my thesis on the Riccati eqution (University of Oslo 2008) a way to build synthetic trace very similar to those of Wang (Figure 1) above:

In my notation the 1-D one-way propagation wave equation has a solution:

(2.1)

Where u is the source waveform which in Wang’s paper is defined by the real valued Ricker-wavelet:

(2.2)

The amplitude compensation operator in (2.1) equivalent to the same in (1.11) is:

(2.3)

And the phase correction is given by the operator:

(2.4)

The condition for phase correction only is achieved very easily by setting Λ(w)=1 in (2.3) and this can be done by setting B(w) =0. (no absorption). Zero-phase is achieved when A(w)=1.

In my theses, Sørsdal (2008), I have done some calculations on delta-pulses on a trace 0,4 sec. These same calculations can be used on a traces up to 2000 ms and then be compared to the synthetic trace of Wang from fig.1 above.

Assuming that the term in (1.11) is equal to unity, I can connect Wang’s theory to my own very easily assuming zero-phase correction. This is the same I have done in my own theory by assuming A(w)=1. Since I have used the attenuation coefficient to express the factor B(w) as a constant B=0.023 in my calculations, we will look at the connection between Q and B by the relation B=π/Q. Thus we get Q=137. This value of Q will put our calculation very close to the instability limit of Wang’s calculations.

When we introduce phase correction by setting A(w)≠1 , the term deviates from unity, and the study of the different values of the term’s parameters compared to the viscoelastic models in my thesis (2008) in this case could be very interesting. In chapter 2 and 3 in his book Wang has discussed most of the models in the literature where he also includes the models I have used in my thesis.

In my calculations in chapter 8 in my thesis I have used B=0.023 (giving Q=137) and A=0,98. Since this values gives results that is close to the instability limit of Wang’s calculations concerning attenuation, and simply by looking at the pulses on fig.8.3 in my thesis I can see that pulses with arrival time from 0,4 sec and more will be causal, this could be a good choice of variables for a first comparison.

To achieve a minimum-phase solution of (2.1) the attenuation and dispersion term in the equation has to be related with Kramers-Krönig dispersion relation using Hilbert transform on the amplitude compensation operator (2.2) and the phase correction operator (2.3). A study of this is done in Wang (2008). I will not go further into this in my notes here, but simply use the values mentioned above.

a) B=0,001 A=0.98 earth Q-filter (blue) inverse Q-filter (red)

b) B=0.003 A=0.98 earth Q-filter (blue) inverse Q-filter (red)

c) B=0.005 A=0.98 undamped earth filter (red) inverse Q-filter (blue)

d) B=0.006 A=0.98 undamped earth-filter(red) inverse Q-filter (blue)

e) B=0.007 A=0.98 undamped earth-filter(red) inverse Q-filter (blue)

Fig.2.The earth Q-filter and the inverse Q-filter for different values of B. A=0.98 introducing phase-shift in the solution of (2.1). Instability occurs for B=0.005
Fig.1 shows graphs of calculations on traces up to 2000 ms. With B=0.005 we see the first signs of instability with energy accumulating at the end of the trace.(c). On (d) and (e) instability is even more dominant and we are not able to recover the original amplitude of the pulses. Fig.3 use an even higher B=0.008 for the earth Q-filter, that should give even more instability, but since we have set B=0 in the inverse Q-filter, we achieve a phase only inverse Q-filter and this introduce no instability. Fig.3 shows that the filter corrects completely for phase-shift, but do nothing to recover amplitude. Wang states that an inverse phase only Q-filter will always be stable. Fig.4-6 shows the calculations for some pulse-forms.

B=0.008 A=0.98 Earth Q-filter (blue) Undamped pulses (red)

B=0.008 A=0.98 Undamped pulses (red) inverse phase only Q-filter (blue)

Fig.3.upper graph: undamped pulses (red) are damped by attenuation and phase delayed by dispersion. lower graph: dispersion is corrected by inverse phase only Q-filter but pulse amplitude is not restored.

Pulse shaping and inverse Q-filtering – delta, Ricker and minimum-phase pulse

Fig.4.upper graph: delta puls initial pulse for the earth Q-filter(original pulse (blue)

Fig.5. Upper graph: Rickerpulse initial pulse for the earth Q-filter. Original pulse (red).

Lower graph: the phase shift is completely corrected by an inverse phase-only Q-filter, but the amplitude is not restored.

Fig.6. Upper graph: Minimum-phase initial pulse for the earth Q-filter. Original pulse (red).

Lower graph: the phase shift is completely corrected by an inverse phase-only Q-filter, but the amplitude is not restored.

3.Standard linear model (Zener model) and minimum phase

Wang has discussed the Kolsky (1956) model in a way Futterman (1962) described it. Wang introduced the equation representing both attenuation and dispersion:

) (3.1)

Where γ=(πQ_r)^-1

According to Kolsky and Futterman we are free to choose w_r following the phenomenological criterion that it be small compared with the lowest measured frequency w. If we use a Rickerpuls with minimum frequency 10 Hz we must choose w_r less than 10π. This gave a solution of the wave equation satisfying laboratory measurement done by them. However, the basic Kolsky model does not satisfy the minimum delay condition (minimum phase) in dispersive earth media. This is easy to see from the discussion above. If w_r is less than all other frequencies in the solution of (3.1), the term A(w) in our model will always be greater than 1 and the solution of the wave equation will give pulses arriving before actual arrival time.

We can show that this will be the case for all damping models that are developed from the standard linear model. The Zener (1948) or standard linear solid model is the most general linear equation in stress a, strain e and their first time derivatives a and s . In the Fourier transform domain, the stress and the strain are linked by the complex modulus defined as M(w) (Ben-Menahem and Singh, 1981):

 (3.2)

Where τ₃ is the strain relaxation time, τ₄ is the stress relaxation time and M_R is a deformation modulus with sub-index _R denoting the relaxed modulus.

The mechanical 'relaxation' means that the strain produced by the sudden application of a fixed stress to a material increases asymptotically with time. Similarly, the stress produced when the material is suddenly strained relaxes asymptotically (Kolsky, 1953). It is found that stress waves whose periods are close to the 'relaxation times' of such a medium are severely attenuated when passing through it.

To get damping we must assume that τ₃ > τ₄. This we can see by substituting them into the equations for attenuation and dispersion given earlier, and we obtain the attenuation (from Wang)

                    (3.3)

If we shall have positive attenuation τ₃ > τ₄ and the phase velocity:

(3.4)

From (3.3) and (3.4) we can draw a very important conclusion about this model. Since τ₃ > τ₄ the paranthesis on right side of expression (3.4) will always be less than one. That means that the wave number k  always will be less than one and A(w)>1. Because of that we will never get a causal solution of the wave equation using this model.