GEO 365N/384S Seismic Data Processing Computational Assignment 5

Linear Radon transform (slant stack)

In the first part of the assignment, we will revisit the groundroll attenuation problem in the Alaska data and will try to solve it using the linear Radon transform (slant stack).

1. Change directory to hw5/radon.
2. Run

   scons -c
   

   to remove (clean) previously generated files.

   

   signal,signal2
   Figure 1. Signal and noise separation using (a) Fourier method from Assignment 2, (b) Radon method from this assignment.

   

   
   

3. To compare the Radon method with the previously used Fourier method, we will extract one shot gather from the Alaska dataset. Figure 1a shows the extracted shot gather and the result of signal-and-noise separation using the Fourier method from Emc-Hammer. To reproduce it on your screen, run

   scons signal.view
   

   

   radon
   Figure 2. From left to right: input shot gather, its Radon transform, data reconstructed by the inverse transform.

   

   
   

4. Figure 2 shows the input shot gather, its Radon transform (slant stack), and the gather reconstructed by the inverse transform. To make sure that the transform is invertible, this implementation uses the Toeplitz structure of the matrix coming from the normal equations in least-squares inversion (Kostov, 1990).

   To reproduce the figure on your screen, run

   scons radon.view
   

5. The asymptotic theory of the generalized Radon transform predicts that a seismic event $T(x)$ will be transformed to an event $\tau(p)$ in the transform domain. The transformation follows from the tangent condition
   

   $$\left\{\begin{array}{rcl} T(x) & = & \theta(x;\tau,p)\;, \\ ... ...c{\partial}{\partial x} \theta(x;\tau,p)\;, \end{array}\right.$$ (1)

   
   
   where $\theta(x;\tau,p)$ describes a family of stacking curves. Prove that, when the stacking curve is a line,
   

   $$\theta(x;\tau,p) = \tau + p\,x\;,$$ (2)

   
   
   a hyperbolic seismic event
   

   $$T(x) = \displaystyle \sqrt{t_0^2 + \frac{x^2}{v_0^2}}$$ (3)

   
   
   will be transformed into an ellipse.

   

   Can you spot ellipses in the Radon transform domain in Figure 2?

   

   radon-cut
   Figure 3. Separating noise from signal in the Radon domain. From left to right: the Radon transform, separated noise components, noise reconstructed by the inverse transform.

   

   
   

6. To isolate dipping noise events associated with the groundroll, we can carve them out from the Radon domain according to their slopes. The selected cut and its inverse Radon transform are shown in Figure 3. To reproduce this figure on your screen, run

   scons radon-cut.view
   

7. The final result of separating the signal and noise using the Radon transform is shown in Figure 1b. You can reproduce it on your screen by running

   scons signal2.view
   

   or compare it with the Fourier-domain result by running

   sfpen Fig/signal.vpl Fig/signal2.vpl
   

8. Name some theoretical advantages and disadvantages of using the Radon transform in comparison with the Fourier transform.

   

9. The separation of signal and noise in Figure 1b is not optimal: there are remaining pieces of noise in the estimated signal and remaining signal in the estimated noise. Your creative task is to try improving this results by changing the parameters of the Radon transform or changing the way the noise is separated from the signal in the transformed domain.

from rsf.proj import *

# Download one shot gather
Fetch('shot.HH','alaska')

Flow('shot','shot.HH','dd form=native')
Plot('shot','grey title="Selected Shot" clip=2')

# Previous method: Fourier transform
####################################

Flow('fft','shot','fft1 | fft3')
Plot('fft',
       '''
       window max1=100 | math output="abs(input)" | real | 
       grey allpos=y title="Fourier Transform" 
       ''')
              
v1=9
v2=11

# Program from Emc-Hammer 

prog = Program('filter.c')
    
Flow('filter','fft %s' % prog[0],
     './${SOURCES[1]} logis=30 v1=%g v2=%g type=1' % (v1,v2))

Result('filter',
       '''
       window max1=100 | math output="abs(input)" | 
       real | grey allpos=y title="Filtered" 
       ''')

Flow('mute','fft %s' % prog[0],
     '''
     math output=1 | ./${SOURCES[1]} v1=%g v2=%g | real
     ''' % (v1,v2))

Result('mute','window max1=100 | grey allpos=y tile="Mute"')

Flow('signal','filter','fft3 inv=y | fft1 inv=y')
Plot('signal','grey title=Signal clip=2')

Flow('noise','shot signal','add scale=1,-1 ${SOURCES[1]}')
Plot('noise','grey title=Noise clip=2')

Result('signal','shot signal noise','SideBySideAniso')

# New method: Radon transform
#############################

# slope sampling
pmax=0.4
np=401
dp=2*pmax/(np-1)

Flow('radon','shot',
     'radon adj=y spk=n p0=%g dp=%g np=%d' % (-pmax,dp,np))
Plot('radon',
     '''
     grey title="Radon Transform" 
     label2=Slope unit2=s/kft clip=0.1
     ''')

Flow('inv','radon','radon adj=n nx=96 dx=0.11 ox=-5.225')
Plot('inv','grey title="Inverse Radon Transform" clip=2')

Result('radon','shot radon inv','SideBySideAniso')

# minimum and maximum slope for the noise
p1=0.15
p2=0.3

Flow('cut','radon',
     '''
     math output=1 |
     cut min2=%g max2=%g | 
     cut max2=%g | cut min2=%g 
     ''' % (-p1,p1,-p2,p2))
Flow('radon-cut','radon cut','mul ${SOURCES[1]}')
Plot('radon-cut',
     '''
     grey title="Noise Radon Transform" 
     label2=Slope unit2=s/kft clip=0.1
     ''')

Flow('noise2','radon-cut','radon adj=n nx=96 dx=0.11 ox=-5.225')
Plot('noise2','grey title=Noise clip=2')

Result('radon-cut','radon radon-cut noise2','SideBySideAniso')

Flow('signal2','shot noise2','add scale=1,-1 ${SOURCES[1]}')
Plot('signal2','grey title=Signal clip=2')

Result('signal2','shot signal2 noise2','SideBySideAniso')

End()

GEO 365N/384S Seismic Data Processing Computational Assignment 5