GETTING THE REFLECTION COEFFICIENTS FROM THE WAVES

The best starting point for inverse problems is the Kunetz equation (19).
 (50)

We need also the expression for the escaping wave (41)
 (51)

We also need to recall that . With this (50) becomes
 (52)

Multiplying through by A(Z) we get
 (53)

Since A(Z) is minimum-phase, A(Z) may be written as 1/B(Z) or A(1/Z) = 1/B(1/Z). Thus (53) becomes
 (54)

Identifying coefficients of zero and positive powers of Z as simultaneous equations, we get a set of equations which for a three-layer model looks like (r0 = 1).
 (55)

In (55) we see our old friend the Toeplitz matrix. It used to work for factoring spectra and predicting time series. Notice that -c3 has been inserted in (55) as the highest coefficient of A(Z). This is justified by reference back to the definition of A(Z) in terms of F(Z) and G(Z) which were in turn defined from the ck. It is by reexamining the Toeplitz simultaneous equations (55) and the Levinson method of solution () that we will learn how to compute the reflection coefficients from the waves.

The first four equations in (55) would normally be thought of as follows: Given the first three reflected pulses r1,r2, and r3 we may solve the equations for A, incidentally getting the reflection coefficient c3. Knowing A, the 5th equation in (55) may be used to compute r4. If the model were truly a three-layer model, it would come out right; if not, the discrepancy would be indicative of another reflector c4 which could be found by expanding equation (55) from 4th order to 5th order. In summary, given the reflected pulses rk, the Levinson recursion successively turns out the reflection coefficients ck.

Now suppose we begin by observation of the escaping wave E(Z). One way to determine the reflection coefficients would be to form 1 + R(Z) + R(1/Z) by the autocorrelation of E(Z); then, the Levinson recursion could be used to solve for the reflection coefficients. The only disadvantage of this method is that E(Z) contains an infinite number of coefficients so that in practice some truncation must be done. The truncation is avoided by an alternative method. Given E(Z) polynomial division will find A(Z). The heart of the Levinson recursion is the building up of A(Z) by Ak(Z) = Ak-1(Z) - ck ZkAk-1(1/Z). In particular, from () we have
 (56)

which shows how to get A3(Z) from A2(Z) and c3. To do it backwards, we see first that c3 is -a3. Then write (56) upside-down
 (57)

Next multiply (56) by 1/(1 - c32) and add the product to (57) multiplied by c3/(1 - c32). Notice that the upside-down vectors on the right-hand side cancel, leaving
 (58)

Equation (58) is the desired result which shows how to reduce Ak+1(Z) to Ak(Z) while learning ck+1. A program to continue this process is

       COMPLEX A,C,AL,BE,TOP,CONJG
C(1)=-1.; R(1)=1.; A(1)=1.; V(1)=1.
300   DO 310 I=1,N
310   C(I)=A(I)
DO 330 K=1,N
J=N-K+2
AL=1./(1.-C(J)*CONJG(C(J)))
BE=C(J)*AL
JH=(J+1)/2
DO 320 I=1,JH
TOP=AL*C(I)-BE*CONJG(C(J-I+1))
C(J-I+1)=AL*C(J-I+1)-BE*CONJG(C(I))
320   C(I)=TOP
330   C(J)=-BE/AL


A program to compute reflection coefficients ck from the prediction-error filter A(Z). The complex arithmetic is optional.

An inverse program to get R and A from c is

       COMPLEX C,R,A,BOT,CONJG
C(1)=-1.; R(1)=1.; A(1)=1.; V(1)=1.
100   DO 120 J=2,N
A(J)=0.
R(J)=C(J)*V(J-1)
DO 110 I=2,J
110   R(J)=R(J)-A(I)*R(J-I+1)
JH=(J+1)/2
DO 120 I=1,JH
BOT=A(J-I+1)-C(J)*CONJG(A(I))
A(I)=A(I)-C(J)*CONJG(A(J-I+1))
120   A(J-I+1)=BOT


This program computes both R and A from c.

 8-12 Figure 12 Earthquake seismogram geometry with white light incident from below. In the top layer, the sum of the waves vanishes representing zero pressure at the free surface. The difference of up- and down- is the observed vertical component of velocity.

Finally, let us see how to do a problem where there are random sources. Figure 12 shows the earthquake geometry." However, in order to introduce a statistical element, the pulse incident from below has been convolved with a white light series wt of random numbers. Consequently, all the waves internal to Figure 12 are given by the convolution of wt with the corresponding wave in the impulse incident model. Now suppose we are given the top-layer waves D = 1 -U = XW and wish to consider downward continuation. We have the layer matrix
 (59)

which can be re-written as
 (60)

The Burg prediction-error scheme can be written in the form
 (61)

which makes it equivalent within a scale factor to downward continuing surface waveforms. The remaining question is whether Burg's estimate of the reflection coefficient, namely,
 (62)

turns out to estimate the reflection coefficient ck in the physical model. To see how Burg's is related to the ck arising in the Levinson recursion, we define and for k=2 as
 (63)

Next from the dot product
 (64)
Now utilize the fact that (1,a1) satisfies the system. Following the Levinson recursion (64) can be written as
 (65)
Likewise we can deduce that .Thus, the Levinson calculation of the reflection coefficient can be written as
 (66)
The Burg treatment differs from the Levinson treatment in that Burg omits end-effect terms on (63). Instead of (63) he has
 (67)
For a sufficiently long data sequence the Burg method and the Levinson technique thus become indistinguishable. For a data sample of finite duration we must make a choice. The Levinson technique with (63) is equivalent to assuming the data sample vanishes off the ends of the interval in which it is observed. In most applications this is untrue, and so the Burg technique is usually preferable.

EXERCISES:

1. An impulse and the first part of a reflection seismogram, that is, 1 + 2R(Z) is .What are the first three reflection coefficients? Assuming there are no more reflectors what is the next point in the reflection seismogram?
2. A seismogram X(Z) = 1/(1 - .1Z + .9Z2) is observed at the surface of some layers over a halfspace. Sketch the time function and indicate its resonance frequency and decay time. Find the reflection coefficients if X(Z) is due to an impulsive source of unknown magnitude in the halfspace below the layers.
3. A source b0 + b1Z deep in the halfspace produces a seismogram .What are the layered structure and the source time function?