Homework 4

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

$$\min \Vert\mathbf{N} \mathbf{f} - \mathbf{d}\Vert^2\;,$$ (4)

where $\mathbf{d}$ is the data, $\mathbf{f}$ is a matching filter, and $\mathbf{N}$ represents convolution of the noise model $\mathbf{n}_0$ with the filter. After minimization, $\mathbf{n} = \mathbf{N} \mathbf{f}$ becomes the new noise model, and $\mathbf{d}-\mathbf{n}$ 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 <rsf.h>

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);
}

Your task:

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