Sylvain Nintcheu Fata

Bojan B Guzina

Marc Bonnet

Computational framework for the BIE solution to inverse scattering problems in elastodynamics *

Keywords: 3D seismic imaging, Cavity detection, Inverse scattering, Elastic waves, Boundary integral equation methods, Adjoint field method

The focus of this paper is a computational platform for the non-intrusive, active seismic imaging of subterranean openings by means of an elastodynamic boundary integral equation (BIE) method. On simulating the ground response to steady-state seismic excitation as that of a uniform, semi-infinite elastic solid, solution to the 3D inverse scattering problem is contrived as a task of minimizing the misfit between experimental observations and BIE predictions of the surface ground motion. The forward elastodynamic solution revolves around the use of the half-space Green's functions, which analytically incorporate the traction-free boundary condition at the ground surface and thus allow the discretization and imaging effort to be focused on the surface of a hidden cavity. For a rigorous approach to the gradient-based minimization employed to resolve the cavity, sensitivities of the trial boundary element model with respect to (geometric) void parameters are evaluated using an adjoint field approach. Details of the computational treatment, including the regularized (i.e.

Cauchy principal value-free) boundary integral equations for the primary and adjoint problem, the necessary evaluation of surface displacement gradients and their implementation into a parallel code, are highlighted. Through a suite of numerical examples involving the identification of an ellipsoidal cavity, a parametric study is presented which illustrates the importance of several key parameters on the imaging procedure including the prior information, "measurement" noise, and the amount of experimental input.

Introduction

Nonintrusive identification of underground openings by elastic waves, which has relevance to diverse areas such as exploration seismology, nondestructive material testing, environmental remediation, and defense applications, has been investigated for a long time. For this class of inverse scattering problems, a variety of computational approaches are available. Most of such imaging solutions are based either on the conventional travel-time analysis (see [START_REF] Aki | Quantitative Seismology[END_REF], or the full waveform interpretation of seismic signals (e.g. Bunks et al., 1995). In the case of three-dimensional seismic exploration, these interpretation techniques necessitate an extensive experimental and computational effort (Sheriff and Geldart, 1995). In many instances, however, there is a need for rapid and inexpensive, yet accurate 3D imaging of major subterranean features on a smaller scale; such is the case with the identification of hidden underground facilities, detection of unexploded ordinances, delineation of abandoned mines, and location of buried wastes. For such problems, imaging algorithms based on boundary integral equation (BIE) methods, which provide a direct link between the observed waveforms and the geometry of a hidden object, could be used to effectively compensate for the limited field data. Although the latter approach has drawn significant attention in acoustic and electromagnetic theories (e.g. [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF][START_REF] Litman | Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set[END_REF] its application to seismic imaging is still in the early stage [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF].

Among numerous reviews dealing with the computational treatment of 3D inverse scattering problems, one may mention Sheriff and Geldart (1995) in exploration seismology, [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF] in acoustics, and Pike and Sabatier (2002) as a general reference spanning applications from elastic to acoustic, electromagnetic, x-ray, and cosmic-scale scattering. Beyond the qualitative application of seismic refraction techniques to underground cavity detection (e.g. [START_REF] Cooper | Geophysical exploration for cavity detection Karst terrain[END_REF], much of the 3D work in exploration seismology revolves around the seismic reflection method [START_REF] Yilmaz | Seismic Data Prospecting[END_REF] which postulates near-vertical wave propagation, i.e. that the source-receiver spacing is much smaller than depth to reflector. To eliminate the effect of multiple reflections discarded (among other phenomena) by the reflection analysis, this method requires significant redundancy of experimental data, as high as 96-fold, in terms of the common midpoint technique [START_REF] Sharma | Environmental and Engineering Geophysics[END_REF]. In the approach, the seismic travel times captured by each in-line array of geophones are interpreted first in a 2D fashion via the far-field approximation (i.e. ray tracing, [START_REF] Aki | Quantitative Seismology[END_REF] that accounts for the in-plane dip and curvature of subsurface reflectors. The resulting images across vertical "slices" are finally migrated (or corrected) for the out-of-plane scattering using 3D ray tracing analysis, often in terms of time-consuming finite difference simulations (Sheriff and Geldart, 1995). An advanced alternative to the foregoing travel-time analysis, which interprets entire seismic records (including multiple reflections) and does not require identification of arrival times is the full waveform analysis.

The computational effort involved with the latter technique, however, is so intensive that most de-velopments so far have been two-dimensional (Bunks et al., 1995;[START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF]. As it stands, consistent 3D seismic imaging of subterranean profiles requires either a major experimental input, high interpretation cost, or both. Moreover, the reflection methods appear to have a limited use in the detection of shallow cavities owing to the interference between surface and reflected body waves [START_REF] Belesky | Seismic and microseismic methods for cavity detection and stability monitoring of near-surface voids[END_REF].

Recently, the use of 3D seismic methods has been advocated for the detection of anti-personnel and anti-tank mines with low metal content, which are almost invisible to ground-penetrating radar and metal detectors. The approach revolves around a concurrent use of elastic and electromagnetic waves via 3D finite difference simulations [START_REF] Scott | Investigation of a technique that uses both elastic and electromagnetic waves to detect buried land mines[END_REF][START_REF] Schröeder | Elastic waves interacting with buried land mines: A study using the FDTD method[END_REF]. In view of its high computational cost, however, the model is developed to run exclusively on large-scale supercomputers such as the Cray T3E.

In contrast to seismic exploration, the developments in radar and sonar technologies have been driven primarily by i) the necessity for rapid object identification, and ii) simplicity of the background (acoustic or electromagnetic) medium where the assumption of homogeneity is often adequate [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]. As a result, the latter 3D imaging solutions are predominantly based on the boundary integral equation (BIE) methods [START_REF] Colton | Recent developments in inverse acoustic scattering theory[END_REF] where the scattering effects of a major object (e.g. a submarine) hidden in a relatively uniform, infinite or semi-infinite domain are intrinsically amenable to the use of fundamental solutions. So far, however, application of the radar and sonar solutions to seismic imaging has been impeded not only by the inherent heterogeneity of geological profiles, but also by the fact that elastic waves, unlike their acoustic counterpart, take many different forms (compressional, shear, Love, Rayleigh and Stoneley waves, see [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]), which renders their BIE interpretation challenging. Few exceptions dealing with the inverse scattering of elastic waves include crack identification in infinite elastic solids by [START_REF] Kress | Inverse elastic scattering from a crack[END_REF] (2D frequency domain study) and [START_REF] Nishimura | Crack determination problems[END_REF] (3D analysis in time domain). To further the relevance of radar and sonar technologies in seismic imaging, the focus of this investigation, which builds on the recent analytical framework in [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF], is the development of a comprehensive computational basis for the identification of cavities via an elastodynamic BIE method, for the case involving three-dimensional elastic wave propagation in a semi-infinite solid. Several key features of this boundary-only imaging technique, which make it suitable for expedient engineering applications, include limited computational effort and full waveform analysis which accounts for both body and surface waves scattered by the cavity.

By means of a well-defined incident seismic field and a set of surface motion sensors used to monitor elastic waves scattered by the cavity, the inverse problem of underground opening identification is reduced to the minimization of a cost function representing the misfit between the field observations and their predictions for an assumed cavity location. For a computationally-effective treatment of the featured body and surface wave fields, the predictive model used in this study [START_REF] Pak | Seismic soil-tructure interaction analysis by direct boundary element methods[END_REF] is based on the fundamental solution for a uniform elastic half-space (see [START_REF] Guzina | On the Analysis of Wave Motions in a Multi-Layered Solid[END_REF], which analytically handles the traction-free boundary condition representing the ground surface. In formulating the gradient search technique employed by the inverse solution, necessary sensitivities of the cost function are evaluated semi-analytically via an adjoint field approach that offers a superior performance relative to conventional finite-difference gradient estimates. Details of the computational treatment, including the regularized boundary integral equations for the primary and adjoint field, featured evaluation of the surface displacement gradients and their implementation into a parallel code, are highlighted. Effectiveness of the proposed imaging method is further elevated via the use of a quasi-Newton minimization algorithm in a parallel computational setting. To investigate the effects of "measurement" noise, prior information, and the amount of experimental input on the imaging procedure, a parametric study is presented for the numerical example involving identification of an ellipsoidal cavity.

Formulation

This investigation deals with an inverse scattering problem for the isotropic, homogeneous elastic half-space containing a hidden cavity, see Fig. 1. With reference to the Cartesian basis {O; ξ 1 , ξ 2 , ξ 3 }, the semi-infinite solid Ω = {(ξ 1 , ξ 2 , ξ 3 )| ξ 3 > 0} is characterized by the Lamé constants λ and µ, mass density ρ, and is bounded on top by the free surface S = {(ξ 1 , ξ 2 , ξ 3 )| ξ 3 = 0}. The void inside the half-space, Ω C ⊂ Ω, is bounded by a smooth closed surface Γ. To aid the ensuing formulation, let Ω - denote the semi-infinite domain surrounding the cavity, i.e. Ω -= Ω\ (Ω C ∪ Γ), and let Γ R = Σ R ∩ Ω be an auxiliary surface, where Σ R is a sphere of radius R centered at the origin. The respective subsets of Ω, Ω -and S that are bounded by Γ R will be denoted as Ω R , Ω - R and S R , with an implicit assumption that R is sufficiently large so that Ω C ⊂ Ω R . With such definitions, n will be identified with the unit normal to S R ∪ Γ ∪ Γ R directed toward the exterior of the "truncated" half-space Ω - R .

Forward problem

To mimic an active seismic imaging configuration, the cavity Γ in Fig. 1 is assumed to be exposed by a vertical vibratory point source, acting at ξ = ζ ∈ S with excitation frequency ω. For brevity, the implicit time-harmonic factor e iωt and the functional dependence of the featured displacement and stress fields on ω will be omitted hereafter. On representing the seismic source by the body force distribution

f (ξ) = P δ(ξ-ζ) e 3 , (1)
where δ stands for the three-dimensional Dirac delta function, e 3 is the unit vector in the ξ 3 -direction, and P is the force magnitude, the total displacement field u(ξ), ξ ∈ Ω -induced in the half-space by ρ , , µ λ

ξ 1 ξ 3 t ω i Pe o Γ n n ζ x m Γ Ω R R Figure 1:
Illumination of an underground cavity by seismic waves f , can be defined as the solution to the elastodynamic boundary value problem

∇•(C : ∇u) + f + ρ ω 2 u = 0, ξ ∈ Ω - t ≡ n • C : ∇u = 0, ξ ∈ Γ ∪ S (2) lim R→∞ Γ R ûk (ξ, x) • t(ξ) - tk (ξ, x) • u(ξ) dΓ ξ = 0, x ∈ Ω - R , k = 1, 2, 3.
In (2), t is the traction vector; C denotes the isotropic elasticity tensor, i.e.

C = λ I 2 ⊗ I 2 + 2 µ I 4 (3)
where I n is the symmetric n-th order identity tensor; vectors ûk (ξ, x) and tk (ξ, x) are respectively the elastodynamic displacement and traction Green's functions for the semi-infinite solid, i.e. the response of the unexcavated half-space Ω at ξ ∈ Ω due to a unit time-harmonic point force acting at x ∈ Ω in the k-th coordinate direction. The far-field condition (2c), implicit to wave propagation problems involving unbounded domains, is known as the generalized radiation condition for the displacement field u (see [START_REF] Pak | Seismic soil-tructure interaction analysis by direct boundary element methods[END_REF].

Inverse problem

On introducing a discrete set S obs ⊂ S of surface control points (ξ = x m , m = 1, 2, . . . M ) where the response of the "excavated" half-space Ω -(due to prescribed excitation f ) is monitored to resolve Γ, the inverse problem of cavity identification can be set forth as a task of finding the void location and shape by minimizing the cost functional

J (Γ) = Q 2 M m=1 u(x m )-u obs (x m ) 2 L 2 (S obs ) + 1 2 q-q p •G• q-q p , (4)
subject to constraint (2) on the trial displacement u (also known as the primary field) (see [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF]. The first term on the right-hand side of (4) can be interpreted as a misfit (over S obs ) between the experimental observations u obs and their elastodynamic predictions u for a trial void geometry Γ. For generality, the cost functional is augmented with a second term to include any prior information on the shape and location of the hidden cavity. In (4), Q is a normalizing factor, G stands for a real symmetric positive definite matrix of weighting coefficients, q p is a constant vector synthesizing the prior information on the geometry of the hidden cavity, and q = q(Γ) is its counterpart evaluated for the trial cavity. Among various geometric parameters that could be used to construct q, one can mention the cavity volume V (Γ), centroid coordinates c i (Γ), and the inertia tensor

I ij (Γ)
defined respectively as

V (Γ) = - 1 3 Γ ξ • n dΓ ξ , c i (Γ) = - 1 4 V Γ ξ i ξ • n dΓ ξ , I ij (Γ) = - 1 5 Γ (ξ i -c i )(ξ j -c j )(ξ -c) • n dΓ ξ , i, j = 1, 2, 3. (5)
One may note that the minus sign in front of surface integrals in (5) accounts for the fact that n is directed toward the interior of the cavity.

Boundary integral framework

With the aid of the elastodynamic reciprocal theorem in frequency domain, the radiating solution u

to the scattering problem (2) can be shown to admit the boundary integral representation

u k (x) = P û3 k (x, ζ) - Γ tk i (ξ, x) u i (ξ)dΓ ξ , x ∈ Ω -. (6)
In ( 6), P û3 k (x, ζ) synthesizes the incident field, i.e. the displacement at x due to seismic source f in the void-free half-space; accordingly, the integral term signifies the scattered field, i.e. the perturbation of the incident field due to presence of the cavity. If the (total) primary field u is assumed to be Hölder continuous, i.e. u ∈ C 0,δ , the limit of (6) as x approaches the boundary Γ can be written in terms of the regularized (i.e. Cauchy principal value-free) boundary integral equation

u k (y) + Γ [ tk i (ξ, y)] 1 {u i (ξ) -u i (y)} dΓ ξ + Γ [ tk i (ξ, y)] 2 u i (ξ) dΓ ξ = P û3 k (y, ζ), k = 1, 2, 3, y ∈ Γ, (7)
where [ tk (ξ, y)] 1 and [ tk (ξ, y)] 2 are respectively the singular and regular parts of the elastodynamic traction Green's function for the semi-infinite solid (see [START_REF] Pak | Seismic soil-tructure interaction analysis by direct boundary element methods[END_REF].

Geometric transformation and material differentiation

To quantify the influence of (geometric) cavity perturbations on the cost functional J (Γ), it is instructive to adopt the framework of continuous deformations and consider small perturbations Ω - R (τ ) ≡ Ω - R (0 ≤ τ < ǫ) of the "truncated" half-space around the reference configuration Ω - R (0), associated with an initial trial cavity Γ 0 . With the assumption that R is sufficiently large so that the trial cavity with boundary Γ(τ ) ≡ Γ remains in Ω R , the foregoing domain evolution with respect to the pseudo-time parameter τ can be formally stated as

ξ = Φ(X, τ ) = X + τ V(X) ∈ Ω - R , τ ∈ [0, ǫ) (8)
where V is the initial transformation velocity, i.e.

V(X) = ∂Φ(X, τ ) ∂τ | τ =0 = V • Φ -1 (ξ, 0) ≡ θ(ξ). (9)
In continuum mechanics, X is usually called the material or Lagrange variable, whereas ξ is referred to as the spatial or Euler coordinate. To facilitate the ensuing developments, the transformation velocity θ is assumed to maintain the external boundaries of Ω - R unchanged, i.e.

θ n ≡ n•θ = 0 on S R θ = 0 on Γ R . (10)
On the basis of ( 8) and ( 9), the material differentiation, ⋆ f , of a given field quantity f (ξ, τ ) at τ = 0 follows directly from the general results of continuum kinematics so that

⋆ f = lim τ →0 1 τ f (ξ, τ ) -f (X, 0) = f ,τ (ξ, 0) + θ(ξ)•∇f (ξ, 0) (11)
where f ,τ = ∂f /∂τ stands for the partial "time" derivative. By virtue of ( 4) and ( 11), the material derivative of the featured cost functional, J (Γ), can be written explicitly as

⋆ J = Q

Sensitivity analysis

A procedure commonly employed in the sensitivity analysis of functionals such as J revolves around the use of suitable finite difference formulae. In the context of ( 12), finite difference computation of ⋆ J translates into numerical estimation of pseudo-time derivative of the primary field, u, τ , due to small cavity perturbations. In view of the well-known drawbacks (in terms of stability, accuracy and computational efficiency) associated with finite difference estimators, however, a closed-form expression for ⋆ J is preferred for an efficient solution to the inverse problem. Adopting the methodology proposed in the context of BIE formulations (e.g. [START_REF] Choi | Boundary integral equation method for optimization of elastic structures[END_REF] and [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF]), a convenient approach to this goal is to incorporate constraints (2) on the primary field into the cost functional via the Lagrangian

L R (u, ũ; Γ) = J (Γ) + Re Ω - R ∇•(C : ∇u) + f + ρω 2 u • ũ dΩ ξ - Γ+S R t• ũ dΓ ξ ( 14
)
where ũ is the Lagrange multiplier. On integrating ( 14) by parts and applying the results for material differentiation of volume and surface integrals [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], it can be shown that the formula for ⋆ L R consists of two distinct components, namely (i) a surface integral containing the initial transformation velocity θ, and (ii) a second, non-homogeneous term. To satisfy the fundamental requirement that the material derivative, ⋆ L R , vanishes for a trivial distribution of the transformation velocity (i.e. no cavity perturbation), the second term must be set to zero. Assuming that ũ ∈

C 2 (Ω -) ∩ C 1 (Ω -∪ Γ ∪ S)
, this requirement can be realized provided that the vector field ũ satisfies the following equations [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF])

∇•(C : ∇ũ) + Q M m= (u -u obs )δ(ξ-x m ) + ρ ω 2 ũ = 0, ξ ∈ Ω - t ≡ n • C : ∇ũ = 0, ξ ∈ Γ ∪ S (15) lim R→∞ Γ R ûk (ξ, x) • t(ξ) - tk (ξ, x) • ũ(ξ) dΓ ξ = 0, x ∈ Ω - R , k = 1, 2, 3.
One may note that (15) define ũ as an elastodynamic state, herein termed the adjoint field, over Ω -; in fact, the only difference between (15) and the primary field ( 2) is that the adjoint body force distribution, f , is defined through

f (ξ) = Q M m= (u -u obs )δ(ξ-x m ), (16)
as a system of point forces, located at observation points x m ∈ S obs , which are proportional to the misfit between experimental observations u obs and theoretical predictions u. On the basis of ( 7) and ( 16), the regularized boundary integral equation for the adjoint field ũ can be directly written as

ũk (y) + Γ [ tk i (ξ, y)] 1 {ũ i (ξ) -ũi (y)} dΓ ξ + Γ [ tk i (ξ, y)] 2 ũi (ξ) dΓ ξ = Q M m=1 ûi k (y, x m )(u i (x m ) -u obs i (x m )), k = 1, 2, 3, y ∈ Γ, ( 17
)
with the right-hand side synthesizing the effect of the fictitious body force f. 2), ( 14) and ( 15), it can be shown that the material derivative of (4),

ρ , , µ λ Ω ξ 1 ξ 3 o Γ n
⋆ J ≡ ⋆ L R , permits the integral representation ⋆ J = Re Γ ρω 2 u• ũ -µ ∇ S u : 2λ λ+2µ I 2 ⊗ I 2 + 2I 4 -n ⊗ I 2 ⊗ n : ∇ T S ũ θ n dΓ ξ + (q-q p )•G• ⋆ q, (18)
where ∇ S u is the surface gradient of u, defined in terms of its normal derivative u, n = n•∇u as

∇ S u = ∇u -n ⊗ u, n (19)
Formula ( 18) is amenable to an effective computational treatment via boundary element techniques as it is written in terms of tangential derivatives that are directly computable from the nodal values of u and ũ on Γ.

Computational platform

For imaging purposes, the location and shape of Γ can be assumed to depend on a finite set of design parameters, p = (p 1 , p 2 , . . . , p D ). With such hypothesis, partial derivatives ∂J /∂p d required for the minimization of J can be obtained directly from (18) by letting τ = p d , d = 1, 2, . . . , D. A usual procedure in numerical evaluation of integrals such as ( 18) is to assume that Γ = K k=1 E k , where {E k } K k=1 are closed and non-overlapping surface elements. With such partition, each boundary element E k can be parametrized by a mapping E → E k that introduces local coordinates, η = (η 1 , η 2 ) ∈ E, over E k ⊂ Γ where E is a polygonal domain in R 2 . As long as Γ(τ ) remains simply connected, the evolving boundary element mesh Γ h = K k=1 E h k approximating Γ(τ ) can be generated by interpolating a suitable set of parameter-dependent nodes ξ q (p) ∈ E k with fixed, i.e. pre-defined mesh connectivity.

To this end, the Q-noded approximation (E h k ) of a generic surface element E k ⊂ Γ can be written as

ξ(η) = Q q=1 N q (η) ξ q (p), ξ ∈ E h k , ξ q ∈ E k η ∈ E, (20)
where N q (η) are the shape functions for the Q-noded element E h k satisfying the fundamental property

Q q=1 N q (η) = 1, η ∈ E (21)
over the parent domain E (e.g. [START_REF] Brebbia | Boundary element techniques[END_REF]. In a manner similar to (20), the primary displacement field u at an arbitrary point ξ ∈ E h k can be expressed as

u = Q q=1 N q (η) u q , (22)
where u q is the nodal displacement at ξ q ∈ E k . With such approximation, the surface gradient of the primary field ∇ S u can be readily computed as

∇ S u = Q q=1 ∇ S N q (η) ⊗ u q . ( 23
)
Eq. ( 23) demonstrates how tangential derivatives of u can be shifted onto surface derivatives of the shape functions N q (η).

For the gradient-based minimization of J , each partial derivative ∂J /∂p d is computed via (18)

with the (prescribed) transformation velocity θ = θ d , associated with the variation of the d-th design parameter, given by

θ d = Q q=1 N q (η) ∂ξ q ∂p d , η ∈ E. ( 24
)
On the basis of (20), ( 22) and the set {y l , l = 1, 2, . . . , L} ⊂ Γ of all nodal points on Γ, the discretized form of the boundary integral equation ( 7) for the primary field can be expressed as

K k=1 Q q=1 E h k [ tj i (ξ(η), y l )] 1 (N q (η) -δ q k ,l ) dΓ ξ + E h k [ tj i (ξ(η), y l )] 2 N q (η)dΓ ξ u q k i +    δ ji - K k=1   1- Q q=1 δ q k ,l   E h k [ tj i (ξ(η), y l )] 1 dΓ ξ    u l i = P û3 j (y l , ζ), j = 1, 2, 3 l = 1, 2, . . . , L (25)
in terms of nodal displacements u l = u(y l ), l = 1, 2, . . . , L, where δ ji is a Kronecker delta and q k is the global number of the q-th element node. The integrals in (25) are approximated by a product

Gaussian quadrature in which a special mapping technique [START_REF] Guzina | Seismic Response of Foundations and Structures in Multilayered Media[END_REF], similar to the Duffy's transformation (e.g. Golberg and [START_REF] Chen | Discrete Projection Methods for integral Equations[END_REF], is used to reduce the weak singularity of the integrand when the integration is performed over an element that contains the collocation point. The system of equations ( 25) can be recast in the matrix form as

HU = F, F = P û3 j (y l , ζ) j = 1, 2, 3 L l=1 (26)
where

U = u 1 1 , u 1 2 , u 1 3 , . . . , u L 1 , u L 2 , u L
3 is a vector containing all nodal displacements, H is a (complex non-symmetric) 3L × 3L displacement influence matrix, and F is a vector whose entries are the values of the right-hand side of (25) evaluated at collocation points y l .

Owing to the fact that the elastodynamic field equations governing the primary and adjoint field (under identical boundary conditions) differ only in the body force term (cf. ( 1) and ( 16)), a linear algebraic system of equations for the values of adjoint field (17) at collocation points can be written as

H Ũ = F, F = Q M m=1 ûi k (y l , x m )(u i (x m ) -u obs i (x m )), k = 1, 2, 3 L l=1 (27)
where the influence matrix H is the same as in (26).

Shape parametrization

In this study, motion of the evolving trial boundary, Γ(τ ), within the half-space Ω = {(ξ 1 , ξ 2 , ξ 3 )| ξ 3 > 0} is restricted to translation and stretch in the Cartesian frame {O; ξ 1 , ξ 2 , ξ 3 }. For problems involving identification of a single cavity, such description entails the use of a six-dimensional parametric space p = (c 1 , c 2 , c 3 , α 1 , α 2 , α 3 ) which incorporates the cavity's centroidal motion (c i , i = 1, 2, 3) and principal stretches (α i > 0, i = 1, 2, 3) in the above Cartesian basis. With such definition, analytical dependence of the nodal coordinates, ξ q = ξ q (p), on the trial boundary Γ is introduced as an affine deformation of a fictitious shape with respect to Lagrange coordinates (X 1 , X 2 , X 3 ) so that

ξ q i = α i X q i + c i , i = 1, 2, 3 (28)
assuming no summation over index i. By virtue of ( 20), ( 21), ( 24)and ( 28), the normal transformation velocity associated with the variation of the d-th design parameter,

θ d n ≡ n • θ d (d = 1, • • • , 6
), can be shown to admit the representation

θ d n = n 1 , n 2 , n 3 , ξ 1 -c 1 α 1 n 1 , ξ 2 -c 2 α 2 n 2 , ξ 3 -c 3 α 3 n 3 . ( 29
)

Surface derivatives

On the basis of ( 18) and ( 29), an explicit expression for the partial derivative ∂J /∂p d directly follows in the form of

∂J ∂p d = Re K k=1 E h k ρω 2 u• ũ - 2λµ λ+2µ div S u div S ũ -µ ∇ S u + ∇ T S u : ∇ T S ũ + µ ∇ S u•n)•(∇ T S ũ•n θ d n dΓ ξ + (q-q p )•G• ⋆ q, d = 1, • • • , 6 (30)
where div S u is the surface divergence of the displacement field u, computable as

div S u = Q q=1 u q •∇ S N q (η), (31)
over a representative boundary element E h k . By use of ( 22), a discretized formula for u• ũ is readily obtained as

u• ũ = Q p=1 Q q=1 N p (η)N q (η) u p • ũq . (32)
Similarly, with the aid of ( 23), the remaining part of the bracketed term in (30) can be written as

∇ S u + ∇ T S u : ∇ T S ũ -∇ S u•n)•(∇ T S ũ•n = Q p=1 Q q=1 u p •∇ S N q (η) ũq •∇ S N p (η) + u p • ũq -(u p •n)(ũ q •n) ∇ S N p (η)•∇ S N q (η) . ( 33
)
To evaluate the tangential derivative ∇ S N q (η) featured in ( 31) and ( 33), it is useful to recall several basic formulae of differential geometry. On assuming that the parametrization ξ = ξ(η 1 , η 2 ) given by ( 20) is regular, a local companion basis {r 1 , r 2 , n} can be established at every point ξ ∈ E h k where As shown in Fig. 3, vectors r 1 and r 2 are tangent to E h k , with the triplet {r 1 , r 2 , n} forming a positiveoriented basis of R 3 . By virtue of (34), the covariant components g IJ of the metric tensor g over the surface element E h k can be expressed as

r 1 = Q q=1 ∂N q ∂η 1 ξ q , r 2 = Q q=1 ∂N q ∂η 2 ξ q , n = r 1 × r 2 ||r 1 × r 2 || . ( 34
g IJ = r I • r J , I, J = 1, 2 (35)
Likewise, the contravariant components g IJ of g are given by

g IL g LJ = δ I J . ( 36
)
where δ I J denotes the Kronecker delta. With the above definitions, the surface gradient of shape functions ∇ S N q takes an explicit form 37) in terms of the local basis r J , J = 1, 2. On the basis of ( 36), the differential area dΓ ξ on E h k can also be expressed in terms of local coordinates as

∇ S N q = ∂N q ∂η I g IJ r J , q = 1, 2, . . . Q (
dΓ ξ = g 11 g 22 -(g 12 ) 2 dη 1 dη 2 . ( 38
)
By use of ( 20), ( 29), ( 31) to (33), and ( 35) to ( 38), the sensitivities ∂J /∂p d of the cost functional (4) with respect to cavity perturbations directly follow from (30) once the primary field u q and adjoint field ũq are computed at nodal points y l ∈ Γ.

Parallel computation

The regularized BIE solution (see respectively ( 7) and ( 17) for the primary and adjoint field) is implemented, together with the material derivative (30), in a data parallel code using the message-passing interface [START_REF] Pacheco | Parallel Programming with MPI[END_REF]. The message-passing interface (MPI) is a freely available specification for the message-passing libraries. Data-type parallelism, of interest in this investigation, typically applies when identical operations can be performed concurrently on different data items. Such is the case with repeated, yet computationally expensive evaluations of the elastodynamic Green's functions, ûk (y l , ζ) and tk (y l , ζ), at collocation points y l (l = 1, 2, . . . L) in ( 26) and ( 27). Thus, the running time of the overall program can be significantly reduced by block-distributing the computation of the coefficient matrix H and vectors F and F featured in ( 26) and ( 27) among all processes participating in the parallel job. On denoting the number of participating processes by n p , this distribution of the computational effort for L collocation points can be summarized as

L = n p × q + R ≡ R(q+1) + (n p -R)q, 0 ≤ R < n p -1 (39)
where the first R processes are each assigned q+1 blocks of data items, so that every of the remaining n p -R processes receives q blocks of computation. Each process works on a different data set; upon completion of their respective computing tasks, they exchange the output data using message-passing to solve the problem. In this study, eight processes are used to solve the inverse scattering problem.

To measure the performance of the parallel code, the speed-up ratio, defined as the ratio of the elapsed time of a serial program over that of its parallel counterpart with n p processes, is adopted. On a Linux cluster with 2 nodes, each with 4 Intel Xeon processors, the speed-up factor of 4.58 is obtained when evaluating the objective function J and its partial derivatives ∂J /∂p d (d = 1, • • • , 6) at a sample location of the trial cavity Γ. This result indicates a considerable gain in running time of the parallel code.

Minimization

In this study, minimization of J (p) is effected using the BFGS quasi-Newton method (e.g. [START_REF] Nocedal | Numerical Optimization[END_REF]) and an inexact line search algorithm

φ(β) ≡ J (p k + β d k ) → min, ( 40
)
where d k is the descent direction at the k-th major iteration, p k . To determine the optimal step length, β β * > 0, the line search uses mixed quadratic and cubic polynomial interpolation that satisfies the

sufficient decrease condition φ(β * ) ≤ φ(0) + γ 1 β * φ ′ (0), (41)
as well as the curvature condition

|φ ′ (β * )| ≤ γ 2 |φ ′ (0)|, ( 42
)
with 0 < γ 1 < γ 2 < 1 [START_REF] More | Line search algorithms with guaranteed sufficient decrease[END_REF][START_REF] Nocedal | Numerical Optimization[END_REF]. In the minimization presented below, the values of γ 1 = 10 -4 and γ 2 = 0.99 are adopted for the BFGS quasi-Newton algorithm. To take into account the physical inequality constraints

c 3 -α 3 > 0, α 1 > 0, α 2 > 0, α 3 > 0 (43)
on the cavity centroid and principal stretches featured in the vector of design parameters, p = (c 1 , c 2 , c 3 , α 1 , α 2 , α 3 ), every trial step length β (including the optimum one) is chosen so that it conforms with (43).

Results

In the following examples, surfaces of the "true" and trial cavities are discretized via eight-node isoparametric quadratic boundary elements [START_REF] Brebbia | Boundary element techniques[END_REF]. The elastodynamic half-space Green's functions, featured in (25) and its counterpart for the adjoint field, are computed semianalytically using the Hankel transform approach aided by numerical contour integration and the method of asymptotic decomposition [START_REF] Guzina | On the Analysis of Wave Motions in a Multi-Layered Solid[END_REF].

Sensitivity estimates

To illustrate the performance of the adjoint field approach as a tool for sensitivity calculation, a numerical experiment was performed with reference to the spherical cavity of known radius so that its relevant vector of design parameters, p = (p 1 , p 2 , p 3 ), represents the position of its center. In this case, motion of the trial spherical cavity reduces to the rigid-body translation in Ω, and the normal transformation velocity θ d n , given more generally by ( 29), scales down to The true and trial cavities are centered respectively at p true = (0, 0, 2a) and p 0 = (-2a, 2a, 6a), where a denotes the radius of the sphere (see Fig. 4). The testing configuration shown in the Figure is comprised of nine control points arranged in a square grid over the free surface S. In sequence, each grid node is taken as a location of the vertical point source (magnitude P = 0.2µa 2 ), with the remaining eight nodes used as receiver locations, so that a total of 9 × 8 × 3 = 216 synthetic observation data are generated for the test problem. The elastic parameters of the half-space and the testing frequency are chosen so that

θ d n = (n 1 , n 2 , n 3 ) . (44
λ = 3 2 µ, ω = ωa µ/ρ = 1. ( 45
)
Table 1 compares the performance of the adjoint field formula (30) with the central difference estimates (step size ∆p d = 0.002a) in terms of sensitivities ∂J /∂p d (d = 1, 2, 3) for an assumed location of the trial cavity. For the clarity of comparison, the cost function J and its derivatives are computed with Q = 10 6 µ a and no prior information. As can be seen from the display, the two methods yield consistent results. It should be noted, however, that the computational time for the adjoint field approach is approximately 1/6 of that for the central difference method. For a general sensitivity problem involving D design parameters, the foregoing efficiency ratio can be extrapolated as 1/(2D) since the central difference method requires 2D computations of the BIE (7), each corresponding to a different (i.e. perturbed) configuration of Γ.

Void mapping

As an example of the minimization procedure, the next problem deals with the elastic-wave imaging of an ellipsoidal void whose semi-axes are aligned with the global coordinate system. The surface of the cavity is discretized via 96 eight-node boundary elements. As specified earlier, the void geometry is parametrized in terms of its centroid coordinates c i and semi-axes lengths α i (i = 1, 2, 3), resulting in a six-dimensional parametric space p = (c 1 , c 2 , c 3 , α 1 , α 2 , α 3 ).

The "true" cavity geometry is given by p true = (-4a,-2a, 4a, 1.8a, 0.9a, 0.6a); its trial counterpart is taken as p 0 = (-1.5a, -0.5a, 5a, a, a, a). In succession, the cavity is exposed using nine vertical point sources with magnitude P = 0.2µa 2 according to the testing grid depicted in Fig. 4. For each source location, surface displacements (u obs ) are monitored over 36 control points uniformly spaced over the square observation area (14a × 14a) as shown in Fig. 5; to maintain the physical relevance, observations coinciding spatially with the seismic source are discarded. To mimic the effect of modeling and measurement uncertainties, synthetic observations obs are "contaminated" via the perturbation factor (1+̺) applied to their scattered component, where ̺ is a random variable uniformly distributed over the interval [-η, η]. In what follows, the noise amplitude is varied from η = 0.05 to η = 0.3.

For consistency with the previous example, the constitutive parameters for the half-space and testing frequency are chosen after (45), resulting in the shear wave length (Λ = 2πa) that is approximately twice the largest diameter of the "true" cavity. A priori information To expose the effect of prior information on the featured imaging procedure, the cost function ( 4) and its sensitivities are computed with η = 0.05 assuming both i) no prior information on the hidden cavity i.e. by setting the quadratic form 1 2 q-q p •G• q-q p to zero, and ii) that the volume of the hidden cavity, V true = 4.072a 3 , is known beforehand. On setting the normalizing factor in (4) to Q = 10 6 µa, the latter prior information is incorporated into the minimization problem by taking q = V (Γ), q p = V true and G = 10 -6 µa -3 . It should be noted that the penalty coefficient G is chosen relative to Q so that the prior knowledge component in ( 4) is an order of magnitude larger than its misfit counterpart at p = p 0 . Figs. 6a and 6b graph respectively the iterative process of mapping the hidden cavity, p true , starting from p = p 0 for the cases involving i) no prior knowledge, and ii) prior information on the cavity volume. As can be seen from the diagram, the optimization procedure converges to the global minimum for both cases, although with somewhat different convergence patterns. It is interesting to note from Fig. 6 that the case with prior information requires more iterations for convergence, primarily because of the restricted "movement" in the parametric space caused by the penalty term 1 2 q-q p •G• q-q p . For clarity, Fig. 7 depicts the above minimization procedure geometrically in a 3D setting. Experimental input As examined before, the imaging process illustrated in Figs. 6 and7 is based on the synthetic experimental input with η = 0.05 that incorporates nine point sources (i.e. 9 incident seismic fields) and 36 control points x m (see Fig. 5). To assess the effect of the density of spatial sampling on 3D seismic imaging, the minimization procedure was repeated by retaining all nine point sources, but varying the number of observation points from 9 (3 × 3 grid) to 16 (4 × 4 grid), 25, 36, 49, and 64 (8 × 8 grid). In each case, the observation points were uniformly spaced over the square observation area 14a × 14a. From numerical simulations for the case that assumes prior information on the cavity volume, it was found that the imaging procedure successfully converges to the global minimum (i.e. the true cavity) for all receiver configurations tested. In contrast, the case with no prior information requires at least 25 observation points (i.e. a 5 × 5 grid) for global convergence. For completeness, it should be mentioned that the ratio Λ/s, where Λ is the shear wave length and s is the smallest inter-receiver spacing, ranges from 3.14 for the 8 × 8 grid, down to 0.90 for the 3 × 3 grid of observation points.

Measurement noise

To examine the susceptibility of the minimization algorithm to measurement errors, the imaging problem involving prior information on the cavity volume (see Figs. 6b and7b) is re-examined under varying levels of random noise (added to the scattered component of synthetic observations) ranging from η = 0.05 to η = 0.30. From the numerical simulations summarized in Table 2, it was found that the minimization procedure fails to reach the global minimum (i.e. to identify the correct cavity location and geometry) for noise levels exceeding η = 0.25. Although this important topic clearly warrants a more systematic approach, e.g. in terms of the Bayesian framework [START_REF] Tarantola | Inverse Problem Theory[END_REF], present results nonetheless provide a point of reference for the signal-to-noise ratios that may be acceptable for the three-dimensional, boundary-only seismic imaging. Estimates of the true cavity, p true = (-4a,-2a, 4a, 1.8a, 0.9a, 0.6a)

(noise) p 1 = c 1 p 2 = c 2 p 3 = c 3 p 4 = α 1 p 5 = α 2 p 6 = α 3
0.05 -3.9999a -1.9998a 4.0008a 1.7923a 0.9018a 0.6013a 0.10 -3.9999a -1.9996a 4.0017a 1.7835a 0.9035a 0.6031a 0.20 -3.9998a -1.9991a 4.0033a 1.7668a 0.9071a 0.6062a 0.25 -3.9999a -1.9990a 4.0041a 1.7584a 0.9091a 0.6076a 0.30 2.2240a -1.0378a 7.9997a 0.6025a 0.6463a 2.4913a

Initial guess It should be emphasized that the success of the foregoing imaging technique is strongly dependent on the choice of a starting point, a pitfall that is common to all gradient-based algorithms. This is illustrated in Fig. 8 where the cost function J (c 1 , c 2 , 4a, 1.8a, 0.9a, 0.6a) is plotted for the reduced case involving only two cavity parameters, i.e. its centroid coordinates c 1 and c 2 , no prior information, and no measurement noise (η = 0). The results are shown for two dimensionless frequencies, namely ω = 1 and ω = 3, and two densities of control points, corresponding respectively to the 4 × 4 and 8 × 8 testing grids. Regardless of the excitation frequency and testing configuration, the cost functions in Fig. 8 are highly non-convex, a problem caused primarily by the oscillatory nature of wave patterns in the "background" half-space medium rather than the complexity of cavity geometry (see, e.g. Pike and Sabatier, 2002, for discussion in the context of conventional seismic imaging). This is substantiated by the observation from Fig. 8 that the number of local extrema increases, and the so-called basin of attraction (i.e. the region around the global minimum where the cost function maintains its convexity) narrows down with increasing frequency, i.e. diminishing wave length in the half-space.

Conclusions

In this paper, a computational treatment is presented for the inverse elastodynamic problem involving boundary-only imaging of cavities hidden in a homogeneous semi-infinite solid. On employing the regularized boundary integral equation (BIE) method to simulate the ground response to surface seismic excitation, the inverse scattering problem is reduced to the gradient-based minimization of a misfit between the observed surface motion and its BIE prediction for an assumed void location.

In the formulation, necessary sensitivities of the predictive model with respect to (geometric) cavity parameters are evaluated semi-analytically via an adjoint problem approach which circumvents the major computational drawbacks of finite difference estimates. Computational details of the sensitivity analysis, including the regularized boundary integral equations for the primary and adjoint problem as well as the necessary evaluation of surface displacement gradients, are highlighted. Due to repeated, yet independent evaluation of the elastodynamic half-space Green's functions required by the BIE solution, the imaging procedure, which revolves around a modified-Newton optimization algorithm, is implemented in a parallel computational setting via the message-passing interface. Numerical results show that the adjoint problem approach produces sensitivity estimates that are consistent with their finite difference counterparts, while reducing the computational effort by a factor of 2D where D is the number of design parameters used to describe the cavity. With the aid of noise-contaminated synthetic observations, performance of the imaging technique is illustrated for the case of a single ellipsoidal cavity hidden in a semi-infinite solid. On varying the input to this numerical example, a limited parametric study is presented which highlights the importance of factors such as prior information, measurement noise, and the amount of experimental input in boundary-only seismic imaging.

Figure 2 :

 2 Figure 2: Seismic sources for the adjoint field

Figure 3 :

 3 Figure 3: Local basis on the surface element

Figure 4 :

 4 Figure 4: Spherical cavity and testing grid in the elastic half-space

Figure 5 :

 5 Figure 5: Testing configuration with nine sources and 36 control points

Figure 6 :

 6 Figure 6: Evolution of design parameters in the minimization process (ω = 1, η = 0.05)

Figure 7 :

 7 Figure 7: Sample of major iterates in the minimization process (ω = 1, η = 0.05)

Figure 8 :

 8 Figure 8: Cost function in the two-dimensional parametric space: a) 4 × 4 testing grid and b) 8 × 8 testing grid (η = 0)

Table 1 :

 1 Sensitivity estimates: central difference versus adjoint field approach Mesh 96 elements, J = 61.27 µa 3 294 elements, J = 61.64 µa 3

	Method	Finite diff.	Adjoint	Finite diff.	Adjoint
	1/(µ a 2 ) ∂J/∂p 1	0.9328	0.9200	0.9366	0.9304
	1/(µ a 2 ) ∂J/∂p 2	-0.9328	-0.9200	-0.9366	-0.9304
	1/(µ a 2 ) ∂J/∂p 3	-4.1679	-4.1483	-4.1871	-4.1849

Table 2 :

 2 Sensitivity of minimization to the level of measurement noise (ω = 1)

	η

M m=1 Re (u(x m , τ ) -u obs (x m ))•u, τ (x m , 0) + q-q p •C• ⋆ q, (12)where ⋆ q is comprised of material derivatives of the prior information entries listed in (5). By virtue of a generic formula for material differentiation of integrals over a regular moving surface (e.g.[START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], these entries can be differentiated as⋆ V = -Γ θ n dΓ ξ , ⋆ c i = -1 V Γ (ξ i -c i ) θ n dΓ ξ , ⋆ I ij = -Γ (ξ i -c i )(ξ j -c j ) θ n dΓ ξ (13)in terms of the (prescribed) normal transformation velocity θ n (ξ).

Acknowledgment The support provided by the National Science Foundation through CAREER Award No. CMS-9875495 to B. Guzina and the University of Minnesota Supercomputing Institute during the course of this investigation is gratefully acknowledged. Special thanks are due to MTS Systems Corporation for providing the opportunity for M. Bonnet to visit the University of Minnesota through the MTS Visiting Professorship of Geomechanics.


B B Guzina

M Bonnet

TOPOLOGICAL DERIVATIVE FOR THE INVERSE SCATTERING OF ELASTIC WAVES

Topological derivative for the inverse scattering of elastic waves

Introduction

Stress-wave identification of cavities and objects embedded in an elastic solid is a longstanding problem in mechanics prompted by its applications in exploration seismology, nondestructive material testing, and underground facility detection. For this class of inverse scattering problems, often associated with the semi-infinite domain assumption for the 'host' elastic solid and free-surface distribution of motion sensors, a variety of approaches are available. Most of such imaging solutions are based either on the far-field approximation of the wave equation (ray theory, see [START_REF] Sheriff | Exploration Seismology[END_REF]), or its finite-difference analogue [START_REF] Plessix | [END_REF]. In order Q. Jl Mech. Appl. Math. (2004) 57,161:179 (??), 1-20 c Oxford University Press 2004

to provide comprehensive three-dimensional subterranean images, however, these domain methods commonly necessitate an extensive experimental and computational effort.

For problems where rapid imaging of underground openings is required, on the other hand, the boundary integral equation (BIE) formulations, which furnish a direct mathematical link between the observed surface waveforms and the geometry of a hidden object, could be used to effectively compensate for the limited field data and expedite the interpretation process (see [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF] for acoustic problems). In the context of elastic-wave imaging, however, an appreciable computational cost of the relevant half-or full-space fundamental solutions precludes the use of BIE methods in conjunction with global search techniques like genetic algorithms which entail a large number of forward simulations. As a result, an expedient BIE identification of subterranean cavities (or material defects) via elastic waves is amenable primarily in terms of the gradient-based optimization techniques, especially with the help of analytical shape sensitivity estimates [START_REF] Guzina | [END_REF]. Unfortunately, the stand-alone use of gradienttype minimization for such purpose is not satisfactory for its success is strongly dependent on a reliable prior information about the geometry of the hidden void.

These considerations led the authors to investigate the usefulness of the concept of topological derivative in connection with the elastodynamic inverse problem. Considering a generic cost functional J of the void shape, size, location, and topology, its topological derivative, T (x 0 ), synthesizes the sensitivity of J with respect to the creation of an infinitesimal cavity at a prescribed location x 0 inside the host solid. The information provided by the topological derivative distribution T (x 0 ) is thus potentially very useful, providing a rational basis for selecting the void topology and its initial location/geometry, both of which are necessary for the application of gradient-based minimization to the inverse scattering problem at hand. The concept of topological derivative first appeared in (5) and [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Kopchpositionierungkrieterien[END_REF] in the context of shape optimization of mechanical structures. Recently, its rigorous mathematical formulation has been established within the framework of elastostatic problems and Laplace equation [START_REF] Sokolowski | [END_REF]8). Beyond its direct application to the structural shape optimization problems, however, the topological derivative may also have potential use in solving inverse problems, a topic in its early stages of investigation to which the present article is intended to contribute (see also [START_REF] Gallego | Proc. IABEM 2002 Conference[END_REF] for two-dimensional elastostatics).

In this communication, the concept of topological derivative is extended to threedimensional inverse scattering of elastic waves involving semi-infinite and infinite domains. In the approach, the formula for topological derivative (explicit in terms of the fundamental solution) is obtained within the framework of BIE methods by an asymptotic expansion of the featured cost functional with respect to the creation of an infinitesimal hole in the intact host solid. Dependent on the prescribed shape of the infinitesimal cavity, the formula involves the solution to six exterior elastostatic problems, and becomes fully explicit if the infinitesimal cavity is spherical. A set of numerical results is included to illustrate the approach, which could be used most effectively as a pre-conditioning tool for more accurate gradient-based imaging algorithms [START_REF] Guzina | [END_REF]. Reliant on the availability of suitable fundamental solutions, the proposed method could be directly extended to elastic-wave imaging of finite solid bodies.

Inverse problem

Consider the imaging problem depicted in Fig. 1 where the semi-infinite solid, probed by elastic waves, houses a hidden void (or a system thereof). With the Cartesian frame {O : ξ 1 , ξ 2 , ξ 3 } set at the top surface S, the homogeneous isotropic half-space Ω = {(ξ 1 , ξ 2 , ξ 3 )|ξ 3 > 0} is characterized by the Lamé constants λ and µ, and mass density ρ. Without loss of generality, it will be assumed that the sought cavity, Ω true C ⊂ Ω, is illuminated by a timeharmonic point source f acting at ξ = x S ∈ Ω with frequency ω. For imaging purposes, the induced surface motion u obs is monitored over a finite set of control points ξ = x m ∈ Ω (m = 1, 2, . . . M ).

Following the conventional approach (4), the cost function for the inverse scattering problem described in Fig. 1 can be defined in a misfit-type fashion as

J (Ω\Ω C ; f ) = 1 2 M m=1 {u(x m , ω)-u obs (x m , ω)}•W m •{u(x m , ω)-u obs (x m , ω)}, (2.1)
where Ω C = Ω C ∪ Γ indicates the closure of a trial cavity Ω C bounded by a closed smooth surface Γ; over-bar symbol denotes complex conjugation; W m (m = 1, 2, . . . , M ) are suitable Hermitian and positive definite matrices, and u is the displacement field which solves the elastodynamic scattering problem for the semi-infinite solid Ω\Ω C : u satisfies the field equations

∇•(C : ∇u) + f = -ρ ω 2 u, ξ ∈ Ω \ Ω C , (2.2)
and Neumann boundary conditions

t ≡ n• C : ∇u = 0, ξ ∈ Γ ∪ S. (2.3)
Here, t denotes the surface traction, n is the unit normal on Γ outward to Ω\Ω C , and

C = λ I 2 ⊗ I 2 + 2µ I sym 4 (2.4)
is the isotropic elasticity tensor with I 2 and I sym 4 symbolizing the second-order and symmetric fourth-order identity tensors, respectively.

In what follows, it is assumed that u satisfies the standard continuity requirements for smooth bounding surfaces,

u ∈ C 2 (Ω \ Ω C ) ∩ C 1 ((Ω \ Ω C ) ∪ S).
For this class of scattering problems, the displacement field u satisfying (2.2) and (2.3) can be shown to admit a Somigliana-type integral representation [START_REF] Pak | [END_REF]

u(x, ω) = e k Γ t(ξ, ω; n)• ûk (ξ, x, ω) dΓ ξ -e k Γ u(ξ, ω)• tk (ξ, x, ω; n) dΓ ξ + e k Ω\ΩC f (ξ, ω)• ûk (ξ, x, ω) dΩ ξ , x ∈ Ω\Ω C , (2.5)
where e k is the unit vector in the ξ k -direction (k = 1, 2, 3); Einstein summation convention is assumed over the spatial coordinate index k, and ûk (ξ, x, ω) and tk (ξ, x, ω; n) constitute the elastodynamic fundamental solution for a uniform semi-infinite solid by denoting the respective displacement and traction vectors at ξ ∈ Ω due to a unit (time-harmonic) point force acting at x ∈ Ω in the kth direction. As shown in (11), these Green's functions can be decomposed into a singular part ([û] 1 , [ t] 1 ) and a residual (regular

) component ([û] 2 , [ t] 2 ) via ûk i (ξ, x, ω) = [û k i (ξ, x, ω)] 1 + [û k i (ξ, x, ω)] 2 , tk i (ξ, x, ω; n) = [ tk i (ξ, x, ω; n)] 1 + [ tk i (ξ, x, ω; n)] 2 , (2.6)
where, for point forces located at a non-zero distance from the free surface (x 3 > 0), [û] 1 and [ t] 1 are given by the (elastostatic) Kelvin solution for an infinite solid.

For the ensuing treatment, the total displacement field u featured in (2.5) can be conveniently decomposed as

u(ξ, ω) = u F (ξ, ω) + u S (ξ, ω), ξ ∈ Ω \ Ω C , (2.7)
where u S denotes the scattered field, and u F is the free (incident) field defined as the response of a cavity-free half-space Ω due to given body force distribution f . On the basis of the traction-free boundary condition (2.3) and the Maxwell-type symmetry of the fundamental solution where ûk i (ξ, x, ω) = ûi k (x, ξ, ω), the first and the third term on the right-hand side of (2.5) reduce respectively to zero and u F , resulting in an integral representation for the scattered field in the form of

u S (x, ω) = -e k Γ u(ξ, ω)• tk (ξ, x, ω; n) dΓ ξ , x ∈ Ω \ Ω C .
(2.8)

Topological derivative

To search the semi-infinite domain Ω for cavities in the context of (2.1), let B a (x O ) = x O +aB define the cavity of size a > 0 and volume a 3 |B|, where B ⊂ R 3 is a fixed and bounded open set of volume |B| containing the origin. For further reference, B a is assumed to be bounded by a simply connected, smooth surface Γ a , with B a = B a ∪ Γ a denoting its closure. Without loss of generality, B is chosen so that B a (x O ) is contained inside the sphere of radius a centred at x O (see Fig. 1). With such definitions, the topological derivative of (2.1) can be defined as

T (x O , f ) = lim a→0 J (Ω\B a ; f ) -J (Ω; f ) |B a | = lim a→0 J (Ω\B a ; f ) -J (Ω; f ) a 3 |B| , x O ∈ B a , (3.1)
which furnishes the information about the variation of J (Ω; f ) if a hole of prescribed shape B and infinitesimal characteristic size is created at x O ∈ Ω. Within the framework of shape optimization, it was shown [START_REF] Sokolowski | [END_REF]8) that the elastostatic equivalent of (3.1) can be used as a powerful tool for the grid-based exploration of a solid for plausible void regions in terms of the chosen shape functional J . In what follows, this concept will be extended to elastic-wave imaging of semi-infinite solids on the basis of the elastodynamic fundamental solution for a homogeneous isotropic half-space. As pointed out earlier, the topological derivative is investigated here as a pre-conditioning tool for selecting the initial 'guess' for an optimization-based approach to the solution of inverse scattering problems, and is therefore considered only in relation to the unperturbed half-space Ω, as suggested by (3.1).

The reader is referred to the end of Section 5 for further discussion on this subject.

Small cavity approximation

In the context of (3.1), the trial cavity featured in (2.5) to (2.8) can be specialized to a small void of size a containing x O so that Ω C = B a (x O ). For this problem, integral representation (2.8) can be conveniently expanded by virtue of (2.7) as

u S (x, ω) = -e k Γa u F (ξ, ω)• tk (ξ, x, ω; n) dΓ ξ u S,F -e k Γa u S (ξ, ω)• tk (ξ, x, ω; n) dΓ ξ u S,S , (4.1)
for x ∈ Ω\B a , where n denotes the unit normal on Γ a outward to the exterior domain. By means of the divergence theorem, the first term can be rewritten as

u S,F (x, ω) = -e k Γa n(ξ)• σk (ξ, x, ω)•u F (ξ, ω) dΓ ξ = e k Ba ∇ ξ • σk (ξ, x, ω) •u F (ξ, ω) + σk (ξ, x, ω) : ∇u F (ξ, ω) dΩ ξ , (4.2)
for x ∈ Ω\B a , where σk = ( σk ) T denotes the elastodynamic stress Green's function,

∇ ξ • σk (ξ, x, ω) + δ(ξ-x)e k = -ρω 2 ûk (ξ, x, ω), x, ξ ∈ Ω, k = 1, 2, 3. (4.3)
To evaluate the limiting form of (4.2) as the cavity size a vanishes, one may invoke the Taylor expansion of ûk inside the void,

ûk (ξ, x, ω) = ûk (x O , x, ω) + (ξ-x O )•∇ ξ ûk (ξ, x, ω)| ξ=x O + . . . , ≡ ûk (x O , x, ω) + ψ(ξ, x O , x, ω), ξ ∈ B a , |ψ| = O(a), a → 0, (4.4)
and write similar expressions for σk , u F , and ∇u F . As a result, on employing (4.3), the limiting behaviour of (4.2) as a → 0 can be reduced to

u S,F = a 3 |B| e k σk (x O , x, ω) : ∇u F (ξ, ω)| ξ=x O -ρω 2 ûk (x O , x, ω)•u F (x O , ω) + η = a 3 |B| e k 1 2µ σk (x O , x, ω) : σ F (x O , ω) - ν 2µ(1+ν) (tr σk (x O , x, ω)) (tr σ F (x O , ω)) -ρ ω 2 ûk (x O , x, ω)•u F (x O , ω) + η(x O , x, ω), x ∈ Ω\B a , (4.5)
where |η| = o(a 3 ) as a → 0, σ F = 1 2 C : (∇u F + ∇ T u F ) is the free-field stress tensor and ν = λ/{2(λ + µ)} denotes the Poisson's ratio of the semi-infinite solid. It should be noted that (4.5) rests on the implicit assumption that x O = x, a hypothesis that will be addressed in the sequel.

Boundary distribution of the scattered field

To elucidate the contribution of the second term u S,S in (4.1) for diminishing a, it is necessary first to determine the asymptotic behaviour of the scattered field, u S , along Γ a in the limit as a → 0. For the latter problem, one has

∇•(C : ∇u S ) = -ρ ω 2 u S , ξ ∈ Ω \ B a , t S = -t F for ξ ∈ Γ a and t S = 0 for ξ ∈ S, (4.6)
where t S and t F are the surface tractions associated respectively with u S and u F . On the basis of (4.6), the regularized BIE for the scattered field (10) can be written as

Γa t F (ξ, ω; n)• ûk (ξ, x, ω) dΓ ξ + Γa (u S (ξ, ω)-u S (x, ω))•[ tk (ξ, x, ω; n)] 1 dΓ ξ + Γa u S (ξ, ω)•[ tk (ξ, x, ω; n)] 2 dΓ ξ = -u S k (x, ω), x ∈ Γ a , k = 1, 2, 3. (4.7)
For the ensuing developments, it is useful to employ (2.6) and rewrite (4.7) as

Γa t F (ξ, ω; n)• [û k (ξ, x, ω)] 1 + [û k (ξ, x, ω)] 2 dΓ ξ + Γa (u S (ξ, ω)-u S (x, ω))• [ tk (ξ, x, ω; n)] 1 + [ tk (ξ, x, ω; n)] 2 dΓ ξ = -u S (x, ω)• e k + Γa [ tk (ξ, x, ω; n)] 2 dΓ ξ , x ∈ Γ a , k = 1, 2, 3. (4.8)
Under the assumption that the sampling point x O is located at a fixed non-zero (but otherwise arbitrary) depth inside the half-space, [û] 1 and [ t] 1 in (4.8) are, for any sufficiently small a, given by the Kelvin's solution so that

[û k (ξ, x, ω)] 1 = O |ξ-x| -1 , [û k (ξ, x, ω)] 2 = O(1) [ tk (ξ, x, ω)] 1 = O |ξ-x| -2 , [ tk (ξ, x, ω)] 2 = O(1), ξ, x ∈ Γ a ⇒ |ξ-x| ≤ 2a, (4.9)
as a → 0. On utilizing (4.9) together with the expansion of the free-field traction

t F (ξ, ω) = n(ξ)• σ F (x O , ω) + (ξ-x O )•∇σ F (ξ, ω)| ξ=x O + . . . , ≡ n(ξ)•σ F (x O , ω) + φ(ξ, x O , ω), |φ| = O(a), ξ ∈ Γ a , a → 0, (4.10)
around ξ = x O , the limiting form of (4.8) as a → 0 can be reduced to

Γa n(ξ)•σ F (x O , ω)•[û k (ξ, x, ω)] 1 dΓ ξ + Γa (v(ξ, ω)-v(x, ω))•[ tk (ξ, x, ω; n)] 1 dΓ ξ = -v k (x, ω), x ∈ Γ a , k = 1, 2, 3, (4.11)
where v denotes the leading asymptotic behaviour of u S along the cavity boundary,

v(ξ, ω) = Asym a→0 u S (ξ, ω), ξ ∈ Γ a . (4.12)
With reference to (2.6) and (4.6), a comparison of (4.7) and (4.11) reveals that v can be interpreted as an elastostatic solution of the exterior problem

∇•(C : ∇v) = 0, ξ ∈ R 3 \ B a , t v = -n•σ F (x O , ω), ξ ∈ Γ a , (4.13)
for the cavity B a in an infinite elastic medium, where t v = n• C : ∇v. In fact, since the prescribed boundary traction t v is defined in terms of a constant stress tensor, the solution to (4.13) can be conveniently recast as

v(ξ, ω) = a σ F k (x 0 , ω) ϑ kl (ζ), ζ = (ξ -x 0 )/a, k, l = 1, 2, 3, (4.14)
in terms of the individual solutions ϑ kl = ϑ lk to six canonical problems

∇ ζ •(C : ∇ ζ ϑ kl ) = 0, ζ ∈ R 3 \ B, t ϑ = -1 2 n•(e k ⊗ e l + e l ⊗ e k ), ζ ∈ ∂B, k, l = 1, 2, 3, (4.15)
for the normalized cavity B in an infinite elastic solid where

t ϑ = n• C : ∇ ζ ϑ kl .
It is useful to note that the reduced problems (4.15) are independent of x 0 and a.

Domain variation of the scattered field

From (4.1), (4.5), (4.12) and (4.14), one finds that the limiting behaviour of (4.1) for vanishing cavity size can be synthesized as

u S (x, ω) = a 3 |B| e k σk (x O , x, ω) : A : σ F (x O , ω) -ρ ω 2 ûk (x O , x, ω)•u F (x O , ω) + ϕ(x O , x, ω), x ∈ Ω\B a , |ϕ| = o(a 3 ), a → 0, (4.16)
with the constant fourth-order tensor A given by

A ijkl = 1 2µ I ijkl - ν 1 + ν δ ij δ kl - 1 |B| ∂B ϑ kl i (ζ) n j (ζ) dΓ ζ . (4.17)

Infinitesimal cavity of spherical shape

When vanishing spherical cavities are considered (that is, when B is a unit ball), the elastostatic exterior problem (4.13) can be recast in terms of spherical harmonics. In this way, an explicit solution to (4.13) can be obtained (see [START_REF] Lurie | Theory of Elasticity[END_REF]) in the form

v(ξ) = a 3 4µ σ F •e r r 2 + r 3 e r × ∇× σ F •e r r 2 - 4-10ν 7-5ν r ∇• σ F •e r r 2 e r + 5ν -1 3(7-5ν) r 2 ∇∇• σ F •e r r 2 + r 2 -a 2 7-5ν ∇∇• σ F •e r r 2 , ξ ∈ R 3 \B a , (4.18)
where σ F ≡ σ F (x O , ω) and e r = e r (ξ) is the unit vector in the r-direction of the spherical coordinate system (r, θ, φ), originated at the centre of the cavity. On employing the identities

∇• r -2 σ F •e r = r -3 [tr σ F -3 σ F : (e r ⊗ e r )] , ∇∇• r -2 σ F •e r = -3r -4 [{tr σ F -5 σ F : (e r ⊗ e r )} e r + 2 σ F •e r ] , (4.19)
e r × ∇× r -2 σ F •e r = -3r -3 [{σ F : (e r ⊗ e r )} e r -σ F •e r ] ,
and the fact that r(ξ) = a and e r (ξ) = -n(ξ), ξ ∈ Γ a , for a spherical cavity, the boundary variation of v in (4.18), and thus that of u S as a → 0 can be reduced to

v(ξ) = - a µ 4-5ν 7-5ν σ F (x O , ω)•n(ξ) - 3-5ν 4(7-5ν) (tr σ F (x O , ω)) n(ξ) , ξ ∈ Γ a . (4.20)
By virtue of (4.20), the normalized solutions ϑ kl to (4.15) can be directly written as which applies when ∂B is a unit sphere, the fourth-order tensor A in (4.17) reduces to

ϑ kl (ζ) = - 1 2µ 4-5ν 7-5ν n l (ζ) e k + n k (ζ) e l - 3-5ν 2(7-5ν) δ kl n(ζ) , ζ ∈ ∂B,
A = 3(1-ν) 2µ(7-5ν) 5 I sym 4 - 1+5ν 2(1+ν) I 2 ⊗ I 2 . (4.21)
On substituting (4.21) into (4.16), the limiting behaviour of (4.1) for a vanishing spherical cavity can be expressed explicitly as

u S (x, ω) = 4πa 3 3 e k 3(1 -ν) 2µ(7-5ν) 5 σk (x O , x, ω) : σ F (x O , ω) - 1+5ν 2(1+ν) (tr σk (x O , x, ω)) (tr σ F (x O , ω)) -ρ ω 2 ûk (x O , x, ω)•u F (x O , ω) + ϕ(x O , x, ω), x ∈ Ω\B a , |ϕ| = o(a 3 ), a → 0. (4.22)
5. Topological derivative for elastic-wave imaging From (2.1) and (2.7), one finds that the perturbation of the cost functional, J , with respect to the creation of a hole in an otherwise intact (semi-infinite) medium can be written as

J (Ω\B a ; f )-J (Ω; f ) = 1 2 M m=1 {u F +u-u obs }•W m •{u F +u-u obs } -{u F -u obs }•W m •{u F -u obs } x=x m = 1 2 M m=1 u•W m •u + 2 Re {u F -u obs }•W m •u x=x m , ( 5.1)
in terms of the scattered field u. By virtue of (3.1), (4.16) and (5.1), the formula for topological derivative that is relevant to elastic wave imaging immediately follows as

T (x O , f ) = lim a→0 1 a 3 |B| M m=1 1 2 u•W m •u + 2 Re {u F -u obs }•W m •u x=x m = M m=1 Re {u F (x m , ω) -u obs (x m , ω)}•W m •e k × σk (x O , x, ω) : A : σ F (x O , ω) -ρ ω 2 ûk (x O , x, ω)•u F (x O , ω) , (5.2)
with the Einstein summation convention assumed over index k = 1, 2, 3. In particular, when B a is a spherical cavity of radius a, expression (5.2) holds with the constant tensor A given by (4.21), and the topological derivative takes the explicit form

T (x O , f ) = M m=1 Re {u F (x m , ω)-u obs (x m , ω)}•W m •e k 3(1 -ν) 2µ(7-5ν) 5 σk (x O , x m , ω) : σ F (x O , ω) - 1 + 5ν 2(1+ν) (tr σk (x O , x m , ω)) (tr σ F (x O , ω)) -ρ ω 2 ûk (x O , x m , ω)•u F (x O , ω) . (5.3)
It is worthwhile noting that on setting ω = 0 and replacing ûk and σk (k = 1, 2, 3) with their elastostatic counterpart, (5.3) can be shown to be in agreement with the results obtained in (8) for the three-dimensional elastostatic case. For testing configurations synthesized in Fig. 1 where the incident seismic field is generated by a point force f acting at x S on the surface of the half-space, (5.3) can be further specialized by writing u F (x, ω) = f j ûj (x, x S , ω), x ∈ Ω\ Ba , (5.4) so that

T (x O , f ) = f l M m=1 Re {f j ûj (x m , x S , ω) -u obs (x m , ω)}•W m •e k × 3(1 -ν) 2µ(7 -5ν) 5 σk (x O , x m , ω) : σl (x O , x S , ω) - 1+5ν 2(1+ν) (tr σk (x O , x m , ω)) (tr σl (x O , x S , ω)) -ρω 2 ûk (x O , x m , ω)• ûl (x O , x S , ω) (5.5)
with the summation convention assumed over indexes j, k, l = 1, 2, 3. Similar to the restriction made in (4.5), here it is assumed that no two of the source point x S , the observation point x m , and the trial point x O are coincident, a hypothesis that is common for this class of imaging problems.

A generalization of (2.1) and (5.2) through (5.5) to multiple seismic sources f q (q = 1, 2, . . . Q) is straightforward and involves external summation in the form of

J f (•) ≡ Q q=1 J (•; f q ) and T f (•) ≡ Q q=1 T (•; f q ).
(5.6)

A noteworthy feature of (5.5) and (5.6) 2 is that they are explicit in terms of the elastodynamic fundamental solution, which makes their computation relatively economical. Moreover, (5.5) and the approach leading to it in fact hold for any geometrical configuration Ω for which the elastodynamic fundamental solutions ûη and ση (η = 1, 2, 3) satisfying suitable boundary conditions on ∂Ω are known. This of course includes the infinite domain Ω = R 3 and the corresponding full-space fundamental solution.

Computational issues and results

For a testing configuration involving Q source points x S , M observation points x m , and a probing grid of N sampling points x O , computation of (5.5) involves M Q + M N + N Q evaluations of the fundamental solution, thus permitting a computationally efficient exploration of the semi-infinite solid Ω for cavities. On evaluating T f (x O ) over a suitable spatial grid spanning the volume of interest, regions in Ω where T f takes the largest negative values constitute possible cavity locations (see [START_REF] Sokolowski | [END_REF] within the framework of shape optimization) and may be used to form an initial guess, in terms of the expected topology, location and shape of the hidden cavity, for a more refined gradient-based minimization. With reference to (5.3) and (5.5), it is also worth noting that T ≡ 0 when u F = u obs , which infers no preferential cavity location.

In the context of shape optimization, an iterative algorithm based entirely on the topological derivative is also proposed (8) where for each iteration, holes are created in the solid wherever values of T fall below a certain (negative) threshold value. In the subsequent iteration, the total field in the cavitated solid serves as an unperturbed, 'free' field. In the present context of elastodynamics, however, evaluation of the equivalent of (5.5) for a cavitated half-space is in principle feasible but requires a computational effort that is several orders of magnitude larger than that associated with (5.5). For this reason, the stand-alone iterative minimization using exclusively topological derivatives is considered to be inferior to the proposed 'hybrid ' method which combines the grid search with gradientbased optimization. In what follows, the usefulness of the topological gradient given by (5.5) and (5.6b) in selecting an initial guess for the optimization-based approach to inverse scattering problems in elastodynamics will be demonstrated through numerical examples. First, the correctness of (5.5) is checked numerically in Section 6.1; with such result, the preliminary imaging of a single spherical cavity and a dual cavity system are presented in Sections 6.2 and 6.3, respectively.

Comparison with BIE approximation

In what follows, a reference will be made to the testing configuration depicted in Fig. 2. The 'true' cavity Ω true C is spherical, of diameter D = 0.4d, and centred at (d, 0, 3d) inside the half-space. In succession, the cavity is illuminated by sixteen point sources f q (q = 1, 2, . . . 16) acting normally on the surface of the semi-infinite solid; for each source location x S , Cartesian components of the ground motion, u obs , are monitored via twenty five sensors x m distributed over the square testing grid. In the absence of physical measurements, experimental observations (u obs ) are simulated using the regularized BIE method [START_REF] Pak | [END_REF] with the surface of the cavity discretized via eight-node quadratic boundary elements [START_REF] Brebbia | Boundary Element Techniques[END_REF].

To provide a focus for the numerical study, all topological gradients computed in the sequel are based on the formula (5.5) which postulates that the infinitesimal cavity is of spherical shape.

For simplicity, the matrices of weighting coefficients, W m (m = 1, 2, . . . M ) are taken as identity operators. In a general imaging situation involving repeated experiments with random measurement errors, W m could be taken as an inverse of the data covariance operator characterizing u obs (x m ) to discriminate between the waveform measurements with strong and poor signal-to-noise ratios [START_REF] Tarantola | Inverse Problem Theory[END_REF].

To validate (5.5) numerically, a use is made of (5.1) so that the multi-source formula for T f (x O ) given by (5.6b) and (5.5) can be compared to its finite-difference approximation,

Tf (x O ; a) = Q q=1 T (x O ; f q , a), where T (x O , f ; a) = 3 8πa 3 M m=1 ũS (x m , ω; x O , a)•W m • ũS (x m , ω; x O , a) + 2 Re {u F (x m , ω)-u obs (x m , ω)}•W m • ũS (x m , ω; x O , a) x=x m . (6.1) In (6.1), u F (x m , ω) = f j ûj (x m , x S , ω) = f 3 û3 (x m , x S , ω
) where x S is the point of action 96 96 of f and ũS is, as an approximation to u, computed via a BIE method by discretizing the surface of the trial spherical cavity Γ a (with fixed small radius a) in terms of quadratic surface elements. As an illustration, a comparison between T f and Tf is made at the interior point x O = (-d, d, 2d) and four excitation frequencies

ω ≡ ωd/c s = k, c s = µ/ρ, k = 1, 2, 4, 8, (6.2)
where c s denotes the shear wave speed in the half-space. From (6.2), one may observe that the ratio between the shear wave length, λ s , and the cavity diameter is approximately equal to λ s /D ≈ 8/k (k = 1, 2, 4, 8), thus characterizing the frequencies selected as those belonging to the so-called resonance region [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] where wave lengths are larger than or of a comparable size to the diameter of the scatterer. It is also worth noting that the resonance ('low-frequency') region is the one most commonly explored in the inverse acoustic and electromagnetic theory. For completeness, the number of boundary elements used for calculating u obs and ũS in (5.5) and (6.1) at each excitation frequency is enclosed in Table 1.

In Fig. 3, the ratio Tf /T f is plotted versus the size of the trial cavity, a. From the display, one may observe that the finite difference results approach the respective 'exact' values given by (5.5) with diminishing a for all four frequencies. It is also worthwhile noting that the difference between Tf and T f drops below 0.2% at the end of each curve. As might be expected, higher excitation frequencies (shorter wave lengths) require in general smaller values of a for an accurate finite difference estimate. The only exception to this rule is a 'faster' convergence of the results for ω = 8 than those for ω = 4; an anomaly apparently induced by the magnitude of T f (x O , ω) that is two decades larger at ω = 8.

Preliminary imaging of a spherical cavity

To examine the effectiveness of (5.5) as a tool for delineating plausible void locations in a semi-infinite solid, the imaging problem in Fig. 2 is taken as an example. For this testing configuration, the values of T f (x O ) are computed over the horizontal surface S h = {ξ ∈ Ω| -5d < ξ 1 < 5d, -3d < ξ 2 < 3d, ξ 3 = 3d} passing through the centroid of the 'true' cavity and plotted in Fig. 4 for the suite of frequencies (6.2). The computational grid is chosen so that the sampling points x O are spaced by 0.25d in both ξ 1 and ξ 2 directions. In for comparison, the true cavity is outlined in white in each of the diagrams. The results indicate the potential of the topological derivative as an approximate, yet computationally efficient tool for exposing the cavity location, with 'higher' frequencies (ω = 2, 4) providing in general better resolution. From the diagram for ω = 8 where λ s /D ≈ 1, however, it is also evident that the topological derivative approach performs best when used in conjunction with wave lengths exceeding the cavity diameter. This is perhaps not surprising, since the assumption of an infinitesimal cavity, implicit to (3.1) and (5.5), is better conformed with by finite cavities that are 'small ' relative to the probing wavelength.

For completeness, the variation of T f (x O ) across the vertical planar region S v = {ξ ∈ Ω|-5d < ξ 1 < 5d, ξ 2 = 0, 0.25d < ξ 3 < 6d} is given in Fig. 5. Similar to the earlier diagram, the sampling points x O are spaced by 0.25d in the ξ 1 and ξ 3 directions. A diminished resolution relative to the previous result reflects the major limitation of the 'experimental' data set, that is, the fact that both source and receiver points are limited to a single planar surface. Following the classical trade-off between the resolution and variance common to this class of inverse problems [START_REF] Menke | Geophysical Data Analysis: Discrete Inverse Theory[END_REF], the contour plots for ω = 2 and 4 exhibit greater accuracy than that for ω = 1, but are also plagued with local minima that are absent in the former diagram. Consistent with the previous assertion, the non-informative distribution of T f for ω = 8 indicates that the use of topological derivative in elastic-wave imaging is most effective at 'low' excitation frequencies, that is, those inside the resonance region.

With diagrams such as those in Figs. 4 and5, an algorithm for identifying plausible cavity ω=1 ω=2 ω=4 ω=8

Fig. 4 Distribution of (µd) -1 T f in the ξ3 = 3d (horizontal) plane: spherical cavity locations could be devised on the basis of the non-zero distribution of an auxiliary function

T f (x O ) = T f (x O ), T f < C, 0, T f ≥ C, (6.3)
where C < 0 denotes a suitable threshold value. With such definition, it is also possible to combine the individual advantages of different probing wavelengths by employing the product of (6.3) at several frequencies. As an illustration of the latter approach, Fig. 6 plots the distribution of the product of T f | ω=1 and T f | ω=2 in the vertical plane, with C set to approximately 40% of the global minima of the respective distributions in Fig. 5. Despite the limited accuracy and multiple minima characterizing respectively the individual solutions for ω = 1 and ω = 2, the combined result stemming from (6.3) points to a single cavity with its centre and size closely approximating the true void configuration.

Multiple void problem

In view of the fundamental assumption underlying (5.5), that is, a single spherical hole (of infinitesimal radius) is created inside the semi-infinite solid, it is important to evaluate the ω=1 ω=4 ω=8 ω=2

Sources Receivers

Fig. 5 Distribution of (µd) -1 T f in the ξ2 = 0 (vertical) plane: spherical cavity performance of the topological derivative approach in situations involving non-spherical or multiple cavities. To this end, consider the constellation given in Fig. 7 where the hidden cavity system (Ω true C ) consists of a spherical void of radius D = 0.4d and an ellipsoidal void whose principal semi-axes, (D 1 , D 2 , D 3 ) = (0.4d, 1.2d, 0.6d), are aligned with the reference Cartesian frame. The two cavities are centred respectively at (-d, d, 2d) and (2d, d, 2d) inside the half-space. Also indicated in the Figure is the boundary element mesh used to simulate the experimental measurements. For comparison with the previous example, the surface testing configuration is assumed to be the same as that in Fig. 2.

Fig. 8 plots the distribution of T f (x O ) across ξ 3 = 2d and ξ 2 = d planes for two representative frequencies, ω = 1 and ω = 2. As a reference, intersection of the 'true' cavity system with the respective cutting planes is outlined in white in each of the graphs. In the top left diagram, the pear-shaped lower frequency contours (ω = 1) for the horizontal plane correctly envelop both cavities, albeit failing to provide a more detailed map of the void system. In contrast, the distribution of T f | ξ3=2d for ω = 2 not only points to the correct number of cavities, but also reasonably approximates their individual shapes. Similar to 0 Fig. 6 Distribution of (µd) -2 T f |ω=1× T f |ω=2 in the ξ2 = 0 plane: spherical cavity T f at each frequency are set to approximately 30% of the global minima of the respective (ξ 3 = 2d) distributions in Fig. 8. Notwithstanding the evident lack of accuracy relative to the images in the horizontal plane, the hybrid contour image in Fig. 9 unequivocally points to two distinct cavities approximating the true void constellation.

Additional considerations

The foregoing examples illustrate the potential of topological derivative as a robust, albeit approximate tool for exposing cavities hidden in a semi-infinite solid from elastic waveform measurements. With such features, the proposed methodology can be regarded as a full complement to the gradient-based BIE imaging algorithm (4), whose high resolution potentials are critically dependent on suitable parametrization, and thus on prior information describing the cavity location, topology, and geometry. To build a comprehensive computational platform for the three-dimensional inverse analysis of voidscattered elastic waves, a hybrid scheme could thus be devised where a computationallyeffective probing tool such as (5.5) is used to explore the volume of interest in a grid-type fashion and thus furnish suitable prior information for the more refined, gradient-based imaging stage.

In situations where an approximate location of the void is unknown beforehand, a construction of its full three-dimensional image in the context of previous examples would entail an inspection of many cutting planes and thus an excessive number of sampling points. To mitigate the problem, an algorithm similar to that proposed for the linear sampling method in acoustics (17) could be adopted. In such iterative scheme, a volume of the host solid under consideration would be initially partitioned into a 2 × 2 × 2 grid of cubes; the topological derivative T f would then be evaluated at their centroids, and each cube where T f falls below a certain threshold value would further be subdivided into eight (2 × 2 × 2) sub-cubes for the next iteration.

Summary

In this communication, the concept of topological derivative that has its origins in elastostatics and shape optimization is extended to three-dimensional elastic-wave imaging of semi-infinite and infinite solids. On taking the limit of the boundary integral representation of the scattered field caused by a spherical cavity with diminishing radius, the topological derivative, which quantifies the sensitivity of the featured cost functional due to the creation of an infinitesimal hole, is formulated explicitly in terms of the elastodynamic fundamental solution. A set of numerical examples involving a planar testing surface on top of the semi-infinite solid is included to validate and illustrate the proposed imaging algorithm wherein plausible void regions are delineated through negative values of the topological derivative at sampling points. The results for both single and multiple cavity configurations indicate that the approach is most effective when used at frequencies inside the resonance region, that is, with wave lengths exceeding the cavity diameter. It is also found that the use of multiple excitation frequencies as a tool to illuminate the cavity may enrich the experimental data set and thus mitigate the geometric limitations of the assumed planar testing surface. To the authors' knowledge, an application of topological derivative to inverse scattering problems has not been attempted before; the approach proposed in this study can be extended to deal with acoustic problems as well.

Fig. 1

 1 Fig. 1 Illumination of an underground cavity (Ω true C ) by elastic waves

  along the boundary of B, where n(ζ) is the normal outward to R 3 \ B. With such result and the identity ∂B n ⊗n dΓ ζ = 4π 3 I 2 ,

Fig. 2

 2 Fig. 2 Sample imaging problem

Fig. 3

 3 Fig. 3 Topological derivative estimates: finite difference versus explicit formula

Fig. 7 3 Fig. 8

 738 Fig. 7 Geometry and discretization of the multiple void domain in the half-space (ξ3 > 0)

Fig. 9

 9 Fig. 9 Distribution of (µd) -2 T f |ω=1× T f |ω=2 in the ξ2 = d plane: dual cavity problem

Table 1

 1 Number of surface elements in the BIE solution used for calculating u obs and ũS

	ω	1.0 2.0 4.0 8.0
	u obs	216 216 384 384
	ũS | a=d/10	96
	ũS | a=d/20	96	96
	ũS | a=d/40	96	96	96	96
	ũS | a=d/80	96	96	96	96
	ũS | a=d/160		96	96	96
	ũS | a=d/320

Acknowledgment

The support provided by the National Science Foundation through CAREER Award No. CMS-9875495 to B. Guzina and the University of Minnesota Supercomputing Institute during the course of this investigation is kindly acknowledged. Special thanks are extended to MTS Systems Corporation for providing the opportunity for M. Bonnet to visit the University of Minnesota through the MTS Visiting Professorship of Geomechanics.