Feature Extraction and Automatic Material Classification of Underground Objects from Ground Penetrating Radar Data

. Ground penetrating radar (GPR) is a powerful tool for detecting objects buried underground. However, the interpretation of the acquired signals remains a challenging task since an experienced user is required to manage the entire operation. Particularly difficult is the classification of the material type of underground objects in noisy environment. This paper proposes a new feature extraction method. First, discrete wavelet transform (DWT) transforms A-Scan data and approximation coefficients are extracted. Then, fractional Fourier transform (FRFT) is used to transform approximation coefficients into fractional domain and we extract features. The features are supplied to the support vector machine (SVM) classifiers to automatically identify underground objects material. Experiment results show that the proposed feature-based SVM system has good performances in classification accuracy compared to statistical and frequency domain feature-based SVM system in noisy environment and the classification accuracy of features proposed in this paper has little relationship with the SVM models.


Introduction
For its rapid, continuous, nondestructive, efficient, convenient, and high-resolution properties, GPR is widely used in national defense, airport construction, railways, and highways [1][2][3][4][5].However, the interpretation of the large amount of acquired and stored GPR data requires a human operator with high skill and experience.In addition, the signals reflected from objects are corrupted with noise that cannot be completely removed.Hence, these problems lead to the growing request for the development of accurate and fast automated subsurface object detection and identification techniques in noisy environment.In recent years, a large number of scholars have proposed a variety of methods of automatic target recognition.
Torrione et al. used histograms of oriented gradients (HOG) features to detect landmines and suggest that other techniques from computer vision might also be successfully applied to target detection in GPR data [6].Phase profile was used in detection and location of targets underground because phase profile is a function of depth [7].A timefrequency feature extraction method based on Wigner distribution is proposed [3] and can detect and recognize diseases which have different compositions and shapes.Wahab et al. proposed a new hyperbola fitting technique to estimate the radius of buried utility (pipes and cables) [8].GPR signal explanation model is established based on the support vector machine and the dyadic wavelet transform (DyWT) theory in [9].It is applied in railway diseases detection and the results are effective.Ko et al. used principal component analysis (PCA) and Fourier coefficients as features to detect and identify landmines in different environments [10].Principal components from PCA, Fourier coefficients, and singular values from singular value decomposition are three features which are used to detect landmines in various burial conditions [11].Most of the above methods concentrate on the detection and localization of buried objects.However, very few works have been completed on the material identification of underground objects [12].
The basic motivation of proposing a novel feature extraction method is to improve automatic material classification accuracy in noisy environment and testify the results have little relationship with the SVM models.In this paper, the novel features, statistical features, and frequency domain features extracted from A-Scan data combined with SVM classifiers for material identification are tested.First of all, The synthetic GPR data after mean filtering is corrupted by different intensity of white Gaussian noise, which denotes different levels of noisy environment.Then, we test them on laboratory data in which we do not know the noise.From the experiments results we can see the good performances of the features proposed in this paper combined with SVM classifiers.

Forward Simulation
The evaluation of different feature-based approaches studied in this paper is firstly performed using different generated GPR data synthesized by using the electromagnetic simulator "GprMax2D" [13].GprMax2D was implemented using the FDTD numerical technique.An approximate solution for Maxwell's equations is directly obtained in the time domain by discretizing it in both space and time through an iterative process.
In our simulations, several parameters have to be set for the transmission and acquisition system.Parameters for FDTD simulation are listed in Table 1.
There are four models with the same parameters except the material of cylinders.As metal, air, stone, and polyvinyl chloride polymer (pvc) are common materials of objects in subsurface, in this paper the cylinders' materials embedded in the four models are the four materials, respectively.The parameters of the four materials of cylinders are listed in Table 2. To save computation time and storage space, only 320 traces that can display the target signal completely are taken and the other 162 traces are removed.The simulation results are shown in Figure 2. The results show that there is a hyperbolic curve in every B-Scan datum.However, the shape of the hyperbolic curve is different more or less.For example, the curvature of the four hyperbolic curves, grey value of hyperbolic curve, and the distance between vertices of hyperbolic curve are different.
If the shapes have big difference we can classify them artificially.On the contrary, we hardly classify them with manpower when the shapes are similar (e.g., Figures 2(c) and 2(d)).In addition, particularly difficult is their classification in noisy environment.Thus, it is important to find an automatic and high accuracy recognition method instead of manpower in noisy environment.

Feature Extraction
3.1.Flow Chart.In the process of feature extraction, we first use DWT to transform A-Scan data and extract approximation coefficients.Then, FRFT is applied to transform approximation coefficients into fractional domain.Finally, fractional domain-envelope curve is extracted and we extract feature to construct feature vector.The flow chart of feature extraction is shown in Figure 3.

Data Preparation.
As is known to all, echoes from ground surface are the strongest signals that can decrease the quality of GPR image and interfere with the identification of objects, but they can be removed by mean filtering [14].Hence, to get rid of the interference of echoes from ground surface we apply mean filter to every B-Scan datum in Figure 2. Figure 4 is the result of Figure 2(a) after mean filtering.From Figure 4 we can see the echoes from ground surface are removed and the hyperbolic curve is more obvious.Besides, the received data in practice contains a lot of other noise which cannot be removed completely and can be regarded as white Gaussian noise to some extent.As the data after the mean filtering hardly contain the noise, to simulate the reality we add different amplitude from 10 dBW to 50 dBW of white Gaussian noise to B-Scan data after mean filtering.Figure 5 shows that the 35 dBW white Gaussian noise is added to Figure 4.Meanwhile, the mean amplitude of object's signal is 30 dBW.

Extraction Approximation Coefficients.
In this paper, we use DWT extraction approximation coefficients.If the function being expanded is a sequence of numbers, like samples of a continuous function (), the resulting coefficients are called the discrete wavelet transform (DWT) of ().The DWT transform pair is defined as where (),   0 , (), and  , () are functions of the discrete variable  = 0, 1, 2, . . .,  − 1, () is the preparation A-Scan data in Section 3.2 and Figure 6 is one A-Scan trace of void that 10 dBW white Gaussian noise is added to Figure 4, () is wavelet function, and () is scaling function.Normally, we let  0 = 0 and select  (the length of the discrete samples of ()) to be a power of 2 (i.e.,  = 2  ) so that the summations are performed over  = 0, 1, 2, . . ., −1,  = 0, 1, 2, . . ., −1, and  = 0, 1, 2, . . ., 2  − 1.The transform itself is composed of  coefficients, the minimum scale is 0, and the maximum scale is −1.The coefficients defined in ( 1) and ( 2) are usually called approximation and detail coefficients, respectively.
In practice we select a wavelet from ready-made wavelets for a particular problem.Here we chose a wavelet called Daubechies wavelet.The approximation coefficients can be extracted from (1).
The fractional Fourier transform, a generalization of the Fourier transform (FT), serves as a useful and powerful analyzing tool in optics, communications, signal processing [16], and so forth.
The FRFT of a signal () ∈  2 (R) is defined as where () is approximation coefficients and the transform kernel is given by where   = √(1 − cot)/2,  = /2, and  is order of FRFT.Although  can take any real number, the period of FRFT is 4. Hence, we usually only consider  ∈ [0, 4].The  axis is regarded as the fractional domain.
In this paper, we only consider the case of  ∈ [0, 2] and set the interval of neighbor  as 0.05, since the definition can easily be extended outside the interval [0, ] by noting that  2 is the identity operator for any integer  and that the FRFT operator is additive in index; that is, The FRFT depends on a parameter  and can be interpreted as a rotation by an angle  in the time-frequency plane [17].Whenever  = 1, that is to say,  = /2, (3) reduces to the FT.
The fractional domain-envelope curve can be extracted:  = max ( (  ())) . (5) The envelope curve of the four materials can be seen from Figure 7.It is indicated that the different materials have big difference.Hence, it is possible to extract features from the envelope curve to classify material of underground objects.
In [17], we have known that FRFT can be interpreted as a rotation by an angle  in the time-frequency plane.That is to say, FRFT provides a signal comprehensive description of the whole process from time domain to frequency domain.In addition, with the  from 0 to 1, the FRFT display signal changes gradually from time domain to frequency domain of all variation characteristics. from 0 to 1 and  from 1 to 2 are symmetrical.The reason why different materials have different envelope curves may be that variation characteristics are different from time domain to frequency domain.Figure 8 shows the DWT-FRFT of air.

Feature Vector Construction.
In this paper, we will extract four features from envelope curve to construct feature vector.These features' calculation in MATLAB is as follows: (a) the maximum of envelope curve: (c) the variance of envelope curve: where   is the  th element in the vector  of length ; from Section 3.3 we can see  = 41 ( ∈ [0, 2] and set the interval of neighbor  as 0.05); (d) the kurtosis of envelope curve: where  4 () and  2 () are the fourth moment and second moment of , respectively.
The kurtosis is used to reflect the envelope curve of the top tip forms or flat level.Due to the sharpness of the envelope curves being different, the kurtosis is different.
According to the above, we can construct feature vector:

Experimental Results
SVM was designed originally for solving binary classification problem.However, to solve a multiclass SVM problem, two approaches can be used, either by creating several binary classifiers or by using a larger multiclass optimization problem.However, it is more expensive computationally to solve a multiclass optimization problem in one step than a binary problem using the same data size.In this paper, the major multiclass SVM methods, for example, one-againstone method, based on constructing binary classifiers, are used.It generates (−1)/2 classifiers; each classifier is trained on data from two classes.Many approaches can then be applied for material classification; the one used in this study is the voting approach using "Max Wins" strategy.The linear SVM works by finding an optimal hyperplane between two classes that maximizes the separating margin.For nonlinearly separable data, a kernel method issued to map the data in a higher dimensional feature space, that is, Φ() ∈ R  , to become linearly separable and apply the linear SVM algorithm.The optimal hyperplane ( * ⋅ Φ() +  * ), characterized by the weight vector  * (normal to the hyperplane), and the bias  * can be identified by solving the following optimization problem: where   is called slack variables used to account for nonseparable data and the parameter  is added to control the tradeoff between the slack variable penalty and the width of the margin.Using a Lagrange functional, this optimization problem can be reformulated in (12), where the Lagrange multipliers () can be calculated by using dual optimization: where  = ∑  =1     Φ(  ).The problem formulated in ( 12) is a convex quadratic optimization problem, and  can be obtained using the quadratic optimization problem solver.The final discriminant function is identified by  *  * and can be known, and hence, the support vector set representing the nonzero Lagrange multipliers can be identified.From the optimization problem formulated in (12), it is clear that the mapped data appear as an inner product only.From Mercer's theorem, it is known that, for any two points   and   and a certain mapping function Φ(), a kernel function can be used to evaluate the inner product of the mapped points without knowing the mapping function; for example, Φ(  )Φ(  ) = (  ,   ).Thus, the linear classifier is changed to a nonlinear classifier by substituting the inner product in the optimization problem by the kernel evaluation.One of the most common kernel functions used is the Gaussian function given where  parameter controls the Gaussian kernel width.

SVM Model Building and Training.
From Figure 2 we can obtain 320 traces (A-Scans) in each B-Scan datum, but not all of them are reflected from the cylinder embedded in the model.Hence, 40 traces from 153 to 192 which are echoes of the cylinder are selected in each model and the total number of useful traces in the four models is 160.All the selected A-Scan data should be prepared in advance in Section 3.2.
In training of the SVM classifier which is LIBSVM-Faruto Ultimate Version [18], we uniformly select 10, 20, 40, and 80 traces, respectively, from the 160 traces as training samples to obtain four SVM models.The others are testing data.First, statistical feature vector [12], frequency domain feature vector [19], and DWT-FRFT feature vector proposed in this paper are extracted, respectively.  denotes a feature vector.Then, a tag   ( = {1, 2, 3, 4}) is assigned to each feature vector to indicate the type of each material.Thus, each A-Scan datum can be expressed by a vector   = {  ,   }.Finally, the different number of training samples is used for training SVM classifier and we get different SVM models.
In addition, the SVM efficiency depends on the kernel's type selection, the kernel's parameters, and soft margin parameter .RBF is selected as the kernel function for the SVM model.Best combination of  and  is selected by a grid search with growing sequences of  and .

Automatic Classification Results and Discussion
. Now the trained SVM models are ready for classification of GPR data of different material objects.In this section, several tests are done to examine the classification accuracy of DWT-FRFT feature-based SVM system of different models for identifying different materials of underground objects.The classification accuracy can be computed by where  is the number of correct classifications and  is total test data.The performances of statistical features and frequency domain features in underground objects material identification using SVM classifiers are also studied and compared.
Figures 9(a)-9(d) are the classification accuracy of the three features using different trained SVM models in noisy environment, respectively.From this comparative study, it is inferred that no matter which one SVM model is used, the DWT-FRFT feature-based SVM system outperforms the other two on the whole.On the other hand, it is indicated that the DWT-FRFT feature is more robust against noise effects than the other two.However, one drawback is that the classification accuracy is lower than statistical feature-based SVM system when the noise intensity is about 25 dBW.This may be because, in feature extraction, only the approximation coefficients are extracted, while the useful information in detail coefficients is removed.Besides, the smaller the noise intensity is, the more the useful information in detail coefficients will be discarded.In addition, when the noise intensity is above 40 dBW, the classification accuracy will be below 50% and the credibility is not enough.Thus, it is not important which accuracy is high when the noise intensity is above 40 dBW.
As is known to all, the number of training samples should be as little as possible if the classification accuracy is similar.Figures 10(a)-10(c) show the results of performances comparison holding the same feature using different SVM models.
As can be seen from Figure 10(a), we can find that the classification accuracy curves of DWT-FRFT feature-based  different SVM models are intensive.It is indicated that the classification accuracy has little relationship with the number of training samples.That is to say, we can select fewer samples to train SVM and can obtain satisfactory results.Although Figure 10(a) is not as intensive as Figure 10(c), the biggest difference of classification accuracy is less than 10% except some one point (e.g., the noise intensity is 40 dBW).

Experiments with Laboratory Data.
In order to verify the performance of DWT-FRFT feature-based SVM system for material classification, experiments are realized on the laboratory data.The GPR system used in this work is LTD-2200 system (900 MHz, 1024 samples per scan, 190 traces, and trace spacing 0.015 m).
There are three different materials objects and they are buried in the same position of a 3 m × 9 m × 1 m bunker, respectively.The depth of objects is 0.3 m.The bunker is filled with uniform sand and the surface is flat.The size of the objects is similar (0.5-meter-long, 0.2-meter-thick).The three objects are copper, stone, and soil, respectively.The raw radargram of copper is shown in Figure 11.30 traces from 80 to 109 are reflected from object.
First of all, the raw data of three objects is preprocessed by mean filter and the result of radargram of copper is shown in Figure 12.From Figure 12 we can see the signal of object is more obvious and the echoes from ground surface are removed.
Then, the traces reflected from the three objects are selected and we can obtain 90 useful traces.Statistical feature, frequency domain feature, and DWT-FRFT feature are extracted from the traces, respectively.We uniformly select 9, 15, and 30 traces as training samples to obtain three SVM models and the others are testing data, respectively.The results of classification are shown in Table 3. From Table 3, it can be seen that compared to the other two featurebased SVM systems DWT-FRFT feature-based SVM system for material classification performs well in classification accuracy.In addition, the classification accuracy has little relationship with the number of training samples.

Conclusion
A novel DWT-FRFT feature-based SVM system for material classification of the underground objects from GPR data proposed in this paper is proved from the experiments of

Figure 1 :
Figure 1: One of the simulation models.
(b) the average value of envelope curve: Average = mean () ;

Figure 7 :Figure 8 :
Figure 7: The fractional domain-envelope curve of the four materials.

Figure 9 :
Figure 9: Classification accuracy of three features combined with SVM model.

Figure 10 :
Figure 10: The influence of SVM mode.

Table 2 :
Parameters of the four materials.

Table 3 :
Results of classification.