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:

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

The basic of inverse Q-filtering is explained in my article on Wikipedia:

1.2   Numerical instability of inverse Q filtering

For the discussion of stabilization of inverse Q-filtering I will refer to Wang’s book that is available as a Google book:

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:

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.

I have studied some more in the field of inverse Q-filtering in this article:

Some aspects of seismic inverse Q-filtering theory by Knut Sørsdal

As seen above the problem is not the phase correction term that is always stable but the amplitude correction term in the theory of inverse Q-filtering that is unstable as shown on figure 2 above. Wang introduced a stabilization scheme for this.

Stabilized inverse Q filter

To further improve the performance of inverse Q filtering, Wang  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.8.a from my article on Wikipedia), he express the wavefield U(τ,w) as

(1.17)

Where

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

To stabilize the implementation, he 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,

Knut Sørsdal: http://en.wikipedia.org/wiki/Seismic_inverse_Q_filtering