Stabilized inverse Q filtering
algorithm
Cand.Real.
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. I have not put my
highly personal research in the subject out on internet, so this
papers here on internet are just a few discussing what is already
published by others. 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 in different ways. I use the Riccati
equation, but other and similar downward continuation algorithms can do the job
as well 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. Also Wangs implementation procedure is based on the theory of wavefield downward continuation. 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:
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,
(1.6)
This summation is referred to as the imaging
condition, as in seismic migration. Equations (1.5) and (1.6) 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
(1.7)
where S(w) is
the Fourier transform of a source waveform s(t), defined as the
real-valued Ricker wavelet,
(1.8)
with the dominant
radial frequency w0=100π (i.e. 50 Hz). The Ricker wavelet is a
symmetric,
Given the travel distance r = τc(w0)
for traveltime τ of the plane wave with dominant frequency w0,
the kr term used in equation (1.7)
can be expressed as
(1.9)
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
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
In equation (1.5), the amplitude compensation operator
is
(1.10)
If we set A(w) = 1,
equation (1.5) becomes an inverse Q filter for phase correction only:
(1.11)
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.12)
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.13)
That is, the (accumulated) exponent of the amplitude
factor is
(1.14)
where τ = ∑∆τ
and wq is a time-varying
frequency limit.
Wang
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.15)
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. Figure 4.2 shows that (b) is 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
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
I
take out from my thesis on the Riccati eqution (University of
In
my
(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. (
In my
theses, Sørsdal (2008), I have done some c
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 c
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
In my c
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
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.1.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 γ=(πQr)-1
According to Kolsky and Futterman
we are free to choose wr 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 wr 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 wr 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 τ3 is the strain relaxation time, τ4 is the
stress relaxation time and MR 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
τ3 > τ4. 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 τ3
> τ4 and the phase
velocity:
(3.4)
From (3.3) and (3.4) we can draw a
very important conclusion about this model. Since τ3 > τ4 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.
3.2.Q models, phase
velocity and minimum phase
Wang set up a list of Q-models in two categories:
A.Phase velocity goes to zero when
frequency approaches zero.
1.The modified Kolsky
model (linear attenuation)
2.The Strick-Azimi model
(power-law attenuation)
3.Kjartansson model (constant Q)
4.Azimi’s second and third models models
(non-linear attenuation)
5.Müller’s model (power-law Q)
B.Phase velocity goes to a finite, non-zero
phase velocity at zero frequency.
6.The Zener model (the
standard linear solid.)
7.The Cole-Cole model (a general linear-solid)
8.A new general linear model
An exception is the Cole-Cole model, which can give either zero or
non-zero phase velocity at zero frequency.
Actually these two categories divides the models
into a category that makes a causual solution and a
non-causual solution of the wave-equation. The
division can be seen already in the unmodified and modified Kolsky-model.
The original Kolsky-model can be represented with (3.5):
The Kolsky model (Kolsky,
1953) assumes the attenuation α(w) to be strictly
linear with frequency over the range of measurement:
(3.5)
And defines the phase velocity as (Kolsky
1956):
(3.6)
wr is a frequency below the lowest frequency in the
frequency band used. Then the second term in the parenthesis in (3.6) will
increase for increasing frequency, and give a nonzero solution when frequency
goes to zero. Then this model falls in the category giving a non-causal
solution.
On the other hand the modified Kolsky model
gives the equation:
(3.7)
In this model wh is the highest frequency in
the range. Then the second term in the equation (3.7) will behave as giving a
casual solution of the wave equation.
Wang has simplified (3.6) and (3.7), and we have the expression for
(3.7):
where
γ=1/(πQr)
We get a similar expression for (3.6). This equation
express the above discussion even simpler. With wh
as the highest frequency, c(w) will increase with
increasing w and go to zero when w approximate zero. With wr
as the lowest frequency c(w) will increase with
increasing frequency, and both wrand w
will go to zero at the same time, giving a nonzero value for c(w). (3.6) is a
good model for Q-models
(6-8) and for non-causual solutions. (3.7) is a model
for Q-models (1.5).
Above, in 3.1, in the discussion of the Zener
model we saw that this model goes in the model type of (3.6) and we got a non-causual solution.
4.1.Using a lowpass-filter
As mentioned in three points in discussion of (1.15), we can test a lowpassfilter in the inverse Q-filtering procedure. Phase
compensation is always stabile, therefore we do not
need to lowpassfilter in this case. However for
numerical convenience we will use same truncation for attenuation and phase
(point 1 from Wang)
We will apply a
Butterworth filter on our seismic model as presented on fig.4.1.
Fig.4.1. Impulse response of our seismic mode, green are
reflectors, red are damped with Q=50.
Inverse Q-filtering the trace will give strong artifacts due to
instability.
Fig.4.2. Model from fig.2.1. invers
Q-filter for Q=50 is applied.
Fig.4.2. shows that
inverse Q-filtering introduce strong unstability with artifacts. If we apply the Butterworth
filter on fig.4.3, fig.4.4 shows that the artifacts are very effectively
removed and we get a very good filtering.
Fig.4.3.
A lowpass filter with cutoff f=120 Hz is applied.
MATlabs sptool
are used. The filteret has both phase (green) and
amplitude (blue) response. It is applied on impulsresponse
fig.4.2 in Direct form 2, Second order section, Zero phase IIR.
The result of filtering
is remarkable. Fig.4.4 shows that all unstable energy are
removede. We have got extra energy before first
reflector, but we can remove non-causal energy after filtering by different
methods (ie. Tapering as mentioned of Wang)
Fig.4.4. filteredt impulsresponse wherer all unstabil energy is removed.
4.2.We remove the Nyquist frequency
The very
effective filtering of our seismic model can easily be explained from looking
at the frequency response. Fig.4.5. shows frequency response for the model in
forward Q-filtered version. Wang stated that a trace would be unstable for a
certain Q-value and time t as presented on fig.1 above. Our calculation shows
that frequencies close to the Nyquist-frequency have
a tendency to be unstable..
Fig.4.5.Frequency response for forward
Q-filtered trace. We will expect Inverse Q-filtering will blow up energy
close to the Nyquist-frquency that is 125 Hz.
Fig.4.6
shows frequency response for inverse Q-filtering of our seismic model. We have
strong unstable energy right before Nyquist-frequency
at 125 Hz. Figure gives us the specter as a functions
of digital frequency (and not seismic frequency as on fig. 4.5). In this plot
we get a good picture of what happens around the Nyquist
frequency. This is in digital representation 0.5 and samplingfrequency
equal to 1. We
clearly see how the Nyquist frequency is unstable. After filtering the artefacts is
completely removed.
Fig.4.6. Frequency response as a
function of digital frequency before filtering.
Fig.4.7. Frequency response as a
function of digital frequency after filtering.
I have done several calculations on inverse Q-filtering using the Riccati-equation and believe my calculations
is a development of the subject of inverse Q.filtering.
In the future I will present these calculations on internet.