Seismic inverse Q filtering

Seismic inverse Q filtering is a data processing technology for enhancing the resolution of reflection seismology images. Q is the anelastic attenuation factor or the seismic quality factor, a measure of the energy loss as the seismic wave moves.

Basics
Seismic inverse Q-filtering employs a wave propagation reversal procedure that compensates for energy absorption and corrects wavelet distortion due to velocity dispersion. By compensating for amplitude attenuation with a model of the visco-elastic attenuation model type, seismic data can provide true relative-amplitude information for amplitude inversion and subsequent reservoir characterization. By correcting the phase distortion due to velocity dispersion, seismic data with enhanced vertical resolution can yield correct timings for lithological identification.

However, Wang's outline of the subject is excellent and to follow his path, inverse Q filtering can be introduced based on the 1-D one-way propagation wave equation. He introduce this equation:.


 * $$\frac { dU(r,w)}{ dr} - ik(w)U(r,w)=0 \quad (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. With this approach Wang assumes that the plane wave U(r,w) has already been attenuated by a Q filter through travel distance r. We must have this in mind when we go to the step of finding a solution of (1.1). It is necessary that the initial U(r,w) either is already created by a forward synthetic Q-filtering process or taken directly from seismic surface data. Wang has introduced this concept in chapter 5 in his book. I think it is necessary to have this in mind also when the inverse theory is developed. Equation (1.1) has an analytical solution given by


 * $$U(r+\bigtriangleup r,w) =U(r,w)\exp (ik(w)\bigtriangleup r) \quad (1.2)$$

Kolsky's attenuation-dispersion model
The wavenumber k(w) is an important variable in the solution (1.2). To obtain a solution that can be applied to seismic k(w) must be connected to a function that represent the way U(r,w) propagates in the seismic media. This functions can be regarded as a Q-model.

The Kolsky Model is used extensively in seismic inverse Q-filtering. The model assumes the attenuation α(w) to be strictly linear with frequency over the range of measurement:


 * $$\alpha=\frac {|w|}{(2 c_r Q_r)} \quad (1.3)$$

And defines the phase velocity as:


 * $$\frac {1}{c(w)} =\frac {1}{c_r} (1-\frac {1}{\pi Q_r} ln |\frac{w}{w_r}|) \quad (1.4)$$

Where cr and Qr  are the phase velocity and the Q value at a reference frequency wr.

For a large value of Qr >>1 the solution (1.4) can be approximated to


 * $$\frac {1}{c(w)} =\frac {1}{c_r} |\frac{w}{w_r}|^{-\gamma} \quad (1.5)$$

where


 * $$ \gamma =(\pi Q_r)^{-1} $$

Kolsky's model was derived from and fitted well with experimental observations. A requirement in the theory for materials satisfying the linear attenuation assumption is that the reference frequency wr is a finite (arbitrarily small but nonzero) cut-off on the absorption. According to Kolsky, we are free to choose wr following the phenomenological criterion that it be small compared with the lowest measured frequency w in the frequency band. Those who want a deeper insight into this concept can go to Futterman (1962)

Now we consider (1.3) as the real part of the wavenumber k(w) and (1.5) as the imaginary part. Considering only the positive frequencies, the complex wavenumber k(w) then becomes


 * $$k(w)=(1-\frac {i}{2Q(w)}) \frac {|w|}{c_r} | \frac{w}{w_r}|^{-\gamma}) \quad (1.6)$$

where cr and Qr are the phase velocity c(w) and the Q(w) value, respectively, at an arbitrary reference frequency. Substituting this complex-valued wavenumber k(w) into solution (1.2) produces the following expression:


 * $$U(r+\bigtriangleup r,w)=U(r,w)\exp[\frac {\bigtriangleup r}{2Q(w)} + i\frac{|w|\bigtriangleup r}{c_r}|\frac{w}{w_r}|^{-\gamma}]    \quad (1.7)$$

Replacing the distance increment ∆r by traveltime increment ∆t = ∆r/cr, equation (1.7) is expressed as


 * $$U(t+\bigtriangleup t,w)=U(t,w)\exp\bigg(|\frac{w}{w_h}|^{-\gamma} \frac{|w|\bigtriangleup t}{2Q(w)}\bigg) \exp \bigg( i|\frac{w}{w_h}|^{-\gamma} w\bigtriangleup t\bigg)    \quad (1.8.a)$$

This is a basic inverse Q filter with the Kolsky model. 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. In order for the inverse filter to compensate for a causal solution Wang modified the original Kolsky model by using wr as the highest frequency in the frequency band and call it wh.

As mentioned above, at the same time we start discussing the design and application of an inverse Q-filter, we need to specify a mathematical Q-model that can compute benchmark data similar to Wang's benchmark data that can be used for inversion later. Then we must regard the forward Q-filtering process as our solution of (1.2) which means we must compute the inversion of (1.2) as a Q-forward filtering process. Wang showed how easily this could be done by simply changing the sign before γ and Q.

Wang writes: We call the removal of the phase correction effect given by a previous phase-only inverse Q-filtering as phase-only forward Q-filtering. It would be rather straightforward if we consider it as an inverse problem of the inverse equation and take an inverse solution. But solving an inverse problem is often time-consuming because it involves the computation of an inverse matrix. Then Wang develops the theory both of amplitude forward Q-filtering and phase Q-filtering simply by changing sign before γ and Q, and introduce the equation:


 * $$U(t+\bigtriangleup t,w)=U(t,w)\exp \bigg(-|\frac{w}{w_r}|^{\gamma} \frac{|w|\bigtriangleup t}{2Q(w)}\bigg) \exp \bigg( i|\frac{w}{w_r}|^{\gamma} w\bigtriangleup t\bigg)    \quad (1.8.b)$$

In (1.8.b) I use the reference frequency in the original form given by Kolsky and use the notation wr. Equation (1.8.a and. b) is the basis for 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,


 * $$U(t+\bigtriangleup t) =\int_0^\infty U(t+\bigtriangleup t,w)dw.\quad (1.9)$$

This summation is referred to as the imaging condition, as in seismic migration. Equations (1.8.a) and (1.8.b) must be applied successively to each time sample with sampling interval ∆t, producing u(t) at each level. This last approach completes the downward continuation theory of q-filtering.

Computations
Fig.1 a and b shows a forward Q-filtering solution of (1.8.a-b integrated with (1.9)) for a single pulse presented for this article. This simple graph illustrates very well some aspects of the theory of Q-filtering. Black dotted graph is zero-phase solution of the equation (1.8.b) where we have amplitude only Q-filtering and no effect from the phase (Q=10). On fig.1.a. blue dotted line is (1.8.b) with fr=0.001 Hz. (We compute from the relation w=2:$$ \pi$$ f.) Red line is the solution of (1.8.a) with fh=10 Hz. We have no amplitude compensation for (1.8.a), looking only at the effect from the phase. Comparing with fig.1.b. we can see that the smaller the reference frequency fr in the forward solution, the more phase effect we get. In the inverse solution (1.8.a) we have the opposite effect, the higher value of fh the more the forward phase effect is compensated, toward a causal solution. The graphs give a good illustration of Q-filtering according to the theory of Wang. However, more calculations should be done and more seismic models should be introduced. The effect of amplitude inverse Q-filtering should also be studied.