Defect reconstruction in a 2D semi-analytical waveguide model via derivative-based optimization

This paper considers the reconstruction of a defect in a two-dimensional waveguide during non-destructive ultrasonic inspection using a derivative-based optimization approach. The propagation of the mechanical waves is simulated by the Scaled Boundary Finite Element Method (SBFEM) that builds on a semi-analytical approach. The simulated data is then fitted to a given set of data describing the reflection of a defect to be reconstructed. For this purpose, we apply an iteratively regularized Gauss-Newton method in combination with algorithmic differentiation to provide the required derivative information accurately and efficiently. We present numerical results for three different kinds of defects, namely a crack, a delamination, and a corrosion. These examples show that the parameterization of the defect can be reconstructed efficiently and robustly in the presence of noise.


Introduction
Ultrasonic guided waves, especially Lamb-waves, are used in Non-Destructive Testing (NDT) and Structural Health Monitoring (SHM) to identify defects. Here the long testing range that can be covered using guided waves is one of the main advantages over conventional scanning techniques. However, this advantage also yields new challenges. The multi-modal nature of guided waves with different wave velocities and the large distance between the defect and the sensor makes the localisation and characterisation of defects more difficult. On the other side, Lamb-waves provide the flexibility to choose the most appropriate mode, which is the most sensitive to a specific type of defect.
Most models for crack identification assume that the data containing the waves scattered by the defect are available at certain sensor positions. Depending on the type of sensor, this data can be an electrical signal that represents a displacement or velocity field. We will refer to this data as measurement data, even though for many publications including this one, it is generated artificially by appropriate simulations. There are two main types of methods to detect a defect in a waveguide. One approach is based on imaging to localise and characterise the defects. Another class of algorithms fits a damage model directly with the measured data. The resulting damage model can then be further analyzed, e.g., using fracture mechanics techniques. This also supports cross-section view z surface view z x x y x y defect Sensor Figure 1: Different views for a 2d-simplification of the 3d-problem other goals like NDT and SHM investigations and may answer the question of whether the damage is severe leading to a failure of the structure. Most imaging algorithms assume a sensor network around the defect and consider a surface view of the waveguide as sketched in Figure 1. Often these approaches are called guided wave tomography or migration in the seismic community. The measurement data is then propagated backwards from the sensors on the surface. As a possible method of backward propagation, each pixel of the surface is analyzed by the direct sound path between this pixel and the transmitter/receiver pairs [IC08]. This approach relies on the assumption that there is only one dominant wave mode to use its group velocity for the backward propagation. Other methods do not need a dominating mode, for example, if a time-reversal approach is utilized as backward propagation method [MBK19]. There are many competing approaches for the backward propagation method. Since a full review is out of the scope of this introduction, we refer the reader to [Zha+11] and [NKM20] for a comparison of the various methods. If the backward propagation method is defined by the transmitter/receiverpath, it depends strongly on the number of sensors and their position [Zha+11]. In other cases, the resolution of these approaches is directly linked to the physical properties of the propagating waves, for example, the wavelength [MBK19]. Most tomography algorithms show only a projection on the surface of the waveguide and do not show any depth of the defect.
As an alternative to the above imaging algorithms, the localization and characterization of defects can be viewed as an inverse problem. The algorithm tries to find the set of parameters of a damage model that yield the smallest distance of the simulated data to the measurement data. These approaches consist of three main components. First, a forward model that generates simulated data representing the damage for the current set of parameters. Second, one uses an objective function that compares the measurement and simulation. Third, an optimization algorithm updates the parameter values to minimize the objective function, i.e., obtain a better fit of simulation and measurement. With regard to the simulation method, damage model, target function and optimization, there are different approaches that are briefly discussed in the following paragraphs.
For a discrete damage model, the extended finite element method (XFEM) is often used, because of its possibility to define mesh independent cracks [RGV07; RGV09; JT14; ACB18; LKZ18]. Rabinovich et al. [RGV09] used the XFEM in the frequency domain to reconstruct a single straight crack. The objective function is based on the least square error between measurement and simulation data. This objective function oscillates and has several local minima. A genetic algorithm optimizes the objective function to overcome the difficulty of getting trapped in a local minimum. The same authors improved the results by switching to the time-domain [RGV09] and changing the objective function. The objective function in the time-domain is based on the "arrival time" of the reflection, which is often used in classic defect detection in an ultrasonic test. Using the arrival time, the objective function oscillates less, but there are still local minima. Livani et al. [LKZ18] investigated multiple defects with particle swarm optimization. Finally, three dimensional problems with multiple cracks are solved by Agathos et al. [ACB18]. They used a two-step GA approach for elliptical cracks. Jung et al. [JT14] consider curved cracks and inclusions together with a gradient descent as the optimization method. Multiple optimization runs are performed with different starting points representing different cracks to circumvent local minima.
A time-reversal approach can also be used to optimize defect parameters. If a sensor signal is reversed in time and transmitted back in the domain, it will be refocused at the source. Amitt et al. [AGT14] defined a objective function based on these refocusing of the wave in a specific part of a membrane 1 . Unfortunately, this objective function is oscillating [AGT14]. A full scan of the parameter space that defined the crack and its position was performed for optimization. [Lop+17] considers an experimental validation of this approach.
Seidle et al. [SR16] used another time-reversal approach for a membrane. Here, an integral damage model is used which lowers the stiffness for certain elements. The connection between XFEM and a stiffness-reduction can be found in [Bro21]. This approach, also known as Full Wavefield Inversion, allows many possible defect geometries and is not limiting their number, but the resolution of the defect depends on the computational grid. A highly refined computational grid can lead to very time-consuming calculations, especially for three-dimensional problems. A theoretical and experimental investigation of a similar approach can be found, e.g., in Rao et al. [RRF16;RRF17;Rao+20], where an acoustic inversion is used for elastic wave propagation to lower the computational demand. A change in both stiffness and density is considered in [Rao+20].
The approaches mentioned above exploit the surface of the waveguide, whereas also a crosssection view may serve as alternative description. A plane strain simplification can be used to derive a two-dimensional problem -see Figure 1. Here, it is important to stress that these alternative dimensions of the cross-sectional view enable the analysis of the depth of the defect and the mode conversion due to defects with a certain depth. Wu et al. [WLH02] analyzed composite plates with the strip element method. The model allows investigating horizontal and vertical cracks in the cross-section. Semi-analytical methods such as the strip element method have the advantage that the forward calculation is efficient. This simulation approach is combined with a genetic algorithm to perform the optimization.
Gravenkamp proposed another semi-analytical method, the Scaled Boundary Finite Element Method (SBFEM), for simulating plate-like structures in the context of NDT and SHM [Gra+12]. The SBFEM was first developed by Song and Wolf to compute bounded and unbounded domains with the same approach [SW00; WS00]. The domain is divided into super elements. Quasi-polar coordinates, called scaled boundary coordinates, construct this first type of super element. For problems in the frequency domain, the super element has only degrees of freedom on the boundary, which leads to a reduction of one dimension in the computational cost. Later, Gravenkamp et al. used another coordinate transformation to construct super elements with a constant cross-section for bounded and unbounded domains in two-dimensions [GSP12;GBS15]. These super elements are very efficient for simulating waveguides. The approach was extended to curved waveguides, threedimensional-waveguides and waveguides with fluid/structure interaction [KG17; KGB17;Was+18]. In addition to the simulation of unbounded domains, a discrete crack model is an essential development in SBFEM [SON18]. Quite recently, the SBFEM was used in a shape optimization scheme for a horn speaker and meta-materials [Kha+21] illustrating the flexibility of this formulation in the meshing process and and with respect to infinite boundary conditions.
In this contribution, we propose a general method for reconstructing a single discrete defect model inside a cross-section model of a waveguide. The SBFEM is used as a efficient forward model for the simulation. The defect reconstruction is formulated as an optimization problem and a derivative-based optimization algorithm is used to solve it in a systematic and efficient way. The technique of Algorithmic Differentiation (AD, [GW08;Nau11]) is used to provide the required gradient information exactly and with low computational costs. We analyse the resulting approach with respect to properties of the objective function, the reconstruction quality and the robustness against noise.
The remainder of this paper is structured in the following way: The forward model and the physical assumptions are provided in Section 2. Section 3 summarizes the optimization scheme. This includes a discussion about the objective function and an estimate for the initial value.

SBFEM forward model
This section introduces the forward model to simulate the guided waves. Throughout, we consider only the two-dimensional case. The extension to three dimensions is subject of our current research. As mentioned in the introduction, the damage is incorporated directly into the waveguide model such that specific parameters model the damage. The forward model incorporates all important aspects of an ultrasound test on a waveguide. However, some aspects can only be designed for a specific experiment. Hence, certain simplifying assumptions have to be made to achieve greater general validity.
The first of these assumptions is that the defects are in a specific area within the waveguide modeled as an infinite domain. It is assumed that the vertical edges are far enough away from this area of interest so that the reflections from the vertical edges of the waveguide do not interfere with the reflections from the damage.
The second assumption is that the excitation by a sensor can be modeled by a traction forcě f (x, t) on one part of the boundary. This traction is mapped into the frequency domain by the Laplace-Transformation. In contrast to the Fourier-Transformation, the imaginary part of the Laplace-Transformation should weaken the resonance and wrap-around effects of the numerical model. In the discrete-time version, this procedure, known as Exponential Window Method (EWM) [Kau17], is given byf where c.c. denotes the complex conjugated of the previous term. The complex angular frequency ω n = (nω ∆ − i ζ) ∈ C is defined by a uniform frequency step ω ∆ and a small parameter ζ. Note that only the real part of the complex angular frequencies changes. Algorithmically, the EWM is implemented in the following three steps: (i) Multiplyτ (t) with the window function exp(−ζt) with ζ = 0.5ω ∆ , where the frequency step ω ∆ results from the time vector, use the Fast-Fourier-Transform (FFT) on the product, and search for the relevant frequencies.
(ii) Approximate the displacement spectrum for all relevant complex frequencies ω n (see Equation (3), (4)). (iii) For evaluation in the time-domain, transform the displacement by the inverse FFT and afterwards multiply the inverse window function exp(ζt).
For all examples, the time dependent partτ (t) of the excitation is a sin-modulated Gaussian pulse, i.e.,τ with the center frequency f c = 200 kHz and a time shift of t c = 5f −1 c . Figure 2a shows the pulse in the time-domain, while 2b illustrates the spectrum. Additionally, the dispersion curves of the waveguide are presented [GSP12]. This figure shows that the relevant spectrum of the excitation lies in the area, where only the two fundamental modes are propagating.  For each discrete angular frequency step, the SBFEM is used to approximate the displacement in the equations of linear elastodynamics with the displacement u, the linear stress σ σ σ, the density ρ and n is the outer normal vector. As for many other numerical methods, when using the SBFEM, the PDE in Equation (3) is converted into systems of linear equations with the global, dynamic stiffness matrix S ωn that depends on the complex angular frequency, the nodal displacement vector u ωn and the nodal traction vector f ωn . To compute the global dynamic stiffness matrix S ωn , the domain has to be sub-divided into super elements, i.e., Ω = K k=1 Ω k . Figure 3 shows such a subdivision into super elements. For each super element Ω k , parts of the boundary are meshed with finite elements that are used to compute a local stiffness matrix S k ωn . These local stiffness matrices are then assembled into the global matrix, while the global nodal traction vector is assembled directly from the finite elements. The finite element approximation uses spectral shape functions based on the Gauss-Lobatto points. Figure 3 depicts elements with shape functions of degree p = 4 with 5 nodes (marked with •). The corner points of the finite elements are dependent on a set of parameters q i that is later modified by the optimization process.
The prismatic super elements are used to approximate the extended parts of the waveguide. Figure 3 shows examples with the super elements Ω 1 , Ω 3 , Ω 4 and Ω 6 . These super elements require a constant cross-section. The cross-section parallel to the y-axis is then meshed with finite elements and scaled in the x-direction. If the super element is finite, as Ω 3 in Figure 3, the finite element meshes must coincide on both sides. Note that the length in x-direction does not influence the computational cost because there are only nodes at the cross-sections. The low number of nodes leads to a very efficient calculation for long parts of the domain. The outer boundaries parallel to the x-direction, illustrated by thin lines in Figure 3, have traction-free conditions. Such super elements can also approximate semi-infinite parts of the model that are used to bypass reflections, see, e.g., Ω 1 and Ω 6 in Figure 3. For the derivation of the local stiffness matrix of the super element, see [GBS15;Gra18]. The main part of the computation to set up the local stiffness matrix is an eigenvalue problem. A direct formula [ATM07] computes the derivative of the eigenvalue problem inside the differentiated algorithm when applying AD.
The super element that forms a star-convex polygon is always required when the prismatic super element can no longer define the geometry or boundary conditions. Here, a finite element line mesh is scaled toward a single point. The single point is called the scaling center. In general, the position of the scaling center is also dependent on the parameter set q i to be determined by the optimization. The super element only needs nodes on the boundary yielding a high level of effectiveness. An additional special feature is that the super element also provides a simple crack model because a double node can introduce it. The crack is illustrated by Ω 4 in Figure 3, where a red node marks the double node. The red line is a traction-free boundary due to the double node. Despite the crack tip singularity, one still observes exponential convergence under prefinement [BG19; TG18; BGB19] using such super elements. For the derivation of the local stiffness matrix of the super element, we followed [Che+14]. The derivation of the local stiffness matrix is quite involved. A continued-fraction-based approach builds on one algebraic Riccati equation and several algebraic Lyapunov equations [SW97;BS08]. For this work, both algebraic equations are solved by an eigenvalue decomposition inside the SBFEM algorithm. The solution for the Riccati equation can be found in [SW97], while the Lyapunov equation is solved in the appendix A.
It is worth mentioning that previous investigations have shown that polygonal super elements are quite robust for small and large finite element sizes inside the same mesh [XST18]. This robustness is important if geometric optimization changes the element sizes on the fly.

The Adapted Optimization Approach The Inverse Problem
In general, the procedure described here to identify the defect properties is an indirect method. The aim is to solve the inverse problem where y meas is given by physical measurements performed to detect the defect and F(q) : q → y sim denotes the Forward Operator. The forward operator is given by a simulation of the same quantity assuming that the defect is characterized by some vector q. Here, the forward operator is a short notation for the EWM together with the SBFEM-model described in the previous section. The space Q and its elements q ∈ Q represent one possible parameterization of the considered defect.
For the examples considered in this article the space Q is scaled either to [1, 2] 2 or [1, 2] 3 . This scaling to the same range of values for all parameters supports the solution of the inverse problem.
Solving the optimization problem (5), one obtains a parameterization q min for the defect which best approximates the measurement data.
In practice, inverse problems are usually hard to solve for several reasons. Firstly, the resulting optimization problem need not (and generally does not) have a unique solution which may not depend continuously on the measurement data y meas . Such inverse problems are called ill-posed problems and are pervailant throughout almost all applications of parameter identification problems. In order to solve ill-posed problems one is usually well-advised to deploy custom-tailored regularization techniques.
On the other hand, a different issue that may arise is that the error function in (5) may be hard to optimize as it may have adverse optimization properties such as not being sufficiently smooth or incorporating many only locally optimal points. In this case, it may be advantageous to modify the forward operator F and the measurements y meas so that the objective function is better suited for optimization purposes.
A simple approach uses the response data measured directly at a point. The y-component is more interesting as the x-component because of its practical relevance. For example, a single laser Doppler vibrometer measures the y-component of displacement or velocity, depending on the type of measurement method. We will assume that the displacement is measured. However, since both quantities have very similar characteristics, all results are transferable to the velocity. For all numerical examples, the first geometric parameter is the global position of the defect. Figure 4a shows in gray the y-displacement response of the model in Figure 3 at the point P for the global defect x-position q 1 , where the other parameters, crack depth and angle, are kept constant at 1.5. Additionally the envelope of the signal is shown in black or dark green. Figure 4b depicts the objective function based on the response, where the green signal in Figure 4a is used as measurement data y meas = F(q * ). However, as the waveguide is excited via a sinusoidal pulse, the measured response also incorporates the oscillating behavior. In fact, it can be seen that the objective function (5) also inherits these features, see Figure 5a for an illustration, which is turned zoom of the blue box in Figure 4b. We note that the resulting objective function has many only locally optimal points near the global optima. Hence, many optimization methods would have difficulties solving this optimization problem to global minimality.
The authors also experimented with different variations of the data such as using the spectrum of the response. However, the problem caused by many local minima remained. Many test revealed (a) Response and envelope at point P .  that the objective function with the best properties from an optimization point of view was the calculation of the envelope of the response via a Hilbert-transform, see Figure 4c and 5b. Analog to the response, Figure 4c shows the objective function for all crack positions between the green envelope in Figure 4b with a forward operator which included the calculation of the envelope. In addition to the global minimum, Figure 4c shows two local minima. These local minima are associated with the mode conversion. The excitation, to the left of the blue line, is a pulse that corresponds to the S0-mode. However, the S0-mode and the mode converted A0-mode are reflected. If the S0-wavepackage of the forward simulation is at the position of A0-wavepackage of the measurement data, this leads to a local minimum. Compare the green line with the line above in Figure 4a. Our approach to circumvent the local minima is to get a sufficient initial guess, which lies near the global minimum. The details of the initial guess can be found below.
For the rest of the paper, the measurement data is the envelope of the y-displacement at a single point, while the forward operator includes the calculation of the envelope.   Figures 4b and 4c, respectively.

Algorithmic Differentiation
In order to use gradient-based optimization methods for the reconstruction of defects, derivatives of solutions of the partial differential equation are required. These could be approximated via finite differences (e.g. [JT14]) or computed by setting up and solving the adjoint equation of the partial differential equation considered (e.g. [SR16;RRF16]). We tested the use of finite differences and found the approximation to be insufficient for the optimization procedure. Furthermore, in the development phase of this project we experimented with different boundary conditions. Hence, an adaptation of the adjoint equation as a time consuming and error-prone process would have been necessary in order to compute derivatives via the adjoint approach. Therefore, we apply algorithmic differentiation, also called automatic differentiation (AD), which provides derivative information of arbitrary order for a given code segment. The derivatives are provably accurate up to working accuracy due to the fact that the computer program evaluating the function is broken down into a sequence of elementary evaluations upon which the chain rule of calculus is systematically applied. There are numerous AD tools available for a wide range of programming languages, see the web site www.autodiff.org for an overview of tools and references. Since our forward operator described in detail in the previous section is implemented in MATLAB, we apply the AD-tool ADiMat [Bis+02] to compute derivatives of solutions of the partial differential equations. Some mild modifications to the code were necessary such as computing derivatives of the Lyapunov equation which originally was not available in ADiMat.

The Resulting Optimization Procedure
Even though we made some effort to reformulate the objective function (5) so that it is easier to optimize, we cannot eliminate all issues entirely. Hence, in this subsection, we present the optimization algorithms that we deployed.
A first simple approach for the solution of the optimization problem (5) is to use standard optimization software naively. This can involve derivative-based approaches such as IPOPT [WB06], WORHP [BW12], or the optimization routine fmincon provided by Matlab. Alternatively, one may use derivative-free approaches such as coordinate descent, genetic optimization or other heuristic-based optimization methods like particle-swarm methods.
Derivative-based optimization algorithms have the great benefit of fast and provable convergence but generally only converge locally. With derivative-free methods often there is the hope associated that these methods converge globally. However, global convergence cannot be proven. Moreover, the whole optimization procedure may require thousands of function evaluations. The required number of evaluations is of particular significance as the simulation in our case involves solving a partial differential equation; thus one objective function evaluation is already quite expensive -solving the partial differential equations many times is simply impractical.
In light of these facts, the approach presented here is a two-step approach: We first generate a refined initial guess for all variables. After a good initial point is established, we deploy a gradientbased optimization with fast convergence properties.
Initial guess for the global position As mentioned above, the first geometric parameter q 1 is the global position of the defect for all numerical examples. To get an initial guess for the defect position q 1 , we consider the time of flight to get a physical position L 1 inside the waveguide. This physical position is then mapped linearly into the parameter space. The physical position of the defect L 1 can be estimated by where T ∆ is the time of flight, c ex g is the group velocity at the center frequency of the excited mode, and c re g is the group velocity of the reflected mode with the highest amplitude. For all our examples prior investigations have shown that the highest reflections for the y-component are always associated with the A0-mode (see Fig. 4a), so the reflection velocity c re g is the A0-group velocity c A0 g at the center frequency -see Figure 2b. The computation of the time of flight T ∆ is based on cross-correlating the excitation pulse, which is separated by a defined threshold T ex from the reflections, with the rest of the signal. If the envelope of the cross-correlation is maximal, the argument is an approximation of the time of flight, i.e., T ∆ = arg max t r, where u y is the y-displacement, env is the envelope function, is the cross-correlation operator and [·] is the Iverson bracket 2 . Figure 4a shows the threshold as a blue line, and the excitation pulse is on the left side while the rest of the signal is on the right side. In the presented examples, the function r has a distinct maximum, which is associated with the A0-mode. However, if there is any ambiguity in the mode or the maximum value, it is simple to run the inverse method with multiple start values for the first parameter.
Refined initial guess for all parameters In order to obtain a refined initial guess for the gradient-based optimization, we can use the approximation Equation (6). The approximation for the first coordinate given there is fairly accurate and is located within the valley of the objective function, see Figure 5b. However, with a simple line-search in the vicinity of the valley, the accuracy of the initial guess can be greatly improved. This simple line-search requires very few function evaluations, and the computational cost is insignificant. The accuracy of the first coordinate is important as inaccurate values lead to locally optimal values (especially in the presence of noise) for the other coordinates with which optimization methods have difficulty escaping. The other coordinates are determined by selecting the vector with the smallest objective function from M random vectors, where our default value of M is 10. As the line-search is computationally very cheap, we also refine the other coordinates in the same fashion, beginning with the least critical or problematic coordinates. Due to many only locally optimal points, see Figure 11a, the amount of randomly selected points is increased to 100 for corrosion-type defects. The whole procedure is sketched in Algorithm 1.

Optimization procedure
We experimented with different optimization methods in order to solve the optimization problem Equation (5). One difficulty we encountered was that the objective function has many locally optimal points that are difficult to escape from. However, optimization approaches aiming at a more global scale, such as genetic algorithms, or particle swarm methods are infeasible due to the computation effort caused by one evaluation of the forward operator.
Algorithm 1 Refined initial guess 1: Evaluate Equation (6) for an initial guess of the first coordinate. 2: Sample M random points q 1 0 , . . . , q M 0 ∈ Q for the remaining second and possibly third coordinate. 3: Evaluate corresponding objective values F(q m 0 ) − y meas 2 , m = 1, . . . , M , and select q m 0 with the smallest objective value. 4: Apply coordinate-wise line-search to q m 0 to decrease the objective function value. The resulting parameter is q 0 .
We were able overcome these issues by applying an Iteratively Regularized Gauss-Newton Method (IRGNM), see e.g. [KNS08]. This method modifies the well known Gauss-Newton approach for nonlinear least-squares fitting problems by adding terms that make sure that matrices occurring in the Newton-step are invertible and adds a Thikonov-type term so that the overall solution stays near the initial guess. The overall procedure is given in Algorithm 2, where F is the Jacobian matrix of the forward operator, I is the Identity matrix and (·) H denotes the adjoint operator.

Numerical Experiments
In this section, we present several numerical examples to show the versatility of the approach. Each model considers another type of defect. The first example is about classic crack identification. The second example shows a similar set-up but tries to find a horizontal crack or rather a delamination defect. The third type of defect is a model of a corrosion defect.
All waveguides in the following examples are made out of steel. The isotropic material parameters are given in Table 1, and for all waveguides a plane strain is assumed. All meshes use shape functions of degree p = 4, i.e., elements with five nodes. The necessary degree is determined in advance by investigating the critical geometric parameter. For example, the smallest and the largest crack length are the critical geometric parameter for the first model. For all our models, a simulation with a degree of p = 4 leads to a difference in the signals below 0.1% compared to a simulation with p = 5. It is worth noting that there are numerous convergence studies for the SBFEM [TG18; GSZ20; BG19; Che+14] and the code is thoroughly validated using commercial FEM tools.     Unit global x-position of Ω 5 q 1 200 2200 mm x-position of the crack tip with 0 as the midpoint of Ω 5 q 2 −0.5 0.5 mm y-position of the crack tip with 0 as the midpoint of Ω 5 q 3 −2.25 0 mm

Waveguide with a Crack
This example considers a steel plate with a single straight crack. Figure 3 shows an overview of the model. Note that the x-axis is interrupted to fit the model into the figure. The first section of this figure is for the excitation. Arrows in Figure 3 show the area and direction where spatially constant traction is applied. At these parts, the excitation takes place by a normal traction force The traction is a model for the excitation of a double transducer set-up [SY04]. A double transducer set-up can produce both fundamental modes by applying a force in the same or opposite direction on the plate. For this traction, the S0-mode is excited, so the group velocity c ex g in equation (6) is S0-group velocity c S0 g . Here, the S0-excitation is appropriate because the S0-mode leads to larger reflections for these kinds of cracks.
The second section of Figure 3 shows the evaluation point P . The envelope of the y-displacement is evaluated as a time series. As described above, the y-component is of practical relevance as it can be measured using a single laser Doppler vibrometer. Figure 4 shows the y-displacement at point P in gray and the envelope as a black curve over the signal of different crack positions.
The third section contains the defect. A double node is inserted into a continued-fraction-based super element. The emerging inner traction-free boundary is marked in red in Figure 3. Note that the continued-fraction approach handles the singular stress at the end of this crack model. In total, the model has 180 degrees of freedom. This low number is possible because the semi-analytical approach approximates the large undamaged part of the waveguide. At this center frequency, the wavelength for the S0 mode is around 25 mm. With other simulation methods, this leads to a much higher computational effort.
location. The first parameter q 1 defines the global position of the crack. This change in position is achieved by shifting the super element Ω 5 along the x-axis and changing the length of the super element Ω 4 accordingly. The second and third parameters control the crack length and angle by moving the crack tip inside the super element Ω 5 relative to the boundary. The maximal and minimal values for the parameters are listed in Table 2. The crack length is between 0.25 mm and ca. 2.6 mm. Here the parameters are restricted to allow a numerically stable solution. Similar parameter restrictions can be found in other publications like [ACB18]. The objective function dependent on q 1 and q 2 with regard to the optimization problem (5) is displayed in Fig. 6a. Here, q 3 is kept constant at 1.5, i.e., only orthogonal cracks are considered. The artificial target measurement corresponds to the midpoint of the parameter space. The surface is relatively smooth and has only a few local minima. The global minimum at (1.5, 1.5, 1.5) T has clearly the lowest value. Using the result of Algorithm 1 as refined initial guess, the defect properties can be obtained by applying the IRGN method Algorithm 2 without any need for using a more in-depth initial guess search such as applying genetic algorithms.
Another important aspect of an inverse method, especially in view of the measurement data, is its performance under noise. To analyze the performance, we create artificial measurement data by a single forward simulation with a parameter vector q * and add random noise to this target signal. The noise is defined in terms of the excitation pulse. The same threshold as for the time of flight is used to identify the excitation pulse. In Figure 4, the blue line shows this threshold. The excitation pulse is to the left of the blue line. A noise level of 10 −2 and 10 −3 means that a normal-distributed noise with vanishing mean value a standard deviation of 1% and 0.1% of the maximum amplitude of the excitation pulse has been added to the target signal, respectively. The envelope of the target signal is calculated afterward. Figure 6b shows the reconstruction error |q min i − q * i | per parameter for the midpoint of the parameter space. Here, the first parameter, the global position of the crack, shows the smallest error. A lower error in the first parameter is also expected because the range in the physical space is much larger compared to the other parameters -see Table 2. Translating the error in the physical space, the error stays under 1 mm for a noise level under 3%, which is a satisfactory level of precision for most applications. The error in the second parameter, on the other hand, shows lower robustness to noise. This parameter, which is linked to the crack angle, is likely to produce incorrect results. However, the more critical third parameter, which is more related to the crack length, has a significantly lower reconstruction error. In general, this investigation is only valid for the current target parameter q * ; as the crack length gets smaller with the third parameter, it is also expected that error increases due to the negative influence of the noise. In conclusion, there are also physical boundaries to an appropriate parameter space depending on the noise inside the experimental data. Figure 7a shows the convergence of the reconstruction error with a small noise level of 10 −5 . This noise level is introduced to prevent an "Inverse Crime" as both the target data and the reconstruction algorithm use the same simulation method. We observe an asymptotic rate that coincides with the rate of the regulation parameter in the Algorithm 2.
To further validate the inverse method, ten different target parameters are randomly drawn. Noise with a standard deviation of 10 −5 is added to the corresponding signals, and the parameters are reconstructed. Table 3 lists the reconstruction error for these points. For all reconstructions, the error is below 1%.

Waveguide with a Delamination
This example considers a plate with a single straight horizontal delamination model. Figure 8 shows an overview of the model. Note that once again, the x-axis is interrupted at some places to fit the model into the figure. The two colors inside the figure indicate that these defects are often present between two different layers. However, the material of these layers in our example is the same steel. Prior investigations have shown that the A0-excitation will lead to larger reflections for delaminations. The traction is chosen accordingly as where e 2 is the second unit vector. Hence, c ex g in equation (6) is the A0-group velocity c A0 g at the center frequency -see Figure 2b.
The second section in Figure 8 highlights again the evaluation point P , where the envelope of the y-component of the displacement is taken as the data.
nodes in two super elements. The semi-analytical solution of the SBFEM once again captures the points with singular stress. There are two parameters for the optimization. The first parameter q 1 denotes the global position of the delamination. The second parameter changes the length of the delamination. The maximal and minimal values for the parameters are listed in Table 4. Figure 9a depicts the objective function dependent on q 1 and q 2 . As in the previous example, the target signal is generated with the midpoint of the parameter space for q * . Figure 9a and Figure 6a show a similar shape, but the two side channels with local minima are not present in this example. These side channels are missing because delaminations in the center of a waveguide do not lead to mode conversion for the A0-mode.
For this example, we could set the regularization parameter α n = 0 and achieve convergence within five iteration steps. In Figure 9b, the error for reconstruction with noise is given. Both parameters show high robustness against the noise. Additionally, Table 5 lists the reconstruction error for ten different randomly drawn target parameters with a noise level of 10 −5 . The error indicates a high level of precision.
In total, the model has 226 degrees of freedom. A single forward calculation is performed in under 4 seconds on a modern desktop computer, while the reconstruction is computed in under 5 minutes, including all preprocessing steps. The time is averaged over 100 reconstructions. description parameter min max Unit global x-position of Ω 5 , . . . , Ω 9 q 1 200 2200 mm corrosion length on the surface q 2 7.5 57.5 mm corrosion depth q 3 0.25 1.00 mm Table 6: Parameters for the waveguide with a corrosion defect.

Waveguide with a Corrosion Defect
(a) Objective function.  The last example considers a steel plate with a simple model for a corrosion damage. The assumption is that parts of the steel are missing, and a tapering has formed in the waveguide. Figure 10 shows an overview of the model. Note that again the x-axis is interrupted. The traction f S0 is the same as in first example (Equation (9)), which leads to an S0-mode excitation. The S0-mode shows greater sensitivity than the A0-mode to changes in thickness. Therefore, the S0excitation is used. The two re-entrant corners [BG19] in sub-domains Ω 5 and Ω 9 are located at the scaling center to handle the possible singular stresses for a deep corrosion. The model has 210 degrees of freedom in total.
There are three parameters to be optimized. Once more, the first parameter defines the global position of the defect. The second parameter defines the length of the corrosion by changing the size  Table 7: Error for reconstruction q min for differently shaped defects q * ∈ [1, 2] 3 for the waveguide with a corrosion defect with a noise level of 10 −5 .
of the super element Ω 7 . The third parameter denotes the depth of the corrosion -see Figure 10. The maximal and minimal values for the parameters are listed in Table 6. Figure 11a shows the objective function for a constant corrosion depth. The surface is more complicated due to the two main reflection points at the scaling centers, leading to mode conversion. As mentioned earlier, we had to increase the number of random vectors drawn for the initial guess because of the more complicated shape of the objective function. In Figure 11b the error for reconstruction with noise is given. Additionally, Figure 7b shows the convergence behavior for the reconstruction error at the midpoint of the parameter space. Despite the more complicated behavior of the objective function, the parameters show robustness against noise. For all ten tests, the reconstruction leads to small errors as listed in Table 7.

Conclusion
In the paper, we showed that a semi-analytical waveguide model given by SBFEM can be combined successfully with derivative-based optimization to reconstruct defects of different nature in a twodimensional waveguide, even in the presence of noise. For this purpose, SBFEM is implemented in MATLAB, coupled with the AD tool ADiMat and embedded in an iterative solution procedure. The reconstruction can be performed on a standard desktop computer within a very moderate time such that our approach offers considerable runtime advantages compared to alternative methods. Three numerical tests were conducted to illustrate our approach: The first example covers a crack identification, the second one a delamination, and the third one a corrosion. In all three cases, the waveguides were made out of steel. In the cross-sectional models, the reconstruction with a minimal amount of data is possible. Only the envelope of the y-displacement at a single point is sufficient if a certain geometry of a single defect is assumed. Due to the limited number of local minima, the envelope curve shows the most promising features for an optimization. In all tests carried out, the approach demonstrates robustness against noise and varying defect parameters. Future work will be dedicated to the extension to the three-dimensional case and anisotropic material behavior. The future work will allow the use of the proposed approach for non-destructive testing and structural health monitoring in real-life applications like carbon-enforced structures.

A Implementation Details
This appendix shows how the Lyapunov equation can be solved with the help of an eigenvalue decomposition. The Lyapunov equation is Define the right eigenvalue decomposition of A with the right eigenvector matrix V and the diagonal eigenvalue matrix D. Note that this leads to the left eigenvalue decomposition of A H After pre-multiplying with V −1 and post-multiplying with V −H , Equation (11) leads to a new Lyapunov equation with a diagonal coefficient matrix The matrixX can be computed element-wise by The solution matrix X of the original Lyapunov equation can be derived by the inverse of Equation (15), i.e., X = V HX V.