seismic convolution
The
convolution process uses the phase and amplitude obtained earlier,
so if you don't want to look at the sweep wavelet, the above steps
are not necessary.
Convolution of a Wavelet with the Reflection Coefficients
With the wavelet and reflection coefficient strings both in the
time domain, we are now able to create a synthetic seismogram.
The process is called convolution. In the time domain, the following
steps are needed:
1. reverse the order of the coefficients of the wavelet.
2. multiply each of the W values of the reversed wavelet with
the first W values of the reflection coefficient string and sum
the results; this is the first value of the output trace. (W is
the number of samples in the wavelet).
3. shift the wavelet down the coefficient string by one sample
and multiply/sum as in step 2; this is the second output sample.
4. repeat until the wavelet has been shifted to the end of the
reflection coefficient string.
5. to prevent unwanted end effects, the reflection coefficient
string must be padded with zeros at both ends before the convolution
is started; the length of the pad is equal to the wavelet length.
NUMERICAL
EXAMPLE:
Wavelet
String: |
0 |
5 |
10 |
0 |
-2 |
-1 |
0 |
|
|
|
|
|
|
|
|
|
|
|
Output
|
Reversal
String: |
0 |
-1 |
-2 |
0 |
10 |
5 |
0 |
|
|
|
|
|
|
|
|
|
|
|
Trace |
Reflection
String: |
0
0 0 0 0 0 0 2 0 0 0 1 0 0 0 0 0 0
|
V |
Shift
0: |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
|
|
|
|
|
|
|
|
0 |
1: |
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
|
|
|
|
|
|
|
0 |
2: |
|
|
0 |
0 |
0 |
0 |
0 |
10 |
0 |
|
|
|
|
|
|
|
|
|
10 |
3: |
|
|
|
0 |
0 |
0 |
0 |
20 |
0 |
0 |
|
|
|
|
|
|
|
|
20 |
4: |
|
|
|
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
|
|
|
|
0 |
5: |
|
|
|
|
|
0 |
0 |
-4 |
0 |
0 |
0 |
0 |
|
|
|
|
|
|
-4 |
6: |
|
|
|
|
|
|
0 |
-2 |
0 |
0 |
0 |
5 |
0 |
|
|
|
|
|
3 |
7: |
|
|
|
|
|
|
|
0 |
0 |
0 |
0 |
10 |
0 |
0 |
|
|
|
|
10 |
8: |
|
|
|
|
|
|
|
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
|
|
|
0 |
9: |
|
|
|
|
|
|
|
|
|
0 |
0 |
-2 |
0 |
0 |
0 |
0 |
|
|
-2 |
10: |
|
|
|
|
|
|
|
|
|
|
0 |
-1 |
0 |
0 |
0 |
0 |
0 |
|
-1 |
11: |
|
|
|
|
|
|
|
|
|
|
|
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
Notice
how the tail end of the first reflection interferes with the beginning
of the second reflection at data point 6.
This
technique, while simple to program, is computer intensive. The
number of multiplications equals N squared, where N is the length
of the padded reflection coefficient string.
A
more efficient technique uses the Fast Fourier Transform, which
converts the time series wavelet and reflection coefficient strings
to their equivalent frequency domain amplitude and phase spectra.
Convolution of the two involves the following:
1. calculate phase and amplitude of wavelet, using FFT algorithm.
2. calculate phase and amplitude of reflection coefficients, using
FFT algorithm.
3. multiply the two amplitude spectra and add the two phase spectra.
4. convert these spectra to the time domain using the inverse
FFT algorithm.
This
method uses N * log2(N) multiplications instead of N squared.
For example, a 5000 point transform takes 61,438 multiplications,
while the direct approach takes 25,000,000, a speed improvement
of over 600 times. Thus a computer which did the FFT in one minute
would take 6.7 hours to compute longhand. For a million data points,
the long method would take four weeks compared to a one minute FFT.
A
sample program used for performing the FFT is shown below. As
written, the program uses single precision numbers, and results
may be a bit noisy. Real programs use double precision for the
important array variables.
1
REM **********************************************************************
3 REM FAST FOURIER TRANSFORM PROGRAM - FORWARD
5 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
23 REM **** CAUTION: GWBASIC sets true = NEGATIVE!!!!!!!!!!!!!!!!
24 REM only on arithmetic operators - line 4070 needs fixes if
true = POSITIVE
25 REM **** CAUTION: GWBASIC expects angles in radians
40 OPTION BASE 1 'array indexes start at 1
41 PI = 3.141593 'single precision
42 WVLTLENGTH = 256 'data points
43 SAMPRATE = 1 'milliseconds
44 SWEEPLENGTH = 1 'seconds
45 INITPHASE = 0 'radians
80 FFTLENGTH = WVLTLENGTH/2
83 L = LOG(FFTLENGTH) / LOG(2) 'number of cycles in FFT
84 M = FFTLENGTH 'abbreviated name for use in FFT
85 FREQINT = 1000 / (2^(L+1)) / SAMPRATE 'output frequency sample
rate from FFT
86 DIM R [M], I[M], WVLT[WVLTLENGTH] 'Real and Imag arrays for
FFT
4430 P1 = 0: IF FFTFLAG$ = "IFT" THEN P1 = 1: GOTO 4800
'Inverse Transform
4440 PRINT "FFT Loop 1": K = 0: FOR J = 1 TO M-1: I
= 2 'FFT Starts Here
4455 IF K < M/I THEN GOTO 4470
4460 K = K - M / I: I = I + I: GOTO 4455
4470 K = K + M / I: IF K <= J THEN GOTO 4510
4480 A = R[J+1]: R[J+1] = R[K+1]: R[K+1] = A 'Binary Sort
4490 A = I[J+1]: I[J+1] = I[K+1]: I[K+1] = A
4510 NEXT J: G = .5: P = 1
4511 PRINT "FFT Loop 2": FOR I = 1 TO L: G = G + G:
C = 1: E = 0
4520 Q = ((1-P)/2)^.5 * (1-2*P1)
4522 IF I = 1 THEN P = -1 ELSE P = ((1+P)/2)^.5
4525 FOR R = 1 TO G
4530 FOR J = R TO M STEP G+G: K = J + G
4540 A = C*R[K] + E*I[K]
4550 B = E*R[K] - C*I[K]
4560 R[K] = R[J] - A
4570 I[K] = I[J] + B
4580 R[J] = R[J] + A
4590 I[J] = I[J] - B
4600 NEXT J: A = E*P + C*Q: C = C*P - E*Q: E = A: NEXT R: NEXT
I
4610 IF FFTFLAG$ = "IFT" THEN RETURN 'End of IFT
4801 REM **********************************************************************
4803 REM FAST FOURIER TRANSFORM PROGRAM - INVERSE
4805 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4810 PRINT "FFT Loop 3" 'FFT Continues
4820 A = 4 * ATN(1) / M: P = COS(A) 'Start IFT
4830 Q = (1 - 2*P1) * SIN(A): A = R[1]: R[1] = A+I[1]: I[1] =
A-I[1]
4835 IF FFTFLAG$ = "IFT" THEN GOTO 4850
4840 R[1] = R[1]/2: I[1] = I[1]/2
4850 C = 1 - 2*P1: E = 0: FOR J = 2 TO M/2: A = E*P + C*Q: C =
C*P - E*Q:
E = A: K = M - J + 2
4855 A = R[J] + R[K]
4860 B = (I[J] + I[K])*C - (R[J] - R[K])*E: U = I[J] - I[K]
4865 V = (I[J] + I[K])*E + (R[J] - R[K])*C: R[J] = (A+B)/2: I[J]
= (U-V)/2
4870 R[K] = (A-B)/2: I[K] = -(U+V)/2: NEXT J
4875 I[M/2+1] = -I[M/2+1]: IF FFTFLAG$ = "IFT" THEN
GOTO 4440 'IFT Continues^
4880 FOR I = 1 TO M: R[I] = R[I]/128: I[I] = I[I]/128: NEXT I
4885 RETURN 'End of FFT
4901 REM **********************************************************************
4903 REM RECREATE WAVELET FROM PHASE AND AMPLITUDE
4905 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4907 PRINT "Recreating Wavelet"
4910 FOR I = 1 TO FFTLENGTH: R[I] = A[I]*COS (P[I]*PI/180): I[I]
= P[I]*SIN(P[I]*PI/180): NEXT I
4915 FFTFLAG$ = "IFT": GOSUB 4400 'Do IFT on Wavelet
4920 PRINT "FFT Loop 4R": RECSCALE = 5000*SAMPRATE
4925 FOR I = 1 TO FFTLENGTH: WVLT[2*I-1] = -I[I]/RECSCALE: NEXT
I
4926 FOR I = 1 TO FFTLENGTH: WVLT[2*I] = -R[I]/RECSCALE: NEXT
I
4927 WVLT[1] = 0
4930 PRINT "Plotting Recreated Wavelet"
4935 VIEW (300,10)-(500,110)
4940 WINDOW (-20,-2)-(128,2)
4945 LINE (0,0)-(0,0),0
4950 FOR I = 1 TO LENGTH: LINE -((I-1) * SAMPRATE, WVLT[I]),7:NEXT
I
5000 REM **********************************************************************
|