*Article* **Classification of Variable Foundation Properties Based on Vehicle–Pavement–Foundation Interaction Dynamics**

#### **Robin Eunju Kim**

Department of Civil and Environmental Engineering, Hanyang University, Seoul 04763, Korea; robinekim@hanyang.ac.kr; Tel.: +82-2220-0413

Received: 15 October 2020; Accepted: 31 October 2020; Published: 3 November 2020

**Abstract:** The dynamic interaction between vehicle, roughness, and foundation is a fundamental problem in road management and also a complex problem, with their coupled and nonlinear behavior. Thus, in this study, the vehicle–pavement–foundation interaction model was formulated to incorporate the mass inertia of the vehicle, stochastic roughness, and non-uniform and deformable foundation. Herein, a quarter-car model was considered, a filtered white noise model was formulated to represent the road roughness, and a two-layered foundation was employed to simulate the road structure. To represent the non-uniform foundation, stiffness and damping coefficients were assumed to vary either in a linear or in a quadratic manner. Subsequently, an augmented state-space representation was formulated for the entire system. The time-varying equation governing the covariance of the response was solved to examine the vehicle response, subject to various foundation properties. Finally, a linear discriminant analysis method was employed for classifying the foundation types. The performance of the classifier was validated by test sets, which contained 100 cases for each foundation type. The results showed an accuracy of over 90%, indicating that the machine learning-based classification of the foundation had the potential of using vehicle responses in road managements.

**Keywords:** machine learning-based classification; non-uniform foundation; stochastic analysis; vehicle–pavement–foundation interaction

#### **1. Introduction**

Road infrastructure forms a basic component in transportation, providing connectivity between local, regional and global value chains. Despite the impacts of road's serviceability on the economy and public safety, maintenance is inadequate, due to its extensive nature. For example, the Federal Highway Administration reported that 26 percent of major urban roads in the U.S. are in a poor condition [1]. FHWA also reported that a capital of 182 billion dollars was spent in 2008 on improvements and maintenance of federal highway, while they are still in shortage [1]. Thus, research that develops tools and methods for assessing road conditions assume greater importance.

Typically, a pavement's condition is assessed by measured information, such as ride comfort, surface defects, and structural adequacy [2]. For example, the Pavement Condition Index (PCI), developed by the U.S. Army Corps of Engineers rates the surface operational conditions including rutting, potholes, crackings, etc. [3]. Recent advances in image processing and deep learning technologies demonstrated that an automatic rating of PCI using a visual platform is available. For example, the automated pavement management system equipped with visual inspecting tools were developed to evaluate the pavement deteriorations and cracking [4,5]. Multimetric sensors including wireless sensing modules are utilized as well to examine the surface condition [6,7], fatigue [8], and certain anomalies on road structures [9]. However, unlike visible surface defects, structural adequacy is associated with the load transfer capability of the subgrade layers.

The falling weight deflectometer (FWD) is an essential non-destructive tool that is widely used for evaluating the structural adequacy of the pavement [10]. In FWD, as the weight falls on the pavement to be examined at some height, the response of the pavement, including deflection, is measured. Then, the responses are related to the strain and elasticity of the pavement, to examine the adequacy of the sublayer [11]. However, due to complexities in the testing method, which is usually performed during post-processing of the collected data, and due to required technical expertise, significant time, and costs, the network-level application is limited [12]. With the advancement of the wireless module's sensing and calculating capability, the inspection and monitoring fields are in the transition from human-oriented inspection to machine-based inspection [13]. The following literature shows some successful examples of monitoring foundation noise excitation [14], decentralized road networks [15], which are known to be complicated, compared to other applications. Thus, so far, predicting the capability of road structure with rather portable and automated devices are of interest, but a challenging task, due to its complicated mechanisms.

Within a road structure, the main excitation source is a moving vehicle. To understand the responses of the moving loads, various foundation types were examined, based on analytical models. One of the simplest model was developed to examine vehicle response due to road roughness on a non-deformable foundation. Roughness was first modeled as Gaussian random signals [13]. Then, the models improved to contain a more realistic input, such as a stationary zero-mean process with a certain power spectral density (PSD) [16–18]. Among various research works, Wedig derived a closed-form expression of the covariance response of a vehicle model, by integrating the PSD of road roughness [19]. These models examined the impact of road roughness on vehicle responses, while the interactions due to pavement deflection were neglected.

To consider the deformable foundation, an Euler-Bernoulli beam resting on the viscoelastic foundation was investigated. Hardy and Cebon (1993) developed a quarter-car model on a smooth beam on a uniform Winkler foundation, to examine the vehicle–pavement–foundation interaction [20]. Similar approaches were adopted by other researchers and they used the models to understand the impact of vehicle parameters (including speed) on foundation responses [18,21]. Kelvin foundation under the Bernoulli beam was also adopted by authors in [22]. In their model, the interaction responses were examined by coupling the solutions of two systems—(1) vehicle on rough road and (2) elastic foundation subject to a single load. To eliminate the boundary condition effects, the frequency domain analysis of the interaction problem was performed on an infinite length beam [23]. Instead of handling infinite length, Kim et al. (2019) formulated the interaction system, based on a moving coordinate system, and examined the second-order stationary response of the interaction problem [24]. In the aforementioned studies, the foundation properties such as stiffness and density were assumed to be uniform, while in reality, those quantities might vary along the length of the road.

The non-uniform foundation on a beam was investigated by several groups of researchers. Early efforts focused on formulating a closed-form equation for varying foundation modulus, targeting statistical analyses. The linearly varying solutions were presented by Franklin and Scott (1979) [25] and higher-order variations were solved by the authors in [26,27]. The free vibration of the beams on the non-uniform foundation was studied by the following authors [28–30]. The authors in [30] compared the impact of nonlinear foundation on the deflection shapes and natural frequencies of the beam. Then, dynamic responses of a beam on the variable Winkler foundation, subject to a moving load, were studied by [31–33], and a moving mass was investigated by [34]. Although previous studies captured the effect of variable foundation on the pavement system, due to computational complexities, studies mostly neglected the inertial force effect from the moving vehicles.

In this study, the impact of the non-uniform foundation on vehicle responses was solved by developing the vehicle–pavement–(non-uniform)-foundation interaction model. In the model, the vehicle was represented with a moving-oscillator (a quarter-car). The pavement roughness was described with a filtered white noise model. The rigid foundation was modeled to have a finite-length Euler–Bernoulli beam on a deformable foundation. The top layer was modeled using the assumed

modes method. The subgrade was modeled with a Winkler-type foundation, in which the stiffness and damping properties varied along the length. The interaction model was then formulated in an augmented state-space representation. To effectively examine the response of the vehicle, the covariance of the response was then solved for the time-varying Lyapunov equation. Then, the equations were solved for various pavement roughness and foundation cases, to construct the covariance responses. Based on the estimated responses, six features that could distinguish the foundation types were selected and employed on a classifier. Subsequently, noise-added responses were employed on a linear classifier and demonstrated that the measured dynamics of a vehicle due to interaction could distinguish the foundation types and variations with an accuracy of over 90%.

#### **2. Model Formulation**

#### *2.1. Overview*

The vehicle–pavement–foundation interaction model considered herein is shown in Figure 1. The rigid pavement was modeled with an Euler–Bernoulli beam that had constant material properties. Elastic modulus (E), the moment of inertia (I), thickness (*tb*), cross-sectional area (*A*), density (ρ), and length (*L*). The vertical displacement of the beam due to interaction was defined as *u*B(*x*, *t*). The elastic foundation was taken as a Winkler-type foundation with varying stiffness (*k*f(*x*)) and viscous damping (*c*f(*x*)), along the length. The roughness of the pavement was modeled as a profile ξ(*x*) and superimposed on top of the beam.

**Figure 1.** Vehicle–pavement–foundation system.

The vehicle was represented with a quarter-car model, consisting of sprung mass (*m*s) and unsprung mass (*m*u). Their vertical movement was defined as *u*<sup>s</sup> and *u*u, while spring stiffness and damping coefficients at suspension and tire were denoted with *k*s, *k*t, *c*s, and *c*t, respectively. The vehicle was assumed to have constant velocity (*V*), as it traveled along the length. In the subsequent sections, a model formulation of the interaction system was introduced. Note that some derivations such as a state-space representation of the roughness were briefly discussed, while more detailed formulations could be found in the related literature [24,35,36].

#### *2.2. Basic Equations*

Employing the assumed modes method, the vertical deflection of the Euler–Bernoulli beam *u*B(*x*, *t*) could be defined in a series of sine functions, assuming a simple support boundary condition:

$$\phi(n) = \sin(\frac{n\pi(\mathbf{x} + L)}{L}) \text{ } n = 1 \text{ } 2 \text{ } 3 \text{ } \dots N \text{ } \tag{1}$$

where *N* is the total number of modes in the shape function. Then, the deflection of the beam could be rewritten as *u*B(*x*, *t*) = **N**B(*x*)*q*B(*t*). **N**B(*x*) is a mode shape vector containing defined mode shapes (φ(*n*)) and *qB*(*t*) is a time-dependent generalized displacement of the beam. The relationships for the first- and second-time derivative are . *u*B(*x*, *t*) = **N**B(*x*) . *<sup>q</sup>*B(*t*) and .. *u*B(η, *t*) = **N**B(*x*) .. *q*B(*t*).

Then, defining the displacement vector as *<sup>x</sup>*c(*t*) <sup>=</sup> *q*B(*t*) *uu zs* T , the equations of motion for the vehicle–pavement–foundation interaction system could be formulated as follows:

$$\mathbf{M}\_{\mathbb{C}}\ddot{\mathbf{x}}\_{\mathbb{C}}(t) + \mathbf{C}\_{\mathbb{C}}\dot{\mathbf{x}}\_{\mathbb{C}}(t) + \mathbf{K}\_{\mathbb{C}}\mathbf{x}\_{\mathbb{C}}(t) = -\mathbf{P}\_{\mathbb{S}} + \mathbf{P}\_{\mathbb{C}}\xi(t) \tag{2}$$

where

$$\mathbf{M}\_{\mathbb{C}} = \begin{bmatrix} \begin{bmatrix} \int\_{0}^{L} \rho A \mathbf{N}\_{\mathbb{B}}^{T} \mathbf{N}\_{\mathbb{B}} dx \\\ \begin{bmatrix} \mathbf{0} \end{bmatrix}\_{1 \times N} & m\_{s} + m\_{u} & m\_{s} \\\ \begin{bmatrix} \mathbf{0} \end{bmatrix}\_{1 \times N} & m\_{s} & m\_{s} \end{bmatrix} \end{bmatrix} \tag{3}$$

$$\mathbf{C}\_{\mathsf{C}} = \begin{bmatrix} c\_{\mathsf{f}}(Vt) \int\_{0}^{L} \mathbf{N}\_{\mathsf{B}}^{\mathsf{T}} \mathbf{N}\_{\mathsf{B}} d\mathbf{x} + c\_{\mathsf{t}} \mathbf{N}\_{\mathsf{B}}^{\mathsf{T}}(Vt) \mathbf{N}\_{\mathsf{B}}(Vt) & -c\_{\mathsf{t}} \mathbf{N}\_{\mathsf{B}}^{\mathsf{T}}(Vt) & [\mathbf{0}]\_{\mathsf{N} \times 1} \\ -c\_{\mathsf{t}} \mathbf{N}\_{\mathsf{B}}(Vt) & c\_{\mathsf{t}} & 0 \\ [\mathbf{0}]\_{1 \times N} & 0 & c\_{\mathsf{s}} \end{bmatrix} \tag{4}$$

$$\mathbf{K}\_{\mathbf{c}} = \begin{bmatrix} \int\_{0}^{L} \mathbf{E} \mathbf{N}\_{\mathbf{B}}^{\mathrm{r}T} \mathbf{N}\_{\mathbf{B}}^{\mathrm{r}} d\mathbf{x} + k\_{\mathrm{f}}(Vt) \int\_{0}^{L} \mathbf{N}\_{\mathbf{B}}^{\mathrm{T}} \mathbf{N}\_{\mathbf{B}} d\mathbf{x} + k\_{\mathrm{i}} \mathbf{N}\_{\mathbf{B}}^{\mathrm{T}}(Vt) \mathbf{N}\_{\mathbf{B}}(Vt) & -k\_{\mathrm{i}} \mathbf{N}\_{\mathbf{B}}^{\mathrm{T}}(Vt) & [\mathbf{0}]\_{\mathrm{N} \times 1} \\ -k\_{\mathrm{i}} \mathbf{N}\_{\mathbf{B}}(Vt) & k\_{\mathrm{i}} & 0 \\ [\mathbf{0}]\_{1 \times \mathcal{N}} & 0 & k\_{\mathrm{s}} \end{bmatrix} \tag{5}$$

$$P\_{\mathfrak{F}} = \begin{bmatrix} [\mathfrak{0}]\_{\mathbf{N}\times\mathbf{1}\prime} & -(m\_{\mathfrak{s}} + m\_{\mathfrak{u}})\_{\mathfrak{F}\prime} & -m\_{\mathfrak{s}}\mathfrak{g} \end{bmatrix}^{\mathrm{T}} \tag{6}$$

$$\mathbf{P}\_{\mathbf{c}} = \begin{bmatrix} -k\_{\mathrm{t}} \mathbf{N}\_{\mathrm{B}}^{\mathrm{T}}(Vt) & -c\_{\mathrm{t}} \mathbf{N}\_{\mathrm{B}}^{\mathrm{T}}(Vt) \\\ k\_{\mathrm{t}} & c\_{\mathrm{t}} \\\ 0 & 0 \end{bmatrix} \tag{7}$$

$$\xi(\mathbf{t}) = \begin{bmatrix} \ \xi(\mathbf{t}) & \ \dot{\xi}(\mathbf{t}) \ \end{bmatrix}^{\mathrm{T}} \tag{8}$$

Note that **N**B(*x*) is short noted as **N**B, except for the case when evaluated at *x* = *Vt*. Additionally, ξˆ(t) = ξ(*x*) *<sup>x</sup>*=*Vt <sup>k</sup>*f(*Vt*) = *<sup>k</sup>*f(*x*) *<sup>x</sup>*=*Vt*, *<sup>c</sup>*f(*Vt*) = *<sup>c</sup>*f(*x*) *<sup>x</sup>*=*Vt* and *g* is the gravity term.

Then, defining a state vector *<sup>x</sup>*<sup>T</sup> <sup>=</sup> *xc* . *xc <sup>T</sup>* and organizing Equation (2) in a state-space representation, yields:

$$\dot{\mathbf{x}}\_{\Gamma}(t) = \mathbf{A}\_{\Gamma}\mathbf{x}\_{\Gamma}(t) + \mathbf{B}\_{\Gamma}\boldsymbol{\xi}(t) = \begin{bmatrix} [\mathbf{0}]\_{N\gamma\times N\_{\Gamma}} & [\mathbf{I}]\_{N\gamma\times N\_{\Gamma}} \\ -\mathbf{M}\_{\mathbf{c}}^{-1}\mathbf{K}\_{\mathbf{c}} & -\mathbf{M}\_{\mathbf{c}}^{-1}\mathbf{C}\_{\mathbf{c}} \end{bmatrix} \mathbf{x}\_{\Gamma}(t) + \begin{bmatrix} [\mathbf{0}]\_{N\gamma\times 2} \\ \mathbf{M}\_{\mathbf{c}}^{-1}\mathbf{P}\_{\mathbf{c}} \end{bmatrix} \dot{\boldsymbol{\xi}}(t) \tag{9}$$

where *NT* is *N* + 2. The output vector *y*T(*t*) could be defined to contain arbitrary information about the system. In this study, the vehicle responses including displacement and velocity of the unsprung and sprung masses were considered as the output, i.e., *<sup>y</sup>*T(*t*) <sup>=</sup> *uu zs* . *uu* . *zs* . Then, the observation and feedthrough matrices yield:

$$\mathbf{y}\_{\Gamma}(t) = \mathbf{C}\_{\Gamma}\mathbf{x}\_{\Gamma}(t) + \mathbf{D}\_{\Gamma}\boldsymbol{\xi}(t) = \begin{bmatrix} [\mathbf{0}]\_{1 \times N} & 1 & 0 & [\mathbf{0}]\_{1 \times N} & 1 & 0 \\ [\mathbf{0}]\_{1 \times N} & 0 & 1 & [\mathbf{0}]\_{1 \times N} & 0 & 1 \end{bmatrix} \mathbf{x}\_{\Gamma}(t) + [\mathbf{0}]\_{2 \mathbf{N}\_{\Gamma}}\boldsymbol{\xi}(t) \tag{10}$$

#### *2.3. Augmented Equations of Motion*

This section further arranges the equations derived for the interaction problem in Equation (9) to yield an augmented system in which the primary input is white noise. The white noise input allows a much simpler calculation and the use of stochastic analyses, when compared with manually inputting the roughness profile to the system. Authors in [35,36] constructed a state-space model for the stochastic roughness, when it follows a specified PSD, (Sζζ(ω)). In their approach, the transfer function (Hζ*w*(ω)) was approximated using polynomial representation, as below:

$$\mathcal{S}\_{\zeta\zeta}(\omega) = \mathcal{H}^2\_{\zeta\nu}(\omega)\mathcal{S}\_0 \tag{11}$$

where *S*<sup>0</sup> is the degree of unevenness and ω is radian per second.

Then, <sup>H</sup>ζ*w*(ω)is realized in a state-space model as below to have output vector *<sup>y</sup>*<sup>f</sup> <sup>=</sup> ξˆ(t) . ξˆ(t) *T* , i.e., .

$$
\dot{\mathbf{x}}\_{\mathbf{f}} = \mathbf{A}\_{\mathbf{f}} \mathbf{x}\_{\mathbf{f}}(t) + \mathbf{B}\_{\mathbf{f}} w(t) \tag{12}
$$

$$y\_{\mathbf{f}}(t) = \mathbf{C}\_{\mathbf{f}} \mathbf{x}\_{\mathbf{f}}(t) \tag{13}$$

where *x*<sup>f</sup> is the state vector, **A**f, **B**f, and **C**<sup>f</sup> are system, input, and observation matrices, respectively. The output vector is defined to contain the roughness and time derivative term of the roughness, i.e., *y*f(*t*) = ξ (*t*) . ξ(*t*) T . An example of designing a pavement filter using a second-order low-pass filter and polynomial approximation approaches is discussed in detail [35].

Finally, by combining Equations (9), (10), (12), and (13), the augmented state vector was defined as *<sup>x</sup>*<sup>a</sup> <sup>=</sup> *x*T <sup>T</sup> *<sup>x</sup>*<sup>T</sup> f T . Then,

$$\dot{\mathbf{x}}\_{\mathbf{a}} = \mathbf{A}\_{\mathbf{a}} \mathbf{x}\_{\mathbf{a}} + \mathbf{B}\_{\mathbf{a}} w(t) = \begin{bmatrix} \mathbf{A}\_{\mathbf{T}} & \mathbf{B}\_{\mathbf{T}1} \mathbf{C}\_{\mathbf{f}1} + \mathbf{B}\_{\mathbf{T}2} \mathbf{C}\_{\mathbf{f}2} \\ 0 & \mathbf{A}\_{\mathbf{f}} \end{bmatrix} \begin{bmatrix} \mathbf{x}\_{\mathbf{T}} \\ \mathbf{x}\_{\mathbf{f}} \end{bmatrix} + \begin{bmatrix} 0 \\ \mathbf{B}\_{\mathbf{f}} \end{bmatrix} w(t) \tag{14}$$

$$\mathbf{y}\_{\mathsf{T}} = \mathbf{C}\_{\mathsf{a}} \mathbf{x}\_{\mathsf{a}} = \left[ \begin{array}{c} \mathbf{C}\_{\mathsf{T}} & \mathbf{D}\_{\mathsf{T}1} \mathbf{C}\_{\mathsf{f}1} + \mathbf{D}\_{\mathsf{T}2} \mathbf{C}\_{\mathsf{f}2} \\\\ \mathbf{x}\_{\mathsf{f}} \end{array} \right] \begin{array}{c} \mathbf{x}\_{\mathsf{T}} \\\\ \mathbf{x}\_{\mathsf{f}} \end{array} \tag{15}$$

Note that **B**T1 and **B**T2 indicate the first and the second columns of **B**T, respectively. Similarly, **D**T1 and **D**T2 correspond to the first and the second columns of **D**T. The equations do not contain the feedthrough terms, implying that the system was strictly proper. In addition, **A**<sup>a</sup> and **C**<sup>a</sup> were time-dependent matrices, due to the variable foundation coefficients, kf(*Vt*) and cf(*Vt*).

#### **3. Stochastic Vehicle Response**

The covariance of the augmented system, **Γ***xa* , could be determined through a linear differential equation, when the input is a white noise process [37]:

$$\dot{\boldsymbol{\Gamma}}\_{\boldsymbol{x}\_{\boldsymbol{a}}}(t) = \mathbf{A}\_{\mathbf{a}}(t)\boldsymbol{\Gamma}\_{\boldsymbol{x}\_{\boldsymbol{a}}} + \boldsymbol{\Gamma}\_{\boldsymbol{x}\_{\boldsymbol{a}}}(t)\mathbf{A}\_{\mathbf{a}}^{\mathrm{T}}(t) + 2\pi\mathbf{S}\_{0}\mathbf{B}\_{\mathbf{a}}\mathbf{B}\_{\mathbf{a}}^{\mathrm{T}} \tag{16}$$

where *S*<sup>0</sup> is the level of the white noise indicating the level of roughness. Solving Equation (16) is beneficial as it does not contain the principal matrix, in which an explicit format of the matrix is generally unknown in time-varying systems.

In the case of uniform foundation, i.e., kf(*x*) = kf, cf(*x*) = cf with an infinite length beam, the system becomes stationary. Then, assuming that the initial conditions could be described by a random vector, *x*a(0) = *x*a0 the initial condition of Equation (16) Γ*x*<sup>a</sup> (0) = **Γ**<sup>0</sup> becomes:

$$\mathbf{T}\_0 = E[(\mathbf{x}\_{\rm a0} - \mu\_{\rm x\_{\rm a0}})(y\_{\rm a0} - \mu\_{\rm x\_{\rm a0}})^\mathrm{T}] \tag{17}$$

where μ*xa*<sup>0</sup> and **Γ**<sup>0</sup> indicates the mean and the covariance, respectively. If the initial conditions are all deterministic, Γ<sup>0</sup> = 0. Then, the covariance of the structure responses, Γy, is given by:

$$
\Gamma\_\mathbf{y} = \mathbf{C}\_\mathbf{a} \Gamma\_{\mathbf{x}\_\mathbf{a}} \mathbf{C}\_\mathbf{a}^\mathbf{T} \tag{18}
$$

Further, with a zero-mean white noise being the input to the augmented system in Equation (14), the stationary covariance responses could be obtained by the solution of

$$0 = \mathbf{A\_{a}}\Gamma\_{\mathbf{x\_{a}}} + \Gamma\_{\mathbf{x\_{a}}}\mathbf{A\_{a}^{T}} + 2\pi\mathbf{B\_{a}}S\_{0}B\_{a}^{T} \tag{19}$$

which is known as the Lyapunov equation [38]. Note that the equation is linear in unknown covariances and can only examine the moments of the responses under the stationary process.

However, the presented study consisted of a non-uniform foundation, in which the quantities varied over the length of the beam. Thus, the basic assumptions made in Equation (20) was no longer valid. Instead of directly integrating Equation (14), the general covariance response in Equation (16) was solved for **Γ***xa* , for which the matrix components are described below:

$$
\Gamma\_{\mathbf{x}a} = \begin{bmatrix}
\Gamma\_{\mathbf{q}\mathbf{q}} & \Gamma\_{\mathbf{q}\mathbf{u}} & \Gamma\_{\mathbf{q}\mathbf{z}} & \Gamma\_{\mathbf{q}\dot{\mathbf{q}}} & \Gamma\_{\mathbf{q}\dot{\mathbf{u}}} & \Gamma\_{\mathbf{q}\dot{\mathbf{z}}} & \Gamma\_{\mathbf{q}\mathbf{f}} \\
& \Gamma\_{\mathbf{u}\mathbf{u}} & \Gamma\_{\mathbf{u}\mathbf{z}} & \Gamma\_{\mathbf{u}\dot{\mathbf{q}}} & \Gamma\_{\mathbf{u}\dot{\mathbf{z}}} & \Gamma\_{\mathbf{u}\dot{\mathbf{z}}} \\
& & \Gamma\_{\mathbf{z}\mathbf{z}} & \Gamma\_{\mathbf{z}\dot{\mathbf{q}}} & \Gamma\_{\mathbf{z}\dot{\mathbf{u}}} & \Gamma\_{\mathbf{z}\dot{\mathbf{z}}} & \Gamma\_{\mathbf{z}\mathbf{f}} \\
& & & \Gamma\_{\mathbf{q}\dot{\mathbf{q}}} & \Gamma\_{\mathbf{q}\dot{\mathbf{u}}} & \Gamma\_{\mathbf{q}\dot{\mathbf{z}}} & \Gamma\_{\mathbf{q}\mathbf{f}} \\
& & & & & & \Gamma\_{\mathbf{z}\dot{\mathbf{z}}} & \Gamma\_{\mathbf{z}\mathbf{f}} \\
& & & & & & & \Gamma\_{\mathbf{f}\mathbf{f}}
\end{bmatrix}
\tag{20}
$$

where **qB**, *uu*, and z*<sup>s</sup>* are short noted as **q**, u, and z, respectively. Then, the symbolic covariance matrix, for which the number of variables is Nvar = (Na × (*N*<sup>a</sup> + 1)/2) was plugged into Equation (15) to construct Nvar distinct differential relationships.

Finally, the desired time-varying covariance responses of the vehicle, **Γ***y***<sup>T</sup>** (*t*) = **Γ**uu(*t*) **Γ**zz(*t*) **Γ . u . <sup>u</sup>**(*t*) **Γ. z . <sup>z</sup>**(*t*) T , was obtained by solving Equation (16) via the time-step integration method embedded in Matlab®(e.g., ode45).

#### **4. Illustrative Examples**

This section demonstrates the proposed approach by examining the covariance response of a vehicle on a non-uniform foundation. To first validate the solution procedure, steady-state covariance responses were compared by slowing the speed of the vehicle. Then, covariance responses were compared for various pavement scenarios. Finally, covariance response features were selected to classify and examine the foundation properties.

#### *4.1. Vehicle and Pavement Model Properties*

The properties of the quarter-car model used in the numerical examples are drawn from [39,40] and summarized in Table 1. Note that *k*<sup>t</sup> is sought using a calibration index to well approximate the empirical model in [41], on a non-deformable rigid foundation with varying roughness [35].


**Table 1.** Vehicle properties.

The transfer function to approximate the PSD of road roughness, as in Equation (11), is considered as follows:

$$\mathcal{S}\_{\xi\xi}(\omega) = \mathcal{S}\_0 \left(\frac{\Omega}{\Omega\_0}\right)^{-\nu} \tag{21}$$

where Ω is the spatial circular frequency (ω/V); Ω<sup>0</sup> = 1 rad/m, ν is the waviness that is taken as 2.45, to match the average roads in the U.S. [42]. Note that *S*<sup>0</sup> is varied to match the International Roughness Index (IRIs), which is a measure of road roughness on ride comfort [43]. A lower IRI value indicates a smooth pavement, while a higher IRI implies a rough pavement. In this study, IRIs are varied from 1 to 5, and the corresponding *S*0's are approximated at those integers, using the golden car approach described in [35].

The typical pavement system was adopted herein, and the uniform properties of the Euler–Bernoulli beam (top layer in Figure 1) are summarized in Table 2. The elasticity of the top-layer used in the study represents the medium soil [44].


**Table 2.** Euler–Bernoulli beam property.

Finally, to accommodate different Winkler type foundations, spring and damping coefficients were varied linearly or quadratically. Thus, the following equations were adopted for each case:

$$\mathbf{k}\_{\mathbf{f}}(\mathbf{x}) = k\_{\mathbf{f}0} \times \mathbf{z}\_{\mathbf{f}}(\mathbf{x}) \mathbf{c}\_{\mathbf{f}}(\mathbf{x}) = \mathbf{c}\_{\mathbf{f}0} \times \mathbf{z}\_{\mathbf{f}}(\mathbf{x}) \tag{22}$$

$$\text{Linear:}\,\mathbf{z}\_{\mathbf{f}}(\mathbf{x}) = 1 - a\mathbf{x}, \qquad 0 \le \mathbf{x} \le \mathbf{L}/2\\\mathbf{z}\_{\mathbf{f}}(\mathbf{x}) = \left(2\mathbf{a} - 1\right) - 2\mathbf{x}(\mathbf{a} - 1)/\mathbf{L}, \qquad \mathbf{L}/2 \le \mathbf{x} \le \mathbf{L} \tag{23}$$

$$\text{Quadratic: }\text{ z}\_{\text{f}}(\mathbf{x}) = 1 + 4\mathbf{x}(a-1)/L - (a-1)4\mathbf{x}^{2}/L^{2}, \text{ } 0 \le \mathbf{x} \le L \tag{24}$$

where the reference parameters for stiffness (*k*f0) and damping (*c*f0) are 30 kPa/mm and 2.4 <sup>×</sup> <sup>10</sup><sup>7</sup> <sup>N</sup>·s/m, respectively. The reduction factor α ≤ 1 was selected such that the soil had the most reduced value at the mid-span of the beam (*L*/2). In this study, α was varied from 0.5, 0.7, and 0.9, which implied 50%, 70%, and 90% of the reference parameters. An illustration of zf(*x*) for each α is shown in Figure 2. Herein, the variation profiles are described with L and Q, for a linear and quadratic shape, respectively, followed by two digits describing α. For example, a dashed-dot plot in Figure 2a (linear with α = 0.5)

was denoted as 'L50'. Similarly, 'Q90' indicates that zf(*x*) varies in a quadratic manner (see a solid line in Figure 2b).

**Figure 2.** Nonlinear variation in the Winkler foundation (**a**) linear variation; (**b**) quadratic variation.

#### *4.2. Validation of the Solution Approach*

In this section, the steady-state responses of the vehicle were compared with that of near-stationary responses. The purpose of the presented study was to validate the time-varying covariance solutions in Equation (16) by slowing the speed of the vehicle, and also to illustrate the difference in the response due to various profiles of the subgrade and roughness, to be used in the following sections.

Figure 3 shows **Γ***y*<sup>T</sup> (*t*) in comparison with steady-state responses. Here, the speed of the vehicle was slowed to 0.5 m/s (1.8 km/h). The total number of modes used in the Euler–Bernoulli beam was 10 sine modes and a roughness of IRI = 3 m/km was used. Regarding subgrade, the uniform foundation was represented by zf(*x*) = 1 in Equation (19), while α = 0.5 was used for linearly and quadratically varying foundations. Steady-state responses were calculated by fixing the location of the vehicle at the mid-span and solving Equation (16) with **. Γ***xa* (*t*) = 0, where the responses were as small as <sup>Γ</sup>*uuuu* = 1.39 <sup>×</sup> 10−<sup>5</sup> [m/s] 2 ; <sup>Γ</sup>*zszs* = 1.40 <sup>×</sup> 10−<sup>5</sup> [m/s] 2 ; Γ . *uu* . *uu* <sup>=</sup> 1.42 <sup>×</sup> <sup>10</sup>−<sup>5</sup> [m/s] 2 ; Γ. *zs* . *zs* <sup>=</sup> 0.189 <sup>×</sup> <sup>10</sup>−<sup>5</sup> [m/s] 2 . Note that the impact of foundation property change in the case of a stationary process was negligible, as reported by [24]. However, non-stationary covariance responses showed large humps within the first 2 s, due to the dynamic effect of boundary conditions. Peaks on velocity covariance responses were much higher than the displacement responses, rapidly converging to stationary responses.

**Figure 3.** Covariance response of the vehicle over time (**a**) Γ*uuuu* ; (**b**) Γ*zszs* ; (**c**) Γ . *uu* . *uu* ; and (**d**) Γ. *zs* . *zs* .

Furthermore, the responses were examined on various road roughness. To compare the results efficiently, Figure 4a plots Γ*uuuu* with linearly varying foundation, as the IRIs varied from 1 to 5 m/km. A time-step integration method, ode45 was used as the vehicle crossed over the length with vehicle speed (V = 20 km/h). As could be seen, the effect of surface roughness on the nonstationary covariance responses were negligible. To better visualize the difference in the responses, the differences were plotted in Figure 4b, in percentage. The difference was estimated for each IRI (Γ*IRIi* ), with respect to the response at IRI = 1 km/m (Γ*IRI*<sup>1</sup> ), as below:

$$
\Delta\Gamma\_{IRI} = \left(\frac{\Gamma\_{IRI\_i} - \Gamma\_{IRI\_1}}{\Gamma\_{IRI\_1}}\right) \ast 100 \left[\% \right] (i = 1, 2, \dots, 5) \tag{25}
$$

ΔΓ*IRI* tended to diverge as the vehicle moved along the beam, indicating that dynamic responses were accumulated. However, within the domain, the maximum difference was less than 0.025% which was negligible, compared to the governing dynamic responses. Although the impact of change in IRI might increase as the beam length gets longer and the speed of the vehicle increases, the result indicated that the responses were governed more by the non-uniform features of the foundation. This fact emphasized the importance of conducting nonstationary response analyses because the stationary response analyses could not capture such differences, as reported by [24].

Based on the study, non-stationary responses under various subgrades converged to steady-state responses with time, while velocity covariance responses showed a faster rate. Additionally, the subgrade variation types affected the vehicle responses while road roughness had a negligible impact. Therefore, in the subsequent sections, road roughness was fixed to IRI = 3 m/km, to examine Γ*uuuu* (*t*) on various α's.

**Figure 4.** (**a**) Covariance response of vehicle under varying IRIs. (**b**) Difference in the response among IRIs in percentage.

#### *4.3. Time-Varying Covariance Responses*

Figure 5 shows the time-varying covariance responses, Γ*uuuu* (*t*), as the vehicle runs over the pavement with different types of foundation properties. The simulations were carried out on the basis of combinations of two variables—(1) foundation profile, implying how the subgrade properties change, i.e., either in a linear or in a quadratic manner, and (2) α, which varied from 0.5, 0.7, to 0.9. Then, the variable time-step was used for the integration and then downsampled to present 100 data points within the duration. After the resampling procedure, one could consider that the responses on the vehicle were measured at about 110 Hz. The chosen sampling rate was low enough to be easily realized, yet could capture the key features of the responses. The response shown in Figure 5 are deterministic, as the foundation variation profiles, α, and the speed of the vehicle are known. However, deterministic identification based on Figure 5 was unrealistic, because the measured signals tended to be contaminated with noises.

**Figure 5.** Covariance response on various foundation properties (without noise).

Now, noises were randomly selected to have a signal to noise level (SNR) between 25–50 dB. Note that the noises were added to the raw signals and then resampling was performed to represent the signal noises. Here, only the measured noises were considered because the zero-mean noises in the responses did not affect the covariance responses. For example, Figure 6 illustrates the covariance response with noises added on linearly varying foundation with α = 50. SNR of Figure 6a is about 45 dB and Figure 6b is about 50 dB.

**Figure 6.** Covariance response with measured noise (**a**) SNR of 45 dB; and (**b**) SNR of 25 dB.

From the outcome of this subsection, the following conclusions could be made:


Thus, to resolve the issue, a machine-learning based classification of subgrades based on Γ*uuuu* (*t*) is discussed in the subsequent section.

#### **5. Machine-Learning Based Classification**

Machine-learning techniques are recognized in the civil engineering field as a promising component for monitoring and inspecting [13]. Machine-learning tools can provide pattern recognition strategies, when a deterministic model is difficult to be identified [45]. With their highlighted importance and computational advances, the Matlab software incorporated the Statics and Machine Learning toolbox containing considerable machine-learning techniques [46,47].

Among classifiers provided in the Classification Learner App in Matlab R2019b, one of the traditional classifier, linear discriminant analysis (LDA) is implemented for identifying the changes in the foundation properties from vehicle responses. The LDA method assumes that the data are distributed in Gaussian and that each attribute has the same variance. Then, the Bayes' theorem is applied to estimate the posterior probability that the observation belongs to a certain class. Then, the costs are evaluated from the maximal difference between the computed sample covariance and the empirical covariance matrix. In LDA, the cost function is linear with respect to the observation [46,48]. With these assumptions, the LDA model attempts to express one dependent variable in terms of a linear combination of other features or measurements [46]. Thus, to enhance the classification accuracy, features in the covariance responses must be selected carefully.

Based on the previously presented results, the following six features were selected—(1) maximum amplitude, A1; (2) time corresponding to A1, T1; (3) minimum tangent occurring between 0.2 s and 0.8 s, A2; (4) time corresponding to A2, T2; (5) slope of the linear regression between 0.2 s and 0.8 s, A3; and (6) y-intercept of the regression found in (5), A4. Then, to incorporate the measured noise in covariance responses, RMS noises were added in Γ*uuuu* (*t*). The extracted features showed some relationships among them. Figure 7 illustrates the distribution of features over the range of A1. As can be seen, some features, such as A2 and A4, show a higher correlation with A1, while other features are more scattered over the range of A1. In addition, L70, L90, and Q90 seem to overlap with each other (as in Figure 7a–e, making it hard to differentiate with classification models that are based on decision trees, etc.

**Figure 7.** The distribution features over the range of A1 (Maximum Amplitude); (**a**) Relationship between A1-T1 (Maximum Time); (**b**) Relationship between A1–A2 (Maximum Tangent); (**c**) Relationship between A1-T2 (Time at Maximum Tangent); (**d**) Relationship between A1–A3 (Linear regression slope); and (**e**) Relationship between A1-A4 (Linear regression y-intercept).

Subsequently, the collected datasets were trained using LDA. The average success rate for using an LDA classifier was over 94%, with at most 10% noise. Table 3 is a confusion matrix when 77 training data were used. The table shows that the foundation was mostly classified, while 18% of the L50 case was misclassified as Q70, and vice versa. Note that LDA showed the highest accuracy when compared with other classification tools; the linear support vector machine showed 82% accuracy, ensemble provided 84% accuracies, while other methods showed over 40% errors.


**Table 3.** Confusion matrix for foundation property test.

Now validation tests were conducted to verify the performance of the developed LDA model. For each case of the foundation, 100 test sets with 1~10% RMS noises added on the responses were generated. The accuracy of the classifier was plotted as a bar chart shown in Figure 8. As could be expected from the confusion matrix, Q70 showed the lowest accuracy, 83%, followed by L50. Overall accuracy was about 94.5% This result supports that by adopting LDA model, the vehicle responses could classify the change in the foundation properties with good accuracies.

**Figure 8.** The success rate for identifying foundation property.

#### **6. Conclusions**

This paper presented a machine-learning-based classification of non-uniform foundation properties using vehicle responses. The dynamics response of the quarter-car model on the stochastic deformable pavement with a finite length was evaluated. A filtered white noise was used to represent the stochastic pavement roughness. The deformable subgrade was modeled by an Euler–Bernoulli beam on a Winkler-type foundation. The non-uniform characteristics were represented with varying stiffness and damping coefficients of the subgrade. Then, the vehicle–pavement–foundation interaction model was combined to yield an augmented state-space representation, which had white noise as the primary input to the system. In this study, the model could accommodate any time of foundation that was describable with a longitudinal axis, although only the impacts of linear and quadratic variations were discussed. A time-varying Lyapunov equation governing the covariance of the responses was solved to effectively obtain the response of the vehicle. From the steady-state Lyapunov solution, the solution approaches were validated. Then, various values of the subgrade's properties, along with surface roughness were compared. The parametric study showed that the stiffness produced some difference in the response profile, while the roughness produced negligible change. This fact opposed the uniform foundation case, indicating the importance of considering the non-uniform foundation. Then, a set of simulations for measuring noise were performed and used for feature extraction for a classifier. Using an application embedded in Matlab®, linear discriminant analysis was employed to show an average accuracy of 94%. Finally, a total of 600 test sets were generated to demonstrate that the estimated foundation properties were mostly correct. Based on the outcome of this study, the contribution of the presented work and the specific conclusions are summarized as follows:


In conclusion, the presented work demonstrated the potentials of monitoring the subgrade anomalies from an inspecting vehicle that is only equipped with a set of accelerometers. Unlike the current approach of a subgrade survey that is only limited to suspicious spots, the successful realization of the presented methodology might allow a complete survey and the construction of a database for road management.

**Author Contributions:** R.E.K. contributed to the conceptualization, methodology, software, funding, and validation. The author has read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the research fund of Hanyang University (HY-2019).

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
