Variable elimination

A seismic inversion is the process of transforming seismic reflection data into a quantitative rock-property description of a reservoir.When we do a seismic inversion with the Riccati equation we perform nonlinear computations to evaluate the density and seismic velocity distribution from the observed reflection seismogram.

1. Introduction

The first attempt with this inversion method was done by Gjevik et al. (1975). In another article Nielsen and Gjevik (1978) made more calculations with the method. The new thing with this second article was that absorption was included in the inversion.

2. Basics

It is an ultimate aim for explorational geophysicists to be able to evaluate the density and seismic velocity distribution from the observed reflection seismograms. This process is generally called the inversion of seismic reflection data. Because of the complexity of the geological structure, we are forced to be rather reserved about a general success of the inversion attempts. But in special favorable cases the geological structures are rather uncomplicated, and the propagation of seismic waves through these structures can be described fairly accurate by simple physical laws. In such cases an attempt at the inversion of reflection data might become successful. Gjevik et al. formulated an iterative inversion method, which was found to be suitable for numerical computations. They used the Riccati equation. They applied the method on synthetic reflection data, and found that it provides a very fast and accurate inversion.

In order to describe the wave field we introduce z as the depth of the seismic medium, and we denote time by t. We will study the reflection of a plane, longitudinal wave propagating through a medium where the density ρ and the velocity are functions of depth z. The displacement and the stress are represented respectively by W and p, and we have the equation of motion:

$\frac{d^{2}W}{d^{2}t}=\frac{1}{\rho }\frac{dp}{dz}(2.1)$

In the elastic case we assume a stress-strain relationship of the form:

$p=({\kappa }_0+\frac{4}{3}{\nu }_0)\frac{dW}{dz}(2.2)$

This is the stress-strain relationship in the one-dimensional case. We have depth dependency both in the compression modulus k₀ and in the shear modulus u₀. Bland (1960) has given a good discussion of this expression, and we have added his theory to that of Nilsen and Gjevik. We will introduce the correspondence principle. This principle says that we can introduce viscoelasticity into the stress-strain relationship equation (2.1) and (2.2). We simply replace the elastic modules k₀ and u₀ with complex visco-elastic modules and introduce the time Fourier transform of a function:

If we assume that the time Fourier transform of W and p exist, we obtain from (2.1) and (2.2):

$\frac{dp}{dz}=-\rho w^{2}W(2.3)$

$p=(\kappa (iw)+\frac{4}{3}\nu (iw))\frac{dW}{dz}(2.4.a)$

The stress-strain relationship in the equations (2.4.a) can be given many interpretations. Nilsen and Gjevik use a stress-strain relation of a 'firmo-viscous' substance known as the Kelvin-Voigt model, and introduced a stress-strain relation of the form:

$p=(\rho c^2+q\frac{d}{dt})\frac{dW}{dz}(2.4.b)$

q represent the rate of energy absorption. The relation can be written on a Fourier-transformed form where elastic modules have been replaced with viscoelastic of the Kelvin-Voigt type:

$p=\rho c_0^2{Q}^2\frac{dW}{dz}(2.5)$

The function Q represents depth and frequency-dependent absorption and has been rather complicated in Nilsen and Gjevik's paper. In this outline we can simply regard it as a single variable.

Now we split the wave field into two components, denoting the stress-amplitude for these components by D and R, seeking a solution for (2.4) and (2.5) of the form:

$P=D+R(2.6)$

$W=\frac{1}{iw{c}_{0}Q}(D-R)(2.6)$

D and R represents waves propagating respectively in down and upward direction in a viscoelastic medium. We assume this interpretation also holds for a weakly inhomogeneous medium. Invoking (2.6) we obtain from (2.3) and (2.5) after some manipulations:

$\frac{dD}{dz}=\frac{iw}{c_{0}Q}D+T(D-R)+(2.7)$

$\frac{dR}{dz}=\frac{iw}{c_{0}Q}RD+T(D-R)(2.8.a)$

Where the function T is defined by:

$T=\frac{1}{2}(\frac{1}{\rho c_0}\frac{d\rho c_0}{dz}+\frac{1}{\rho c_0}\frac{dQ}{dz})(2.8.b)$

In cases where absorption is small, the first term on the right hand side of this expression is always much larger than the second term, and we have:

$T=\frac{1}{2}(\frac{1}{\rho c_0}\frac{d\rho c_0}{dz})=\gamma (2.8.c)$

y is the well known relative impedance variation in a seismic medium. For waves used in reflection studies of rocks and sediments we can almost always use this expression. Multiplying the first of the equations (2.7) by the factor –R/D², the second factor with 1/D and adding the resulting equations, we find that the ratio K=R/D satisfies the Riccati equation:

$\frac{dR}{dz}=\frac{2iw}{c_{0}Q}K-\gamma (i-K^2)(2.8.d)$

K is a complex function of z and w. By the relation:

$t=2\underset{0}{\overset{\mathrm{\infty }}{\int }}dz,$

we introduce the two-way travel-time and from (2.8) we get the equation:

$\frac{dR}{dz}=\frac{2iw}{Q}K-\gamma (i-K^2)(2.9)$

where γ is a function of t.

If we assume that γ and K vanishes for t > T, equation (2.9) is equivalent to the integral:

$K(t,iw)=exp(\varphi (t,iw))\underset{t}{\overset{T}{\int }}\gamma (s)exp(-\varphi (s,iw))(1-K^2(s,iw))ds(2.10)$

where the phase function is defined by.

$\varphi (t,iw)=iw\underset{0}{\overset{t}{\int }}\frac{1}{Q}ds(2.11.a)$

Nilsen and Gjevik replaced the variable Q with a function that included a depth-dependent function h, and got the equation:

$\varphi (t,iw)=iw(\tau -i\underset{0}{\overset{\tau }{\int }}h(s)ds)(2.11.b)$

A more deep insight into the development of the Kelvin-Voigt model used in Nilsen and Gjevik's article can be found in Jaeger, J. (1962)

The reflection response of the reflecting layer is:

$K(0,iw)=\underset{0}{\overset{T}{\int }}\gamma (s)exp(-\varphi (s,iw))(1-K^2(s,iw))ds(2.12)$

Equation (2.10-11) are on the same form as equation (12-13) in Nilsen and Gjevik’s paper. And our phase function represents a viscoelastic absorption model of the Kelvin Voigt type. We now have a straight forward method to find the reflection response including absorption for a great variety of viscoelastic models when γ and the viscoelastic model is given.

The inversion procedure closely follows Nilsen and Gjevik’s paper. If the reflection response K and the viscoelastic properties of the medium are known γ can be computed by an inversion iteration procedure. We simply multiply (2.12) by 1/(2π) exp(iwt) and integrate with respect to w. By interchanging the sequence of integration, we obtain:

$\underset{t}{\overset{T}{\int }}\gamma (s)H(t-s,s)ds=\frac{1}{2\pi }\underset{0}{\overset{\mathrm{\infty }}{\int }}K(0,iw)exp(iwt)dw+\frac{1}{2\pi }\underset{t}{\overset{T}{\int }}\gamma (s)\underset{0}{\overset{\mathrm{\infty }}{\int }}{K}^{2}exp(iwt)exp(-\varphi (s,iw))ds(2.13)$

The function H in equation (2.13) represents attenuation and is equivalent with function H in Nilsen and Gjevik’s paper. In the case of no energy absorption, the function H approaches the delta function and the left hand side of equation (2.13) is simply γ .

Computations

The iterative inversion procedure is started by neglecting the second term on the right hand side of (2.13). Hence the integral on the right hand side can be computed, and if h(t) can be estimated, we may obtain γ(t) by deconvolution or seismic inverse Q filtering. The second order approximation to γ(t) is obtained by substituting the first order approximation to γ(t) and K in the second term on the right hand side of (2.13). The process may be continued as long as there are significant changes in the computed values of γ(t) . In cases where h=0 (no energy absorption) it has been shown (Gjevik et al. 1976) that this iterative procedure often provides a convergent sequence whereby the product of density and velocity (acoustic impedance) can be computed. Fig.1. shows the impedance that can be computed from γ (2.13) after 3 iterations with no absorption. Gjevik et al. (1976) explains very well how the impedance is computed from γ.

• 
  Fig.1.Impedance model computed With Riccati equation

Fig.2. shows a solution of inversion with the Riccati-equation and no absorption where we have included more data. Blue graph is the impedance computed from γ. Blue and red bars represents seismic layers corresponding to the solution of (2.12)-(2.13). Red represents a positive seismic event - blue is negative. There are four times more samplepoints in fig.2 than fig.1 so samplepoint 100 is the same time as samplepoint 400 in fig.2. Fig.2.a (left) presents 1 iteration for (2.13). This solution gives us the reflectors and can also be called the reflectivity. Fig.2.b. gives us the solution of (2.13) after 6 iterations and represents the 'synthetics' or the synthetic seismogram. Fig.2.c. gives us the solution of (2.13) after 6 iterations and is the actual seismic inversion.

A simple illustration of the procedure can be given by considering left graph as the starting position. Data from left graph is iterated so we got the solution in the middle. In the right graph we have operated with inverse computations on the data in the middle, trying to get a solution equal to the one to the left. We can see that we got a sequence of bars more like the graph to the left than the sequence of bars in the middle. So we have got a useful solution for our inversion procedure. The reason why our solution (right) is not identical to the solution (left) is partly due to machine calculation errors in the inversion procedure, and Gibbs phenomenon. A better solution could be found with correct filtering to correct this. The stability of the inversion will also be of importance for a correct solution. The last problem must be dealt with by correct scaling of the seismogram.

• 
  Fig.2.a) Reflectors/Reflectivity b) Synthetics c) Inversion.

Fig.3. shows a more complicated impedance that can be computed from γ (2.13) after 2 iterations with no absorption. Of special interest is the event immediately after the second barrier at 300 samplepoints on fig.4. In the synthetics the polarity of the bar representing the seismic layer has changed so we have two blue bars. It is restored back to one blue and one red in the inversion so we can compare with the reflectivity. Much of the energy caused by multiples in the synthetics for arrivals after 400 sample points are removed, and we have to some degree restored events from 300 to 400 sample points. So we have got a good inversion.

• 
  Fig.3.Impedance model computed with Riccati equation

• 
  Fig.4) Reflectors/Reflectivity b) Synthetics c) Inversion.

References

43 year old Petroleum Engineer Harry from Deep River, usually spends time with hobbies and interests like renting movies, property developers in singapore new condominium and vehicle racing. Constantly enjoys going to destinations like Camino Real de Tierra Adentro.

External links

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