A Hybrid Method Applied to Improve the Efficiency of Full-Waveform Inversion for Pavement Characterization

In the hybrid method, the MLPs are trained with the GPR data simulated by the GRP model proposed by Lambot et al. In their method, the GPR in the far-field region can be modeled by a linear equation system under the assumption that the electric field measured by the antenna tends to be a plane wave. The GPR signal can accordingly be calculated by [7]: where, b and a are the backscattered waves and the incident waves, respectively, $H_i$ is the mutual coupling coefficient between the transmitting and receiving antenna, $H$ is the transmitting-receiving transfer function accounting for antenna gain and phase delay, $H_f$ is the feedback loss accounting for the multiple reflections occurring between the antenna and pavement, and $G_{xx}^0$ is the zero-offset Green’s function of the pavement layers. The coefficients of the GPR model can be determined by solving the system of Equation (1) for different model configurations (denoted n, ranging from 1 to K). Same as that in [7], the well-defined model configurations with the antenna situated at different heights above a large copper plate is used. Therefore, the Green’s functions $G_{xx,n}^0(\omega )$ can be calculated, while the functions $S_{11,\mathrm{n}}(\omega )$ can be readily measured. Then, the system of Equation (1) can be rewritten as

The linear system of (2) can be written in matrix form as [7] where

Using the least squares approach, the vector $x$ can be easily calculated as where the superscript H denotes the Hermitian operation. More detailed information of the GPR modeling method is in [7].

When the GPR model is calibrated, it can be used to simulate the received GPR data under various pavement models to train the multilayer perceptrons (MLPs). MLPs network is one type of multilayer feedforward neural networks. It consists of an input layer, one or more hidden layers, and an output layer. MLPs have the ability to approximate any continuous function. The number of neurons in each hidden layer of MLPs is heuristically set according to the problem to be processed. The more neurons there are, the stronger approximation ability of the network, and the greater risk of overfitting [20]. However, the number of neurons in the output layer is determined by the dimensions of the problem being processed. For example, for a multi-class classification problem, the number of neurons in the output layer is equal to the number of categories to be classified. Figure 1 shows the MLPs with one hidden layer.

In Figure 1, $a_j^l$ denotes the activation of the $j\mathrm{th}$ neuron of the $l\mathrm{th}$ layer, $b_j^l$ denotes the bias of the $j\mathrm{th}$ neuron of the $l\mathrm{th}$ layer, $w_{jk}^l$ denotes the weight connecting the $j\mathrm{th}$ neuron of the $l\mathrm{th}$ layer and the $k\mathrm{th}$ neuron of the $\left(l-1\right)\mathrm{th}$ layer. With these notations, the activation value of the $j\mathrm{th}$ neuron of the $l\mathrm{th}$ layer can be associated with the activation values of the $\left(l-1\right)\mathrm{th}$ layer by the following equation where the summation is performed on all k neurons in the $\left(l-1\right)\mathrm{th}$ layer, and $\sigma$ denotes the activation function. In the classification scenario, the MLPs has to perform a mapping from the input space into a finite set of class $Y=\left\{1,\dots ,Lab\right\}$, where $Lab$ is the number of classes. In the training phase, the parameters of the MLPs are determined from a finite training set here each training data $V^{\gamma }$ is labelled with its class membership $y^{\gamma }$, and $\mathrm{TR}$ is the number of training data. To determine the parameters of the network, MLPs are trained in a supervised manner with a highly popular algorithm known as the error back-propagation algorithm [21]. This algorithm is based on the error-correction learning rule. Most often the mean square error (MSE) is minimized [22] where $y_j^{\gamma }$ means the desired value of the $j\mathrm{th}$ network output for the $\gamma \mathrm{th}$ instance, $a_j^{out,\gamma }$ denotes the actual value of the $j\mathrm{th}$ network output for the $\gamma \mathrm{th}$ instance. In the training process of MLPs, the weights and biases are updated in the manner of gradient descent where ${\widehat{w}}_{jk}^l$ and $w_{jk}^l$ are respectively the updated weight value and the weight value before the update in each iteration, ${\widehat{b}}_j^l$ and $b_j^l$ are respectively the updated bias and the bias before the update in each iteration, $\eta$ is the learning rate, and $\frac{\partial \cdot }{\partial \cdot }$ denotes the partial derivative operation. The training process of MLPs is sensitive to the learning rate. If the learning rate is too large, the error function will oscillate during the descent process. If the learning rate is too small, the error function will drop slowly which reducing the efficiency of training. A highly popular way is to make the learning rate exponentially decrease from a large value with the number of iterations of training. In order to improve the training efficiency of MLPs, some other techniques are often introduced in training process, such as dropout, regularization, bath normalization et al., and more detail information about MLPs can be found in [23,24]. In the inference phase, an instance without label is fed to the network, and the categorization is performed by assigning the input vector $V$ the class of the output neuron with maximum activation:

The proposed hybrid approach is a method based on supervised ML algorithm. In order for the ML algorithm to have a good effectiveness in practical applications, the training data need to cover all possible situations as much as possible. In pavement inspection, the pavement is often assumed to be horizontally stratified under the antenna footprint. Under this condition, the GPR signals emulated by the Lambot’s far-field model have a high degree agreement with the measured data [7]. Therefore, in the proposed hybrid method, the far field model is used to simulate the training data. In the simulation process, each simulated GPR signal is determined by the height of the antenna and the characteristic parameters of the pavement model. In pavement inspections, the antenna is always located at a height with a small variation due to the unevenness of the pavement surface. Hence, when simulating the training data, the height parameter of the antenna is set to a random value within a small variation range. In addition, two kinds of pavement models are used in the simulation (Figure 2). The first kind pavement model consists of two layers, the first layer of it being a thin layer, and the second layer being an infinite one. The second pavement model consists of three layers, the first layer of it being a thin layer, the second being a thick layer with respect to the probing pulse width, and the third one being an infinite layer. In order to estimate the thickness of the first layer and the permittivities of the first two layers, only the echoes from the first two interfaces need to be matched in the full-waveform inversion. Therefore, this three-layer pavement model can cover all situations where the number of layers is greater than two. Here, the two kinds of pavement models are used to generate training data set, which can ensure that the ML algorithm learns the essential features in these data and enhances the robustness of the prediction model.

According to the above descriptions, it is assumed that the thickness of the first layer of the pavement is within the range of $[{\mathrm{h}}_{1,\mathrm{Low}},{\mathrm{h}}_{1,\mathrm{Up}}]$, and the permittivity of the pavement is in the range of $\left[{\epsilon }_{\mathrm{Low}},{\epsilon }_{\mathrm{Up}}\right]$. In order to improve the efficiency of parameter search in the process of the full-waveform inversion, these parameter ranges are divided into small intervals. For example, the permittivity range is divided into $\left[{\epsilon }_{Low},{\epsilon }_1\right)$, $\left[{\epsilon }_1,{\epsilon }_2\right)$, …, $\left[{\epsilon }_{i-1},{\epsilon }_i\right)$, …, $\left[{\epsilon }_{\mathrm{N}-1},{\epsilon }_{\mathrm{Up}}\right]$ as N intervals, where $\Delta \epsilon ={\epsilon }_i-{\epsilon }_{i-1}$. Same as the previous operation, the thickness range of the first layer is divided into $\left[h_{1,Low},h_{1,1}\right)$, …, $\left[h_{1,i-1},{\mathrm{h}}_{1,i}\right)$, …, $\left[h_{1,\mathrm{M}-1},{\mathrm{h}}_{1,\mathrm{Up}}\right]$ as M intervals, where $\Delta h_1=h_{1,i}-h_{1,i-1}$. With these small parameter intervals, all possible properties combinations of the first two layers of the pavement can be divided into a limited number of categories, which corresponds to a multi-class classification problem. For each properties interval combination, a label $y$ can be assigned according to $y_{(l1,l2,l3)}=Label({\epsilon }_{1,l1},{\mathrm{h}}_{1,l2},{\epsilon }_{2,l3})$, where, l1, l2, l3 are respectively the index number of intervals in original ranges, $Label(\cdot )$ denotes a mapping operation. For the class $y_{(l1,l2,l3)}$, a number of random values are generated according to the uniform distribution over the parameter intervals belonging to it. Then, these randomly generated parameter combinations are used to simulate GPR signals, and these simulated GPR data belong to the class $y_{(l1,l2,l3)}$. It should be noted that for the three-layer pavement model, the thickness of the second layer can be set to a random value within a meaningful range, and the permittivity of the third layer can be simply set to a fixed value. However, the thickness of the second layer and the permittivity of the third layer do not affect the classification results. Since the GPR signal obtained by the simulation is in frequency domain form, it is rearranged as follow: where $\Re \left[S\left({\mathrm{f}}_k\right)\right]$ and $\Im \left[S\left({\mathrm{f}}_k\right)\right]$ are respectively the real part and imaginary part of the radar data at frequency ${\mathrm{f}}_k$, Q is the number of frequency points of the GPR data, and vector $V$ is the instance used to train the MLPs. Figure 3 shows the MLPs training process starting from the simulated data. During the phase of inference, an input GPR signal is predicted to a class by MLPs. Then, the midpoint values of all parameter intervals under this class are used as the initial solution for the full-waveform inversion of this GPR signal. If the classification result is accurate, the maximum distance between each component of the initial solution and its corresponding global optimal solution is $\delta /2$, where $\delta$ is the length of the interval. Therefore, if $\delta$ is small, the initial solution is close to the global optimal solution, so an efficient local optimization algorithm can meet the requirement of the full-waveform inversion, i.e., the Levenberg-Marquardt algorithm.

When the MLPs network is well trained, it can be used to classify the measured GPR data into a specific class $y_{(\mathrm{p}1,\mathrm{p}2,\mathrm{p}3)}$. Since the class $y_{(\mathrm{p}1,\mathrm{p}2,\mathrm{p}3)}$ is associated with the parameter intervals combination $\left\{\left[{{\displaystyle \epsilon }}_{1,p1-1},{{\displaystyle \epsilon }}_{1,p1}\right),\left[{{\displaystyle h}}_{1,p2-1},{{\displaystyle h}}_{1,p2}\right),\left[{{\displaystyle \epsilon }}_{2,p3-1},{{\displaystyle \epsilon }}_{2,p3}\right)\right\}$, the midpoint values of the these intervals under this combination are used as the initial solution for the full-waveform inversion. The full-waveform inversion is a nonlinear optimization problem which consists in finding the parameter vector $z=\left[{\epsilon }_u,{\sigma }_u,h_u\right]$, u = 1, …, R, so that an objective function $\mathrm{obj}(z)$ is minimized [7]. Since the GPR only receives signals in a limited time range, the full-waveform inversion is generally performed in the time domain. Therefore, the objective function is accordingly defined as follows: where $s^{*}(\mathrm{t})$ is the vector containing the time domain target signal for inversion, $s(\mathrm{t})$ is the vector containing the time domain signal simulated by the unknown parameter vector $z$ during the inversion process, $\mathrm{IFFT}[\cdot ]$ denotes the inverse Fourier transform, $S_{11}^{*}(\omega )$ denotes the measured GPR data, $S_{11}(\omega ,z)$ denotes the simulated GPR data, $\Vert \cdot \Vert$ denotes the norm of the vector, and the superscript T denotes the transpose operation. Since the inversion is only performed on the thickness of the first layer and the dielectric constants of the first two layers, it is necessary to properly select the time range t of the target signal so that only the echoes from the first two interfaces are included in the target signal. Within the selected time range, the target signal can be approximated as the GPR response of the two-layer pavement. Therefore, in the inversion process, the simulated GPR data $S_{11}(\omega ,z)$ is obtained by using the two-layer pavement model. In most electromagnetic inverse problems, this objective function is highly nonlinear and is characterized by an oscillatory behavior and a multitude of local minima [7]. This complex topography necessitates the use of a robust global optimization algorithm. However, the global optimization algorithms are often inefficient. In our method, the initial solution is provided by the MLPs instead of randomly selected. Since the initial solution provided by the MLPs is close to the global optimal solution, the inversion can be performed by using an efficient local optimization algorithm. Figure 4 shows the process of the inversion starting from the initial solution provided by the MLPs.

This section aims at illustrating that the proposed hybrid method can improve the efficiency of the full-waveform inversion. First of all, a mono-static GPR composed of a homemade horn antenna and a portable vector network analyzer (VNA) is modeled with the method introduced in Section 2.1. The portable VNA model is produced by COPPER MOUNTAIN TECHNOLOGIES (S5048) [25], and its operating frequency ranges from 20 kHz to 4.8 GHz. The dimensions of the homemade horn antenna are 22 cm length and 16 cm $\times$ 30 cm aperture area. Its nominal frequency is 0.75–4.5 GHz, and its isotropic gain ranges from 6–14 dBi. In the experiment, the mono-static GPR operates on the frequency range from 0.9 GHz to 3.5 GHz with 108 frequency points. To calibrate the mono-static GPR model, measurements were performed with antenna at 41 different heights from 0.4 m to 0.6 m over a 2.5 m $\times$ 2.5 m copper sheet. During the preparation of the training data, white noise is added to the simulated GPR data with a moderate signal to noise ratio (SNR) i.e., 20 dB. The antenna height is in the range from 0.45 m to 0.47 m. The thickness of the first layer is within the range from 2 cm to 5 cm, the relative permittivity of the first layer is within the range from 4 to 7, and the relative permittivity of the second layer is within the range from 5 to 8. The dielectric constant of the second layer is set larger than that of the first layer because the second layer has a larger degree of compaction, which is proportional to the dielectric constant [26]. According to the descriptions in Section 2.4, the thickness range of the first layer ($h_1\in [2,5]$) is divided into $\left[2,2.5\right)$, $\left[2.5,3\right)$, $\left[3,3.5\right)$, $\left[3.5,4\right)$, $\left[4,4.5\right)$, $\left[4.5,5\right]$; the permittivity range of the first layer (${\epsilon }_1\in \left[4,7\right]$) is divided into $\left[4,4.6\right)$, $\left[4.6,5.2\right)$, $\left[5.2,5.8\right)$, $\left[5.8,6.4\right)$, $\left[6.4,7\right]$; and the permittivity range of the second layer (${\epsilon }_2\in \left[5,8\right]$) is divided into $\left[5,5.6\right)$, $\left[5.6,6.2\right)$, $\left[6.2,6.8\right)$, $\left[6.8,7.4\right)$, $\left[7.4,8\right]$. In particular, for the three-layer pavement model, the thickness of the second layer is set to a random value between 8 cm and 12 cm, and the relative permittivity of the third layer is fixed at 4. It can be seen that the permittivity ranges of the first layer and the second layer are partially overlapped. In order to avoid insignificant parameter combinations (i.e., the permittivity of the first layer and second layer are equal), the following conditions must be met when combining these parameter intervals: where, l1, l2, l3 are respectively the index numbers of the intervals in original ranges, and $y_{(l1,l2,l3)}$ is the number of the class. There are a total of 90 classes. For each class, 150 GPR data are respectively simulated by using the two-layer pavement models and the three-layer pavement models. Therefore, there are a total of 27,000 training data. To validate the classification performance of the MLPs, fifteen percent of the training data is used as the test set, and the test set would not be used to train the MLPs. In the numerical experiment, MLPs with two hidden layers are trained. There are 200 and 150 neurons in the two hidden layers, respectively. The number of neurons in the input layer and the output layer is respectively 216 and 90. The learning rate is in the form of exponential decay, with a base of 0.9 and an attenuation factor of 0.99. The training process of the MLPs is finished in the framework of TensorFlow. In the experiment, the MLPs are trained for 400 iterations. After each iteration, the trained MLPs are applied to the test data, and the classification accuracy is plotted versus the number of iterations in Figure 5. It can be found that the maximum classification accuracy is 98.6%, which means that the well trained MLPs can accurately classify the test data into corresponding classes according to the preset rule with a significantly high probability.

To verify the performance of the hybrid method at the respect of full-waveform inversion, three GPR data are selected from the test set to perform full-waveform inversion, and two of which are correctly classified and one is misclassified. The two correctly classified GPR signals are respectively simulated by the two-layer pavement model and the three-layer pavement model. Since only the thickness of the first layer and the permittivities of the first two layers are estimated, the signals before 6.5 ns are matched when performing the full-waveform inversion. The correctly classified GPR data generated by the two-layer pavement model has the simulation parameter combination of ${\epsilon }_1=4.493$, $h_1=2.41$, ${\epsilon }_2=5.532$, and the corresponding class number of 1. According to the preset rule, the initial solution of the full-waveform inversion is ${\epsilon }_1^{\mathrm{int}}=4.300$, $h_1^{\mathrm{int}}=2.25$, ${\epsilon }_2^{\mathrm{int}}=5.300$. In addition, the antenna height is 0.462 m, which is considered as known in the full-waveform inversion. The antenna height is considered as known because it could be determined independently using other methods such as laser measurements in practical applications [8]. For the full-waveform inversion, the Levenberg-Marquardt (LM) algorithm is set with a maximum number of iterations of 70 and a step tolerance of $1\times {10}^{-5}$. The inversion results are ${\epsilon }_1^{\mathrm{inv}}=4.494$, $h_1^{inv}=2.41$, ${\epsilon }_2^{inv}=5.535$, which agree well with the simulation parameter combination. Figure 6a shows the GPR data simulated by the inversion results and the target signal, and significant agreement can be seen between them. Figure 6b shows the trajectory of the objective function during the inversion process. It can be seen that the algorithm converges after only thirty iterations. According to the previous study on the topology structure of the inversion objective function [8], the function gradient is small in some ranges of the inversion parameters and large in other ranges. Therefore, in the inversion process, the objective function will show a tendency that it sometimes changes slowly and sometimes changes dramatically (Figure 6b).

The correctly classified GPR data generated by the three-layer pavement model has the simulation parameter combination of ${\epsilon }_1=5.115$, $h_1=4.42$, ${\epsilon }_2=7.992$, $h_2=11$, ${\epsilon }_3=4$, and the corresponding class number of 53. Since only the thickness of the first layer and the permittivities of the first two layers are estimated, the time range for the inversion is within 6.5 ns, and the initial solution of the full-waveform inversion is ${\epsilon }_1^{\mathrm{int}}=4.900$, $h_1^{\mathrm{int}}=4.25$, ${\epsilon }_2^{\mathrm{int}}=7.700$. The antenna height is 0.460 m, and the configuration of the LM algorithm remains the same as before. The inversion results are ${\epsilon }_1^{\mathrm{inv}}=5.112$, $h_1^{inv}=4.42$, ${\epsilon }_2^{inv}=8.011$, which agree well with the simulation parameters of the target signal. Figure 7a shows the GPR signal simulated by the inversion results and the target signal, and significant agreement can be seen between them within 6.5 ns. Figure 7b shows the trajectory of the objective function during the inversion process. It can be found that the algorithm also converges after about thirty iterations. For the misclassified GPR data, it has the simulation parameter combination of ${\epsilon }_1=5.680$, $h_1=3.92$, ${\epsilon }_2=6.850$, $h_2=10$, ${\epsilon }_3=4$, and the corresponding class number of 64. However, it is misclassified into class 58. According to the preset rule, the initial solution derived from the misclassification result is ${\epsilon }_1^{\mathrm{int}}=5.500$, $h_1^{\mathrm{int}}=3.75$, ${\epsilon }_2^{\mathrm{int}}=6.500$. Then, the initial solution is used for the full-waveform inversion, and the trajectory of the objective function during the inversion process is plotted in Figure 8. If this GPR data is correctly classified, its initial solution would be ${\epsilon }_1^{\mathrm{int}}=5.500$, $h_1^{\mathrm{int}}=3.75$, ${\epsilon }_2^{\mathrm{int}}=7.100$ according to the preset rule. This initial solution deduced by the correct classification result is also used for full-waveform inversion, and the trajectory of the objective function during the inversion process is also plotted in Figure 8. It can be seen from Figure 8 that the LM algorithm converges faster when the initial solution derived by the correct classification result is used. Nevertheless, the LM algorithm converges fast for both of the two initial solutions. The inversion results for both the two initial solutions are ${\epsilon }_1^{\mathrm{inv}}=5.682$, $h_1^{inv}=3.92$, ${\epsilon }_2^{inv}=6.859$, which agree well with the simulation parameters of the target signal. In order to understand the phenomenon in Figure 8, the 2-D response surface topography of the objective function is computed. The objective function is computed with $h_1$ and ${\epsilon }_2$ as input. The value of ${\epsilon }_1$ and the antenna height are considered as known. The misclassified test GPR data is the target signal of inversion. To better highlight the topography, the objective function values are express in a logarithmic scale (Figure 9). Figure 9 shows that several local minima are distributed around the global minimum, which is why the full-waveform inversion of the GPR data in the conventional methods require global optimization algorithms. The global minimum is in the center of a basin with a radius of 0.5 along the ${\epsilon }_2$ direction and a radius of 1 cm along the $h_1$ direction. These radii are greater than half the lengths of the corresponding parameter intervals. This fact indicates that if a GPR data is correctly classified by the MLPs, the corresponding initial solution falls into the basin centered on the global optimal solution. Therefore, the gradient-based local optimization algorithm LM can quickly search from the initial solution to the global optimal solution. Since the MLPs extract a set of features from the training GPR data for each class, a newly input GPR signal is predicted to a class based on the degree of the similarity to these features. Because the features of the GPR data are determined by the characteristic parameters of the pavement, when a parameter of the GPR data is close to the parameter interval boundary, it may be misclassified to the nearest neighbor interval with a small probability. In this case, the distance between the initial solution and the global optimal solution would be slightly larger than the length of the half interval. Nevertheless, since the length of the parameter interval is relatively small, the initial solution of the misclassified GPR data will still falls within the basin centered on the global optimal solution. For example, the global optimal solution of the misclassified test GPR data is $opt=\left\{h_1=3.92,{\epsilon }_2=6.850\right\}$, and the initial solution derived from the misclassification result is $\mathrm{mis}\_\mathrm{cls}=\left\{h_1^{\mathrm{int}}=3.75,{\epsilon }_2^{\mathrm{int}}=6.500\right\}$. It can be seen from Figure 9 that the initial solution derived from the misclassification result still falls in the basin centered on the global optimal solution. If the above test data is correctly classified according to the experimental settings, the initial solution would be $cor\_cls=\left\{h_1^{\mathrm{int}}=3.75,{\epsilon }_2^{\mathrm{int}}=7.100\right\}$. The distance between the initial solution derived from the misclassification result and the global optimal solution in the ${\epsilon }_2$ direction is 0.35, and the distance between the initial solution derived from the correct classification result and the global optimal solution in the ${\epsilon }_2$ direction is 0.25. Consequently, the misclassification result increases the distance between the initial solution and the optimal solution by 40%, and also increases the number of iterations of the inversion process.

The full-waveform inversion of these three test data have also been performed by conventional DE algorithm, which is the state-of-the-art method addressing the “thin layer problem”. The parameter search is performed on the corresponding entire range (${\epsilon }_1\in \left[4,7\right]$, $h_1\in [2,5]$, ${\epsilon }_2\in \left[5,8\right]$), the number of the population is 10, and the maximum generation is 80 [10]. Table 1 shows the number of calculations of the objective function to find the global optimal solution. It can be seen from Table 1 that since the initial solution provided by the MLPs is close to the optimal solution, the optimal solution can be found after a small number of searches. The DE algorithm encounters many local minima in the inversion process, so a lot of trials are needed to find the global optimal value. Therefore, the hybrid method can improve the efficiency of the full waveform inversion.

A field experiment is conducted with the GPR system at TIAN BEI road in Beijing (China). TIAN BEI road is an asphalt pavement which consists of three layers. The first layer of the pavement is the wear layer, its design thickness is 4 cm, and its aggregate size is less than 10 mm (AC-10). Due to the small size of the aggregate, the first layer can be regarded as a uniform layer, which is much in line with the pavement model assumption of the Lambot’s method. The design thicknesses of the second and the third layers are 10 cm and 12 cm, respectively. The second layer and the third layer have a similar aggregate ratio, which are composed with aggregate AC-16 and aggregate AC-23. Therefore, the echo from the interface between the second layer and the third layer is weak. Figure 10a shows the experimental setup. We use the GPR to perform A-scan above the road, and the antenna height is 0.46 cm. After the A scan, a field core of the pavement was drilled. Figure 10b is a core of the pavement. The red line in Figure 10b denotes the interface between the first layer and the second layer, and the measured thickness of the first layer was 29 mm. The waveform inversion of the corresponding GPR was performed by the means of DE algorithm and the hybrid method, respectively. Since the thickness of the first layer and the permittivities of the first two layers are estimated, only the signal before 6.5 ns is matched. The inversion results for both the methods are ${\epsilon }_1=6.06$, ${\epsilon }_2=6.89$, $h_1=30.2\mathrm{mm}$, and the thickness estimation error is 1.2 mm. Figure 11 shows the measured signal and the signal simulated by the inversion results. Significant agreement between the measured signal and the simulated signal can be seen before 6.5 ns. Signals from 7 ns to 8 ns primarily consist of the multiple reflections from the pavement surface and the echo from the interface between the second and the third layer. Since the echo from the interface between the second layer and the third layer is weak, the simulated signal and the measured signal agree well between 7 ns to 8 ns. The difference between the measured signal and the simulated signal at 9 ns is due to the reflection between the asphalt pavement and the lime-ash soil. Figure 12a shows the number of calculations of the objective function during the inversion process when using the DE algorithm. For the DE algorithm, the number of population is 10, so objective function would be calculated by 10 times in each generation. Figure 12b shows the number of calculations of the objective function of the hybrid method. It can be seen that the objective function is respectively calculated by 500 times and 70 times for the DE algorithm and the hybrid method to find the global optimal solution. Therefore, the hybrid method greatly improve the efficiency of the full-waveform inversion.

As shown by the results in the Section 3.1, the performance of the proposed hybrid method depends on whether the MLPs can accurately classify the input GPR data. When the classification result of the MLPs is accurate, the initial solution derived from the classification result will be in the basin centered on the global optimal solution. In this case, the global optimal solution can be quickly searched by using the gradient-based local optimization algorithm. If GPR signal has a corresponding layer thickness and permittivity that exceeds the ranges of parameters used to train the MLPs, the MLPs cannot accurately classify the GPR signal. Therefore, the performance of the hybrid method will be affected. In the numerical experiment, the ranges of the layer thickness and permittivities used to train the MLPs are respectively ${\epsilon }_1\in \left[4,7\right]$, $h_1\in [2\mathrm{cm},5\mathrm{cm}]$, ${\epsilon }_2\in \left[5,8\right]$. We use two parameter combinations $S1=\left\{{\epsilon }_1=3.40,h_1=2.30,{\epsilon }_2=5.30\right\}$, $S2=\left\{{\epsilon }_1=4.40,h_1=1.40,{\epsilon }_2=5.70\right\}$ to simulate GPR data. For $S1$, the value of ${\epsilon }_1$ exceeds 15% of the specified range, For $S2$, the value of $h_1$ exceeds 30% of the specified range. The GPR signals simulated by $S1$ and $S2$ are sent to the hybrid method. Table 2 shows the inversion results of these two GPR signals.

It can be seen from Table 2 that since the parameters of $S1$ and $S2$ are beyond the specified parameter ranges to a large extent (more than 15%), the initial solution derived from the classification result of MLPs is far from the global optimal solution, so the estimation errors of the hybrid method are great. In order to accurately invert the parameters of $S1$ and $S2$, we adjust the training parameter ranges of the MLPs to ${\epsilon }_1\in \left[3,6\right]$, $h_1\in [1\mathrm{cm},4\mathrm{cm}]$, ${\epsilon }_2\in \left[4,7\right]$. These parameter ranges are divided into small intervals according to the same operations in Section 3.1. The interval lengths of ${\epsilon }_1$ and ${\epsilon }_2$ are 0.6, and the interval length of $h_1$ is 0.5 cm. The MLPs uses the 27,000 GPR data generated over these parameter intervals for training, with a highest classification of 98.7%. The GPR signals simulated by $S1$ and $S2$ are again sent to the hybrid method. Table 3 shows the inversion results of $S1$ and $S2$ after the training parameter ranges adjustment. It can be seen from the results in Table 3 that when the ranges of the training parameters of the MLPs is adjusted to cover the parameters of $S1$ and $S2$, the parameters of $S1$ and $S2$ can be accurately inverted using the hybrid method. Combining the results of Table 2 and Table 3, it can be seen that the minimum layer thickness and permittivity that can be successfully detected by the hybrid method depend on the ranges of training parameters of the MLPs. In pavement inspection, the training parameter ranges of MLPs can be set to cover the possible properties values, so the hybrid method will have good performance.

Just like the hybrid method, the state-of-the-art full-waveform inversion method based on DE algorithm is also affected by the ranges of search parameters set. We compared the performance of the two methods in case where the characteristics corresponding to the GPR data are slightly outside the training parameter ranges of the hybrid method and the parameter search ranges of the DE inversion method, for example, within 5%. We used other two parameter combinations $S3=\left\{{\epsilon }_1=3.80,h_1=2.30,{\epsilon }_2=5.20\right\}$ and $S4=\left\{{\epsilon }_1=4.40,h_1=1.90,{\epsilon }_2=5.70\right\}$ to simulate GPR data, and the training parameter ranges of the hybrid method and parameter search ranges of the DE method are same (${\epsilon }_1\in \left[4,7\right]$, $h_1\in [2\mathrm{cm},5\mathrm{cm}]$, ${\epsilon }_2\in \left[5,8\right]$). For $S3$ and $S4$, the values of ${\epsilon }_1$ and $h_1$ exceed 5% of the specified ranges. Table 4 shows the inversion results of the hybrid method and the DE inversion method for $S3$ and $S4$ in case where the parameters values are slightly outside the specified ranges. The results show that the hybrid method can accurately invert these parameters, and the DE inversion method has a maximum thickness estimation error of 0.12 cm and permittivity estimation error of 0.43.