GPRlab: Ground Penetrating Radar

Data Processing and Analysis

Software User Manual

Hongqiang Xiong， Zhiyu

2023-08-03

1
 Table of contents

1 Introduction .................................................................................................................4
2 Software Installation ...................................................................................................5
2.1 Installing GPRlab in MATLAB ........................................................................... 5
2.2 Standalone desktop software" .............................................................................. 6
2.3 Uninstall GPRlab ................................................................................................. 6
3 software interface ........................................................................................................7
4 Data module ................................................................................................................8
4.1 Reading in data .....................................................................................................8
4.2 Data parameter setting ..........................................................................................9
5 Image processing module—GPR data processing ....................................................12
5.1 Data flow table ................................................................................................... 12
5.2 Process flow operations ......................................................................................13
5.3 Image Setting and Data Saving .......................................................................... 14
6 Processing algorithm in menu—Function ................................................................ 16
6.1 R_DC ..................................................................................................................16
6.2 R_background .................................................................................................... 17
6.2.1 Remove background .................................................................................... 17
6.2.2 SVD ............................................................................................................. 17
6.2.3 Dynamic window background .....................................................................18
6.2.4 Remove file ..................................................................................................18
6.3 Gain .................................................................................................................... 19
6.3.1 Exponent gain .............................................................................................. 19
6.3.2 Automatic gain .............................................................................................20
6.3.3 Manual gain ................................................................................................. 21
6.4 Math operations ..................................................................................................21
6.4.1 Plus .............................................................................................................. 22
6.4.2 Multiplication .............................................................................................. 22
6.4.3 Absolute ....................................................................................................... 22
6.4.4 Gradient ....................................................................................................... 22
6.4.5 Square .......................................................................................................... 23
6.5 Filter ................................................................................................................... 23
6.5.1 Vertical filtering ........................................................................................... 24
6.5.2 Horizontal filtering ...................................................................................... 26

2
 6.5.3 Predictive deconvolution ............................................................................. 26
6.6 TwoD filter ......................................................................................................... 26
6.6.1 mean filter ....................................................................................................28
6.6.2 F-k filter ....................................................................................................... 28
6.6.3 Empirical mode decomposition filter .......................................................... 29
6.6.4 Variational modal decomposition Filter .......................................................29
6.7 Hilbert Analysis ..................................................................................................30
6.8 Migration ............................................................................................................31
6.8.1 SAR ..............................................................................................................31
6.8.2 F-K migration .............................................................................................. 32
6.9 Wavelet ...............................................................................................................32
6.9.1 TwoD-dwt .................................................................................................... 33
6.9.2 OneD-cwt .....................................................................................................33
7 Wave view—waveform analysis view ......................................................................34
7.1 waveform analysis ..............................................................................................34
7.2 Waveform operation ........................................................................................... 35
8 Time-frequency analysis view—time-frequency analysis ........................................37
9 Examples ...................................................................................................................40
9.1 Case1 .................................................................................................................. 40
9.2 Case2 .................................................................................................................. 43
9.3 Case3 .................................................................................................................. 47
10 Version update instructions .....................................................................................51

3
 1 Introduction

1 Introduction

The software is an open-source and free software for Ground Penetrating Radar
(GPR) data analysis and research. The software was developed using the MATLAB
programming language and compiled with MATLAB R2020b. The software is named
GPRlab because it is designed for GPR technology research, and it provides
researchers with a laboratory to process, distribute, and visualize GPR data freely.
The software can read commercial software formats, such as dzt, rd3, DT1, and
matrix arrangement GPR data. It can save processing workflows and processed data
and transfer data to MATLAB's workspace. The software can display radar and
waveform graphs, and it has complete data processing functions, such as DC removal,
background removal, signal gain, mathematical operations, 1D filtering, 2D filtering,
waveform analysis, and other algorithms. The software's image functions include
saving images, zooming in and out of images, deleting data points, observing data
points, and saving data points.
The software is free and cannot be used for commercial purposes. We also hope
that you will not use it for research on offensive weapons. If you use the software in
your research, please acknowledge that you used GPRlab to facilitate software sharing
and dissemination. The ownership of GPRlab is the responsibility of Hongqiang
Xiong, and if you encounter any difficulties in using the software, please contact the
software author via email at 1014007697@qq.com.

4
 2 Software Installation

2 Software Installation

We recommend installing this software as an APP in MATLAB, although

GPRlab can also be installed as a standalone desktop software.

2.1 Installing GPRlab in MATLAB

(1) First, make sure that your computer has MATLAB R2020b or later versions
installed.
(2) As shown in Figure 2.1, open MATLAB, enter the APP module, click Install
App, and select the GPRlab.mlappinstall we provide.
(3) After installation, GPRlab will appear as an app in your MATLAB (as shown
in Figure 2.2), click to use it.

Figure 2.1 Installing GPRlab in MATLAB

Figure 2.2 GPRlab installed in MATLAB

5
 2 Software Installation

2.2 Standalone desktop software"

if you already have MATLAB R2020b or later version installed on your

computer, you can directly open the exe/for_redistribution_files_only and click on
GPRlab.exe to run GPRlab or create a shortcut of GPRlab.exe on your desktop.
If you do not have MATLAB R2020b or later version installed on your computer,
you can install GPRlab using the installation program MyAppInstaller_mcr.exe
located in the exe/for_redistribution.
Note: We recommend that you use the installation method provided in
section 2.1. The standalone desktop software is not a recommended option
as it may encounter some unknown issues.

2.3 Uninstall GPRlab

If GPRlab is installed as an APP in MATLAB, you can right-click on it in the

MATLAB APP module and select "Uninstall". If GPRlab is installed as a standalone
desktop software, the uninstallation method is the same as any other desktop software.

6
 3 software interface

3 software interface

The software interface is shown in Figure 3.1. The top of the software is the
menu bar, which contains all the algorithms provided by GPRlab. Below the menu bar
are three module windows provided by this software:
(1) The Data module, which is used for GPR data reading and parameter settings.
(2) The Image Processing module, which is used for GPR data processing.
(3) The Wave module, which is used for GPR waveform analysis, which can be
used to analyze the spectrum of a single signal.
(4) Time-frequency analysis module, which is used for analyzing the
time-frequency spectrum and synchronous compression spectrum of a A-Scan.
After reading the data and setting the parameters in the Data module, you can
switch to the Image Processing and Wave modules to process and analyze the data.

Figure 3.1 GPRlab interface

7
 4 Data module

4 Data module

The Data module is used to load data and set parameters.

4.1 Reading in data

When GPRlab reads in single-channel DZT or dzt files (GSSI devices), the time
window in the file header will be loaded.

When GPRlab reads in single-channel rd3 files (MALA devices), the header file
rd needs to be placed in the same directory as the data, and the time window in the
header file will be loaded.

When GPRlab reads in single-channel DT1 files (EKKO devices), the header file
HD needs to be placed in the same directory as the data, and the time window in the
header file will be loaded.

For matrix-arranged data (such as a 512x1600 B-Scan shown in Figure 4.1), it

can be directly read in. Supported formats include: txt, dat or csv (for text files with
delimiters); xls, xlsb, xlsm, xlsx, xltm, xltx or ods (for spreadsheet files). The speed of
reading data matrices in different formats varies, with txt, dat, and csv being the
fastest, and electronic spreadsheet formats being generally slower. It is recommended
to use csv files because they can be opened with Excel and the data can be directly
viewed.

8
 4 Data module

Figure 4.1 Matrix arranged data

When loading data, click the "Load" button. Loading large amounts of data may
take some time. Once the "success!" dialog box pops up, the data loading is complete.
After successful data loading, the values of the data can be viewed in the Raw Data
panel on the right-hand side. As shown in Figure 4.2, we have loaded a file named
"test.dzt".

Figure 4.2 Read in data

4.2 Data parameter setting

Setting parameters for GPR data is important because subsequent data processing
algorithms may use the parameters set. After reading in the data, GPRlab must set the
parameters in the Data information panel. The Data information panel includes the
9
 4 Data module

following parameters:

Antenna frequency (required);

Relative permittivity (required);
Time window (required);
Sampling interval (if the data is evenly spaced, fill in the spacing; otherwise, set it
to 0);
Start point (default is 0);
End point (default is 0).

If it is free acquisition (or time acquisition mode), please fill in the Start and
End point positions, and set the spacing to 0. The software will automatically
calculate the spacing after confirmation. If it is evenly spaced acquisition (or distance
acquisition mode), set End to 0, and the software will automatically calculate the
position of the end point.

Time zero position (the position of the zero point on the time window, default is
0).

If the data is in commercial software format, the software will automatically import
the time window, and the user needs to set the other parameters. After filling in the
parameters, click the "Yes" button to confirm. As shown in Figure 4.3, we have filled
in the parameters for "test.dzt" (which was acquired in free acquisition mode).

10
 4 Data module

Figure 4.3 Parameter settings. a) after filling in the parameters and b) after confirming the parameters

When the parameter data is successfully set, a "Succeed!" dialog box will pop up,
and the software's table below will display the size information of the loaded data,
including: Rows, Columns, Mileage, and Depth. Figure 4.4 shows the results of
setting the data loading and parameters for "test.DZT" in the Data module.

Figure 4.4 Data module is set successfully

11
 5 Image processing module—GPR data processing

5 Image processing module—GPR data processing

The Image processing module (as shown in Figure 5.1) is used to process GPR
data. After reading the data and corresponding parameters in the Data module, the
right-hand side of the Image processing module will display an image of the original
data, with Time and Depth on the left and right y-axes, and Mileage on the x-axis,
showing the entire profile. The lower part of the image area displays the display
settings and data saving options. On the left side of the image are the data processing
flowchart and data processing adjustment buttons.

Figure 5.1 Image processing module

5.1 Data flow table

The data processing flowchart in the Image processing module saves the
processing steps and corresponding parameters. Each processing flow will occupy one
row in the chart. The flowchart will display the algorithms added by the user.

The first column of each row in the flowchart is "Activation," which will be
12
 5 Image processing module—GPR data processing

assigned a value of 1 when a processing flow is added, indicating that the flow is
activated. The second column of each row is "Processing" (the processing algorithm),
and the subsequent columns are the parameters corresponding to the processing
algorithm. Users can select a certain algorithm and use the "up," "down," and "delete"
buttons to change the order of the processing algorithm in the flowchart. After setting
up the flowchart, click the "Yes" button, and the data will be processed. After waiting
for a period of time, the processed data will be displayed in the image area on the
right-hand side. As shown in Figure 5.2, we have simply added three processing flows
to process "test.DZT".

Figure 5.2 Processing test.DZT

5.2 Process flow operations

The operation of the processing flow mainly includes moving the position of the
processing flow, and saving and loading the processing flow. After selecting a certain
flow, the algorithm execution steps can be adjusted by moving the flow up or down.

The "up" button moves the selected flow up;

The "down" button moves the selected flow down;
The "delete" button deletes the selected flow;
"Save processing" saves the processing flow in the current processing chart, with
13
 5 Image processing module—GPR data processing

formats such as txt, csv, xlsx, xls, and dat. We recommend saving it in csv format
for easy opening with office software. An example of a saved flow is shown in
Figure 5.3.
"Load processing" loads an existing processing flow.

The "Reset_rowdata" button resets the processed data to the original data.
After clicking this button, the data in the image area will be reset to the entire
profile of the original data. This button makes it very convenient for using
GPRlab.

The "Yes" button applies the algorithms in the processing flowchart to process
the original data.

Figure 5.3 Saved process

Note: If the user installs the app in Matlab, the processed data and original
data will be directly transferred to the Matlab workspace and named
Processing_data and Raw_data.

5.3 Image Setting and Data Saving

The image defaults to a grayscale image displaying the entire profile. The
buttons in the lower part of the image area can be used to set the color and profile
length of the image display, where:

"Trace Length" determines how many traces are displayed in one image and is
associated with "Distance."
"Distance (m)" determines the length of the profile displayed in one image and is
associated with "Trace Length."
"Yes" confirms the size of the displayed image.

14
 5 Image processing module—GPR data processing

"Color" sets the color mode. GPRlab provides five color modes: Gray, Parula,
Turbo, Jet, and HSV.
"Save Data" saves the processed data in txt, csv, xlsx, or dat format.

After clicking the "Yes" button, the slider below the image will be activated. In
addition, clicking the image area once allows you to use the left and right arrow keys
on the keyboard to scroll the image, with one click scrolling one image distance.
Figure 5.4 shows processed data using the HSV color mode, and Figure 5.5 shows an
image with 1000 traces displayed.

Figure 5.4 HSV display data

Figure 5.5 1000 channels show the results in each picture

Note: You can use the "Reset_rowdata" button to reset the image display to
15
 6 Processing algorithm in menu—Function

show the entire profile.

6 Processing algorithm in menu—Function

6.1 R_DC

In general, the amplitude of the GPR echoes should be symmetrically distributed

about the mean value. However, due to the instability of the receiver circuit, some DC
offsets may be displayed. The DC offset is a physical phenomenon related to the
structure and composition of the sampling gate. DC removal means removing the DC
component of the signal. When the waveform is shifted from the zero baseline, it is
adjusted to the zero baseline, as shown in Figure 6.1. The algorithm formula for DC
removal is:

1 N 1
u *x ,t  u x ,t 
N
u
i 0
x ,t i (0-1)

In the formula, u *x ,t represents the signal after DC removal, and u x ,t represents the

signal before DC removal.

Figure 6.1 a) before removing DC and b) after removing DC

16
 6 Processing algorithm in menu—Function

6.2 R_background

As shown in Figure 6.2, GPRlab provides four ways to remove the background:
Remove background, SVD, Dynamic window background, and Remove file.

Figure 6.2 R_background

6.2.1 Remove background

The "Remove background" method selects the average of a portion of the profile
as the background value for each trace. Set the starting and ending traces for the
background. This is the most commonly used method. If the average of the entire
profile is used as the background, set Start = 1 and End to the total number of traces in
the entire profile. Note: End must be greater than Start.

6.2.2 SVD

The SVD method uses singular value decomposition to remove the background.
This is a background removal method that uses principal component analysis (PCA).
SVD decomposes GPR data into N (sample number, for example 512 or 1024)

17
 6 Processing algorithm in menu—Function

components, with the first few components being the main information, decreasing in
magnitude. The components before Start are discarded, and the components after End
are discarded. Generally, Start should be greater than or equal to 2 but not exceed 5.
End is recommended to be set to 0 by default.
In simple terms, the main component is the strong background component.
Therefore, discarding the components is equivalent to removing the background of the
data. The specific mathematical process can be found by the user, and is not detailed
in this manual.

6.2.3 Dynamic window background

The Dynamic window background method removes the background by selecting

the average value of the surrounding traces as the background for each trace. This is
equivalent to using a sliding window to remove the background of the data. This
method is mainly used for high-frequency antennas or when there are many horizontal
target anomalies and other background removal methods are not suitable. For example,
when detecting steel bars, this method can be used. Note that this method will also
eliminate stratigraphic information.
Window length: The size of the sliding window.
Start point: The starting sample point to be processed. If this method is applied to
the entire trace, Start point = 1.
End point: The ending sample point to be processed. If this method is applied to
the entire trace, End point = N (number of sample points).

6.2.4 Remove file

The Remove file method removes the background by loading a background data
file. The background can be read in the same format as the data, and is a single trace
18
 6 Processing algorithm in menu—Function

data. This processing method is suitable for special cases such as antenna testing,
numerical simulation, and air-coupled antennas.
Note that when using this method, an external file is required. Therefore, if
the saved processing flow contains this method, the saved processing flow will be
invalid.

6.3 Gain

Gain is the process of amplifying data using a gain function curve. As shown in
Figure 6.3, GPRlab provides five gain methods.

Figure 6.3 Gain

6.3.1 Exponent gain

Exponent gain is the process of multiplying the original signal by the

corresponding exponential gain function. Therefore, the parameters are set to specify

this gain function: y  scale * x.exp onent

19
 6 Processing algorithm in menu—Function

scale: The multiplier before the exponential gain, cannot be 0.

exponent: The exponent, cannot be 0.

start: The multiplier for the starting sample point gain, cannot be 0.

end: The multiplier for the ending sample point gain, cannot be 0.

Note: The algorithm uses start and end to obtain the value of x.

6.3.2 Automatic gain

1）Liner Automatic gain

The Linear Automatic Gain method first calculates the average of the absolute
amplitude value within a certain window range, and then calculates the gain value
within the window based on the set "Average value", so that the average absolute
amplitude value within the window reaches the set value. The size of the window is
the parameter "length", and the smaller the window is selected, the smoother the gain
curve will be. The selected window overlaps with the adjacent window as each trace
data slides from the beginning to the end, and the overlap rate is 50%. The value of
"length" is determined according to the number of sample points, and in general, it
can be set to 1/7 of the total number of sample points.

2）RMS Automatic gain

The Root Mean Square Automatic Gain method calculates the root mean square
value based on the current sample point and the gain window (with the current sample
point as the center), and then divides the amplitude value of the current sample point
by this root mean square value. The size of the window is the parameter "length", and
the smaller the window, the stronger the constraint on high amplitude values, while
the larger the window, the stronger the constraint on low amplitude values. "Length"
is the number of sample points included in the window.

20
 6 Processing algorithm in menu—Function

3）GRMS Automatic gain

The Gaussian-Root Mean Square Automatic Gain method calculates the square
of the amplitude value based on the current sample point and the gain window (with
the current sample point as the center), and then uses a Gaussian window for
weighting to calculate the mean square value. Finally, the amplitude value of the
current sample point is divided by this mean square value. The size of the window is
the parameter "length", and the smaller the window, the stronger the gain, while the
larger the window, the weaker the gain. "Length" is the number of sample points
included in the window.

6.3.3 Manual gain

The gain function of the Manual Gain method is a piecewise linear function. The
breakpoint positions "points" and the corresponding gain multipliers "value" can be
freely set by the user. The setting method is to use commas or spaces to separate them.

Note that the first and last sample points must be set in the array.

6.4 Math operations

As shown in Figure 6.4, the Math Operations function provides basic

mathematical operations for each trace data, including addition, multiplication,
absolute value, gradient, and square root.

21
 6 Processing algorithm in menu—Function

Figure 6.4 Math operations

6.4.1 Plus

The Plus function adds a specified value "value" to each sample data. This
method can be used to adjust the size and contrast of the data.

6.4.2 Multiplication

The Multiplication function multiplies each sample data by a specified value

"times" (which cannot be set as 0). This algorithm belongs to a type of gain method
and adds a constant gain to the original data. This method can be used in conjunction
with other mathematical operations.

6.4.3 Absolute

The Absolute function calculates the absolute value of each sample data and then
multiplies it by "times" (which cannot be set as 0).

6.4.4 Gradient

In a trace data, the Subtract function subtracts the sample data at a certain
distance ("points") from each sample data and then multiplies it by a coefficient
"scale". The distance and coefficient can be freely set by the user and are denoted as r
22
 6 Processing algorithm in menu—Function

and a, respectively. The operation can be expressed as:

F ( x )  scale( f ( x )  f ( x  po int s ))

GPRlab provides gradient operations that can be performed in the upward or

downward direction.

6.4.5 Square

The Square function calculates the square root of each sample data, and then
multiplies it by a "times" parameter (which cannot be set as 0). This operation can
weaken strong reflections and to some extent enhance weak reflections. After the
operation, the amplitude differences of various reflection waves in the data are
reduced, which can achieve the purpose of data balance.

6.5 Filter

As shown in Figure 6.5, GPRlab provides one-dimensional filtering, including

vertical filtering, horizontal filtering, and deconvolution. The filtering methods mainly
used are frequency domain filtering, FIR wavelet filtering, and IIR Chebyshev Type II
filtering. The deconvolution method uses the predicted error deconvolution algorithm.

Figure 6.5 Filter

23
 6 Processing algorithm in menu—Function

6.5.1 Vertical filtering

Vertical filtering refers to filtering the data in the time domain, which is a
necessary step in GPR data processing. Filtering is the process of changing the
relative sizes of various frequency components of a signal or completely eliminating
certain frequency components. Commonly used filters are designed to pass certain
frequencies without distortion, while significantly attenuating or eliminating other
frequencies. The frequency range that passes without distortion is called the passband,
and the frequency range that needs to be attenuated or eliminated is called the
stopband. The lowest frequency in the passband is the low-pass frequency, and the
highest frequency in the passband is the high-pass frequency.

In actual field measurements, in order to preserve as much information as

possible, a full-band recording method is often used, so that the effective and
interfering waves are recorded simultaneously. In order to remove interference signals
from the data, filtering methods need to be used to eliminate interference waves based
on the different frequency ranges of effective signals and interference signals in the
data. If there is a clear boundary between the frequency spectra of the effective signal
and the interference signal, the passband of the filter can be designed as the main
frequency spectrum of the effective signal, and the interference noise band outside the
signal band can be filtered out to obtain the filtered result. In data acquisition, for each
trace data, a Fourier transform is first performed on it, and the spectral range between
the low-pass frequency and the high-pass frequency, that is, the passband spectrum, is
selected for inverse Fourier transform to obtain the filtered data.

Frequency domain filtering, also known as ideal frequency domain filtering,

transforms the data to the frequency domain by FFT, multiplies it by a rectangular
window, and then IFFT is used to obtain the filtered data. However, it should be noted
that the rectangular window corresponds to a sinc function in the time domain, which
extends infinitely. After adding the window in the frequency domain, the IDFT
24
 6 Processing algorithm in menu—Function

actually only takes a part of the main value. Therefore, when filtering, oscillations
(high-frequency signals cannot directly become 0) will be formed at the beginning
and end of the signal, causing distortion at the beginning and end positions of the
signal. This has little effect on the vertical direction, but for the horizontal direction in
the wave number domain, it will cause discontinuity and abrupt changes between the
front and back images.

FIR, also known as finite impulse response filtering, is a filter type with a ripple
equal to a wavelet, with a stopband attenuation of 40dB, a passband gain of 0dB, and
a passband ripple limit of 1dB.

IIR, also known as infinite impulse response filtering, is a filter type of

Chebyshev Type II, with a stopband attenuation of 40dB, a passband gain of 0dB, and
a passband ripple limit of 1dB.

In FIR and IIR filters, the stopband low frequency (stop1), passband low
frequency (pass1), passband high frequency (pass2), and stopband high frequency
(stop2) need to be set, as shown in Figure 6.6.

Figure 6.6 Schematic diagram of filtering parameters

25
 6 Processing algorithm in menu—Function

6.5.2 Horizontal filtering

The filtering method is consistent with the previous section, with the direction
being horizontal. The frequency in the wave number domain refers to the number of
wavelengths per meter, with a maximum frequency value of 0.5/sampling interval. It
should be noted that horizontal filtering often eliminates layer information and
highlights horizontal changes. Users need to understand the meaning of the frequency
in the horizontal direction in order to use this method proficiently.

6.5.3 Predictive deconvolution

The prediction error deconvolution is a vertical prediction filtering method

provided by the Wiener filtering method. An operator length, Length, and a prediction
step, Lag, need to be set. The operator length, Length, should be less than half of the
time window, and the Lag should be less than the operator length. When performing
multiple wavelet filtering, Lag should be greater than the position of the first
reflection in the time window. In general, Length can be set to three times the
antenna's main frequency cycle. If there is strong horizontal interference in the data
but some layer information needs to be preserved, this method can be used.

6.6 TwoD filter

Two-dimensional filtering is an operation performed on the entire

two-dimensional data. GPRlab provides three filtering methods: mean filtering,
median filtering, and f_k two-dimensional bandpass filtering, EMD and VMD
filtering as shown in Figures 6.7.

26
 6 Processing algorithm in menu—Function

Figure 6.7 TwoD filter median filter

This method views the GPR profile as a two-dimensional image and uses image
processing ideas for filtering. In median filtering, the sampling value of each point is
sorted among all sampling points within a certain size neighborhood window, and the
value in the middle is taken as the new sampling value for that point. Median filtering
is a nonlinear signal processing technique based on sorted statistical theory, which is
particularly effective for isolated noise points. The specific method is to use a "rows"
x "columns" two-dimensional sliding template. For each data point, the center of the
two-dimensional sliding template is placed on that point, and the sampling values
within the template are sorted by size to generate a monotonically increasing (or
decreasing) data sequence. The value at the middle position of the sequence is taken
as the value for that point, and the final filtered data is obtained. This algorithm has a
good effect on high-frequency random noise in the data, while preserving the details
of the image to some extent.

27
 6 Processing algorithm in menu—Function

6.6.1 mean filter

This method views the GPR profile as a two-dimensional image and uses image
processing ideas for filtering. Mean filtering takes the average of all the sampling
values within a certain size neighborhood window for each point's sampling value,
and takes the average as the new sampling value for that point. The specific method is
to use a "rows" x "columns" two-dimensional sliding template. For each data point,
the center of the two-dimensional sliding template is placed on that point, and the
sampling values within the template are averaged to obtain the average as the value
for that point, and the final filtered data is obtained. This algorithm has a certain
smoothing effect on the data, but it can also make the image blurry and lose details
and contours. The larger the parameter setting, the better the smoothing effect, but the
more blurred the image becomes.

6.6.2 F-k filter

F-k filtering is a typical two-dimensional filtering method that decomposes the

image into the two-dimensional frequency domain through two-dimensional Fourier
transform and performs bandpass filtering in the two-dimensional frequency domain.
Due to the characteristics of GPR echoes, clutter mainly manifests as changes in the
horizontal direction. Therefore, bandpass filtering in the k-direction is generally used
for filtering settings.

The filtering parameters are Lg (m) for the length of the low-frequency signal in
the horizontal direction and Ls (m) for the maximum displacement of the system's
micro-vibration. Lg represents the low-frequency part in the horizontal direction,
which is generally selected near 1/2 of the profile length. When Ls is unknown, a
small value, such as 10 times the trace spacing, can be selected.

The F-k filtering method is similar to the horizontal filtering method and will
suppress low-frequency information in the horizontal direction. This means that this
28
 6 Processing algorithm in menu—Function

algorithm will weaken multiple reflections and suppress low-frequency interference in

the horizontal direction of the layer.

6.6.3 Empirical mode decomposition filter

EMD is a novel time-frequency analysis method that decomposes a signal into

different hierarchical modal classifications called Intrinsic Mode Functions (IMFs).
IMFs differ from the traditional trigonometric-based decomposition approach as they
are not constructed using base functions. Instead, they represent the envelope
components containing different time-scale features of the original signal. EMD first
extracts IMF1 and then decomposes it layer by layer until IMF_n, leaving a residual
component. Therefore, shallow IMFs represent high-frequency variations, while
deeper IMFs represent low-frequency variations. In the context of Ground Penetrating
Radar (GPR), EMD can be applied in both horizontal and vertical directions for signal
decomposition. In practical applications, we have observed that EMD performs better
in horizontally filtering out multiple reflections compared to F-k filtering. This
method requires the user to set three parameters:
Firstly, the user needs to choose either the horizontal (Horizontal) or the vertical
(Vertical) direction.
Start represents the starting IMF component to be retained. The value of Start
cannot exceed 10.
-End represents the position of the ending IMF component to be retained,
counting from the last IMF.

6.6.4 Variational modal decomposition Filter

VMD is an advanced method similar to EMD but with some notable differences.
The IMFs obtained from VMD are signals with finite bandwidth, making the
decomposition more well-defined compared to EMD, and it allows users to specify
the number of IMFs. However, VMD's computational complexity is higher. When

29
 6 Processing algorithm in menu—Function

applying EMD and VMD methods, users need to understand the basic principles of
both approaches. VMD requires the setting of four parameters:
Firstly, the user needs to choose either the horizontal (Horizontal) or the vertical
(Vertical) direction.
Start represents the starting IMF component to be retained. The value of Start
cannot exceed 10.
End represents the ending IMF component to be retained.
Num represents the number of IMFs to be obtained through VMD
decomposition.
Note: VMD is typically more time-consuming, so it is recommended to first use
EMD to test the signal processing performance and then replace it with the VMD
method if necessary.

6.7 Hilbert Analysis

The A-scan of GPR is a real signal x(t ) that can be represented in polar

coordinates:

x(t )  A(t ) cos  (t ) (0-2)

Where A(t ) represents the instantaneous change of the signal's envelope

(instantaneous amplitude),  (t ) represents the phase change of the signal
(instantaneous phase), and  (t ) represents the frequency change of the signal
(instantaneous frequency). The Hilbert Analysis is a method for obtaining these three
instantaneous attributes. To obtain the instantaneous attributes, a Hilbert Analytic
function needs to be constructed:

X (t )  x (t )  ixˆ(t ) (0-3)

Among them, xˆ (t ) represents the Hilbert transform of the signal. Based on the

equation above, the instantaneous attributes can be calculated. GPRlab provides the

30
 6 Processing algorithm in menu—Function

calculation of Hilbert instantaneous amplitude, instantaneous phase, and instantaneous

frequency without the need for input parameters.

6.8 Migration

As shown in Figure 6 8, Migration currently provides Synthetic Aperture Radar

(SAR) and FK migration algorithms.

Figure 6.8 Migration

6.8.1 SAR

GRPlab provides synthetic aperture focusing algorithm (SAR) for offset

processing, which can be used to collect data from both ground-coupled antennas and
air-coupled antennas. This algorithm is similar to the time domain BP algorithm. In
this method, the medium is divided into two layers for synthetic aperture focusing.
The theoretical part of this method can be found in "Research on Synthetic Aperture
Focusing Imaging of Vehicle-mounted GPR Data" by Hongqiang Xiong (2018). This
method has good imaging effects and is not sensitive to parameters. This algorithm
requires four input parameters.

Point: the number of sampling points at the layer interface. If you want to
process ground-coupled antenna data with a homogeneous model, you can set it to the
maximum sampling point number. For air-coupled antenna data, it should be set to the
longitudinal sampling points corresponding to the target interface position. This

31
 6 Processing algorithm in menu—Function

parameter cannot be negative.

Alpha: antenna beam angle (in degrees). This value is the most important control
parameter for the aperture length and should be set based on the antenna performance.
With the performance requirements met, the smaller the value of this parameter, the
smaller the edge effect. If the antenna beam angle is unknown, you can try different
values to get the best image effect that can be used to infer the antenna beam angle.
This parameter has a value range of (0, 360).

Permittivity_1: relative permittivity of the first layer medium. If it is air, it should

be set to 1. This parameter has a value range of [1, 81].

Permittivity_2: relative permittivity of the second layer medium. If it is air, it

should be set to 1. This parameter has a value range of [1, 81].

6.8.2 F-K migration

The algorithm allows for setting a uniform permittivity for the space. It can be
used as a conventional method for offset processing.

6.9 Wavelet

According to Figure 6-9, this software includes 1D and 2D wavelet analysis tools,
specifically including 1D discrete wavelet denoising, 2D discrete wavelet denoising,
and 1D continuous wavelet filtering.

32
 6 Processing algorithm in menu—Function

Figure 6.9 Wavelet OneD-dwt

The software uses the Sym4 wavelet for 1D discrete wavelet denoising.

6.9.1 TwoD-dwt

The software uses the Sym4 wavelet for 2D discrete wavelet denoising.

6.9.2 OneD-cwt

The software uses the Morse wavelet for 1D continuous wavelet (CWT) analysis
and applies it to partial time windows and partial profiles within a given frequency
range. This method is the opposite of bandpass filtering and removes signal
components within the set frequency range, which is equivalent to removing in-band
interference.

LF (MHz): the low-frequency cutoff of the frequency band to be removed;

HF (MHz): the high-frequency cutoff of the frequency band to be removed;

LT (ns): the start time;

HT (ns): the end time;

Start track: the trace position where the profile starts;

End track: the trace position where the profile ends.

33
 7 Wave view—waveform analysis view

7 Wave view—waveform analysis view

The Waveform Analysis view allows users to perform time-domain waveform

analysis and spectrum analysis on the A-Scan of the ground-penetrating radar. This
provides users with a more intuitive understanding of the changes in the raw data
waveform and processed waveform. The waveform graph has functions such as direct
zooming in and out, saving data, and displaying data point values.

7.1 waveform analysis

As shown in Figure 7-1, the Waveform Analysis view displays the time-domain
signal waveform in the upper part and the half-spectrum (with a maximum frequency
of fs/2) in the lower part. The reason for displaying the half-spectrum is that the
discrete Fourier transform (DFT) is symmetric. There are four types of waveform
analysis supported:

Raw data: the original data, displayed as a green line;

Processed data: the processed data, displayed as a red line;

Background: the background data, which uses the mean background and is
displayed as a black line;

Difference: the difference between the original data and the background data,
displayed as a blue line.

Before displaying the waveform, the following parameters need to be set:

Track number: which track to observe specifically;

Background start: the start track of the background data;

34
 7 Wave view—waveform analysis view

Background end: the end track of the background data;

DC+YES: whether to keep the DC component in the spectrum graph;

DC+No: whether to remove the DC component in the spectrum graph.

Note: considering whether to keep the DC component is because the original

data and background data may have a large relative DC value, which makes the other
parts of the spectrum weak in comparison. Here, it is recommended to choose "No".
GPRlab can not only analyze the waveform data but also save the waveform and
spectrum data in the same format as other data

Figure 7.1 Wave module

7.2 Waveform operation

GPRlab provides various interactive operations for the waveform graph. When
the mouse is placed on the corresponding waveform or spectrum graph, the operation
options will appear in the upper right corner, which can be used to select saving,
brightening, reducing, and zooming. Brightening means that a portion of the points on
the image can be brightened, as shown in Figure 7-2:

35
 7 Wave view—waveform analysis view

Figure 7.2 Brush data

After brightening the image, right-clicking on the image will bring up the
Properties dialog box, and you can also perform operations such as copying, replacing,
and deleting data. Clicking on a small gray dot directly on the waveform graph with
the mouse will display the coordinates of that point in a small box, which can be
moved up, down, left, and right with the mouse, as shown in Figure 7-3.

Figure 7.3 Display data

36
 8 Time-frequency analysis view—time-frequency analysis

8 Time-frequency analysis view—time-frequency analysis

Time-frequency analysis view can be used for conducting time-frequency

analysis on GPR A-Scan data. The basic settings for this view are consistent with the
Wave view. In the view, we provide two methods for time-frequency analysis:
Generalized Phase-shifting S-Transform (GPST) as the conventional approach and
Synchrosqueezing Generalized Phase-shifting S-Transform (SSGPST) as its
high-resolution counterpart in the synchronous compressed domain.

GPST and SSGPST are newly proposed time-frequency analysis methods that
can be regarded as variations of ST (S Transform) and SSST (Synchrosqueezing S
Transform), or CWT (Continuous Wavelet Transform) and SS-CWT
(Synchrosqueezing Wavelet Transform), respectively. Users can make use of their
experience with publicly available ST, SSST, CWT and SS-CWT to adjust the settings
for the Time-frequency analysis view. The parameters that need to be set include:

Num_Fpoints: The number of frequency points to be retained (<sampling

points/2).

General parameter: This parameter allows users to adjust the behavior of the
time-frequency analysis method based on their specific requirements.

Figure 8.1 shows the Time–frequency analysis results of the 500th trace data.

37
 8 Time-frequency analysis view—time-frequency analysis

图 8.1 Time–frequency analysis

Figure 8.2 shows the time-frequency plot when the reserved frequency point is
200. The data is sampled 512 times.

图 8.2 Num_Fpoints = 200

As shown in the time-frequency plot Figure 8.3 when General parameter=0.6. It

can be observed that the temporal resolution increases while the frequency resolution
decreases.

38
8 Time-frequency analysis view—time-frequency analysis

图 8.3 General parameter = 0.6

39
 9 Examples

9 Examples

9.1 Case1

ase 1 data comes from testing a 30m arch bridge with four arches using the free
acquisition mode of a GSSI device.

(1) To read in the data, as shown in Figure 8-1, we can see that the data is
sampled at 512 and has 5607 traces.

Figure 9.1 Case1 loading data

(2) Setting the parameters, as shown in Figure 8-2a, since free acquisition is
used, the Sampling Interval is set to 0, and the End is set to 30. After clicking Yes,
GPRlab automatically calculates the Sampling Interval value as 0.0053, as shown in
Figure 8-2b. After that, as shown in Figure 8-3, we can see that the Mileage is set to
30m, and the Depth is calculated based on the dielectric constant as 3.182m.

40
 9 Examples

Figure 8.2 Case1 setting parameters. a) After entering the parameters and b) after confirming the parameters

Figure 8.3 Data module completes the setup

(3) Switching to the Image Processing module, we can observe the raw data as
shown in Figure 8-4.

41
 9 Examples

Figure 8.4 Raw data

(4) First, we applied the R_DC algorithm to prevent DC interference from the
device. Then, we applied the "Remove background" algorithm in R_background to set
the background data to 1-800 traces. This is because this data comes from the
instrument in a stationary state at the bridge head. Finally, through our observation,
we found that the strong multiple reflections between the air-coupled antenna and the
ground had interfered with the entire profile. We selected the f_k filter method in
TwoD filter to remove the multiple waves. Lg is set to 10 (1/3 of the profile), and Ls
is kept at the default value. The processing flow is set as shown in Figure 8-5.

Figure 8.5 Case1 processing flow settings

42
 9 Examples

(5) Click Yes to complete the data processing and adjust the color to obtain the
processing result as shown in Figure 8-6.

Figure 8.6 Processing results of Case1

The data and processing flow of this case are saved in examples/Case1. Users
can also directly load or view the processing flow we provided.

9.2 Case2

Case 2 data comes from a train-mounted GPR device in distance acquisition

mode, obtaining partial data in a high-speed rail tunnel.

(1) To read in the data, as shown in Figure 8-7, we can see that the data is
sampled at 512 and has 1600 traces.

43
 9 Examples

Figure 8.7 Case2 loading data

(2) Setting the parameters, as shown in Figure 8-8a, since distance acquisition is
used, the Sampling Interval is set to 0.02, and the End is set to 0. After clicking Yes,
GPRlab automatically calculates the End value as 31.98, as shown in Figure 8-8b.
After that, as shown in Figure 8-9, we can see that the Mileage is set to 31.98m, and
the Depth is calculated based on the dielectric constant as 3m.

Figure 8.8 Case2 setting parameters. a) After entering the parameters and b) after confirming the parameters

44
 9 Examples

Figure 8.9 The Data module completes the setup

(3) Switching to the Image Processing module, we can observe the raw data as
shown in Figure 8-10.

Figure 8.10 Case2 raw data

(4) First, we applied the "Remove background" algorithm in R_background to

set the background data to all 1-1600 traces. Then, we applied the "GRMS automatic
gain" method in Gain to increase the gain. The length is set to 20. Then we used the
f_k filter and media filter in TwoD filter to remove some multiple reflections and
equipment vibrations, and the latter is to smooth the profile. Finally, we applied an
FIR filter to remove some noise. The processing flow is set as shown in Figure 8-11.
45
 9 Examples

Figure 8.11 Case2 processing flow settings

(5) Click Yes to complete the data processing, and obtain the processing result as
shown in Figure 8-12.

Fig. 8.12 Processing results of Case2

46
 9 Examples

The data and processing flow of this case are saved in examples/Case2. Users
can also directly load or view the processing flow we provided.

9.3 Case3

Case 3 data comes from the 500MHz high-frequency channel GPR of Chang'e 4,
which collected data in front of the Von Karman Crater on the moon during the first
two lunar days. The data is about 105m, and the acquisition mode is unknown (here
we assume it is free acquisition).

(1) To read in the data, as shown in Figure 8-13, we can see that the data is
sampled at 2048 and has 9032 traces.

Figure 8.13 Case3 loading data

(2) Setting the parameters, as shown in Figure 8-14a, since free acquisition is
used, the Sampling Interval is set to 0, and the End is set to 105. After clicking Yes,
GPRlab automatically calculates the Sampling Interval value as 0.0116, as shown in
Figure 8-14b. After that, as shown in Figure 8-15, we can see that the Mileage is set to
105m, and the Depth is calculated based on the dielectric constant as 55.43m.

47
 9 Examples

Figure 8.14 Case3 setting parameters. a) After entering the parameters and b) after confirming the parameters

Figure 8.15 Data module completed setting

(3) Switching to the Image Processing module, we can observe the raw data as
shown in Figure 8-16.

48
 9 Examples

Figure 8.16 Case3 raw data

(4) Unlike Case 1 and Case 2, in Case 3, a provided processing flow (containing
7 algorithms) is loaded, and the parameters can be displayed in the Processing panel,
as shown in Figure 8-11.

Figure 8.17 Case3 loading processing flow

(5) Click Yes to complete the data processing, and then set the Color to Jet to
obtain the processing result as shown in Figure 8-18.

49
 9 Examples

Figure 8.18 Case3 processing results

50
 10 Version update instructions

10 Version update instructions

Version 9.0 updates:

1. Added EMD and VMD algorithm.

2. Added Time-frequency analysis view

3. Some changes to the default settings and prompts.

This update mainly provides some state-of-the-art time-frequency analysis

algorithms, modal decomposition methods, and synchronous compression techniques
in the field of signal analysis. It's important to note that these methods have not yet
been widely applied in practical work within the GPR.

51