GPR DATA PROCESSING IN AUTOMATIC MODE
Authors: R.R. Denisov, Ph.D. Sciences V. V. Kapustin
Reviewer: Doctor of Physics and Mathematics Of Sciences M.L. Vladov
The article was published in
Scientific and technical journal
"Geophysics" No 4 for 2010, pages 7680
ISSN 168145688
Abstract. Authors consider the way of approximate solution of the inversion task for continuous georadar profiling. The backscattering electromagnetic field is used to obtain the section of electrical and physical subsoil properties. Method of stable solution is based on building nonlinear spatial filters that are applied to number set of subsurface permittivity (or electromagnetic waves velocity). The suggested approach allows building on electronic system of georadar data processing that works efficiently in the complex diffusing medium.
Key words: georadar, electronic data processing system, algorithm, subsurface permittivity.
The application of the georadar method for the solution of engineering geological and geotechnical problems is now acquiring an ever wider scale. This is explained by the simplicity and convenience of fieldwork and high productivity method. At the same time, a number of problems arise with processing large data sets and interpreting the resulting material. When using georadar for the study of highways and railways, aerodromes, while working in water areas, the length of georadar profiles can be tens or even hundreds of kilometers. In most cases, rapid receipt of the final material is a necessary condition for conducting georadar work [2].
The existing software for processing georadar data allows for a number of standard tasks to use readymade processing graphs that improve the appearance of raw data [3]. When solving problems of detection in the plan of local areas of electrophysical inhomogeneities (watering areas, clay lenses, voids, cracks on asphalt or concrete cover, etc.), the automation of this process can be partially solved by the methods of attribute analysis. This problem is analogous to the problem of profiling used in various geophysical methods. As the attributes, amplitude, frequency, phase recording characteristics, attenuation, and a number of others can be used.
The problem of determining the spatial position of electrophysical boundaries and inhomogeneities (the probing problem) is more complicated in the sense of using automation tools. This is due to the fact that the automated determination (picking) of electrophysical boundaries, in most cases, requires the interactive participation of the geophysicistinterpreter or the use of complex selflearning algorithms. It is quite natural that to solve most typical problems related to the study of engineeringgeological conditions, it would be desirable to have the possibility of an operative automated processing that makes it possible to obtain material convenient for interpretation.
The problem of constructing a section of a subsurface medium, as a rule, is solved by the method of visual analysis of the wave pattern of the radargram and tracing the inphase axes, which correspond to the boundaries of the layers with different electrophysical parameters [5]. This method does not cause difficulties if the thickness consists of regions with sufficient dielectric contrast, at the boundaries of which the dielectric constant varies discontinuously, and these boundaries are confidently observed on the radargram.
However, it often happens that the characteristics of components, for example, the ground, change smoothly, without a jump, due to the diffuse nature of the contact of neighboring layers or for some other reason. In this case, the inphase axes that correspond to the desired boundaries on the radar picture are difficult to detect, or these boundaries are not visible at all.
In such conditions, the only solution is to manually determine the velocity characteristics of the sections of the section from the diffracted waves found on the radargram, which were formed as a result of reflections from point subsurface objects, and to combine regions with similar velocities into layers. The reliability of the cut thus obtained is determined by the number of detected reflections, the uniformity of their distribution over the area of the radargram, the degree of preparation of the geophysicistinterpreter, and the human factor (inattention, fatigue, etc.). The laboriousness of this method does not allow processing large volumes of georadar data for an acceptable time. Thus, the advantage of the georadiolocation method in the field of high productivity is reduced to nothing.
Taking into account the existing features of the subsurface radar method: the relatively low depth of research, the high resolving power of the method, the complexity and heterogeneity of the upper part of the section, it is quite natural that the field of backscattering of electromagnetic waves (BSEF) can be used to determine the electrophysical structure of the medium [4]. Automation of the process of recognition of diffracted reflections and determination of kinematic characteristics of waves with subsequent construction of a section allows to expand the possibilities of georadiolocation.
Let's consider the basic principles that should be incorporated into the algorithm for analyzing the backscattering field of the GPR waves. The wave field on the radargram is of a rather complex nature and usually consists of the following components: single reflections; multiple reflections; backscattering (diffraction); air reflections; interference. In this connection, the first problem is the separation of the backscattering field from the set of observed data. This problem can be partially solved using standard processing means  frequency and spatial filtering [5]. However, complete separation of the fields can not be achieved. Therefore, after preliminary processing of the data, the resulting numerical array is conventionally considered as a backscattering field that can be subjected to a transformation (such as hyperbolic migration) to calculate the permittivity or velocity V.
During the calculation, the algorithm excludes values that do not fall into a predetermined range 1<<<< 2 and computes in the vicinity of each determination point of the value a number of dynamic parameters  frequency f, damping d, amplitude A, etc. As a result, the vector field [, f, d, A, ...... Mn] of dimension n equal to the number of dynamic parameters can be determined. Further, using a priori or statistically estimated connections between the dielectric constant and one or more dynamic parameters, a nonlinear spatial filter mask can be constructed to suppress the values determined by multiple reflections, air reflections and interference.
Automated processing is performed according to the following algorithm (for example, the dielectric permittivity attribute):
 Detection on the radargram of points through which the diffraction curves pass;
 Determination of kinematic and dynamic characteristics of signals in the vicinity of selected points;
 Culling points by the speed criterion. In the future, the points remaining after the culling will serve as nodal points when creating the section;
 Recalculation of the wave velocities into effective values of the real part of the complex relative permittivity ;
 Correction of values by frequency attribute or group of attributes;
 Creation of a section of the permittivity of the medium in the time domain by interpolation of values by reference points;
 Recalculating the time scale of the cut into a depth scale in accordance with the values
As a result of the first two steps of the algorithm, a set of points that have coordinates and attributes in the form of values of speed and frequency is formed (the option of using one attribute is considered). Further, there is a rejection of points according to the speed criterion  for example, points with velocities above 30 cm / ns (electromagnetic wave velocity in vacuum) and below 3.33 cm / ns (minimum electromagnetic wave velocity in water) are not consideration as nodal points of the section and are deleted.
In Fig. 1 is a radargram showing the positions of the nodal points after the rejection by red markers. High density of points allows to obtain a sufficiently detailed section.
Fig.1 An example of a radargram with the position of nodal points, after the rejection according to the speed criterion.
After sorting the points, the values of the velocity are converted to the values of the dielectric constant, since in the future the kinetic characteristics of the waves and the electrophysical properties of the medium will be investigated. The most important stage of processing is the correction of the values by the frequency (or some other) attribute. It is required to correct possible errors in determining the speed of electromagnetic waves and to avoid distortions caused by the effect of multiple reflections, air reflections and other disturbances.
This procedure is performed, as noted above, using adaptive spatial filtering. The filter mask (decision rule) can be constructed on the basis of one or several empirically established connections between the dielectric permittivity matrix point of the medium and one or more dynamic attributes. Consider a variant of constructing such a filter based on one of the properties of the medium under study.
As is known, the attenuation of the probing pulses of a georadar is frequency dependent  highfrequency components of electromagnetic waves decay more rapidly than lowfrequency components. Along with this, the humidity, on which the dielectric permittivity of the soil depends first of all, also varies with depth, as a rule, it increases, and accordingly the value also increases. If the mean values from the set {} are determined for the values of the central frequencies of the pulses determined in the process of processing, and the dependence is plotted as a point cloud, then the regression equation can be used as a correction function (as a filter mask). As a result of the application of this function, the reliability of interpretation of the section significantly increases. Figure 2 shows the correcting function for the radargram obtained on the surface of frozen soils.
Fig. 2 Correction function obtained at the development site of the seasonal freezing layer.
The horizontal axis is the central frequency in the vicinity of the node point, the vertical axis is the mean values .
It should be noted that the geophysicistinterpreter, based on a priori information on the structure of the cut, can introduce his own corrective function, which more accurately reflects the relationship of and frequency for a given type of section.
The final stage of creating a section of dielectric permittivity is the interpolation of values over nodal points by methods known in computational mathematics (Delaunay triangulation, for example). Thus, the input of the automated processing program is fed with a numerical matrix of amplitudes of the reflected signals obtained by georadar profiling, an array of the same size is formed at the output, but with values for each point of the twodimensional space of the section. It is most informative to visualize the matrix of values in the form of a contour plot with areas filled in according to the chosen color map.
In Fig. 3 shows the result of the automated processing of the profile obtained by a georadar with a central frequency of 300 MHz probe pulses (studies of the characteristics of the ground in the permafrost region).
Fig. 3 GPR profile obtained by probing the soil in the permafrost region and the result of automated processing in the form of a section of the permittivity. Unlike the raw radargram, the degree of freezing of various sections of the subsurface medium can be observed on the section.
The automated creation of a section of dielectric permittivity has the following advantages over traditional methods of processing GPR data:
 The depth of georadiolocation studies is increasing  the method of analysis of BSEF has a high noise immunity and works well enough in the noise area of the radargram [Fig. 4];
 The information content of the research is increasing  the boundaries of the section are fixed where there is not enough jump in dielectric permittivity to form the commonmode axes of signals, which are reflections from the boundaries of the layers. A change in the permittivity within the layer is also clearly visible [Fig. 5];
 The processing speed of the field material is significantly increased, which is important for the constantly increasing volumes of georadar work, especially in the road and railway industries;
 The scope of georadar is expanding. Based on the permittivity profiles, it is possible to construct moisture distributions along the section [Fig. 6], frequency distributions or use, depending on the task, other attributes;
 Minimized the impact of the socalled human factor.
In conclusion, it should be noted that the proposed method is not a panacea for solving all problems of georadiolocation, but is one of the variants of the approximate solution of the inversion problem, which has applicability limits.
Fig. 4. GPR profile that crosses the landslide area and the result of automated processing in the form of a section of the dielectric constant. As a result of processing, the body of the landslide and the possible slip boundary are determined.
Fig. 5. The result of the automated processing of the GPR profile crossing the road with asphalt concrete. Based on the analysis of changes in the values of dielectric permeability within the structural layers, it becomes possible to diagnose, at a higher level, the state of road clothes and underlying soils.
Fig. 6. GPR profile along the sandy embankment of the road. As a result of the automated processing, a section has been constructed, according to which it is possible to judge the distribution of humidity in the area of the culvert that is at a horizontal mark of 234 meters.
References
 Gonzalez R., Woods R., Eddins S., 2006, Digital image processing in the MATLAB environment: Technosphere.
 Kapustin VV, Strochkov Yu.A., 2008, Some features of processing georadar data in the study of building structures: Exploration and protection of mineral resources, 1, 2225.
 Kapustin VV, 2005, Additional possibilities of computer processing of georadar and seismic data: Exploration and conservation of bowels, 12, 2631.
 V.B. Levyant, I.B. Petrov, S.A. Pankratov, 2009, Study of the characteristics of longitudinal and exchange waves of backscatter response from fractured reservoir zones: Seismic survey technologies, 2, 311.
 Starovoitov AV, 2008, Interpretation of georadar location data: Moscow, Izdvo, Moscow State University.
