Homework 6

Next: Completing the assignment Up: Homework 6 Previous: Homework 6

# Computational part

1. We start by returning to the model from Homeworks 1 and 5 with a hyperbolic reflector under a constant velocity layer. The model is shown in Figure 1. Now we will approach the imaging task using reverse-time migration with a two-way wave extrapolation. Figure 2a shows synthetic zero-offset data generated by Kirchhoff modeling. Figure 8 shows an image generated by zero-offset reverse-time migration using an explicit finite-difference wave extrapolation in time.

Your task: Change the program for reverse-time migration to implement forward-time modeling using the ``exploding reflector'' approach.

1. Change directory
```cd hw6/hyper
```
2. Run
```scons view
```
to generate figures and display them on your screen.
3. Run
```scons wave.vpl
```
to observe a movie of reverse-time wave extrapolation.
4. Edit the program in the rtm.c file to implement a process opposite to migration: starting from the reflectivity image like the one in Figure 8 and generating zero-offset data like the one in Figure 2a.
5. Run
```scons view
```
again to observe the differences.

model
Figure 1.
Synthetic velocity model with a hyperbolic reflector.

data,image
Figure 2.
(a) Synthetic zero-offset data corresponding to the model in Figure 1. (b) Image generated by reverse-time exploding-reflector migration. The location of the exact reflector is indicated by a curve.

 ```from rsf.proj import * # Make a reflector model Flow('refl',None, ''' math n1=10001 o1=-4 d1=0.001 output="sqrt(1+x1*x1)" ''') Flow('model','refl', ''' window min1=-3 max1=3 j1=5 | unif2 n1=401 d1=0.005 v00=1,2 | put label1=Depth unit1=km label2=Lateral unit2=km | smooth rect1=3 ''') Plot('model', ''' window min2=-1 max2=2 | grey allpos=y bias=1 title=Model ''') Plot('refl', ''' graph wanttitle=n wantaxis=n yreverse=y min1=-1 max1=2 min2=0 max2=2 plotfat=3 plotcol=4 ''') Result('model','model refl','Overlay') # Kirchoff zero-offset modeling Flow('dip','refl','math output="x1/input" ') Flow('data','refl dip', ''' kirmod nt=1501 dt=0.004 freq=10 ns=1201 s0=-3 ds=0.005 nh=1 dh=0.01 h0=0 twod=y vel=1 dip=\${SOURCES[1]} verb=y | window | costaper nw2=200 ''') Result('data', ''' window min2=-2 max2=2 min1=1 max1=4 | grey title="Zero-Offset Data" grid=y label1=Time unit1=s label2=Distance unit2=km ''') # Reverse-time migration proj = Project() prog = proj.Program('rtm.c') Flow('image wave','model %s data' % prog[0], ''' pad beg1=200 end1=200 | math output=0.5 | ./\${SOURCES[1]} data=\${SOURCES[2]} wave=\${TARGETS[1]} jt=20 n0=200 ''') # Display wave propagation movie Plot('wave', ''' window min2=-1 max2=2 min1=0 max1=2 f3=33 | grey wanttitle=n gainpanel=all ''',view=1) # Display image Plot('image', ''' window min2=-1 max2=2 min1=0 max1=2 | grey title=Image ''') Result('image','image refl','Overlay') End() ```

 ```/* 2-D zero-offset reverse-time migration */ #include #ifdef _OPENMP #include #endif static int n1, n2; static float c0, c11, c21, c12, c22; static void laplacian(float **uin /* [n2][n1] */, float **uout /* [n2][n1] */) /* Laplacian operator, 4th-order finite-difference */ { int i1, i2; #ifdef _OPENMP #pragma omp parallel for private(i2,i1) shared(n2,n1,uout,uin,c11,c12,c21,c22,c0) #endif for (i2=2; i2 < n2-2; i2++) { for (i1=2; i1 < n1-2; i1++) { uout[i2][i1] = c11*(uin[i2][i1-1]+uin[i2][i1+1]) + c12*(uin[i2][i1-2]+uin[i2][i1+2]) + c21*(uin[i2-1][i1]+uin[i2+1][i1]) + c22*(uin[i2-2][i1]+uin[i2+2][i1]) + c0*uin[i2][i1]; } } } int main(int argc, char* argv[]) { int it,i1,i2; /* index variables */ int nt,n12,n0,nx, jt; float dt,dx,dz, dt2,d1,d2; float **vv, **dd; float **u0,**u1,u2,**ud; /* tmp arrays */ sf_file data, imag, modl, wave; /* I/O files */ /* initialize Madagascar */ sf_init(argc,argv); /* initialize OpenMP support */ omp_init(); /* setup I/O files */ modl = sf_input ("in"); /* velocity model */ imag = sf_output("out"); /* output image */ data = sf_input ("data"); /* seismic data */ wave = sf_output("wave"); /* wavefield */ /* Dimensions */ if (!sf_histint(modl,"n1",&n1)) sf_error("n1"); if (!sf_histint(modl,"n2",&n2)) sf_error("n2"); if (!sf_histfloat(modl,"d1",&dz)) sf_error("d1"); if (!sf_histfloat(modl,"d2",&dx)) sf_error("d2"); if (!sf_histint (data,"n1",&nt)) sf_error("n1"); if (!sf_histfloat(data,"d1",&dt)) sf_error("d1"); if (!sf_histint(data,"n2",&nx) || nx != n2) sf_error("Need n2=%d in data",n2); n12 = n1*n2; if (!sf_getint("n0",&n0)) n0=0; /* surface */ if (!sf_getint("jt",&jt)) jt=1; /* time interval */ sf_putint(wave,"n3",1+(nt-1)/jt); sf_putfloat(wave,"d3",-jt*dt); sf_putfloat(wave,"o3",(nt-1)*dt); dt2 = dt*dt; /* set Laplacian coefficients */ d1 = 1.0/(dz*dz); d2 = 1.0/(dx*dx); c11 = 4.0*d1/3.0; c12= -d1/12.0; c21 = 4.0*d2/3.0; c22= -d2/12.0; c0 = -2.0 * (c11+c12+c21+c22); /* read data and velocity */ dd = sf_floatalloc2(nt,n2); sf_floatread(dd[0],nt*n2,data); vv = sf_floatalloc2(n1,n2); sf_floatread(vv[0],n12,modl); /* allocate temporary arrays */ u0=sf_floatalloc2(n1,n2); u1=sf_floatalloc2(n1,n2); ud=sf_floatalloc2(n1,n2); for (i2=0; i2= 0; it-) { sf_warning("%d;",it); laplacian(u1,ud); #ifdef _OPENMP #pragma omp parallel for private(i2,i1,u2) shared(ud,vv,it,u1,u0,dd,n0) #endif for (i2=0; i2

2. Figure 3 shows the Sigsbee velocity model (Paffenholz et al., 2002) and its an approximate filtered reflectivity.

zvel,zref
Figure 3.
(a) Sigsbee velocity model. (b) Approximate reflectivity of the Sigsbee model (an ideal image).

1. Change directory
```cd hw6/sigsbee
```
2. Run
```scons view
```
to generate figures and display them on your screen.
3. Modify the SConstruct file to implement the modeling and migration experiment.
4. Include your results in the paper by editing the hw6/paper.tex file.
5. For EXTRA CREDIT, make sure that your modeling and migration are adjoint to each other and apply least-squares exploding-reflector reverse-time migration to the Sigsbee model.

 ```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('zvel','zvstr', ''' put d1=0.025 o2=10.025 d2=0.025 label1=Depth unit1=kft label2=Lateral unit2=kft | scale dscale=0.001 ''') Result('zvel', ''' window f1=200 | grey title=Velocity titlesz=7 color=j screenratio=0.3125 screenht=4 labelsz=5 mean=y ''') # Compute approximate reflectivity ################################## Flow('zref','zvel', ''' depth2time velocity=\$SOURCE nt=2501 dt=0.004 | ai2refl | ricker1 frequency=10 | time2depth velocity=\$SOURCE ''') Result('zref', ''' window f1=200 | grey title=Reflectivity titlesz=7 screenratio=0.3125 screenht=4 labelsz=5 ''') End() ```

3. In the last part of the homework, you will work with a field dataset: a 2-D line from the Blake Outer Ridge area offshore Florida and Georgia (Figure 4). It was collected by USGS in order to study the occurrence of methane hydrates. The presence of gas hydrates is manifested by a so-called BSR (bottom-simulating reflector). The dataset and its analysis for gas hydrate detection are described by Ecker et al. (2000,1998).

map2
Figure 4.
Map of the Blake Outer Ridge area. An region of known gas hydrate distribution as mapped from seismic bottom simulating reflectors (BSR) is highlighted. The seismic survey area is marked by a rectangle.

The following figures show the dataset at different stages of seismic data processing: from initial data to an image in depth.

Your task: Modify the data processing sequence to create a justifiably better image.

nmo
Figure 5.
Common midpoint gather (left), velocity analysis panel using normal moveout (middle), and common midpoint gather after normal moveout (right). A curve in the middle plot indicates an automatically picked velocity trend.

noff,stack
Figure 6.
(a) Near-offset section. (b) Normal moveout stack.

picks,vel
Figure 7.
(a) Picked stacking velocity. (b) Interval velocity estimated by Dix inversion.

image
Figure 8.
Seismic image created by converting stacked data from time to depth. Can you identify geological features that are not properly imaged?

 ```from rsf.proj import * # get data Fetch('cmps-tp.HH','blake') # CMP (common midpoint) gathers Flow('cmps','cmps-tp.HH', 'dd form=native | reverse which=2') # one CMP Flow('cmp','cmps', ''' window f3=950 n3=1 max1=6 | put o2=0.0 d2=1 ''') Plot('cmp', ''' grey title="CMP gather" unit1=s label2=Offset unit2=km ''') # Near-offset section Result('noff','cmps', ''' window n2=1 n3=1024 | grey title="Near Offset Section" unit1=s label2=Distance unit2=km ''') # offset maps Flow('off1',None,'math n1=24 o1=0.4 d1=0.1 output=x1') Flow('off2',None,'math n1=24 o1=2.8 d1=0.05 output=x1') Flow('off','off1 off2','cat axis=1 \${SOURCES[1]}') Flow('offs','off','spray n=111') # Velocity analysis ################### vscan = ''' vscan offset=\${SOURCES[1]} v0=1.4 nv=61 dv=0.005 half=n semblance=y ''' pick = 'pick rect1=20 rect2=3 vel0=1.5' # compute semblance for one CMP Flow('vscan','cmp off',vscan) Plot('vscan', ''' grey color=j allpos=y title="Velocity Scan" unit1=s label2=Velocity unit2=km/s pclip=100 ''') # pick maximum semblance for one CMP Flow('pick','vscan',pick) Plot('pick0','pick', ''' graph transp=y yreverse=y min2=1.4 max2=1.7 plotcol=7 plotfat=10 pad=n wanttitle=n wantaxis=n ''') Plot('pick1','pick', ''' graph transp=y yreverse=y min2=1.4 max2=1.7 plotcol=0 plotfat=1 pad=n wanttitle=n wantaxis=n ''') Plot('vscan2','vscan pick0 pick1','Overlay') # compute semblance for every 10th CMP Flow('vscans','cmps offs','window j3=10 | ' + vscan) # pick max semblance for every 10th CMP Flow('picks0','vscans',pick) # interpolate picks on the original grid Flow('picks','picks0', 'transp | remap1 n1=1105 d1=0.05 o1=0 | transp') Result('picks', ''' grey color=j scalebar=y barreverse=y allpos=y bias=1.4 title="Stacking Velocity" label1=Time unit1=s label2=Distance unit2=km ''') # Normal moveout and stack ########################## nmo = ''' nmo offset=\${SOURCES[1]} velocity=\${SOURCES[2]} half=n ''' Flow('nmo','cmp off pick',nmo) Plot('nmo', ''' grey title="Normal Moveout" unit1=s label2=Offset unit2=km ''') Result('nmo','cmp vscan2 nmo','SideBySideAniso') Flow('nmos','cmps off picks',nmo) Flow('stack','nmos','stack') Result('stack', ''' window n2=1024 | grey title=Stack unit1=s label2=Distance unit2=km ''') # Convert from time to depth ############################ # Dix inversion Flow('semb','vscans picks0','slice pick=\${SOURCES[1]}') Flow('vel0','picks0 semb', 'dix rect1=20 rect2=2 weight=\${SOURCES[1]}') # interpolate on the original grid Flow('vel','vel0', 'transp | remap1 n1=1105 d1=0.05 o1=0 | transp') Result('vel', ''' grey color=j scalebar=y barreverse=y allpos=y bias=1.4 title="Interval Velocity" label1=Time unit1=s label2=Distance unit2=km ''') Flow('image','stack vel', ''' time2depth velocity=\${SOURCES[1]} intime=y dz=0.005 nz=1001 ''') Result('image', ''' window n2=1024 min1=3 | grey title=Image label1=Depth unit1=km label2=Distance unit2=km ''') End() ```

1. Change directory
```cd hw6/blake
```
2. Run
```scons view
```
to generate figures and display them on your screen.
3. Modify the SConstruct file. Check your results by running
```scons view
```
again.

 Homework 6

Next: Completing the assignment Up: Homework 6 Previous: Homework 6

2013-11-21