Homework 4

Next: Spatial interpolation contest Up: Homework 4 Previous: Theory

# Match Filtering for Attenuation of Surface Seismic Waves

data
Figure 1.
Seismic shot record from sand dunes in the Middle East. The data are contaminated by ground roll propagating in the sand.

Figure 1 shows a section out of a seismic shot record collected over sand dunes in the Middle East. The data are contaminated by ground roll propagating in the sand. A major data analysis task is to separate the signal (reflection waves) from the noise (surface waves).

spec0
Figure 2.
Data spectrum. Solid line - original data. Dashed line - initial noise model and signal model.

A quick look at the data spectrum (Figure 2) shows that the noise is mostly concentrated at low frequencies. We can use this fact to create a noise model by low-pass filtering.

noise0
Figure 3.
(a) Noise model created by low-pass filtering of the original data. (b) Result of subtraction of the noise model from the data.

Figure 3 shows the noise model from low-pass filtering and inner muting and the result of subtracting this model from the data. Our next task is to match the model to the true noise by solving the least-squares optimization problem

 (4)

where is the data, is a matching filter, and represents convolution of the noise model with the filter. After minimization, becomes the new noise model, and becomes the estimated signal. Match filtering is implemented in program match.c. Some parts of this program are left out for you to fill.

 /* Match filtering */ #include int main(int argc, char* argv[]) { bool adj; int n1, n2, i1, i2, i, j, nf; float *data, *noiz, *filt; sf_file inp, out, oth; sf_init(argc,argv); inp = sf_input("in"); out = sf_output("out"); oth = sf_input("other"); if (!sf_getbool("adj",&adj)) adj=false; /* adjoint flag */ if (adj) { /* input data, output filter */ if (!sf_histint(inp,"n1",&n1)) sf_error("No n1="); if (!sf_histint(inp,"n2",&n2)) sf_error("No n2="); if (!sf_getint("nf",&nf)) sf_error("Need nf="); /* filter size */ sf_putint(out,"n1",nf); sf_putint(out,"n2",1); } else { /* input filter, output data */ if (!sf_histint(inp,"n1",&nf)) sf_error("No n1="); if (!sf_histint(oth,"n1",&n1)) sf_error("No n1="); if (!sf_histint(oth,"n2",&n2)) sf_error("No n2="); sf_fileflush(out,oth); /* copy data dimensions */ } filt = sf_floatalloc(nf); data = sf_floatalloc(n1); noiz = sf_floatalloc(n1); if (adj) { for (i=0; i < nf; i++) /* !!! COMPLETE LINE !!! */ } else { sf_floatread(filt,nf,inp); } for (i2=0; i2 < n2; i2++) { sf_floatread(noiz,n1,oth); if (adj) { sf_floatread /* !!! COMPLETE LINE !!! */ } else { for (i1=0; i1 < n1; i1++) data[i1] = 0.; } for (i=0; i < nf; i++) { for (i1=0; i1 < n1; i1++) { j=i1-i+nf/2; /* symmetric filter */ /* zero value boundary conditions */ if (j < 0 || j >= n1) continue; if (adj) { filt[i] += /* !!! COMPLETE LINE !!! */ } else { data[i1] += noiz[j]*filt[i]; } } } if (!adj) sf_floatwrite(data,n1,out); } if (adj) sf_floatwrite /* !!! COMPLETE LINE !!! */ exit(0); } 

1. Change directory to hw4/match
2. Run
scons view

to reproduce the figures on your screen.
3. Modify the match.c file to fill in missing parts.
4. Test your modifications by running the dot product test.
scons dot.test

Repeating this several times, make sure that the numbers in the test match.
5. Modify the SConstruct file to display the results of match filtering and include them in your assignment. Try improving the results by finding better parameters.
6. EXTRA CREDIT The seismic record in the example above is a section from 3-D data record (shot gather). Try using match filtering to predict one 2-D section out of 3-D data from another section. You may need to modify the program to make the filter two-dimensional.

 from rsf.proj import * # Critical parameters ########################### cut = 12 # cutoff frequency nf = 11 # filter length ########################### # Download data Fetch('dune3D.H','mideast') # Plotting macro def grey(title): return ''' window n1=490 | grey clip=2.5 title="%s" label1=Time unit1=s label2=Offset unit2=m ''' % title # Select one 2-D slice Flow('data','dune3D.H', ''' dd form=native | window n3=1 f3=2 n1=500 f1=100 | scale dscale=100 ''') Result('data',grey('Data')) # Create noise model by low-pass filtering Flow('noise0','data', ''' bandpass fhi=%g | mutter half=n v0=1500 t0=0.8 hyper=y tp=0.12 | cut n1=90 ''' % cut) Plot('noise0',grey('Noise Model')) # Signal = Data - Noise Flow('signal0','data noise0','add scale=1,-1 ${SOURCES[1]}') Plot('signal0',grey('Signal Model')) Result('noise0','noise0 signal0','SideBySideIso') # Plot spectrum Plot('spec','data', ''' spectra all=y | graph title="Average Spectrum" max2=3 label2= ''') Plot('nspec0','noise0', ''' spectra all=y | graph wanttitle=n wantaxis=n max2=3 plotcol=5 dash=1 ''') Plot('sspec0','signal0', ''' spectra all=y | graph wanttitle=n wantaxis=n max2=3 plotcol=4 dash=2 ''') Result('spec0','spec nspec0 sspec0','Overlay') # Matching filter program match = Program('match.c')[0] # Dot product test Flow('filt0',None,'spike n1=%d' % nf) Flow('dot.test','%s data noise0 filt0' % match, ''' dottest ./${SOURCES[0]} nf=%d dat=${SOURCES[1]} other=${SOURCES[2]} mod=${SOURCES[3]} ''' % nf,stdin=0,stdout=-1) # Conjugate-gradient optimization Flow('filt','data %s noise0 filt0' % match, ''' conjgrad ./${SOURCES[1]} nf=%d niter=%d other=${SOURCES[2]} mod=${SOURCES[3]} ''' % (nf,2*nf)) # Extract new noise and signal Flow('noise','filt %s noise0' % match, './${SOURCES[1]} other=${SOURCES[2]}') Flow('signal','data noise','add scale=1,-1 \${SOURCES[1]}') End() 

 Homework 4

Next: Spatial interpolation contest Up: Homework 4 Previous: Theory

2014-10-21