Homework 2

Computational part 2

In the computational part, we begin working with field data. The left panel in Figure 3 shows a CMP (common midpoint) gather from the Gulf of Mexico (Claerbout, 2006).

cmp
Figure 3. From left to right: (a) CMP (common midpoint) gather with overlaid traveltime curves. (b) Interval velocity. (c) RMS (root-mean-square) velocity overlaid on the semblance scan. (d) CMP gather after normal moveout.

We will assume a $V(z)$ medium and will use a very simple model of the interval velocity to explain the geometry of the observed data. The model involves two parameters: the initial gradient of velocity and the maximum velocity. The velocity function starts at the water velocity of 1.5 km/s and grows linearly with vertical time until it reaches the maximum velocity, after which point it remains flat. The panels in Figure 3 show the interval velocity, the corresponding RMS velocity (overlaid on the semblance scan), and the CMP gather after NMO (normal moveout).

Your task is to find the best values of the two model parameters for optimal prediction of the traveltime geometry and for flattening the CMP gather after NMO.

1. Change directory

   cd hw2/cmp
   

2. Run

   scons cmps.vpl
   

   to generate and display a movie looping through different values of the maximum velocity. If you are on a computer with multiple CPUs, you can also try

   pscons cmps.vpl
   

   to generate different movie frames faster by running computations in parallel.
3. Edit the SConstruct file to modify the velocity gradient. Check your result by running

   pscons cmps.vpl
   

   again.
4. Edit the SConstruct file to select the best frame of the movie (corresponding to the best maximum velocity). Display it by running

   scons view
   

from rsf.proj import *
from rsf.prog import RSFROOT

# Program compilation
#####################

proj = Project()

# COMMENT THE NEXT LINE FOR FORTRAN OR PYTHON
prog = proj.Program('traveltime.c')

# UNCOMMENT BELOW IF YOU WANT TO USE FORTRAN
#prog = proj.Program('traveltime.f90',
#                   F90PATH=os.path.join(RSFROOT,'include'),
#                   LIBS=['rsff90']+proj.get('LIBS'))

# UNCOMMENT BELOW IF YOU WANT TO USE PYTHON
#prog = proj.Command('traveltime.exe','traveltime.py',
#                    'cp $SOURCE $TARGET')
#AddPostAction(prog,Chmod(prog,0o755))

exe = str(prog[0])

# Donwload data
Fetch('midpts.hh','midpts')

# Select a CMP gather, mute
Flow('cmp','midpts.hh',
     '''
     window n3=1 | dd form=native | 
     mutter half=n v0=1.5 |
     put label1=Time unit1=s label2=Offset unit2=km
     ''')
Plot('cmp','grey title="Common Midpoint Gather" ')

# Velocity scan
Flow('vscan','cmp',
     'vscan half=n v0=1.4 nv=111 dv=0.01 semblance=y')
Plot('vscan','grey color=j allpos=y title="Semblance Scan" ')

##############################
grad = 0.25 # Velocity gradient
##############################

cmps = []
for iv in range(21):
    vmax = 1.5+0.2*grad*iv

    # Interval velocity
    vint = 'vint%d' % iv

    Flow(vint,None,
         '''
         math n1=1000 d1=0.004 
         label1=Time unit1=s 
         output="1.5+%g*x1" | clip clip=%g 
         ''' % (grad,vmax))
    Plot(vint,
         '''
         graph yreverse=y transp=y pad=n plotfat=15 
         title="Interval Velocity" min2=1.4 max2=%g 
         wheretitle=b wherexlabel=t
         label2=Velocity unit2=km/s 
         ''' % (1.6+4*grad))

    # Traveltimes
    time = 'time%d' % iv
    Flow(time,[vint,exe],
         '''
         ./${SOURCES[1]} nr=5 r=285,509,648,728,906
         nh=24 dh=0.134 h0=0.264 type=hyperbolic
         ''')
    Plot(time+'g',time,
         '''
         graph yreverse=y pad=n min2=0 max2=3.996
         wantaxis=n wanttitle=n plotfat=10
         ''')
    Plot(time,['cmp',time+'g'],'Overlay')

    # RMS velocity
    vrms = 'vrms%d' % iv

    Flow(vrms,vint,
         '''
         add mode=p $SOURCE | causint | 
         math output="sqrt(input*0.004/(x1+0.004))" 
         ''')
    Plot(vrms+'w',vrms,
         '''
         graph yreverse=y transp=y pad=n 
         wanttitle=n wantaxis=n min2=1.4 max2=2.5 
         plotcol=7 plotfat=15
         ''')
    Plot(vrms+'b',vrms,
         '''
         graph yreverse=y transp=y pad=n 
         wanttitle=n wantaxis=n min2=1.4 max2=2.5 
         plotcol=0 plotfat=3
         ''')
    Plot(vrms,['vscan',vrms+'w',vrms+'b'],'Overlay')

    # Normal moveout
    nmo = 'nmo%d' % iv

    Flow(nmo,['cmp',vrms],'nmo velocity=${SOURCES[1]} half=n')
    Plot(nmo,
         '''
         grey title="Normal Moveout"
         grid2=y gridcol=6 gridfat=10
         ''')

    # Display it together
    allp = 'cmp%d' % iv
    Plot(allp,[time,vint,vrms,nmo],
         'SideBySideAniso',vppen='txscale=1.5')

    cmps.append(allp)
Plot('cmps',cmps,'Movie',view=1)

###############################
frame = 15
###############################
Result('cmp','cmp%d' % frame,'Overlay')

Flow('time',['vint%d' % frame,exe],
     '''
     ./${SOURCES[1]} nr=1 r=500
     nh=1001 dh=0.01 h0=0 type=hyperbolic
     ''')
Result('time',
       '''
       graph title=Traveltime
       label2=Time unit2=s yreverse=y
       label1=Offset unit1=km
       ''')

End()

Homework 2