Particle Center Supported Plane for Subsurface Target Classiﬁcation based on Full Polarimetric Ground Penetrating Radar

: The subsurface target classiﬁcation of ground penetrating radar (GPR) is a popular topic in the ﬁeld of geophysics. Among the existing classiﬁcation methods, geometrical features and polarimetric attributes of targets are primarily used. As polarimetric attributes contain more information of targets, polarimetric decomposition methods, such as H-Alpha decomposition, have been developed for target classiﬁcation of GPR in recent years. However, the classiﬁcation template used in H-Alpha classiﬁcation is preset depending on the experience of synthetic aperture radar (SAR); therefore, it may not be suitable for GPR. Moreover, many existing classiﬁcation methods require excessive human operation, particularly when outliers exist in the sample (the data set containing the features of targets); therefore, they are not efﬁcient or intelligent. We herein propose a new machine learning method based on sample centers, i.e., particle center supported plane (PCSP). The sample center is deﬁned as the point with the smallest sum of distances from all points in the same sample, which is considered as a better representation of the sample without signiﬁcant effect of the outliers. In this proposed method, particle swarm optimization (PSO) is performed to obtain the sample centers; the new criterion for subsurface target classiﬁcation is achieved. We applied this algorithm to full polarimetric GPR data measured in the laboratory and outdoors. The results indicate that, comparing with support vector machine (SVM) and classical H-Alpha classiﬁcation, this new method is more efﬁcient and the accuracy is relatively high.


Introduction
Ground penetrating radar (GPR) is a type of electromagnetic technique for the detection of subsurface targets.It has been applied widely to various fields such as engineering, archaeology, hydrology, and land mine identification.As underground object detection is always involved in these fields, the subsurface target classification of the GPR becomes a popular topic in the geophysical field.It is aimed to classify different types of targets based on their attributes.Among the existing methods, geometrical features and polarimetric attributes are primarily used in classification.For example, Ikechukwu K. Ukaegbu et al. combined GPR and a gamma ray detector to estimate the nonintrusive depth of buried radioactive wastes [1]; Wentao Li et al. applied a randomized Hough transform to achieve the automatic recognition of tree roots [2]; Byeongjin Park et al. used instantaneous phase analysis of GPR data to perform underground object classification [3]; Xuan Feng et al. applied migration for the detection of underground objects [4][5][6].However, traditional GPR measurements can only yield a small part of information for underground objects.As the full polarimetric data contains much more information of targets [7], full polarimetric GPR has been developed to detect underground objects in recent years [8,9].Polarimetric decomposition [10,11] is a method for obtaining polarimetric attributes from full polarimetric data, in which H-Alpha decomposition [5] and Freeman decomposition [12,13] have been applied to GPR.
H-Alpha decomposition is a method based on the Kennaugh matrix, which was proposed by Cloude and Pottier in 1997 [14].This method has been applied widely to synthetic aperture radar (SAR), as well as to target classification in GPR [5] in recent years.It uses entropy H and angle Alpha to perform the classification.The entropy H is a measure of randomness of scattering mechanisms, and the angle Alpha characterizes the scattering mechanism [15].This classification method depends on the classification template called the H-Alpha plane [5] which contains nine zones, and the physical scattering characteristic associated with each zone provides information for identification.This distinctive advantage, unfortunately, is offset by preset zone boundaries in the H-Alpha plane [15].The determination of classification boundaries depends on the experience of SAR.While some differences in measurement exist between GPR and SAR, such as the medium and distance of propagation for electromagnetic waves, observation methods, and noise characteristics, the template used in SAR may not be suitable for GPR.In this study, we removed the classical H-Alpha template and studied the characteristics of H and Alpha data; subsequently, we propose a new machine learning method to obtain the new classification criterion.
Particle swarm optimization (PSO) is a new evolutionary algorithm proposed by J. Kennedy and R. Eberhart in 1995 [16].This algorithm, which is similar to simulated annealing (SA), starts with a random solution and applies iterations to obtain the optimal solution.Unlike other methods, this method uses particles to represent the random solution.In each iteration, every particle will renew its position based on the optimal position that it has bypassed (local optimal solution) and the optimal position all the particles have bypassed (global optimal solution) [17,18].It is in fact a method to solve the extreme value of an objective function; therefore, we can change the objective function called the fitness function to solve many different problems.The advantages of this algorithm is fast convergence, high-quality solutions, and fewer parameters required for calculation.In recent years, PSO has been applied to full waveform inversion [19,20], magnetotelluric sounding data inversion [21], and unsupervised clustering analysis [22][23][24][25].
Machine learning has been applied widely in many types of disciplines.It primarily focuses on designing algorithms that can obtain laws from a number of samples which are data sets containing the features of targets, and the laws can be used to predict the characteristics of unknown data.In recent years, many machine learning methods have been applied in the GPR field.For example, Xavier Núñez-Nieto et al. applied logistic regression and neural network (NN) techniques to automated landmine and UXO detection [26]; Tao Liu et al. used neural networks to inverse GPR data [27]; Haoqiu Zhou et al. combined full polarimetric GPR and support vector machine (SVM) data for subsurface target classification [28]; Minghe Zhang et al. used freeman decomposition and random forest (RF) to perform underground object detection [29].
In this article, we define the sample center as the point with the smallest sum of distances from all points in the same sample, the sample center is considered a better representation of the sample without significant effects of the outliers.Subsequently, particle center supported plane (PCSP) is proposed to perform the classification.
This article is organized as follows.Section 2 introduces the theory of H-Alpha decomposition, and details the proposed the technique, i.e., particle center supported plane (PCSP).In Section 3, the PCSP is applied to full polarimetric GPR data of four types of targets measured in the laboratory.Section 4 verifies the feasibility of PCSP with real full polarimetric GPR data of underground pipes measured outdoors.In Section 5, an additional experiment is performed for three pipes of different depths, diameters, and materials; subsequently, classical H-Alpha classification and support vector machine (SVM) analysis are performed to compare with the PCSP method.

H-Alpha Decomposition
In full polarimetric GPR, the Sinclair matrix [S] is applied to describe the nature of target scattering.The [S] matrix [13] can be represented as follows: [S] = S HH S HV S V H S VV (1) where S XY represents the GPR data that is collected with an X type polarimetric receiving antenna and Y type polarimetric transmitting antenna, H represents horizontal polarimetric antenna, V represents vertical polarimetric antenna.According to the reciprocity theorem, S HV is equal to S VH .The scattering vector k s is defined as follows [4]: Coherency matrix [T] and its parameterization [4] is shown as (3), where λ 1 , λ 2 , λ 3 , represent eigenvalues, and [U 3 ] is composed with eigenvectors [4], [U 3 ] is shown as ( 4): where e 1 e 2 e 3 are the eigenvectors.The parameterization of a 3 × 3 unitary [U 3 ] matrix in terms of column vectors with different parameters α, β, δ and γ, which are the parameters of the dominant scattering mechanism, is made so as to enable a probabilistic interpretation of the scattering process.Finally, the mean scattering angle α and entropy H [4] can be calculated as follows: where p 1 , p 2 , p 3 are the false probabilities [4] calculated as follows: Herein, H and α are the parameters used to perform the classification.

Particle Center Supported Plane
We define the sample center as the point with the smallest sum of distances from all points in the same sample.Typically, some outliers will be present in the sample.These outliers will affect the training result of the classification significantly (Figure 1a), but the effect on the sample center is small.For example, in the quadrilateral ABCD of Figure 1b, the diagonal intersection O is the point with the smallest sum of distances from four points A, B, C, D; therefore, point O is the sample center of points A, B, C, D. However, wherever the outlier C is at the extension line of AC, the sample center O will never change.Therefore, the sample center is considered a better representation of the sample without the significant effects of the outliers.In this research, we propose a new machine learning method based on sample centers, i.e., particle center supported plane (PCSP), to perform the classification.The sample centers are obtained by particle swarm optimization (PSO).For PSO, particles are used to represent the random solutions of the optimization problem.The coordinate and velocity of mth particle at kth iteration can be represented as follows [22]: where M represents the number of particles, and k is the time of iteration, x where ztni represents the nth training sample point in ith type; t represents "training"; Ni is the number of sample points in ztni; P represents the number of types.
We want the sum of the distances between the particles and all the points in the same sample to be minimum.Therefore, the objective function (fitness function) is defined as follows: From the first iteration to the kth iteration, the minimum value of the fitness function for mth particle is known as the local optimal solution of mth particle at kth iteration, i.e., B k m , m=1,2,...,M; the minimum value of the fitness function for all the particles is known as the global optimal solution at kth iteration, i.e., B k g [22].For all particles, we use the formulas below to renew their positions and velocities [22]: where k represents the time of iteration; ω is the inertia weight; r1, r2 are the random numbers between 0 to 1; c1, c2 are the learning factors.In fact, for the right part of equation (11), the first item represents the impact of the current position, the second item represents the impact of the current optimal position of this one particle, and the third item represents the influence of current optimal position of all particles.When the fitness function becomes minimum, the particles are at the center In this research, we propose a new machine learning method based on sample centers, i.e., particle center supported plane (PCSP), to perform the classification.The sample centers are obtained by particle swarm optimization (PSO).For PSO, particles are used to represent the random solutions of the optimization problem.The coordinate and velocity of mth particle at kth iteration can be represented as follows [22]: where M represents the number of particles, and k is the time of iteration, x k mH and x k mα correspond to H and α coordinates, respectively.The initial position x 0 m and initial velocity v 0 m are random.In this problem, the sample contains H and α, therefore, the training sample can be represented as follows: where z tni represents the nth training sample point in ith type; t represents "training"; N i is the number of sample points in z tni ; P represents the number of types.We want the sum of the distances between the particles and all the points in the same sample to be minimum.Therefore, the objective function (fitness function) is defined as follows: From the first iteration to the kth iteration, the minimum value of the fitness function for mth particle is known as the local optimal solution of mth particle at kth iteration, i.e., B k m , m = 1,2, . . .,M; the minimum value of the fitness function for all the particles is known as the global optimal solution at kth iteration, i.e., B k g [22].For all particles, we use the formulas below to renew their positions and velocities [22]: where k represents the time of iteration; ω is the inertia weight; r 1 , r 2 are the random numbers between 0 to 1; c 1 , c 2 are the learning factors.In fact, for the right part of Equation (11), the first item represents the impact of the current position, the second item represents the impact of the current optimal position of this one particle, and the third item represents the influence of current optimal position of all particles.When the fitness function becomes minimum, the particles are at the center of the sample, the coordinate of sample center is B g , subsequently, we can obtain the sample centers of all types of samples: B g = B g1 , B g2 , . . ., B gP (14) where P represents the number of classes.
After the sample centers of different samples are achieved from PSO, a classification law can be constructed based on these sample centers.Figure 1a indicates that the outliers can affect the slope of the classification boundary significantly.In this method, straight lines are used as the classification boundary, and the slope of every line is determined by the sample centers.Subsequently, the intercepts are calculated by the classification accuracy for training samples.
For all the training samples {z tn1 , z tn2 , . . .,z tnP } and their sample centers {B g1 , B g2 , . . ., B gP }, where P is the number of classes, we use {z tni , z tnj } and their sample centers {B gi , B gj } as examples.In Figure 2a, red points and purple points belong to the ith type, green points and blue points belong to the jth type.A critical point Q ij between two sample centers is calculated based on the classification accuracy of the boundary Line ij .The Line ij is through the critical point Q ij and is perpendicular to the connection line of B gi and B gj ; the equation of Line ij is as follows: where z = (H, α) is the point on the flat.where P represents the number of classes.
After the sample centers of different samples are achieved from PSO, a classification law can be constructed based on these sample centers.Figure 1a indicates that the outliers can affect the slope of the classification boundary significantly.In this method, straight lines are used as the classification boundary, and the slope of every line is determined by the sample centers.Subsequently, the intercepts are calculated by the classification accuracy for training samples.
For all the training samples {ztn1, ztn2,...,ztnP} and their sample centers {Bg1, Bg2,..., BgP}, where P is the number of classes, we use {ztni, ztnj} and their sample centers {Bgi, Bgj} as examples.In Figure 2a, red points and purple points belong to the ith type, green points and blue points belong to the jth type.A critical point Qij between two sample centers is calculated based on the classification accuracy of the boundary Lineij.The Lineij is through the critical point Qij and is perpendicular to the connection line of Bgi and Bgj; the equation of Lineij is as follows: where z = (H, α) is the point on the flat.The points that are on the same side of Lineij with Bgi are considered to belong to the ith type.The accuracy of Lineij for classifying ith type and jth type are as follows, respectively: where Ni ((ztni -Qij)•(Bgi -Qij) > 0) represents the number of points in the ith sample which is on the same side of Lineij with Bgi.Subsequently, the points in the ith sample that are far away from the sample center Bgi are considered as outliers; therefore, Accuracyi and Accuracyj do not have to be 1.We chose Accuracyi,j = 0.8 as the threshold.The objective function is defined as follows: The points that are on the same side of Line ij with B gi are considered to belong to the ith type.The accuracy of Line ij for classifying ith type and jth type are as follows, respectively: ) where represents the number of points in the ith sample which is on the same side of Line ij with B gi .Subsequently, the points in the ith sample that are far away from the sample center B gi are considered as outliers; therefore, Accuracy i and Accuracy j do not have to be 1.We chose Accuracy i,j = 0.8 as the threshold.The objective function is defined as follows: Accuracy i,j ≥ 0.8 Subsequently, the optimal Q ij as well as the particle center supported plane Line ij can be determined by solving the maximum value of ( 18) under the condition of (19).
By repeating the above steps, all the boundary lines between every two samples of {z tn1 , z tn2 , . . .,z tnP } can be determined.Subsequently, classification can be performed for the new data z cn where c represents "classified".The discriminant function for the ith type with the other jth type is as follows: where j = 1,2, . . .,P, j = i.The complete discriminant function of the ith type is: If F i (z cn ) = 0, z cn belongs to the ith type.As Figure 2b shows, if the new sample point is on the same side of the red line with B g1 , it belongs to the type of B g1 .

Data Processing Flow Chart
The Data processing flow chart of the proposed method is in Figure 3. Firstly, the H-Alpha decomposition is performed to obtain the H and α data for training; subsequently the particle center supported plane (PCSP) is performed to obtain the classification boundaries; finally, the new H-Alpha data is applied to the PCSP to obtain the classification result.Subsequently, the optimal Qij as well as the particle center supported plane Lineij can be determined by solving the maximum value of ( 18) under the condition of (19).
By repeating the above steps, all the boundary lines between every two samples of {ztn1, ztn2,...,ztnP} can be determined.Subsequently, classification can be performed for the new data zcn where c represents "classified".The discriminant function for the ith type with the other jth type is as follows: where j=1,2,...,P, j ≠ i.The complete discriminant function of the ith type is: If Fi (zcn) = 0, zcn belongs to the ith type.As Figure 2b shows, if the new sample point is on the same side of the red line with Bg1, it belongs to the type of Bg1.

Data Processing Flow Chart
The Data processing flow chart of the proposed method is in Figure 3. Firstly, the H-Alpha decomposition is performed to obtain the H and α data for training; subsequently the particle center supported plane (PCSP) is performed to obtain the classification boundaries; finally, the new H-Alpha data is applied to the PCSP to obtain the classification result.

Data Analysis of Typical Targets in the Laboratory and Outdoors
In this section, full polarimetric GPR data measured in the laboratory was applied to the proposed method to obtain the classification boundaries and verify its feasibility.

Full polarimetric GPR Measurement
The full polarimetric measurements of four types of targets (Figure 4) were performed in the laboratory.The survey lines were right on the targets, the number of measurement points was 101, the measurement point interval was 1 cm, the frequency ranges between 800 MHz and 4000 MHz, and the number of sampling points was 1024.The four typical targets are as follows: metallic sphere, representing Bragg surface scattering; metallic cylinder, representing linear body scattering; metallic dihedral, representing double-bounce scattering; and metallic multibranch, representing multiple scattering.The diameter of the sphere was 15 cm.The length of the cylinder was 31.5 cm, the

Data Analysis of Typical Targets in the Laboratory and Outdoors
In this section, full polarimetric GPR data measured in the laboratory was applied to the proposed method to obtain the classification boundaries and verify its feasibility.

Full polarimetric GPR Measurement
The full polarimetric measurements of four types of targets (Figure 4) were performed in the laboratory.The survey lines were right on the targets, the number of measurement points was 101, the measurement point interval was 1 cm, the frequency ranges between 800 MHz and 4000 MHz, and the number of sampling points was 1024.The four typical targets are as follows: metallic sphere, representing Bragg surface scattering; metallic cylinder, representing linear body scattering; metallic dihedral, representing double-bounce scattering; and metallic multibranch, representing multiple scattering.The diameter of the sphere was 15 cm.The length of the cylinder was 31.5 cm, the diameter was 5 cm.The dihedral was made of two plates at a 90 • angle whose length and width were 35 cm and 20 cm, respectively.The length of the multibranch scatterer was approximately 40 cm.The four targets were immersed in a dry sand trough whose length, width and depth were 2.53 m, 2.53 m, and 0.85 m, respectively; the wave velocity of the sand was approximately 0.2 m/ns.The buried depths for the top of the sphere and cylinder were 8 cm; the buried depth for the intersection of two plates was 20 cm; the buried depth for the trunk of the multibranch scatterer was 10 cm.
In the measurement, a vector network analyzer and a Cartesian coordinate robot (Figure 5) with different antenna combinations (Figure 6) were used for data collection, the interval between feeding points of transmitting and receiving antennas was 8 cm for all three polarimetric modes with the same axis direction.The data of the sphere, cylinder, dihedral, and multibranch are shown in cm.The four targets were immersed in a dry sand trough whose length, width and depth were 2.53 m, 2.53 m, and 0.85 m, respectively; the wave velocity of the sand was approximately 0.2 m/ns.The buried depths for the top of the sphere and cylinder were 8 cm; the buried depth for the intersection of two plates was 20 cm; the buried depth for the trunk of the multibranch scatterer was 10 cm.
In the measurement, a vector network analyzer and a Cartesian coordinate robot (Figure 5) with different antenna combinations (Figure 6) were used for data collection, the interval between feeding points of transmitting and receiving antennas was 8 cm for all three polarimetric modes with the same axis direction.The data of the sphere, cylinder, dihedral, and multibranch are shown in Figures 7-10, respectively.cm.The four targets were immersed in a dry sand trough whose length, width and depth were 2.53 m, 2.53 m, and 0.85 m, respectively; the wave velocity of the sand was approximately 0.2 m/ns.The buried depths for the top of the sphere and cylinder were 8 cm; the buried depth for the intersection of two plates was 20 cm; the buried depth for the trunk of the multibranch scatterer was 10 cm.
In the measurement, a vector network analyzer and a Cartesian coordinate robot (Figure 5) with different antenna combinations (Figure 6) were used for data collection, the interval between feeding points of transmitting and receiving antennas was 8 cm for all three polarimetric modes with the same axis direction.The data of the sphere, cylinder, dihedral, and multibranch are shown in Figures 7-10, respectively.

H-Alpha Decomposition for GPR Data
After the collection for the full polarimetric GPR data of four types of targets, H-Alpha decomposition is applied to obtain the H and α data.Subsequently, the data of four targets are divided into two equal parts: the training data and the testing data.Figure 11a shows the spatial distribution of the training H and α data of four types of targets.

Particle Center Supported Plane
Subsequently, PSO is performed for the H and α data to obtain the sample centers of the four types of targets.Twenty particles are used in the calculation (M = 20), and we chose w = 0.5, c1 = 1.5, c2 = 2.5; when the relative error between the kth iteration and (k+1)th iteration is under 0.01, the calculation will stop.Figure 11b shows the sample centers of the training data.

H-Alpha Decomposition for GPR Data
After the collection for the full polarimetric GPR data of four types of targets, H-Alpha decomposition is applied to obtain the H and α data.Subsequently, the data of four targets are divided into two equal parts: the training data and the testing data.Figure 11a shows the spatial distribution of the training H and α data of four types of targets.

Particle Center Supported Plane
Subsequently, PSO is performed for the H and α data to obtain the sample centers of the four types of targets.Twenty particles are used in the calculation (M = 20), and we chose w = 0.5, c1 = 1.5, c2 = 2.5; when the relative error between the kth iteration and (k+1)th iteration is under 0.01, the calculation will stop.Figure 11b shows the sample centers of the training data.

H-Alpha Decomposition for GPR Data
After the collection for the full polarimetric GPR data of four types of targets, H-Alpha decomposition is applied to obtain the H and α data.Subsequently, the data of four targets are divided into two equal parts: the training data and the testing data.Figure 11a shows the spatial distribution of the training H and α data of four types of targets.

Particle Center Supported Plane
Subsequently, PSO is performed for the H and α data to obtain the sample centers of the four types of targets.Twenty particles are used in the calculation (M = 20), and we chose w = 0.5, c 1 = 1.5, c 2 = 2.5; when the relative error between the kth iteration and (k+1)th iteration is under 0.01, the calculation will stop.Figure 11b shows the sample centers of the training data.

H-Alpha Decomposition for GPR Data
After the collection for the full polarimetric GPR data of four types of targets, H-Alpha decomposition is applied to obtain the H and α data.Subsequently, the data of four targets are divided into two equal parts: the training data and the testing data.Figure 11a shows the spatial distribution of the training H and α data of four types of targets.

Particle Center Supported Plane
Subsequently, PSO is performed for the H and α data to obtain the sample centers of the four types of targets.Twenty particles are used in the calculation (M = 20), and we chose w = 0.5, c1 = 1.5, c2 = 2.5; when the relative error between the kth iteration and (k+1)th iteration is under 0.01, the calculation will stop.Figure 11b shows the sample centers of the training data.In Figure 12, the objective functions of every two sample centers are presented, the right part of the blue line is the range where Accuracy i > 0.8, and the left part of the red line is the range where Accuracy j > 0.8.Therefore, we must obtain the maximum value of fcost (shown in (18)) between the blue line and red line.The "Distance" coordinate corresponding to the maximum value of fcost indicates the best position of Q ij at the connection line of B gi and B gj .The obtained boundaries are shown in Figure 13 and the classification types of six boundaries are shown in Table 1.
In Figure 12, the objective functions of every two sample centers are presented, the right part of the blue line is the range where Accuracyi > 0.8, and the left part of the red line is the range where Accuracyj > 0.8.Therefore, we must obtain the maximum value of fcost (shown in (18)) between the blue line and red line.The "Distance" coordinate corresponding to the maximum value of fcost indicates the best position of Qij at the connection line of Bgi and Bgj.The obtained boundaries are shown in Figure 13 and the classification types of six boundaries are shown in Table 1.In Figure 12, the objective functions of every two sample centers are presented, the right part of the blue line is the range where Accuracyi > 0.8, and the left part of the red line is the range where Accuracyj > 0.8.Therefore, we must obtain the maximum value of fcost (shown in (18)) between the blue line and red line.The "Distance" coordinate corresponding to the maximum value of fcost indicates the best position of Qij at the connection line of Bgi and Bgj.The obtained boundaries are shown in Figure 13 and the classification types of six boundaries are shown in Table 1.L1 is the line between the multibranch and dihedral; L2 is the line between the multibranch and cylinder; L3 is the line between the multibranch and sphere; L4 is the line between the dihedral and cylinder; L5 is the line between the dihedral and sphere; L6 is the line between the cylinder and sphere.We used these six lines to classify the testing data; the classification results and accuracy are shown in Figure 14  L1 is the line between the multibranch and dihedral; L2 is the line between the multibranch and cylinder; L3 is the line between the multibranch and sphere; L4 is the line between the dihedral and cylinder; L5 is the line between the dihedral and sphere; L6 is the line between the cylinder and sphere.We used these six lines to classify the testing data; the classification results and accuracy are shown in Figure 14 and Table 2, respectively.The Figure 14 and Table 2 indicate that PCSP can classify these four types of targets with good accuracy.All the correct rates are greater than 80%.In fact, not all the lines contributed to the classification.After some analyses, the real plane partition is shown in Figure 15.The red lines in Figure 15a are the lines which contributed to the classification, the real classification plane is as Figure 15b.The Figure 14 and Table 2 indicate that PCSP can classify these four types of targets with good accuracy.All the correct rates are greater than 80%.In fact, not all the lines contributed to the classification.After some analyses, the real plane partition is shown in Figure 15.The red lines in Figure 15a are the lines which contributed to the classification, the real classification plane is as Figure 15b.

Identification for Real Data of Subsurface Pipes
According to Yue Yu [30], we can use the single polarimetric GPR to obtain the full polarimetric GPR data.Three time measurements with different measuring angles of antenna directions and survey direction are applied in this method as Figure 16 shows.We assume that the data collected by using this method are M0°, M45°, and M90°, respectively, the real matrix [S] can be represented as ( 22) [30].We applied this method to the data collection of subsurface pipes using 500 MHz shielded MALA GPR.The buried depths of two pipes were approximately 0.5 m and 0.55 m; the pipes extend two hundred meters in the north-south direction; the estimated pipe diameters were both approximately 20 cm; the applied GPR system had two horizontal polarimetric antennas whose interval is 18 cm.The number of measurement points was 98; the measurement point distance was 2 cm; the number of sampling points was 600.The applied GPR system and three types of measuring methods are shown in Figure 17.The full-polarimetric data obtained by using (22) are shown in Figure 18; the red rectangles indicate the positions of two subsurface pipes.The classification result of PCSP is shown in Figure 19.The accuracy is shown in Table 3.

Identification for Real Data of Subsurface Pipes
According to Yue Yu [30], we can use the single polarimetric GPR to obtain the full polarimetric GPR data.Three time measurements with different measuring angles of antenna directions and survey direction are applied in this method as Figure 16 shows.We assume that the data collected by using this method are M 0 • , M 45 • , and M 90 • , respectively, the real matrix [S] can be represented as ( 22) [30].

Identification for Real Data of Subsurface Pipes
According to Yue Yu [30], we can use the single polarimetric GPR to obtain the full polarimetric GPR data.Three time measurements with different measuring angles of antenna directions and survey direction are applied in this method as Figure 16 shows.We assume that the data collected by using this method are M0°, M45°, and M90°, respectively, the real matrix [S] can be represented as ( 22) [30].We applied this method to the data collection of subsurface pipes using 500 MHz shielded MALA GPR.The buried depths of two pipes were approximately 0.5 m and 0.55 m; the pipes extend two hundred meters in the north-south direction; the estimated pipe diameters were both approximately 20 cm; the applied GPR system had two horizontal polarimetric antennas whose interval is 18 cm.The number of measurement points was 98; the measurement point distance was 2 cm; the number of sampling points was 600.The applied GPR system and three types of measuring methods are shown in Figure 17.The full-polarimetric data obtained by using (22) are shown in Figure 18; the red rectangles indicate the positions of two subsurface pipes.The classification result of PCSP is shown in Figure 19.The accuracy is shown in Table 3.We applied this method to the data collection of subsurface pipes using 500 MHz shielded MALA GPR.The buried depths of two pipes were approximately 0.5 m and 0.55 m; the pipes extend two hundred meters in the north-south direction; the estimated pipe diameters were both approximately 20 cm; the applied GPR system had two horizontal polarimetric antennas whose interval is 18 cm.The number of measurement points was 98; the measurement point distance was 2 cm; the number of sampling points was 600.The applied GPR system and three types of measuring methods are shown in Figure 17.The full-polarimetric data obtained by using (22) are shown in Figure 18; the red rectangles indicate the positions of two subsurface pipes.The classification result of PCSP is shown in Figure 19.The accuracy is shown in Table 3.    Figure 19 shows that majority of the data points of three pipes focused on the blue zone.Table 3 shows that for the two pipes, the probabilities that they belonged to the type of cylinder were far higher than the other three types, the probabilities were 68.91% and 72.31%, respectively.Figure 19 shows that majority of the data points of three pipes focused on the blue zone.Table 3 shows that for the two pipes, the probabilities that they belonged to the type of cylinder were far higher than the other three types, the probabilities were 68.91% and 72.31%, respectively.Figure 19 shows that majority of the data points of three pipes focused on the blue zone.Table 3 shows that for the two pipes, the probabilities that they belonged to the type of cylinder were far higher than the other three types, the probabilities were 68.91% and 72.31%, respectively.Figure 19 shows that majority of the data points of three pipes focused on the blue zone.Table 3 shows that for the two pipes, the probabilities that they belonged to the type of cylinder were far higher than the other three types, the probabilities were 68.91% and 72.31%, respectively.

Identification for Three Pipes of Different Depths, Diameters and Materials.
To verify the feasibility of the proposed method.An additional experiment was performed.In this experiment, three pipes (Figure 20) of different depths, diameters, and materials were buried in the dry sand trough.The measurement environment and instruments were the same with that of Section 3.1.The parameters and buried states of three pipes are in Table 4 and Figure 21, respectively.The distance between the three pipes was 45 cm.The survey lines were right on the targets, the number of measurement points was 89, the measurement point interval was 2 cm, the frequency ranged between 800 MHz and 4000 MHz, and the number of sampling points was 1024.To verify the feasibility of the proposed method.An additional experiment was performed.In this experiment, three pipes (Figure 20) of different depths, diameters, and materials were buried in the dry sand trough.The measurement environment and instruments were the same with that of section 3.1.The parameters and buried states of three pipes are in Table 4 and Figure 21, respectively.The distance between the three pipes was 45 cm.The survey lines were right on the targets, the number of measurement points was 89, the measurement point interval was 2 cm, the frequency ranged between 800 MHz and 4000 MHz, and the number of sampling points was 1024.The measurement data are shown in Figure 22.Subsequently, we perform the PCSP to identify these three pipes.The identification results and accuracy are shown in Figure 23 and Table 5, respectively.To verify the feasibility of the proposed method.An additional experiment was performed.In this experiment, three pipes (Figure 20) of different depths, diameters, and materials were buried in the dry sand trough.The measurement environment and instruments were the same with that of section 3.1.The parameters and buried states of three pipes are in Table 4 and Figure 21, respectively.The distance between the three pipes was 45 cm.The survey lines were right on the targets, the number of measurement points was 89, the measurement point interval was 2 cm, the frequency ranged between 800 MHz and 4000 MHz, and the number of sampling points was 1024.The measurement data are shown in Figure 22.Subsequently, we perform the PCSP to identify these three pipes.The identification results and accuracy are shown in Figure 23 and Table 5, respectively.The measurement data are shown in Figure 22.Subsequently, we perform the PCSP to identify these three pipes.The identification results and accuracy are shown in Figure 23 and Table 5, respectively.The Figure 23 and Table 5 indicate that for the three pipes of different depths, diameters, and materials, most of the data points focused on the blue zone; the probabilities that they belonged to the type of cylinder are far higher than that of other three types, the probabilities were 88.69%, 80.56%, and 88.24%, respectively.The results prove that, for the pipes of different depths, diameters, and materials, the PCSP method can identify them with good accuracy, therefore the proposed method is appropriate.

Comparison with Classical H-Alpha Classification and Support Vector Machine (SVM)
In this section, classical H-Alpha classification and support vector machine (SVM) were performed to compare with the PCSP method.Classical H-Alpha classification is based on the classical template of SAR with nine zones (Figure 24), and the physical scattering characteristic associated with each zone provides information for classification [14].The Figure 23 and Table 5 indicate that for the three pipes of different depths, diameters, and materials, most of the data points focused on the blue zone; the probabilities that they belonged to the type of cylinder are far higher than that of other three types, the probabilities were 88.69%, 80.56%, and 88.24%, respectively.The results prove that, for the pipes of different depths, diameters, and materials, the PCSP method can identify them with good accuracy, therefore the proposed method is appropriate.

Comparison with Classical H-Alpha Classification and Support Vector Machine (SVM)
In this section, classical H-Alpha classification and support vector machine (SVM) were performed to compare with the PCSP method.Classical H-Alpha classification is based on the classical template of SAR with nine zones (Figure 24), and the physical scattering characteristic associated with each zone provides information for classification [14].The Figure 23 and Table 5 indicate that for the three pipes of different depths, diameters, and materials, most of the data points focused on the blue zone; the probabilities that they belonged to the type of cylinder are far higher than that of other three types, the probabilities were 88.69%, 80.56%, and 88.24%, respectively.The results prove that, for the pipes of different depths, diameters, and materials, the PCSP method can identify them with good accuracy, therefore the proposed method is appropriate.

Comparison with Classical H-Alpha Classification and Support Vector Machine (SVM)
In this section, classical H-Alpha classification and support vector machine (SVM) were performed to compare with the PCSP method.Classical H-Alpha classification is based on the classical template of SAR with nine zones (Figure 24), and the physical scattering characteristic associated with each zone provides information for classification [14].Support Vector Machine (SVM) is a new method of pattern recognition which was put forward by Corinna Cortes and Vapnik in 1995 [31].It has been applied to target classification for GPR [27].This method uses a hyper-plane to classify two different samples.The support vector is defined as the points in the two samples which are nearest to the hyper-plane, they are the most significant points to determine the hyper-plane and hard to be classified.The key idea is to minimize the spaces between support vectors and the plane to search for the optimal hyper-plane (Figure 25).The classical H-Alpha classification is an unsupervised method which does not need a training process; SVM is a supervised method which needs a training process.Therefore, the training data used in the SVM was the same as that of PCSP.The same testing data was applied in three methods.The sample points of testing data are projected into the classical H-Alpha classification space as Figure 26 presents.
From Figure 26, we can find that most sample points were in the zones they belong to, but the proportions were not high, particularly for cylinder and multibranch.Subsequently, the proportions of the sample points of the four targets in their zones were calculated; the SVM with linear kernel (L-SVM) was performed.The comparison of accuracy is shown in Table 6; the comparisons of operation time and RAM usage for L-SVM and PCSP are shown in Table 7. Support Vector Machine (SVM) is a new method of pattern recognition which was put forward by Corinna Cortes and Vapnik in 1995 [31].It has been applied to target classification for GPR [27].This method uses a hyper-plane to classify two different samples.The support vector is defined as the points in the two samples which are nearest to the hyper-plane, they are the most significant points to determine the hyper-plane and hard to be classified.The key idea is to minimize the spaces between support vectors and the plane to search for the optimal hyper-plane (Figure 25).Support Vector Machine (SVM) is a new method of pattern recognition which was put forward by Corinna Cortes and Vapnik in 1995 [31].It has been applied to target classification for GPR [27].This method uses a hyper-plane to classify two different samples.The support vector is defined as the points in the two samples which are nearest to the hyper-plane, they are the most significant points to determine the hyper-plane and hard to be classified.The key idea is to minimize the spaces between support vectors and the plane to search for the optimal hyper-plane (Figure 25).The classical H-Alpha classification is an unsupervised method which does not need a training process; SVM is a supervised method which needs a training process.Therefore, the training data used in the SVM was the same as that of PCSP.The same testing data was applied in three methods.The sample points of testing data are projected into the classical H-Alpha classification space as Figure 26 presents.
From Figure 26, we can find that most sample points were in the zones they belong to, but the proportions were not high, particularly for cylinder and multibranch.Subsequently, the proportions of the sample points of the four targets in their zones were calculated; the SVM with linear kernel (L-SVM) was performed.The comparison of accuracy is shown in Table 6; the comparisons of operation time and RAM usage for L-SVM and PCSP are shown in Table 7.The classical H-Alpha classification is an unsupervised method which does not need a training process; SVM is a supervised method which needs a training process.Therefore, the training data used in the SVM was the same as that of PCSP.The same testing data was applied in three methods.The sample points of testing data are projected into the classical H-Alpha classification space as Figure 26 presents.
From Figure 26, we can find that most sample points were in the zones they belong to, but the proportions were not high, particularly for cylinder and multibranch.Subsequently, the proportions of the sample points of the four targets in their zones were calculated; the SVM with linear kernel (L-SVM) was performed.The comparison of accuracy is shown in Table 6; the comparisons of operation time and RAM usage for L-SVM and PCSP are shown in Table 7.  Table 6 shows that the accuracy of L-SVM and PCSP were in the acceptable range, and the accuracy of PCSP was slightly higher than that of L-SVM.The correct rates were different with different targets.However, for the H-Alpha method, the accuracy of cylinder classification was under 50%, which means the cylinder can not be identified by using this method.
Table 7 shows that the RAM usage of PCSP is slightly higher than that of L-SVM.However, the operation time of PCSP was obviously less than that of L-SVM.In general, when RAM is sufficient, the accuracy of PCSP is slightly higher than that of L-SVM, and the calculation speed of PCSP is obviously faster than that of L-SVM.

Limitations of the Proposed Technique
However, the proposed technique for the underground object classification of GPR data are based on the classical parameters (H, α) obtained by H-Alpha decomposition; the type of objects we classified in this study was limited.Follow-up studies are underway to perform the classification for more types of objects using different polarimetric attributes and machine learning methods.Table 6 shows that the accuracy of L-SVM and PCSP were in the acceptable range, and the accuracy of PCSP was slightly higher than that of L-SVM.The correct rates were different with different targets.However, for the H-Alpha method, the accuracy of cylinder classification was under 50%, which means the cylinder can not be identified by using this method.

Conclusions
Table 7 shows that the RAM usage of PCSP is slightly higher than that of L-SVM.However, the operation time of PCSP was obviously less than that of L-SVM.In general, when RAM is sufficient, the accuracy of PCSP is slightly higher than that of L-SVM, and the calculation speed of PCSP is obviously faster than that of L-SVM.

Limitations of the Proposed Technique
However, the proposed technique for the underground object classification of GPR data are based on the classical parameters (H, α) obtained by H-Alpha decomposition; the type of objects we classified in this study was limited.Follow-up studies are underway to perform the classification for more types of objects using different polarimetric attributes and machine learning methods.

Conclusions
In this article, we proposed the particle center supported plane (PCSP) as a classification method based on sample centers to classify the H-Alpha data of ground penetrating radar.To verify the method, the PCSP method was performed to classify both data measured in the laboratory and real data measured on the road outdoors.The results indicate that this method could classify the four types of targets and subsurface pipes with good accuracy.To verify the feasibility of the proposed method, the identification for three pipes of different depths, diameters, and materials were performed, the results indicate that the PCSP method can identify the pipes of different depths, diameters, and materials with good accuracy.To verify the superiority of the proposed method, comparisons between the new method and classical H-Alpha classification, as well as L-SVM, were performed.The comparison of accuracy indicates that the accuracy of L-SVM and PCSP are in the acceptable range, and the accuracy of PCSP is slightly higher than that of L-SVM.However, for the H-Alpha method, the accuracy of cylinder classification was under 50%, which means the cylinder can not be identified by using the H-Alpha method.Moreover, the comparison of operation time and RAM usage indicate that if RAM usage is sufficient, the calculation speed of PCSP is obviously faster than that of L-SVM.

Figure 1 .
Figure 1.Outlier and sample center.(a) The influence of outlier, where Class 1 and Class 2 represent the samples of two different type of targets.(b) Sample center of four points.
H and α coordinates, respectively.The initial position x 0 m and initial velocity v 0 m are random.In this problem, the sample contains H and α, therefore, the training sample can be represented as follows:

Figure 1 .
Figure 1.Outlier and sample center.(a) The influence of outlier, where Class 1 and Class 2 represent the samples of two different type of targets.(b) Sample center of four points.

Figure 15 .
Figure 15.The real classification plane.(a) Lines that make contribution to the classification.(b) The real classification plane.

Figure 16 .
Figure 16.Three measuring methods with different measuring angles.

Figure 15 .
Figure 15.The real classification plane.(a) Lines that make contribution to the classification.(b) The real classification plane.

) 19 Figure 15 .
Figure 15.The real classification plane.(a) Lines that make contribution to the classification.(b) The real classification plane.

Figure 16 .
Figure 16.Three measuring methods with different measuring angles.

Figure 16 .
Figure 16.Three measuring methods with different measuring angles.

1 .
Identification for Three Pipes of Different Depths, Diameters and Materials.

Figure 20 .
Figure 20.Three different pipes in the try sand trough.

Table 4 .Figure 21 .
Figure 21.Three different pipes in the try sand trough.

Figure 20 .
Figure 20.Three different pipes in the try sand trough.

Figure 20 .
Figure 20.Three different pipes in the try sand trough.

Table 4 .Figure 21 .
Figure 21.Three different pipes in the try sand trough.

Figure 21 .
Figure 21.Three different pipes in the try sand trough.

Figure 22 .
Figure 22.The measurement data of three different pipes in the try sand trough.(a) HH polarization.(b) HV polarization.(c) VV polarization.

Figure 23 .
Figure 23.The identification results of three different pipes in the try sand trough.(a) Pipe 1.(b) Pipe 2. (c) Pipe 3.

Figure 22 .
Figure 22.The measurement data of three different pipes in the try sand trough.(a) HH polarization.(b) HV polarization.(c) VV polarization.

Figure 22 .
Figure 22.The measurement data of three different pipes in the try sand trough.(a) HH polarization.(b) HV polarization.(c) VV polarization.

Figure 23 .
Figure 23.The identification results of three different pipes in the try sand trough.(a) Pipe 1.(b) Pipe 2. (c) Pipe 3.

Figure 23 .
Figure 23.The identification results of three different pipes in the try sand trough.(a) Pipe 1.(b) Pipe 2. (c) Pipe 3.

Figure 25 .
Figure 25.Support vector machine, where Class 1 and Class 2 represent the samples of two different type of targets.

Figure 25 .
Figure 25.Support vector machine, where Class 1 and Class 2 represent the samples of two different type of targets.

Figure 25 .
Figure 25.Support vector machine, where Class 1 and Class 2 represent the samples of two different type of targets.

Table 1 .
The types of six boundaries.

Table 1 .
and Table2, respectively.The types of six boundaries.

Table 2 .
Accuracy of particle center supported plane.

Table 2 .
Accuracy of particle center supported plane.

Table 3 .
Classification probability of subsurface pipes.

Table 3 .
Classification probability of subsurface pipes.

Table 3 .
Classification probability of subsurface pipes.

Table 3 .
Classification probability of subsurface pipes.

Table 4 .
The parameters of three pipes.Identification for Three Pipes of Different Depths, Diameters and Materials.

Table 5 .
Identification accuracy of three pipes.

Table 5 .
Identification accuracy of three pipes.

Table 5 .
Identification accuracy of three pipes.

Table 6 .
Comparison of accuracy.

Table 7 .
Comparisons of operation time and RAM usage.

Table 6 .
Comparison of accuracy.

Table 7 .
Comparisons of operation time and RAM usage.