An overview of a unified theory of dynamics of vehicle–pavement interaction under moving and stochastic load

This article lays out a unified theory for dynamics of vehicle–pavement interaction under moving and stochastic loads. It covers three major aspects of the subject: pavement surface, tire–pavement contact forces, and response of continuum media under moving and stochastic vehicular loads. Under the subject of pavement surface, the spectrum of thermal joints is analyzed using Fourier analysis of periodic function. One-dimensional and two-dimensional random field models of pavement surface are discussed given three different assumptions. Under the subject of tire–pavement contact forces, a vehicle is modeled as a linear system. At a constant speed of travel, random field of pavement surface serves as a stationary stochastic process exciting vehicle vibration, which, in turn, generates contact force at the interface of tire and pavement. The contact forces are analyzed in the time domain and the frequency domains using random vibration theory. It is shown that the contact force can be treated as a nonzero mean stationary process with a normal distribution. Power spectral density of the contact force of a vehicle with walking-beam suspension is simulated as an illustration. Under the subject of response of continuum media under moving and stochastic vehicular loads, both time-domain and frequency-domain analyses are presented for analytic treatment of moving load problem. It is shown that stochastic response of linear continuum media subject to a moving stationary load is a nonstationary process. Such a nonstationary stochastic process can be converted to a stationary stochastic process in a follow-up moving coordinate.

Abstract This article lays out a unified theory for dynamics of vehicle-pavement interaction under moving and stochastic loads. It covers three major aspects of the subject: pavement surface, tire-pavement contact forces, and response of continuum media under moving and stochastic vehicular loads. Under the subject of pavement surface, the spectrum of thermal joints is analyzed using Fourier analysis of periodic function. One-dimensional and two-dimensional random field models of pavement surface are discussed given three different assumptions. Under the subject of tire-pavement contact forces, a vehicle is modeled as a linear system. At a constant speed of travel, random field of pavement surface serves as a stationary stochastic process exciting vehicle vibration, which, in turn, generates contact force at the interface of tire and pavement. The contact forces are analyzed in the time domain and the frequency domains using random vibration theory. It is shown that the contact force can be treated as a nonzero mean stationary process with a normal distribution. Power spectral density of the contact force of a vehicle with walking-beam suspension is simulated as an illustration. Under the subject of response of continuum media under moving and stochastic vehicular loads, both time-domain and frequency-domain analyses are presented for analytic treatment of moving load problem. It is shown that stochastic response of linear continuum media subject to a moving stationary load is a nonstationary process. Such a nonstationary stochastic process can be converted to a stationary stochastic process in a follow-up moving coordinate.
Keywords Vehicle-pavement interaction Á Random field Á Continuum medium Á Spectral analysis Á Green's function Á Linear system The investment of the United States in the nation's transportation infrastructure alone (highways, bridges, railways, and airports) amounted to $7 trillion by 1999. To preserve infrastructure longevity in a cost-effective manner, the research in pavement design and infrastructure management has been growing rapidly in recent years. From a pavement design point of view, pavement response, damage, and performance are essentially the result of long-term vehicle-pavement interaction. When vehicle speed is low, the dynamic effect of vehicular loads on pavements is insignificant. However, with the promotion of high-speed surface transportation in the world, this dynamic effect must be taken into account to develop more rational pavement design methods. For instance, real causative mechanisms that lead to fatigue damage of pavement material might be frequency dependent. From an infrastructure management point of view, vehicle-pavement interaction has a profound impact on the way that existing technologies of structural health monitoring, environmental vibration mitigation, nondestructive testing and evaluation, and vehicle weight-in-motion are to be improved, innovated and implemented. For instance, modern high-speed surface transportation systems are normally accompanied by rises in levels of noise and vibration that may cause a significant detrimental effect to the ecology. Vehiclepavement interaction-induced structure-borne and groundborne vibrations emit and propagate toward some extent. Residents may experience hardship from uncomfortable vibration, and high-precision equipment may suffer from malfunctioning to irreparable damage. To mitigate noise and vibration in surrounding areas of roadway, it is necessary to investigate predominant frequencies of vehiclepavement interaction, the source of vibration, so as to develop effective noise and vibration countermeasures. The study of vehicle-pavement interaction also plays a critical role in developing better inversion algorithms for nondestructive testing and evaluation of transportation infrastructure. In addition, taking into account dynamic effects of vehicle vibration caused by rough surface may considerably improve accuracy and reliability of weigh-in-motion systems, which measure a vehicle's weight as it travels at a normal speed. Vehicle-pavement interaction is not only a central problem to pavement design, but also has a profound impact on infrastructure management, vehicle suspension design, and transportation economy. Figure 1 shows a central role of vehicle-pavement interaction played in a wide variety of applications. From a vehicle-design point of view, a designer needs to consider vehicle vibration and controllability, which affect ride quality and vehicle maneuver. Vehicle-pavement interaction also has a huge economic impact. As pavement performance deteriorates and roadway surface gets rougher, both operation costs (fuel, tire wear, and routine maintenance) and roadway maintenance cost of the vehicle increase dramatically, accompanied by decreasing transportation productivity. There is no doubt that there is a great and urgent need for fundamental research on vehicle-pavement interaction due to rapid deterioration of huge highway infrastructure nationwide, tight maintenance budget, and the key role played by vehicle-pavement interaction.
Since the American Association of State Highway Officials (AASHO) road test, the fourth power law has been widely used by pavement engineers to design highway and airfield pavement and to predict the remained life and cumulated damage of pavements [1][2][3]. Besides the damage caused by static loads, dynamic loads may lead to additional pavement damage. A consequence of a high power in the damage law is that any fluctuation of pavement loading may cause a significant increase in the damage suffered by pavement structures. A number of recent filed measurements and theoretical investigations showed that vehicle vibration-induced pavement loads are moving stochastic loads [4][5][6][7][8]. The researchers concluded that vibrations of vehicles were related primarily to pavement surface roughness, vehicle velocity, and suspension types [9][10][11][12][13].
Estimation of pavement damage caused by dynamic loads varies anywhere from 20 % to 400 % [14]. The smaller estimates are based on the assumption that peak dynamic loads (and hence the resulting pavement damage) are distributed randomly over the pavement surface. The larger estimates are based on the assumption that vehicles consistently apply their peak wheel loads in the same areas of the pavement. Theoretical studies by Cole [15] and Hardy and Cebon [16] confirmed that for typical highway traffic, certain areas of the pavement surface always suffered the largest wheel forces, even when the vehicles had a wide range of different suspension systems, payloads, and speeds.
Complicated relationships exist among vehicle suspension, dynamic wheel loads, pavement response, and damage [17][18][19]. On one hand, it has been known for years on how to manufacture automobiles operating properly on a variety of pavement surfaces. On the other hand, however, the effect of vehicle design on pavement has not been thoroughly studied. For instance, Orr [20] stated in his study that comparatively little was known about the influence of suspension design on pavement in the automobile industry yet.
A number of obstacles exist in revealing vehicle-pavement interaction. A theoretical foundation universally applicable to the involved specific problems is no doubt very attractive. It will not only provide a guide for experimental study and validation, but will also enable better design and maintenance of vehicles and pavements. As such, this article provides an overview of a unified theory for dynamics of vehicle-pavement interaction under moving and stochastic loads. This article covers three major aspects of the subject: pavement surface, tire-pavement contact forces, and response of continuum media under moving and stochastic vehicular loads developed by Sun and his associates (Fig. 2). The remainder of this article is organized as follows. Section 2 addresses mathematical description of pavement surface roughness. Section 3 studies contact force generated due to vehicle-pavement interaction, which is the source to both vehicle dynamics and pavement dynamics. Section 4 investigates pavement response under moving stochastic loads. Section 5 discusses difficulties and deficiencies in the existing research, and projects further research and applications. Section 6 makes concluding remarks.

Pavement surface
Irregularities in pavement surface are often called roughness. All pavement surfaces have some irregularities. The National Cooperative Highway Research Program defines ''roughness'' as ''the deviations of a pavement surface from a true planar surface with characteristic dimensions that affect vehicle dynamics, ride quality, and dynamic pavement loads [21].'' Factors contributing to roughness, in a general sense, include vertical alignment, cracks, joints, potholes, patches and other surface distresses.
Newly constructed pavements may be poorly finished or have design features such as construction joints and thermal expansion joints, which can be main sources of vehicle vibration [22,23]. Measurements show that pavement roughness can be modeled as a random field consisting of different wavelength. Considerable effort has been devoted to describing pavement surface roughness [24][25][26]. The peak-to-valley measurements, the average deviation from a straight edge, and the cockpit acceleration are several distinct approaches that have been suggested for pavement surface characterization [27]. Practical limitations of these three approaches could be found in the study of Hsueh and Penzien [28].
Spectrum of a deterministic function referring to its Fourier transform reveals sinusoid components that a deterministic periodic or nonperiodic function is composed of. Spectrum analysis of a deterministic function f ðtÞ can be obtained if and only if this deterministic function satisfies Dirichlet condition and absolute integrability condition [29]: Pavement surface as a random field of elevation does not decay as spatial coordinates extend to infinity; therefore, it does not satisfy inequality (1). For this reason, taking Fourier analysis to a sampled pavement surface does not make too much sense in theory. However, the correlation vector of random field does decay as spatial coordinates extend to infinity. Hence, applying Fourier transform to correlation vector of pavement surface, commonly known as power spectral density, does serve the purpose of spectrum analysis in theory [30]. As such, mathematical description of pavement surface takes place under the theoretical framework of spectral analysis of stochastic process. The following subsections consist of a general mathematical framework for describing pavement surface roughness, a review of statistical description of pavement surface roughness, a description of periodic joints of cement concrete pavement, and three different hypotheses when projecting a one-dimensional road profile to a twodimensional pavement surface.

One-dimensional description
Given the manifest complexity exhibited in pavement surface as a random field, making some assumptions becomes indispensable to simplify mathematical description of pavement surface. Commonly used assumptions on surface roughness are that surface roughness is an ergodic and homogeneous random field with elevation obeying Gaussian distribution [31][32][33]. The assumption of ergodicity makes sure that the temporal average of a sample of stochastic process equals to the statistical mean of stochastic process, which enable one to obtain the statistical characteristics of a stochastic process by measuring only a few samples. The assumption of homogeneity ensures random property of surface roughness is independent to the sites measured. When a vehicle travels at a constant speed, a homogeneous random filed is converted to a stationary stochastic process. The assumption of Gaussian distribution ensures that transformation of such randomness through a linear system is still Gaussian.
Road profile along longitudinal direction (i.e., direction of travel) is often an exchangeable term of pavement roughness. The simplest model describes the pavement surface as a cylindrical surface defined by a single longitudinal profile n(x). Figure 3 shows an illustrative example of road profile. Instead of varying with time, the elevation n of the surface varies with respect to longitudinal distance x along the direction of travel. In the spatial domain, lowfrequency components correspond to long wavelengths while high frequency components correspond to short wavelengths. Let n(x) be a zero mean, homogeneous, Gaussian random field. Its probabilistic structure can be completely described by either the autocorrelation function or the power spectral density (PSD).
According to the Winner-Khintchine theory [31], the following expressions constitute a pair of Fourier transform: where X represents the distance between any two points along the road. Wavenumber spectrum, S nn ðXÞ, is the direct PSD function of wavenumber X , which represents spatial frequency defined by X = 2p/k where k is the wavelength of roughness. Under the assumption of homogeneity, spatial autocorrelation function R n ðXÞ is defined by for any x 1 and X. Here E½Á is the expectation operator and can be calculated by This expectation, differing from statistical mean of a random process, is temporal average of a sample of stochastic process.

PSD roughness
As a one-dimensional random field, roughness can be characterized in the spatial domain and the wavenumber domain. An analytic form of PSD roughness is often desired for theoretical treatment [34,35]. Spectral analysis of pavement roughness has been the subject of considerable research for years. Many function forms of PSD roughness have been proposed to characterize bridge deck surface [36], rail-track surface [33,37], and airfield runway surface [24,38].
Sayers et al. [39,40] suggested a power function for PSD roughness: Early studies on PSD of longitudinal profiles of runways and roads are a special case of (5) for a ¼ À2 [41][42][43][44]. Although power function (5) is convenient for parameter estimation and design purpose, it creates mathematical difficulties at X ¼ 0, where S nn ðXÞ becomes infinite. For this reason, two distinct function forms have been proposed in later studies.
Another form is to use piecewise functions. For instance, (11) was suggested by the International Organization for Standardization (ISO) to cover different frequency ranges [49][50][51][52].
where X a represents reference frequency, and X b represents cut-off frequency. Iyengar and Jaiswal [33] provided a similar split power law in the form of (11) for rail-track in which C sp = 0.001(m -2 cycles/m) and w 1 = w 2 = 3.2.
Marcondes et al. [22,23] proposed an exponential piecewise function: In all aforementioned PSD roughness, A, A 1 , A 2 , p, q, k, w 1 , w 2 , C sp , a b, S 0 , and X 0 are the real and positive parameters estimated from field test.
Although various PSD functions have been proposed for fitting the measured PSD roughness, it has been found that most pavement profiles including road surface, runway surface, and rail-track surface exhibit very similar trends of PSD curves. Some researchers have even reduced the PSD curves to only two families: one for rigid (cement concrete pavements), and the other for flexible (asphalt concrete pavements) [25,26].

Effect of periodic joints of rigid pavement
Thermal expansion joints are commonly set on rigid pavement surface to reduce thermal stress in pavement. Thermal joints of rigid pavement can be a significant source of excitation of vehicle vibration, especially when the joint sealing material gets lost. This could cause undesirable riding qualities at a high speed [43]. Figure 4 shows an idealized rigid pavement surface with periodic joints. An actual rigid pavement surface might be considered as a combination of this periodic profile and a random field.
In the spatial domain, pavement roughness n(x) in Fig. 4 can be described by a periodic function within a period [-d/2,d/2]. nðxÞ ¼ 0; x 2 ½Àd=2; ÀD=2Þ cðxÞ; x 2 ½ÀD=2; D=2: 0; x 2 ðD=2; d=2 Here spatial period along the longitudinal direction, d, represents slab length of rigid pavement, D and h, respectively, represent maximum width and depth of joint, c(x) is named shape function of joint, which is a symmetrical function and describes the shape of the joint. With this description, pavement roughness n(x) becomes an even function with period d.
A period function can be legitimately expanded using Fourier series, which is sometimes called frequency spectrum analysis of periodic function [53]. As such, coefficients of the expanded Fourier series of n(x) a n and b n (n = 0,1,2,)ÁÁÁ are given by cðxÞ cos 2npx d dx; ð14aÞ b n ¼ 0; n ¼ 0; 1; 2; . . . (for nðxÞ is an even function): Wavenumber spectrum of roughness is thus expressed as In which X n represents discrete wavenumber and S nn ðX n Þ is discrete wavenumber spectrum. It is evident from (18) that different shapes of joints determine different discrete spatial spectrum. Four types of shape function, namely rectangular curve, parabolic curve, cosine curve, and triangular curve have been observed and investigated by Sun [54]. Table 1 gives the discrete spatial spectrums of these four types of shape functions computed from (15a). The comparison between these discrete spectrums for specified rigid pavement with slab length 5 m and joint width 0.03 m is plotted in Fig. 5 where S ¼ hD=d, X n ¼ 2pn=d, n ¼ 0; 1; 2; Á Á Á, and D=d ¼ 0:006. From this figure it can be seen that the spectrum of rectangular joint is the greatest in four types of joints, while the spectrum of triangular joints is the smallest. The ratio of the magnitude of four types of joints based on S nn ðX 0 Þ approximates to 4:1.78:1.62:1. Because S 2 = (hD/d) 2 , the effect of joint of rigid pavement with long slab length d is much less than that of the pavement with short slab length.

Two-dimensional description
One-dimensional random field model of pavement surface is adequate for two-wheel vehicles, such as bicycles and motor-cycles but inadequate for cars and trucks having two or more wheel per axle. Actual highway and airfield pavement consists of a two-dimensional surface of finite width, a nominal camber and grader. The elevation of pavement surface exhibit random fluctuations about the nominal geometry and therefore should be more accurately treated as a two-dimensional random field, n(x, y), with space coordinates x and y as the indexing parameters as shown in Fig. 6.
When ignoring those isolated large fluctuations such as potholes, fluctuations of pavement surface can be approximate by a homogenous, Gaussian random field with a zero mean [49,55]. Probabilistic structure of a two-dimensional random field can be completely defined either by twodimensional autocorrelation function or by two-dimensional PSD, S nn ðX; KÞ, where X = x 2x 1 , Y = y 1 -y 2 , X and K represent wavenumber in x-axis and y-axis directions, respectively. Here S nn ðX; KÞ is defined as a double-sided Fourier transform of autocorrelation function (16).
RðX; YÞe ÀiðXXþKYÞ dXdY: 1 À ðDX n =pÞ 2 2 S nn ðX 0 Þ: The determination of the two-dimensional autocorrelation function or two-dimensional PSD relies on the measured elevation of the entire pavement surface, which is a possible but formidable task in data acquisition, numerical computation and data storage for hundreds of thousands of miles of roadways. Efforts have, therefore, been made to simplify the two-dimensional random field such that it can be uniquely generated from the onedimensional random field models by imposing certain hypothetic properties to pavement surface [49,50,55].

The hypothesis of isotropy
Consider a two-dimensional random field model of pavement surface as shown in Fig. 7. There are two parallel wheel paths separated by a constant distance Y along the transverse direction. Longitudinal profiles along each wheel path can be derived from the two-dimensional random field n(x, y) as follows.
From (16) the autocorrelation functions between A and B and between D and C are, respectively, given by and Since nðx; yÞ is homogenous, Crosscorrelation functions, R AC ðX; YÞ and R CA ðX; YÞ, are even functions of X and Y: Now, assume that n(x, y) be an isotropic random field [49,56]. The property of isotropy requires that for any profile making an angle h with x-axis, the following condition holds: Rðq cos h; q sin hÞ ¼ Rðq; 0Þ: From (16), (19), and (21), it follows that Since autocorrelation function and PSD form a Fourier pair, cross-PSD is given by Equations (19) and (23) are general results for spectral analysis under the hypothesis of isotropy. When using the definition of a normalized cross-PSD developed by Dodds and Robson [49], cross-PSD can also be written as where gðXÞ is coherence function between direct-PSD and cross-PSD of two parallel wheel paths, which is always not  greater than 1 [57]. Heath [58] further presented a closedform gðXÞ where J 1 ðÁÞ is the first-order Bessel function and X ! 0.

The hypothesis of uncorrelation
Parkhilovskii [55] proposed the hypothesis of uncorrelation. It assumes that the two profiles of two parallel wheel paths, n(x, y 1 ) and n(x, y 2 ) nðx; y 2 Þ, can be derived from two uncorrelated random fields, that is It then follows that Thus, direct-and cross-PSD can be expressed in terms of direct PSDs of fðxÞ and cðxÞ. Robson [59] has examined the physical and mathematical basis of Parkhilovskii model. He concluded that Parkhilovskii model can be made compatible with isotropic model for a profile-pair description of pavement (i.e., two parallel railway tracks), and may be used where isotropy assumption is not valid.

The Hypothesis of Shift
Sun and Su [60] proposed the hypothesis of shift to construct two-dimensional PSD using one-dimensional PSD. According to homogeneity in (19), autocorrelations of n(x, y 1 ) and n(x, y) equal to each other. Sun and Su [60] assumed nðx 1 þ X; y 1 þ YÞ ¼ nðx 1 þ X þ s; y 1 þ YÞ. In other words, the elevation of one wheel path is equal to the elevation of a parallel but shifted wheel path with a spatial lag s. With where spatial lag s is a parameter that can be estimated from field test of autocorrelation functions of parallel wheel path. Some properties about the spatial lag can be Under the hypothesis of shift, since PSD is the Fourier transform of correlation function, we have Replacing R nn ðX þ sÞ in (29) by (30) we obtain Reversing the order of integration and using the convergence concept of a generalized function we have the inner integral Here the following integral is used Substituting (32) into (31) we get Here the property of Dirac delta function (Lighthill 1958) is used in the derivation of (34). In summary, several used assumptions related to pavement surface are homogeneity, ergodicity and normal distribution. These have been supported by a number of measurements of pavement surface. However, there are still counterexamples demonstrating the existence of pavement surfaces that do not satisfy the homogeneous and ergodic properties [62]. In those situations, pavement surfaces are perceived by a vehicle as nonstationary stochastic processes [47,63,64]. It is worthwhile noting that nowadays with technology advancement, one can measure and obtain entire two-dimensional topology of pavement surface using remote sensing, line scanning laser, and synthetic aperture radar. For more accurate applications, it is no longer necessary to use hypothetical random field model of pavement surface.

Contact force on vehicle-pavement interface
Wheel loads in specifications of highway and airport pavement design are presently treated as a static load pressure with uniform distribution [65] PðrÞ where r 0 is the radium of circle distribution of wheel loads and P 0 represents the total loads applied. Such a simplification of actual vehicle/aircraft load in design can be adequate when the speed of travel is low and pavement surface is smooth. However, when the speed of travel is high and pavement surface is uneven due to deterioration, influence of dynamic load generated by vehicle-pavement interaction can become very significant. As a matter of fact, pavement damage is directly resulting from a long-term effect of dynamic traffic loading. Therefore, contact force on vehicle-pavement interface needs to be considered in pavement design to more closely reflect actual loading condition.
All pavement surfaces exhibit irregularities that cause vehicle vibration, which, in turn, results in moving stochastic contact forces on pavement [66]. At low speed of travel vehicle vibration is insignificant but at high speed of travel vehicle may vibrate significantly due to poor pavement surface conditions. It does not only reduce the ride quality but also generates additional damage to vehicle and pavement structures beyond static load.
The interaction between vehicles and pavements has been investigated extensively by the automobile industry for improving ride quality and reducing the mechanical fatigue of vehicle [5,25,50,[67][68][69]. Hedric [9], Abbo et al. [10], and Markow et al. [70] identified some critical factors that affect dynamic loads on pavement. These factors include vehicle and axle configuration, vehicle load, suspension characteristics (stiffness, damping), speed of travel, pavement roughness, faults, joint spacing and slab warping. From the perspective of pavement design and maintenance, however, little has been known about the effects of vehicle vibration on response, damage and performance of pavement structures [20,54]. Although some efforts have been devoted to the measurement and prediction of dynamic wheel loads [4,6,11,12,71,72], few of them provides a complete theoretical foundation for describing contact force induced by vehicle-pavement interaction. This can be a crucial element leading to a more precise dynamic analysis of pavement structures [54,73,74].
Of three approaches for studying vibration-generated contact forces, namely, analytic, experimental, and numerical simulation methods, only the last two have been used widely. Experimental method offers real results, it is however very costly and limited by safety requirements. Moreover, results from experimental method only suits in some degree for specific test conditions (e.g., vehicle and pavement types, speed of travel) [54]. In many circumstances, numerical simulation associated with a limited field tests serves as the most prevailing approach.
Numerical simulation has the advantage of capable of extrapolating experimental results over a range of test conditions where the experiment would be too dangerous or too expensive. To study the contact forces between vehicle and pavement using numerical simulation, a vehicle must be simplified to a vehicle model so as to simulate the real operation conditions of the vehicle. Then, based on the vehicle model associated with measured or simulated pavement surface roughness, dynamic contact forces can be generated using a validated vehicle simulation program.
3.1 Pavement surface as a source of excitation to vehicle As far as the vehicle vibration is concerned, pavement surface serves as a source of excitation when the vehicle travels along the road. As such, spatial fluctuation of pavement surface gets converted to temporal random excitation of vehicle suspension system. In this article, we assume that the speed of travel is constant. Different stochastic process of excitation can be resulting from different hypotheses on two-dimensional random field of pavement surface, which are present in this section.

General two-dimensional random field
Let x ¼ ðx; yÞ be the fixed coordinates in which the spatial random field is defined. Let x 0 ¼ ðx 0 ; y 0 Þ be the moving coordinates connected with the moving vehicle and travel at the constant speed v (vector quantity) with respect to the fixed coordinates x. Denote R n ðx 1 ; x 2 Þ the correlation function that describes the spatial random field nðxÞ. Denote R g ðx 0 1 ; x 0 2 ; t 1 ; t 2 Þ the correlation function that describes the temporal random excitation gðx 0 ; tÞ. Since we have From the assumption mentioned above, nðxÞ is a homogeneous random field. Accordingly, its correlation function R n ðx 1 ; Since vehicle velocity is a constant vector, the homogenous random field gðx 0 ; tÞ in spatial domain is converted to a stationary stochastic process in time domain. In other words, (38) can be replaced by where X 0 ¼ x 0 2 À x 0 1 and s ¼ t 2 À t 1 . Now we define S n ðXÞ as the spatial PSD of two-dimensional random field, i.e., The PSD of random excitation gðx 0 ; tÞ can be expressed as An overview of a unified theory of dynamics 143 S g ðX; xÞ ¼ S n ðXÞdðx À XvÞ: Equation (41) can be proved as follows. Since Applying Eqs. (41) to (42) gives This is the same as Eq. (39). It is necessary to represent the PSD of random excitation gðx 0 ; tÞ in terms of angular frequency x and spatial lag X 0 , i.e., If the direction of vehicle velocity vector is the same as that of the x-axis, then the following transformation applies: Substituting (45) into (44), (9) is converted to where X ¼ ðX; YÞ and X 0 ¼ ðX 0 ; Y 0 Þ. Corresponding to (40), spatial correlation function R n ðXÞ is the inverse Fourier transform of the spatial PSD S n ðXÞ, i.e., where X ¼ ðX; KÞ. Combining Eqs. (47) and (46), we obtain In the derivation of (13), the following equation is used.
Equation (48) is a general result applicable to any form of two-dimensional random field. Under the restriction of the specified vehicle systems, we can represent (48) furthermore. Figure 8 shows a plan view of a vehicle with N tires. We now consider the cross-PSD of between i th and j th tires. It is clear that if i th and j th tires are located in the same side of the vehicle, either in the left side or in the right side, then we should have X 0 = X ij and Y = 0. If i th tire is located in right side and j th tire is located in left side, or vise verse, then we should have X 0 = X ij and Y = 0.
Given Fig. 7 condition, (48) can be written in terms of spatial PSD roughness: Since the correlation function R n ðXÞ and PSD function S n ðXÞ of the two-dimensional random field are both known from (40) and (47), there is no obstacle to compute the excitation PSD of (40) and (41).

Two-dimensional random field under hypothesis of isotropy
When pavement surface is treated as a two-dimensional random field under hypothesis of isotropy [56,58], its twodimensional PSD can be represented in terms of onedimensional PSD. In other words, in light of (26) and the definition of two-dimensional PSD (40), the two-dimensional PSD can be constructed from the following Fourier transform: Notice that Here, the general integration and the property of Dirac delta function [75] Z are used in the derivation of Eq. (53).
Comparing (53) with (27), it is clear that the left side of (53) is just the PSD given by S AC ðXÞ. Since it has been proved that S AC ðXÞ ¼ gðXÞS nn ðXÞ where gðXÞ is an ordinary coherence function and S nn ðXÞ is one-dimensional PSD roughness. Hence, (40) and (41) can be further expressed as

Two-dimensional random field under hypothesis of uncorrelation
When pavement surface is treated as a two-dimensional random field under hypothesis of uncorrelation [55], the integration of the left side of (53) becomes either S AB ðXÞ if i and j are located on the same side of the vehicle, or S AC ðXÞ if i and j are located on the different side of the vehicle, which are further given by (30a) and (30b), respectively. Without any difficulty, it is straightforward to see that (40) and (41) can be further written as

Two-dimensional random field under hypothesis of shift
When pavement surface is constructed from a twodimensional random field under hypothesis of shift, (40) and (41) can be written as

Vehicle models
Here, we consider vehicle models for simulating dynamic contact forces. Three kinds of vehicle models are commonly used nowadays, i.e., quarter-vehicle, half-vehicle, An overview of a unified theory of dynamics 145 and full-vehicle models. Making a reasonable choice between a complex vehicle model and a simple vehicle model relies on the nature of the problem. The advantage of the complex vehicle model is that it can accurately predict the acceleration response at different locations of the vehicle. The complex vehicle model is thus suitable for studying influences of vehicle vibrations on human body and the fatigue life of vehicles caused by vehicle vibrations. However, the complex vehicle model is highly sophisticated and requires detailed input and long execution times even for simple problems. Because vehicle sizes and loads vary greatly, it is more difficult to select representative parameter values for the complex vehicle model than for the simple vehicle model, which increases the difficulty to compare computer simulation results conducted by different research groups. This should be always kept in mind. In Sect. 3.1, pavement surface is converted into a stationary stochastic process as input to excite vehicle vibration when the vehicle travels at a constant speed. To proceed further, we also assume that vehicle's suspension system is linear, under which the equations of motion of vehicle suspension systems for small oscillations can be derived using the Lagrange equations.
Let Y be a set of M independent generalized coordinates (e.g., the absolute displacement of components of vehicle systems) that completely specify the configuration of the system measured from the equilibrium position. Then, the kinetic energy T, potential energy U, and dissipative energy D, can be expressed as where M abs , K abs and C abs are mass, stiffness, and viscous damping matrices with respect to absolute displacement vector Y ¼ fY 1 ; Y 2 ; . . .; Y M g T , respectively. The Lagarange equations of motion are Under the assumption that all components of the vehicle system be linear, (61) and (62) yield a set of simultaneous second-order differential equations with constant coefficients: where g and _ g are, respectively, the pavement surface elevation vector and its derivative process. Clearly, here we assume that the contact between tires and pavement surface is in the form of point distribution.
Without loss of generality, let components of the vehicle system directly contacting with ground in the moving coordinates x 0 be numbered from 1 to N (see Fig. 7), where N represents the total number of tires. Let Equation (63) can be rewrote as where M, K and C are mass, stiffness, and viscous damping matrices with respect to relative displacement Z ¼ fZ 1 ; Z 2 ; . . .; Z M g T , respectively. Since the total number of degree of freedom is assumed to be M, and g is Ndimensional temporal excitation vector representing pavement surface-induced displacements in coordinates in whichZ andg are, respectively, the Fourier transform of Z and g, i.e., Multiplying the inverse matrix of ½Àx 2 M þ ixC þ K, it is straightforward to see The primary aim of establishing vehicle models is to find the scalar response of ZðtÞ as well asZðxÞ. This can be accomplished by using stochastic process theory of linear systems. The dynamic characteristics of the vehicle at angular frequency x are defined by M 9 N frequency response function (FRF) matrix, H(x), of the vehicle system, where H(x) = [H(x) ij ], i ¼ 1; 2; . . .; M for the M outputs and j = 1,2,….,N for the N inputs of pavement excitation to tires. These functions are defined as follows. If a vertical displacement g j ðtÞ ¼ e ixt is applied to the vehicle at jth tire, with all the other inputs being zero, then the response of the vehicle is given by ZðtÞ ij ¼ HðxÞ ij e ixt . Based on Eq. (68), it is straightforward to see that FRF matrix can be given by and where AR and AI; respectively, represent real and imaginary parts of matrix A, QR and QI; respectively, represent real and imaginary parts of matrix Q. It is convenient to express FRF matrix in terms of real part HR and imaginary part HI where HR ¼ AR Á QR À AI Á QI; From the stochastic process theory, PSD response of a linear system and PSD excitation satisfy the following relationship [76]: or equivalently, where S Z ðX ij ; xÞ ¼ ½S Z i Z j ðX ij ; xÞ MÂM ,S g ðX ij ; xÞ ¼ ½S g i g j ðX ij ; xÞ NÂN , H Ã ðxÞ and H T ðxÞ are, respectively, the conjugate matrix and the transposed matrix of FRF matrix HðxÞ. It may be noticed that H Ã ðxÞ and H T ðxÞ can be, respectively, expressed in the form of real and imaginary parts by means of Eqs. (72a) and (72b) Similarly, we may write S Z ðX ij ; xÞ as in which SR Z ðX ij ; xÞ and SI Z ðX ij ; xÞ are, respectively, the real and imaginary parts of response spectral density S Z ðxÞ. If the pavement excitation spectral density S g ðxÞ is a real spectrum matrix, then by expanding (72b) using (73a) and (73b), we have SR Z ðX ij ; xÞ ¼ HR Á S g ðX ij ; xÞ Á HR T þ HI Á S g ðX ij ; xÞ Á HI T : ð75aÞ and SI Z ðX ij ; xÞ ¼ HR Á S g ðX ij ; xÞ Á HI T À HI Á S g ðX ij ; xÞ Á HR T : If the pavement excitation spectral density S g ðxÞ is a complex spectrum matrix, say S g ðX ij ; xÞ ¼ SR g ðX ij ; xÞ þ iSI g ðX ij ; xÞ ð 76Þ then by expanding Eqs. (72b) using (73a), (73b), and (76), we have SR Z ðX ij ; xÞ ¼ ½HR Á SR g ðX ij ; xÞ þ HI Á SI g ðX ij ; xÞ Á HR T ; À ½HR Á SI g ðX ij ; xÞ À HI Á SR g ðX ij ; xÞ Á HI T ; and SI Z ðX ij ; xÞ ¼ ½HR Á SI g ðX ij ; xÞ À HI Á SR g ðX ij ; xÞ Á HR T þ ½HR Á SI g ðX ij ; xÞ À HI Á SR g ðX ij ; xÞ Á HI T : Figure 9 shows a sketch of the ith tire contacted with rough pavement surface. From Newton's second law of motion it is direct to know that the contact force between ith tire and pavement, P i ðtÞ, can be given by

Time-domain analysis
where k i and c i , respectively, the ith tire spring stiffness and viscous damping, both being constant. According to stochastic process theory, if the input of a linear timeinvariable system is a stationary random process, then the output of the system is also a stationary random process [76]. As was noted, since temporal random excitation of pavement surface Gaussian ergodic random process with zero mean, the response Z i (t) and its derivative process _ Z t ðtÞ are both zero mean Gaussian stationary stochastic processes. Taking expectation to both sides of (78), we get the mean function of dynamic contact forces, m P m P i ðtÞ ¼ E½P i ðtÞ ¼ k i E½Z i ðtÞ þ c i E½ _ Z i ðtÞ ¼ 0: It should be pointed out that the mean function m P (t) here, not containing the static load of vehicle, only represents the statistical average of dynamic effect. If the static load is considered, then the complete mean function of dynamic contact forces, P i , is given by An overview of a unified theory of dynamics where g is the acceleration due to gravity, P 0i and m i are, respectively, the effective weight and effective mass distributed by the vehicle on ith tire. Without special statements, the derivation of formulation in the context is based on (79).
It is of interest to know the correlation function of dynamic contact forces between vehicle and pavement because this information is critical for analyzing the dynamic response of pavement structure under vehicle loads. According to the definition of correlation function and using Eq. (78), we have It is known from stochastic differential principles that there exists exchangeability between expectation and square differential [31]. In other words, if a random process Z(t) is differentiable for any order, it is proved that o nþm ot n os m R ZZ ðt; sÞ ¼ R Z ðnÞ Z ðmÞ ðt; sÞ; where Z ðnÞ ¼ d n ZðtÞ=dt n . In addition, if ZðtÞ is a stationary process, then (73a) becomes ðÀ1Þ n o nþm os nþm R ZZ ðsÞ ¼ R Z ðnÞ Z ðmÞ ðsÞ: ð83Þ Applying Eq. (83) into Eq. (81) gives It is clearly seen that autocorrelation function is obtained as i ¼ j and crosscorrelation function is obtained as i 6 ¼ j.

Frequency-domain analysis
Having obtained correlation function, it is not difficult to represent the power spectral density (PSD), which is defined as the Fourier transform of correlation function. From Eq. (81) we have where S P i P j ðX ij ; xÞ ¼ 1 2p R P i P j ðX ij ; sÞe Àixt ds; ð86aÞ Correlation function R P i P j ðsÞ is the Fourier inverse transform of PSD function S Z i Z j ðxÞ: Taking derivative to both sides of Eq. (87) gives [31,182] d n ds n R Z i Z j ðX ij ; xÞ ¼ i n By combining Eqs. (83) and (88), we find that Substituting (88) into (84), we eventually get the PSD forces in terms of PSD of relative displacement response: Equation (89) applies to both direct-spectrum and crossspectrum. In addition, as far as direct-spectrum is concerned, we can rewrite PSD forces in more concise form: Standard deviation of contact forces can be expressed in terms of PSD forces. Applying Eq. (90) to the definition of standard deviation directly shows that Furthermore, if the standard deviations of displace response and velocity response, r Z i and r _ Z i , are deployed, Eq. (91) can be rewritten as in which

A Case Study of Walking-Beam Suspension System
As a case study, a walking-beam suspension system is shown in Fig. 10. Governing equations of this vehicle's suspension system are also given by (65) in matrix form, where parameters and their values related to this walkingbeam system are listed in Table 2 A two-dimensional isotropic random field is used for numerical simulation. The one-dimensional PSD roughness proposed by the International Organization for Standardization (ISO) in the power form [49,52] with parameters S 0 = 3.37 9 10 -6 m 3 /cycle and c = 2.0 is adopted for numerical computation. Figure 11 shows the PSD contact force of the right tire at the speed of 20 m/s. In aforementioned study, the contact area between the tire and the pavement surface is assumed as a point contact. There is no difficulty to extend this point contact to a distributed contact by considering the footprint as a weighted integration of contacting points [77][78][79]. It should be noted that the effect of nonlinearity in vehicle suspension and variable speed of travel (e.g., acceleration, deceleration, etc.), and inhomogeneity in pavement surface on contact force between vehicle and pavement have not been addressed here, which have been considered in various studies [80].

Background
A large class of time-dependent sources such as vehicle, submarines, aircraft, and explosion-induced waves belongs to moving sources. The study of response of media (pavement, runway, rail-track, bridge, air, and sea) to moving sources is named moving source problem (MSP), which is of particular interest to structural design, noise assessment, target detection, etc. [81]. A number of studies have been addressed to MSP in various fields of physics. For instance, the response of an ice plate of finite thickness caused by moving loads was discussed by Strathdee et al. [82]. Wells and Han [83] analyzed the noise generated by a moving source in a moving medium. As for the aspect of Fig. 9 A sketch of the ith tire on rough pavement surface An overview of a unified theory of dynamics 149 elastodynamics, even more investigations are being done [84][85][86][87][88][89][90][91][92][93][94]. One may refer to Sun [95] for detailed review of MSP. Methods for solving MSP primarily include integral transformation, characteristic curve, and modal analysis [96][97][98][99][100]. A common characteristic of previous studies is to utilize Galilean transform. The advantage of using Galilean transform is that the governing field equations, usually partial differential equations, can be reconstructed in a moving coordination so that the effect of source velocity may be reflected in parametric ordinary differential equations. However, since Galilean transform requires that the source is steadily moving, the methods based on Galilean transform only apply to the steady-state solution. As for transient solution, a feasible way is to directly apply high-order integral transformation to the field equations. Since velocity parameter is included in boundary or initial condition, there is intractable obstacle when integrating the field equations. This can be the reason why most of the MSPs are solved only for steady-state solution.
As far as pavement response is concerned, some pioneers have dedicated their research to the subject of effects of moving loads on pavement structure. To perform the dynamic analysis of pavement structures, one of the indispensable considerations is the dynamic traffic loading, which is the excitation source of pavement structures. Since pavement loads are caused by vehicle vibration, it is necessary to include the vehicle in the investigations of pavement loads. Employing MIT heavy truck simulation program Hedric and his associates, Abbo and Markow [9][10][11], examined the influence of joint spacing, step faulting, vehicle suspension characteristics, and vehicle velocity on pavement damage as defined by fatigue cracking in concrete slab. The response of pavement under their consideration is achieved using PMARP, a static finite element program for pavement analysis developed by the Purdue. The conclusions obtain provide a valuable approximation of pavement response.While PMARP is a static finite element program, the properties of inertia and damping of the pavement structure were not included in consideration. An advanced pavement model called MOVE program was developed by Chen [101] and Monismith et al. [73] at UC-Berkeley. This model can take the motion of load into account by means of finite element methods and is an important dedication to pavement structural analysis. The vehicle load is assumed to be an infinite line load to simplify the analysis from a three-dimensional problem to a two-dimensional problem. In addition, the model is based on the deterministic elastodynamics; therefore, neither pavement surface roughness nor vehicle suspension can be considered in the presented framework.
Since dynamic effects have been increasingly important in the prediction of pavement response, damage, and performance [14,16,18], it has become necessary to develop better mathematical models to account the effects of motion and fluctuation of contact forces caused by various types of vehicle suspensions [73]. The author has carried out a number of studies on deterministic and stochastic MSPs over the last two decades by providing a general approach that can include surface roughness and vehicle suspensions in the response of continuum media [54,74].
The complexity in MSP is the variable load position, which not only makes the representation of the field  Pavement structures are traditionally classified into two categories: rigid pavement and flexible pavement. The former often refers to Portland cement concrete (PCC) pavement, the latter often refers to asphalt concrete (AC) pavement. This classification is not very strict because PCC pavements are not ''rigid,'' or are AC pavements ''flexible'' in many situations. From the standpoint of structural nature PCC pavements behave more like a slab, while AC pavement acts more like a multi-layered system. Therefore, it is more reasonable to classify pavement on the base of mechanic properties. There have been developed many models to simulate the behavior of pavement structures. For instance, Hardy and Cebon [16,18] simplified AC pavement into an Euler-Bernoulli beam. From the theoretical perspective, most of the proposed pavement models belong to continuum media [102,[183][184][185][186][187]. x' x Consider a linear medium with the region R or the boundary B being at rest initially. A load is applied when the medium is at rest and moves according to a given law of motion (Fig. 12). The MSP is to solve the response of the medium to the moving load. The governing equation of such a system thus belongs to linear partial differential equation. According to the theory of linear equation, the solution of the equation can be constructed by the integration of the fundamental solution of the equation or the so-called Green's function. It is usually described as the superposition principle or, equivalently, the reciprocal principle [103].
Besides the above-mentioned description of the continuum media, the following assumptions are made here. One assumption is that, comparing to the mass of vehicle, the mass of pavement structures (including surface layer, base, subgrade, and soil foundation) is large enough such that pavement vibration is much smaller than vehicle vibration. Another assumption is that wave velocity excited by a dynamic load in pavement structure is much faster than vehicle's speed of travel. With these two assumptions, the couple effect of vehicle-pavement interaction can be negligible.

Moving stochastic vehicle loads
Contact forces applied on pavements by vehicles are moving stochastic loads. The statistical characteristics of the dynamic contact forces have been discussed in detail previously. Dynamic contact forces follow with two aspects of meanings. One is that the location of the force is changing continuously with traveling of vehicle, and another is that the amplitude of the force is varying due to vehicle vibrations.
Without any loss of generality, a vehicle is assumed to be a linear system traveling along the x-axis at a constant speed. The contact forces between the vehicle and pavement can be given by a concentrated moving load: where t and v are time variable and vehicle velocity, respectively, d() is the Dirac-delta function, which is defined by In addition, P sto ðtÞ, a function of time, represents amplitudes of stochastic contact forces and can be expressed as two terms below: Here P 0 , representing the static load applied at a tire when the vehicle is at rest, is a constant quantity, and P(t) represents the dynamic portion of the stationary stochastic contact forces with a zero mean. Its correlation function, power spectral density and standard deviation are, respectively, denoted as R PP ðsÞ, S PP ðxÞ, and r P .

Representation theory of MSP
It is convenient to assume a three-dimensional configuration with observation variable x = (x,y,z), source variable n ¼ ðn; g; fÞ, and time t ! 0. Suppose a linear differential operator O describes the dynamic property of a physical system and appropriate interface and boundary conditions relate the field quantities of specified problems. Obviously, for a concrete elastodynamic problem, the linear differential operator O is given by the well-known Navier-Stoke's field equations. The Green's function is then defined as the fundamental solution of the system. In other words, for the problem discussed here, the Green's function corresponds to the solution of the governing equations as the point source takes the form of a Dirac delta function in both spatial and temporal domains.

General Formulation
Without loss of generality, vanishing initial conditions are assumed here. According to the causality of a physical system, the Green's function Gðx; tÞ ¼ 0 for t\0. We may then write O½Gðx À n; t À sÞ ¼ dðx À nÞdðt À sÞ: Here, the initiation of the source is delayed by s. The causality of a physical system requires that for Green's function t s: Since the initial condition of the linear medium is zero, the dynamic response could be expressed as where x ¼ ðx; y; zÞ, n ¼ ðn; g; fÞ, v ¼ ðv x ; v y ; v z Þ, dn ¼ dndgdf, and S is the region R or the boundary B. It is know that the Green's function Gðx; t; nÞ corresponds to the solution of the equation when a unit point impulse is applied at position n. Assume that the medium is infinite in those dimensions of interest. It is known that the governing equations are linear. According to the reciprocal principle, the response of media at field point x in the fixed coordinates when the source lies at n in the fixed coordinates is equal to the response of media at field point x À n in the fixed coordinates when the source lies at 0 in the same coordinate.
Define the impulse response function (IRF) hðx; tÞ as the solution of O½hðx; t À sÞ ¼ dðxÞdðt À sÞ: According to the above mentioned analysis, we have hðx À n; tÞ ¼ Gðx; t; nÞ ¼ Gðx À n; t; 0Þ: Substituting (101) into (99), we get Furthermore, if the concrete form Fðn; sÞ in (103) is considered, we may rewrite (102) as in which vs ¼ ðv x s; v y s; v z sÞ. Also, the property of Dirac delta function (102) is used here. If the transformation h ¼ t À s is used, then (103) can be expressed as Equations (103) and (104) are general results for MSP and are named generalized Duhamel's integral (GDI). From the point of view of time history, we might regard a moving load as a series of impact on continuum media during a number of tiny time intervals. The integration of the response of the medium excited by each impulse is thus equal to the cumulative effect of the moving load. Although solving for the Green's function is still a nontrivial task, convolution (103) does provide a sound theoretical representation, which can be very powerful when combining with numerical computation such as finite element method.

Deterministic analysis for a moving constant load
The solution of the problem described here can be constituted using the solutions of two individual problems because of the superposition principle of the solution of linear equations. One problem deals with the deterministic response of the medium under the moving constant load P 0 dðx À vtÞdðyÞ. The other problem deals with the random response of the medium under the moving stochastic load PðtÞdðx À vtÞdðyÞ. In this section, the first problem is analyzed. In the next section, the second problem is analyzed. The summary of the solution is provided in the sections followed. If the load is a constant with amplitude P 0 , the response are given by It is obvious that the response is no longer a constant independent on the time t. If the source is a fixed load, according to our example, the solution should be a static quantity. Actually, this idea is straightforward demonstrated by putting the velocity variable v = 0 into (105). The result gives that uðx; tÞ ¼ P 0 Clearly, the response is a constant without depending on the time variable.

Stochastic analysis for a moving stochastic vehicle load
In this section, a stochastic moving source is analyzed. In the derivation of the GDI in (103), we require no special assumptions on P sto ðtÞ. Therefore, if P sto ðtÞ is a stochastic process, and (103) becomes an integral in the sense of Stieltjes integration. Meanwhile, the response uðx; tÞ becomes a stochastic process. We now consider the response of media to PðtÞdðx À vtÞdðyÞ.
As mentioned before, P(t) is a zero mean stationary process with autocorrelation function R PP ðsÞ, PSD S PP ðxÞ and standard deviation r P .Taking the expectation of both sides of (103) and using the exchangeability of expectation and integration, we obtain the mean function of the response, i.e., uðx; tÞ ¼ It is not difficult to obtain the spatial-time correlation functions for the response, i.e.,

À1
E½Pðs 1 ÞPðs 2 Þhðx 1 À vs 1 ; y 1 ; z 1 ; t 1 À s 1 Þ Â hðx 2 À vs 2 ; y 2 ; z 2 ; t 2 À s 2 Þds 1 ds 2 ; ð108Þ where R u ðx 1 ; x 2 ; t 1 ; t 2 Þ is the correlation function of response uðx; tÞ. Let x 1 ¼ x 2 ¼ x; then we obtain the time autocorrelation function where h j ¼ t j À s j ðj ¼ 1; 2Þ. By substituting t 1 ¼ t 2 ¼ t into (108) and (109), it is straightforward to find second moment functions, i.e., the mean square functions of the random response. It has been known that for a linear system with a stationary stochastic excitation at a fixed position, the response of that system is still a stationary stochastic process [76]. However, this conclusion only applies to the fixed source problem. For a linear system with a moving stochastic source, the random response of that system is a nonstationary stochastic process even if the random excitation PðtÞ is a stationary stochastic process [54,74]. To show this, we inspect (107) and (109), which are apparently not stationary because time variable t is contained in the kernel function of the solution. In other words, in a situation where a source is moving with respect to a receiver a nonstationary signal will be recorded at the observer position, even when the source produces a stationary output. This effect is known as the Doppler shift. It is also useful to realize that although the system is a linear system, it essentially becomes a time varying system when a moving source is applied.
In many circumstances, there may be a demand to push the analysis further into the frequency domain. For instance, response information of amplitude distributions and frequency components of media to moving vehicle loads is needed for vehicle optimum control and pavement performance prediction. To fulfill this purpose, it is necessary to use spectral analysis techniques to obtain the required information. As such, one may encounter difficulties when performing Fourier spectral analysis technique because this technique has been devised primarily for stationary signals. Although some variations of Fourier spectral analysis technique have been introduced to dealing with nonstationary stochastic processes, Fourier spectral analysis are not ideal tools for nonstationary signals induced by MSP. This shortcoming has led Sun [104]to develop the so-called follow-up spectral analysis, by which the commonly used spectral analysis technique is still applicable with sound theoretical foundation.
Let coordinates oxyz and o 0 x 0 y 0 z 0 be, respectively, fixed coordinates and follow-up coordinates moving with the moving source. The relationship between the two coordinates is Therefore, a moving source x 0 þ vt in the fixed coordinate oxyz becomes a fixed source x 0 ¼ x 0 in the follow-up coordinate o 0 x 0 y 0 z 0 . Here x 0 is a constant vector. Now consider the response of the medium at a moving field point x þ vt þ x 0 in the fixed coordinates. Utilizing (104), we have uðx þ vt þ x 0 ; tÞ ¼ The mean function and autocorrelation function of the response described by (111) are, respectively, given by and R u ðx þ vt þ x 0 ; sÞ According to the definition of a stationary stochastic process [31], (112), and (113) indicate that the response at a moving field point x þ vt þ x 0 in the fixed coordinates becomes a stationary process. It is evident that the moving field point x þ vt þ x 0 in the fixed coordinates becomes a fixed field point x þ x 0 in the follow-up coordinates as illustrated in Fig. 12. To show this, we replace x in (110) In other words, the response of a fixed position x ? x 0 in the follow-up coordinates to a moving stationary stochastic load possesses the stationary property. Thus, Fourier spectral analysis technique is still applicable here. It should be pointed out that the explanation of the stationary process uðx þ vt þ x 0 ; tÞ is essentially different from the commonly described stationary process. In general, the mean function of a commonly described stationary process refers to the time average of the random process, while the mean in (112) is indeed interpreted as a spatial average of the random response. The same explanation applies to the autocorrelation function shown in (113).
In light of the aforementioned explanation, there is no difficulty to define PSD in the follow-up coordinates. The following spectral analysis is performed in the follow-up coordinates. We may rewrite the autocorrelation function in the follow-up coordinates as R u ðx 0 ; sÞ ¼ where x 0 is expressed by (114). Let s ¼ 0 and put it into Eq. (115), we obtain the mean square function of the random response in the follow-up coordinate In addition, since PðtÞ is a zero mean stationary process, we have where r p and Var½PðtÞ are, respectively, the standard deviation and variance of PðtÞ. So we can rewrite (117) as where r u and w 2 u are the standard deviation and variance of the response field, respectively. Define the relationship between the frequency response function and the impulse unit response function in the follow-up coordinates as According to Wiener-Khintchine theory, the PSD and the autocorrelation function form a Fourier transform pair. Taking Fourier transform to both sides of (115) and noticing the (119), we obtain the expression of PSD in follow-up coordinates, i.e., Similarly, an expression for the time autocorrelation function can be obtained by taking the Fourier inverse transform of (120) Hence, the mean square function is also given by

Pavement Models
There are mainly two types of pavement structures: Portlant cement concrete pavement and asphalt concrete pavement (Fig. 13). These pavement structures can be modeled by a beam, a slab, a layered medium on a halfspace or rigid bedrock. A number of studies have been carried out lately by Sun and his associates using analytic method and analytic-numerical method . Beskou and Theodorakopoulos [164] provided a recent review on numerical methods for studying dynamic effects of moving loads on road pavements. Sun and Greenberg [74], Sun and Luo [125][126][127][128], Sun et al. [153], Luo et al. [165] and Sun et al. [123] provide concrete examples of pavement models subject to moving loads. The response of pavement systems under dynamic loads may be expressed in partial-differential equations. A generic description of the governing equations is: where u½Á is a partial-differential operator, x ¼ ðx; y; zÞ a spatial vector in Cartesian coordinates, t is the time variable, Pðx; tÞ is the applied dynamic load (i.e., the input), which can be recorded by data acquisition system during laboratory experiments and field tests, uðx; t; hÞ the pavement response vector (i.e., the output in the form of displacements, stresses, and strains), and h ¼ ðh 1 ; h 2 ; . . .; h n Þ is the parameter vector to be identified. A pavement structure usually consists of a surface course, base courses and subgrade. Within each course, there may be several sub-layers made up of different materials. A two-dimensional Kirchhoff thin slab resting on a Winkler foundation is the common model for PCC pavements. The operator u½Á for a Kirchhoff thin slab is where D ¼ Eh 3 =½12ð1 À l 2 Þ, and q, l and h are the density, Poisson's ratio and thickness of the slab, E is the Young's elastic modulus, K is the modulus of subgrade reaction, C is the radiation damping coefficient, and r 2 ¼ o 2 =ox 2 þ o 2 =oy 2 is the Laplace operator. The parameter vector is h ¼ ðE; h; l; K; C; qÞ. For a multilayered flexible pavement system, the governing equation is controlled by a three-dimensional Navier-Stokes's equation for each layer where r ¼ o=ox þ o=oy þ o=oz, f the body force vector, the Laplace operator , g d is the hysteretic damping coefficient, Lamb constants k and G the bulk modulus and shear modulus, respectively, k Ã and G Ã the complex counterparts of k and G, respectively. The subgrade may be artificially divided into a number of thin layers. Within each layer the soil is characterized to be isotropic, homogenous and have the same structural and material properties, while these properties vary for different layers. Furthermore, physical nonlinearity may possibly be presented in asphalt surface layer and soil subgrade using nonlinear constitutive models involving viscoelasticity-viscoplasticity. Eq. (125) adopts the simplest model to account for viscoelasticity. A more generic model is the generalized viscoelastic model, which includes the Burgers model, Maxwell model and Kelvin model as its special cases. Fig. 14 presents a schematic plot of a list of viscoelastic models.
To solve a viscoelastic problem, elastic solutions is sought first and then the correspondence principle is applied to convert the elastic solution into a viscoelastic solution. Two elastic/viscoelastic subgrade models will be studied: a half-space and a layer resting on bedrock. Clearly, the parameter vector h ¼ ðE; h; k; G; qÞ varies from layer to layer. When viscoelasticity is considered, more parameters will appear in the parameter vector. In principle, the adoption of a generalized viscoelastic model in forward dynamic analysis introduces no significant difficulty. The number of Kelvin components in this model can be estimated through a thorough investigation. With properly specified initial and boundary conditions, Eqs. (123)-(125) constitute a complete mathematical description of the forward dynamic problem. Forward analysis aims to solve for the response uðx; t; hÞ provided that the excitation Pðx; tÞ and the parameter vector h are known. These mathematical, physical models describe the behavior of different types of pavement systems and they will be studied in great depth.
Equation (123) belongs to a wave equation from a mathematical physics point of view. Its solution can be obtained in the form of a Lebesgue-Stieltjes integral using proper integral transformation, depending upon the nature of the problem (e.g., steady-state vs. transient) and upon how the problem is formulated (e.g., in Cartesian or cylindrical coordinates). Let the Green's function (the fundamental solution) of Eq. (1) be Gðx; t; hÞ ¼ u À1 ½dðx; tÞ, in which dðÁÞ is the Dirac-delta function and u À1 ½Á denotes the inverse operator of u½Á. Let pavements be at rest prior to the NDE test, leading to a vanishing initial condition. The solution of Eq. (123) under loading condition Pðx; tÞ can be constructed as where S is the region where Pðx; tÞ is defined. Equation (126) can also be equivalently represented in the transformed domaiñ where TfÁg is a transformation operator,ũðn; x; hÞ and Gðn; x; hÞ are the response and the Green's function in the transformed domain, respectively,S is the region in the transformed domain where the transformed dynamic load Pðn; xÞ is defined; n ¼ ðn; g; fÞ and x are the counterpart of spatial vector x ¼ ðx; y; zÞ and time variable t. The Thomas-Haskell method relates a transformed response at the bottom of a layer, in the form of a transfer matrix, to a corresponding quantity at the top of a lower layer. For the last half century, this method has served as the cornerstone for numerous studies in multilayered elastic analysis. Another benchmark proposed by Kausel and Roesset contributes an alternative, in which the dynamic stiffness matrix is expanded in terms of wave number and approximated by taking terms only up to the second order of the wave number. The Thomas-Haskell method is more accurate than the Kausel-Roesset method but demands more computational effort. Forward dynamic analysis for multilayered viscoelastic media will not add significant difficulty because both methods are still applicable, though it is a demanding task. Since each method has its own strengths and weaknesses, both methods will be adopted in the project to tackle wave propagation through multilayered viscoelastic media. Equations (126) and (127) are amenable to numerical evaluation of the dynamic response in the timespace domain and in the transformed domain, respectively, provided that Pðx; tÞ and h are known. Without doubt, the computation here involves intensive numerical evaluation of multifold integration of complex functions with unstable characteristics in time and space.

Discussion and Future Research
In this section, discussion and future research of the unified theory of dynamics of vehicle-pavement interaction under moving and stochastic load is carried out from the following aspects: (a) nature of the problem, (b) modeling, (c) methodology, and (d) further extension and engineering application.
(a) Nature of the problem. The setting of vehiclepavement interaction can be categorized either in a deterministic framework or a stochastic framework. Stochasticity may be presented in speed, magnitude and position of the loading, structural models of vehicle-pavement system and constitutive models. A stochastic framework provides a more realistic setting but exhibit more complexity. The setting of vehicle-pavement interaction can also be categorized either as a steady-state problem or as a transient problem. Except very few studies [104,145], almost all existing literatures belong to a steady-state problem described in a deterministic framework. The modeling of vehicle-pavement interaction involves the load model, the vehicle model and the pavement model. The load model addresses spatial distribution of the load (e.g., concentrated load, distributed load, multiple load, etc.), speed of the load (e.g., constant speed, varying speed, etc.), trajectory of the load (e.g., straight line, curve, etc.), and magnitude of the load (e.g., constant load, impact load, sinusoidal load, varying load, and random load). The vehicle model addresses vehicle suspension system, which involves mass distribution, dimension, and configuration of the vehicle. Quartervehicle model, half-vehicle model, and whole vehicle model have all been established with increasing complexity involving different number of spring and dashpot elements. Regardless of how many spring and dashpot elements are being used in the built vehicle model, almost all the literature studies tend to use linear suspension system due to its ease of computation when integrated with a pavement model. The pavement model addresses the simplification of structural system of pavement using beam model, slab model, and layered medium on a half space or bedrock. In terms of constitutive models of the material, almost all studies only considered linear elastic and linear viscoelastic materials when tackling dynamics of vehicle-pavement interaction because of the complexity of the problem. Asphalt concrete pavement actually exhibits complex viscoelastic-viscoplastic-damage properties [166][167][168][169][170][171][172], which should be integrated into the subject. Also, it is rare to consider the coupling effect of vehicle-pavement interaction, while this is not the case in train-railway interaction because of the mechanism of the interaction and the relative mass difference between vehicles (e.g., car, truck, train, and airplane) and transportation infrastructure (e.g., highway, bridge, and railway). None of the studies has addressed vehicle-pavement interaction using deteriorated pavement models as well as long-term pavement damage and failure due to vehicle-pavement interaction.
(c) Methodology. The methods for solving vehiclepavement interaction problem can be classified into analytic approach and numerical application, though the implementation of analytic approach still requires numerical computation. For vehicle dynamics, equation of motion is first established as a set of linear differential equation system and solved in the frequency domain using frequency response function or in the time domain using numerical sequential integration. For pavement dynamics, analytic approach makes use of integral transform to treat wave propagation in continuum media, while numerical approach such as finite difference, finite element and boundary element methods makes use of discretization in time and in space. The advantage of analytic approach is that it provides insights for revealing physics of the wave propagation in continuum media and can be highly efficient in terms of computation implementation, particularly when the spatial scale of the problem involves hundreds of kilometers (e.g., seismology) or the speed of the load is very high and close to various critical speeds of waves in media. The advantage of numerical approach is that they can deal with pavement having complex geometric structure as well as nonlinear constitutive model of the material.
(d) Further extension and engineering application. The study of dynamics of vehicle-pavement interaction provides a deep understanding for improving vehicle design (e.g., road-friendly vehicle suspension system), road transportation safety [173][174][175][176][177], long-lasting pavement structures design [178][179][180][181], ride quality and infrastructure asset management. An improved quantitative understanding on effects and mechanisms of various factors on dynamics of vehicle-pavement interaction is the fundamentals for increased application and accuracy and reliability of structural health monitoring, nondestructive testing and evaluation, environmental vibration mitigation and weight-in-motion. It will also benefit the nation's transportation economy by reducing operation and maintenance costs of vehicles and transportation infrastructure as well as increasing transportation productivities.

Conclusions
Irregularities of pavement surface, from the small-scale unevenness of material on pavement surface to the largescale undulating of vertical curve of a highway or an airport, all belong to spatial fluctuation of pavement surface at different scale. Extensive study has been accomplished both domestically and internationally, toward measuring the physical aspects of pavement roughness, analyzing the resulting data, and evaluating the riding performance of pavements. Instrumentation and analysis technique including the spectral analysis approaches have been summarized in the article. A number of the PSD functions including the effect of thermal joints on PSD roughness have been presented here. These PSD functions are similar in their shapes and only different in their mathematical descriptions. Under the condition of constant speed of travel and linear vehicle suspension, it is proven that dynamic contact force between vehicle and pavement is a stationary stochastic process. Its mean function is given by (79) without the consideration of static loads or by (80) with the consideration of static loads. The correlation function, PSD forces, and standard deviation are, respectively, given by (84), (89) and (92). The concept and the methodology present in the article are not restricted to specific surface roughness and/or vehicle models. They are generally applicable to all kinds of linear vehicle models and measured pavement surface conditions. The response of linear continuum media to moving stochastic vehicle loads is analyzed herein. We show that there exists predicable relation among surface roughness, vehicle suspensions and speed, and the response of continuum media. The theory developed here is widely applicable to moving vehicle loads.