PSO-Based Multimodal Inversion of Rayleigh-Wave Dispersion Curves for the Geotechnical Characterization of an Embankment Profile

1. Introduction

Small-strain soil stiffness is one of the fundamental parameters in geotechnical engineering because it directly affects the reliability of numerical analyses, the assessment of deformations, and the interpretation of the behaviour of geotechnical structures within the working load range [1,2]. In the very small strain range, soil behaviour is commonly described by the initial shear modulus, G₀, which is particularly important in soil–structure interaction analyses, dynamic response assessment, and the definition of representative soil deformation parameters [1,2,3]. Reliable determination of small-strain stiffness is therefore important not only for general geotechnical site characterization, but also for the analysis of the behaviour of embankments and other earth structures [2,3].

The initial shear modulus is directly related to the shear-wave velocity, V_S, through the relationship G₀ = ρV_S² [1,2]. For this reason, the shear-wave velocity profile with depth is considered one of the most important indicators of stiffness variation within layered soil deposits [2,3,4]. In geotechnical practice, the determination of the Vs profile has become an important part of modern soil characterization, especially when it is necessary to investigate a larger volume of soil than can be covered solely by laboratory testing or by a limited number of point in situ tests [3,4].

Among the available approaches, methods based on the analysis of surface waves occupy a special place, as they enable efficient non-destructive testing and the assessment of changes in shear-wave velocity with depth [3,4,5]. Their main advantage lies in the dispersive behaviour of Rayleigh waves in a layered medium, where waves of different frequencies and wavelengths carry information from different depths within the investigated profile [4,5]. As a result, these methods are now widely accepted in geotechnical practice for the characterization of natural deposits, embankments, transport infrastructure, and other shallow to intermediate-depth soil profiles [3,4,5,6]. Within this group of methods, multichannel analysis of surface waves is of particular interest because, with appropriate processing, it also enables the identification of higher modes of the dispersion curve [4,5,6].

However, the transition from the measured dispersion curve to the shear-wave velocity profile is not unique. The non-uniqueness of the inversion process is one of the key challenges in the interpretation of surface-wave data, because different layering models may lead to very similar agreement between theoretical and experimental dispersion curves [5,6,7,8]. This problem is particularly pronounced under more complex layering conditions, such as velocity inversions, marked stiffness contrasts, or the presence of thin layers [6,7,8]. For this reason, the literature increasingly emphasizes the importance of robust global search procedures, as well as approaches that include higher modes of the dispersion curve in the inversion process [7,8,9,10]. The inclusion of higher modes can increase the sensitivity of the inversion procedure to individual model parameters, improve the resolution of the investigated profile, and contribute to a more stable determination of soil stiffness properties [8,9,10].

In this context, multimodal inversion and global optimization algorithms are of particular interest for geotechnical applications. Previous studies have shown that the inclusion of higher modes can improve the reliability of V_S-profile estimation compared with approaches based solely on the fundamental mode [8,9]. At the same time, particle swarm optimization algorithms have proven to be suitable for nonlinear and multimodal dispersion-curve inversion problems, enabling efficient global exploration of the space of possible solutions [9,10]. This makes them suitable for the development of practical procedures for the geotechnical characterization of layered soil profiles.

Based on the above, the aim of this study is to develop a procedure for determining small-strain soil stiffness properties through the measurement and inversion of Rayleigh-wave dispersion curves in layered soil. The basic hypothesis is that the reliability of shear-wave velocity profile determination can be improved by incorporating higher modes into the inversion process. Accordingly, this paper presents a multimodal inversion procedure based on the particle swarm optimization algorithm, its verification on synthetic soil profiles, and its validation on an embankment field profile through comparison with seismic CPT results.

2. Methodology

2.1. Surface-Wave Method Framework

Methods based on the analysis of surface waves are among the most commonly applied non-destructive geophysical techniques for determining small-strain soil stiffness. Their application is based on the fact that, in the very small strain range, the initial shear modulus is directly related to shear-wave velocity. Therefore, determining the shear-wave velocity profile with depth enables the assessment of stiffness changes within layered soil deposits [4,5].

The basis of all these methods is the dispersive behaviour of Rayleigh waves in a vertically heterogeneous medium. Waves of different frequencies, or wavelengths, sample different depth zones, so the variation of phase velocity with frequency contains information on the distribution of stiffness with depth. In a horizontally layered medium, the propagation of surface waves is generally a multimodal phenomenon, meaning that, in addition to the fundamental mode, higher vibration modes may also exist. This is of particular importance for the interpretation of experimental dispersion curves and the subsequent inversion procedure [4,5,15].

In geotechnical practice, the most commonly used methods are continuous surface-wave analysis, or the continuous surface-wave system (CSW/CSWS), spectral analysis of surface waves (SASW), and multichannel analysis of surface waves (MASW). All of these methods are based on the generation and recording of Rayleigh waves at the ground surface, but they differ in terms of the source type, the number of receivers, and the signal-processing procedure [3,4,5,11,12,13,14,15].

CSW or CSWS uses a controlled sinusoidal source, most commonly a vibrator, to generate excitation over a prescribed frequency range. Such an approach enables good control of the frequency content of the excitation and allows the investigation to be targeted toward the depth range of interest. In addition, owing to the controlled excitation, the useful signal can be extracted more effectively under conditions of increased noise, which is one of the reasons why this method has found application in geotechnical investigations aimed at assessing soil stiffness [3,11].

SASW is an older and widely used active surface-wave method developed for determining shear-wave velocity profiles solely from the ground surface. In its classical configuration, it uses two receivers, and the analysis is based on the phase relationship between signals recorded at a known spacing between the receivers. The main advantage of the method is its simple and economical field implementation, whereas its sensitivity to the signal-to-noise ratio and to the quality of phase-difference estimation is one of the reasons why more robust multichannel approaches were later developed [12,15].

MASW is based on multichannel recording, most commonly with uniformly spaced receivers arranged in a linear geometry and an impulsive source. The main advantage of this method is that multichannel acquisition and wavefield transformation allow a clearer representation of dispersion trends, more effective separation of coherent noise, and easier identification of the fundamental and higher modes. For this reason, MASW enables not only the determination of one-dimensional shear-wave velocity profiles, but also the construction of two-dimensional representations of soil stiffness variation along the profile, with relatively efficient and economical field implementation [13,14,15].

For the purposes of this study, particular attention is focused on the MASW method. The reason is that, compared with CSWS and SASW, it is better suited for the identification of higher modes of the dispersion curve, which is of direct importance for the multimodal inversion procedure developed in this work. Accordingly, in the following sections, surface-wave methods are considered primarily as the experimental basis for determining Rayleigh-wave dispersion curves, while the emphasis of the subsequent analysis is placed on the theoretical calculation of dispersion curves and their multimodal inversion [5,13,14,15].

2.2. Calculation of Theoretical Dispersion Curves

For the theoretical calculation of dispersion curves, the soil is considered as a horizontally layered elastic medium consisting of a finite number of homogeneous, isotropic, and linearly elastic layers overlying an elastic half-space [16,17,18,19,20]. Each layer is described by its thickness $h_i$, density ${\rho }_i$, shear-wave velocity $V_{Si}$, and compression-wave velocity $V_{Pi}$. An example of such a model, for the case of two layers overlying an elastic half-space, is shown in Figure 1.

The objective of the theoretical calculation is to determine, for a given soil profile, the relationship between frequency, or wavelength, and the phase velocity of Rayleigh waves. The resulting theoretical dispersion curve provides the basis for the subsequent inversion procedure. The phase velocity of a Rayleigh wave can be expressed as

$$c_R={\displaystyle \frac{\omega }{k}}=f\lambda$$

(1)

where $c_R$ is the phase velocity of the Rayleigh wave, $\omega$ is the angular frequency, $k$ is the wavenumber, $f$ is the frequency, and $\lambda$ is the wavelength. In a layered medium, the phase velocity is not constant, but depends on frequency or wavelength, which gives rise to the dispersion curve.

In this paper, the theoretical calculation is based on the formulation using the dynamic stiffness matrix of a horizontally layered medium [18,19]. Within this approach, the condition for free-wave propagation can be written as

$$\mathbf{K}(\omega ,c_R)\mathbf{u}=0$$

(2)

where $\mathbf{K}$is the global dynamic stiffness matrix of the system and $\mathbf{u}$is the displacement vector. A non-trivial displacement solution exists only when the matrix $\mathbf{K}$ is singular, that is, when the system satisfies the free-wave condition for the considered combination of frequency and phase velocity.

In the literature, two numerical approaches are most commonly considered for identifying the modes of the dispersion curve. The first is the DET method, in which the modes are identified by tracking the zeros of the determinant of the global stiffness matrix,

$$\mathrm{d}\mathrm{e}\mathrm{t}\mathbf{K}(\omega ,c_R)=0$$

(3)

The second is the SAE method, in which matrix singularity is assessed through the smallest absolute eigenvalue of the global stiffness matrix. In this case, a mode is identified for the phase velocity value satisfying

$$\mathrm{m}\mathrm{i}\mathrm{n}\left(\mid {\mu }_n\left(\mathbf{K}\right)\mid \right)\to 0$$

(4)

where ${\mu }_n\left(\mathbf{K}\right)$ is an eigenvalue of the global dynamic stiffness matrix. Both approaches are based on the same physical condition, namely system singularity, and yield the same results, but they differ in terms of numerical implementation and numerical stability [18,19,21].

In this way, for each considered frequency or wavelength, the points of the fundamental and higher modes of the dispersion curve can be determined, i.e., multiple possible surface-wave velocities for the same wavelength. Figure 2 shows an example of a theoretical multimodal dispersion curve, including the fundamental mode and the first two higher modes of Rayleigh waves.

In this paper, the SAE approach is adopted because it provides smooth and numerically stable theoretical dispersion curves and is suitable for the determination of higher modes, which is particularly important for the multimodal inversion procedure. Therefore, in the remainder of the paper, the theoretical dispersion curves are calculated using the SAE method [19,21].

Since the aim of this paper is multimodal inversion, the theoretical calculations consider not only the fundamental mode, but also the higher modes of Rayleigh waves. This provides an appropriate theoretical basis for the subsequent inversion procedure, in which experimental and theoretical dispersion curves are compared across multiple modes in order to achieve a more reliable determination of shear-wave velocity profiles in layered soil.

2.3. PSO-Based Inversion Procedure

The inversion of dispersion curves is formulated as a global optimization problem in which the objective is to determine the set of soil model parameters that provides the best agreement between the theoretical and experimental dispersion curves. Since this is a nonlinear and multimodal problem, the inversion process is carried out using the particle swarm optimization (PSO) algorithm, which belongs to the class of population-based evolutionary algorithms [22,23]. In this approach, each particle represents one possible solution, i.e., one set of soil model parameters, and the optimization is performed through an iterative search of the solution space.

The vector of a single particle is defined by the soil model parameters to be determined through the inversion procedure. These parameters generally include the thicknesses of individual layers and the corresponding shear-wave velocities, while the theoretical dispersion curves for each trial model are calculated using the procedure described in the previous subsection. For each particle, a theoretical multimodal dispersion curve is obtained from the prescribed parameter vector, and its deviation from the experimental dispersion curve is then evaluated. The optimization process is directed toward minimizing this mismatch, i.e., the objective function.

The objective function may be written as the sum of normalized differences between the theoretical and experimental phase velocities for all included modes and all considered points of the dispersion curve:

$$F(\mathbf{m})={\displaystyle \frac{1}{N}}{\sum }_{j=1}^M{\sum }_{i=1}^{N_j}{\left({\displaystyle \frac{c_{ij}^{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}}\left(\mathbf{m}\right)-c_{ij}^{\mathrm{e}\mathrm{x}\mathrm{p}}}{c_{ij}^{\mathrm{e}\mathrm{x}\mathrm{p}}}}\right)}^2$$

(5)

where $\mathbf{m}$ is the vector of unknown soil model parameters, $M$ is the number of considered modes, $N_j$is the number of points on the dispersion curve for the $j$-th mode, $N=\sum N_j$ is the total number of considered points, and $c_{ij}^{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}}$ and $c_{ij}^{\mathrm{e}\mathrm{x}\mathrm{p}}$ are the theoretical and experimental phase velocities for the $i$-th point and the $j$-th mode, respectively. In this way, the optimization problem is reduced to finding the set of model parameters that minimizes the value of the function $F$.

The basic principle of the particle swarm optimization algorithm is that each particle changes its velocity and position in the search space during the iterations, taking into account its own best solution found so far and the best solution found by the entire swarm. The basic principle of particle movement in the search space is illustrated in Figure 3.

The updating of the velocity and position of the $i$-th particle at step $k+1$ can be written as [22,23]

$${\mathbf{v}}_i^{k+1}=w{\mathbf{v}}_i^k+c_1{r}_1\left({\mathbf{p}}_i^k−{\mathbf{x}}_i^k\right)+c_2{r}_2\left({\mathbf{g}}^k−{\mathbf{x}}_i^k\right)$$

(6)

$${\mathbf{x}}_i^{k+1}={\mathbf{x}}_i^k+{\mathbf{v}}_i^{k+1}$$

(7)

where ${\mathbf{x}}_i^k$ is the position of the particle, ${\mathbf{v}}_i^k$ is its velocity, $w$ is the inertia factor, $c_1$ is the cognitive factor, $c_2$is the social factor, $r_1$and $r_2$are random numbers in the interval from 0 to 1, ${\mathbf{p}}_i^k$is the best position found so far by the particle under consideration, and ${\mathbf{g}}^k$is the global best solution found by the entire swarm.

At the beginning of the optimization, the initial position and velocity of each particle are randomly assigned within the predefined lower and upper bounds of the parameters. The value of the objective function is then calculated for all particles, their best individual positions are determined, and the global best solution of the swarm is identified. The procedure requires the definition of lower and upper parameter bounds, the number of particles, the values of the cognitive and social factors, the maximum number of iterations, and the stopping criteria.

After each iteration, it is checked whether the new particle positions fall outside the prescribed search-space bounds. If this occurs, their values are corrected to remain within the allowable range. The objective function is then recalculated, the best individual positions are updated, and a new global best solution is determined. The procedure is repeated until the prescribed maximum number of iterations is reached or until the change in the best particle position, or the change in the optimal value of the objective function, becomes smaller than a predefined threshold.

In addition to the basic PSO approach, a modified variant with periodic resetting of a portion of the particles is also applied. Compared with the classical algorithm, this variant allows a certain fraction of particles at a given iteration step to re-enter the search process from a newly assigned random position. This increases the possibility of escaping local minima and improves the exploration of the solution space. As the optimization progresses, the number of particles whose initial positions are reset is gradually reduced, so that in the final stages of the iteration this number tends toward a single particle [24,25,26].

Since the objective is multimodal inversion, the objective function is not based solely on the fundamental mode, but also includes the higher modes of the dispersion curve. This increases the amount of information incorporated into the optimization process and reduces the possibility that different soil profiles may provide equally good agreement in only one mode. Such an approach is particularly important under more complex layering conditions, where the inclusion of higher modes may lead to a more reliable determination of the shear-wave velocity profile [9,10].