Detection of Underground Anomalies Using Analysis of Ground Penetrating Radar Attribute

Need of specifying underground construction works for supporting further tasks as maintenance, repairing, or setting up new underground structures. For these needs, ground penetrating radar, one of the efficient geophysical methods, can bring high-resolution and quick underground image revealing existence of both natural and artificial anomalies. Its fixed receiver-transmitter antennas setting as constant offset is commonly used in urban areas. Conventionally, hyperbolae events are crucial indicator for scattering objects as kinds of pipes, water drainage system, and concrete building structures as well as sink holes. Calculation of their depths and sizes requires migration analysis with the environment velocity. Migrated sections with different velocity show different chaos degrees of transformation from a hyperbola diffraction curve to its focused area. We have researched diagrams of different Ground Penetrating Radar attributes as energy, entropy, and varimax dependent on two variables, velocity and window zone covering diffraction events from a set of synthetic data and real data, in specifying the environment velocity. We have developed a novel technique for evaluation of the ground velocity and object’s size by combination of the new varimax diagram and the Kirchhoff migration method. The technique can define contribution of diffracted ground penetrating radar waves for building the diagram after removing the reflection contribution. The synthetic datasets consist of different random background noise levels and expressions of different-sized circular and rectangular pipes. The real data is measured for detecting two underground gas pipes in Ba Ria – Vung Tau province, Vietnam.


Introduction
Electromagnetic waves propagate through the medium and bounce back to receiver antenna after hitting the boundary of two zones of different electric permittivity. The electromagnetic waves characteristics lead to ability of investigation depth depending on their frequencies and electrically conductivity distribution in shallow surface environment (Doolittle andCollins, 1995, Smith andJol, 1995). Specifically, "skin effect" shows that the higher frequency Ground Penetrating Radar (GPR) section can illustrate higher resolution images rather than lower frequency GPR ones. However, the lower frequency GPR ones can provide images of deep structures better than the higher frequency ones.
Like seismic, the conventionally processed GPR data and their attributes could illuminate underground hidden geophysical or artificial structures and provide tool for estimation of the environment velocity ( The wave propagation velocity can help to correctly map the underground anomalies or structures through a specific processing step, migration. Known electric permittivity parameters are valuable for calculating electromagnetic velocity but it is unavailable where no drill hole is provided. To compensate the limitation, different techniques by analyzing GPR travel-time and amplitude sections are developed. Common Mid-Point (CMP) gather or a prior known object's depth can provide tools to evaluate velocity (Yilmaz, 2001 Zhao et al., 2015). Migration technique can collapse the diffraction hyperbolae into highly energy focus points in which the chosen migrated velocity responds to the environment velocity. In the Kirchhoff migration, summation of seismic amplitudes along a diffraction hyperbola is the secondary source amplitude positioning at the peak of the diffraction hyperbola. The GPR data can take an advantage of the migration technique when sharing the same dynamic characteristics with seismic wave. The hyperbola curves are also function of two-way travel time, velocity, size of diffractor and its depth, and antennas distance in which the velocity is analyzed (Sham and Lai, 2016).
Relationship between migrated velocity and a function of migrated data points over a specified window can form a useful GPR attribute diagram for estimation of environment velocity. The diagram can show its extreme value corresponding to the environment velocity. Entropy and its inverse, varimax, are of great indicator for velocity analysis in both geophysical data, seismic and GPR data (De Vries and Berkhout, 1984, Wiggins, 1978, Prego et al., 2017, Fomel et al., 2007, Levy and Oldenburg, 1987, Clair and Holbrook, 2017 can show that maximum energy difference can occur in the peak of diffraction hyperbolae with environment velocity applying to the GPR data In this research, we develop a novel approach by combining the Kirchhoff migration and the diagram of GPR attributes to evaluate the environment velocity and underground anomalies' properties (i.e., their locations and their sizes). We will apply the approach in the GPR synthetic data and real data in Ba Ria-Vung Tau province, Vietnam.

Method
Analyzing expressions of GPR amplitudes in migrated sections with a velocity band can provide tools of velocity estimation. According to research works (Yilmaz, 2001, Nguyen et al., 2017 when GPR data is migrated with correct velocity, the diffraction hyperbola turns into highly energy point at its peak. If the migrated velocity is not correct, over-migration or under-migration effects can occur with higher velocity or lower velocity than the environment velocity, respectively ( Figure 1). In the wrong velocity case, the GPR migrated amplitudes are still stretched in curved shapes upward or downward, causing highly entropy data (De Vries and Berkhout, 1984, Wiggins, 1978, Prego et al., 2017, Fomel et al., 2007, Levy and Oldenburg, 1987, Clair and Holbrook, 2017. For these research works, the varimax and its inverse, entropy, within window zone containing the peak of diffraction hyperbola can also reach maximum or minimum, respectively. Three functions of GPR migrated amplitude, varimax, entropy, and energy, are our interest. Varimax and entropy can relate to chaos degree of the GPR amplitude. For the energy variable, we are interested in investigating its strength over the reflection or diffraction events. Relationship between the GPR attributes (i.e., varimax, entropy, and energy) over two variables as velocity and a window can be of great explanation of the environment velocity. Also, we can investigate how useful the energy attribute can contribute to the velocity estimation from the small diffractor. The raw data is migrated with different velocities. The migrated data can fall into one of three categories: under-migration, correct migration, and over-migration. If correct migration is, no "smiling" or "frown eyes" are recognized (Yilmaz, 2001). c) The varimax section is formed by representation of varimax depending two variables, window length and velocity.

Analysis of GPR attributes
Building GPR attribute diagram needs a set of migrated GPR sections from a velocity band in advanced. The extremes of the GPR attribute (i.e., entropy, energy, and varimax) can reflect to the suitable environment velocity (Fomel et al., 2007, Clair andHolbrook, 2017). Each GPR attribute is formed using two variables, velocity and window zone through its equation as followed: (i) For the chaos of the GPR migrated data, entropy definition is introduced (Sava et al., 2004): (iii) For the inverse of the entropy, varimax term is used for minimum the order of "spikiness" or detecting the smallest number of the largest spikes as the equation (Wiggins, 1978): (v) We tempted to add one more energy term for researching the chaos level of migrated data as the equation (dGB Earth Sciences, 2015): where, s i is an amplitude in a GPR migrated section with a defined velocity. i is the location of the data point within the window zone. The flowchart for one example of a GPR attribute versus two variables, velocity, and window length (i.e., varimax) is expressed in the Figure 1.
The rectangular window zone has width as time distance in nanosecond and length as distance in meter covering the peak of the diffracted hyperbola. The maximum varimax can show the chosen velocity in which its migrated data show the smalless number of spikiness (Wiggins, 1978).

The new technique: varimax diagrams made by from the diffraction contribution
We have developed a new workflow of calculating varimax section for defining environment velocity and size of an object. For small objects, the conventional workflow (Fomel et al., 2007, Clair andHolbrook, 2017) expresses that maximum varimax value in a full rectangular window can relate to environment velocity. Apparently, maximum varimax values can lead to smallest number of spikiness in the window area of the GPR migrated data. That is, the diffracted hyperbola can converge into the focus point. However, the focus data point cannot answer to the full size of the big object. For the big object, our technique can calculate contribution of diffracted GPR waves in varimax diagram in which the maximum varimax responds to the environment velocity.
We apply the idea of achieving smallest number of spikiness of the GPR data after mitigating diffraction effects. For the big object, there can be two separated zones, (i) zone being responsible for reflection effect and (ii) zone for diffraction effect. Then, after removal of reflection zone, the diffraction zone can converge into the stage of being the smallest number of spikiness of the GPR data with the suitable velocity responding to the max varimax parameter. Note that exact separation of reflection and diffraction is a big challenge (Fomel et al., 2007). Our workflow described in Figure 2 shows that two different varimax sections with two different window zones can produce two different velocities. Our technique needs three main processing factors: (i) calculation of a varimax diagram within a designated window zone from different migrated sections, (ii) the migrated section with velocity corresponding to the maximum varimax can give the object size, and (iii) condition for stopping the workflow depends on how dynamic characteristics of migrated events appear.
If correct migration occurs, the workflow stops. In case, over-migration case occurs, the new designated window zone does not include the reflection size of the object and a new varimax diagram is re-calculated. The full window zone shown in Figure  With the new varimax section, the migrated section ( Figure  2f) shows the more focused object image.

Numerical model 3.1.1. Building synthetic data
Setup model: The model includes six anomalies, three rectangular pipes and three circle pipes (Figure 3). In the model, the background velocity is 0.075 m/ns (equally, 0.75 x108 m/s) and its anomalies' ones is 0.122 m/ns (equally, 1.22 x108 m/s). The top of all the anomalies are in the same depth, 2 meter. Each anomaly type, circle or rectangular, has sizes as 0.1 m, 0.5 m, and 1 m.
Forward modelling: The modelling tool of the source code MATGPR (Tzanis, 2006, Tzanis, 2010) is used to build constant offset (CO) GPR data from the model (Figure 3). The tool utilized the modelling work of Bitri and Grandjean (1998) in which a phase shift technique in frequency-wavenumber domain and solution of 2D Maxwell 's equations is used for wavefield extrapolation.
We have added different noise levels into the synthetic GPR data for testing robustness and effectiveness of our workflow. The newdata combining modelled signal and noise contributions is formed by the equation: New_data=synthetic_data + constant_value * white_gausian_ noise

Processing and result
Our workflow using varimax section analysis is applied into three cases for the synthetic dataset (i) small objects including circle or rectangular pipes having diameter size as 0.1 m, (ii) big-sized circle (i.e., its diameter as 1 m), and (iii) bigsized rectangle (i.e., its diameter as 1 m). We would apply two routines for calculating environment velocity, (i) conventional and (ii) our suggest ones. For the other medium-sized objects as 0.5 m, we would compare their migrated results from the two routines with the discussed model in Figure 3.
Preparation for our varimax analysis, migrated sections with different velocities and survey settings in the Opend-Tect software (dGB Earth Sciences, 2015, Huck, 2012) are set up. All the GPR CO datasets are migrated with the velocity band ranging from 0.05 to 0.14 m/ns. For speeding migration process, parallel computing with four cores is applied in the calculation. The processor is Intel (R) Core™ i7-6700HQ CPU at 2.6 GHz in baseline running. Window length ranges from 0.05 m to 2.9 m horizontally and 6.35 ns vertically. The datasets consist of the pure synthetic data and different white Gaussian noise addition levels. Therefore, each 3D varimax is built up with three dimensions, velocity, window length, and level of noise (white Gaussian noise addition level) defined as constant_value. For 3D view, inline and crossline are nonlinearly defined as window length and level of noise, respectively. The velocity is linearly expressed as Z direction.
We applied analysis of attribute diagrams for three cases in the model ( Figure 3); (i) the smallest objects with the first object for the smallest circular pipe and the fourth for smallest rectangular pipe, (ii) the largest circle object for the third, and (iii) The largest rectangular object for the sixth.

Conventional routine
Case 1: the small objects The interest objects as small circles and rectangles with their diameters as 0.1 m, locate at the distances 5 m and 35 m, respectively. The three kinds of GPR attributes as energy, entropy and varimax are input for velocity estimation in an object. The window area for the attributes are full rectangles with the constant time gate (6.35 ns) and varied window length from 0.05 m to 2.9 m. Figure 5 provides the diagrams of entropy, energy, and varimax for the small objects from the synthetic GPR data with different noise levels. In energy analysis, the maximum energy values responding to velocity changes versus the window length (i.e., see the dashed white curve) although good estimation for environment velocity just works well with the  Figure 5).
Case 2: the largest circle object According to the representations of different GPR attributes for the largest circle object in Figure 6, the environment velocity is not correctly picked with the extreme values of the GPR attributes. It shows that the extreme values for the GPR energy section change due to the window length; meanwhile, the entropy or varimax cannot answer the environment velocity with their extreme values although they refer a constant velocity.
For analysis of entropy and varimax in the biggest circle object, the invariant velocity responding to the extreme values, minimum entropy or maximum varimax, is not equal to the environment velocity when the conventional workflow of full window zone covering the diffraction is applied. Figure 7 expresses three velocities corresponding to the extreme values shown by the three 3D entropy or varimax bars of the three circles objects. It illustrates that the chosen velocity relating to the extremes increases if the objects' sizes increase.
Case 3: The largest rectangular object In Figure 8, all the three GPR attributes (energy, entropy, and varimax) do not show any visible extreme values responding to the environment velocity. It shows that the full rectangular window zone is not helpful for illuminating the differences of diffractions or reflection signal. Figures 7 and 8 inspires us to try a new strategy in finding environment velocity from GPR diffracted signals of the biggest objects. For velocity estimation in traditionally window zone, some remarks are made; (i) energy sections cannot provide a good tool for velocity estimation, (ii) entropy or varimax for the circle type can show the constant velocity with different axis parameters but it is not the environment velocity if its object size is big, and (iii) entropy or varimax for the rectangular object type does not reflect any extreme values that connect to environment velocity.

Our new technique
Our new flowchart for defining environment velocity is provided (see Section 2.2). For analyzing velocity from the synthetic datasets, GPR data zone connecting to the biggest objects as circles and rectangles are discussed. Taking advantages of the migrated sections from the velocity band, several rounds of calculating 3D varimax distribution with contribution of diffraction are applied. The idea of the flowchart is that reflection contribution is removed out of the calculation of varimax and the extreme varimax will respond to the contribution of the diffraction.
Case 2: the largest circle object We have worked with four rounds of calculating varimax sections (Figures 9 and 10) for achieving the environment velocity and its sizes. The first round comes up with the velocity 0.082 m/ns showing the maximum varimax and its migrated section shows over-migration effect. Many over-migration effect or upward curves in the migrated section (Figure 9c) shows that the chosen velocity bigger than the environment one. In the second round, its window zone does not include the reflection zone as 0.15 m which is extracted from the migrated section of the first round (Figure 9b). The second-round migrated section with the chosen velocity, 0.080 m/ns, also has over-migration effect (Figure 9d) although the object size increase to 0.36 m better than the previous result 0.15 m.
In the third round, the migrated section also has over-migration effect with the better results, the object size 0.51 m and the chosen velocity 0.078 m/ns. The best-chosen velocity is 0.076 m/ns and the object size as 0.95 m in the migrated section for the fourth round shows the acceptable errors, 1.3% and 5%, respectively. The migrated section (Figure 10f) show the focused result image.
Case 3: The largest rectangular object Two rounds of calculating varimax sections is done for calculating the biggest rectangular size, 1 m and the environment velocity, 0.75 m/ns ( Figure 11). The first round ( Figure  11a) does not show any meaningful max varimax value that could help to define velocity or its migrated section. In the second round, reflection contribution of the rectangular pipe as 1.05 m is calculated from the synthetic GPR data (Figure 4) in advance. After removing the reflection contribution, the new varimax in the second round (Figure 11b) shows the great chosen velocity

Real data
The data was recorded in Ba Ria Vung Tau province, Vietnam for checking present map of two underground metal gas Some brief interpretation can be extracted from the processed data (Figure 12, right image). A significant layer lo-cates at the time 20 ns. Two strong symmetrical and unsymmetrical hyperbolae events related to diffraction effects at the two locations (x = 2 m, t = 30 ns) and (x = 9 m, t = 35 ns) can respond to the two underground metal gas pipes, NCS and BH, respectively. Our workflow of velocity estimation by the varimax diagrams (see Section 2.2) is applied to the first location where the symmetrical hyperbola peak is captured in distance-time domain, 2 meter and 30ns.
We have used our workflow to the real GPR data with different rounds (Figures 13 and 14).
• The first round: In the first round (Figure 13), the varimax is established with velocity band from 0.08 m/ns to 0.16 m/ns and the window length from 0.05 to 2.9 m. Obviously, the largest window zone can cover the biggest pipe. The varimax section ( Figure  13a) firstly shows the velocity corresponding to the max varimax (see white arrow) with different full window lengths. The zoom image extracted from the whole migrated section (Figure 13c) can illustrate several "smile curves" as over-migra-

Discussion
The synthetic datasets with different added error levels and real data are tested using varimax diagrams. The testing shows that using varimax sections can define the environment velocity.
For small objects, the traditional entropy or varimax diagram shows their great applications in velocity estimation. The reason is that diffraction mainly comes from the small objects.
In big circle object's cases, traditional rectangular window for a zone of hyperbola can lead to the over-migrated effect with its resulted velocity. It could be explained when Kirchhoff migration sum all the amplitude locating in the hyperbola curve of diffraction, it assumes the circular shape of the object as a part of its diffracted hyperbola and turns circle the shape into the focused point. For the big rectangular pipe, the tradition full rectangular window cannot be well applied. The energy, entropy and varimax do show no extreme values that could relate to the velocity.
In our suggested workflow, the varimax diagram without reflection contribution can reveal the strongest varimax responding to the environment velocity. The reflection contribution can be defined from migrated GPR sections.

Conclusion
We have applied the traditional workflow and our workflow of velocity estimation by using varimax diagrams to the synthetic and real data. For small objects, the both workflows could provide the environment velocity. However, in the big object case our workflow can work better than the tradition one when our new varimax sections with a suitable window setting can catch diffraction contribution after removal of reflection contribution. Moreover, the categories of over-migration, correct migration, and under-migration is the stop condition of our workflow. Finally, energy diagram does not show its effectiveness in detecting environment velocity because the instability of energy extremes occurs with different window lengths.

Acknowledgments
This research is funded by Vietnam National University HoChiMinh City (VNU-HCM) under grant number C2019-18-08. I would like to thank Anh L. T. Ha for her support.