Processing Ground Penetrating Radar Data in MATLAB

Ground Penetrating Radar (GPR) is a non-destructive geophysical method that uses radar pulses to image subsurface structures. The raw data, however, is noisy and requires significant processing before it tells you anything useful. MATLAB is my go-to tool for this workflow.

Understanding GPR Data

A GPR survey produces a 2D matrix (B-scan) where:

• Rows = time samples (representing depth via two-way travel time)
• Columns = traces (representing horizontal position along the survey line)

Each column is a single radar trace — the reflected signal recorded at one position. The goal is to enhance reflections from subsurface features while suppressing noise.

Loading and Visualizing Raw Data

% Load GPR data (assuming .mat format)
load('gpr_survey.mat');  % variable 'data' is the B-scan matrix

% Display raw radargram
figure;
imagesc(data);
colormap(gray);
colorbar;
xlabel('Trace Number');
ylabel('Sample Number');
title('Raw GPR B-Scan');

The raw image is typically dominated by the direct wave (surface coupling) and horizontal banding. Processing is essential.

Standard Processing Pipeline

1. DC Shift Removal (Dewow)

Remove low-frequency "wow" caused by signal saturation:

function data_dewow = dewow(data, window)
    % Moving average subtraction
    [nsamples, ntraces] = size(data);
    data_dewow = zeros(size(data));
    
    for i = 1:ntraces
        trace = data(:, i);
        avg = movmean(trace, window);
        data_dewow(:, i) = trace - avg;
    end
end

data = dewow(data, 50);  % window size depends on data

2. Time-Zero Correction

Align all traces to the same starting time (the air/ground interface):

function [data_aligned, t0_indices] = time_zero(data, threshold)
    [nsamples, ntraces] = size(data);
    t0_indices = zeros(1, ntraces);
    
    for i = 1:ntraces
        trace = abs(data(:, i));
        idx = find(trace > threshold * max(trace), 1, 'first');
        t0_indices(i) = idx;
    end
    
    % Shift all traces
    min_t0 = min(t0_indices);
    max_shift = max(t0_indices) - min_t0;
    data_aligned = zeros(nsamples - max_shift, ntraces);
    
    for i = 1:ntraces
        shift = t0_indices(i) - min_t0;
        data_aligned(:, i) = data(shift+1:shift+size(data_aligned,1), i);
    end
end

3. Background Removal

Remove horizontal banding (the average trace):

% Simple mean subtraction
mean_trace = mean(data, 2);
data_bgr = data - mean_trace;

This is effective but can also remove real horizontal reflectors. A sliding-window approach is more conservative:

function data_bgr = background_removal_sliding(data, window)
    [nsamples, ntraces] = size(data);
    data_bgr = zeros(size(data));
    
    for i = 1:ntraces
        start_idx = max(1, i - floor(window/2));
        end_idx = min(ntraces, i + floor(window/2));
        local_mean = mean(data(:, start_idx:end_idx), 2);
        data_bgr(:, i) = data(:, i) - local_mean;
    end
end

4. Gain Application

Compensate for signal attenuation with depth:

% Exponential gain
function data_gained = apply_gain(data, alpha)
    [nsamples, ~] = size(data);
    gain = exp(alpha * (1:nsamples)');
    data_gained = data .* gain;
end

% SEC (Spreading and Exponential Compensation) gain
function data_sec = sec_gain(data, alpha, t_pow)
    [nsamples, ~] = size(data);
    t = (1:nsamples)';
    gain = (t.^t_pow) .* exp(alpha * t);
    data_sec = data .* gain;
end

5. Bandpass Filtering

Remove frequencies outside the GPR antenna's useful range:

function data_filt = bandpass_filter(data, fs, flo, fhi)
    [nsamples, ntraces] = size(data);
    data_filt = zeros(size(data));
    
    [b, a] = butter(4, [flo fhi] / (fs/2), 'bandpass');
    
    for i = 1:ntraces
        data_filt(:, i) = filtfilt(b, a, data(:, i));
    end
end

% For a 400 MHz antenna, typical bandpass: 100-800 MHz
data = bandpass_filter(data, sampling_freq, 100e6, 800e6);

Depth Conversion

Converting from time to depth requires knowing the signal velocity:

% v = velocity in m/ns (typical: 0.1 for wet soil, 0.15 for dry soil)
velocity = 0.12;  % m/ns
dt = 1 / sampling_freq;  % time step in ns
depth = (0:nsamples-1) * dt * velocity / 2;  % divide by 2 for two-way travel

Final Visualization

figure;
imagesc(distance_axis, depth, processed_data);
colormap(bone);
caxis([-std2 std2] * 3);  % clip at 3 standard deviations
xlabel('Distance (m)');
ylabel('Depth (m)');
title('Processed GPR Profile');
set(gca, 'YDir', 'reverse');
colorbar;

Conclusion

GPR data processing is a methodical pipeline where each step builds on the previous one. MATLAB's matrix operations and signal processing toolbox make it an excellent choice for this work. The key is understanding what each processing step does physically — don't just apply filters blindly.