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/mx/inverseQfiltering.html

 

 

In this paper I discuss a stabilizarion 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 γ=(πQr)-1 , and cr and Qr are the phase velocity c(w) and the Q(w) value, respectively, at an arbitrary reference frequency. The tuning parameter wh 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/cr, 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 cr 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 w0=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(w0) for traveltime τ of the plane wave with dominant frequency w0, the kr term used in equation (4.8) can be expressed as

 

 

                                    (1.10)

 

so that it is independent of the velocity c(w0). 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:

 

                                               (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 wq 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 wq, with cosine-square tapering;

 

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

 

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 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. This aspect is discussed m the following stabilization scheme.

 

1.3 Stabilized inverse Q filter

 

To further improve the performance of inverse Q filtering, we now propose a stabilized approach to the wavefield downward continuation, where Q is a 1-D function, Q(r), varying with depth-time τ.

 

Considering downward continuation from the surface τ0 = 0 to the depth-time level τ using equation (1.7), we may express the wavefield U(τ,w) as

 

          (1.17)

Where

 

γ(τ) = 1/πQ(τ)                                                                                                (1.18)

 

To stabilize the implementation, we rewrite equation (1.17) as

 

                                         (1.19)

 

where

                                                          (1.20)

 

Solving equation (1.19) as an inverse problem with stabilization, we derive the following amplitude-stabilized formula:

 

                                            (1.21)

Where

 

                                                                   (1.22)

 

and σ2  is  a stabilization factor, a real positive constant used to stabilize the

solution. Equation (1.21) is the basis for an inverse Q filter in which

downward continuation is performed on all plane waves in the frequency

domain. Finally, we sum these plane waves to produce the time-domain

seismic signal,

 

                   (1.23)

 

We refer to this expression as stabilized inverse Q filtering.

 

Stabilized  inverse  Q  filtering,   equation  (1.23),   must be  performed successively to each time sample τ. Therefore, we may discretize it to

                                        (1.24)

 

 

 or present it in a vector-matrix form as

 

x = Az,                                    (1.25)

 

where x ={(u(τi)} is the time-domain output data vector, z = {U(wj )} is the frequency-domain input data vector and A ( M x N ) is the inverse Q filter with elements defined as

 

    (1.26)

 

in which A(τi,wj) is the stabilized amplitude-compensation coefficient, and the exponential term is the phase-correction term of the inverse Q operator.

 

References: Yanghua Wang:Seismic inverse Q Filtering, Blackwell Publishing 2008,

Stolt, R.H. (1978) Migration by Fourier Transform, Geophysics 43, 23-48

Hargreaves and Calvert (1991): Inverse Q-filtering by Fourier Transform, Geophysics 56,519-27

Bano, (1996): Q-phase compensation og seismic records in te frequency domain, Bulletin of the seismological  Society of America 86, 1179-86

Claerbout, (1976): Fundamentals of Geophysical Data Processing, McGrawHill Book Co, New York

 

 

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.