Prediction of Frequency for Simulation of Asphalt Mix Fatigue Tests Using MARS and ANN

Fatigue life of asphalt mixes in laboratory tests is commonly determined by applying a sinusoidal or haversine waveform with specific frequency. The pavement structure and loading conditions affect the shape and the frequency of tensile response pulses at the bottom of asphalt layer. This paper introduces two methods for predicting the loading frequency in laboratory asphalt fatigue tests for better simulation of field conditions. Five thousand (5000) four-layered pavement sections were analyzed and stress and strain response pulses in both longitudinal and transverse directions was determined. After fitting the haversine function to the response pulses by the concept of equal-energy pulse, the effective length of the response pulses were determined. Two methods including Multivariate Adaptive Regression Splines (MARS) and Artificial Neural Network (ANN) methods were then employed to predict the effective length (i.e., frequency) of tensile stress and strain pulses in longitudinal and transverse directions based on haversine waveform. It is indicated that, under controlled stress and strain modes, both methods (MARS and ANN) are capable of predicting the frequency of loading in HMA fatigue tests with very good accuracy. The accuracy of ANN method is, however, more than MARS method. It is furthermore shown that the results of the present study can be generalized to sinusoidal waveform by a simple equation.


Introduction
Fatigue cracking, due to repeated traffic load, is among the most common mode of failures that flexible pavements experience. The fatigue life of asphalt mixes in laboratory tests is commonly determined by applying a dynamic load having sinusoidal or haversine waveform with a specific frequency. The frequency of horizontal tensile stress and strain pulses in both longitudinal and transverse directions depends on several factors such as vehicle speed, loading properties, environmental conditions, and pavement structure.
Several relationships were proposed to predict the duration of vertical stress pulse at different depth of asphalt layers [1][2][3][4][5][6]. Nevertheless, the knowledge on the frequency of stress and strain pulses in horizontal directions at the bottom of the asphalt layer is still poorly understood. Brown (1973) derived an equation to calculate the loading time as a function of vehicle speed and depth beneath the pavement surface. The loading time was considered as the average of the pulse times of the stresses in the three directions as obtained from the elastic layered theory. The relationship between the loading time (s), depth (m), and vehicle speed V (km/h) was as follows [7]: log ( ) = 0.5 + 0.2 − 0.94 log ( ) . (1) When (1) is plotted for different speeds and thicknesses between 150 and 400 mm, it can be seen that the approximation = 1/ ( = average speed in km/h) is a reasonable fit for the range of thicknesses studied [8].
In the development of the Mathematical Model of Pavement Performance (MMOPP), Ullidtz (2005) used the loading time corresponding to the middle depth of the asphalt layer. It was calculated based on the simplified assumption that the load at that depth is uniformly distributed over a circular area with the radius of + ℎ 1 /2, where is the radius of the contact area and ℎ 1 is the thickness of the asphalt layer. The Scientific World Journal Based on this assumption, the time of loading is defined as follows [5]: where is the time of loading, a is the radius of the contact area, ℎ 1 is the thickness of the asphalt layer, and is the vehicle speed. According to Ullidtz (2005), since no reductions are made for the influence of dual tires or for lateral distribution of the loads, the results should be on the conservative side [5]. Garcia and Thompson (2008) measured the durations of longitudinal and transverse tensile strain pulses in four sections tested with the accelerated pavement testing machine (ATLAS) [9]. They found a very strong relationship between the longitudinal and transverse strain pulse durations. In general, the transverse pulse durations were significantly, about three times, higher than those in the longitudinal direction [9]. Robbins and Timm (2009) used the instrumented sections of the National Center of Asphalt Technology (NCAT) to validate the procedure proposed by the Mechanistic Empirical Pavement Design Guide (MEPDG) to calculate the longitudinal pulse duration [10]. They developed a regression model for the duration of longitudinal strain pulse at the bottom of asphalt layer with three variables as follows [10]: where is the strain pulse duration; ℎ is the thickness of asphalt layer; is the vehicle speed, and is the mid-depth temperature with , , , and as regression coefficients. Hernandez (2010) performed an experimental testing program at the Accelerated Pavement Load Facility (APLF) of Ohio University on four pavement test sections [11]. He studied the influence of load, temperature, offset, and thickness of asphalt layer on the amplitude and duration of strain pulses at the bottom of the asphalt layer in both longitudinal and transverse directions. These experimental results showed that the load amplitude does not affect the longitudinal pulse duration [11]. The influence of the offset was furthermore dependent on the magnitude of the applied load. According to Hernandez (2010), the procedures recommended by MEPDG and Hu et al. (2010) [6] are not suitable to predict the magnitude of the longitudinal pulse duration, overestimating the desired variable by more than 2.5 times [11].
Restrepo-Velez (2011) evaluated the effect of several factors on the duration of tensile strains at the bottom of asphalt layer [12]. The pavement responses were measured on the perpetual section AC 664, of the WAY-30 project. Similar to the observation by Hernandez [11], Restrepo-Velez (2011) also concluded that the pulse durations were higher at lower values of speed and temperature. Additionally, Restrepo-Velez (2011) compared the observed responses with pavement responses predicted using the MEPDG method and the multilayer elastic analysis software, JULEA. Based on this comparison, MEPDG procedure led to an overprediction of the strain pulse durations of around 80% compared to those measured in the field [12]. It is worth nothing that, some studies have shown that the loading times calculated from the longitudinal strain gauges fall within the method suggested by MEPDG method [9,13]. Based on the results of these studies there is no agreement in applicability of MEPDG method for prediction of loading times in longitudinal direction. This is due to the fact that the MEPDG method was originally developed to compute the loading time resulted from vertical stress, not horizontal responses.
According to the viscoelastic analysis of 112 flexible pavement sections,  proposed regression equations for determining the duration of tensile stress and strain pulses at the bottom of asphalt layer in both longitudinal and transverse directions [14]. Proposed equations were developed based on haversine waveform by means of weighted nonlinear regression. Using this method of regression, the fitted haversine waveform only depends on three parameters including the speed of moving wheel, the thickness of asphalt layer, and asphalt layer temperature [14].
A recent study showed that several parameters such as thickness of different layers, ratio of resilient modulus for two succeeding layers, and contact radius of tire affect shape and frequency of longitudinal and transverse response pulses at the bottom of asphalt layer [15]. All of these parameters should be therefore considered to build up a comprehensive method for predicting loading time (i.e., frequency) of longitudinal and transverse response pulses. To the best of our knowledge, no general relationship or method has been proposed for determining of frequency of tensile horizontal stress and strain pulses at the bottom of asphalt layers. Existing relations are basically only applicable to vertical stress durations and limited number of equations which were developed based on field measurement of horizontal strain pulses cannot be used generally in other situations, where the pavement section and loading characteristics are very different from the desired sections.
The objective of this paper is to propose two methods for predicting the frequency of tensile stress and strain pulses at the bottom of asphalt layer in both longitudinal and transverse directions based on haversine and sinusoidal waveforms. By applying these methods, the frequency of loading in HMA fatigue laboratory tests such as four-points bending beam and indirect tensile (IDT) fatigue tests can be determined more realistically based on pavement design and loading characteristics. Results of this study can be used for more realistic simulation of dynamic loading of asphalt mix fatigue tests in both stress control mode and strain control mode.

Analysis of Pavement Sections.
In order to establish the dataset related to the normalized stress and strain pulses at the bottom of the asphalt layer, 5000 flexible pavement sections were analyzed using layered elastic theory (LET). The longitudinal and transverse stress and strain values were calculated at different radial distances from the center of the contact area. Moving load was assumed to be a single wheel The Scientific World Journal 3 having the contact pressure of 700 kPa. Since the contact pressure has no effect on the normalized shape and duration of response pulses, only the contact radius was considered as a variable. In each case, the pavement structure was considered as a four-layered system and all layers were treated as linear elastic. Interface of two succeeding layers was considered as fully-bounded. Pavement structure characteristics as well as tensile responses considered in this study are illustrated in Figure 1. The range of thickness of layers, the ratio of resilient modulus of each layer to the immediate succeeding layer, and the radius of the contact area are given in Table 1. Each pavement section was analyzed using layered elastic analysis program, Non-PAS, which enables the analysis of a pavement system consisting of a maximum of ten linear or nonlinear elastic layers subjected to the maximum of ten circular contact loads. Detailed verification of Non-PAS program using Kenlayer program confirms that the Non-PAS program can accurately predict the pavement responses subjected to single and multiple loading [16]. Because of the limitation of the available pavement analysis programs such as Kenlayer and ELSYM 5 in respect to the number of response points, developing such a program was necessary. The developed code facilitates the calculation of responses (stresses, strains, and deflections) for an unlimited number of points in radial direction.
Stress and strain values in the longitudinal and transverse directions were computed at different radial distances with one centimeter interval, as long as the amount of tensile stress or tensile strain reduces to 1% of the maximum stress or strain, respectively.
The longitudinal and transverse stress and strain pulses at the bottom of HMA layer according to LET analysis are presented in Figure 2 for a typical pavement section.
As can be seen, the response pulse in longitudinal direction generally consists of two compression zones and one tension zone, while in case of transverse response pulse, the HMA layer commonly experiences tensile stress or strain. In general, the shape of strain pulses computed by NonPAS program was very similar to those measured in full-scale pavement tests [9][10][11][12]. In this research, with respect to prediction of loading frequency for simulation of asphalt fatigue damage, only the tension zone of longitudinal and transverse pulses was considered for modeling.

Fitting of the Haversine Function to the Analytical
Response Pulse. Effect of loading waveform on the fatigue life of asphalt mixes can be explained by the energy put into the system per loading cycle [17]. It has been suggested that the energy is proportional to the area occupied by a load waveform in a stress (or strain) against time coordinate system [18]. On the other hand, full-scale test results showed that haversine function is a good representation of strain pulses in the longitudinal and transverse directions [9]. According to these observations, for each record of dataset (stress and stain pulses in both longitudinal and transverse directions at the bottom of asphalt layer for 5000 pavement sections), the haversine function was fitted to analytical response pulse such that the area under the haversine function was equal to the area under the analytical response pulse obtained by LET analysis. The normalized haversine function can be expressed as follows: where eff is the effective length and ( ) is the normalized value of the pulse at distance of . Since the analytical pulses have been determined in distance domain, duration of haversine pulse is also determined in distance domain which is called effective length ( eff ). Given the speed of moving wheel, the duration of haversine pulse as well as loading frequency can be obtained using where is the duration in time domain in second, eff is the effective length based on fitted haversine waveform to analytical pulse in centimeter, is the speed of moving wheel in km/h, and is the frequency of haversine waveform in Hertz.
Before fitting the haversine function to response pulse, the response pulse was normalized by dividing all data points to the maximum value. The area under the normalized haversine function ( ) can be obtained using the following integral: Given the area under the analytical response pulse ( ), the effective length ( eff ) of fitted haversine function can be determined by eff = 2 .
As a criterion for goodness-of-fit, the coefficient of determination ( 2 ) between fitted haversine function and analytical response pulse was determined. The frequency histogram of coefficients of determination ( 2 ) is given in Figure 3. As evident, the maximum and minimum values of 2 are 0.99 and 0.70, respectively. Figure 3 also indicates that the haversine function is fitted better to stress and strain pulses in longitudinal direction than transverse direction. This observation is in agreement with full-scale tests results [9].
Although the haversine waveform is more appropriate for representing the tensile responses at the bottom of asphalt layer (especially in transverse direction) and has been recommended by ASTM D7460, some standards such as AASHTO T321 and EN 12697-26 recommend the sinusoidal waveform to perform fatigue tests of asphalt mixes. In such cases, the where eff is the effective length of a sinusoidal waveform which has the area equal to a haversine pulse with effective length of eff . In the next sections of this paper, the haversine waveform is considered for modeling of tensile response pulses at the bottom of asphalt layer. It is obvious that all results can be generalized to sinusoidal waveform by using (8).

Effective Length ( ) of Horizontal Stress and Strain
Pulse. The statistical data relating to computed effective length ( eff ) for stress and strain pulses in both longitudinal and transverse directions is given in Table 2. As can be seen, in the case of both stress and strain pulses, the effective length in transverse direction is larger than that in longitudinal direction. The frequency of stress and strain pulses in both directions at two different speeds of 10 and 80 km/h are given in Table 3. Compared with the frequency of 5-10 Hz, which is commonly used in the asphalt mix laboratory test, the real frequency of loading may range from 0.30 to 77.5 Hz based on the vehicle speed, load specifications, and pavement design.
Relations between effective length of stress and strain pulses in longitudinal and transverse directions are illustrated in Figure 4. As can be seen, linear relationships exist between effective length of stress and strain pulses in both longitudinal and transverse directions. These relations are very useful and practical for computation of the effective length of different response pulses, where the effective length of one response pulse is known.
The ratio of effective length between transverse and longitudinal directions of strain pulse is about 3.14 ( Figure 4(b)). This value is in agreement with the findings of previous researchers who observed that the duration of transverse tensile strain was almost three times of what was measured in longitudinal direction [9,11].

Multivariate Adaptive Regression Splines (MARS)
3.1. Theory of MARS. The theory of Multivariate Adaptive Regression Splines (MARS) was developed by Friedman (1991) for solving regression-type problems [19]. The MARS technique has become particularly popular in the area of 6 The Scientific World Journal  The Scientific World Journal 7 data mining because it is categorized into nonparametric regression procedures that make no assumption about the form of functional relationship (e.g., linear and logistic) between the dependent and predictor variables. Instead, useful models (i.e., models that yield accurate predictions) can be derived even in situations where the relationship between the predictors and the dependent variables is nonmonotone and difficult to approximate with parametric models [20]. MARS divides the whole space of input variable into various subregions. It defines a different mathematical equation for each subregion. The fundamental idea of MARS is to use the combination of the linear truncated basis functions to approximate the model. Thus, the functions of MARS consist of single spline functions or the product of two or more of the truncated power functions to allow for the interactions. MARS model can be written as follows: where 0 is the coefficient of the constant basis function 0 ( ) = 1, ( ) is the th basis function, which may be a single spline function or product of two or more, is the coefficient of the basis function, and is the number of basis functions in the model. Each basis function, ( ), takes one of the three forms of the following: (1) a constant, (2), a hinge function ( − ) + or ( − ) − and (3) a product of two or more hinge functions. A product of two or more hinge functions can model interaction between two or more variables. The hinge functions have the following form: where is a constant, called the knot. . The forward phase is executed until one of the following conditions occurs: change in residual error is smaller than threshold or the maximum number of terms is reached. These parameters are specified by the user beforehand. Due to the nature of hinge functions, this search can be done more rapidly by using fast least-squares update technique or by using a heuristic method that reduces the number of parent terms considered at each step [19].

Learning in
The Backward Pass. The forward pass usually builds an overfitted model. So, a backward deletion phase is engaged to build a model with better generalization ability. In this phase the model is pruned by removing one least effective to find the best submodel. Model subsets are compared using the generalized cross-validation criterion (GCV) which is defined as follows: where is the number of basis functions in the model, denotes the fitted values of the current MARS model, denotes the number of data points and is the penalizing parameter. The numerator is the common residual sum of squares, which is penalized by the denominator, which accounts for the increasing variance in the case of increasing model complexity. The penalizing parameter can be chosen arbitrarily. A conventional value is = 4. A smaller generates a larger model with more basis functions; a larger creates a smaller model with less basis functions [21].
At the end of the backward phase, from those "best" models of each size, a model with lowest GCV value is selected as the final one. One of the backward phase's benefits over the forward phase is that at any step it can select any term to delete, whereas the forward phase, at each step, can only see the next pair of terms. The forward phase adds terms in pairs; however, the backward pass typically discards one side of the pair. Thus, terms are often not seen in pairs in the final model [19].

Prediction of by Means of MARS.
For developing some equations to predict the effective length ( eff ) of stress and strain pulses at the bottom of asphalt layer in both longitudinal and transverse directions, on the basis of MARS method, STATISTICA program was employed. The degree of interactions between predictors was assumed to be 3, since higher values had no effect on improving the model. Maximum number of basis function and penalizing parameter were set to 25 and 4, respectively. The where 1 is thickness of asphalt layer (cm), 2 is thickness of base layer (cm), 3 is thickness of subbase layer (cm), 1 / 2 is the ratio of asphalt resilient modulus to base resilient modulus, 2 / 3 is the ratio of base resilient modulus to subbase resilient modulus, 3 / 4 is the ratio of subbase resilient modulus to subgrade resilient modulus.
Regression statistics for (12) to (15) are given in Table 4. Capability of (12) to (15) to predict the effective length of tensile stress and strain pulses in both longitudinal and transverse directions at the bottom of asphalt layer is illustrated in Figure 5. As evident, these equations give good accuracy to predict effective length. For all pulses, the error percentage decreases with increasing the effective length. It can be therefore concluded that for high values of effective length the proposed equations have sufficient accuracy; however, for low values of effective length (pavements with thick asphalt layer), proposed equations are of limited accuracy.

Theory of ANN.
In order to predict the effective length accurately, the Artificial Neural Network (ANN) method was employed. The ANN approach is a computer methodology which attempts to simulate some important features of the human nervous system; in other words, the ability to solve problems by applying the information gained from the past experiences to new problems or case scenarios. Analogous to a human brain, an ANN uses many simple computational elements, named artificial neurons, connected by variable weights [22]. A typical artificial neuron is illustrated in Figure 6. A neural network can be trained to perform a particular function by adjusting the values of the connections (weights) between the elements. Neural networks are trained so that a particular input leads to a specific target output. The network is adjusted based on a comparison of the output and the target until the network output matches the target.
Typically many such input/target output pairs are used to train a network.
A feed-forward backpropagation neural network ( Figure 7) is a kind of ANN which is useful in addressing problems requiring recognition of complex patterns and performing nontrivial mapping function [23].
The training of a feed forward neural network using a backpropagation algorithm involves the following two phases [24,25].
(i) Forward Phase. During this phase, the free parameters of the network are fixed, and the input signal is propagated through the network layer by layer. The forward phase ends with the computation of an error signal using the following: where is the desired response and is the actual output produced by the network in response to the input .
(ii) Backward Phase. During this second phase, the error signal is propagated through the network in the backward direction, hence, the name of the algorithm. It is during this phase that adjustments are applied to the free parameters of the network so as to minimize the error in a statistical sense.

Prediction of by Means of ANN.
Several implementations of the backpropagation algorithm are possible. In the present study, the Levenberg-Marquardt algorithm was adopted due to its efficiency in training networks. This implementation is readily available in Matlab software within its neural network toolbox. The testing (30%), cross validating (10%), and training (60%) sets for ANN training procedure were selected randomly from the established dataset. Several networks with one hidden layer and different numbers of neurons in each hidden layer (5 to 20) were explored to determine the optimal architecture of BPN. The optimal architecture of ANN with the minimum possible size and acceptable accuracy was found to be 7-15-3 (15 hidden neurons) architecture. Hyperbolic tangent sigmoid and linear transfer functions were used for the hidden layer and output layer, respectively. More details to implement the proposed ANN for prediction of the effective length of different tensile responses at the bottom of asphalt layer are given in Appendix.
The performance of ANN modeling to predict the effective length ( eff ) of longitudinal and transverse stress and strain pulses related to training and testing sets is demonstrated in Figure 8. The frequency histogram of ANN output residuals is presented in Figure 9 and is compared to normal distribution. For prediction of frequency of different responses at the bottom of asphalt layer based on ANN, a program was developed using the macrocapability of Microsoft Excel. This program enables user to determine the frequency of stress and strain pulses at the bottom of asphalt layer with respect to input parameters (moving speed, thickness of layers, resilient modulus of different layers, and contact radius of moving load). User interface of this program is represented in Figure 10.

Parametric Analysis.
A standard pavement section composed of a HMA layer, an aggregate base, an aggregate subbase, and subgrade soil was assumed to investigate the effect of various parameters on the effective length of the tensile response pulses. These controlling parameters include the radius of contact area, thickness of different layers, and ratio between the modules of each layer and the immediate succeeding layer below. The thickness and the resilient modulus of different layers are represented in Table 5. The standard contact radius for parametric analysis was considered to be 10 cm.
In order to study the effect of different parameters on the effective length of tensile stress and strain pulses, the ANN  x pi x pN · · · · · · · · · · · · · · · · · · h j o k w Figure 7: A three layer feed-forward backpropagation network architecture. method was employed. Using ANN, the effective length of the pulse was computed by changing the desired parameter. The results of the parametric analysis are illustrated in Figure 11. The effect of subbase thickness ( 3 ) is not shown, because, similar to 2 , it has a very slight effect on effective length. As can be seen in Figure 11, some parameters including the ratio of resilient modulus of subbase to subgrade ( 3 / 4 ) and the thickness of base and subbase layers have a slight effect on the effective length of response pulses, especially in longitudinal direction. The most prominent factors that affect the effective length of tensile response pulses are as follows: (i) contact radius of wheel, (ii) thickness of asphalt layer, (iii) the ratio of resilient modulus of asphalt layer to base layer, (iv) the ratio of resilient modulus of base to subbase layer.
All of these factors have a direct effect on the effective length of tensile response pulses at the bottom of asphalt layer. The effective length of stress and strain pulses in both longitudinal and transverse directions increase with the increase of these parameters.
It is worth nothing that the effect of temperatures is similar to the reciprocal of 1 / 2 ratio. In fact when the temperature increases, the 1 / 2 ratio decreases and vice versa. Full-scale tests results show that the duration of tensile response pulses increases with decreasing the asphalt layer temperature which is in agreement with the effect of 1 / 2 ratio on effective length [10][11][12].

Conclusions
The following conclusions can be drawn from this study.
(1) Haversine function was fitted to stress and strain pulses in longitudinal direction better than transverse direction.
(2) Response pulse in longitudinal direction generally consists of two compression zones and one tension zone, while in case of transverse response pulse, the HMA layer commonly experiences tensile stress or  (3) There are strong correlations between the effective length of stress and strain pulses in both longitudinal and transverse directions. Using these relations and measuring/computing the effective length for one of the stress or strain pulses, the effective length for other responses can be obtained easily.
(4) MARS approach predicts the effective length with a good accuracy when the effective length is large (e.g., wide-base tires or pavement sections with thick   asphalt layer). A more accurate method is therefore needed when the effective length is small.
The most prominent factors that affect the effective length of tensile response pulses are the contact radius of wheel, the thickness of asphalt layer, the ratio of resilient modulus of asphalt layer to the base layer, and the ratio of resilient modulus of base to sub-base layer. The effective length of tensile response pulses in both longitudinal and transverse directions increases with the increase of these factors.