基于多种群遗传算法考虑浅层基岩及粘弹性的路基模量反演方法

张军辉; 刘杰; 范海山; 张石平; 丁乐; ZHANG Junhui; LIU Jie; FAN Haishan; ZHANG Shiping; DING Le

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Back-calculation of subgrade modulus considering shallow bedrock and viscoelasticity based on multi population genetic algorithm PDF

- ORCID：
ZHANG Junhui
✉
- ORCID：
LIU Jie
- ORCID：
FAN Haishan
- ORCID：
ZHANG Shiping
- ORCID：
DING Le

School of Traffic and Transportation Engineering, Changsha University of Science & Technology, Changsha 410114, P. R. China

CLC： U416.1

Updated：2023-03-07

DOI：10.11835/j.issn.2096-6717.2022.021

Abstract

This article proposes a new idea for back-calculation of the subgrade modulus considering shallow bedrock and the viscoelastic characteristics. For the subgrade model, displacement boundary conditions and the Kelvin model are adopted to describe the depth of the shallow bedrock and the viscoelasticity, respectively. The portable falling weight deflectometer (PFWD) field test is simulated by ABAQUS general finite element (FE) software, and the optimal value of the modulus is iterated by a multi population genetic algorithm (MPGA). Based on the new method, back-calculation results from FE simulation tests show that the modulus average error of the forward model considering shallow bedrock is 7.0%, while that of the forward model considering half space is as high as 16.2%, indicating that negletct of the shallow bedrock in the forward model of the back-calculation program may cause a significant error in the inversion modulus, but its influence decreases with the increase of the depth of the shallow bedrock, and the depth limit is 3 m. Similarly, for the FE model considering viscoelasticity, the maximum error for neglect of this attribute in the forward model reaches 27.9%, compared with the error of only 7.4% when considering viscoelasticity in the forward model. Due to the difficulty exploring the depth of shallow bedrock, examinations are only conducted from the theoretical aspect.

Keywords

back-calculation; subgrade modulus; shallow bedrock; viscoelasticity; genetic algorithm

1 Introduction

As the foundation of the pavement structure, durable and stabilized subgrade is critical for the whole road service life. The elastic modulus is the key parameter for pavement design, evaluation of subgrade characteristics^[

1], dynamic response calculation^{[Reference 2

Baidu Scholar}2] and rapid detection of subgrade strength^{[Reference 3

Baidu Scholar}3]. Generally, there are three kinds of methods used to obtain the subgrade modulus: modulus back-calculation, the laboratory dynamic triaxial test ^{[Reference 4-5}4-5] and the modulus prediction model ^{[Reference 6-7}6-7]. The first method is used to gain the subgrade's structure modulus that can reflect the anti-deformation ability, while the latter two can only achieve the material modulus of subgrade soil, and the conversion between material modulus and structure modulus requires empirical relationships.

Before back-calculating the subgrade modulus, the real load and displacement time history curve should be measured. Presently, the falling weight deflectometer (FWD)^[

8] and the portable falling weight deflectometer (PFWD)^{[Reference 1

Baidu Scholar}1] tests are the measuring instruments most widely used for subgrade dynamic deflection detection. However, due to the high density of the FWD load, it will enlarge the maximum strain of the subgrade top, which makes the FWD method have poor adaptability for the detection of subgrade quality. With its lightweight design, not only has the PFWD method made the detection process more convenient, it can simulate the real stress state of the subgrade, which is more beneficial for rapid detection of the bearing capacity of the subgrade.

Many researchers have carried out a lot of exploration work in the modulus back-calculation domain. The methods most commonly used are the peak value method ^[

9], the least square method ^{[Reference 10

Baidu Scholar}10], the shuffled complex evolution (SCE) algorithm ^{[Reference 11

Baidu Scholar}11], the data mining approach ^{[Reference 12-13}12-13], the genetic algorithm (GA) ^{[Reference 14-15}14-15] and the artificial neural network (ANN) method ^{[Reference 16-17}16-17]. Among them, the peak value method has been traditionally applied in the engineering field to back-calculate the elastic modulus of the subgrade by considering the measured peak load and its corresponding deflection, which may cause great systematic errors in the inversion modulus when the subgrade material shows obvious viscoelastic properties ^{[Reference 18

Baidu Scholar}18]. Although the least square method, the SCE algorithm, the data mining approach and the ANN are fast in calculation speed, the accuracy of their results depends on the use of a certain number of data samples ^{[Reference 19

Baidu Scholar}19]. GA, which does not rely on the gradient and number of samples, has a good convergence rate and calculation speed, but poor local search ability ^{[Reference 20-21}20-21]. However, a multi population genetic algorithm (MPGA) can make up for the shortcoming of the standard genetic algorithm well ^{[Reference 22

Baidu Scholar}22].

By using a MPGA, modulus back-calculation analysis can be divided into two parts, selecting a mechanical subgrade model, and carrying out the iteration procedure of matching the measured displacement time history curve. Generally, the subgrade is regarded as a semi-infinite space body ^[

23], and the shallow bedrock existing in the subgrade is difficult to consider in the mechanical model for the subgrade in mountainous highways and areas with a high degree of soil consolidation ^{[Reference 24

Baidu Scholar}24]. Several studies have been conducted on this topic. Uddin et al. ^{[Reference 25

Baidu Scholar}25] developed a back-calculation program for an in-situ modulus considering a shallow bedrock layer in the pavement-subgrade system. The findings showed that neglect of the shallow bedrock may cause gross errors in the subgrade inversion modulus. Briggs et al.^{[Reference 26

Baidu Scholar}26] suggested a four layers pavement-base-subgrade-bedrock system to back-calculate the subgrade modulus. The calculation results indicated that compared with the modulus of the rigid layer, the depth has more influence on the inversion modulus. Hudson et al. ^{[Reference 27

Baidu Scholar}27] proposed that the measurement depth of a rigid layer should be used in the analysis of FWD deflection data, but this method is both costly and impractical because a stiff layer is not only depth-unknown, but kind-various ^{[Reference 28

Baidu Scholar}28]. Seok et al. ^{[Reference 29

Baidu Scholar}29] back-calculated the elastic modulus of the subgrade when investigating the depth of the shallow bedrock. The findings showed that when the depth of the bedrock exceeds 4.0 m, its effects on the subgrade back-calculation modulus are negligible. Uzan et al.^{[Reference 30

Baidu Scholar}30] found that the depth of the rigid layer is 5 m.

The above scholars based their results on the peak value method and linear elastic theory in their research on the influence of the rigid bedrock layer on the subgrade inversion modulus. There are little researches that consider the viscoelastic characteristics of the subgrade material. However, due to the wide distribution of cohesive soil in South China ^[

31], only considering the linear elastic model of the subgrade will cause significant errors to accrue in the modulus inversion value. In addition, most of the existing subgrade modulus back-calculation studies considering shallow bedrock focus on the deflection basin analysis of the overall pavement structure, while few researches take the displacement and load time history curves on the subgrade surface as the iterative parameters for the MPGA. For the depth of the bedrock layer, because its influence degree on the subgrade inversion modulus is great, this factor is selected in this study for modulus inversion research. Finally, the researches on the rigid bedrock layer were mainly in the 1980s and 1990s. Based on the above descriptions, in order to explore the influence of the shallow bedrock layer on the subgrade inversion modulus, a new subgrade model considering the shallow bedrock layer and viscoelastic characteristics of subgrade soil that combines the MPGA and ABAQUS was put forward. For this model, on the premise of ignoring the modulus of the shallow bedrock layer, the depth is simplified as a boundary condition and used for the ABAQUS modeling. It is proved to be of high accuracy and is a reasonable simplification of the existing methods. The integral numerical expressions of mechanical model response under the time domain are derived. This method can be adapted to different subgrades considering the shallow bedrock layer by changing the physical parameters of the subgrade. Although this paper only explores the influence of the rigid bedrock depth and subgrade viscoelastic characteristics on the modulus back-calculation value, it can help follow-up researchers understand the necessity of considering the above two factors in the mechanical model. To optimize the calculation accuracy and convergence speed, back-calculation of the subgrade modulus programmed by MATLAB has been developed, which can provide a powerful reference for further practical application.

2 Descriptions of mechanical model and back-calculation method

For the dynamic mechanical model considering the rigid bedrock of the subgrade, when the PFWD test is conducted on the top of the subgrade, the dynamic wave generated will be reflected due to the existence of the shallow rigid bedrock layer, which will affect the accuracy of the inversion modulus. Sirithepmontree et al. ^[

32] found that the existence of a stiff layer had a great influence on the soil dynamic response, and the back-calculation modulus of the subgrade mainly depended on its mechanical model. Therefore, it is necessary to consider this factor in the mechanical model when shallow bedrock exists in the real situation.

2.1 Mechanics forward model

In this study, the subgrade is regarded as a single-layer structure with a shallow rigid bedrock layer. A diagram of the subgrade mechanical model and load distribution is presented in Fig. 1. Under the axisymmetric coordinate system, the vertical displacement expression under the impact load based on the rigid plate is obtained by combining the corresponding boundary conditions after solving the dynamic governing equations and constitutive equations.

Fig. 1 Schematic diagram of the subgrade mechanical model and load distribution

Note: E and η refer to the elastic modulus of the spring and the viscosity coefficient of the dashpot in the Kelvin model under the time domain; μ is the Poisson's ratio; ρ is the dynamic density; h₀ is the depth of the rigid bedrock layer.

Subgrade soil is a kind of visco-elastoplastic material, and its plastic deformation will be eliminated by multiple drop hammer impacts. Due to the short impact time of the PFWD load ^[

33], the subgrade can be regarded as a viscoelastic body. The Kelvin model is employed in this study to describe the viscoelastic properties of the subgrade. At the beginning, the external constant force is borne by the damper, and the stress on the spring is zero. However, with the increase of time, the system will reach a stable level, and the external force will be gradually transferred from the damper to the spring. In this case, it takes more than a single spring system for the Kelvin model to reach a stable stress-strain level, which causes the "hysteresis" phenomenon of dynamic response. On the premise of ensuring the accuracy of calculation, the Kelvin model can improve the efficiency of the calculation and reduce the complexity of the computation. The expression of the modulus E(s) under the frequency domain is as follows.

E (s) = E + η s

（1）

where s is the parameter of the Laplace transform.

For a long time, during the pavement design process, in order to simulate the vertical effect of automobile tires on the pavement structure, researchers have put forward various load distribution forms, such as circular, hemispherical and "bowl shape". Based on these traditional load forms, Gerrard et al. ^[

34] proposed the general form of the arbitrary vertical axisymmetric load （Eq. （2））, but it can only express the horizontal distribution form of the load. Actually, for PFWD, its input load is changing over time, according to the least square method, and it can be fitted in the form of a half sine wave （Fig. 2）.

p (r, t) = \{\begin{array}{l} m p (t) {(1 - \frac{r^{2}}{δ^{2}})}^{m - 1}, r < δ \\ 0, r > δ \end{array}

（2）

Fig. 2 Schematic diagram of half sine curve fitting

Note: According to the principle that the area enclosed by the curve and the time axis are equal, the measured curve is fitted by a half sine curve, and the fitting parameters p_max and T₀ are obtained.

where p（t）=p_maxsin（πt/T₀）（H（t）-H（t-T₀））, T₀, p_max and H（∙） mean the action time of the PFWD impact load, the peak pressures of the distributed load concentration and the unit step function, respectively; p（r,t） is the expression of load distribution, r is the distance from the calculation point to the center of the loading plate; m is the load type coefficient. For the load under the rigid bearing plate, m=0.5; δ is the radius of loading plate; p（t） is the distributed load concentration on the upper surface of the loading plate. For a better understanding of the foregoing, Fig. 2 is given below.

Several common load expressions are included in this general expression of axisymmetric vertical loads （Eq. （2））. For example, when m=0.5, 1, 1.5, it represents the stress distribution forms under the rigid bearing plate, of the circular uniform vertical load and of the hemispherical vertical load, respectively. For the circular uniform distributed load, the load concentration changes abruptly at the edge of the plate （r=δ）, which makes some theoretical results inconsistent with the actual physical phenomena, referred to as the "singular point" problem in mathematics. And for the hemispherical load, although it is continuous around the periphery, the disadvantage is that the stress is too concentrated at the midpoint of the loading plate, which easily leads to inconsistency between theory and practice. For the stress distribution form under the rigid loading plate, its shape is a parabola, which is suitable to simulate the interaction between the subgrade and the loading plate. The relevant load expression after the Hankel-Laplace transform is as follows.

\begin{array}{l} \tilde{\bar{p}} (ξ, s) = \frac{π T_{0}}{s^{2} {T_{0}}^{2} + π^{2}} (1 + e^{- s T_{0}}) \frac{p_{m a x} δ s i n (ξ δ)}{2 ξ} \end{array}

（3）

where $ξ$ is the parameter of the Hankel integral transform. Thus, the boundary condition of the subgrade mechanical model under the Hankel-Laplace domain can be shown below.

\{\begin{array}{l} {\tilde{\bar{u}}}_{r 1} (ξ, h_{0}, s) = 0 \\ {\tilde{\bar{ω}}}_{z} (ξ, h_{0}, s) = 0 \\ {\tilde{\bar{σ}}}_{z} (ξ, 0, s) = - \tilde{\bar{p}} (ξ, s) \\ {\tilde{\bar{σ}}}_{z r} (ξ, 0, s) = 0 \end{array}

（4）

For the dynamic mechanical model of the subgrade, differential equations of motion, geometric equations and physical equations under the axisymmetric problem are shown as Eq. (5) to Eq. (7), respectively.

\{\begin{array}{l} \frac{\partial σ_{r} (r, z, t)}{\partial r} + \frac{\partial τ_{z r} (r, z, t)}{\partial z} + \frac{σ_{r} (r, z, t) - σ_{θ} (r, z, t)}{r} = ρ \frac{\partial^{2} u_{r} (r, z, t)}{\partial t^{2}} \\ \frac{\partial σ_{z} (r, z, t)}{\partial z} + \frac{\partial τ_{z r} (r, z, t)}{\partial r} + \frac{τ_{z r} (r, z, t)}{r} = ρ \frac{\partial^{2} w (r, z, t)}{\partial t^{2}} \end{array}

（5）

A = {(\begin{matrix} \frac{\partial}{\partial r} & \frac{1}{r} & 0 & \frac{\partial}{\partial z} \\ 0 & 0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial r} \end{matrix})}^{T} (\begin{matrix} u_{r} (r, z, t) \\ w (r, z, t) \end{matrix})

（6）

(\begin{matrix} σ_{r} (r, z, t) \\ σ_{θ} (r, z, t) \\ σ_{z} (r, z, t) \\ τ_{z r} (r, z, t) \end{matrix}) = (\begin{matrix} λ + 2 G & λ & λ & 0 \\ λ & λ + 2 G & λ & 0 \\ λ & λ & λ + 2 G & 0 \\ 0 & 0 & 0 & G \end{matrix}) A

（7）

where $σ_{r} (r, z, t)$ , $σ_{θ} (r, z, t)$ and $σ_{z} (r, z, t)$ are the normal stresses in the r, θ and z direction, respectively; $τ_{z r} (r, z, t)$ is the shear stress in the z-r plane; $u_{r} (r, z, t)$ and $w (r, z, t)$ are the displacements in the r and z direction, respectively; $A = {(\begin{matrix} ε_{r} & ε_{θ} & ε_{z} & γ_{z r} \end{matrix})}^{T}$ , in which $ε_{r} (r, z, t)$ , $ε_{θ} (r, z, t)$ and $ε_{z} (r, z, t)$ are the normal strains in the r, θ and z direction, respectively, and $γ_{z r} (r, z, t)$ is the shear stress in the z-r plane; λ=μE/(1+μ)(1+2μ); G=2E/(1+μ). Dynamic basic equations under the frequency domain can be obtained by performing a Laplace-Hankel transform on Eq. （5） to Eq. （7）. Combining them, the following equations can be obtained.

\{\begin{array}{l} (\frac{d^{2}}{d z^{2}} - m_{1}) {\tilde{\bar{u}}}_{r 1} (ξ, z, s) = m_{2} \frac{d \tilde{\bar{ω}} (ξ, z, s)}{d z} \\ (\frac{d^{2}}{d z^{2}} - m_{3}) \tilde{\bar{ω}} (ξ, z, s) = - m_{4} \frac{d {\tilde{\bar{u}}}_{r 1} (ξ, z, s)}{d z} \end{array}

（8）

where m₁=（λ（s）+2G（s）ξ²+ρs²）/G（s）; m₂=ξ（λ（s）+G（s））/G（s）; m₃=（（ρs²+G（s）） ξ²）/（λ（s）+ 2G（s））; m₄=ξ（λ（s）+G（s））/（λ（s）+2G（s））; λ（s）=μE（s）/（（1+μ）（1+2μ））; G（s）=2E（s）/（1+μ）; ${\tilde{\bar{u}}}_{r 1} (ξ, z, s)$ and $\tilde{\bar{ω}} (ξ, z, s)$ are the displacements in the r and z direction, respectively. Its expressions are shown below.

\{\begin{array}{l} {\tilde{\bar{u}}}_{r 1} (ξ, z, s) = \int_{0}^{+ \infty} \int_{0}^{+ \infty} r u_{r} (r, z, t) e^{- s t} J_{1} (ξ r) d t d r \\ \tilde{\bar{ω}} (ξ, z, s) = \int_{0}^{+ \infty} \int_{0}^{+ \infty} r ω (r, z, t) e^{- s t} J_{0} (ξ r) d t d r \end{array}

（9）

By solving Eq. （8）, the following equations can be obtained.

\{\begin{array}{l} {\tilde{\bar{u}}}_{r 1} (ξ, z, s) = A e^{α z} + B e^{- α z} + C e^{δ z} + D e^{- δ z} \\ \tilde{\bar{ω}} (ξ, z, s) = k_{1} A e^{α z} - k_{1} B e^{- α z} + k_{2} C e^{δ z} - k_{2} D e^{- δ z} \end{array}

（10）

where k₁=（α²-m₁）/（m₂α）; k₂=（δ²-m₁）/（m₂δ）. The expressions of α and δ are as follows, respectively.

α, δ = \sqrt[]{((m_{3} + m_{1} - m_{2} m_{4}) \pm \sqrt[]{{(m_{3} + m_{1} - m_{2} m_{4})}^{2} - 4 m_{1} m_{3}}) / 2}

（11）

Combining Eq. （10） with the physical equation and the geometric equation under the Hankel-Laplace domain, the expressions of ${\tilde{\bar{σ}}}_{z} (ξ, z, s)$ and ${\tilde{\bar{σ}}}_{z r} (ξ, z, s)$ are as follows.

\{\begin{array}{l} {\tilde{\bar{σ}}}_{z} (ξ, z, s) = k_{3} A e^{α z} + k_{3} B e^{- α z} + k_{4} C e^{δ z} + k_{4} D e^{- δ z} \\ {\tilde{\bar{σ}}}_{z r} (ξ, z, s) = k_{5} A e^{α z} - k_{5} B e^{- α z} + k_{6} C e^{δ z} - k_{6} C e^{- δ z} \end{array}

（12）

where k₃=G（s）ξ+（λ（s）+G（s））（（α²-m₁）/m₂）;k₄=G（s）ξ+（λ（s）+G（s））（（δ²-m₁）/m₂）; k₅=（α- ξ（α²-m₁）/（αm₂））G（s）; k₆=（δ-ξ（δ²-m₁）/（δm₂））G（s）. The undetermined parameters A, B, C and D can be solved by combining Eq. （4）, Eq. （10） and Eq. （12）.

\{\begin{array}{l} A = - \frac{\tilde{\bar{p}} (ξ, s) (e^{- (α - δ) h_{0}} k_{6} (k_{1} + k_{2}) + e^{- (α + δ) h_{0}} k_{6} (k_{1} - k_{2}) - 2 k_{2} k_{5})}{Δ} \\ B = - \frac{\tilde{\bar{p}} (ξ, s) (e^{(α + δ) h_{0}} k_{6} (k_{1} - k_{2}) + e^{(α - δ) h_{0}} k_{6} (k_{1} + k_{2}) - 2 k_{2} k_{5})}{Δ} \\ C = - \frac{\tilde{\bar{p}} (ξ, s) (e^{(α - δ) h_{0}} k_{5} (k_{1} + k_{2}) - e^{- (α + δ) h_{0}} k_{5} (k_{1} - k_{2}) - 2 k_{1} k_{6})}{Δ} \\ D = - \frac{\tilde{\bar{p}} (ξ, s) (e^{(α + δ) h_{0}} k_{5} (k_{2} - k_{1}) + e^{- (α - δ) h_{0}} k_{5} (k_{1} + k_{2}) - 2 k_{1} k_{6})}{Δ} \end{array}

（13）

where ∆=(e⁽^α+δ⁾^h⁰+e^-(α+δ)^h⁰)(k₁-k₂)(k₃k₆-k₄k₅)+(e^(α-δ)^h⁰·(k₁+k₂)+e^-(α-δ)^h⁰k₁)(k₃k₆+k₄k₅)-4(k₁k₄k₆+k₂k₃k₅).

Due to the complexity of the displacement expression, it can only be solved by the numerical integration method. First, substitute Eq. （13） into Eq. (10), and then perform inverse Hankel integral transform on it. Finally, the numerical solution $\bar{ω} (ξ, 0, s)$ can be calculated using the Gaussian interpolation algorithm.

\begin{array}{l} \bar{ω} (r, 0, s) = \int_{0}^{\infty} ξ \tilde{\bar{ω}} (ξ, 0, s) J_{0} (ξ r) d ξ \approx \\ \int_{0}^{x_{s}} ξ \tilde{\bar{ω}} (ξ, 0, s) J_{0} (ξ r) d ξ \approx \\ \sum_{i = 1}^{N} \sum_{j = 0}^{20} A_{j} [x_{i j} \bar{ω} (x_{i j}, 0, s) J_{0} (x_{i j} r)] \end{array}

（14）

where $x_{s}$ is the upper limit of the numerical integration function; N is the total number of integration segments and ΔL is the length of each segment, N= $x_{s}$ /ΔL; A_j is the weight coefficient of the quadrature formula of the Gauss type, in order to ensure the accuracy of the calculation, the 20 node Gauss interpolation formula is adopted; x_j is the integral node of the Gaussian interpolation algorithm, x_ij=(i-1)ΔL+ x_j.

In order to meet different calculation accuracies, many numerical solution methods of the inverse. Laplace transform have been developed for the dynamic layered method, the fast Fourier transform method and the Dehin method. Due to the complexity of the displacement components under the frequency domain and the high precision requirement of the solutions, the numerical method with complex expression is generally used for the inverse Laplace transform. In this study, the Debin method is adopted, and the expression is shown below.

\begin{array}{l} ω (r, 0, t_{j}) = \frac{2 e^{a j Δ t}}{T} R e \{- \frac{1}{2} [\bar{ω} (r, 0, a)] + \\ \sum_{k = 0}^{N - 1} A_{k} {(c o s \frac{2 π}{N} + j s i n \frac{2 π}{N})}^{j k}\} \end{array}

（15）

where $A_{k} = \sum_{m = 0}^{L} \{\bar{f} (r, 0, a + i (k + m N) \frac{2 π}{T})\}$ ; T is the total calculation time; N is the calculating steps; t_j is the calculation time required for j time steps （j=0,1,2 $\dots$ N-1）, t_j=j∙T/N; L and a are the corresponding calculation parameters of the Durbin method, L∙N=50-5 000, a∙T=5-10; i is the imaginary unit, i²=-1; Re（·） is a function that takes the real part of a complex number.

2.2 Modulus back-calculation method

Genetic algorithm （GA） is a random search algorithm based on the biological natural selection mechanism. The optimization process has strong robustness and global search ability because it does not depend on gradient. Due to the complexity of the subgrade mechanical response expression, in the field of modulus back calculation, the use of a GA has a great advantage. However, although the GA has good global search ability, it has poor local search ability and is prone to premature convergence ^[

35], but the MPGA can make up for its shortcomings ^{[Reference 36

Baidu Scholar}36]. Guan et al. ^{[Reference 37

Baidu Scholar}37] proved theoretically that the GA based on elitist strategy has global convergence, but it is easy to make the population lose diversity in the later stage of the optimization process. Hussain et al. ^{[Reference 38

Baidu Scholar}38] and Yu et al. ^{[Reference 39

Baidu Scholar}39] showed that the fitness function based on linear transformation is helpful to solve the problem of premature convergence as it can improve individual differences, and the real number coding method can expand the value range of the parameters and eliminate the decoding time, thus improving the accuracy of the calculation and the speed of the convergence. Therefore, on the basis of multi-population strategy, this paper considers the real coding method, elitist preservation strategy and linear transformation of the fitness function to achieve the purpose of finding the optimal global solution.

For the process of modulus back-calculation, to find the optimal value, an MPGA is used to obtain the minimum variance between the real displacement time history curve series $A_{1} = [\bar{ω} (t_{1}), \bar{ω} (t_{2}) \dots \bar{ω} (t_{n})]$ and the theoretical displacement time history curve series $A_{2} = [ω_{i} (t_{1}), ω_{i} (t_{2}) \dots ω_{i} (t_{n})]$ . The influence of the peak displacement difference is also considered. Based on this goal, the specific expression to evaluate the fitness of each generation is shown in Eq. （16）. The fitness function based on linear transformation is shown in Eq. （17）.

\begin{array}{l} Q_{i} = \frac{1}{\sqrt[]{\frac{1}{n} \sum_{j = 1}^{n} {(\bar{ω} (t_{j}) - ω_{i} (t_{j}))}^{2}}} + \\ \frac{1}{|m a x (\bar{ω} (t)) - m a x (ω_{i} (t))|} \end{array}

（16）

ξ_{i} = \frac{100 \times (Q_{i} - Q_{m i n})}{Q_{m a x} - Q_{m i n}}

（17）

where $Q_{i}$ is the individual fitness of the i-th generation offspring; $Q_{m a x}$ and $Q_{m i n}$ are the maximum and minimum value, respectively; $\bar{ω} (t_{j})$ is the real displacement value at the j-th time point of the generation offspring; $ω_{i} (t_{j})$ is the theoretical displacement value at the j-th time point of the i-th generation offspring; $m a x (\cdot)$ is the maximum value function; $ξ_{i}$ is the individual score of the i-th generation offspring based on the linear transformation method.

Specifically, the steps of MPGA are as follows:

1） Data acquisition: load and displacement time history curves are measured by PFWD, and displacement series A₁ and A₂ are obtained according to a certain time series （t₁, t₂⋯t_n） and Eq. （15）, respectively. The load peak value p_max and impact time T₀ are achieved by fitting the load time history curve with the half sine load.

2） Setting of the initialization parameters of the decimal population: among the optimization range, four populations are randomly generated with 80 individuals in total. The cross probability of each population is p_c=0.95+0.05×rand $(\cdot)$ and the variation probability is p_m=0.02+0.03×rand $(\cdot)$ , where the function $r a n d (\cdot)$ can generate random numbers between 0 and 1; the total genetic algebra is 100 generations; E, η and ρ are a group of genes of the population.

3） Genetic manipulation: the individual fitness of each population is calculated according to Eq. （16） and Eq. （17）, and the optimal individual is recorded. Based on the fitness function values, the new generation of subpopulations is obtained by genetic operations such as selection, crossover and mutation. Then, two populations are randomly selected from these subpopulations, and new populations are generated by immigration operations.

4） Algorithm termination condition: when the genetic algebra reaches 100 generations, the algorithm stops and selects the optimal value of the recorded "optimal individual" as the output result. Otherwise, it returns to step （3）. The whole reverse calculation process is summarized in Fig. 3.

Fig. 3 Back-calculation flow chart

Note:The back-calculation process based on a MPGA is presented on the right side of the flow chart, and the specific process of obtaining individual fitness is shown on the left side of the figure. In order to connect the two sides, a group of genetic genes are randomly generated on the right side and input into the left side for fitness calculation.

3 Theoretical analysis and verification

Because the depth of shallow bedrock is difficult to ascertain, this paper only explores the influence of this factor and subgrade viscoelasticity on the modulus inversion process theoretically, and the shallow bedrock is represented in the form of a rigid constraint in ABAQUS. Based on the above assumptions, a model is developed to simulate the field detection process of the PFWD, in which the computing time and the initial time step of the input load are 0.02 s and 0.000 5 s, respectively. Bedrock depths of 0, 0.5, 1, 2 and 3 m, viscosity coefficients of 0, 54 and 134 kPa∙s, and theoretical subgrade moduli of 30, 40, 50 and 60 MPa are selected for combined analysis.

3.1 Finite element modeling

In order to verify the accuracy of the above subgrade modulus back-calculation method, this study uses the ABAQUS general FE software to simulate the test process of PFWD. As shown in Fig. 4, the subgrade model is composed of a drop hammer, a spring, a loading plate and a subgrade. During the PFWD test, an impact load with a drop weight of 10 kg is applied to the loading plate with a radius of 0.15 m sitting on the subgrade surface. A hard spring is used to connect the falling weight and the loading plate, and its stiffness K is 560×10⁵N/m. At the initial velocity of 4.2 km/h, the drop hammer compresses the spring downward, and applies the impact load to the subgrade by converting kinetic energy into elastic potential energy. At the same time, rigid constraints are set in the bottom of the subgrade to simulate the shallow bedrock. Because the load form, subgrade shape and mechanical response are symmetrical to the same axis, the axisymmetric model is selected in this paper. The axisymmetric FE model of the PFWD and its local enlarged drawing are shown in Fig. 4.

Fig. 4 Axisymmetric FE model of the PFWD and its local enlarged drawings

Note: The whole model is shown in the top left of the figure; the contact behavior is shown in the lower left figure; the figure on the right side shows how the impact load is applied in the model.

In the case of a certain bedrock depth, when the load is applied to the subgrade, a dynamic wave will transmit to the boundary of the FE model to produce a reflection, which will greatly influence the accuracy of the calculation results, thus affecting the modulus inversion results ^[

40]. In order to further determine the size of the model, as shown in Table 1, the model selects the widths of 1, 2, 3, 4, 5, 6, 10, 15 and 20 m for mechanical response analysis. On the premise of ensuring the accuracy of the calculation, this paper introduces the "finite element and infinite element boundary" subgrade model^{[Reference 41

Baidu Scholar}41], and an infinite element boundary is set on the right side of the model along the radial direction to eliminate the influence of the dynamic wave reflection on the subgrade deformation^{[Reference 42

Baidu Scholar}42]. By comparing the results in Table 1, it can be found that when the width of the model reaches 20 m, the mechanical response is consistent with the 5 m model （infinite element boundary）, which confirms the accuracy of the modeling method.

Table 1 Verification of mesh convergence (1 m, 30 MPa, 54 kPa·s)

Width of subgrade/m	Peak value at the center point of the loading plate
Width of subgrade/m	Stress/kPa	Displacement/μm
1	148.02	835.85
2	146.51	837.89
3	146.50	837.82
4	147.73	838.01
5	147.73	837.62
6	146.44	837.64
10	147.75	837.47
15	146.45	837.74
20	146.34	837.48
5 (infinite element boundary)	146.71	837.63

Note: The above results are calculated when the bedrock depth, the subgrade modulus and the viscosity coefficient of the subgrade are 1 m, 30 MPa and 54 kPa·s, respectively.

For the model element types, the drop weight and the loading plate are CAX4R. The element types of the finite element region and the infinite element boundary of the subgrade are CAX4R and CINAX4, respectively. As for the mesh generation rules, the drop weight, the loading plate, and the horizontal grid of the subgrade contacting the loading plate are refined with a form of uniform division; there are 10 meshes. In addition, the rest of the horizontal and vertical meshes of the subgrade are gradually coarsened from right to left and from top to bottom ^[

43], with a minimum and maximum size of 0.015 m and 0.05 m, respectively. In terms of material parameter setting, elastic parameters are used for the loading plate and the drop weight, while elastic and viscoelastic parameters are used for the subgrade. The specific conditions of the model parameters and grid properties are shown in Table 2.

Table 2 Model parameters and mesh properties

Name	Element type	Material property	Meshing type	Modulus/MPa	Poisson's ratio	Density/(kg·m^-3)
Drop weight	CAX4R	Elastic	By number	2.1×10⁶	0.25	7.5×10³
Loading plate	CAX4R	Elastic	By number	2.1×10⁶	0.25	7.5×10³
Subgrade	CAX4R CINAX4	Elastic Viscoelastic	By number By size	30,40,50,60	0.35	2.0×10³
Spring	Spring stiffness coefficient: 5.6×10⁵N/m

For the contact setting, face-to-face contact is selected for the model, and the penalty contact method is chosen for the mechanical contact formula, in which the upper contact surface is the bottom of the bearing plate and the lower surface is the upper part of the subgrade. As for the frictional behaviour between the contact surfaces, the tangential friction coefficient is 10 000, and the normal friction behaviour is hard contact. In addition, the viscoelastic property of the subgrade is modelled as a Kelvin model in terms of the Prony series （Eq. （18））, but it is inconsistent with the generalized Maxwell model built in the ABAQUS general FE software. In order to achieve the equivalent result on the existing basis, as shown in Fig. 5, an equivalent viscoelastic parameters setup method is presented, in which the modulus of the series spring is 99 times than that of the single spring. The back-calculation results show that the instantaneous elasticity generated in the FE model has no adverse effect on the accuracy of the calculation.

Fig. 5 A setup method of equivalent viscoelastic parameters

Note: The graph is the approximate representation of the Kelvin model by giving appropriate parameters to the built-in generalized Maxwell of ABAQUS. G^*is the elastic modulus of a single spring; G₁ and η are the spring modulus and the viscous element parameter in the first Maxwell model, respectively.

G_{R} (t) = G_{0} (1 - \sum_{i = 1}^{N} {\bar{g}}_{i}^{P}) + \sum_{i = 1}^{N} {\bar{g}}_{i}^{P} G_{0} e^{- t / τ_{i}^{G}}

（18）

where G_R（t） is the shear modulus of the spring; ${\bar{g}}_{i}^{p}$ , $τ_{i}^{G}$ are the parameters of the Prony series, in which ${\bar{g}}_{i}^{P} = G_{1} / G_{0}$ ; $G_{0} = G^{*} + G_{1}$ ; $τ_{i}^{G} = η / G_{1}$ .

Fig. 6 shows the stress and strain nephogram of the subgrade model, with the corresponding subgrade modulus being 30 MPa, the viscosity coefficient being 54 kPa·s and the depth of the shallow bedrock being 3 m. Fig. 6（a） is the stress nephogram under the peak load. Its stress values in the legend increase from bottom to top, and the stress is small in the center and large on both sides along the radial direction of the subgrade, which is consistent with the stress distribution law of the subgrade top under a rigid slab. In addition, the influence depth of the stress nephogram under peak load is about 62 cm, which shows the necessity of finishing the roadbed part of the subgrade. Fig. 6（b） is a strain nephogram corresponding to the peak load time. The strain values gradually decrease along the depth direction and the radial direction. The effective influence depth is about 76 cm, which is deeper than the stress nephogram under the peak load, indicating that the propagation of a dynamic wave has a diffusion effect. Fig. 6（c） is the strain nephogram under the peak displacement; its effective influence depth is 102 cm. It shows that the viscoelastic properties of the subgrade make the peak displacement lag behind the peak load, with a lag time of about 3 ms. It also indicates that when the bedrock is shallow, the reflection effect on the downward transmission of the dynamic wave is large, which affects the accuracy of the inversion modulus.

（a） Stress nephogram under the peak load

（b） Strain nephogram corresponding to the peak load time

（c） Strain nephogram under the peak displacement

Fig. 6 Stress and strain nephogram of the calculated model

3.2 Parameter sensitivity analysis

After calculating the above FE model, load and displacement time history curves at the center point of the loading plate are read and then imported into the inverse calculation program in which the forward model is the viscoelastic dynamic model considering shallow bedrock. Then, the subgrade modulus inversion calculation is carried out according to the thought of curve matching. In the mechanical model, the existence of shallow bedrock will cause stress reflection on the subgrade top, which will affect its deformation ^[

44]. The viscoelastic properties of the subgrade are characterized by the viscosity coefficient. As a viscoelastic material, the maximum displacement value of the subgrade will lag behind its maximum load value, just like the asphalt concrete mixtures ^{[Reference 45

Baidu Scholar}45].

If the existence of shallow bedrock or subgrade viscoelasticity is ignored when selecting the forward model, the theoretical time history curve will deviate from the actual condition, which will cause the back- calculation analysis of the subgrade modulus to have a larger error. The influence degree of shallow bedrock depth and subgrade viscoelasticity on the inversion modulus of the subgrade is discussed below, and the back-calculation results of the peak value method are compared with the results considering the shallow bedrock depth.

3.2.1　Analysis of the influence of shallow bedrock.

In order to explore the influence degree of the depth of shallow bedrock on the modulus back-calculation results, two subgrade forward models are considered. One is the viscoelastic dynamic model considering shallow bedrock, and the other is the viscoelastic dynamic model considering half space. The related parameters and results are shown in Table 3（a）-（c）, and the error comparison diagram and curve comparison diagrams of the two models are drawn as Fig. 7 and Fig. 8, respectively.

Table 3 (a) Comparison between the back-calculated values of the model considering shallow bedrock and half space(0 kPa∙s)

Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation with shallow bedrock			Back-calculation with half space
Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
0	0.5	30	31.98	0.997	6.6	33.37	0.978	11.2
		40	42.67	0.997	6.7	45.79	0.970	14.5
		50	55.03	0.997	10.1	58.97	0.995	17.9
		60	65.11	0.997	8.5	72.00	0.970	20.0
	1	30	33.25	0.999	10.8	35.12	0.999	17.1
		40	42.22	0.999	5.6	45.50	0.996	13.7
		50	53.27	0.995	6.5	58.90	0.994	17.8
		60	66.00	0.999	10.0	68.98	0.995	15.0
	2	30	31.29	0.998	4.3	33.44	0.999	11.5
		40	41.92	0.999	4.8	45.02	0.999	12.6
		50	52.26	0.999	4.5	55.98	0.999	12.0
		60	63.26	0.999	5.4	67.75	0.999	12.9
	3	30	31.98	0.999	6.6	31.96	0.999	6.5
		40	42.70	0.999	6.8	42.83	0.999	7.1
		50	53.98	0.999	8.0	53.78	0.999	7.6
		60	64.83	0.999	8.1	64.78	0.999	11.2
Average value				0.998	7.1		0.993	13.0

（a） Average error of modulus between the model considering bedrock and considering half space

（b） Average error of modulus at different viscoelastic

Fig. 7 Comparison charts of the average error of the modulus

coefficients under two models

;

Note: Fig. 7（a） shows the variation of the average error of the subgrade inversion modulus with the depth of shallow bedrock under the two forward models, in which the legends "back-calculation with bedrock" and "back-calculation with half space" represent the average error of the back-calculation modulus of the dynamic viscoelastic model considering the bedrock and the half space, respectively. Fig. 7（b） shows the variation of the average error of the subgrade inversion modulus of the dynamic viscoelastic model with the bedrock depth under different subgrade viscosity coefficients, in which the legends "0 kPa·s, considering bedrock" and "0 kPa·s, considering half space" are the average error of the back-calculation modulus of the dynamic viscoelastic model considering bedrock and half space, respectively.

（a） 0.5 m, 0 kPa·s, 30 MPa

（b） 1 m, 0 kPa·s, 30 MPa

（c） 2 m, 0 kPa·s, 30 MPa

（d） 3 m, 0 kPa·s, 30 MPa

Fig. 8 Curve matching diagrams

Note: Legends "Simulated curve" and "Curve considering bedrock" are the displacement time history curve read in ABAQUS and displacement time history curve iterated in MATLAB, respectively, in which the subgrade model is the viscoelastic dynamic model considering shallow bedrock. In addition, "curve considering half space" represents the displacement time history curve iterated in MATLAB when the forward model is the viscoelastic dynamic model considering half space. The parameters of the title from left to right are bedrock depth, subgrade viscosity coefficient and subgrade modulus, respectively.

Table 3 (b) Comparison between the back-calculated values of the model considering shallow bedrock and the model considering half space (54 kPa∙s)

Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation with shallow bedrock			Back-calculation with half space
Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
54	0.5	30	32.24	0.999	7.5	40.67	0.997	35.6
		40	42.32	0.999	5.8	53.52	0.997	33.8
		50	53.38	0.999	6.8	53.51	0.997	7.0
		60	64.09	0.999	6.8	80.11	0.996	33.5
	1	30	30.97	0.999	3.2	32.15	0.998	7.2
		40	41.29	0.999	3.2	44.81	0.996	12.0
		50	51.94	0.999	3.9	55.48	0.995	11.0
		60	62.81	0.999	4.7	64.15	0.995	6.9
	2	30	32.47	1.000	8.2	32.71	1.000	9.0
		40	43.74	1.000	9.4	43.32	1.000	8.3
		50	53.72	1.000	7.4	53.93	1.000	7.9
		60	65.61	1.000	9.4	65.66	1.000	9.4
	3	30	31.19	0.997	4.0	31.30	0.999	4.3
		40	44.67	0.996	11.7	41.63	0.999	4.1
		50	52.37	0.999	4.7	52.19	0.999	4.4
		60	62.87	0.999	4.8	63.28	0.999	5.5
Average value				0.999	6.3		0.998	12.5

Table 3 (c) Comparison between the back-calculated values between the model considering shallow bedrock and themodel considering half space (134 kPa∙s)

Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation with shallow bedrock			Back-calculation with half space
Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
134	0.5	30	33.60	0.999	12.0	45.22	0.999	50.7
		40	43.76	0.999	9.4	58.77	0.999	46.9
		50	54.94	0.999	9.9	73.53	0.999	47.1
		60	65.43	0.999	9.1	87.56	0.999	45.9
	1	30	32.67	0.999	8.9	37.05	0.999	23.5
		40	43.06	1.000	7.7	48.71	0.999	21.8
		50	53.86	0.999	7.7	60.87	0.998	21.7
		60	64.70	0.999	7.8	71.82	0.998	19.7
	2	30	31.67	0.999	5.6	34.99	0.999	16.6
		40	41.84	0.999	4.6	45.51	1.000	13.8
		50	50.23	0.999	0.5	56.83	1.000	13.7
		60	61.99	0.999	3.3	65.77	0.999	9.6
	3	30	32.71	0.999	9.0	33.14	0.999	10.5
		40	44.23	0.999	10.6	45.51	0.999	13.8
		50	54.16	0.999	8.3	54.26	0.999	8.5
		60	65.16	0.999	8.6	65.16	0.999	8.6
Average value				0.999	7.7		0.999	23.3

Note: "Back-calculation with shallow bedrock" and "Back-calculation with half space" represent that the subgrade forward models are the dynamic viscoelastic model considering shallow bedrock and that considering half space, respectively. For the two calculation modes, the viscosity coefficients of the subgrade are assumed to be 0 54 and 134 kPa∙s, respectively. Every viscosity coefficient corresponds to four sets of bedrock depth （0.5, 1, 2 and 3 m）. At the same time, four sets of subgrade modulus E₀ (30, 40, 50 and 60 MPa) are selected in the FE model under each bedrock depth for calculation. In addition, E₁, R² and Error are the inversion modulus of the subgrade, the curve matching correlation coefficient and the average error of the inversion modulus, respectively.

As shown in Table 3（a）-（c） and Fig. 7（a）, when the forward model of the back-calculation program is the viscoelastic dynamic model considering shallow bedrock, which is the same as the FE model, the accuracy of the subgrade back-calculation modulus is better than that of the viscoelastic dynamic model considering half space. Specifically, the average error of the former is 7.0%, while that of the latter is as high as 16.2%, and this trend is not altered with the change of the subgrade soil's viscosity coefficient. As shown in Fig. 7（b）, when the subgrade viscosity coefficients are 0, 54 and 134 kPa∙s, the average modulus errors of the former forward model are 7.1%, 6.3% and 7.7%, while those of the latter model are 13.0%, 12.5% and 23.3%, respectively. It can be seen from the above data that the existence of shallow bedrock has a significant influence on the modulus inversion of the subgrade, therefore, it is necessary to take the shallow bedrock into account in the forward model during the modulus back-calculation process. In addition, as can be seen from Fig. 7（b）, as the depth of shallow bedrock increases, the gaps in the average modulus error between the two mechanical models gradually narrows. When the bedrock depths are 0.5, 1, 2 and 3 m, the disparities between the models are 35.1%, 8.9%, 6% and -0.2%, respectively, showing that when the depth of the bedrock reaches 3 m, it has little effect on the modulus inversion results. Furthermore, it can be seen from Fig. 8 that on the basis of the dynamic viscoelastic subgrade FE model with shallow bedrock, when the forward model of the back-calculation program is consistent with it, the "curve considering bedrock" matches the "simulated curve" well. However, when the forward model is the viscoelastic dynamic model considering half space, the matching degree between "curve considering half space" and "simulated curve" becomes poor. However, the situation varies with the increase of bedrock depth. When the bedrock depth reaches 3 m, three displacement curves basically coincide, as shown in Fig. 8（d）, which shows that the bedrock at enough depth has little effect on either the subgrade deformation or the subgrade modulus inversion value.

In order to further illustrate the necessity of considering shallow bedrock in the subgrade forward model when it exists in the FE model, this paper compares the back-calculation results of the viscoelastic dynamic model considering half space with that of the peak value method. The corresponding calculation results are summarized in Table 4 and Fig. 9.

Table 4 Comparison of the back-calculated values between the dynamic viscoelastic model with half space and the peak method

Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation with half space			The peak value method
Viscosity coefficient/(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₂/MPa	Error/%
0	0.5	30	33.37	0.978	11.2	31.88	6.3
		40	45.79	0.970	14.5	43.57	8.9
		50	58.97	0.995	17.9	56.50	13.0
		60	72.00	0.970	20.0	68.05	13.4
	1	30	35.12	0.999	17.1	34.35	14.5
		40	45.50	0.996	13.7	45.28	13.2
		50	58.90	0.994	17.8	54.70	9.4
		60	68.98	0.995	15.0	65.89	9.8
	2	30	33.44	0.999	11.5	32.57	8.6
		40	45.02	0.999	12.6	43.56	8.9
		50	55.98	0.999	12.0	54.15	8.3
		60	67.75	0.999	12.9	64.52	7.5
	3	30	31.96	0.999	6.5	31.60	5.3
		40	42.83	0.999	7.1	41.83	4.6
		50	53.78	0.999	7.6	52.03	4.1
		60	64.78	0.999	11.2	62.31	3.9
Average value				0.997	13.0		8.7

Note: "Back-calculation with half space" and "The peak value method" mean that the back-calculation results in this table are based on the dynamic viscoelastic model considering shallow bedrock and the peak value method, respectively. For the calculation mode of "Back-calculation with half space", the viscosity coefficient of the subgrade is assumed to be 0 kPa·s; the viscosity coefficient corresponds to four sets of bedrock depth (0.5, 1, 2 and 3 m); and four sets of subgrade modulus E₀ (30, 40, 50 and 60 MPa) are selected in the FE model under each bedrock depth for calculation.

Fig. 9 Comparison chart of the average relative error of the modulus

Note: This figure shows the variation of the average error of subgrade inversion modulus with the depth of shallow bedrock under the two forward models, in which the legends "the peak value method" and "back-calculation with half space" represent the average error of back-calculation modulus of the two methods.

As shown in Table 4 and Fig. 9, the modulus inversion accuracy of the peak value method is better than that of the dynamic viscoelastic model considering half space. Specifically, the average error of modulus for the peak value method is 8.7%, while that of the latter reaches 13%, and this trend is not altered at different bedrock depths. When the depths of shallow bedrock in the FE model are 0.5, 1, 2 and 3 m, the average errors of the modulus for the peak value method are 10.4%, 11.7%, 8.3% and 4.4%, respectively, while those of the dynamic viscoelastic model with half space are 15.9%, 15.9%, 12.3% and 8.1%, respectively, which shows the calculation accuracy of the dynamic viscoelastic model considering half space is not as good as that of the peak method. Furthermore, the curve matching correlation coefficient under various parameter conditions is greater than 0.990 when the existence of the bedrock is ignored in the forward model, but the accuracy of the inversion modulus is not that high, which indicates the necessity of considering the shallow bedrock in the mechanical forward model from another aspect.

3.2.2　Analysis of the influence of viscoelastic property.

To explore the influence degree of viscoelasticity on the back-calculation accuracy of the subgrade model considering shallow bedrock, the above-mentioned FE subgrade models considering viscoelasticity and shallow bedrock are established, and two forward subgrade models with shallow bedrock considering viscoelasticity and elasticity are adopted in the back-calculation program for calculation. The related parameters and results are shown in Table 5（a）-（c）, and the average error comparison charts and curve comparison diagrams of the two models are drawn as Fig. 10 and Fig. 11, respectively.

Table 5 (a) Comparison between the back-calculated values of the model considering viscoelasticity and elasticity (0 kPa∙s)

Viscosity coefficient/ (kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation considered viscoelasticity			Back-calculation considered elasticity
Viscosity coefficient/ (kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
0	0.5	30	31.98	0.997	6.6	32.56	0.996	8.5
		40	42.67	0.997	6.7	43.41	0.991	8.5
		50	55.03	0.997	10.1	52.20	0.834	4.4
		60	65.11	0.997	8.5	66.13	0.996	10.2
	1	30	33.25	0.999	10.8	34.30	0.995	14.3
		40	42.22	0.999	5.6	45.40	0.999	13.5
		50	53.27	0.995	6.5	54.50	0.998	9.0
		60	66.00	0.999	10.0	66.53	0.999	10.9
	2	30	31.98	0.998	6.6	33.87	0.993	12.9
		40	42.67	0.999	6.7	45.26	0.999	13.2
		50	55.03	0.999	10.1	56.23	0.998	12.5
		60	65.11	0.999	8.5	68.13	0.998	13.6
	3	30	33.25	0.999	10.8	32.68	0.998	8.9
		40	42.22	0.999	5.6	43.54	0.998	8.9
		50	53.27	0.999	6.5	54.62	0.993	9.2
		60	66.00	0.999	10.0	65.49	0.997	9.1
Average value				0.998	8.1		0.986	10.5

（a） Average error of modulus between the model considering viscoelasticity and considering elasticity

（b） Average error of modulus at different bedrock depths under two models

Fig. 10 Comparison charts of the average error of modulus

Note: In Fig. 10（b）, legends "0.5 m V" and "0.5 m E" represent the subgrade forward model considering shallow bedrock, with the depth being 0.5 m, having viscoelasticity and elasticity, respectively, and other legends have the same meaning.

（a） 0.5 m, 0 kPa·s, 30 MPa

（b） 0.5 m, 54 kPa·s, 30 MPa

（c） 0.5 m, 134 kPa·s, 30 MPa

Fig. 11 Curve matching diagrams

Note: Legends "Simulated curve" and "Curve considering viscoelasticity" are the displacement time history curve read in ABAQUS and the displacement time history iterated in MATLAB, respectively, in which the subgrade model is the viscoelastic dynamic model considering shallow bedrock. In addition, the legend "Curve considering elasticity" represents the displacement time history curve iterated in MATLAB when the forward model is the elastic dynamic model considering shallow bedrock.

Table 5 (b) Comparison between the back-calculated values of the model considering viscoelasticity and elasticity (54 kPa∙s)

Viscosity coefficient /(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation considered viscoelasticity			Back-calculation considered elasticity
Viscosity coefficient /(kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
54	0.5	30	32.24	0.999	7.5	43.12	0.911	43.7
		40	42.32	0.999	5.8	53.73	0.921	34.3
		50	53.38	0.999	6.8	67.66	0.907	35.3
		60	64.09	0.999	6.8	76.71	0.941	27.9
	1	30	30.97	0.999	3.2	35.60	0.948	18.7
		40	41.29	0.999	3.2	46.49	0.973	16.2
		50	51.94	0.999	3.9	58.76	0.962	17.5
		60	62.81	0.999	4.7	69.36	0.972	15.6
	2	30	32.47	1.000	8.2	41.08	0.896	36.9
		40	43.74	1.000	9.4	49.77	0.957	24.4
		50	53.72	1.000	7.4	61.44	0.938	22.9
		60	65.61	1.000	9.4	72.73	0.960	21.2
	3	30	31.19	0.997	4.0	39.20	0.894	30.7
		40	44.67	0.996	11.7	45.35	0.998	13.4
		50	52.37	0.999	4.7	57.95	0.949	15.9
		60	62.87	0.999	4.8	69.68	0.960	16.1
Average value				0.999	6.3		0.943	24.4

Table 5 (c) Comparison between the back-calculated values of the model considering viscoelasticity and elasticity (134 kPa∙s)

Viscosity coefficient/ (kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	Back-calculation considered viscoelasticity			Back-calculation considered elasticity
Viscosity coefficient/ (kPa∙s)	Thickness of bedrock/m	Modulus E₀/MPa	E₁/MPa	R²	Error/%	E₁/MPa	R²	Error/%
134	0.5	30	33.60	0.999	12.0	48.72	0.592	62.4
		40	43.76	0.999	9.4	57.43	0.694	43.6
		50	54.94	0.999	9.9	68.21	0.764	36.4
		60	65.43	0.999	9.1	78.97	0.813	31.6
	1	30	32.67	0.999	8.9	50.73	0.764	69.1
		40	43.06	1.000	7.7	61.03	0.815	52.6
		50	53.86	0.999	7.7	71.41	0.870	42.8
		60	64.70	0.999	7.8	81.28	0.909	35.5
	2	30	31.67	0.999	5.6	52.64	0.678	75.5
		40	41.84	0.999	4.6	62.93	0.770	57.3
		50	50.23	0.999	0.5	73.00	0.833	46.0
		60	61.99	0.999	3.3	83.33	0.874	38.9
	3	30	32.71	0.999	9.0	49.61	0.744	65.4
		40	44.23	0.999	10.6	59.91	0.770	49.8
		50	54.16	0.999	8.3	69.63	0.834	39.3
		60	65.16	0.999	8.6	79.41	0.877	32.4
Average value				0.999	7.4		0.788	48.7

Note: "Back-calculation considering viscoelasticity" and "Back calculation considering elasticity" show that the subgrade forward models are the dynamic viscoelastic model considering shallow bedrock and the dynamic elastic model considering shallow bedrock, respectively. For the two calculation modes, the viscosity coefficients of the subgrade are assumed to be 0, 54 and 134 kPa·s, respectively. Each viscosity coefficient corresponds to four sets of bedrock depth (0.5, 1, 2 and 3 m); and four sets of subgrade modulus E₀ (30, 40, 50 and 60 MPa) are selected in the FE model under each bedrock depth for calculation.

As shown in Table 5（a）-（c）, when the forward model is consistent with the FE model, the inversion accuracy of the subgrade modulus is less affected by the viscoelastic property, and things are opposite when it is different from the FE model. Specifically, the average modulus error of the former is 7.4%, while that of the latter reaches 27.9%, and this trend is not altered by changing the viscosity coefficient of the subgrade. As shown in Fig. 10（a）, when the subgrade viscosity coefficients are 0, 54 and 134 kPa·s, the average relative modulus error of the former is 8.1%, 6.3% and 7.7%, while that of the latter is 10.5%, 24.4% and 48.7%, respectively. In addition, it can be seen from Fig. 10（b） and the above data that when the viscosity coefficient of the two forward models is 0 kPa·s, the difference in the back-calculation modulus errors is small. As the viscosity coefficient increases, the back-calculation modulus error is almost the same when the forward subgrade model is the viscoelastic dynamic model considering shallow bedrock, while the back-calculation modulus error increases almost linearly when the forward subgrade model is the elastic dynamic model considering shallow bedrock, which indicates that the viscoelasticity of the subgrade soil cannot be ignored in the forward model. Furthermore, it can be seen from Table 5（a）-（c） and Fig. 11 that on the basis of the FE dynamic viscoelastic subgrade model considering shallow bedrock, when the forward model of the back-calculation program is consistent with it, "Curve considering viscoelasticity" matches "Simulated curve" well, and its correlation coefficients R² are greater than 0.998. Moreover, with the increase of the viscosity coefficient of subgrade material, the phenomenon that the peak displacement value lags behind the peak load value becomes more obvious; the corresponding lag time are 0, 1.5 and 3 ms, respectively. However, when the forward model of the back-calculation program is the elastic viscoelastic dynamic model considering shallow bedrock, the matching degree between the "Curve considering elasticity" and "Simulated curve" becomes poor, and the value of R² is 0.986, 0.943 and 0.788, respectively. All in all, it shows that when the subgrade model considering shallow bedrock has viscoelasticity, this property must be considered in the forward model of the back-calculation program.

4 Conclusion

In this article, on the basis of the viscoelastic theory and MPGA, a new model for the elastic modulus back-calculation of subgrade considering shallow bedrock is proposed. The new model is a good supplement to the existing back-calculation method of the subgrade modulus, and is a simplification of the existing FE simulation method for shallow bedrock. To comprehensively consider shallow bedrock and subgrade viscoelasticity, the axisymmetric viscoelastic model considering shallow bedrock and the implicit analysis method are employed to derive the displacement function, and the displacement time history curves are computed according to the half-sine fitting method. Through the verification of PFWD FE simulation tests, the validity and applicability of the new way are proved. The specific results are as follows:

1） For the FE subgrade model considering shallow bedrock，when the subgrade forward model is the model considering half space, its accuracy of the inversion modulus is worse than that of the dynamic viscoelastic subgrade model considering shallow bedrock, and even worse than that of the peak value method. As the bedrock depth increases, its influence on the accuracy of the inversion modulus decreases, and the subgrade model is basically equivalent to the model considering half space when the depth reaches 3 m. Additionally, the back-calculation error of the model considering half space becomes greater when the viscosity coefficient of the subgrade soil is large.

2） When the subgrade forward model is the same as the FE subgrade model considering viscoelasticity, the error of the inversion modulus is scarcely influenced while the displacement hysteresis is greatly influenced. When the subgrade forward model is the model considering elasticity, the average relative error of the inversion modulus increases almost linearly with the increase of the viscosity coefficient. The depth of the shallow bedrock has little impact on the error of the inversion modulus when changing the viscoelastic property of the subgrade.

It is worth noting that this paper only theoretically analyzes the influence of the shallow bedrock and corresponding subgrade viscoelastic properties on the results of modulus back analysis, and further practical promotion will be explained in the follow-up study.

References

JIANG Z W, GUO X Y, LI W T, et al. Self-shrinkage behaviors of waste paper fiber reinforced cement paste considering its self-curing effect at early-ages [J]. International Journal of Polymer Science, 2016, 2016: 8690967. [Baidu Scholar]

NOIKE T, GOO I S, MATSUMOTO H, et al. Development of a new type of anaerobic digestion process equipped with the function of nitrogen removal [J]. Water Science and Technology, 2004, 49(5/6): 173-179. [Baidu Scholar]

HEIDARABADIZADEH N, GHANIZADEH A R, BEHNOOD A. Prediction of the resilient modulus of non-cohesive subgrade soils and unbound subbase materials using a hybrid support vector machine method and colliding bodies optimization algorithm [J]. Construction and Building Materials, 2021, 275: 122140. [Baidu Scholar]

LIN B, ZHANG F, FENG D C. Long-term resilient behaviour of thawed saturated silty clay under repeated cyclic loading: Experimental evidence and evolution model [J]. Road Materials and Pavement Design, 2019, 20(3): 608-622. [Baidu Scholar]

CHEN S H, DU D Y. Degradation of n-butyl xanthate using fly ash as heterogeneous Fenton-like catalyst [J]. Journal of Central South University, 2014, 21(4): 1448-1452. [Baidu Scholar]

AI X M, YI J Y, ZHAO H, et al. An empirical predictive model for the dynamic resilient modulus based on the static resilient modulus and California bearing ratio of cement- and lime-stabilised subgrade soils [J]. Road Materials and Pavement Design, 2021, 22(12): 2818-2837. [Baidu Scholar]

YAO Y S, ZHENG J L, ZHANG J H, et al. Model for predicting resilient modulus of unsaturated subgrade soils in South China [J]. KSCE Journal of Civil Engineering, 2018, 22(6): 2089-2098. [Baidu Scholar]

SAHA S, GU F, LUO X, et al. Development of a modulus of subgrade reaction model to improve slab-base interface bond sensitivity [J]. International Journal of Pavement Engineering, 2020, 21(14): 1794-1805. [Baidu Scholar]

ASLI C, FENG Z Q, PORCHER G, et al. Back-calculation of elastic modulus of soil and subgrade from portable falling weight deflectometer measurements [J]. Engineering Structures, 2012, 34: 1-7. [Baidu Scholar]

MA Q, SUN Z. Direct inversion of Young's modulus and Poisson's ratio using exact Zoeppritz equations based on double constraints [J]. Journal of Seismic Exploration, 2019, 28(2): 175-204. [Baidu Scholar]

GOPALAKRISHNAN K, KIM S. Global optimization of pavement structural parameters during back-calculation using hybrid shuffled complex evolution algorithm [J]. Journal of Computing in Civil Engineering, 2010, 24(5): 441-451. [Baidu Scholar]

SALTAN M, TERZI S, KÜÇÜKSILLE E U. Backcalculation of pavement layer moduli and Poisson's ratio using data mining [J]. Expert Systems with Applications, 2011, 38(3): 2600-2608. [Baidu Scholar]

TERZI S, SALTAN M, KÜÇÜKSILLE E U, et al. Backcalculation of pavement layer thickness using data mining [J]. Neural Computing and Applications, 2013, 23(5): 1369-1379. [Baidu Scholar]

PARK S W, PARK H M, HWANG J J. Application of genetic algorithm and finite element method for backcalculating layer moduli of flexible pavements [J]. KSCE Journal of Civil Engineering, 2010, 14(2): 183-190. [Baidu Scholar]

WANG M H, CHI S C, XIE Y F, et al. Dynamic parameters inversion analysis of rockfill materials considering interaction effects based on weak earthquakes [J]. Soil Dynamics and Earthquake Engineering, 2020, 130: 105968. [Baidu Scholar]

SHI X M, QUILTY S M, LONG T, et al. Managing airport stormwater containing deicers: Challenges and opportunities [J]. Frontiers of Structural and Civil Engineering, 2017, 11(1): 35-46. [Baidu Scholar]

ZHANG X R, OTTO F, OESER M. Pavement moduli back-calculation using artificial neural network and genetic algorithms [J]. Construction and Building Materials, 2021, 287: 123026. [Baidu Scholar]

VARMA S, EMIN KUTAY M. Backcalculation of viscoelastic and nonlinear flexible pavement layer properties from falling weight deflections [J]. International Journal of Pavement Engineering, 2016, 17(5): 388-402. [Baidu Scholar]

ZHANG J H, FAN H S, ZHANG S P, et al. Back-calculation of elastic modulus of high liquid limit clay subgrades based on viscoelastic theory and multipopulation genetic algorithm [J]. International Journal of Geomechanics, 2020, 20(10): 04020194. [Baidu Scholar]

DENG W, XU J J, ZHAO H M. An improved ant colony optimization algorithm based on hybrid strategies for scheduling problem [J]. IEEE Access, 2019, 7: 20281-20292. [Baidu Scholar]

ZAABAR I, CHATTI K, LEE H S, et al. Backcalculation of asphalt concrete modulus master curve from field-measured falling weight deflectometer data [J]. Transportation Research Record: Journal of the Transportation Research Board, 2014, 2457(1): 80-92. [Baidu Scholar]

JIAO Y L, XING X C, ZHANG P, et al. Multi-objective storage location allocation optimization and simulation analysis of automated warehouse based on multi-population genetic algorithm [J]. Concurrent Engineering, 2018, 26(4): 367-377. [Baidu Scholar]

YUN K K. Comparison of static and dynamic strains of asphalt concrete pavement from FWD tests [J]. KSCE Journal of Civil Engineering, 1997, 1(1): 39-48. [Baidu Scholar]

LIU G S, LI L, YANG X C, et al. Numerical analysis of stress distribution in backfilled stopes considering interfaces between the backfill and rock walls [J]. International Journal of Geomechanics, 2017, 17(2): 06016014. [Baidu Scholar]

UDDIN W, MEYER A, HUDSON W R. Rigid bottom considerations for nondestructive evaluation of pavements (discussion) [M]. Transportation Research Board, 1986. [Baidu Scholar]

BRIGGS R C, NAZARIAN S. Effects of unknown rigid subgrade layers on backcalculation of pavement moduli and projections of pavement performance [M]. Transportation Research Board, 1989. [Baidu Scholar]

HUDSON W R, ELKINS G E, UDDIN W, et al. Evaluation of pavement deflection measuring equipment [M]. Composite Pavements, 1987. [Baidu Scholar]

ROHDE G T, SMITH R E. Determining depth to apparent stiff layer from FWD data [R]. Federal Highway Administration, 1991. [Baidu Scholar]

SEOK Y T, SUNG S J. Effect of structural geometry of jointed concrete pavement on backcalculation using AREA method [J]. International Journal of Highway Engineering, 2007, 9(1): 39-46. [Baidu Scholar]

UZAN J, LYTTON R L, GERMANN F P. General procedure for backcalculating layer moduli [M]//Nondestructive Testing of Pavements and Backcalculation of Moduli, 1989, 1026: 217-228. [Baidu Scholar]

ZHANG J H, PENG J H, ZHANG A S, et al. Prediction of permanent deformation for subgrade soils under traffic loading in Southern China [J]. International Journal of Pavement Engineering, 2022, 23(3): 673-682. [Baidu Scholar]

SIRITHEPMONTREE H, SAPSATHIARN Y. Dynamic soil models for backcalculation of material properties from falling weight deflectometer deflection data [J]. Procedia Engineering, 2017, 189: 152-157. [Baidu Scholar]

ZHA X D. Subgrade dynamic backcalculation based on PFWD [J]. Advanced Materials Research, 2013, 723: 220-229. [Baidu Scholar]

GERRARD C M, WARDLE L J. Rational design of surface pavement layers [J]. Australian Road Research, 1980， 10(2): 3-15. [Baidu Scholar]

SENSENEY C T, KRAHENBUHL R A, MOONEY M A. Genetic algorithm to optimize layer parameters in light weight deflectometer backcalculation [J]. International Journal of Geomechanics, 2013, 13(4): 473-476. [Baidu Scholar]

PARK J, PARK M W, KIM D W, et al. Multi-population genetic algorithm for multilabel feature selection based on label complementary communication [J]. Entropy, 2020, 22(8): 876. [Baidu Scholar]

GUAN B X, ZHANG C S, NING J X. Genetic algorithm with a crossover elitist preservation mechanism for protein-ligand docking [J]. AMB Express, 2017, 7: 174. [Baidu Scholar]

HUSSAIN A, CHEEMA S A. A new selection operator for genetic algorithms that balances between premature convergence and population diversity [J]. Croatian Operational Research Review, 2020, 11(1): 107-119. [Baidu Scholar]

YU F, XU X Z. A short-term load forecasting model of natural gas based on optimized genetic algorithm and improved BP neural network [J]. Applied Energy, 2014, 134: 102-113. [Baidu Scholar]

BORJA R I, WU W H, SMITH H A. Nonlinear response of vertically oscillating rigid foundations [J]. Journal of Geotechnical Engineering, 1993, 119(5): 893-911. [Baidu Scholar]

LI M Y, WANG H. Prediction of asphalt pavement responses from FWD surface deflections using soft computing methods [J]. Journal of Transportation Engineering, Part B: Pavements, 2018, 144(2): 04018014. [Baidu Scholar]

WANG H, AL-QADI I L. Importance of nonlinear anisotropic modeling of granular base for predicting maximum viscoelastic pavement responses under moving vehicular loading [J]. Journal of Engineering Mechanics, 2013, 139(1): 29-38. [Baidu Scholar]

WANG H, LI M Y. Comparative study of asphalt pavement responses under FWD and moving vehicular loading [J]. Journal of Transportation Engineering, 2016, 142(12): 04016069. [Baidu Scholar]

NAM B H, KEE S H, YOUN H, et al. Methodology to improve the AASHTO subgrade resilient modulus equation for network-level use [J]. Journal of Transportation Engineering, 2015, 141(12): 04015027. [Baidu Scholar]

LEI L G. Backcalculation of asphalt concrete complex modulus curve by layered viscoelastic solution [D]. Michigan, USA: Michigan State University, 2011. [Baidu Scholar]