A divide-alternate-and-conquer approach for localization and shape identification of multiple scatterers in heterogeneous media using dynamic XFEM

A numerical method for localization and identification of 
multiple arbitrarily-shaped scatterers (cracks, voids, or 
inclusions) embedded within heterogeneous linear 
elastic media is described. An elastodynamic implementation of the extended finite element 
method (XFEM), which is endowed with a spline-based parameterization of the scatterer boundaries, is employed to solve the forward (wave propagation) problem. This particular combination enables direct, sensitivity-based, and computationally efficient manipulation of the scatterers' boundaries over a stationary background mesh during the inversion process. The inverse problem is cast as a formal optimization problem 
whereby the discrepancy between the measured wave responses and those from the estimated 
scatterers is minimized. The solution is achieved through a gradient-based 
 procedure that is steered by a divide-alternate-and-conquer strategy. The divide-and-conquer 
segment of the search algorithm seeks the global minimizer among potentially multiple 
solutions, whereas the alternate-and-conquer segment adaptively refines the shapes of identified scatterers. The results of 
several synthetic experiments with various types of scatterers are provided. These experiments verify the overall approach, and demonstrate that it is robust, accurate, and effective even at high levels of measurement noise.


(Communicated by Jun Zou)
Abstract. A numerical method for localization and identification of multiple arbitrarily-shaped scatterers (cracks, voids, or inclusions) embedded within heterogeneous linear elastic media is described. An elastodynamic implementation of the extended finite element method (XFEM), which is endowed with a spline-based parameterization of the scatterer boundaries, is employed to solve the forward (wave propagation) problem. This particular combination enables direct, sensitivity-based, and computationally efficient manipulation of the scatterers' boundaries over a stationary background mesh during the inversion process. The inverse problem is cast as a formal optimization problem whereby the discrepancy between the measured wave responses and those from the estimated scatterers is minimized. The solution is achieved through a gradient-based procedure that is steered by a divide-alternate-and-conquer strategy. The divide-and-conquer segment of the search algorithm seeks the global minimizer among potentially multiple solutions, whereas the alternateand-conquer segment adaptively refines the shapes of identified scatterers. The results of several synthetic experiments with various types of scatterers are provided. These experiments verify the overall approach, and demonstrate that it is robust, accurate, and effective even at high levels of measurement noise.
1. Introduction. Identification and quantification of hidden scatterers (e.g., damage or defects such as cracks, voids, and inclusions) at an incipient stage before fracture in structures is an integral part in predictive maintenance as well as structural health monitoring (SHM) [24,16]. Ultrasonic inspection-one of the wellestablished nondestructive evaluation (NDE) methods for assessing the current state of a structural system-has been widely applied in various areas such as SHM [24], medical imaging [19,54] and geophysical prospecting [28,57]. In a typical ultrasonic inspection, a single actuator and multiple sensors are positioned along the boundary of the host medium. Mechanical waves produced by the actuators travel through the medium and are reflected/refracted when encountering internal scatterers. From the wave signals returning back to each sensor, their locations and sizes are identified. Such an approach utilizing ultrasonic waves, however, has difficulties in the precise reconstruction of multiple scatterers within a heterogeneous medium [46]. This is primarily due to the interference/superposition of multiply scattered waves from the scatterers [57] and the wave reflections from the internal interfaces between two media with different types of property [59]. For accurate estimation of buried scatterers, more sophisticated and systematic identification approaches based on the finite element method (FEM) and the boundary element method (BEM) have recently been proposed for different application areas (see, for example, [58,8,9,31]).
Another class of methods to localization and identification of scatterers/inclusions comprises the sampling-type approaches [42], which include the Linear Sampling Method (LSM) [15,11,10], the Multiple Signal Classification (MUSIC) algorithm [14,2,1,13], and the Direct Sampling Method (DSM) [25,26,34,27]. LSM [15] is a computationally efficient approach (by being based on linear inversions) and needs no a priori information about the target scatterer. The MUSIC algorithm [14] is a non-iterative imaging method, which has been developed and widely utilized in numerous applications (e.g., to identify internal corrosion of a pipe [2], and inclusions embedded in half-space [1]). DSM is a more recent technique [25] that has been devised primary for electromagnetic and acoustic waves. It offers computational efficiency over both LSM and MUSIC, in that it does not require any matrix operations; and it can accurately reconstruct scatterers/heterogeneities using limited information with very few incident waves. Despite their advantages, these methods are generally semi-analytical in nature and require exact solutions (or Green's functions) of wave propagation; as such, handling arbitrary domain geometries and boundary conditions as well as heterogeneous/anisotropic host media can become significant challenges.
In general, parameter identification (here, inverse-scattering) problems can be cast as minimization problems, wherein an objective (aka cost or error ) functional that quantifies the discrepancy between predictions of measurements from a forward simulation run with a parametric model, and the measurements-e.g., displacements obtained via ultrasonic testing-themselves [35]. The said identification process (here, seeking first the location, and subsequently, the optimal shape parameters of the scatterers) requires two main ingredients: (i) an accurate forward modeling method; and (ii) a robust optimization algorithm. The forward modeling predictions of measurements can be obtained for a given set of iteratively estimated scatterer (aka updating) parameters. In the inversion process, these parametersdescribing the shape, size, and location of each scatterer-are updated using a proper optimization method that steers them towards the minimizer(s) of the objective functional.
A classical challenge in shape optimization and identification methods is the need to employ an adequately general method for forward simulations-so that, for example, arbitrary loading and boundary conditions, material heterogeneities, and different sensor-actuator deployment schemes can be handled. As such, the classical finite element method, at first, appears to be suitable choice for generating forward simulations. However, iterative refinement of the location and the shape of the scatterer(s) pose a significant challenge, in that the classical FEM would require re-meshing of the entire domain for every iteration. This issue is not only a challenge because of potentially intractable computational costs, but also because of lack of robust and adequately general automatic re-meshing methods for arbitrary domain geometries and topologies, as well as element types [17]. In this thesis, this critical issue will be circumvented through the adoption of the so-called eXtended Finite Element Method (XFEM), which obviates re-meshing during the minimization procedure.
Since its inception by Belytschko and Black in [4], XFEM has been continuously improved, and utilized in a variety of applications, including the modeling of crack propagation [45], material interfaces and dislocations [5], multi-phase and free-surface flows [48], bone fracture [20], and biofilms [49] to name a few. As stated above, XFEM is very well-suited for shape identification and optimization problems because (i) re-meshing procedures are not required to render the background mesh to conform perfectly to the arbitrarily evolving boundary of a scatterer during an iterative identification procedure [6,30]; and (ii) the host medium can be arbitrarily heterogeneous and anisotropic [50].
A key issue in any generic inverse problem is to filter out multiple erroneous solutions (local minima), which inherently originate due to temporal and spatial sparseness of the measurements. A remarkably fruitful avenue in that regard has been the use of Genetic Algorithms (GAs) [23] and other nature-inspired techniques such as the Artificial Bee Colony (ABC) algorithm [32]. This is because such heuristic algorithms facilitate an almost exhaustive search 1 ; and they can be applied to complex (nonlinear, multi-dimensional) optimization problems-including inverse scattering [18]-with ease.
In recent years, and pertinent to the present study, [43] proposed a combined XFEM-GA approach for single or multiple scatterers with regular geometry (i.e., a straight line crack, a circle, or an ellipse). The same approach of combining XFEM with heuristic minimization algorithms continued to the present day with various algorithmic advances and some experimental validation [44,56,12,51].
While successful, the aforementioned XFEM-GA algorithms entail a large number of forward solves in order to converge on the optimum values; and the computational cost exponentially increases with the number of updating parameters. Moreover, existing literature either has only dealt with simple inclusions shapes, or elasto-static data.
In the present paper, an identification methodology of internal multiple arbitraryshaped scatterers with several different types of geometries (e.g., a crack, a void or an inclusion) is proposed, which carries out an extension of early our studies [30,29]. Herein, the methodology to the parametric identification problems can be largely divided into two main components: (i) the dynamic XFEM coupled with a cubic spline-based modeling which is capable of describing the arbitrary shape of scatterers; and (ii) a gradient-based optimization under a divide-alternate-and conquer approach for efficiently eliminating multiple local minima and precisely delineating complex geometries of multiple scatterers. Noteworthy features in this study compared to previous work [43,44,56,12,51]-wherein XFEM served as a forward solver to inverse scattering problems-are that: (i) a pointwise-parametric cubic spline method is applied to a dynamic XFEM framework to accurately represent open/closed forms of scatterers (e.g., curves/concave shapes); (ii) a divide-conquerand-alternate approach is adopted to enhance the effectiveness of the gradient-based optimization; and (iii) the host media can be heterogeneous; and (iv) in order to avoid local minima, a set of independent inversions are carried out first for initially The reminder of the paper is organized as follows. In §2, we briefly explain the dynamic XFEM formulation for a crack, a void, and an inclusion as well as the cubic spline modeling for arbitrary geometries. §3 provides the description of the gradient-based optimization method, and discusses the divide-alternate-andconquer strategy. §4 gives several numerical experiments to testify the performance of the proposed methodology. Finally, concluding remarks are presented in §5.

2.
Dynamic XFEM with cubic splines for solving the forward problem.

2.1.
Equations governing wave motion. We consider an inverse scattering problem in which a linear elastic heterogeneous host medium occupies the open set Ω ⊂ R 2 bounded by Γ, as shown in Figure 1, such that where Γ u and Γ t denote the specified displacement and traction along the outer boundary of Ω, respectively; Γ d , the crack surfaces; Γ h , the void surfaces; and Γ I , the material interfaces of inclusions.
The strong form of the governing equations describing elastic wave motion and the initial conditions are u(x, 0) = 0 in Ω, where ∇· denotes the divergence operator; σ := σ(x, t), the Cauchy stress tensor; x, the spatial position; t, time; b := b(x, t), the body force per unit volume; ρ := ρ(x), the mass density; u := [u 1 (x, t) u 2 (x, t)] T , the displacement vector; and( ) and( ) indicate the first and second partial derivatives of the state variables with respect to time t, respectively. The associated boundary conditions are given as σ · n = 0 on Γ d and/or Γ h , [|σ · n|] = 0 on Γ I , where n is the outward unit vector normal on the boundary of Ω;ū andt are, respectively, the prescribed displacements and tractions on Γ u and Γ t ; the displacement and traction boundary conditions are, respectively, held along Γ u and Γ t ; the traction-free conditions are imposed on Γ d and Γ h ; and the interface traction continuity conditions are assumed over Γ I . The constitutive relationships for linear elastic media are defined as where C is the material elasticity tensor; ǫ is the strain tensor; and ∇ is the gradient operator.
2.2. The weak form of the governing equations and its spatial discretization. Multiplied by δu ∈ H 1 (Ω), the weak form corresponding to Eq. (2) is given as where δu are the test functions that belong to H 1 (Ω)-i.e., the first-order Sobolevspace continuous functions-except for the essential boundary Γ u in Ω. The first term of the above expression can be integrated by parts to give By the divergence theorem in Ω and the boundary conditions in Eqs (6,7,8), the first term on the left-hand side of Eq. (11) becomes Substituting Eqs. (11,12) into Eq. (10), the discretized weak form of linear elastodynamics for u h ∈ U h and δu h ∈ V h therefore leads to where H 1h ⊆ H 1 denote finite dimensional Hilbert spaces consisting of the shape functions; u h and δu h are the finite element approximations of the trial and test functions.
We note here that this discretized weak form of the linear elastodynamic problem (a.k.a., the initial-boundary value problem) embodies a well-posed set of boundary conditions in Eqs. (5)(6)(7)(8), so that it yields a unique solution with respect to the unknown displacement field u(x, t) [7]. A general treatment of well-posedness-of second-order damped systems-may be found in [3].
2.3. Time-stepping. Neglecting damping, the semi-discrete form of the linear dynamic equations (13) at the n-th time step is where M and K are the global mass and stiffness matrices; U n ,Ü n , f n are the nodal displacement, acceleration, and external load vectors at the n-th time step. The Newmark time-stepping method, originated by N.M. Newmark [41], is adopted to integrate Eq. (16). Newmark directly derived and expressed the state of velocities and displacements at the n-th time step, based on Taylor's serieṡ By the assumption of linear acceleration, i.e., ... Un = (Ü n −Ü n−1 )/∆t , the standard Newmark formula, which specify the state of variables at the n-th time step, given the prior state, are derived aṡ By substitution of the above two-state approximation into Eq. (16), a system of equations for the acceleration is obtained as below Here, the Newmark parameters are set at α = 0.5 and β = 0.25 in subsequent calculations, which is generally preferred for long time simulations because it produces unconditionally stable solution regardless of a time step size ∆t. Moreover, this method offers second-order accuracy.
2.4. The XFEM approximation for cracks. Based on the concept of partition of unity [39], the XFEM-approximated displacement field for general cracks is decomposed into two parts: (i) a conventional finite element approximation part (denoted by the summation with respect to i ∈ I), and (ii) an enrichment approximation part (denoted by the summations with respect to j ∈ I, and k ∈ I Λ ) [22] as Here, I denotes the set of all nodes over computational domain (including the cracks, voids, and inclusions); I Γ , the set of enriched element nodes that are fully cut/traversed by the cracks Γ d (marked with red-filled circles in Figure 2(a)); I Λ , the set of enriched nodes of the elements containing the crack tips (marked with blue-filled squares in Figure 2(a)); N (i, j or k) (x), the standard finite element shape functions associated with the (i, j or k)-th node at the location of x (i, j or k) ; u i , the standard nodal displacement vectors at node i; b j and c α k , the enriched nodal displacement vectors involving the discontinuous function ψ sign (x) and the asymptotic crack-tip function F α (x) at node j and k, respectively.
The discontinuity in the interior of the cracks is described through the signed distance level set function ψ sign (x) (see, Figure 3), whose value is defined as with where · denotes the Euclidean norm.
The structural behavior around the crack tips in an isotropic elastic medium can be characterized by using crack-tip enrichment functions, which are defined by the asymptotic displacements given as [21], where r and θ are the polar coordinates in the local crack-tip coordinate system shown in Figure 4.
2.5. The XFEM approximation for voids. The XFEM approximation for voids can be simply made out using the enrichment function V (x) [50], where As illustrated in Figure 2 elements inside voids (♦). Thus, the XFEM solution can be obtained simply by numerical integration of element matrices in the same way as in conventional FEM over Ω, but by omitting integration points the inside voids (♦).
2.6. The XFEM approximation for inclusions. For the modeling of inclusions, the XFEM displacement approximation function suggested by Moës et al. [40] can be used, which employs the level set function ζ, as in, where  Note that for voids, the level set function ζ is used as well, in order to distinguish among those nodes that are inside and outside the voids.
Here, the level set function ζ (see, Figures 2(c), 5) play a crucial role in the representation of closed geometries (i.e., inclusions and voids). As this is closely related to the discussion of parametric cubic splines, the determination of the sign and the magnitude of ζ is deferred to §3.3 and §3.4, respectively.
3. Spline-based representation of arbitrary shapes in XFEM. In this section, we briefly explain how to construct the open/closed geometries of arbitrarilyshaped scatterers (cracks/voids or inclusions) using cubic splines, based on our earlier works [30,36,37,29]. In particular, we describe a parametric cubic spline method that allows a closed geometry to be constructed with a few number of control points and shape parameters, which is a highly favorable attribute for the optimization problems at hand.

3.1.
Shape representation for open curves. Cubic splines are used for representing curved cracks with a minimal number of discrete control points/segments. As shown in Figure 6, each discrete curve segment (marked with the numbered squares) is defined through control points (marked with purple-filled circles). Each control point carries three shape parameters-viz., the x and y coordinates and the curve slopes (θ). For the generic i-th curve segment, the cubic interpolating spline function is given by  where the coefficients can be determined as 3.2. Shape representation for closed curves. Parametric cubic splines are used for the modeling the boundary of an arbitrarily shaped void or inclusion. This approach offers the following two advantages: (i) the variables x and y are uncoupled from each other (through the use of a curvilinear coordinate system), which makes it easier to construct complicated geometries compared to conventional cubic spline methods; (ii) each individual control point of the cubic spline segments has two shape parameters, x, y coordinates, excluding θ.
The parametric cubic spline functions are defined with respect to an independent variable s ∈ [0, 1] for x i and y i along the i-th curve segment as where s is equally spaced by (1/n) over s ∈ [0, 1] for simplification; s 0 indicates the value of s at the starting point of each segment; n is the total number of control points. In the closed-loop case (see, Figure 7), an n + 1-th control point is added in the sequence as in i = {1, 2, . . . , n, 1}, and P n+1 = P 1 , where P i denotes x i and y i in each set of cubic spline functions, respectively. The coefficients corresponding to the cubic spline functions x i (s) and y i (s) are given by, where h is the spacing between segments, which is chosen as h = 1/n for convenience. Also, the second derivatives (P ′′ i ) are determined by inverting the following (n + 1) × (n + 1) matrix Here, the sub-matrix for rows i = 2 to m takes a tridiagonal matrix form.

3.3.
Determination of the sign of the level set function ζ. In order to determine the sign (+ or −) of the level set function ζ, which differentiates between the outside and inside of voids/inclusions with closed curves, we devised a modified crossing number method. In the conventional crossing number method, the number of times that a ray-drawn from a point, previously known to be either inside or outside of a closed shape, to a target point-crosses the boundaries of the geometry is counted [53]. After that, from the total crossing number, the sign at the target location is determined. If the crossing number is even (odd), then the target point is on the same (opposite) side of the geometry. However, in this approach, we must recognize a priori that the given target point's location is either inside or outside the geometry. Hence, a novel modified crossing number method is proposed herein, making up for the weakness of the conventional approach, as well as localizing the calculation area near the object. Figure 8 illustrates how the signs at points ζ 1−3 , located on the inside or outside of the geometry, are determined, which requires the following steps: (i) Set the centroid o c of the closed geometry from the given informationi.e., control points (purple-filled circles) and cubic curve segments (numbered squares). (ii) Draw a ray passing through each target point in question (here, ζ 1−3 ).
It is useful to note here that (a) the crossing number is not counted, as an exception for the case of ζ 4 , wherein the ray goes through the tangent point o 3 lying on the curve segment 9 ; and (b) the centroid o c is assumed not to be situated on the boundary of the object, which can be checked by using the parametric cubic spline functions that represent each curve segment. Table 1. Determining the sign at points ζ 1−3 . Figure 9. A plot of the sign of ζ for arbitrarily shaped geometries using the modified ray-crossing number method; bluecolored/empty areas indicate the negative/positive sign of ζ. Figure 9 displays examples of the sign of ζ for arbitrarily shaped geometries using the modified ray-crossing number method for which the blue-and white-colored areas correspond to negative and positive ζ values, respectively.

3.4.
Determination of the magnitude of the level set function ζ. In the present study, the magnitude of the level set function ζ is defined by the shortest distance from any point within the computational domain to the surface of the targeted two-dimensional void/inclusion formed by cubic spline curves.
From §3.2, the parametric cubic spline function for any point (x Γ h/I , y Γ h/I ) within a certain curve segment along the boundary of Γ h or Γ I is given by The distance between (x Γ h/I , y Γ h/I ) and any given point (x, y) is Here, the shortest distance is the optimal value of s obtained by minimizing d(s). To that end, Eq. (35) is substituted into Eq. (36), and the resulting equation is differentiating with respect to s.
Note that the correct answer can be sifted out from the multiple solutions to Eq. (37) by using the specific range of s in each curve segment, which is bracketed by each s 0 , which denote the starting points of the curve segments.

4.
A gradient-based minimization method enhanced with robust search algorithms for solving the inverse problem.
4.1. The objective functional. Identification problems can be cast as minimization problems, wherein properly constructed objective functionals map morphological differences (size, shape, and location) between the estimated and actual (target) scatterers into quantitative values. The objective functional adopted here is defined as (see, also [30]), where u D (x s i , t; ξ) and u M (x s i , t), respectively, indicate the computed and measured wave responses-here, the displacements recorded at the i-th sensor's location x s icorresponding to the estimated and actual 2 scatterers. The term ξ denotes the set of unknown shape parameters; N s , the number of sensors; T , the observation duration; · , the Euclidean norm of a vector.
The algorithm for minimizing the objective functional given in Eq.(38) is a gradient-based search that is enhanced with a back-tracking line-search method, as well as a divide-and-conquer strategy. A detailed exposition of these algorithms is omitted here for brevity, and may be found in [30,29] wherein those methods were utilized for localization and shape identification of a single scatterer. The following section describes a more comprehensive search algorithm to localize and, subsequently, to identify the shapes of multiple scatterers.

4.2.
A divide-alternate-and conquer strategy for multiple inclusions. For a single arbitrarily shaped scatterer, Jung et al. [30] suggested a two-phase divideand-conquer methodology. In the first phase, the approximate location of the single scatterer was determined by using a minimal number of unknown shape parameters; assuming that the type of the scatterer is known a priori, a straight line or a circle was used initially for an arbitrarily shaped crack or void/inclusion, respectively. This first phase was carried out with the aforementioned simple geometries with multiple initial estimates that were distributed uniformly within the host domain. In the second phase, only the successful localization result(s) from the first phase were used, and the simple shapes were replaced with adaptively refined spline-based shapes for shape identification. The adaptive refinement (i.e., the insertion of new control points between existing ones) was carried out until either a desired error norm was attained, or the total number of control points exceeded a pre-determined maximum.
In the present study, the aforementioned two-phase divide-and-conquer strategy is modified and extended to identify multiple scatterers with arbitrary shapes. We retain the divide-and-conquer approach of the first phase-here, to localize multiple scatterers. In the second phase, we employ an alternate-and-conquer approach to delineate the boundary of each scatterer using cubic splines. Details of the new X Y Figure 10.
It should be noted here that the finite element mesh size must be smaller than the shortest wavelength that can be calculated using the material properties of the host media and the frequency range of the excitation source. Naturally, the mesh Figure 11. Flowchart for identifying multiple arbitrarily shaped scatterers: the left and right columns on the flowchart are involved in the divide-and-conquer and alternate-and-conquer strategy, respectively.
size and the frequency-content of the source excitation have effects on the sizes of detectable target scatterers.

5.2.
The divide-alternate-and conquer strategy. The flowchart for the dividealternate-and-conquer approach is shown in Figure 11. As seen, it comprises two main procedures: the left portion of the flowchart describes the divide-and-conquer procedure, which is used for roughly localizing multiple scatterers. The alternateand-conquer procedure (shown on the right) follows this first phase, and sequentially and adaptively identifies/refines the shapes of the now localized scatterers. It is noted that the convergence criteria C I , C II , and C III ) can be devised depending on material systems, inclusions types, as well as computational cost, and are defined as,  Figure 12. An illustration of the divide-and-conquer approach: eighteen uniformly-spaced initial estimates (circles) are employed to seek multiple voids (marked by black dashed lines); some of the initial estimates (blue dashed-dotted lines marked with numbers 3, 6, 10, 13 and 16) converged near the global minima; the others (orange solid lines) fell into local minima. Also, each arrow (→) indicates their converged locations.
with P (n) = 2P (n−1) − 1 with P 0 = 2 for a curved crack, 2P (n− 1) with P 0 = 1, P 1 = 4 for a void/inclusion. (41) where P (n) denotes the total number of shape control points at the n-th step (n = n+1 as passing by the process of adding shape control points to the current selected scatterer's set in the flowchart). The new control points are positioned along the cubic spline curves connecting between the existing control points.
In order to localize the multiple unknown scatterers, the divide-and-conquer strategy is used as the first phase. This procedure is similar to that we used in our earlier study for detecting a single scatterer [30]. The validity of the said approach for the multiple scatterers is evaluated here through the a generic (yet adequately general) example below, and subsequently applied to several numerical experiments on various types of scatterers.
As shown in Figure 12, multiple initial estimates 1 -18 (orange/blue-colored circles) are uniformly assigned within the computational domain first; and the minimization process for each estimate is performed independently until the convergence criterion I (C I , see Eq. (40)) is satisfied. In this first step, the circle shapes of the initial estimates are kept, and only their radii and centroids are used as minimization parameters. The converged/final locations of the initial estimates for each independent inversion are displayed with arrows (→) in Figure 12 where the blue dashed-dotted estimates (i.e., 3 , 6 , 10 , 13 and 16 ) converged to the same (nearly identical) position near each target scatterer (global minima) while the others (orange-colored circles in the figure) reached various other locations (local minima).
At this point, the objective functional values of all converged estimates can be compared. The histogram in Figure 13 displays those values, and visualizes the stark differences of the normalized objective functional values between the cases that converged to the target locations (colored blue) and those that fell into local minima (colored orange). The initial estimates that converged near the target scatterers (the black dashed lines) have distinctly lower values of L than the rest. Therefore, we can effectively select the candidate cases for the second phase.
In the second phase (the right column of the flowchart on Figure 11), we optimize the shape of each candidate scatterer using cubic splines under an alternate-andconquer strategy. At the start of this phase, total number of shape control points (for the case of voids/inclusions) is quadrupled from P 0 = 1 (i.e., the centre point of the circle) to P 1 = 4 (i.e., 4 points equally spaced points along the perimeter of the circle) through the procedure in the flowchart that adds shape control points to the current chosen scatterer set. In subsequent iterations, the previous number of control points are doubled (i.e., P 2 = 8). Figure 14 illustrates progression of the alternate-and-conquer procedure for the example problem. The variation of the value of normalized objective functional with iterations is shown in Figure 15, where it can be seen that at the beginning of each alternation, this value drops dramatically, and reaches a state of diminishing returns (that is, without adding new points, or improving the-presently frozen-estimates of the other scatterers, no significant improvement on the objective functional can be achieved). Ultimately, it is seen that the alternate-and-conquer strategy works well and Figure 14 shows that multiple voids are successfully identified, even though small peaks occur during the alternation step (e.g., between (a) and (b) or (b) and (c) in Figure 14) due to  Remarks. It is useful to note that each of "iteration" that is referred to in all of the figures above includes multiple of forward wave propagation analyses. Specifically, each iteration requires forward analyses twice the current total number of shape parameters plus (ς) due to the finite difference approximations of the gradient and the backtracking line-searches, respectively. For example, in the case of an void/inclusion with 4 control points (i.e., 8 shape parameters), (2 × 8 + ς)-times forward analyses are carried out in a single iteration.

5.3.
Numerical experiments with various types of multiple scatterers in homogeneous and bi-material hosts. Numerical experiments are performed to investigate the performance and effectiveness of the proposed identification approach for various types of multiple arbitrary-shaped scatterers (cracks, voids, or inclusions) within a homogeneous or heterogeneous medium, as shown in Figures 16,  17, and 18. These figures visualize the sequence of estimations of the scatterers' shapes: the figures-(a) display the localization of the estimated scatterers obtained using the divide-and-conquer approach (phase I); and the other figures-(b)-(d)/(e) display the iteratively refined shape estimates (approximated using cubic splines) through the alternate-and-conquer approach.
As explained in previous section for the divide-and-conquer approach, the starting geometries of initial estimates (e.g., a straight line for a crack and a circle for a void/inclusion) are chosen based on the anticipated type of the target scatterers that can be determined through, for example, bulk ultrasonic testing [33,52,38].  is slightly altered by adding the new ones. However, convergence is eventually achieved in all cases.

5.4.
Analysis of the effects of noise on the identification results. We used synthetic data (u M ) instead of real measured responses in the preceding numerical studies, and measurement and modeling errors were not taken into account. In this section, we investigate the effect of measurement errors (noise) on the identification results. In the present study, we use scatterers that are smaller than those before to make it more challenging for the inverse estimation method described above. Figure 20 displays the test problem, which has multiple voids with arbitrary shapes. Again, they are localized first as seen in Figure 20  (c) (d) Figure 20. Evolutionary identification process for multiple voids (dashed lines) from (a) to (d): (a) the initial-estimated scatterers (solid lines) are localized by using the divide-and-conquer strategy; (b)-(d) the scatterers's shapes are refined by using cubic splines under the alternate-and-conquer strategy. Note that a random noise with 10% of the root-mean-square value of the pure/noiseless signal response is injected into the measurements (u M ).
inout data (u M ) is first polluted with a white Gaussian random noise with the standard deviation whose amplitude is 10% of the root mean square of pure/noisefree wave responses. This level of measurement error is generally deemed high, and thus adequate, to test the robustness of estimation algorithms. Figure 21 shows the convergence behavior of Λ (i.e., Λ i = L i /L 0 ) in each subprocess ((b)-(d)) of the alternate-and-conquer approach, wherein the value of Λ during the stage between (b) and (c) sharply declines, whereas its value in stage (d) remains almost constant, which is due to the effect of the injected noise. Nevertheless, in spite of the relatively high level of measurement noise, the final converged geometries are fairly accurate as seen in Figure 20-(d). (b) (c) (d) Figure 21. The variation of Λ with the iteration number, i, during the detection process (b) to (d) in Figure 20 for multiple voids.

Iteration number
Here, Λ i = L i /L 0 .
6. Summary and conclusions. An effective algorithm to localize multiple scatters within a heterogeneous medium and to identify their shapes was presented. The computational engine for the said algorithm is cubic-spline-based XFEM, which allows accurate representation of the iteratively varied geometries of the unknown/estimated scatterers without re-meshing. Parameters of the cubic spline serve as the updating parameters for the inverse problem, which is set up as an optimization problem whose objective functional measures the discrepancy between measured and predicted wave responses in the time domain. The proposed search algorithm comprises two phases. In the first (i.e., divideand-conquer) phase, the approximate locations and shapes of the multiple scatterers are obtained. In the second (alternate-and-conquer) phase, the estimated shapes of the now localized scatterers are iteratively refined one at a time. The minimization problem is solved with a gradient-based method that is safeguarded with a backtracking line-search technique.
Numerical experiments on various types of multiple arbitrarily shaped scatterers clearly demonstrated the effectiveness and accuracy of proposed approach, which was also shown to be robust against relatively high levels of measurement noise.
The primary novel aspects of this study compared to earlier works are the development of (i) a modified ray-crossing method applied to parametric cubic splines, which effectively delineates arbitrarily varied shapes of the scatterers during the iterative optimization process; and (ii) an enhanced divide-conquer-and-alternate approach, which consists of a (sub-) gradient-based optimization method for localization and shape identification of multiple unknown scatterers within a heterogeneous host.
A potential extension of the approach presented here would be to three-dimensional inverse scattering problems. In such an extension, the core issue would be to devise a method to parametrically describe a three-dimensional scatterer's surface. While numerous alternatives exist, the so-called NURBS appear to offer a viable pathway toward such an extension [47].