Georadar Data Processing by an Automatic Algorithm

B-algorithm was used for georadar data processing. The algorithm detects heterogeneities correlated with region boundaries and “regime” shifts in georadar data paths and constructs their distribution. The use of the algorithm in the MATLAB system simplifies and accelerates calculations, bringing the estimates obtained closer to optimum solutions. The algorithm, tested on model medium and winter lake georadar data, was shown to be sensitive to “regime” shifts and to be able to process data automatically. Key-Words: Georadar, boundary, heterogeneity, path, algorithm, “regime” shift.


Introduction
The current trend in the development of automated shallow-depth geophysical systems is an attempt to construct algorithms that simplify data processing and interpretation [1], e.g. a georadiolocation method to study subsurface media. This simple and convenient method is productively used for survey but it becomes less operative at a data processing stage, when an interpreter is used. Therefore, attempts should be made to find automated data processing algorithms, at least for some typical applications. When conducting engineeringsurveying and geological-geophysical studies, one common task is to delineate the boundaries of bodies or medium heterogeneity regions. The wave forms of georadar data contain the data required. Analyzing their image, the interpreter looks for the equiphase condition axes of the impulses reflected from layer boundaries, from which he calculates the time taken by a sounding impulse to travel from the source to the boundary and back to the signal receiver. Some calculation difficulties arise because the contact of the layers is diffuse, the parameters are uncertain, the signals are affected by noise, the impulses are overlapped and the data obtained are processed visually [2, 3]. The automatic algorithm single-way processes the paths of georadar data and displays the results, improves the performance of research, and reduces the impact of the human factor. The application of B-algorithm, which automatically reveals "regime" shifts in time data series, to process the georadar data paths of Balgorithm, is discussed in the present paper [4]. The B-algorithm gives a simplified schematic representation of the process by the model of "regimes", reveals the structure of process and the duration of the regimes, helps to establish causeeffect relations. Its application to georadar data processing is based on the correspondence of path "regime" shifts to medium heterogeneity region boundaries. To analyze the results of algorithm application, let us focus on estimating the shift moments of georadar data paths without recognizing the heterogeneities identified and estimating their electrophysical parameters, which is a separate problem [5]. To adapt B-algorithm constructed in the MATLAB system [6] to the problem with regard for its distinctive features, it was modified [7]. The operation of the algorithm was demonstrated on model data and the fragments of a georadar section for Lake Mustajärvi, Suojärvi District, Republic of Karelia, Russian Federation, obtained by profiling from the ice surface in winter.

Materials and Methods
An Oko-2 georadar with a central frequency of 400 MHz and a resolution of 0.15 m for free space was used for surveying the Mustajärvi Lake sequence. A similar georadar was used by L.L. Fedorova and K.O. Sokolov for physical modelling of an artificial medium [8]. The plot of the path of this work, digitized by the Digitizer software [9] and interpolated linearly, was used in this paper as an example of data processing. Georadar OKO is a commercially available geophysical device with a line of replaceable antenna units operating in the range of 35 to 2000 MHz, designed for ground, underwater and above-ground shooting. Mobile and compact OKO georadars allow to carry out non-destructive surveys of environment with high detail to the depth of up to several tens of meters. Georadar motion sensors, connected to a GPS navigator, provide an accurate reference to the terrain. Optical isolation by signal and information circuits is implemented in the instrument for noise immunity of measurements. Independent operation of the device in field conditions is provided with additional power supply system. Professional software contains basic options for georadar signal processing and results visualization. The algorithms of the Georadar-Expert package are similar to B-algorithm [10]. The BSEF algorithm analyses a backscattering field, representing a medium as a cross-section through the attributes of electrophysical characteristics. Controlling representation contrast, the algorithm makes it possible to vary the vertical resolution and depth of study. The B-detector algorithm detects medium layer boundaries using Canny's image recognition technology [11]. The technology uses procedures for signal smoothing, search for maximum gradients, noise suppression, double threshold filtration and tracing of heterogeneity regions. Bdetector combines high-speed data processing with high-precision boundary positioning. Wavelets, understood as functions localized in frequency and time, are used to analyze a signal on various scales and to reveal the moments of rapid changes and impulse overlapping regions [12]. It is easy to use Hilbert transform for the energy representation of signals in the form of an envelope [13]. A variety of methods, functions and scales, a large number of parameters and visual interpretation of results make it difficult to process, control and adjust georadar systems. R-algorithm, which recognizes "regimes" in climatological data series, is easier to use [14]. Operating a small number of parameters, the algorithm estimates a noise constituent, calculates a mean signal sampling value in the sliding window, compares it with the previous value and decides to shift the "regime». The wide-scale application of the аlgorithm is favoured by its free use in the Visual Basic system [15]. A wish to avoid false operation decreases the sensitivity of the algorithm to a "regime" shift. The boundary "regimes" it identifies are often poorly supported by available data. B-algorithm has a simple mathematical apparatus and a small number of parameters. The algorithm consecutively forms a model of path "regimes" at clusterization, comparison and optimization stages.
by varying their lengths. ni is the number of series members in a section and Аi,Сi is their mean value. Using the algorithm in the MATLAB computer mathematics system, one can use its functions, simplify and accelerate calculations, brings the solutions obtained closer to optimum and present the results in well-defined form. To automatically perform MATLAB functions, the parameters required are preset. The "mining" clusterization parameter, r ("subclust"), is preset by a number in the interval [0.1; 0.5]. p-parameter is the confidence level in the command «ttest2» of Student's test, which compares mean signal values on adjacent path sections. Default, its value is taken to be 0.05. The natural parameter b with + signs, added to section lengths, is responsible for the length variation range in the "patternsearch" command -a rapidly converging variety of a genetic algorithm. As B-algorithm with n-number of iterations is iterated, the smallest total deviation variance for the sections is achieved. Extra algorithm procedures perform the enumeration of georadar data paths and the construction of a "regime" shift image in the form of a binary matrix, in which 1 corresponds to the boundary and 0 to its absence. Comparison of B-and R-algorithms on the examples showed that B-algorithm reveals stable "regimes" of the processes, has high sensitivity to "regimes" shifts, distinguishes the shifts and the trends of the "regimes" [4].

Example 1
The example explains the correspondence of a "regime" shift and heterogeneity regions on a model medium consisting of a soil, frozen sand, concrete sequence (Fig.1а) [8]. When processing georadar data with wavelets, no reflections from the lower boundary of the concrete layer were detected because of impulse damping in the electroconducting frozen sand layer. B-and R-algorithms ( Fig. 1c and Fig.1d, respectively) were applied to one of the paths [8], the plot of which (Fig.1b) was digitized and interpolated linearly. The image of the "regime" shift of R-algorithm is not consistent with the model. Consistency is not achieved by varying parameters over a wide range. As window length decreases, the algorithm emits intensive impulses; as window length increases, the algorithm describes the path using one "regime". The inconsistency seems to be due to the low sensitivity of the algorithm to "regime" shift. B-algorithm recognized a "regime" shift at all model medium boundaries, including the lower boundary of the concrete layer, representing the "regime" shift by one count of the path. The last "regime" shift was estimated less precisely than other shifts presumably because of signal damping. However, this example shows that one can work with dying-out signals.

Example 2
The thickness and structure of ice sheets are controlled when ships move in ice-bound water bodies, when ice passages are constructed and when building operations in permafrost areas are carried out [16]. It is promising to survey huge ice areas with a georadar installed on an aircraft and to process data using algorithms that identify ice types, water and air lenses and fractures [17]. Ice sheets are suitable for georadar studies because of  the high quality of a signal, a small number of images to be identified and the simplicity of checking the results by drilling. B-algorithm was used for processing fragment I of the ice-covered georadar section of Lake Mustajärvi (Fig. 2). Fig. 3 shows the results of processing of the first path of fragment I by B-algorithm at various rparameter values and the structure of the ice cover "regime" of the fragment. Fig. 3 d shows the structure as consisting of unconsolidated surface ice layers, consolidated ice and lower more poorly-ordered water-ice mixture. Fig. 4а shows near-shore fragment II of the georadar section (Fig. 2). Fig 4b and 4c show the results of its processing by Canny's technology and Hilbert transform using MATLAB «edge» and «hilbert» commands. A scheme of the fragment's «regimes» recognized by B-algorithm is shown in Fig.4d. Canny's technology was used to reveal ice sheets, water column heterogeneities and the mode of occurrence of bottom rocks. Hilbert transform was used to estimate ice cover thickness and the positions of water column and bedrock heterogeneities. B-algorithm also revealed ice layers, the water column pattern and the gentle descent of the lake bottom bedrock. The results of the algorithms are hard to compare in detail, but the pattern of the "regime" shift of B-algorithm correctly shows the positions of the major heterogeneities of the medium studied. The positions of "regime" shifts in the layers revealed by B-algorithm are blurred with depth. Fig.3 The path of the first picket of fragment I and its ten-fold image of the "regime" of B-algorithm at b=4, n=10 and r = 0.24, 0.36, 0.3, respectively (a,b,c) d -shift structure of ice cover "regimes" for parameters (с). This could be due to impulse damping and the algorithm parameters chosen. To better show "regime" shifts on a long path, the r-parameter value should be small enough to be consistent with the smaller section length at an algorithm clusterization stage. At preset b-parameter values this can decrease the estimation accuracy of "regime" shift moments. This problem could be solved using aposterior control of strengthening along the path, analyzing the enveloping path and comparing variances in its sections. Fig. 4. Fragment II of the georadar section of Lake Mustajärvi (а) and the results of its processing by Canny's technology (b) and Hilbert transform (c), an image of its "regime" shift provided by B-algorithm with the following parameters: r=0.1, b=4, n=15 (d).

Example 3
Another feature, revealed while processing long paths, is time-consuming. The common way of making them shorter is to reject consecutive processing and to modify the algorithm by making calculations parallel.