Homework 4

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

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.

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.

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

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

 Homework 4

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

2013-10-29