next up previous [pdf]

Next: Completing the assignment Up: Homework 3 Previous: Theoretical part

Computational part

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

    data
    data
    Figure 2.
    2-D synthetic data.
    [pdf] [png] [scons]

    model
    Figure 3.
    (a) Synthetic model: curved reflectors in a $V(z)$ velocity.
    model
    [pdf] [png] [scons]

    vscan
    vscan
    Figure 4.
    Velocity semblance scan.
    [pdf] [png] [scons]

    vnmo
    vnmo
    Figure 5.
    RMS velocity (a) and picked NMO velocity (b).
    [pdf] [png] [scons]

    dstack zoff
    dstack,zoff
    Figure 6.
    (a) DMO stack. (b) Zero-offset section.
    [pdf] [pdf] [png] [png] [scons]

    tmig
    tmig
    Figure 7.
    Kirchhoff poststack time migration.
    [pdf] [png] [scons]

    dmig2
    Figure 8.
    Time migration converted to depth, with reflectors overlaid.
    dmig2
    [pdf] [png] [scons]

    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.25 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
    gulf
    Figure 9.
    2-D field dataset from the Gulf of Mexico.
    [pdf] [png] [scons]

  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()
    


next up previous [pdf]

Next: Completing the assignment Up: Homework 3 Previous: Theoretical part

2013-10-10