Multilevel Genetic Algorithm–Based Acoustic–Elastodynamic Imaging of Coupled Fluid–Solid Media to Detect an Underground Cavity

Detecting underground cavities (e.g., karstic cavities, caves, tunnels, and so forth) is a challenging task for geotechnical engineering projects due to the geological and hydrogeological complexity of the subsurface environment. Geological hazards, such as collapsing soil or urban ground collapse due to subsurface voids, could induce significant damage to infrastructure. Such a hazard is one of the major issues for land planning, infrastructure operation and maintenance, and disaster management. In addition, cavity evolution that occurs due to hydrogeological process may cause sinkhole collapse. For example, Florida’s karst environment involves active groundwater recharge to the Floridan Aquifer that erodes overburden soils; thus, a cavity could grow gradually, leading to a sudden collapse (cover-collapse type) or gradual subsidence (cover-subsidence type) (Fig. 1) (Beck 1986; Tihansky 1999; Xiao et al. 2016; Kim et al. 2020; Nam et al. 2020). These large underground water-filled cavities hidden below building and transportation infrastructure should be predetected so that prevention and mitigation measures are applied before catastrophic collapse or excessive ground settlement takes place.

Cone penetration tests (CPTs) and standard penetration tests (SPTs) have been used to detect and characterize subsurface cavities in a cover soil layer, referred to as a raveled zone in karst areas (Nam et al. 2018; Nam and Shamet 2020). Although SPTs and CPTs provide a continuous subsurface profile—revealing soil type, resistance, and stratigraphy—these invasive methods collect only pointwise profiles of the subsurface environment. It is time-consuming and labor-intensive to apply CPTs and SPTs to wide areas multiple times.

Numerous geophysical studies over the last decades have introduced nondestructive evaluation (NDE) methods for detecting and estimating the geometry of subsurface cavities. For example, ground-penetrating radar (GPR), microgravimetry, and resistivity imaging and their combinations can characterize an underground cavity based on its size and depth position (Fehdi et al. 2014; Kiflu et al. 2016). GPR is inefficient in highly heterogeneous media—such as backfills, moist clays with high soil conductivity, and saturated soils below the groundwater table. The microgravimetric method is effective in detecting shallow cavities (Butler 1984; Bishop et al. 1997), but its density contrast cannot be obtained clearly if the size of a void is relatively small compared with its depth. An electrical resistivity method has been applied to detect air-filled or water-filled cavities (Van Schoor 2002; Coşkun 2012), but it requires a large areal space for surveying and cannot detect small-scale irregularities in the geologic interfaces because it measures average resistivity values and is easily ruined by noise sources, such as piping, power lines, and house structures.

On the other hand, elastodynamic imaging has been used widely in site characterization (Brown et al. 2002; Kallivokas et al. 2013). It uses elastodynamic waves that are initiated by wave sources on the ground surface and/or borehole sources, and then are reflected or refracted due to material heterogeneity within a medium. The resulting wave motions can be measured on the ground surface or at borehole seismic arrays. Inverse modeling, associated with elastic wave propagation analyses, provides promising results for identifying the material properties of the media. Inverse modeling—employing a partial differential equation (PDE)-constrained optimization method and the finite-element method (FEM) or spectral-element method (SEM)—has been utilized to infer the spatial distribution of material properties of host media at fine resolution (Kang and Kallivokas 2011; Fathi et al. 2015a). For example the geotechnical site characterization methodologies have been investigated in a soil domain that is truncated by perfectly matched layers (PML), in which waves decay and are not reflected off surrounding boundaries (Fathi et al. 2015b), using state-adjoint-control equation-based full-waveform inversion (FWI) approaches (Kang and Kallivokas 2010, 2011; Kallivokas et al. 2013; Fathi et al. 2015a, 2016; Kucukcoban et al. 2019). In particular, there has been progressive development of the material inversion method for site characterization from a simplified one-dimensional (1D) setting to a full three-dimensional (3D) setting followed by field experimental validation. Kang and Kallivokas (2010) investigated a new material inversion method for identifying the vertical distribution of the shear wave speed within a PML-truncated 1D soil column with a hypothesis of horizontal layering without considering geometrical damping. Kang and Kallivokas (2011) further developed their material inversion method for estimating the spatial distribution of the shear wave speed in a two-dimensional (2D) scalar wave setting at fine resolution without considering the realistic behaviors of elastic waves (i.e., vector waves). Kallivokas et al. (2013) continued the investigation of the new material inversion method to estimate the spatial distributions of both P- and S-wave speeds within a 2D linear, elastic, undamped PML-truncated solid setting in consideration of vector wave behaviors. They used an assumption that the material properties of the soils are uniform in the antiplane direction of the considered 2D plane. They validated their numerical method using field experimental data that used a line loading–like application of wave sources, which replicated the 2D plane-strain setting of the computational study. Fathi et al. (2015a) developed the material inversion method for inverting the spatial distribution of both P- and S-wave speeds in a 3D linear, elastic, undamped PML-truncated solid. Fathi et al. (2016) validated the result of the inverse modeling of Fathi et al. (2015a) using field experimental data in a 3D setting. In the experiment, a point wave source was employed on the ground surface, and its corresponding wave responses were measured on the surface and used as measurement data in the inverse simulations. Fathi et al. (2016) compared the inverted spatial distribution of the wave speeds with that of a reference solution from a conventional method, such as spectral analysis of surface waves (SASW).

In addition, the Gauss–Newton-based FWI method, which is based on a gradient vector and a Hessian matrix, has been studied for characterizing the material profiles in a truncated two-dimensional solid domain (Tran and McVay 2012; Pakravan et al. 2016) for the application of site characterization. The PDE-constrained optimization also has been used in consideration of the boundary-element method (BEM), the computational efficiency of which is much greater than that of the FEM due to the reduction of the dimensionality, to detect the geometries of wave-scattering objects in host media. There have been studies of inverse scattering algorithms using BEM wave solvers, hinging on the moving boundary concept and the total derivative (Petryk and Mroz 1986), which allows computing the derivative of an objective functional with respect to the geometry variables of scattering objects (Guzina et al. 2003; Jeong et al. 2009). The extended finite-element method (XFEM) has been used as a wave solver in studies to identify cracks or air-filled voids in solid media because it can avoid expensive remeshing processes in inverse iterations while using the flexibility of the FEM for heterogeneous materials (Jung et al. 2013). In addition, a frequency-continuation scheme was devised to help the inversion solver address the multiplicity of solutions of subsurface imaging problems (Bunks et al. 1995; Kang and Kallivokas 2010; Fathi et al. 2015a). That is, when the convergence rate of the inversion solver is decreased due to a local minimum of an objective functional that is composed of measurement data corresponding to a given dominant frequency of excitation, another set of measurement data that are induced by excitation of a different dominant frequency is employed. Such a frequency shift can change the curvature of the objective functional so that a minimizer can escape a local minimum. In general, one may consider increasing the frequency from one inversion set to another because it has been reported that using a lower frequency leads to the inversion result of a lower resolution (due to a larger wavelength); then one can fine-tune it using a higher frequency (due to a smaller wavelength) (Bunks et al. 1995; Kang and Kallivokas 2010; Fathi et al. 2015a).

Despite such recent developments in elastodynamic wave imaging techniques, few papers have demonstrated the feasibility of identifying the properties of a solid system that includes fluid-filled voids. Finite-difference method (FDM)-based FWI approaches have been presented, which employ all solid elements in an entire computational domain and detect the areas with smaller values of shear wave speeds $V_s$ (e.g., about $100\mathrm{m}/\mathrm{s}$) than those of typical soils and rocks to indicate air-filled voids (Tran and McVay 2012; Tran et al. 2013; Mirzanejad et al. 2020). However, those works have not explicitly modeled the interfaces between air voids and the solid domain. Thus, the authors are concerned that, first, modeling a fluid domain as solid elements with a small value of $V_s$ could introduce numerical error in forward wave solutions. Second, using a small value of $V_s$ requires a small solid element for wave simulations. However, the aforementioned FDM works did not use such a small element to model a void, so they could suffer from the additional error of wave solutions. To address this issue, the authors suggest explicitly using fluid elements to model voids filled with water (or air) so that accurate fluid–solid coupling (Everstine 1997; Lloyd et al. 2016) can be incorporated into the wave solver. Then, using such a wave solver, inverse modeling could identify accurately the interfaces between fluid and solids as well as the material profiles of solids. As a prototype work of the suggested method, this paper investigates the feasibility of estimating the interfaces between fluid and solid layers and the material properties of solid layers in a multilayered system using acoustic and elastodynamic waves generated from the ground surface. This investigation proved the feasibility of detecting a water-filled void and identifying its location and geometry in soils, due to which overlaying ground surface potentially could sink.

This work used the FEM to study the wave responses of the coupled fluid–solid media in a prototype one-dimensional setting. The solid formation of a semi-infinite extent was truncated by using an absorbing boundary, and a proper fluid–solid interface condition was considered. The inverse problem was cast into a minimization problem, in which the value of an error between the synthetic measured wave response due to target profiles and the computed counterpart due to guessed profiles was minimized. The solution multiplicity of the minimization process was addressed using a multilevel genetic algorithm (GA) in this work. Under this new scheme, the authors suggest the following steps: (1) the optimizer identifies an inversion solution around a local minimum after a number of GA iterations in the first-level GA; (2) a misfit function then can be r-calculated for a different measurement signal induced by a pulse signal of a different central frequency in the next-level GA; (3) by doing so, the GA optimizer is able to escape the local minimum of the previous-level GA; and (4) the same technique is employed to address the solution multiplicity of the next-level GA. The numerical results show that the numerical simulations using the multiple-frequency-level GA accurately detects the interfaces between fluid and solid layers and the material properties of solid layers. The authors also demonstrate that a preliminary GA inversion can determine whether the subsurface media include a fluid layer. Specifically, the preliminary inversion can indicate the existence of a fluid layer by detecting a very large value of Young’s modulus (with a value similar to that of the bulk modulus of water) and a mass density with a value similar to that of water. When the existence of a fluid layer is detected, a primary GA inversion using the fluid–solid model can estimate further the locations of a fluid layer and the material properties of surrounding solids.

This study identified the depths of the interfaces between a fluid layer and its surrounding solid layers and the stiffness of the solid layers using a dynamic test, which generates elastodynamic waves into the subsurface media from the top surface. A 1D fluid–solid model, which consists of three different layers, i.e., two solid layers and one fluid layer (Fig. 2), was considered. The governing differential equations for the compressional waves in the three layers arewhere $L_{s1}$, $L_f$, and $L_{s2}$ = depth of interface between top solid layer and fluid layer, depth of interface between fluid and bottom solid layer, and depth of truncation wave absorbing boundary, respectively; $u_{s1}(x,t)$ and $u_{s2}(x,t)$ = displacement fields of wave motions of solid particles in top and bottom solid layers; and $P$ = fluid pressure field of wave motions in fluid. The wave speeds of the three layers are given by $c_{s1}=\sqrt{E_{s1}/{\rho }_{s1}}$, $c_f=\sqrt{{\kappa }_f/{\rho }_f}$, and $c_{s2}=\sqrt{E_{s2}/{\rho }_{s2}}$, respectively, where $E$, $\rho$, and $\kappa$ are the modulus of elasticity, the density, and the bulk modulus. The problem of interest also includes the following boundary conditions:

A pulse signal $f(t)$ is applied at the top of the solid layer at $x=0$ [Eq. (4)]. The depth of the bottom solid layer is truncated by using the absorbing boundary condition [Eq. (5)].

The solid and fluid layers are coupled via the following interface conditions:

The preceding interface conditions originate from the dynamic interface condition between a solid and fluid in a multidimensional setting, presented by Everstine (1997), which is composed of the following. First, the effect of fluid pressure on the structure is imposed as a load. This corresponds to Eqs. (6) and (8). Second, the effect of structural motion on the fluid should be modeled as $\nabla P·\mathbf{n}=-{\rho }_f[({\partial }^2{\mathbf{u}}_s)/\partial t^2]·\mathbf{n}$, where $\mathbf{n}$ is the normal vector at the interface. This corresponds to Eqs. (7) and (9).

The system initially is at rest

To derive the weak form, the governing wave equations in the strong form Eqs. (1)–(3) are multiplied by the test functions $v_1(x)$, $w(x)$, and $v_2(x)$. By integrating them by parts and using the boundary and interface conditions of the strong form, the following weak forms are obtained:

Then the functions are approximated aswhere $\mathit{\varphi }(x)$, $\mathbf{\Psi }(x)$, and $\mathbf{\Omega }(x)$ denote vectors of global basis functions constructed by local shape functions; and ${\mathbf{u}}_{s1}$, $\mathbf{P}$, and ${\mathbf{u}}_{s2}$ are vectors of nodal solutions. Introducing the finite-element approximations Eqs. (15)–(17) into the weak forms Eqs. (12)–(14) provides the following discrete form:where the matrixes are defined asand the force vector is $\mathbf{f}=\mathit{\varphi }(0)f(t)$. The discrete Eqs. (18)–(20) lead to the following coupled discrete form:which, at a discrete time $t_n$, can be writtenwhere ${\mathbf{q}}_n={[{\mathbf{u}}_{s1}^T,{\mathbf{P}}^T,{\mathbf{u}}_{s2}^T]}^T$ is the vector containing the nodal solutions of solid displacements and acoustic pressures at the $n$th time step; and ${\dot{\mathbf{q}}}_n$ and ${\ddot{\mathbf{q}}}_n$ denote the first- and second-order derivatives of ${\mathbf{q}}_n$ with respect to time. The second-order ordinary differential [Eq. (28)] is solved using Newmark’s unconditionally stable time-integration scheme (Newmark 1959).

This study identified the four control parameters ($E_{s1}$, $E_{s2}$, $L_{s1}$, and $L_f$) by using measured wave responses initiated by a known excitation signal $f(t)$. The values of the other variables ($L_{s2}$, ${\rho }_{s1}$, ${\rho }_{s2}$, ${\kappa }_f$, and ${\rho }_f$) were set to be known during the inversion process, and the excitation signal $f(t)$ also was known. This problem was formulated into a minimization problem to identify estimated control parameters that could lead to the minimum of an objective functionalwhere $u_m(0,t)$ = wave signal recorded by sensor at $x=0$ due to target profile of control parameters; and $u(0,t)$ = computed counterpart at the same measurement location due to guessed profile. In this computational study, $u_m(0,t)$ was obtained synthetically by running the presented FEM solver using a target profile as an input. To prevent an inverse crime from taking place, smaller values of element sizes were used when $u_m(0,t)$ was computed than when $u(0,t)$ was computed.

This work used a GA, which is a heuristic algorithm, to reconstruct the control parameters that are the fittest to the objective of the problem. In each generation of GA there are individuals, each of which has a set of different values of control parameters. In this work, the GA runs the FEM wave solver for the estimated values of control parameters of each individual to compute $u(0,t)$ and evaluate the misfit [Eq. (29)]. Then the GA evaluates the fitness of each individual using the value of its misfit. The GA then weeds out less-fitting individuals when it creates the next generation of individuals, each of which is characterized by its own set of estimated parameters. In this process, by virtue of the mutation process, in which the GA learns the better-fitting characteristics of control parameters, the GA gives rise to new individuals in the next generation. The GA calculates the sensitivity of the fitness function with respect to the changes of the values of the control parameters and determines how the next generation of individuals is created. Namely, the GA repeats the tasks of weeding out, selecting the best individual, and mutating at each generation. Therefore, the suggested inverse modeling solver could identify targeted values of the parameters using the iterative process of the GA, in which the values of the estimated parameters of the best-fit individual evolve over the progress of generations.

In the overall GA process in this problem, constraints were imposed on the allowable ranges of each control parameter. For example, $L_{s1}$ (the depth of the upper face of a fluid layer) was set to be always smaller than $L_f$ (the depth of the lower face) so that the thickness of the fluid layer always was positive. In addition, the values of elastic moduli $E_{s1}$ and $E_{s2}$ were set to be always positive.

It is known that in the minimization problem of a small number (e.g., less than 20) of control parameters, a GA is likely to find a set of control parameters that are close to the global minimum (Jeong et al. 2017; Guidio and Jeong 2021). Because there are only four control parameters, the authors originally hypothesized that the GA is a suitable solution approach to finding the global minimum solution of this inverse problem. However, the numerical experiments showed that the GA suffers from the solution multiplicity, and, thus, the authors tested a new multiple frequency–level GA approach. Namely, after the first level of the GA is finished, the presented inversion solver uses the reconstructed values from the first level to define the upper and lower limits of the control parameters in the next GA level, which uses a pulse signal of a different central frequency from its predecessor, creating a new synthetic $u_m$. In the final-level GA, the final best-fit estimated control parameters are obtained as the inversion solution.

This section presents a set of numerical experiments studying the performance of the presented multilevel GA–based parameter-estimation method. In the first two examples, three targeted fluid–solid models, which differed from each other in terms of $L_{s1}$, $L_f$, and $L_{s2}$, were considered as follows. Model 1 had a total length $L_{s2}$ of 60 m and its fluid was located between $x$ of $L_{s1}=20\mathrm{m}$ and $x$ of $L_f=25\mathrm{m}$. The corresponding geometry parameters of Model 2 were $L_{s2}$ of 60 m, $L_{s1}$ of 40 m, and $L_f$ of 45 m; and those of Model 3 were $L_{s2}$ of 120 m, $L_{s1}$ of 50 m, and $L_f$ of 80 m. Parameters $L_{s1}$ and $L_f$ are unknown parameters in the inversion. In all the models, the mass density ($\rho$) of the solid layers was $\mathrm{2,000}\mathrm{kg}/{\mathrm{m}}^3$; the Young’s moduli were $E_{s1}=2\times {10}^8\mathrm{Pa}$ and $E_{s2}=5\times {10}^8\mathrm{Pa}$; and the bulk modulus ($k_f$) and mass density (${\rho }_f$) of the fluid layer were $2.34\times {10}^9\mathrm{N}/{\mathrm{m}}^2$ and $\mathrm{1,021}\mathrm{kg}/{\mathrm{m}}^3$, respectively. Whereas $E_{s1}$ and $E_{s2}$ were set to be unknown, $\rho$, $k_f$, and ${\rho }_f$ are set to be known in the inversion.

In the presented numerical experiments, a three-level GA–based inversion method was used with a frequency-continuation scheme. Namely, in each GA set, a dynamic force with a different frequency content was used as follows:

The computational procedure of the presented multilevel GA–based inversion method is summarized in Algorithm 1.

To avoid an inverse crime, to compute $u_m$ induced by targeted control parameters, the fluid–solid domain was discretized using an element size of 0.1 m, whereas an element size of 0.2 m was used to compute $u$ due to estimated control parameters. In the forward and inverse modeling, the time step was 0.001 s, and the total observation duration $T$ was set to 1 s.

Three examples of numerical experiments are presented. The first example tested the inversion performance of the three-level GA approach by using a dynamic force of a Ricker pulse signal with a central frequency of 5, 10, and 15 Hz, respectively, in each level. In the second example, the inversion performance was examined by using a decreasing frequency-continuation counterpart of 15, 10, and 5 Hz. The last example investigated the utilization of preliminary inversion using all-solid layers for detecting the presence of a fluid–filled cavity by identifying $E$ and $\rho$ of all the solid layers, among which one layer had a very large value of $E$ (with a value similar to that of the bulk modulus of water) and $\rho$ of the value of the mass density of water. Once the preliminary inversion is completed, the optimizer uses the presented three-level GA approach based on the fluid–solid wave model as a primary inversion.

To analyze the inversion results, the error between each target control parameter ($L_{s1}$, $L_f$, $E_{s1}$, and $E_{s2}$) and its estimated counterpart of the fittest individual at each generation is calculated as

The average error of all the control parameters of the best-fit individual at each generation is calculated aswhere ${\mathcal{E}}_k$ = error [Eq. (30)] of reconstruction of $k$th control parameter.

This work presents a new multilevel GA–based inversion method combined with a frequency-continuation scheme. Three levels of GA are used, and in each level, a pulse signal with a different dominant frequency is utilized to tackle the solution multiplicity of the presented inverse problem. Numerical results show the following findings about the proposed inversion method. First, it is feasible to estimate the depths of solid–fluid interfaces and Young’s moduli of the surrounding solid layers by employing an acoustic-elastodynamic wave model and a measured dynamic response at a sensor on the top surface. Second, combining the multilevel GA–based optimizer with a frequency-continuation scheme improves the performance of identifying targeted control parameters. Third, the frequency-continuation scheme in which the frequency is decreased progressively for each GA level is as effective as the conventional frequency-continuation scheme (i.e., increasing the frequency for each level) in the presented 1D setting. Lastly, a preliminary inversion using a solid-only model can address the uncertain presence of a fluid layer in a domain by determining that one of the layers has a very large value of Young’s modulus (with a value similar to that of the bulk modulus of water) and a mass density similar to that of water. After the preliminary inversion, a primary inversion—based on a multiple frequency–level GA using a fluid–solid model—can adjust the fluid layer’s depth information and the material properties of the soil layers. The authors also discussed how to extend the presented method in realistic 2D or 3D settings. Such an extension effectively may detect water-filled underground cavities and estimate the risk of potential urban cave-in and karst sinkholes. In typical engineering projects, engineers do not rely on a single geophysical method no matter how advanced its theories or computation capacities are. Thus, the presented method and its extension can be an important supplement to existing testing methods rather than a standalone solution in practice.

Some or all data, models, or code generated or used during the study are available from the corresponding author by request, including MATLAB code (.m format) of the presented inverse modeling, and MATLAB datasets (.mat format) of the presented numerical results.

This material is based upon work supported by the National Science Foundation under Awards CMMI-2044887 and CMMI-2053694. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation. The authors are grateful for the support by the Faculty Research and Creative Endeavors (FRCE) Research Grant at Central Michigan University. This paper is a contribution of the Central Michigan University Institute for Great Lakes Research. The authors also greatly appreciate the reviewers’ very constructive comments, which substantially helped to improve the paper.