GPR data migration

Note:

• This R-package is still in development, and therefore some of the functions may change in a near future.
• If you have any questions, comments or suggestions, feel free to contact me (in english, french or german): emanuel.huber@alumni.ethz.ch.

Objectives of this tutorial

Learn how to migrate GPR data.

Note that his tutorial will not explain you the math/algorithms behind the different processing methods.

Preliminary

File organisation

I suggest to organise your files and directories as follows:

/2014_04_25_frenke   (project directory with date and location)
/processing      (here you will save the processed GPR files)
/rawGPR          (the raw GPR data, never modify them!)
RGPR_tutorial.R  (this is you R script for this tutorial)


Install/load RGPR and set the working directory

• Install and load the RGPR-package

# install "devtools" if not already done
if(!require("devtools")) install.packages("devtools")
devtools::install_github("emanuelhuber/RGPR")
library(RGPR)       # load RGPR in the current R session

• Set the working directory:

myDir <- "~/2014_04_25_frenke"  # adapt that to your directory structure
setwd(myDir)    # set the working directory
getwd()         # Return the current working directory (just to check)


A <- readGPR(fPath = "rawGPR/LINE00.DT1")   # the filepath is case sensitive!

## Warning in readGPR(fPath = "rawGPR/LINE00.DT1"): Use argument 'dsn' instead
## of 'fPath' because argument 'fPath' is deprecated.


Pre-processing

We assume that for each GPR record there is a file containing the (x, y, z) coordinates of every traces. The header of these files is “E”, “N”, “Z” instead of “x”, “y”, “z” because in topography “x” sometimes designates the North (“N”) and not the East (“E”) as we would expect. The designation “E”, “N”, “Z” is less prone to confusion and therefore we chose it!

1. Define the filepaths to the topo files:
TOPO <- file.path(getwd(), "coord/topo/LINE00.txt")

1. Read all the files with the funciton readTopo() that creates a list whose elements correspond to the GPR record and contain all the trace coordinates:
TOPOList <- readTopo(TOPO, sep = "\t")

## read /mnt/data/huber/Documents/WORKNEW/GPR_Project/RGPR-gh-pages/2014_04_25_frenke/coord/topo/LINE00.txt...

1. Set the list of coordinates as the new coordinates to the GPRsurvey object:
coord(A) <- TOPOList[[1]]

## [1] 23


DC shift removal

Remove the DC-offset estimated on the first n samples usind the function dcshift(). This function takes as argument the GPR object and the sample index used to estimate the DC shift (in this case, the first 110 samples):

A1 <- dcshift(A, 1:110)   # new object A1

## [1] 21


First wave break estimation and set time-zero

The first wave break time, t_{}, is estimated for each traces

tfb <- firstBreak(A1)   # take some time


Convert the first wave break time tf**b into time-zero t0 with firstBreakToTime0().

Here we define t0 = tf**b − a/c0, where a is the distance between the transmitter and receiver and c0 is the wave velocity in the media between the transmitter and receiver (in our case, air). The value a/c0 corresponds to the wave travel time from the transmitter to the receiver.

t0 <- firstBreakToTime0(tfb, A1)
time0(A1) <- t0     # set time0 to A1


To shift the traces to time-zero, use the function time0Cor.

A2 <- time0Cor(A1, method = "spline")

## [1] 21


Dewow

Remove the low-frequency components (the so-called “wow”) of the GPR record with:

A3 <- dewow(A2, type = "runmed", w = 50)

## [1] 21


Frequency filter

Eliminate the high-frequency (noise) component of the GPR record with a bandpass filter. We define as corner frequencies at 150 MHz and 260 MHz, and set plotSpec = TRUE to plot the spectrum with the signal, the filtered signal and the filter.

A4 <- fFilter(A3, f = c(150, 260), type = "low", plotSpec = TRUE)


## [1] 21


Time gain

Apply a power gain and a spherical gain to compensate for geometric wave spreading and attenuation (Kruse and Jol, 2003; Grimm et al., 2006).

A5 <- gain(A4, type = "power", alpha = 1, te = 150, tcst = 20)

## [1] 21

A6 <- gain(A5, type = "exp", alpha = 0.11, t0 = 0, te = 125)

## [1] 21


Topographic Kirchhoff migration

See Dujardin & Bano (2013, Topographic migration of GPR data: Examples from Chad and Mongolia, Comptes Rendus Géoscience, 345(2):73-80. Doi: 10.1016/j.crte.2013.01.003)

Pre-processing

Constant offset correction

Time correction for each trace to compensate the offset between transmitter and receiver antennae (it converts the trace time of the data acquired with a bistatic antenna system into trace time data virtually acquiered with a monostatic system)

A7 <- timeCorOffset(A6)

## [1] 21

plot(A7)


Time upsampling (sinc-interpolation) of the GPR data to reduce the aliasing risk.

A8 <- upsample(A7, n = c(3,1))

## [1] 21


Topographic Kirchhoff migration.

Vertical resolution of the migrated data: dz = 0.01m.

Dominant frequency from spec(A9): fdo = 80 MHz (used to estimate the Fresnel zone).

For the moment the algorithm works only with a constant radar wave velocity. In this example the velocity is:

vel(A8)         # velocity

## [1] 0.1

depthunit(A8)   # units: nano-second (ns)

## [1] "ns"


To change the velocity, simply do:

vel(A8)  <- 0.09        # velocity in ns

A9 <- migration(A8, type="kirchhoff", max_depth = 20,
dz = 0.01, fdo = 80)

## [1] 21

plot(A9)


You don’t see so much: we need some post-processing!

Post-processing

Trace smoothing with a Gaussian filter

A10 <- filter1D(A9, type = "Gaussian", w = 2.5)

## [1] 21


Automatic gain control

A11 <- gain(A10, type = "agc", w = 0.55)

## [1] 21


inverse normal transformations

A12 <- traceScaling(A11, type = "invNormal")

## [1] 21


Comparison before/after migration

Before migration

plot(traceScaling(A8, type = "invNormal"))

## [1] 22


After migration

plot(A12)