Imaging the underworld_GPR edition: All in One View

Last updated on 2025-04-14 | Edit this page

Overview

Questions

How do you write a lesson using Markdown and sandpaper?

Objectives

Explain how to use markdown with The Carpentries Workbench
Demonstrate how to include pieces of code, figures, and nested challenge blocks

Introduction

This is a lesson created via The Carpentries Workbench. It is written in Pandoc-flavored Markdown for static files and R Markdown for dynamic files that can render code into output. Please refer to the Introduction to The Carpentries Workbench for full documentation.

What you need to know is that there are three sections required for a valid Carpentries lesson:

questions are displayed at the beginning of the episode to prime the learner for the content.
objectives are the learning objectives for an episode displayed with the questions.
keypoints are displayed at the end of the episode to reinforce the objectives.

Challenge 1: Can you do it?

What is the output of this command?

R

paste("This", "new", "lesson", "looks", "good")

Output

OUTPUT

[1] "This new lesson looks good"

Challenge 2: how do you nest solutions within challenge blocks?

Show me the solution

You can add a line with at least three colons and a solution tag.

Figures

You can use standard markdown for static figures with the following syntax:

![optional caption that appears below the figure](figure url){alt='alt text for accessibility purposes'}

You belong in The Carpentries!

Callout

Callout sections can highlight information.

They are sometimes used to emphasise particularly important points but are also used in some lessons to present “asides”: content that is not central to the narrative of the lesson, e.g. by providing the answer to a commonly-asked question.

Math

One of our episodes contains $\LaTeX$ equations when describing how to create dynamic reports with {knitr}, so we now use mathjax to describe this:

$\alpha = \dfrac{1}{(1 - \beta)^2}$ becomes: $\alpha = \dfrac{1}{(1 - \beta)^2}$

Cool, right?

Key Points

Use .md files for episodes when you want static content
Use .Rmd files for episodes when you need to generate output
Run sandpaper::check_lesson() to identify any issues with your lesson
Run sandpaper::build_lesson() to preview your lesson locally

Content from Intro and Theory Refresher

Last updated on 2025-04-17 | Edit this page

Overview

Questions

What are the core components of a GPR system?
How does GPR work to detect subsurface changes?
What does a GPR trace represent?
How can we process and visualize SEG-Y GPR data using Python?

Objectives

Understand the structure and functionality of a GPR system.
Explain the concept of GPR signal reflection and trace acquisition.
Interpret what a single GPR trace and a 2D radargram represent.
Demonstrate loading and visualizing GPR SEG-Y data using Python and obspy.

Key Points

GPR uses EM waves to detect subsurface features.
A trace is a time series of reflections from a point.
SEG-Y files can be visualized using ObsPy.

What is geophysics?

In the broadest sense, the science of geophysics is the application of physics to investigations of the Earth, Moon and planets. The subject is thus related to astronomy. To avoid confusion, the use of physics to study the interior of the Earth, from land surface to the inner core, is known as solid earth geophysics. Applied geophysics’ covers everything from experiments to determine the thickness of the crust (which is important in hydrocarbon exploration) to studies of shallow structures for engineering site investigations, exploring for groundwater and for minerals and other economic resources, to trying to locate narrow mine shafts or other forms of buried cavities, or the mapping of archaeological remains, or locating buried pipes and cables – but where in general the total depth of investigation is usually less than 100 m.

Ground penetrating radar (GPR) is now a well-accepted geophysical technique. The method uses radio waves to probe “the ground”. In its earliest inception, GPR was primarily applied to natural geologic materials. Now GPR is equally well applied to a host of other media such as wood, concrete, and asphalt. The most common form of GPR measurements deploys a transmitter and a receiver in a fixed geometry, which are moved over the surface to detect reflections from subsurface features. Ground penetrating radar systems are conceptually simple; the objective is to measure field amplitude versus time after excitation. The heart of a GPR system is the timing unit, which controls the generation and detection of signals. Most GPRs operate in the time domain. GPRs are normally characterized by their center frequency. The common range varies between 50MHz to 2.5 GHz.

Most GPR consist of three components: a console, an antenna, and an encoder. The first two are mandatory. The third is typical for 99 percent of GPR users that work in the utility locating, civil engineering, concrete scanning, and archaeology industries. The console is the brains of the system. This data logger communicates with both the encoder and the antenna to initiate a signal and record the responses. The antenna is where the GPR signal is produced. The encoder, also referred to as a survey wheel, controls for locational accuracy of the GPR data. As the wheel turns a fraction of a turn, the console will initialize the antenna. Most GPR antenna actually consist of a transmitter and a receiver, and once it is initialized by the console, the transmitter will produce the signal. The receiver will record responses for a defined amount of time and once that “time window” has reached its maximum time, the console will terminate the receiver from recording any more data.

The GPR antenna produces an electromagnetic pulse at the ground surface. This pulse, which is in the form of a wave, travels into the subsurface and will continue until one of two things happen: the signal strength totally decays, or the wave encounters a change in the physical properties of the subsurface material. This change can be a small target such as a buried pipe or a large geological layer such a bedrock. When the wave encounters a change in material properties, some of the wave’s energy will reflect back to the ground surface and some will transmit further into the ground if any energy is left. This produces a one-dimensional view just below the antenna. If multiple of these one-dimensional views are collected next to each other (i.e. the GPR is pushed along the ground surface in a line), then a twodimensional view of the subsurface will be generated by the GPR. This two-dimensional view is equivalent to looking at the wall of an excavation trench. It is a vertical profile of the subsurface directly beneath the path that the GPR was pushed.

Diagram of a GPR system showing console, antenna, and encoder.

In this episode, we’ll introduce how to read and visualize a GPR trace from SEG-Y data using Python and the obspy library.

To visualize Ground Penetrating Radar (GPR) data, we can use the obspy library, which supports the SEG-Y format.

Below is a simple example of how to load and display a single GPR trace:

PYTHON

from obspy.io.segy.segy import _read_segy
import matplotlib.pyplot as plt

# Load SEG-Y using ObsPy (handles variable-length traces)
stream = _read_segy("LINE01.sgy", headonly=False)

# Access trace data
trace0 = stream.traces[0].data
n_traces = len(stream.traces)

# Plot first trace
plt.figure(figsize=(6, 4))
plt.plot(trace0)
plt.title("First Trace from LINE01.sgy")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()

PYTHON

import matplotlib.pyplot as plt
from obspy.io.segy.segy import _read_segy
import numpy as np

# Load SEG-Y file
segy_file = "LINE01.sgy"
stream = _read_segy(segy_file, headonly=False)

# Number of traces and max trace length
n_traces = len(stream.traces)
max_len = max(len(tr.data) for tr in stream.traces)

# Fill 2D array with seismic data
data = np.zeros((max_len, n_traces))
for i, tr in enumerate(stream.traces):
    data[:len(tr.data), i] = tr.data

# Plot
plt.figure(figsize=(12, 6))
plt.imshow(data, cmap="gray", aspect="auto", origin="upper", interpolation="none")
plt.title("Seismic Section (LINE01.sgy)")
plt.xlabel("Trace Number")
plt.ylabel("Time Sample")
plt.colorbar(label="Amplitude")
plt.show()

Content from AGC

Last updated on 2025-04-16 | Edit this page

Overview

Questions

What isAGC

Objectives

Describe AGC

Key Points

GPR AGC

Ground Penetrating Radar (GPR) write about AGC:

PYTHON

import numpy as np
import matplotlib.pyplot as plt
from obspy.io.segy.segy import _read_segy

def apply_agc(data, window_len):
    """
    Apply automatic gain control (AGC) on a 2D seismic array.

    Parameters:
        data (ndarray): [samples x traces]
        window_len (int): Window length in samples

    Returns:
        ndarray: AGC-applied seismic section
    """
    agc_data = np.zeros_like(data)
    half_window = window_len // 2

    for i in range(data.shape[1]):  # loop over traces
        trace = data[:, i]
        agc_trace = np.zeros_like(trace)

        for j in range(len(trace)):
            start = max(j - half_window, 0)
            end = min(j + half_window, len(trace))
            window = trace[start:end]
            rms = np.sqrt(np.mean(window ** 2)) + 1e-12  # avoid divide by zero
            agc_trace[j] = trace[j] / rms

        agc_data[:, i] = agc_trace

    return agc_data

# --- Load SEG-Y ---
segy_file = "LINE01.sgy"
stream = _read_segy(segy_file, headonly=False)

# --- Build data array ---
n_traces = len(stream.traces)
max_len = max(len(tr.data) for tr in stream.traces)
data = np.zeros((max_len, n_traces))
for i, tr in enumerate(stream.traces):
    data[:len(tr.data), i] = tr.data

# --- Apply AGC ---
window_length = 50  # samples (adjust based on data)
agc_result = apply_agc(data, window_length)

# --- Plot ---
plt.figure(figsize=(12, 6))
plt.imshow(agc_result, cmap="gray", aspect="auto", origin="upper")
plt.title("Seismic Section with AGC (window = {} samples)".format(window_length))
plt.xlabel("Trace Number")
plt.ylabel("Time Sample")
plt.colorbar(label="Amplitude (AGC-normalized)")
plt.show()

Content from Background Removal

Last updated on 2025-04-16 | Edit this page

Overview

Questions

What is bakground removal

Objectives

Describe background removal

Key Points

GPR bkr

Ground Penetrating Radar (GPR) write about bkr:

PYTHON

# Background removal by subtracting mean trace
def remove_background(data):
    background = np.mean(data, axis=1, keepdims=True)
    return data - background

# Apply it on your GPR section
data_bg_removed = remove_background(agc_result)

# Plot
plt.figure(figsize=(12, 6))
plt.imshow(data_bg_removed, cmap="gray", aspect="auto", origin="upper")
plt.title("GPR Section with Background Removed")
plt.xlabel("Trace Number")
plt.ylabel("Time Sample")
plt.colorbar(label="Amplitude")
plt.show()

Content from Time-Zero Correction

Last updated on 2025-04-16 | Edit this page

Overview

Questions

What is bakground removal

Objectives

Describe background removal

Key Points

GPR bkr

Ground Penetrating Radar (GPR) write about bkr:

PYTHON

import numpy as np
import matplotlib.pyplot as plt
from obspy.io.segy.segy import _read_segy
from scipy.signal import hilbert
import seaborn as sns

def detect_first_break(trace, threshold_ratio=0.05, min_sample=5):
    """
    Detect the first break index using an amplitude threshold on the envelope,
    starting the search from a specified minimum sample.
    """
    envelope = np.abs(hilbert(trace))
    threshold = threshold_ratio * np.max(envelope)
    for i in range(min_sample, len(envelope)):
        if envelope[i] > threshold:
            return i
    return min_sample  # fallback

def time_zero_correction_auto(data, threshold_ratio=0.05, min_sample=5):
    """
    Automatically estimate the target sample (median of all first breaks),
    align all traces to it, and trim the section to re-zero time.
    """
    n_samples, n_traces = data.shape
    shifts = [detect_first_break(data[:, i], threshold_ratio, min_sample) for i in range(n_traces)]

    # Estimate target sample using median of detected first breaks
    target_sample = int(np.median(shifts))
    print(f"Estimated target sample for alignment: {target_sample}")

    corrected = np.zeros_like(data)
    for i in range(n_traces):
        shift = shifts[i] - target_sample
        trace = data[:, i]
        shifted_trace = np.zeros_like(trace)
        if shift >= 0:
            shifted_trace[:len(trace)-shift] = trace[shift:]
        else:
            shifted_trace[-shift:] = trace[:len(trace)+shift]
        corrected[:, i] = shifted_trace

    # Trim the top 'target_sample' rows so time-zero starts from there
    corrected_trimmed = corrected[target_sample:, :]
    return corrected_trimmed, shifts, target_sample

# --- Load GPR (SEGY) data ---
segy_file = "LINE01.sgy"
stream = _read_segy(segy_file, headonly=False)

# Build 2D data array [samples x traces]
n_traces = len(stream.traces)
max_len = max(len(tr.data) for tr in stream.traces)
data = np.zeros((max_len, n_traces))
for i, tr in enumerate(stream.traces):
    data[:len(tr.data), i] = tr.data

# --- Auto Time-Zero Correction ---
threshold_ratio = 0.05
min_sample = 5
corrected_data, shifts, target_sample = time_zero_correction_auto(data, threshold_ratio, min_sample)

# --- First-Break Shift Histogram ---
plt.figure(figsize=(10, 3))
sns.histplot(shifts, bins=20, kde=True)
plt.title("Distribution of First-Break Shifts (Samples)")
plt.xlabel("Sample")
plt.ylabel("Frequency")
plt.tight_layout()
plt.show()

# --- Overlay Original vs Corrected for One Trace ---
plt.figure(figsize=(10, 4))
idx = 1000
plt.plot(data[:, idx], label="Original", alpha=0.6)
aligned = np.zeros_like(data[:, idx])
shift = shifts[idx] - target_sample
if shift >= 0:
    aligned[:len(aligned) - shift] = data[:, idx][shift:]
else:
    aligned[-shift:] = data[:, idx][:len(aligned) + shift]
plt.plot(aligned[target_sample:], label="Corrected", linestyle="--")
plt.title(f"Trace {idx} - Original vs Auto-Aligned")
plt.xlabel("Sample")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- Plot First-Break Picks on Sample Traces ---
plt.figure(figsize=(10, 6))
for i in range(5):
    idx = i * (n_traces // 5)
    plt.plot(data[:, idx], label=f"Trace {idx}")
    plt.axvline(x=shifts[idx], color='r', linestyle='--')
plt.axhline(target_sample, color='k', linestyle=':', label='Auto Target Sample')
plt.title("First-Break Picks on Sample Traces")
plt.xlabel("Sample")
plt.ylabel("Amplitude")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# --- Full Section: Before vs After Correction ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharey=False)

axs[0].imshow(data, cmap="gray", aspect="auto", origin="upper")
axs[0].set_title("Before Time-Zero Correction")
axs[0].set_xlabel("Trace")
axs[0].set_ylabel("Sample")

axs[1].imshow(corrected_data, cmap="gray", aspect="auto", origin="upper")
axs[1].set_title(f"After Auto Time-Zero Correction (Aligned to Sample {target_sample})")
axs[1].set_xlabel("Trace")

plt.tight_layout()
plt.show()

Content from Bandpass Filtering

Last updated on 2025-04-16 | Edit this page

Overview

Questions

What isAGC

Objectives

Describe AGC

Key Points

GPR AGC

Ground Penetrating Radar (GPR) write about AGC:

PYTHON

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt
from obspy.io.segy.segy import _read_segy

def bandpass_filter(data, lowcut, highcut, fs, order=4):
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist
    if not 0 < low < 1 or not 0 < high < 1:
        raise ValueError(f"Cutoff frequencies must be within (0, {nyquist:.2f} Hz). Got low={lowcut}, high={highcut}")
    b, a = butter(order, [low, high], btype='band')
    
    filtered = np.zeros_like(data)
    for i in range(data.shape[1]):
        filtered[:, i] = filtfilt(b, a, data[:, i])
    return filtered

def plot_frequency_spectrum(trace, fs, title):
    """
    Plot frequency amplitude spectrum of a single trace.
    """
    n = len(trace)
    freqs = np.fft.rfftfreq(n, d=1/fs)
    spectrum = np.abs(np.fft.rfft(trace))
    
    plt.plot(freqs / 1e6, spectrum, label=title)  # MHz
    plt.xlabel("Frequency (MHz)")
    plt.ylabel("Amplitude")
    plt.grid(True)

# --- Load SEG-Y GPR data ---
segy_file = "LINE01.sgy"
stream = _read_segy(segy_file, headonly=False)

# Build 2D data array [samples x traces]
n_traces = len(stream.traces)
n_samples = max(len(tr.data) for tr in stream.traces)
data = np.zeros((n_samples, n_traces))
for i, tr in enumerate(stream.traces):
    data[:len(tr.data), i] = tr.data

# --- Sampling rate ---
dt_us = stream.binary_file_header.sample_interval_in_microseconds  # microseconds
if dt_us == 0:
    raise ValueError("Sample interval is missing or zero in SEG-Y header. Set manually if needed.")

fs = 1e6 / dt_us  # Hz
nyquist = fs / 2
print(f"Sampling frequency: {fs/1e6:.2f} MHz, Nyquist: {nyquist/1e6:.2f} MHz")

# --- Dynamically determine bandpass range ---
lowcut = 0.05 * nyquist   # 5% of Nyquist
highcut = 0.95 * nyquist  # 95% of Nyquist
print(f"Applying bandpass filter: {lowcut/1e6:.2f}–{highcut/1e6:.2f} MHz")

# --- Apply bandpass filter ---
filtered_data = bandpass_filter(data, lowcut, highcut, fs)

# --- Plot before and after ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharey=True)

axs[0].imshow(data, cmap="gray", aspect="auto", origin="upper")
axs[0].set_title("Original GPR Data")
axs[0].set_xlabel("Trace")
axs[0].set_ylabel("Sample")

axs[1].imshow(filtered_data, cmap="gray", aspect="auto", origin="upper")
axs[1].set_title(f"Bandpass Filtered GPR\n({lowcut/1e6:.1f}–{highcut/1e6:.1f} MHz)")
axs[1].set_xlabel("Trace")

plt.tight_layout()
plt.show()

# --- Plot Frequency Spectra ---
middle_trace_idx = n_traces // 2
plt.figure(figsize=(10, 5))
plot_frequency_spectrum(data[:, middle_trace_idx], fs, title="Before Filter")
plot_frequency_spectrum(filtered_data[:, middle_trace_idx], fs, title="After Filter")
plt.title("Frequency Spectrum of Middle Trace")
plt.legend()
plt.tight_layout()
plt.show()

Content from Deconvolution

Last updated on 2025-04-16 | Edit this page

Overview

Questions

What isAGC

Objectives

Describe AGC

Key Points

GPR AGC

Ground Penetrating Radar (GPR) write about AGC:

PYTHON


import numpy as np
import matplotlib.pyplot as plt
from obspy.io.segy.segy import _read_segy
from scipy.fft import rfft, irfft, rfftfreq

def spectral_decon(trace, fs, stab=0.01):
    """
    Apply frequency-domain (spectral) deconvolution to a trace.

    Parameters:
        trace (1D array): input trace
        fs (float): sampling frequency in Hz
        stab (float): stabilization constant

    Returns:
        deconvolved trace (1D array)
    """
    n = len(trace)
    f = rfftfreq(n, d=1/fs)
    spectrum = rfft(trace)
    amp = np.abs(spectrum)
    phase = spectrum / (amp + 1e-12)  # preserve phase

    # Invert the amplitude spectrum with stabilization
    inverse_amp = 1 / (amp + stab * np.max(amp))
    decon_spectrum = inverse_amp * spectrum

    # Reconstruct the time-domain signal
    decon_trace = irfft(decon_spectrum, n=n)
    return decon_trace

# --- Load SEG-Y GPR data ---
segy_file = "LINE01.sgy"
stream = _read_segy(segy_file, headonly=False)

# Build 2D data array [samples x traces]
n_traces = len(stream.traces)
n_samples = max(len(tr.data) for tr in stream.traces)
data = np.zeros((n_samples, n_traces))
for i, tr in enumerate(stream.traces):
    data[:len(tr.data), i] = tr.data

# --- Sampling frequency from header ---
try:
    dt_us = stream.binary_file_header.sample_interval_in_microseconds
    if dt_us == 0 or dt_us > 10:
        print(f"Header sample interval {dt_us} µs seems invalid. Overriding to 0.5 µs.")
        dt_us = 0.5
except:
    print("Sample interval not found. Using default of 0.5 µs.")
    dt_us = 0.5

fs = 1e6 / dt_us  # sampling frequency in Hz
print(f"Sampling rate: {fs/1e6:.2f} MHz")

# --- Apply spectral deconvolution ---
decon_data = np.zeros_like(data)
for i in range(n_traces):
    decon_data[:, i] = spectral_decon(data[:, i], fs=fs, stab=0.5)

# --- Plot comparison ---
fig, axs = plt.subplots(1, 2, figsize=(14, 6), sharey=True)

axs[0].imshow(data, cmap="gray", aspect="auto", origin="upper")
axs[0].set_title("Original GPR Data")
axs[0].set_xlabel("Trace")
axs[0].set_ylabel("Sample")

axs[1].imshow(decon_data, cmap="gray", aspect="auto", origin="upper")
axs[1].set_title("Spectral Deconvolved GPR Data")
axs[1].set_xlabel("Trace")

plt.tight_layout()
plt.show()

Content from Wrap-up and Discussion

Last updated on 2025-04-14 | Edit this page

Questions

questions are displayed at the beginning of the episode to prime the learner for the content.
objectives are the learning objectives for an episode displayed with the questions.
keypoints are displayed at the end of the episode to reinforce the objectives.

Objectives

scsc.sdl ## Keypoints

dcs