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: **

**http://en.wikipedia.org/wiki/Seismic_inverse_Q_filtering**

**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:

http://books.google.no/books?id=IpwAjT-F_TgC

**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.*

**Further reading by Knut Sørsdal:**

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

http://bki.net/ricc/inverseQfilter

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(w_{j} )} 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},w_{j}) *is the stabilized
amplitude-compensation coefficient, and the exponential term is the
phase-correction term of the inverse

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

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