Per Halvorsen

pmhalvor@uio.no

Task 1 - Convolution and frequency spectra¶

An input signal $x(n)$ convolved by the FIR filter $h(n)$ is given by

$$ y(n) = x(n) * h(n) = \sum_{k=-\infty}^{\infty} h(k)x(n-k) $$

conv_in3190¶

Write a function conv_in3190 that convolves an input signal $x(n)$ with a FIR filter $h(n)$ and returns $y(n)$. For an input with length $x_{length}$ the "support" for the output signal $y(n)$ will be longer than the input signal. The length of this extension depends on the length of the filter, $h_{length}$.

In addition to the two input parameters h and x, there should also be an input parameter ylen, where one can specify whether the length of the output signal $y_{length}$ should be the same as the input signal without a time shift compared to the input signal, or if it should be $y_{length} = h_{length} + x_{length} − 1$ (e.g., by allowing ylen to take the values 0 or 1). The function definition with help text is shown below.

Answer:

In [1]:
def conv_in3190(x, h, ylen=0):
    N_x = len(x)
    N_h = len(h)

    if ylen == 0:
        N_y = N_x   # for n > N_x, y[n] = 0
    elif ylen == 1:
        N_y = N_x + N_h - 1 # non-truncated support

    y = [0] * N_y

    for n in range(N_y):
        for k in range(
            max(0, n-N_x+1),  # k starts from max(0, n-N_x+1) to avoid negative indices in x-array
            min(n+1, N_h)     # k goes up to min(n+1, N_h) to stay within the bounds of h-array
        ):
            y[n] +=  h[k] * x[n-k] # convolution sum, throws error if n-k < 0 or n-k >= N_x or k >= N_h
    return y

freqspec_in3190()¶

For a signal $x(n)$, the frequency spectrum (discrete-time Fourier transform) is defined as: $$ X(e^{j\omega}) = \sum_{k=-\infty}^{\infty} x(k)e^{-j\omega k} $$

Write a Matlab function freq_spec_in3190 that calculates the complex frequency spectrum $X(e^{j\omega})$ of a signal $x(n)$. In addition to the input parameter x, the function should also have the following two input parameters:

  1. N, which specifies the number of points on the unit circle at which the frequency spectrum should be sampled.
  2. fs, which specifies the sampling frequency of the signal, so that the function returns the frequency spectrum with physical frequency (in Hertz).

In addition to the frequency spectrum, your Matlab function should also return the frequency axis. The function definition with help text is shown below

Answer:

In [2]:
import math
import cmath

def freq_spec_in3190(x, N, sample_freq):
    X = [0] * N

    for k in range(N):
        for n in range(len(x)):
            # Calculate the DFT of x
            X[k] += x[n] * cmath.exp(-1j * 2 * math.pi * k * n / N)

    # Calculate the frequency axis
    freqs = [k * sample_freq / N for k in range(N)]

    return X, freqs

$x(n) = \sin(2\pi f_1 t) + \sin(2\pi f_2 t)$¶

We will now test these functions using the signal $x(n)$, the sum of two sinusoids, with the following parameters:

  • Length in time: $t_{length} = \text{1 second}$
  • Sampling frequency: $f_s = 100 \text{ Hz}$
  • Frequencies: $f_1 = 10 \text{ Hz}$ and $f_2 = 20 \text{ Hz}$

The frequencies should be filtering using FIR filter:

$$ h(n) = \frac{1}{5} \sum_{k=0}^4 \delta(n-k) $$

Plot the frequency spectrum $| H(e^{j\omega})|$ of $h(n)$. Then plot the frequency spectra $|X(e^{j\omega})|$ and $|Y(e^{j\omega})|$ of $x(n)$ and $y(n) = x(n) * h(n)$ in the same plot. Explain what you see. Include your code for generating $x(n)$, $y(n)$ and the frequency spectra.

Answer:

In [3]:
import matplotlib.pyplot as plt

t_length = 5
f_s = 100
f_1 = 10
f_2 = 20 

t = [i/f_s for i in range(t_length*f_s)]
len(t) 

x = [math.sin(2*math.pi*f_1*i) + math.sin(2*math.pi*f_2*i) for i in t]

# delta = lambda n: 1 if n == 0 else 0
# h_n = lambda n: (1/5) * sum(delta(n-k) for k in range(5))
h = [1/5] * 5
In [4]:
# plot frequency spectrum for h
H, f_H = freq_spec_in3190(h, 100, f_s)
plt.plot(f_H, [abs(H_i) for H_i in H])
plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of $h(n)$")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()
No description has been provided for this image
In [5]:
# plot frequency spectrum for x
X, f_X = freq_spec_in3190(x, 444, f_s) # N \in [25, 50, 100, 125, 250, 500] give perfect alignment, otherwise bumpy plots 
plt.plot(f_X, [abs(X_i) for X_i in X])
plt.title(r"$|X(e^{jw})|$ , Frequency spectrum of $x(n)$")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()
No description has been provided for this image
In [6]:
Y = conv_in3190(x, h)
Y, f_Y = freq_spec_in3190(Y, 444, f_s)
plt.plot(f_Y, [abs(Y_i) for Y_i in Y])
plt.title(r"$|Y(e^{jw})|$ , Frequency spectrum of $y(n)$")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()
No description has been provided for this image

In the frequency spectra of the double sinusoid, we saw 4 peaks in frequency magnitude, at 10, 20, 80, and 90. For certian sample sizes (25, 50, 100, 125, 250, 500) we saw perfect alignment on these peaks, meaning our frequency spectrum showed clear prominate peaks at these frequencies.

The frequency spectra for our filter shows zeros at 20, 40, 60 and 80 Hz, and peaks at 0 and 100 (though these peaks are likely due to distrotion at the edges of our filter).

The convolution of the two signals maintained the peaks at 10 and 90 Hz, ince these were not zeros in our filter. However, the prominate peaks at 20 and 80 Hz were removed, as these were zeros in our filter.

Task 2 – Noise removal¶

image.png

We have the following data:

  • seismogram1: Gather 1.
  • seismogram2: Gather 2.
  • offset1: Offset for gather 1 [m].
  • offset2: Offset for gather 2 [m].
  • t: Recording time for both gathers [s].

We will initially work only with seismogram1, which is a gather with a six-kilometer offset. Since we assume that our velocity model is the same regardless of location, and we have limited noise, a single gather is sufficient for our analysis in this task.

In an ideal and noise-free world, our collected seismic data $x(t)$ would look like the gather on the left in the above figure. Due to noise from the ship, waves, towed streamers, recording equipment, and many other factors, seismic data will always be affected by noise. In this task, unfortunately, high-frequency noise has crept into the data, as we can see in the gather on the right in the above image. We aim to reduce the noise in our data, and for this, we wish to apply an FIR filter.

In our seismic processing toolbox, we have the two lowpass FIR filters $h_1$ and $h_2$. These are also downloadable here: https://www-int.uio.no/studier/emner/matnat/ifi/IN3190/h24/teaching_material/mandatory/data/filters.m

In [7]:
h_1 = [0.0002, 0.0001, -0.0001, -0.0005, -0.0011, -0.0017, -0.0019,
    -0.0016, -0.0005, 0.0015, 0.0040, 0.0064, 0.0079, 0.0075, 0.0046,
    -0.0009, -0.0084, -0.0164, -0.0227, -0.0248, -0.0203, -0.0079,
    0.0127, 0.0400, 0.0712, 0.1021, 0.1284, 0.1461, 0.1523, 0.1461,
    0.1284, 0.1021, 0.0712, 0.0400, 0.0127, -0.0079, -0.0203, -0.0248,
    -0.0227, -0.0164, -0.0084, -0.0009, 0.0046, 0.0075, 0.0079, 0.0064,
    0.0040, 0.0015, -0.0005, -0.0016, -0.0019, -0.0017, -0.0011,
    -0.0005, -0.0001, 0.0001, 0.0002]

h_2 = [-0.0002, -0.0001, 0.0003, 0.0005, -0.0001, -0.0009, -0.0007,
    0.0007, 0.0018, 0.0005, -0.0021, -0.0027, 0.0004, 0.0042, 0.0031,
    -0.0028, -0.0067, -0.0023, 0.0069, 0.0091, -0.0010, -0.0127,
    -0.0100, 0.0077, 0.0198, 0.0075, -0.0193, -0.0272, 0.0014, 0.0386,
    0.0338, -0.0246, -0.0771, -0.0384, 0.1128, 0.2929, 0.3734, 0.2929,
    0.1128, -0.0384, -0.0771, -0.0246, 0.0338, 0.0386, 0.0014, -0.0272,
    -0.0193, 0.0075, 0.0198, 0.0077, -0.0100, -0.0127, -0.0010, 0.0091,
    0.0069, -0.0023, -0.0067, -0.0028, 0.0031, 0.0042, 0.0004, -0.0027,
    -0.0021, 0.0005, 0.0018, 0.0007, -0.0007, -0.0009, -0.0001, 0.0005,
    0.0003, -0.0001, -0.0002]

Task 2a¶

Plot the two filters $h_1(t)$ and $h_2(t)$ in the same figure. Then, use your function freq_spec_in3190 to calculate the frequency spectrum $|H_1 (e^{j\omega})|$ and $|H_2 (e^{j\omega})|$. Plot the two frequency spectra in the same figure in decibels, i.e., with values.

Answer:

In [8]:
# h_1 and h_2 
plt.plot(h_1, label=r"$h_1$", color='blue')
plt.plot(h_2, label=r"$h_2$", color='red')
plt.title(r"$h_1$ and $h_2$")
plt.legend()
plt.xlabel("t")
plt.ylabel("h(t)")
plt.show()
No description has been provided for this image
In [9]:
# define decibel conversion
dB = lambda H: 20 * math.log10(abs(H))
In [10]:
H_1, f_H_1 = freq_spec_in3190(h_1, 444, f_s)
H_2, f_H_2 = freq_spec_in3190(h_2, 444, f_s)

plt.figure(figsize=(15, 5)) 
plt.plot(f_H_1, [dB(H_i) for H_i in H_1], label=r"$|H_1(e^{jw})|$", color='blue')
plt.plot(f_H_2, [dB(H_i) for H_i in H_2], label=r"$|H_2(e^{jw})|$", color='red')
plt.title(r"$|H_1(e^{jw})|$ and $|H_2(e^{jw})|$")
plt.legend()
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image

$h_1(t)$ is a stricter filter, allowing only frequencies less than 10 Hz to pass or greater than 90 Hz to pass. $h_2(t)$ is a less strict filter, allowing frequencies less than 20 Hz to pass or greater than 80 Hz to pass.

Task 2b¶

For a seismic wavefield propagating through the ground, the absorption will depend on the frequency, with the high frequencies being absorbed the most. This means that the wavefield reflected from deeper layers will have "lost" the high frequencies. Therefore, when examining the frequency content of a seismic signal $x(t)$, it is common to look at the portion of the signal that contains the shallowest reflection(s).

Plot some near traces (those with low offset) over a time window that includes the shallowest reflection(s) and the corresponding frequency spectrum. Since the time window has a limited length, this will result in edge effects in the frequency spectrum. You should apply a window function to your signal before calculating the frequency spectrum to mitigate these effects; for example, see the function tukeywin in Matlab. The window probably only needs to modify a small number of samples at the beginning and end of the time window. Document which window function you use and its parameters, and create the following plots:

  • The near traces with and without a window function applied, in the same figure.
  • The frequency spectra of these, with and without a window function applied, in the same figure. Use a Decibel scale.

Explain what part of the frequency spectra that are due to the signal, and what part that is due to noise.

Answer:

In [11]:
import seisplot
from scipy.io import loadmat

filename = "37.mat"
data = loadmat(f"../data/active_seismics/{filename}")
print("Time:", data["t"].shape)
print("Offset:", data["offset1"].shape)
print("Seismogram:", data["seismogram1"].shape)

# calculate sampling frequency from time data
elapsed_t = (max(data["t"]) - min(data["t"]))[0]
sample_freq = math.floor(len(data["t"]) / elapsed_t)  # samples per second
print("Sample frequency:", sample_freq)

# calculate offset units from offset data
offset_units = (max(data["offset1"][0]) - min(data["offset1"][0])) / len(data["offset1"][0])
print("Offset step size:", offset_units)

# manually adjust time frame 
N = 50
t_start = 0
t_end = 1500


truncated_seismogram = data["seismogram1"][t_start:t_end, :N]
truncated_t = data["t"][t_start:t_end].flatten()
truncated_offset = data["offset1"][0][:N]

seisplot.plot(
    truncated_seismogram.T, 
    plottype="wiggle", 
    vaxis=truncated_t,
    haxis=truncated_offset, 
    vlabel=f"Time (seconds)", 
    hlabel="Offset (meters)", 
    width=10, 
    height=8, 
    title=f"Seismogram 1 ({filename}) - no window"
)
plt.show()

Time: (1501, 1)
Offset: (1, 601)
Seismogram: (1501, 601)
Sample frequency: 250
Offset step size: 9.983361064891847
No description has been provided for this image

Explore window options: https://docs.scipy.org/doc/scipy/reference/signal.windows.html

In [12]:
from scipy.signal import get_window, windows

# reccomended 
plt.plot(windows.tukey(t_end - t_start, 0.75), label="Tukey window (alpha=0.75)")
plt.plot(windows.tukey(t_end - t_start, 0.5), label="Tukey window (alpha=0.5)")  # good enough 
plt.plot(windows.tukey(t_end - t_start, 0.25), label="Tukey window (alpha=0.25)")

# others
plt.plot(windows.hann(t_end - t_start), label="Hann window")
plt.plot(windows.hamming(t_end - t_start), label="Hamming window")
plt.plot(windows.blackman(t_end - t_start), label="Blackman window")

plt.title("Scipy window options")
plt.legend()
Out[12]:
<matplotlib.legend.Legend at 0x115b339d0>
No description has been provided for this image
In [13]:
tukey_50 = windows.tukey(t_end-t_start, 0.50)
In [14]:
seismogram_tukey_50 = truncated_seismogram.T * tukey_50

seisplot.plot(
    seismogram_tukey_50, 
    plottype="wiggle", 
    vaxis=truncated_t,
    haxis=truncated_offset, 
    vlabel=f"Time (s)", 
    hlabel="Offset (m)", width=10, height=8, title=f"Seismogram 1 ({filename}) - " + r"Tukey ($\alpha=0.50$)"
)
Out[14]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Seismogram 1 (37.mat) - Tukey ($\\alpha=0.50$)'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image
In [15]:
tukey_75 = windows.tukey(t_end-t_start, 0.75)
seismogram_tukey_75 = truncated_seismogram.T * tukey_75
seisplot.plot(
    seismogram_tukey_75, 
    plottype="wiggle", 
    vaxis=truncated_t,
    haxis=truncated_offset, 
    vlabel=f"Time (s)", 
    hlabel="Offset (m)", 
    width=10, height=8, 
    title=f"Seismogram 1 ({filename}) - " + r"Tukey ($\alpha=0.75$)"
    )
Out[15]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Seismogram 1 (37.mat) - Tukey ($\\alpha=0.75$)'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image
In [16]:
hann = windows.hann(t_end-t_start)
seismogram_hann = truncated_seismogram.T * hann
seisplot.plot(
    seismogram_hann, 
    plottype="wiggle",
    vaxis=truncated_t,
    haxis=truncated_offset, 
    vlabel=f"Time (s)", 
    hlabel="Offset (m)", 
    width=10, height=8, 
    title=f"Seismogram 1 ({filename}) - Hann"
)
Out[16]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Seismogram 1 (37.mat) - Hann'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image

Notice how the edges flattened out, but the middle of the signal looks the same. This shows how windowing reduces the effects of edges.

Frequency spectra with no-windowing

In [18]:
H_seismogram_freq_no_window = [freq_spec_in3190(truncated_seismogram[i], 512, sample_freq) for i in range(N)]
f_seismogram_no_window = H_seismogram_freq_no_window[0][1]
H_seismogram_no_window = [H for H, _ in H_seismogram_freq_no_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_seismogram_no_window, [dB(H) for H in H_seismogram_no_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of seismogram 1 (no window)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_seismogram_no_window, [abs(H) for H in H_seismogram_no_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of seismogram 1 (no window)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (absolute)")
plt.show()
No description has been provided for this image
No description has been provided for this image

Frequency spectra w/ windowing

In [17]:
alpha = 0.25
seismogram_tukey = truncated_seismogram.T * windows.tukey(t_end-t_start, alpha)

H_seismogram_freq_window = [freq_spec_in3190(seismogram_tukey[i], 1024, sample_freq) for i in range(N)]
f_seismogram_window = [f for _, f in H_seismogram_freq_window]
H_seismogram_window = [H for H, _ in H_seismogram_freq_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_seismogram_window[i], [dB(H) for H in H_seismogram_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of seismogram 1 (Tukey windowed $\alpha$"+f"={alpha})")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image
In [20]:
seismogram_hann = truncated_seismogram.T * windows.hann(t_end-t_start)

H_seismogram_freq_window = [freq_spec_in3190(seismogram_hann[i], 512, sample_freq) for i in range(N)]
f_seismogram_window = [f for _, f in H_seismogram_freq_window]
H_seismogram_window = [H for H, _ in H_seismogram_freq_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_seismogram_window[i], [dB(H) for H in H_seismogram_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of seismogram 1 (Hann windowed)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image
In [21]:
# plot only absolute magnitude to see what it looks like 
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_seismogram_window[i], [abs(H) for H in H_seismogram_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of seismogram 1 (Hann windowed)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (absolute)")
plt.show()
No description has been provided for this image

Looks like the true signal lies between 0 - 40 Hz and 210 - 250 Hz, with the rest being noise. This is becuase the magnitudes for frequencies in this region are compact and have a clear peak. The other frequency ranges have large variations in magnitude, and no clear peaks.

At around 125 Hz, we see a mirror effect for our frequencies. Though, through some research, I found its common to only look at the first 125 Hz, and disregard the higher frequencies.

Task 2c¶

Apply each of the two filters given by h1 and h2 to the seismic data, i.e., the entire gather, using the convolution function you have already implemented. The filtered gathers should have the same time axis as the unfiltered gather. Generate plots for the filtered gathers $y_1(t) = h_1(t)*x(t)$ and $y_2(t) = h_2(t)*x(t)$.

Based on your conclusion from Tasks 2a and 2b, which of the two filters is the best fit for noise reduction in our gather? Explain! Also, plot individual traces for both $y_1(t)$ and $y_2(t)$ in the same figure over an appropriate time window, which you should use to substantiate your explanation. For the remaining tasks, you will use your filtered gather.

In [18]:
# apply to single trace to see whay it looks like TODO remove
y_1 = conv_in3190(seismogram_tukey[0], h_1)
y_2 = conv_in3190(seismogram_tukey[0], h_2)

# plot y directly to see what it looks like TODO remoe
plt.figure(figsize=(15, 5))
plt.plot(y_1, label=r"$y_1$", color='blue', alpha=0.6)
plt.plot(y_2, label=r"$y_2$", color='red', alpha=0.6)
plt.title(r"$y_1$ and $y_2$")
plt.legend()
plt.xlabel("t")
plt.ylabel("y(t)")
plt.show()
No description has been provided for this image
In [20]:
import numpy as np

N = 600
y_full_1 = np.array([conv_in3190(data["seismogram1"].T[i], h_1) for i in range(N)])
y_full_2 = np.array([conv_in3190(data["seismogram1"].T[i], h_2) for i in range(N)])

seisplot.plot(
    y_full_1, 
    plottype="wiggle", 
    vaxis=data["t"].flatten(),
    haxis=data["offset1"][0][:N], 
    vlabel=f"Time (s)", 
    hlabel="Offset (m)", 
    width=10, height=8, 
    title=f"Filtered traces - {r'$h_1$'} ({filename})"
)
seisplot.plot(
    y_full_2,
    plottype="wiggle",
    vaxis=data["t"].flatten(),
    haxis=data["offset1"][0][:N], 
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Filtered traces - {r'$h_2$'} ({filename})"
)
Out[20]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Filtered traces - $h_2$ (37.mat)'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image
No description has been provided for this image
In [22]:
y_full_1.shape[1]
Out[22]:
1501
In [25]:
y_full_1_window = y_full_1 * windows.hann(y_full_1.shape[1])

H_1_freq_window = [freq_spec_in3190(y_full_1_window[i], 512, sample_freq) for i in range(600)]
f_1_window = [f for _, f in H_1_freq_window]
H_1_window = [H for H, _ in H_1_freq_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_1_window[i], [dB(H) for H in H_1_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of $y_1$ (Hann windowed)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image
In [26]:
y_full_2_window = y_full_2 * windows.hann(y_full_2.shape[1])

H_2_freq_window = [freq_spec_in3190(y_full_2_window[i], 512, sample_freq) for i in range(600)]
f_2_window = [f for _, f in H_2_freq_window]
H_2_window = [H for H, _ in H_2_freq_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_2_window[i], [dB(H) for H in H_2_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H(e^{jw})|$ , Frequency spectrum of $y_2$ (Hann windowed)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image

In 2a, I mentioned that $h_1$ was a stricter filter, only allowing frequencies less than 10 Hz to pass or greater than 90 Hz to fully pass though. When this filter is applied to our signal, which we previously inferred has a true signal between 0 - 40 Hz, we get a clear peak of magnitude around 10-15Hz. I think the clear peak in magnitude tells us this is the optimal filter.

(In a sense I "cheated" here by looking at the results in the next tasks for both filters. The stricter filter seems to have preserved more of the "interesting" part of our signal, while getting rid of a lot fo the wiggle noise seen in the trace plots.)

Task 3 - Far-field signature¶

In a perfect world, our source (the air gun) would emit a sound signal that approximates a pulse (a spike). The spectrum of a spike has infinite bandwidth, meaning it has a constant value. Such a spectrum is technically impossible to achieve in practice, and the sound signal from the air gun is instead a (very) band-limited version of a pulse. An example of such a band-limited pulse can be seen below.

image.png

Task 3a¶

We will now try to estimate the frequency spectrum of the pulse for the air gun used in our example. To do this, we use the direct arrival of the wavefield through the water, i.e., without reflections. As we saw in Figure 2, we can illustrate the propagation of a wavefield using so-called sound rays. In Figure 8, the wave energy that propagates directly between the source and the receivers is shown. This is called the direct arrival and is illustrated in Figure 8. image.png

This direct arrival for a noise-free shot is shown in Figure 9. By isolating this part of the signal, we can recreate the pulse and find the dominant frequency from the spectrum of the pulse. Plot the selected trace over a time window that contains the direct arrival, as well as the frequency spectrum of this signal in decibels.

Answer:

In [95]:
t_start = 40
t_end = 220
N = 120

# only plot the traces for y_full_1 for the first 200 samples, and 100 traces
seisplot.plot(
    y_full_1[0:N, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"][t_start:t_end].flatten(),
    haxis=data["offset1"][0][:N],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Filtered traces - {r'$h_1$'} ({filename})"
)
seisplot.plot(
    y_full_2[0:N, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"][t_start:t_end].flatten(),
    haxis=data["offset1"][0][:N],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Filtered traces - {r'$h_2$'} ({filename})"
)
Out[95]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Filtered traces - $h_2$ (37.mat)'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image
No description has been provided for this image
In [97]:
direct_arrival = y_full_1[:N, t_start:t_end]

H_direct_arrival_freq = [freq_spec_in3190(direct_arrival[i], 512, sample_freq) for i in range(N)]
f_direct_arrival = [f for _, f in H_direct_arrival_freq]
H_direct_arrival = [H for H, _ in H_direct_arrival_freq]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_direct_arrival[i], [dB(H) for H in H_direct_arrival[i]], label=f"Trace {i+1}")

plt.title(r"$|H_1(e^{jw})|$ , Frequency spectrum of direct arrival")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image
In [98]:
y_2_direct_arrival = y_full_2[:N, t_start:t_end]

H_2_direct_arrival_freq = [freq_spec_in3190(y_2_direct_arrival[i], 512, sample_freq) for i in range(N)]
f_2_direct_arrival = [f for _, f in H_2_direct_arrival_freq]
H_2_direct_arrival = [H for H, _ in H_2_direct_arrival_freq]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_2_direct_arrival[i], [dB(H) for H in H_2_direct_arrival[i]], label=f"Trace {i+1}")

plt.title(r"$|H_2(e^{jw})|$ , Frequency spectrum of direct arrival")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image

Some reassurance that these plots are supposed to look this wild: https://duckduckgo.com/?q=frequency+spectra+of+seismic+data&t=brave&iar=images&iax=images&ia=images

(Also, Cut off at 125 is perfectly reasonable.)

Task 3b¶

To avoid edge effects, a window should be applied to your signal before calculating the frequency response, similar to what you did in Task 2b. Plot the frequency spectrum of the pulse in decibels with and without the window function in the same plot. What is the dominant frequency?

Answer:

In [99]:
# without window 
direct_arrival_truncated = y_full_1[:N, t_start:t_end]

H_direct_arrival_freq_truncated = [freq_spec_in3190(direct_arrival_truncated[i], 512, sample_freq) for i in range(N)]
f_direct_arrival_truncated = [f for _, f in H_direct_arrival_freq_truncated]
H_direct_arrival_truncated = [H for H, _ in H_direct_arrival_freq_truncated]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_direct_arrival_truncated[i], [dB(H) for H in H_direct_arrival_truncated[i]], label=f"Trace {i+1}")

plt.title(r"$|H_1(e^{jw})|$ , Frequency spectrum of direct arrival")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image
In [100]:
# with window 
direct_arrival_window = y_full_1[:N, t_start:t_end] * windows.hann(t_end - t_start)

H_direct_arrival_freq_window = [freq_spec_in3190(direct_arrival_window[i], 512, sample_freq) for i in range(N)]
f_direct_arrival_window = [f for _, f in H_direct_arrival_freq_window]
H_direct_arrival_window = [H for H, _ in H_direct_arrival_freq_window]

# plot
plt.figure(figsize=(15, 5))
for i in range(50):
    plt.plot(f_direct_arrival_window[i], [dB(H) for H in H_direct_arrival_window[i]], label=f"Trace {i+1}")

plt.title(r"$|H_1(e^{jw})|$ , Frequency spectrum of direct arrival (Hann windowed)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude (dB)")
plt.show()
No description has been provided for this image

Dominating frequency is around 10 Hz. Though, using a window doesn't seem to have much an effect on the signal. The noise seems a bit more contained when using a window.

Task 3c¶

When we image a medium using wave energy, the resolution of the image is determined by the frequency of the wavefield and the properties of the medium (primarily velocity in our case). The vertical resolution indicates how accurately (vertically) we can differentiate structures in the image.

We start with Figure 10, where the boundary between two layers in the ground has a step with height h, causing it to have two different depths. The figure also shows two sources (red circles) fired over the deep and shallow parts of the step. To distinguish the step in the figure, the difference in the arrival time must cause a delay at the receiver (green triangles) of approximately 1/8 of the wavelength.

image.png

The formula for frequency is given by $f = c/\lambda$, where c is the speed of sound and $\lambda$ is the wavelength. If we assume that the layer causing the step has a velocity of 3000 m/s, what is the vertical resolution if we assume that the frequency in the wavefield is the same as the dominant frequency in the previous task?

Answer: $$ f = \frac{c}{\lambda} \rightarrow \lambda = \frac{c}{f} = \frac{3000 \text{ m/s}}{10 \text{ Hz}} = 300 \text{m} $$

The vertical resolution is then $300/8= 37.5 m$.

Task 4 – Reflections and multiples¶

Let’s take a closer look at the different events in our gather. Seismic reflections are simply wave energy that has traveled down into the ground and been reflected back from layers at various depths. However, this does not explain everything we see in our gathers.

Primary reflections as wave energy that has been reflected only once. In Figure 11 other wave paths that occur when wave energy propagates through the ground are illustrated. image.png In addition to downward-traveling wave energy being reflected back up, we also see upward-traveling wave energy being reflected downward. These reflections are called multiples when we observe them in the data.

We see three examples of wave energy that has been reflected more than once between the source and the receiver:

  1. The wave energy has been reflected twice in the water layer.
  2. The wave energy has been reflected once in the shallowest sedimentary layer and once in the water layer.
  3. The wave energy has been reflected twice in the shallowest sedimentary layer.

Task 4a¶

Let’s first look at seabed multiples. This is wave energy that has been reflected between the seabed and the sea surface at least once between the source and the receiver. If the time it takes for the wave energy to propagate vertically to the seabed and back to the sea surface is $t_w$, and the reflection coefficient for the seabed is $c_w$, then the time series describing the primary reflection and the multiples will be given by:

$$ x(t) = [1, 0, . . . , 0, −c_w, 0, . . . , 0, c^2_w, 0, . . . , 0, −c^3_w, 0, . . . ]  % = \sum_{n=0}^{\infty} (-1)^n c^n_w \delta(t − nt_w) $$

The separation between the impulses in this time series is $t_w$. Identify $t_w$ in your gather by examining the near trace, and create a plot of the near trace with the primary reflection and the first three seabed multiples. Also, plot the impulses in your gather at the first offset, with the main reflection as a cross and the reflections as circles, similar to what was done in the plot on the right in Figure 4.

Answer:

Start by looking at the fitlered signal again, to get an idea of the time window.

In [24]:
plt.figure(figsize=(15, 5))
plt.plot(y_full_1[0])
plt.xlabel("sample")
plt.ylabel("amplitude")
plt.title(r"First trace of $y_1$")
plt.show()
No description has been provided for this image

The first and second big peak occur in isolation, at roughly samples 250 and 475. The thrid peak seems to be a combination of two impulses, with different amplitudes/frequencies, around sample numbers 650-700. The trailing peaks also seem to be combinations of different multiples, getting progressively weaker as the sample count increases.

Now let's replot the full gather, to see what ranges we should zoom in on to plot the start of each reflection and multiple.

In [28]:
seisplot.plot(
    y_full_1,
    plottype="wiggle",
    vaxis=data["t"].flatten(),
    haxis=data["offset1"][0][:y_full_1.shape[0]],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8, 
    title=f"Full gather - {r'$y_1$'} ({filename})"
)
Out[28]:
(<Figure size 1000x800 with 1 Axes>,
 <Axes: title={'center': 'Full gather - $y_1$ (37.mat)'}, xlabel='Offset (m)', ylabel='Time (s)'>)
No description has been provided for this image

We'll label reflection off the first layer with a red x, then its corresponding (and trailing) multiples with red o's.

In [35]:
t_start = 0
t_end = 3 * sample_freq
N=50

seisplot.plot(
    y_full_1[: N, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"].flatten()[t_start:t_end],
    haxis=data["offset1"][0][:N],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8, 
    title=f"1st layer reflection and multiples - {r'$y_1$'} - Zoomed in ({filename})"
)

plt.plot([100], [0.95], 'rx', markersize=15, label="reflection")
plt.plot([100], [1.82], 'ro', markersize=15, fillstyle='none', label="multiple #1")
plt.plot([100], [2.68], 'ro', markersize=15, fillstyle='none', label="multiple #2")
plt.legend(loc='upper right')
plt.show()
No description has been provided for this image

Notice our placement of the last multiple. When looking at the full gather, we can see that the first part of this section actually corresponds to another reflection, following a different trajectory than the reflections coming from the first layer. We choose roughly the middle of this sequence to show that we are interested in the bottom part of these combined signals.

Task 4b¶

Create a similar time series as in Task 4a for the reflection from the deepest sedimentary layer that has then been reflected in the water layer. Make a plot of the time series with the primary reflection and the first two seabed multiples. Plot the impulses in your gather at the first offset, similar to what was done in the plot on the right in Figure 4.

Answer: We'll label the first sediment layer similar as above, but in blue.

In [38]:
t_start = 0
t_end = 5 * sample_freq
N=50

seisplot.plot(
    y_full_1[:N, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"].flatten()[t_start:t_end],
    haxis=data["offset1"][0][:N],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8, 
    title=f"2nd layer reflection and multiples - {r'$y_1$'} - zoomed in ({filename})"
)

plt.plot([100], [2.55], 'bx', markersize=15, label="reflection")
plt.plot([100], [3.45], 'bo', markersize=15, fillstyle='none', label="multiple #1")
plt.plot([100], [4.33], 'bo', markersize=15, fillstyle='none', label="multiple #2")
plt.legend(loc='upper right')
plt.show()
No description has been provided for this image

Task 5 - Water velocity¶

In Figure 9, the direct arrival was indicated in our data. Since the water layer has an approximately constant sound velocity, we see that this takes the form of a straight line.

We know that velocity $v$ is defined as

$$ v = \frac{s}{t} $$

where $s$ is distance and $t$ is time. Create plots showing the absolute value of two traces at different offsets within a time window that contains the direct arrival. Use the formula for velocity to estimate the speed of sound in the water layer using the direct arrival.

Hint: We have previously mentioned that water typically has a sound velocity of around 1500 m/s. Consider how the choice of offset might affect the estimate of the sound velocity.

Answer:

In [111]:
N = 120
seisplot.plot(
    y_full_1[0:N, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"][t_start:t_end].flatten(),
    haxis=data["offset1"][0][:N],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Filtered traces - {r'$h_1$'} ({filename})"
)
plt.plot(
    [100,     1000],
    [0.275,    0.88],
    "g",
    label="direct arrival (impulse start)"
)
plt.plot(
    [100,    1000],
    [0.17,   0.78],
    "r",
    label="direct arrival (impulse end)"
)
Out[111]:
[<matplotlib.lines.Line2D at 0x144b62590>]
No description has been provided for this image

We drew lines for the start and end of the direct arrival, in green and red respectively. From the points we drew, and can caluate the slope of the line, which is the velocity of the sound wave in the water layer.

Green $$ v = \frac{\Delta \text{s}}{\Delta \text{t}} = \frac{(1000-100)m}{(0.88 - 0.275)s} = \frac{900m}{0.605s} \approx 1487.6 m/s $$

Red $$ v = \frac{\Delta \text{s}}{\Delta \text{t}} = \frac{(1000-100)m}{(0.78 - 0.17)s} = \frac{900}{0.61s} \approx 1475.4 m/s $$

From our approximated slopes, we can estimate the sound velocity in the water layer to be around 1480 m/s.

Task 6 – Sediment velocity I¶

We see from the gathers earlier in the task that reflections as a function of offset are hyperbolas. It can be shown that the vertical reflection time $t_0$ for a horizontal reflector in the medium is given by

$$ t_0^2 = t^2 - \frac{x^2}{v^2_{\text{NMO}}} $$

where $t$ is the reflection time, $x$ is the offset, and $v_{\text{NMO}}$ is the velocity that flattens the hyperbolic-shaped reflection. The length of the $v_{\text{NMO}}$ vector should be the same as the traces, i.e., one value for each time sample. If there are multiple layers, this formula will be an approximation.

We will try to estimate the correct sound velocities by testing values for $v_\text{NMO}$ and seeing which values flatten the reflections.

In [112]:
# doc-strings found at: github.com/pinga-lab/nmo-tutorial/
from scipy.interpolate import CubicSpline
import numpy as np

def nmo_correction(cmp, dt, offsets, velocities):
    """
    Performs NMO correction on the given CMP.
    
    The units must be consistent. E.g., if dt is seconds and
    offsets is meters, velocities must be m/s.
    
    Parameters
    ----------
    
    cmp : 2D array
        The CMP gather that we want to correct.
    dt : float
        The sampling interval.
    offsets : 1D array
        An array with the offset of each trace in the CMP.
    velocities : 1D array
        An array with the NMO velocity for each time. Should
        have the same number of elements as the CMP has samples.
        
    Returns
    -------
    nmo : 2D array
        The NMO corrected gather.
    """
        
    nmo = np.zeros_like(cmp)
    nsamples = cmp.shape[0]
    times = np.arange(0, nsamples*dt, dt)
    for i, t0 in enumerate(times):
        for j, x in enumerate(offsets):
            t = reflection_time(t0, x, velocities[i])
            amplitude = sample_trace(cmp[:, j], t, dt)
            if amplitude is not None:
                nmo[i, j] = amplitude
    return nmo

def reflection_time(t0, x, vnmo):
    """
    Calculate the travel-time of a reflected wave.
    
    Doesn't consider refractions or changes in velocity.
        
    The units must be consistent. E.g., if t0 is seconds and
    x is meters, vnmo must be m/s.
    
    Parameters
    ----------
    
    t0 : float
        The 0-offset (normal incidence) travel-time.
    x : float
        The offset of the receiver.
    vnmo : float
        The NMO velocity.
        
    Returns
    -------
    
    t : float
        The reflection travel-time.
        
    """
    t = np.sqrt(t0**2 + x**2/vnmo**2)
    return t

def sample_trace(trace, time, dt):
    before = int(np.floor(time/dt))
    N = trace.size
    samples = np.arange(before - 1, before + 3)
    if any(samples < 0) or any(samples >= N):
        amplitude = None
    else:
        times = dt*samples
        amps = trace[samples]
        interpolator = CubicSpline(times, amps)
        amplitude = interpolator(time)
    return amplitude

Task 6a¶

We already have an estimate for the water velocity from Task 5. Using the velocity formula, you can then estimate the sea depth.

Use a constant $v_\text{NMO}$ function to see how this value flattens the seabed reflection using the provided Python function nmo_correction.

Try varying the water velocity and observe how this affects the flatness of the seabed reflection. How sensitive would you say the NMO correction of the seabed is to varying water velocities? Use your best estimate of the water velocity for the remainder of the task.

Answer:

For implementation assistance, I found this notebook: https://github.com/pinga-lab/nmo-tutorial/blob/master/step-by-step-nmo.ipynb.

The first reflection looks to arrive at the first offset hydrophone at roughly 950 ms, or 0.95 s. Using the estimated water velocity of 1480 m/s, the reflected signal from our air gun traveled $0.95 * 1480 = 1410$ m. Given that our air gun and receiver are assumed to be at the surface, we can draw a right triangle to fine our depth, with a hypotenuse of half the traveled distance, or 705 m. The distance to this first hydrophone is 100m, the short of of our right triangle will be 50m (helft the distnace to build a right triangle with the depth as the long side). Using the Pythagorean theorem, we can find the depth of the sea floor to be $d = \sqrt{705^2 - 50^2} \approx 703.22 m$.

In [42]:
np.sqrt(705**2 - 50**2)
Out[42]:
np.float64(703.2247151515652)

After producing many plots w/ varying $v_\text{NMO}$, I would say this method is very sensitive and prone to error. See the plots below for each of the ranges of $v_\text{NMO}$ I tested. Or just skip past to task 6b. (Shortcut to Task 6b)

Between 1700 and 2000 m/s¶

In [403]:
v_nmo_txt = r"$v_\text{NMO}=$"

y_1_T = y_full_1.T
time = data["t"].flatten()
dt = time[1] - time[0]
offsets = data["offset1"][0][:y_1_T.shape[1]].astype(np.float64)
nsamples = y_1_T.shape[0]
noffsets = y_1_T.shape[1]

times = np.arange(y_1_T.shape[0])*dt
v1, t1 = 100, 0.275  # used from task 5
v2, t2 = 1000, 0.88

# equation for a straight line between two points.
v_nmo = v1 + ((v2 - v1)/(t2 - t1))*(times - t1)
print("v_nmo", v_nmo.shape, v_nmo[0], v_nmo[-1])

plt.imshow(y_1_T, aspect="auto", extent=[offsets[0], offsets[-1], time[-1], time[0]],cmap='seismic')
plt.title("Seismogram (non-NMO corrected)") 
plt.xlabel("Offset (m)")
plt.ylabel("Time (s)")
plt.show()

for v_nmo in np.arange(1750, 2500, 75):
    velocities = (v_nmo)*(times - t1)
    nmo = nmo_correction(y_1_T, dt, offsets, velocities)
    plt.imshow(nmo, aspect="auto", extent=[offsets[0], offsets[-1], time[-1], time[0]],cmap='seismic')
    plt.title(f"{v_nmo_txt}{v_nmo} ({filename})")
    plt.xlabel("Offset (m)")
    plt.ylabel("Time (s)")
    plt.show()

    seisplot.plot(
        nmo.T,
        plottype="wiggle",
        vaxis=time,
        haxis=offsets,
        vlabel=f"Time (s)",
        hlabel="Offset (m)",
        width=7,
        height=5, 
        title=f"{v_nmo_txt}{v_nmo} ({filename})"
    )
    plt.show()

v_nmo (1501,) -309.0909090909092 8616.528925619836
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

Task 6b¶

Now that we have estimated the sound velocity in the water layer, we will try to find the sound velocity in the shallowest sedimentary layer. We will do this through trial and error by inserting different values for $v_\text{NMO}$ below the seabed and finding the velocity that gives the flattest gather. Include a plot that shows the NMO-corrected gather with constant velocities in the water layer and the sedimentary layer below. The result will be better if we smooth the transition between the two layers below the seabed. In addition to a plot that shows the NMO-corrected gather with this smoothed velocity, also include the code that shows the smoothing you have chosen with comments.

Answer:

Trying to achieve the "flattest" gather through NMO-correction feels ambiguous considering the results found in 6a. There were quite a few extra artifacts (i.e. mirroring hyperbolas) that appear when applying an NMO correction. These make it hard to determine when the NMO-corrected gather is "flattest".

Instead, I think its easiest here to rather plot reflection times that closest follow the first reflection trajectory representing the seabed. This elimiates a computation step, by avoiding generating NMO corrected gathers for each test-velocity. Additionally, we can directly see the velocity that gives the best matching curve, instead of having to interpret the noisy NMO corrected gathers.

In solving this task in my proposed manner, I am still getting the same result as asked for, but just more efficiently.

In [432]:
seisplot.plot(
    y_full_1,
    plottype="wiggle",
    vaxis=data["t"].flatten(),
    haxis=offsets[:y_full_1.shape[0]],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"1st sedimentary layer velocity estimation ({filename})"
)
trace_numbers = list(np.arange(1, noffsets)*10)

for vnmo in [2200, 2400, 2600, 2800]:
    reflect_times = reflection_time(t0=2.6, x=offsets[:y_full_1.shape[0]], vnmo=vnmo)
    plt.plot(trace_numbers, reflect_times, linestyle=(0, (1, 10)), label=f"v_nmo = {vnmo}", alpha=0.75)

trace_steps = list(range(1, noffsets + 1, 25))
ax.set_xticks(trace_steps)
ax.set_ylim(nsamples*dt, 0)
ax.set_xlim(0.5, noffsets + 0.5)
fig.tight_layout()
plt.legend()
plt.show()

/opt/homebrew/Caskroom/miniconda/base/envs/in3190/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)
No description has been provided for this image

It seems the first sedimentary layer has a velocity of around 2400 m/s.

Task 7 – Sediment velocity II¶

So far, we have looked at sound waves that have been reflected back to the surface. However, this represents only part of the wave energy that returns to the surface from where there is an impedance contrast between two layers. Some of the wave energy will also bend (refract), as shown in Figure 13. If the wavefield hits the interface (layer boundary) between two layers at the critical angle (the threshold angle for total reflection), the wavefield will follow the layer boundary. Some of this wave energy will then return back to the surface. While reflections appear as hyperbolas in the data, we see refractions as linear events since they primarily depend on the velocity in the layer they are refracted through.

image.png

As we saw earlier, the reflection time forms a hyperbola. In our shot gathers, refractions appear in the far traces (those with large offsets) as straight lines. Therefore, we need a large distance between the source and the receivers to effectively utilize refractions. In areas with complex geology, such as the Gulf of Mexico, seismic survey configurations often employ multiple seismic vessels to achieve distances of tens of kilometers between sources and receivers. This is done to measure refractions that return to the surface from depths of several kilometers.

Task 7a¶

In your gather, try to identify the refraction that returns from the boundary layer between the water layer and the shallowest sedimentary layer. By using the velocity formula from Task 5, what sound velocity do you find for the shallowest sedimentary layer? Compare your answer with what you found in Task 6b.

Answer:

In [23]:
offsets = data["offset1"][0][:y_full_1.shape[1]].astype(np.float64)

seisplot.plot(
    y_full_1,
    plottype="wiggle",
    vaxis=data["t"].flatten(),
    haxis=offsets[:y_full_1.shape[0]],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"1st sedimentary layer velocity estimation ({filename})"
)
plt.plot(
    [1500, 4000],
    [1.4, 2.42],
    linestyle=(0, (1, 10)),
    label="refraction estimation",
    alpha=0.75
)

plt.legend(loc="upper right")
plt.show()
No description has been provided for this image
In [24]:
# find the slope of the estimated line representing the refraction
(4000 - 1500)/(2.42 - 1.4)
Out[24]:
2450.9803921568628

$$ v = \frac{s}{t} = \frac{4000 - 1500}{2.42 - 1.4} = \frac{2500}{1.02} \approx 2451 \text{ m/s} $$

This is very close to the 2400 m/s found in Task 6b.

Task 7b¶

In your gather, try to identify the refraction that returns from the boundary layer between the shallowest and the deepest sedimentary layer. By once again using the velocity formula from Task 5, what sound velocity do you find for the deepest sedimentary layer?

In [39]:
time = data["t"].flatten()
dt = time[1] - time[0]
In [40]:
t_start = int(2.8/dt)
t_end = int(3.8/dt)
N_start = 500
N_end = 600

print(t_start, t_end, N_start, N_end)

seisplot.plot(
    y_full_1[N_start:N_end, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"].flatten()[t_start:t_end],
    haxis=offsets[:y_full_1.shape[0]][N_start:N_end],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Zoomed plot ({filename})"
)

plt.plot(
    [5000, 6000],
    [3.375, 3.65],
    # linestyle=(0, (1, 10)),
    "r-",
    linewidth=3,
    label="refraction estimation",
    alpha=0.75

)

plt.show()

699 949 500 600
No description has been provided for this image

In this zoomed in plot, we can see parts of the first-layer refraction, found previously in 7a (top right). The prominate reflection across the middle to lower half of the plot is the first reflction of this first layer of sediment. There doesn't seem to be any clear refraction tangents in this signal, so we'll just have to use the slope of this first reflection at high offsets to estimate the velocity of the second sedimentary layer.

In [41]:
(6000 - 5000)/(3.65 - 3.375)
Out[41]:
3636.3636363636374

From our estimated plot, we can predict the velocity of the second sedimentary layer to be approximately 3636.3636 m/s.

Task 7c¶

The data we have used so far is from an old seismic survey. Recently, a new seismic survey has been conducted over the same area, but now with eight-kilometer-long streamers instead of six-kilometer-long streamers. The gather for the new dataset corresponding to what we have worked with so far can be found in the variable seismogram2. Repeat the exercise from Task 7b, but with the new dataset. Do you see any change in the estimated sound velocity? If so, how would you explain this difference?

Answer:

In [42]:
data["seismogram2"].shape, data["offset2"].shape, data["t"].shape
Out[42]:
((1501, 801), (1, 801), (1501, 1))
In [43]:
# filter the second seismogram
y_s2 = np.array([conv_in3190(data["seismogram2"].T[i], h_1) for i in range(801)])
y_s2.shape
Out[43]:
(801, 1501)
In [44]:
offsets_2 = data["offset2"][0][:y_s2.shape[1]]
offsets_2.shape
Out[44]:
(801,)
In [45]:
t_start = int(3/dt)
t_end = int(4/dt)
N_start = 500
N_end = 800

print(t_start, t_end, N_start, N_end)

seisplot.plot(
    y_s2[N_start:N_end, t_start:t_end],
    plottype="wiggle",
    vaxis=data["t"].flatten()[t_start:t_end],
    haxis=offsets_2[N_start:N_end],
    vlabel=f"Time (s)",
    hlabel="Offset (m)",
    width=10,
    height=8,
    title=f"Seismogram 2 - Zoomed plot ({filename})"
)

plt.plot(
    [5000, 6000, 7000],
    [3.375, 3.65, 3.94],
    # linestyle=(0, (1, 10)),
    "r-",
    linewidth=3,
    label="refraction estimation",
    alpha=0.75

)

plt.show()

750 1000 500 800
No description has been provided for this image
In [46]:
(7000 - 6000)/(3.94 - 3.65), (6000 - 5000)/(3.65 - 3.375), (7000 - 5000)/(3.94 - 3.375)
Out[46]:
(3448.275862068965, 3636.3636363636374, 3539.823008849558)
In [47]:
np.mean([(7000 - 6000)/(3.94 - 3.65), (6000 - 5000)/(3.65 - 3.375)])
Out[47]:
np.float64(3542.3197492163013)

Using the larger offsets, we get slightly smaller velocity estimates for the second layer. This is expected, since the further out we look, the close the slope become to the true velocity of the layer (less influence of the hyperbolic shape of the reflection). Though, the difference in velocity estimates is very small. We can use an average of the two segments draw to represent our estimated velocity for the second sedimentary layer.

Task 7d¶

Use your estimated velocities for the water layer and the two sedimentary layers, as well as the depth of the two sedimentary layers, to fill out a table with the depths and velocities of each sedimentary layer, as well as the velocity of the water layer.

Answer:

Description Value
Depth sedimentary layer 1 705 m
Depth sedimentary layer 2 1954.44 m
Velocity water layer 1480 m/s
Velocity sedimentary layer 1 2450 m/s
Velocity sedimentary layer 2 3542 m/s

The last piece of information we needed to calculate was the depth to the start of sedimentary layer 2.

We know that the first reflection traveled a total of 1410 m, taking 0.95 seconds to reach the first hydrophone.

On our plot in 4b, we found that the first reflection of the second sedimentary layer arrived at roughly 2.55 seconds.

So through subtraction, we find how long the signal was traveling through the first layer of sediment, $2.55 - 0.95 = 1.6$ seconds. Given the estimated velocity of the first sedimentary layer was approximately 2450 m/s, we can conclude the signal traveled a total of $1.6 * 2450 = 3920$ m through the first layer of sediment, 1960 m down and 1960 m back up again.

For simplicity, we let's assume the angle of reflection is the same in the second layer as in the first layer. We can then draw similar right triangles as we did in 6a to find the depth using a trigonomic ratio:

$$ \cos \theta = \frac{703}{705} = \frac{d}{1960} \quad \Rightarrow \quad d = 1954.44 m. $$

In [48]:
(703/705)*1960
Out[48]:
1954.4397163120568

In reality, the angle of reflection could likely be different in the second layer. Though I think my calculations provide a good enough approximation for this exercise.