GPR data migration


Table of Contents

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.


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

Read GPR data

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.


Add topographic data (coordinates)

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/RGPR/CODE/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]]

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

First wave break estimation and set time-zero

The first wave break time, tfb, is estimated for each traces

tfb <- firstBreak(A1, w = 20, method = "coppens")   # take some time

Convert the first wave break time tfb into time-zero t0 with firstBreakToTime0().

Here we define $t_0 = t_{\text{fb}} -\frac{a}{c_0}$, 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 $\frac{a}{c_0}$ 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")


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

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

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)

frequency filter

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)
A6 <- gain(A5, type = "exp", alpha = 0.11, t0 = 0, te = 125)

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)


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)

plot data

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

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

Topographic Kirchhoff migration.

Vertical resolution of the migrated data, we set dz = 0.01m.

The dominant return frequency can be estimated by visual inspection of the spectrum of A8 (to display the frequency spectrum of A8 use the function spec(), e.g. spec(A8)). Here, the dominant return frequency is estimated to be $80 MHz$ and therefore we set fdo = 80 (the dominant frequency is 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 <- migrate(A8, type="kirchhoff", max_depth = 20,
                 dz = 0.01, fdo = 80)

plot migrated data

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


Trace smoothing with a Gaussian filter

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

Automatic gain control

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

inverse normal transformations

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

Comparison before/after migration

Before migration

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

before migration

After migration


after migration