Homework 3

Computational part

In the computational part, we begin working with a real data set. The dataset is a 2-D line from the Blake Outer Ridge area offshore Florida. It was collected by USGS in order to study the occurrence of methane hydrates. The dataset and its analysis for gas hydrate detection are described by Ecker et al. (2000,1998).

Figure 1 shows an example CMP (common midpoint) gather from the dataset. Note the change in the trace spacing caused by using a non-linear cable in the acquisition. The presence of gas hydrates is manifested by a so-called BSR (bottom-simulating reflector). You can get an idea of the spatial extent of BSR from the near-offset section in Figure 2.

cmp
Figure 1. CMP (common midpoint) gather from the Blake Outer Ridge dataset. The presence of gas hydrates is manifested by BSR (bottom-simulating reflector).

noff Figure 2. Near offset section of the Blake Outer Ridge data.

1. We will use a very simple model to try explaining the geometry of the observed events in the data. The left plot in Figure 3 shows the model and rays from two-point ray tracing. The model consists of a constant velocity water layer and a sedimentary layer that contains BSR with a slowly variable velocity. The right plot in Figure 3 shows the CMP gather and traveltime curves obtained by ray tracing.

   

   ray
   Figure 3. Left: A two-layer model and ray trajectories. Right: A CMP gather and traveltime curves.

   

   
   

   Your task is to modify the model so that the predicted traveltime curves match the geometry of the sea bottom and the BSR reflection observed in the data.

   1. Change directory

      > cd ~/geo384w/hw3/blake
      

   2. Run

      > scons view
      

      to generate figures and display them on your screen.
   3. Edit the top of the SConstruct file to modify the model parameters. Check your result by running

      > scons ray.view
      

   #############################################################
   #
   # Adjustable parameters (modify for the assignment)
   #
   #############################################################
   
   velw=1  # water velocity (km/s)
   vels=2  # sediment velocity (km/s)
   water=2 # water depth (km)
   bsr=1   # bsr depth (km)
   
   #############################################################
   #
   # You don't need to change anything below 
   #
   #############################################################
   
   from rsfproj import *
   
   #############################################################
   #
   # Get data
   #
   #############################################################
   
   Fetch('cmps-tp.HH','blake')
   
   #############################################################
   #
   # Plot CMP Gather
   #
   #############################################################
   
   Flow('cmp','cmps-tp.HH',
        '''
        window f3=950 n3=1 max1=6 |
        dd form=native |
        reverse which=2 |
        put o2=0.0 d2=1
        ''')
   Plot('cmp',
        '''
        grey title="CMP Gather" wheretitle=t label1=Time unit1=s
        label2="Trace Number" wherexlabel=b max1=6
        ''')
   Plot('floorbox',None,
        'box x0=6.153 y0=6.565 label="Seafloor" xt=1 yt=1')
   Plot('bsrbox',None,
        'box x0=4.722 y0=4.778 label="BSR" xt=-1 yt=-1')
   Plot('gcmp','cmp floorbox bsrbox','Overlay')
   
   Flow('off1',None,'math n1=24 o1=0.4 d1=0.1  output=x1')
   Flow('off2',None,'math n1=24 o1=2.8 d1=0.05 output=x1')
   Flow('off','off1 off2','cat axis=1 ${SOURCES[1]}')
   Plot('wcmp','cmp off',
        '''
        wiggle xpos=${SOURCES[1]} transp=y yreverse=y 
        label1=Time unit1=s label2=Offset unit2=km
        poly=y title="CMP Gather"
        ''')
   
   Result('cmp','wcmp gcmp','SideBySideAniso')
   
   #############################################################
   #
   # Create a model
   #
   #############################################################
   
   depth = {'floor': water,
            'bsr': water+bsr}
   
   for ref in depth.keys():
       Flow(ref,None,
            'spike n1=501 d1=0.005 o1=-0.25 mag=%g' % depth[ref])
       Plot(ref,
            '''
            graph min2=0 max2=5 yreverse=y plotcol=0 plotfat=10 
            wanttitle=n wantaxis=n pad=n scalebar=y
            ''')
       
   Flow('vel','floor',
        '''
        unif2 d1=0.005 n1=1001 v00=%g,%g dvdz=0,0.1
        ''' % (velw,vels))
   Plot('vel',
        '''
        grey color=j bias=%g scalebar=y wanttitle=n
        label1=Depth unit1=km label2=Lateral unit2=km
        barlabel=Velocity barunit="km/s" barreverse=y
        wherexlabel=b
        ''' % (0.5*(velw+vels)))
   
   #############################################################
   #
   # Ray tracing
   #
   #############################################################
   
   rays=Split('vel floor bsr')
   times=['cmp',]
   for ref in depth.keys():
       ray=ref+'ray'
       for case in range(2):
           Flow(ray+str(case),'vel',
                '''
                shoot2 yshot=0 tol=0.001
                zshot=%g nr=24 r0=%g dr=%g
                ''' % ((depth[ref],0.2,0.05),
                       (depth[ref],1.4,0.025))[case])
       Flow(ray,[ray+'0',ray+'1'],'cat axis=1 ${SOURCES[1]}')
       Flow(['t'+ray,ray+'.ray'],['vel',ray],
            '''
            cell2 yshot=0 zshot=%g
            anglefile=${SOURCES[1]} traj=${TARGETS[1]}
            ''' % depth[ref])
       Plot(ray,[ray+'.ray','vel'],
            '''
            plotrays frame=${SOURCES[1]} plotcol=7 scalebar=y
            ''')
       rays.append(ray)
       Plot('t'+ray,
            '''
            scale dscale=2 |
            graph yreverse=y min2=4 max2=6 plotfat=5
            wantaxis=n wanttitle=n
            ''')
       times.append('t'+ray)
       
   Plot('ray',rays,'Overlay')
   Plot('time',times,'Overlay')
   
   Result('ray','ray time','SideBySideAniso')
   
   End()
   

2. In the second part of the computational assignment, you will use the water velocity that you obtained in the previous part to produce a seismic image of the near-offset section from Figure 2. We will use the zero-offset assumption that states that the reflector surface is obtained by backward extrapolating the recorded zero-offset reflection time into the subsurface using a velocity that is half of the true velocity and stopping at zero time. Using the intermediate time steps, you will create a movie of the backward traveltime extrapolation.
   1. Change directory

      > cd ~/geo384w/hw3/image
      

   2. Run

      > scons view
      

      to generate figures and display them on your screen.
   3. Edit the SConstruct file to change the water velocity to the one you found in the previous part. Check your result by running

      > scons image.view
      

   4. Modify the last part of the SConstruct file to generate a movie of wavefront extrapolation from the surface to the seafloor.

   from rsfproj import *
   
   #############################################################
   #
   # Get data
   #
   #############################################################
   
   Fetch('cmps-tp.HH','blake')
   
   #############################################################
   #
   # Plot near offset section
   #
   #############################################################
   
   Flow('noff','cmps-tp.HH','window f2=47 | dd form=native')
   Plot('noff',
        '''
        grey title="Near Offset Section" 
        label1=Time unit1=s label2=Distance unit2=km
        wheretitle=t wherexlabel=b
        ''')
   Plot('floor',None,
        'box x0=6.238 y0=6.705 label="Seafloor" xt=-1 yt=1')
   Plot('bsr',None,
        'box x0=10.918 y0=5.922 label="BSR" xt=1 yt=1')
   
   #############################################################
   #
   # Pick seafloor reflection
   #
   #############################################################
   
   Flow('sfr','noff',
        '''
        math output="input*input" | scale axis=2 |
        window max1=5.5 |
        transp | pick rect1=30 vel0=5.148 | window
        ''')
   Plot('sfr',
        '''
        graph yreverse=y min2=4 max2=6.496 pad=n
        wantaxis=n wanttitle=n plotcol=3 plotfat=5
        ''')
   
   Result('noff','noff floor bsr sfr','Overlay')
   
   #############################################################
   #
   # Water layer
   #
   # MODIFY BELOW !!!
   #
   #############################################################
   
   velw=1  # water velocity (km/s)
   
   Flow('water','noff',
        '''
        put d1=0.01 o1=0 |
        window max1=%g |
        math output=%g
        ''' % (velw*3,velw*0.5))
   
   #############################################################
   #
   # Backward time extrapolation
   #
   #############################################################
   
   Flow('time','water sfr','timecont surf=${SOURCES[1]}')
   Plot('time',
        '''
        contour allpos=y nc=1 c0=0
        wanttitle=n wantaxis=n plotfat=5
        ''')
   
   #############################################################
   #
   # Wave extrapolation
   #
   #############################################################
   
   Flow('water2','water','scale dscale=2')
   Flow('image','noff water2',
        '''
        cosft sign2=1 | stolt vel=%g | cosft sign2=-1 |
        time2depth velocity=${SOURCES[1]}
        ''' % velw)
   Plot('image',
        '''
        grey label1=Depth unit1=km label2=Distance unit2=km
        title=Image
        ''')
   
   Result('image','image time','Overlay')
   
   #############################################################
   #
   # Movie
   #
   #############################################################
   
   ts = range(5)
   ts.reverse()
   times = []
   
   for t in ts:
       time = 'time%d' % t
       times.append(time)
   
       ##########################
       #
       # ADD COMMANDS HERE
       #
       ##########################
   
   Plot('times',times,'Movie')
   
   End()
   

Homework 3