The asphalt pavement behavior under the FWD was modelled as an axisymmetric problem. The FWD loading was assumed to be uniformly distributed in a circular area with the radius of 0.15 m. The multilayered pavement consists of layers with a finite thickness above the infinite half-space. All layers are homogeneous, isotropic, and horizontally infinite. The asphalt layer is modelled as linear viscoelastic and other layers as linear elastic. Considering the linear viscoelastic property of AC and the dynamic effect of FWD loading, the spectral element method (SEM) was utilized as the forward analysis engine to simulate the structural behavior of asphalt pavement (Zhao et al., 2015). The governing equation for SEM in terms of displacement is shown in Equation 1. λ + μ ∇ ∇ · U + 2 μ ∇ 2 U = ρ U , ¨ (1) where λ = ν E / 1 + ν 1 − 2 ν and μ = E / 2 1 + ν are Lamé constants for materials related to the modulus E and the Poisson’s ratio ν ; U represents the displacement vector of the materials; ∇ is the differential operator; ∇ · U and ∇ 2 U are the divergence and the Laplacian of the displacement vector U , respectively; dots above the symbol represent the second differentiation of the vector with respect to time; and ρ is the density of materials.

In SEM, the linear viscoelastic property of AC was characterized by the modified MHN model (Zhao et al., 2013), as shown in Equation 2. Five parameters of the MHN in addition to the temperature-related factor could accurately simulate the linear viscoelastic behavior of AC. E * ω = E 0 * + E ∞ * − E 0 * / 1 + 1 + ω 0 i ω α β , (2) where E * ω is the complex modulus related to load frequency ω ; E 0 * and E ∞ * are complex moduli as ω approaches 0 and ∞ , respectively; ω 0 , α , and β are model coefficients. The parameters E ∞ * and E 0 * characterize the equilibrium and instantaneous behavior of the linear viscoelastic property. The coefficients α and β control the width and the skewness of the master curve, which are related to the relaxation mechanism behavior. The parameter ω 0 controls the horizontal position of the master curve along the frequency axis and comes into effect in the horizontal shifting between different test temperatures together with the time–temperature shift factor. It should be noted that the MHN model is essentially a complex modulus model, which means that the dynamic modulus and phase angle could be derived through the linear viscoelastic theory using the same coefficients. Only five coefficients of the asphalt layer are required to be determined, which can significantly enhance the accuracy of the layer property backcalculation compared to the Prony series, which requires more than 13 coefficients.

The governing equation was transferred into the frequency–wavenumber domain via Fourier transform. The dynamic effect and the inertia distribution of mass can be described accurately without using sublayers. One spectral element is used to simulate the whole layer of the pavement, and the total number of SEM elements is equal to the number of layers. The spectral element stiffness matrix for each layer is constructed, and the global stiffness matrix is assembled in each frequency–wavenumber domain following the same procedure as the finite-element method. The displacement is obtained based on the boundary conditions and the equilibrium equation. The solution is summed over the wavenumber domain to capture the frequency-dependent behavior of the system, and the inverse Fourier transform is used to obtain the solution in the time domain. The superiority of the SEM is owing to the fact that it integrates the analysis efficiency of the spectral analysis and the organization of the finite element method. The formulated spectral element can precisely characterize the wave propagation dynamics, and solving is carried out in the frequency domain without model discretization. The SEM significantly reduces the total number of degrees-of-freedom and enhances the computational efficiency. Considering the dynamic effect of the FWD impulse load and the viscoelastic properties of AC, the analysis model could objectively simulate the structural behavior under the FWD load. The detailed information and validation of SEM are not introduced herein and can be found elsewhere (Zhao et al., 2015; Cao et al., 2025; Cao et al., 2019).

Intelligent optimization has become a valuable tool to solve problems across various research fields, especially those related to engineering optimization. The ANN is the typical fundamental technique that enables simulating human behavior to discover the knowledge or the underlying patterns within the datasets. The unit of the ANN, known as a neuron, receives the information from various excitation sources and outputs the result through the active function in addition to a bias value, as shown in Equation 3. The primary tasks in constructing an ANN include designing the network architecture and determining the model’s hyperparameters (Li et al., 2022). The commonly utilized ANN structure consists of input layers, hidden layers, and output layers. The input layer is used to receive information from input variables, and the number of neurons corresponds to the number of input parameters. The output layer is the final layer of the ANN to return the prediction, and the number of neurons is related to the number of optimized parameters. The hidden layers connect the input layer and the output layer and perform the computation as a black box. The number of the hidden layer and the neurons in each layer varies and was related to the complexity of the problem. o u t p u t = f i n p u t + B i a s i s . (3)

The layer property backcalculation that identifies a reasonable combination of layer properties that accurately simulates the structural behavior under the FWD load can be formulated as a nonlinear optimization problem. The traditional method treated the difference between the measured and calculated peak deflections as the objective function, and the seed layer property was iteratively updated until the difference fell within a specified tolerance. Important characteristics related to the dynamic FWD load and damping materials, such as the time lag between the load peak and the peak deflection at various sensor locations, along with the deflection rate, are considered. However, only the peak value of the deflection time history is used in the static procedure, which limits the utilization of the full information contained in the measured data. Since the forward analysis engine must be executed numerous times, both the accuracy and computational efficiency are critical concerns for the iterative approach.

We used the ANN framework to build the mapping between the measured deflection and the layer property. In order to capture the influence of dynamic loading and damping materials, the measured time history of the FWD load and the deflection at various locations were considered as input variables in the ANN architecture. The length and complexity of the data posed a significant consideration for constructing the ANN input layer. The geometric characteristics were carefully extracted from the measured data to form the input variables. For the load time-history, four features were selected to represent the load properties, namely, peak load, peak time when the peak load occurs, load area, and load duration. For the deflection time history, the hysteresis curve was introduced to characterize the geometric property. Since the head and tail portions of the data are prone to environmental noise, the measured deflection greater than 30% of the magnitude and the corresponding loading were utilized to form the hysteresis curves, as shown in Figure 1. Six characteristics were extracted from each hysteresis curve. The peak deflection, the area under the hysteresis curve, and the ratio of the long axis to the short axis were considered the fundamental descriptors of each deflection curve. The peak deflection represents the deflection basin, while the area under the hysteresis curve and the ratio of the long axis to the short axis describe the relationship between deflection and loading. The enclosure degree was represented by the vertical distance between the initial and final deflection points, which were denoted as the length of the segment AB. The curve density was defined as the coordinate of the intersection point between the long and short axes, labelled as point O, which represented the center of the deflection peaks and the corresponding load. The inclination degree was defined as the tangent of the angle between the long axis and the tangent to segment AB. These characteristics were considered geometric indices to characterize the properties of the load–deflection hysteresis curve. The input variable characteristics are summarized in Table 1. It should be noted that the hysteresis curve was formed using the deflection and load time histories. The extracted geometric indices for deflections capture time domain characteristics, which allows the dynamic effects and damping influences to be reasonably represented. The number of neurons in the input layer was determined based on the features of the load time history and the load–deflection hysteresis curves. The output variable contains the mechanical property of each layer, including the linear viscoelastic coefficient of the asphalt layer and the elastic moduli for other layers.

This study employed a four-layer ANN model to establish the mapping between the measured deflection and layer properties (Wang and Zhao, 2022). The number of neurons in the input and the output layers was determined based on the structural information and the FWD configuration, while the number of neurons in the hidden layer was based on the relationships between adjacent layers. The weights and the biases of each neuron were considered as hyperparameters to be optimized. The weight connects the neurons and adjusts the importance of the contribution of each input excitation. The bias provides the flexibility of the neural network by acting as an offset or threshold. The ANN hyperparameters were initially predetermined and subsequently optimized during the training process (Marcelino et al., 2021; Mohammadi and Sheikholeslam, 2023). The built-in algorithm for the determination of hyperparameters of the ANN is the gradient-based method, which is prone to premature convergence and local optima (Abdolrasol et al., 2021). Various techniques have been applied to optimize the hyperparameters of the ANN. The widely used swarm intelligent optimization algorithms, namely, PSO, MA, and HPO, were utilized in this study to optimize the hyperparameters. Thus, a brief introduction of the algorithms mentioned above is provided in the following section.

The PSO is a well-known swarm-based optimization algorithm inspired by the social behavior observed in bird and fish flocks (Kennedy and Eberhart, 1995). The number of particles and the dimension of the PSO are related to the complexity of the optimization problem. The position and the velocity of the individual particles are important characteristics for the performance of the PSO algorithm. The position of each particle represents a candidate solution to the optimization objective, while the velocity is the critical parameter related to optimization efficiency. Initially, the position of each particle is uniformly distributed within the constrained search space and is updated by adding the current velocity vector. The position with the best fitness performance is recorded as the personal best ( p b e s t ) of the current particle and the global best ( g b e s t ) across all particles and generations. The exploration direction and velocity are updated based on the distance from the current position to the p b e s t and g b e s t positions, respectively, as shown in the following Equation 4. v i t + 1 = w v i t + c 1 r 1 p b e s t i − x i t + c 2 r 2 g b e s t − x i t (4) where, v i t and v i t + 1 are the velocity of the i ^th particle at the current iteration t and the next iteration t + 1 , respectively. p b e s t i and g b e s t are the optimal solutions of the current and the historical generations, respectively. w is the inertia factor governing the contribution of the current speed; c 1 and c 2 are the acceleration factors related to the contribution of the p b e s t and g b e s t , respectively. r 1 and r 2 are random numbers distributed between [0,1]. The kernel point of the PSO is that the global optimal solution is updated via shared information among the generations. The individual particle adjusts its exploration direction by updating the velocity using Equation 4 in each iteration until the objective error meets the requirement or the maximum iteration is reached.

Inspired by the social behavior of mayfly insects, the Mayfly algorithm was proposed for single-objective or multi-objective optimization problems (Zervoudakis and Tsafarakis, 2020). The position of each individual mayfly is considered the potential candidate for the optimization objective, while the speed is the optimization direction related to computational efficiency and precision, which is the same as that in PSO. The advantage of MA is that the mayfly in each generation is grouped into female and male sets, and the offspring are generated through crossover and mutation operations, similar to that in the GA. The MA combined the advantages of the PSO and the GA to enhance the optimization performance. The male MA will gather in swarms and dance up and down around the water level. The velocity of the position is updated following the Equation 5, considering the distance between from individual position to the p_best and g_best , respectively. v i j t + 1 = g v i j t + a 1 e − β r p 2 p b e s t i j − x i j t + a 2 e − β r g 2 g b e s t i j − x i j t , (5) where v i j t and v i j t + 1 are the velocities of the j ^th element in the i ^th male mayfly at iteration t and t+1, respectively; a 1 and a 2 are the positive attraction coefficients of social function; r p is the Cartesian distance between x i j t and p b e s t i j ; r g is the Cartesian distance between x i j t and g b e s t i j ; and g is the gravity parameter.

Both the female and male mayfliesare ranked based on the fitness error. The female mayfly flies toward the male mayfly with the same rank order following the attraction function. The position of the female mayfly is updated using the following Equation 6. v f i j t + 1 = g v f i j t + a 2 e − β r r m f 2 x i j t − x f i j t , i f f x f i > f x i g v f i j t + f l t × r , i f   f x f i < f x i , (6) where v f i j t and v f i j t + 1 are the velocities of the j ^th element of the i ^th female mayfly at time step t , and t + 1 , x i j t and x f i j t are the positions of the male and female mayflies with the same fitness rank. f l is a random walk coefficient, and r is a random value in the range [−1, 1]. Subsequently, the position of the female mayfly is updated by adding the velocity of the female at the current step, similar to the male mayfly.

After the optimized male and female mayflies are obtained at the end of the generation, the female will be paired with the male mayfly based on the attraction rules in MA. The crossover and the mutation operations are conducted to generate the offspring. The random selection strategy is utilized to extract the information from the parent generation. The inherited element from the parent generation is combined to form two offspring using the following Equation 7. o f f s p r i n g 1 = L × m a l e + 1 − L × f e m a l e o f f s p r i n g 2 = 1 − L × m a l e + L × f e m a l e , (7) where L is a random value within a certain range. It is noted that the offspring inherits both the advantage and disadvantage from the parents; a premature convergence will occur during the optimization process. A Gaussian distribution with the average of 0 and standard deviation of 1 is added into the offspring to lead the offspring to visit around the search space to obtain the global optimization, as shown in Equation 8. o f f s p r i n g n = o f f s p r i n g n + σ N n 0 , 1 , (8) where σ is the standard deviation of the normal distribution and N n 0 , 1 represents the normal Gaussian distribution. Thus, the generated offspring will be treated as initial members of the next generation. The optimization process will be performed following the routine until the accuracy meets the requirements.

The HPO method (Naruei et al., 2022) is another population-based method to simulate the dynamic interaction behavior of the animal hunting process. The core principle of the HPO is that the hunter will dynamically adjust the position based on the average position of the prey to capture the marginal prey, while the prey moves toward designated safety zones. The position search procedure is classified into two steps: exploration and exploitation. The exploration step contains a highly random behavior to discover the promising areas by updating the position significantly. The exploitation step contains a less random behavior to update the position around the promising regions. The position of the hunter will update following Equation 9. H i j t + 1 = H i j t + 0.5 * 2 C Z P p o s − H i j t + 2 1 − C Z μ j − H i j t , (9) where H i j t + 1 and H i j t are the positions of the hunter of the current and next iteration, respectively; P p o s is the location of the prey; μ j = 1 / N ∑ i = 1 N H i j is the average position of the hunter; C is the optimization parameter that controls the balance between exploration and exploitation; Z is the adaptive parameter. Once the prey is attacked, the prey will update its position to reach the safe zone, and the hunter may choose another prey. The global optimization is considered a safe place for the prey with a better chance of survival. The prey updates the position using the following Equation 10. P i , j t + 1 = T p o s j + C Z ⁡ cos 2 π R × T p o s j − P i , j t , (10) where P i , j t and P i , j t + 1 are the positions of the prey of the current and next iteration, respectively. T p o s j is the global optimization position. R is a random number in the range [−1, 1], and the COS function is to improve the prey exploitation performance by changing the radials and angles from the optimum position. According to the natural behavior, the trapped prey will die, and the hunter will move to the remaining prey swarm. A decreasing mechanism is adopted to simulate the variation in the swarm number following the Equation 11. k b e s t = r o u n d C × n , (11) where n is the number of search agents. The number of the search agent will decrease during the optimization process, in accordance with the nature that the hunter attacks the prey at the farthest distance from the average position. The optimization will stop when the termination criteria are met and k b e s t is equal to the requirements of the search agent.

The key point of ANN-based optimization methods is to establish the mapping between input and output variables, of which the preliminary task is to determine the hyperparameter. Since the hyperparameters have a crucial influence on the accuracy and robustness of the method, various swarm optimization algorithms have been applied to perform hyperparameter optimization. Due to the complexity and variability of engineering problems, a fixed framework cannot be appropriate for all field problems. We, thus, propose a universal backcalculation framework and evaluate the performance of intelligent optimizations at a project level. Thus, the procedure of the intelligent optimization combined with a swarm optimization algorithm should be a universal framework, and the implementation of the intelligent optimization in layer property backcalculation was briefly introduced.

Step 1

Generating the synthetic dataset using forward analysis

For the ANN-based method, the primary task is to establish the dataset for hyperparameter training and model validation. In this study, the analysis of an asphalt pavement under FWD loading is treated as a dynamic problem. In order to simulate the pavement behavior more accurately, the asphalt layer was treated as a linear viscoelastic material characterized by the MHN model, and others were treated as damped elastic materials. The SEM mentioned above was utilized to calculate the deflection time history. For a given pavement structure, the layer thickness and the physical properties of the paving materials were treated as known parameters. The mechanical properties of each layer, including the MHN model coefficients and elastic moduli for other layers, varied within a reasonable range. The deflection time history at locations with different offsets from the loading center was collected in the database.

Once the deflection time history was computed, the geometric characteristics of the applied load and deflection were extracted from the individual time-history and load–deflection hysteresis curves of each sensor. In total, four characteristics of the applied load time-history and six characteristics of the deflection–load hysteresis curve of each sensor were extracted to generate the synthesis database.

Step 2

Constructing the ANN structure and determining the hyperparameters of the ANN

The ANN shows excellent performance for the nonlinear mapping problem, where the complexity of the ANN plays an important role in the accuracy and efficiency. In this study, the geometric characteristics of the measured deflection and the applied load were treated as input variables, while the layer properties for each layer were output variables. Thus, the number of the input and output variables was determined based on the structural information and the field measurement configuration. Since there is no definitive method to determine the exact number of neurons in each hidden layer, this study sets the number of neurons per layer. Specifically, the number of neurons in the first hidden layer is equal to the size of the input layer, while the number of neurons in subsequent hidden layers is approximately half of that of the preceding hidden layer (Kim et al., 1995; Chou et al., 2017).

The unit of the ANN, denoted as a neuron, performs the function of the human brain to deal with the information source from the preceding layer, including the multiplication operation by a weight and the addition operation by a bias. The information is subsequently transferred into the next layer through the activation function. The activation function governs the transformation of the input data, while the weight and bias determine the influence of each input variable. We treated the weights and biases of each neuron as critical hyperparameters. The intrinsic optimization method in the ANN is typically gradient-based, which is susceptible to becoming trapped in the local optima. To address this, various population-based optimization algorithms were employed to initialize the weights and biases. The generated synthetic dataset was used to determine the initial hyperparameters. A piecewise chaos mapping function with the threshold of 0.4 (Varol Altay and Alatas, 2020) was used to pre-set the parameter with a uniform distribution within the given range following the Equation 12. f j + 1 = f j / 0.4 , 0 ≤ f j < 0.4 f j − 0.4 / 0.5 − 0.4 , 0.4 ≤ f j < 0.5 1 − 0.4 − f j / 0.5 − 0.4 , 0.5 ≤ f j < 1 − 0.4 , 1 − f j / 0.4 , 1 − 0.4 ≤ f j < 1 (12) where f j is a chaotic variable.

The pre-set weights and biases were assigned to the ANN to establish the relationship between the geometric characteristics of the measured data and layer properties. The normalized difference between the predicted and the actual layer property was considered the fitness function, as shown in the following equation. The weight and bias of each neuron were updated through the population-based optimization mentioned above until the termination condition was achieved. F = 1 n ∑ i = 1 n P a r a A c t u a l , i − P a r a B a c k , i / P a r a A c t u a l , i , (13) where F is the fitness function, n is the total number of predicted parameters, and P a r a A c t u a l , i and P a r a B a c k , i are the actual and backcalculated layer properties of analyzed pavement structures, respectively.

The swarm optimization algorithm shows excellent performance for global searching due to a highly random behavior, and the built-in gradient descent algorithm facilitates local optimization. This study used the swarm optimization to find the optimal ranges, leveraging its strong exploratory capabilities. The optimized results were then used as initial weights and biases for the ANN. The mature backcalculation framework was developed through the ANN training process, capitalizing on the advantages of local optimization and the potential for nonlinear mapping. The termination criterion was controlled by the number of iterations of 1,000, and the learning ratio was set as 0.015. Once the backcalculation framework was available, the layer property of various stations could directly be obtained based on the measured deflection time history.

Step 3

Optimizing the temperature-related parameter of the AC

The MHN coefficients of the asphalt layer and the elastic moduli of other layers could be obtained from the backcalculation model. Once the MHN model was available, the complex modulus, i.e., dynamic modulus and phase angle, can be computed based on the linear viscoelastic theory. The viscoelastic material exhibits temperature-dependent behavior; however, temperature-related parameters could not be determined because only one temperature data point was available from the single test station. The temperature correction for the backcalculated layer property should be performed for comparison and assessment at the project level. The AC is a typical linear viscoelastic material exhibiting temperature-related behavior. The WLF equation (Equation 14) was used to characterize the temperature-dependent property of AC. log ⁡ α T = − C 1 T − T 0 C 2 + T − T 0 , (14) where α T is the horizontal shift factor; C 1 and C 2 are the model coefficients; T 0 is the reference temperature; T is the field temperature. The air temperature and the surface temperature could automatically be collected while the FWD test was conducted. The temperature at the middle of the AC layer was treated as the field temperature and calculated based on the measured surface temperature and average air temperature using a regression function, as shown in Equation 15 (Asefzadeh et al., 2017). T z = 0.9511 × T s u r f + 0.1462 × Z sin 7.5442 × t + 6.6341 , (15) where T z is the temperature at depth Z , which is the middle depth of the AC layer in this study; T s u r f is the temperature at the surface and is collected from the FWD apparatus; t is the time of FWD test fraction of the day.

The WLF function coefficient was determined through an optimization process. The backcalculated layer property and the field temperature at a specified test station were selected as the reference conditions. The complex modulus of the asphalt layer at other stations could be computed using the WLF function via the linear viscoelastic theory, which was denoted as the corrected complex modulus. The difference between the corrected and the backcalculated layer property was adopted as the objective function, as shown in Equation 16. The optimization routine was performed using the HPO algorithm, which showed an excellent performance in the backcalculation of the asphalt layer property. R M S E = 1 M N ∑ i = 1 M ∑ j = 1 N E B a c k , i , j − E C o r r , i , j / E B a c k , i , j 2 + 1 M N ∑ i = 1 M ∑ j = 1 N ϕ B a c k , i , j − ϕ c o r r , i , j / ϕ B a c k , i , j 2 , (16) where R M S E is the root-mean-squared error between the backcalculated and transferred complex modulus at the i ^th station and j ^th frequency; the symbols E and ϕ represent the dynamic modulus and the phase angle, respectively. Thus, the mechanical behavior of the AC for a project could be comprehensively characterized using the backcalculated layer property and the coefficients of the WLF function.

A typical three-layer asphalt pavement was used in theoretical evaluation of the effectiveness of intelligent optimizations. The pavement information is listed in Table 2. The layer thickness was determined according to typical asphalt pavement structures. The Poisson’s ratio, damping ratio, and density were assigned according to the material type (Cao et al., 2025). During the forward analysis, the range for layer properties, including the MHN model coefficient of the AC layer and the elastic modulus of other layers, was determined according to the material type. The dynamic modulus range for viscoelastic materials was set between 5,000 MPa and 45,000 MPa, the elastic modulus for cement-treated aggregate ranged from 2,000 MPa to 25,000 MPa, and the elastic modulus for soil ranged from 30 MPa to 400 MPa. In order to improve the accuracy of the proposed method by extending the cover content, the layer property was carefully selected within the recommended range and randomly combined to create analysis scenarios. In total, 10,000 scenarios were generated to perform forward analysis.

Typical loading time history from various test equipment was collected from the long-term pavement performance (LTPP) program dataset, with the amplitude of 0.7 MPa. The surface deflection at the fixed offsets from the load center of 0, 210, 330, 500, 800, 1,100, 1,400, and 1,700 mm was computed to perform the layer property backcalculation. Four characteristics of the load and six characteristics of the deflection–load hysteresis curve of each offset were utilized as input variables, and the MHN model coefficient of the AC layer and the elastic moduli of other layers were used as output variables. A four-layer ANN (52-26-13-7) structure was established to build the mapping between input and output parameters. The weights and biases were optimized using the procedure mentioned above and the built-in algorithm.

A pure backpropagation neural network (BPNN) was used to validate the feasibility of the ANN in layer property backcalculation. The pure BPNN utilizes the gradient-based algorithm to optimize the hyperparameter. Figure 2 plots the comparison between the calculated and the backcalculated deflection time history, showing a good agreement between the measured and backcalculated deflection time history, thus indicating that the backcalculated layer property could be considered an appropriate solution for pavement behavior characterization. The backcalculated layer property and the corresponding actual values are listed in Table 3. Once MHN was available, the complex modulus could be derived based on the linear viscoelastic theory. The comparison between the actual and back-calculated dynamic modulus and phase angle master curves is plotted in Figure 3. The phase angle shows a discrepancy at lower frequencies, and a slight difference in the dynamic modulus could be observed at higher frequencies. Both the backcalculated dynamic modulus and phase angle master curves exhibit a good agreement during the effective frequency range (Fu et al., 2020). The normalized difference between the backcalculated elastic moduli of the base and subgrade was less than 10%. The negligible difference demonstrates the feasibility of the ANN-based layer backcalculation.

However, a crucial concern is the robustness of the ANN-based method, the performance of which depends on the optimization strategy. The results of the two analysis scenarios under duplicated trials are listed in Table 4. Both the backcalculated deflection histories matched well with the calculated data, which is not shown here for brevity. The dynamic modulus at the frequency of 10 Hz was calculated for comparison. The difference between the actual and the backcalculated layer modulus is listed in parentheses. The positive value denotes that the backcalculated modulus was larger than the actual one, and vice versa. A significant difference between the layer properties from different trials could be observed in the results. For scenario #1, the backcalculated surface layer property from trial 1 was larger than the actual one, while the backcalculated modulus of the base layer was lower than the actual one, which could be interpreted using the concept of layer compensation (Cao et al., 2019; Stubstad et al., 2006; Nakhaei and Timm, 2023). However, the backcalculated layer properties of both the surface and base layer from trial 2 were less than the actual ones. The normalized difference (ND) between the results obtained from the different trials was calculated and is listed in the last row. The ND for the surface layer between the two trials was 14.72%. For scenario #2, the difference between the backcalculated and actual base layer property could be more than 3,000 MPa, with an ND of approximately 20%. The result indicated that the robustness of the ANN-based method may produce an unreasonable result for pavement assessment.

In order to compare the performance of the backcalculation method with that of various intelligent optimization algorithms, the layer properties of scenarios #1 and #2 were backcalculated using the ANN model with different swarm optimization algorithms. The result was marked by the algorithms used in the hyperparameter optimization. The ND between the actual and the backcalculated parameter was also calculated and is listed in Table 5. The backcalculated layer properties from different models were reasonable, indicating the effectiveness of the backcalculation methods. The ND of the backcalculated layer property from the ANN model with swarm optimization algorithms was lower than that from the BPNN model. For the asphalt layers, the layer properties obtained from ANN + MA and ANN + HPO were closer to the actual values than those from BPNN and ANN + PSO. Especially for the base layer of scenario 2, the ND from the BPNN was 15.89%, while those from ANN + PSO, ANN + MA, and ANN + HPO were 12.04%, 4.34%, and 7.87%, respectively. The ND from the ANN + MA and ANN + HPO was almost 60% lower than that of the BPNN. The result indicated that the ANN + MA and ANN + HPO models exhibit a better overall performance for layer property backcalculation compared to BPNN and ANN + PSO.

In order to evaluate the performance of the proposed ANN-based model at the project level, a four-layer field asphalt pavement was selected to perform layer property backcalculation. The four-layer highway is located in Liaoning Province, China. The layer material and the thickness were determined by the as-designed document and are shown in Table 6. Three asphalt surface lifts were combined into one AC layer to reduce the complexity of the optimization. The field FWD test was performed with an interval of 50 m along the vehicle path. The magnitude of the test was set as 0.7 MPa, and the deflection time-history at various offsets of 0, 200, 300, 450, 600, 900, 1,200, 1,500, 1,800, and 2,100 mm was collected.

The mechanical property of each layer varied within the reasonable range according to the material type. A total of 30,000 scenarios were generated, and the deflection time history was calculated using the SEM. The characteristics of the loading and the deflection–load hysteresis curve were calculated to generate the database for training the ANN model with various swarm optimization algorithms. The input variable consisted of four characteristics of the load and six characteristics of the deflection–load hysteresis curve of 10 sensors. The output variable was the mechanical property of each layer. Based on the structural information and the test configuration, a four-layer ANN (64-32-16-8) architecture was developed to map the input to the output parameters. The ANN structure was trained using the theoretical deflection database following the procedure mentioned above. The trained ANN model was used to backcalculate the layer property based on field-measured deflection. The typical comparison between the measured and the backcalculated deflection from the BPNN model is shown in Figure 4, yielding a good agreement, indicating that the backcalculated layer property can accurately simulate the pavement behavior under FWD load conditions. The backcalculation results from other models showed a similar performance, which is not plotted here for brevity. The backcalculated layer properties from various models are listed in Table 7.

It can be seen that the layer property from various models fell within reasonable ranges and matched well with each other, indicating that the trained ANN model can effectively capture the essential relationship between the layer property and pavement performance. The complex modulus of the AC layer could be derived from the backcalculated MHN model and is shown in Figure 5. Within the effective frequency range of 5 Hz–65 Hz (Li and Wang, 2019), the master curve from all the models matched well with each other. Notably, the master curve was calculated using the backcalculated MHN model, for which the reference temperature was the field-measured data. The backcalculated layer property should be primarily corrected into the same temperature condition.

The field-measured deflection was carefully selected based on the criterion that no macroscopic crack appeared at the test station. The backcalculated layer property and the measured temperature of 27.81 °C of station 1 were treated as the reference scenario. The backcalculated layer properties of other stations were used to obtain the optimal coefficient of the WLF function within the effective range. The optimized WLF coefficients were 10.6 and 145.25 for C 1 and C 2 , respectively. Once the coefficient of the MHN model and the WLF function was available, the dynamic modulus at the artificial condition could be derived via the LVE theory. The dynamic modulus of the test stations could be obtained based on the layer property of station 1 and the WLF function coefficient. The comparison of the corrected dynamic modulus and the backcalculated dynamic modulus was plotted and is shown in Figure 6. A good agreement could be observed, indicating the feasibility of the temperature modification procedure.

In order to evaluate the performance of the proposed backcalculation method, the dynamic modulus at the frequency of 10 Hz and the temperature of 20 °C was calculated using the backcalculated LVE parameter and the measured temperature of each station. The box-and-whisker plot was adopted to illustrate the statistical characteristics of the backcalculated layer properties obtained from different methods, as shown in Figure 7. The length of the box-and-whisker plot represents the distribution performance of the result. The maximum, minimum, and average values of the backcalculated layer property were also inserted in the figure.

The backcalculated layer property from various methods varied within a reasonable range, indicating the effectiveness of the trained ANN model. For the asphalt layer, the difference between the maximum and minimum of the dynamic modulus from the BPNN model was greater than that from the ANN + PSO, ANN + MA, and ANN + HPO models. The interquartile range of the box-and-whisker from BPNN, ANN + PSO, and ANN + MA was greater than that from the ANN + HPO, indicating the ANN + HPO exhibits a better performance in layer property backcalculation for viscoelastic materials. The phase angle showed a similar trend to the dynamic modulus; however, the phase angle obtained from ANN + MA exhibited a more scattered distribution. For the elastic material, all backcalculation models can obtain moduli with a reasonable range in accordance with the corresponding material types. There were only slight differences among the average moduli obtained from the different methods. In order to quantitatively analyze the performance, the average, standard deviation (STD), and the coefficient of variation (COV) of the layer property from various models were calculated and are listed in Table 8.

As shown in Table 8, the average moduli for all the layers from various ANN-based models had a reasonable range in accordance with the material type. The difference of the average backcalculated layer property from different models was less than 10%. The STD from the BPNN model for all the layers was significantly greater than that from the ANN model with a swarm optimization algorithm. For the asphalt layer, the STD from the ANN + HPO model was 140.33, which was obviously less than 341.24 observed from the BPNN model. The COV of the backcalculated layer property from ANN + HPO was only 50% of that from the BPNN model. For the elastic layers, both the STD and COV from the BPNN and ANN + PSO were greater than those from ANN + MA and ANN + HPO. The STD and COV for the base layer from the ANN + MA were 38% and 45% less than those from the BPNN, respectively. The statistical characteristic of the ANN + MA was slightly less than that of the ANN + HPO model. Among the swarm optimization algorithm models, the ANN + HPO shows a better performance in the backcalculation of the asphalt layer, while the ANN + MA has a better performance for the elastic layers. The computational time for layer property backcalculation is a critical factor for engineering applications. The time consumption of the backcalculation procedure included forward analysis for generating the synthetic dataset and the hyperparameter optimization process. The backcalculation was performed on a computer with Intel i7 CPU of 3.20 GHz. The number of generations was set to 100, and the number of iterations was set to 150 for all procedures. The computational time for the SEM was approximately 3 s per scenario, while the optimization processes required 3,830, 3,780, 3,215, and 2,797 s for BPNN, ANN + PSO, ANN + MA, and ANN + HPO, respectively. The time consumption of the ANN + HPO model was significantly less than that of the other models.

The superiority of intelligent models lies in the ability of the swarm population to extend the candidate solution range and utilize information from each generation to approach the optimal result. The MA algorithm combines the advantages of PSO and the inheritance mechanisms of the GA; crossover and mutation operations leverage individual information to escape the local optima and find the global optimum. The HPO algorithm divides particles into hunters and prey, both of which update their positions toward the best solution. The interaction between these two groups can enhance the performance of the ANN model. The number of the search agents decreases during the optimization process, which could reduce the optimization amount and significantly enhance the computer efficiency. Integrating the ANN framework with swarm optimization algorithms enables the construction of a mapping between the measured deflection and layer properties, thereby improving pavement assessment at the project level. The influence of environmental noise, pavement deterioration on the accuracy of layer property backcalculation, and the application of the backcalculated layer property in the maintenance decision-making will be analyzed in subsequent research to improve the applicability of the proposed procedure.