next up previous [pdf]

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

Computational part

  1. In the first part of the computational assignment, you will experiment with imaging a synthetic seismic reflection dataset from Homework 3 using prestack velocity continuation.

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

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

    mig
    mig
    Figure 3.
    Initial constant-velocity migration.
    [pdf] [png] [scons]

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

    Figure 6 shows a synthetic reflection dataset computed from a reflector model shown in Figure 2 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 data.

    Figure 3 shows an initial prestack common-offset time migration using a constant velocity of 1.5 km/s. Figure 4 shows the result of prestack time migration after velocity continuation, extraction of a velocity slice, and conversion from time to depth.

    1. Change directory
      cd hw4/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. Run
      pscons velcon.vpl
      
      to display a movie of the velocity continuation process.
    4. Run
      pscons semb.vpl
      
      to display a movie slicing through a semblance cube computed from velocity continuation.
    5. The processing flow in the SConstruct file involves some cheating: the exact RMS velocity is used to extract the final image. Modify the processing flow so that only properties estimated from the data get used.

    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.01 n1=1001'
    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.25*x1"
         d2=0.05 n2=201 d1=0.01 n1=501
         label1=Depth unit1=km
         label2=Distance unit2=km
         label=Velocity unit=km/s
         ''')
    Plot('vofz',
         '''
         window min2=2.75 max2=7.25 |
         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 dscale=100')
    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.25 verb=y |
         tpow tpow=1 |
         put d2=0.05 label3=Midpoint unit3=km 
         ''',split=[1,1001], reduce='add')
    
    # Add random noise
    Flow('data','modl','noise var=1e-6 seed=101811')
    
    Result('data',
           '''
           byte |
           transp plane=23 |
           grey3 flat=n frame1=750 frame2=100 frame3=10 
           label1=Time unit1=s 
           label3=Half-Offset unit3=km 
           title=Data point1=0.8 point2=0.8
           ''')
    
    # Initial constant-velocity migration
    #####################################
    Flow('mig','data',
         '''
         transp plane=23 |
         spray axis=3 n=1 d=0.1 o=0 |
         preconstkirch vel=1.5 |
         halfint inv=1 adj=1
         ''',split=[2,51],reduce='cat axis=4')
    
    Result('mig',
           '''
           byte | window |
           grey3 flat=n frame1=750 frame2=100 frame3=10 
           label1=Time unit1=s 
           label3=Half-Offset unit3=km 
           title="Initial Migration" point1=0.8 point2=0.8
           ''')
    
    # Velocity continuation
    #######################
    
    Flow('thk','mig',
         'window | transp plane=23 | cosft sign3=1')
    Flow('velconk','thk',
         'fourvc nv=81 dv=0.01 v0=1.5 verb=y',
         split=[3,201])
    Flow('velcon','velconk','cosft sign3=-1')
    
    Plot('velcon',
         '''
         transp plane=23 memsize=1000 |
         window min2=2.5 max2=7.5 |
         grey title="Velocity Continuation"
         ''',view=1)
    
    # Continue data squared
    Flow('thk2','mig',
         '''
         mul $SOURCE |
         window | transp plane=23 | cosft sign3=1
         ''')
    Flow('velconk2','thk2',
         'fourvc nv=81 dv=0.01 v0=1.5 verb=y',
         split=[3,201])
    Flow('velcon2','velconk2','cosft sign3=-1')
    
    # Compute semblance
    Flow('semb','velcon velcon2',
         '''
         mul $SOURCE | divn den=${SOURCES[1]} rect1=25
         ''',split=[3,201])
    
    Plot('semb',
         '''
         byte gainpanel=all allpos=y |
         transp plane=23 |
         grey3 flat=n frame1=750 frame2=0 frame3=48 
         label1=Time unit1=s color=j 
         label3=Velocity unit3=km/s movie=2 dframe=5
         title=Semblance point1=0.8 point2=0.8
         ''',view=1)
    
    # Extracting images
    ###################
    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))" 
         ''')
    
    # Using vrms is CHEATING
    ########################
    Flow('slice','velcon vrms','slice pick=${SOURCES[1]}')
    
    # Using vofz is CHEATING
    ########################
    Flow('dmig','slice 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. In the second part of the computational assignment, we will use velocity continuation again but this time on a synthetic zero-offset section containing diffraction events.

    Figure 5 shows a famous Sigsbee synthetic velocity model. We will focus on the left part of the model, which is appropriate for time-domain imaging. A synthetically generated zero-offset section is shown in Figure 6.

    Our processing strategy is to extract diffractions from the data (Figure 7) and to image them using zero-offset velocity continuation (Figure 8). In addition, we are going to analyze the image by expanding it in dip angles by using dip-angle migration (Figure 9).

    vel
    Figure 5.
    Sigsbee velocity model.
    vel
    [pdf] [png] [scons]

    data
    data
    Figure 6.
    Zero-offset synthetic data.
    [pdf] [png] [scons]

    dif
    dif
    Figure 7.
    Diffractions extracted from the data by plane-wave destruction.
    [pdf] [png] [scons]

    dimage
    dimage
    Figure 8.
    Time-migrated image of diffractions.
    [pdf] [png] [scons]

    anglemig
    anglemig
    Figure 9.
    Dip angle gathers from constant-velocity angle-domain migration.
    [pdf] [png] [scons]

    1. Change directory
      cd hw4/sigsbee
      
    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. Generate a movie displaying the velocity continuation process. Is it possible to detect velocities from focusing zero-offset diffractions?
    4. Modify the program in the anglemig.c file to input a variable migration velocity instead of using a constant velocity. Regenerate Figure 9 using a variable velocity
      pscons anglemig.view
      
      Do you notice a difference?
    5. For EXTRA CREDIT, implement a method for estimating migration velocity from the input data.

    from rsf.proj import *
    
    # Download velocity model from the data server
    ##############################################
    vstr = 'sigsbee2a_stratigraphy.sgy'
    Fetch(vstr,'sigsbee')
    Flow('zvstr',vstr,'segyread read=data')
    
    Flow('vel','zvstr',
         '''
         put d1=0.00762 o2=3.048 d2=0.00762
         label1=Depth unit1=km label2=Distance unit2=km
         label=Velocity unit=km/s |
         scale dscale=0.0003048
         ''')
    Result('vel',
           '''
           grey wanttitle=n color=j allpos=y 
           screenratio=0.3125 screenht=4 labelsz=4
           scalebar=y barreverse=y
           ''')
    
    # Window a portion
    Flow('vel2','vel','window max2=9.5')
    
    dt = 0.002
    nt = 5001
    
    # Convert to RMS
    Flow('voft','vel2',
         'depth2time velocity=$SOURCE dt=%g nt=%d' % (dt,nt))
    Flow('vrms','voft',
         '''
         add mode=p $SOURCE | causint | 
         math output="sqrt(input*%g/(x1+%g))" |
         smooth rect2=10 | window j2=2
         ''' % (dt,dt))
    
    # Download zero-offset from the data server
    ###########################################
    Fetch('sigexp_ns.rsf','sigsbee')
    Flow('data','sigexp_ns.rsf',
         '''
         dd form=native |
         bandpass flo=2 fhi=60 |
         window max2=9.5 j1=2 j2=2 n1=%d |
         costaper nw1=50 nw2=50
         ''' % nt)
    
    Result('data',
           'window min1=3 | grey title="Zero-Offset Data" ')
    
    # Slope estimation
    Flow('dip','data','fdip rect1=100 rect2=10')
    Result('dip',
           '''
           grey color=j scalebar=y
           title="Dominant Slope"
           barlabel=Slope barunit=samples
           ''')
    
    # Plane-wave destruction
    Flow('dif','data dip','pwd dip=${SOURCES[1]}')
    Result('dif',
           '''
           window min1=3 | 
           grey title="Separated Diffractions" 
           ''')
    
    # Velocity continuation
    Flow('fourier','dif','pad n2=1025 | cosft sign2=1')
    Flow('velconf','fourier',
         '''
         put o3=0 | stolt vel=1.4 | 
         spray axis=2 n=1 o=0 d=1 | 
         fourvc pad2=8192 nv=61 dv=0.02 v0=1.4 verb=y
         ''',split=[2,1025],reduce='cat axis=3')
    Flow('velcon','velconf',
         '''
         transp plane=23 memsize=1000 | 
         cosft sign2=-1  | window n2=424 | 
         transp plane=23 memsize=1000
         ''')
    
    # Picking a slice
    #################
    Flow('dimage','velcon vrms',
         'slice pick=${SOURCES[1]}')
    Result('dimage',
           '''
           window min1=3 |
           grey title="Imaged Diffractions"
           ''')
    
    # Angle-gather migration
    ########################
    proj = Project()
    prog = proj.Program('anglemig.c')
    
    Flow('anglemig','dif %s' % prog[0],
         './${SOURCES[1]} vel=2 na=90 a0=-45 da=1')
    
    Result('anglemig',
           '''
           window min2=2 |
           transp | transp plane=23 memsize=1000 | 
           byte gainpanel=all | grey3
           frame1=1000 frame2=200 frame3=60 unit3="ô_"
           title="Dip Angle Gathers" point1=0.7 point2=0.7
           ''')
    
    End()
    

    /* 2-D angle-domain zero-offset migration. */
    #include <rsf.h>
    
    static float get_sample (float **dat, 
    			 float t, float y,
    			 float t0, float y0, 
    			 float dt, float dy, 
    			 int nt, int ny) 
    /* extract data sample by linear interpolation */
    {
        int it, iy;
    
        y = (y - y0)/dy; iy = floorf (y);
        y -= (float)iy;
        if (iy < 0 || iy >= (ny - 1)) return 0.0;
        t = (t - t0)/dt; it = floorf (t);
        t -= (float)it;
        if (it < 0 || it >= (nt - 1)) return 0.0;
    
        return (dat[iy][it]*(1.0 - y)*(1.0 - t) +
    	    dat[iy][it + 1]*(1.0 - y)*t +
    	    dat[iy + 1][it]*y*(1.0 - t) +
    	    dat[iy + 1][it + 1]*y*t);
    }
    
    int main (int argc, char* argv[]) 
    {
        int it, nt, ix, nx, ia, na;
        float dt, vel, da, a0, dx, z, t, x, y, a;
        float **dat, *img;
        sf_file data, imag;
    
        sf_init (argc, argv);
    
        data = sf_input ("in");
        imag = sf_output ("out");
    
        /* get dimensions */
        if (!sf_histint (data, "n1", &nt))   sf_error ("n1");
        if (!sf_histint (data, "n2", &nx))   sf_error ("n2");
        if (!sf_histfloat (data, "d1", &dt)) sf_error ("d1");
        if (!sf_histfloat (data, "d2", &dx)) sf_error ("d2");
    
        if (!sf_getint("na",&na)) sf_error("Need na="); 
        /* number of angles */
        if (!sf_getfloat("da",&da)) sf_error("Need da="); 
        /* angle increment */
        if (!sf_getfloat("a0",&a0)) sf_error("Need a0="); 
        /* initial angle */
    
        sf_shiftdim(data, imag, 1);
    
        sf_putint(imag,"n1",na);
        sf_putfloat(imag,"d1",da);
        sf_putfloat(imag,"o1",a0);
        sf_putstring(imag,"label1","Angle");
    
        /* degrees to radians */
        a0 *= SF_PI/180.;
        da *= SF_PI/180.;
    
        if (!sf_getfloat("vel",&vel)) vel=1.5; 
        /* constant velocity */
    
        dat = sf_floatalloc2(nt,nx);
        sf_floatread (dat[0],nt*nx, data);
    
        img = sf_floatalloc (na);
     
        /* loop over image location */
        for (ix = 0; ix < nx; ix++) {
    	x = ix*dx;
            sf_warning ("CMP %d of %d;", ix, nx);
       
    	/* loop over image time */
    	for (it = 0; it < nt; it++) {
    	    z = it*dt;
    
    	    /* loop over angle */
                for (ia = 0; ia < na; ia++) { 
    		a = a0+ia*da;
    		
    		t = z/cosf(a);           
                    /* escape time */
    		y = x+0.5*vel*t*sinf(a); 
                    /* escape location */
    
    		img[ia] = get_sample (dat,t,y,0.,0.,
    				      dt,dx,nt,nx);
    	    } /* ia */
    
                sf_floatwrite (img, na, imag);
            } /* it */
        } /* ix */
     
        exit(0);
    }
    
    

  3. In the final part of the computational assignment, we return to the 2-D field dataset from the Gulf of Mexico. The zero-offset data after a DMO stack are shown in Figure 10.

    stack
    stack
    Figure 10.
    2-D field dataset from the Gulf of Mexico after DMO stack.
    [pdf] [png] [scons]

    1. Change directory
      cd hw4/gulf
      
    2. Run
      scons view
      
      to generate Figure 10 and display it on your screen.
    3. Edit the SConstruct file to implement a processing flow involving velocity continuation and angle-gather migration. Make sure to select appropriate processing parameters.

    from rsf.proj import *
    
    Fetch('bei-stack.rsf','midpts')
    Flow('stack','bei-stack',
         '''
         dd form=native |
         put label2=Distance unit2=km label1=Time unit1=s
         ''')
    
    Result('stack','grey title="DMO Stack" ')
    
    End()
    


next up previous [pdf]

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

2013-10-29