Homework 3

Computational part

1. In the first part, you will experiment with creating and imaging a synthetic seismic reflection dataset.

   

   data
   Figure 2. 2-D synthetic data.

   

   
   

   model Figure 3. (a) Synthetic model: curved reflectors in a $V(z)$ velocity.

   

   vscan
   Figure 4. Velocity semblance scan.

   

   
   

   

   vnmo
   Figure 5. RMS velocity (a) and picked NMO velocity (b).

   

   
   

   

   dstack,zoff
   Figure 6. (a) DMO stack. (b) Zero-offset section.

   

   
   

   

   tmig
   Figure 7. Kirchhoff poststack time migration.

   

   
   

   dmig2 Figure 8. Time migration converted to depth, with reflectors overlaid.

   Figure 2 shows a synthetic reflection dataset computed from a reflector model shown in Figure 3 and assuming a velocity model with a constant vertical gradient $V(z) = 1.5 + 0.36\,z$. A small amount of random noise is added to the synthetic data.

   Figure 4 shows a conventional velocity semblance scan. Figure 5 compares the RMS velocity with the NMO velocity picked from the scan. Figure 6 compares a DMO stack and the zero-offset section from the data. Finally, Figures 7 and 8 show the result of Kirchhoff time migration before and after conversion from time to depth.

   1. Change directory

      cd hw3/synth
      

   2. Run

      scons view
      

      to generate the figures and display them on your screen. If you are on a computer with multiple CPUs, you can also try

      pscons view
      

      to run certain computations in parallel.
   3. Explain the cause of the difference between the RMS velocity and the NMO velocity in Figure 5.
   4. Edit the SConstruct file to switch on the antialiasing parameter in Kirchhoff migration (by changing it from 0 to 1). Generate the migration figures again. What differences do you observe?
   5. Edit the kirchhoff.c file to improve the result of migration. This task is open-ended. The change is up to you, as long as you can achieve an improvement. Possible options:
      • Use a more accurate traveltime computation.
      • Introduce an amplitude weight.
      • Change from time to depth migration (you can assume a locally $V(z)$ medium and use results from Homework 2.)
      • Change from poststack to prestack migration.
   6. The processing flow in the SConstruct file involves some cheating: the exact depth velocity and the exact RMS velocity are used without being estimated from the data. Modify the processing flow so that only properties estimated from the data get used. This task is open-ended as well, different data processing strategies are possible.

   from rsf.proj import *
   
   # Generate a reflector model
   
   layers = (
        ((0,2),(3.5,2),(4.5,2.5),(5,2.25),
          (5.5,2),(6.5,2.5),(10,2.5)),
        ((0,2.5),(10,3.5)),
        ((0,3.2),(3.5,3.2),(5,3.7),
          (6.5,4.2),(10,4.2)),
        ((0,4.5),(10,4.5))
   )
   
   nlays = len(layers)
   for i in range(nlays):
        inp = 'inp%d' % i
        Flow(inp+'.asc',None,
             '''
             echo %s in=$TARGET
             data_format=ascii_float n1=2 n2=%d
             ''' %            (' '.join(map(lambda x: ' '.join(map(str,x)),
                              layers[i])),len(layers[i])))
   
   dim1 = 'o1=0 d1=0.001 n1=10001'
   
   Flow('lay1','inp0.asc',
        'dd form=native | spline %s fp=0,0' % dim1)
   Flow('lay2',None  ,
        'math %s output="2.5+x1*0.1" '      % dim1)
   Flow('lay3','inp2.asc',
        'dd form=native | spline %s fp=0,0' % dim1)
   Flow('lay4',None  ,'math %s output=4.5'  % dim1)
   
   Flow('lays','lay1 lay2 lay3 lay4',
        'cat axis=2 ${SOURCES[1:4]}')
   
   graph = '''
   graph min1=2.5 max1=7.5 min2=0 max2=5
   yreverse=y wantaxis=n wanttitle=n screenratio=1
   '''
   Plot('lays0','lays',graph + ' plotfat=10 plotcol=0')
   Plot('lays1','lays',graph + ' plotfat=2 plotcol=7')
   Plot('lays2','lays',graph + ' plotfat=2')
   
   # Velocity
   
   Flow('vofz',None,
        '''
        math output="1.5+0.36*x1"
        d2=0.01 n2=1001 d1=0.01 n1=501
        label1=Depth unit1=km
        label2=Distance unit2=km
        label=Velocity unit=km/s
        ''')
   Plot('vofz',
        '''
        window min2=2.5 max2=7.5 |
        grey color=j allpos=y bias=1.5
        title=Model screenratio=1
        ''')
   
   Result('model','vofz lays0 lays1','Overlay')
   
   # Model data
   
   Flow('dips','lays','deriv scale=y')
   Flow('modl','lays dips',
        '''
        kirmod cmp=y dip=${SOURCES[1]} 
        nh=51  dh=0.1  h0=0
        ns=201 ds=0.05 s0=0
        freq=10 dt=0.004 nt=1501
        vel=1.5 gradz=0.36 verb=y |
        pow pow1=1 |
        put d2=0.05 
        label1=Time        unit1=s
        label2=Half-Offset unit2=km 
        label3=Midpoint    unit3=km 
        ''',split=[1,10001], reduce='add')
   
   # Add random noise
   Flow('data','modl','noise var=1e-6 seed=101013')
   
   Result('data',
          '''
          byte |
          transp plane=23 |
          grey3 flat=n frame1=750 frame2=100 frame3=10 
          title=Data point1=0.8 point2=0.8
          ''')
   
   # Velocity estimation
   #####################
   Flow('voft','vofz',
        'depth2time velocity=$SOURCE dt=0.004 nt=1501')
   Flow('vrms','voft',
        '''
        add mode=p $SOURCE | causint | 
        math output="sqrt(input*0.004/(x1+0.004))" 
        ''')
   
   # Velocity scan
   Flow('vscan','data',
        'vscan v0=1.5 dv=0.02 nv=51 semblance=y',
        split=[3,201], reduce='cat')
   
   Result('vscan',
          '''
          byte allpos=y gainpanel=all |
          transp plane=23 |
          grey3 flat=n frame1=750 frame2=100 frame3=25 
          label1=Time unit1=s color=j
          label3=Velocity unit3=km/s 
          label2=Midpoint unit2=km
          title="Velocity Scan" point1=0.8 point2=0.8
          ''')
   
   # Velocity picking
   Flow('vnmo','vscan','pick rect1=100 rect2=10')
   
   for vel in ('vrms','vnmo'):
       Plot(vel,
        '''
        grey color=j allpos=y bias=1.5 clip=0.7
        scalebar=y barreverse=y barunit=km/s
        label2=Midpoint unit2=km label1=Time unit1=s
        title="%s Velocity" 
        ''' % vel[1:].upper())
   Result('vnmo','vrms vnmo','SideBySideIso')
   
   # Stacking
   ##########
   
   Flow('nmo','data vnmo','nmo velocity=${SOURCES[1]}')
   Flow('stack','nmo','stack')
   
   # Using vrms is CHEATING
   ########################
   Flow('nmo0','data vrms','nmo velocity=${SOURCES[1]}')
   Flow('dstack','nmo0',
        '''
        window f1=250 | 
        logstretch | fft1 | 
        transp plane=13 memsize=1000 |
        finstack | 
        transp memsize=1000 |
        fft1 inv=y | logstretch inv=y | 
        pad beg1=250 | put unit1=s
        ''')
   
   Flow('zoff','data','window n2=1')
   
   stacks = {
       'stack': 'Stack with NMO Velocity',
       'dstack': 'DMO Stack',
       'zoff': 'Zero Offset'
       }
   
   for stack in stacks.keys():
       Result(stack,
              '''
              window min1=1.5 max1=5 min2=1 max2=9 | 
              grey title="%s" 
              ''' % stacks[stack])
   
   # Kirchhoff Migration
   #####################
   
   proj = Project()
   
   prog = proj.Program('kirchhoff',['kirchhoff.c','aal.c'])
   exe = str(prog[0])
   
   # Using vrms is CHEATING
   ########################
   Flow('tmig','dstack %s vrms' % prog[0],
        './${SOURCES[1]} vel=${SOURCES[2]} antialias=0')
   
   Result('tmig',
          '''
          window min2=2.5 max2=7.5 |
          grey title="Time Migration"
          ''')
   
   # Using vofz is CHEATING
   ########################
   Flow('dmig','tmig vofz',
        'time2depth velocity=${SOURCES[1]}')
   
   Plot('dmig',
        '''
        window max1=5 min2=2.5 max2=7.5 | 
        grey title="Time -> Depth" screenratio=1
        label2=Distance label1=Depth unit1=km
        ''')
   
   Result('dmig','Overlay')
   Result('dmig2','dmig lays2','Overlay')
   
   End()
   

   /* 2-D Post-stack Kirchhoff time migration. */
   #include <rsf.h>
   
   #include "aal.h" /* antialising routines */
   
   int main(int argc, char* argv[])
   {
       int nt, nx, nz, ix, iz, iy, i;
       float *trace, **out, **v;
       float x,z, dx, ti, tx, t0,dt, z0,dz, vi,aal;
       sf_file inp, mig, vel;
   
       sf_init (argc,argv);
       inp = sf_input("in");
       vel = sf_input("vel");
       mig = sf_output("out");
   
       if (!sf_histint(inp,"n1",&nt)) sf_error("No n1=");
       if (!sf_histint(inp,"n2",&nx)) sf_error("No n2=");
   
       if (!sf_histfloat(inp,"o1",&t0)) sf_error("No o1=");
       if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d1=");
       if (!sf_histfloat(inp,"d2",&dx)) sf_error("No d2=");
   
       if (!sf_getint("nz",&nz)) nz=nt;
       if (!sf_getfloat("dz",&dz)) dz=dt;
       if (!sf_getfloat("z0",&z0)) z0=t0;
   
       if (!sf_getfloat("antialias",&aal)) aal=1.0;
       /* antialiasing */
   
       v = sf_floatalloc2(nz,nx);
       sf_floatread(v[0],nz*nx,vel);
   
       trace = sf_floatalloc(nt);
       out = sf_floatalloc2(nz,nx);
   
       for (i=0; i < nz*nx; i++) {
   	out[0][i] = 0.;
       }
   
       /* loop over input traces */
       for (iy=0; iy < nx; iy++) { 
   	sf_floatread (trace,nt,inp);
   	aal_doubint(nt,trace);
           
   	/* loop over output traces */
   	for (ix=0; ix < nx; ix++) { 
   	    x = (ix-iy)*dx;
   
   	    /* loop over output time */
   	    for (iz=0; iz < nz; iz++) {
   		z = z0 + iz*dz;  
   		vi = v[ix][iz];
   
   		/* hypot(a,b) = sqrt(a*a+b*b) */
   		ti = hypotf(z,2.0*x/vi);		
   
   		/* tx = |dt/dx| */
   		tx = 4.0*fabsf(x)/(vi*vi*(ti+dt));
   
   		out[ix][iz] += 
   		    aal_pick(ti,tx*dx*aal,trace,nt,dt,t0);
   	    } 
   	} 
       } 
       
       sf_floatwrite(out[0],nz*nx,mig);        
   
       exit(0);
   }
   

   

   gulf
   Figure 9. 2-D field dataset from the Gulf of Mexico.

   

   
   

2. In the second part of the computational assignment, you will use the processing strategy developed for synthetic data to process a 2-D field dataset from the Gulf of Mexico, shown in Figure 9.
   1. Change directory

      cd hw3/gulf
      

   2. Run

      scons view
      

      to generate the figures and display them on your screen.
   3. Edit the SConstruct file to implement your processing strategy and to generate a seismic image. Make sure to select appropriate processing parameters.

   from rsf.proj import *
   
   # Download data
   Fetch('beinew.HH','midpts')
   
   # Set dimensions
   Flow('gulf','beinew.HH',
        '''
        dd form=native |
        put
        label1=Time unit1=s
        label2=Half-Offset unit2=km
        label3=Midpoint unit3=km
        ''')
   
   # Display
   Result('gulf',
          '''
          byte |
          transp plane=23 |
          grey3 flat=n frame1=500 frame2=160 frame3=10 
          title=Data point1=0.8 point2=0.8
          ''')
   
   End()
   

Homework 3