Computational part

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 Figure 1. 2-D synthetic data.

model Figure 2. (a) Synthetic model: curved reflectors in a velocity.


mig Figure 3. Initial constant-velocity migration.

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

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 . 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.

Change directory
```
cd hw4/synth
```
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.
Run
```
pscons velcon.vpl
```
to display a movie of the velocity continuation process.
Run
```
pscons semb.vpl
```
to display a movie slicing through a semblance cube computed from velocity continuation.
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()

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.


data Figure 6. Zero-offset synthetic data.


dif Figure 7. Diffractions extracted from the data by plane-wave destruction.


dimage Figure 8. Time-migrated image of diffractions.


anglemig Figure 9. Dip angle gathers from constant-velocity angle-domain migration.

Change directory
```
cd hw4/sigsbee
```
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.
Generate a movie displaying the velocity continuation process. Is it possible to detect velocities from focusing zero-offset diffractions?
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?
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);
}

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 Figure 10. 2-D field dataset from the Gulf of Mexico after DMO stack.

Change directory
```
cd hw4/gulf
```
Run
```
scons view
```
to generate Figure 10 and display it on your screen.
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()